Hyper Suprime-Cam view of the CMASS galaxy sample. Halo mass as a function of stellar mass, size, and Sérsic index

November 14, 2018

An HSC view of the CMASS galaxy sample. Halo mass as a

function of stellar mass, size and Sérsic index.

Alessandro Sonnenfeld

, Wenting Wang

, and Neta Bahcall

1 _{Leiden Observatory, Leiden University, Niels Bohrweg 2, 2333 CA Leiden, the Netherlands}

e-mail: sonnenfeld@strw.leidenuniv.nl

2 _{Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan}

3 _{Department of Astrophysical Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ 08544, USA}

ABSTRACT

Aims.We wish to determine the distribution of dark matter halo masses as a function of the stellar mass and the stellar mass profile, for massive galaxies in the BOSS CMASS sample.

Methods. We use grizy photometry from HSC to obtain Sérsic fits and stellar masses of CMASS galaxies for which HSC weak lensing data is available, visually selected to have spheroidal morphology. We apply a cut in stellar mass, log M∗/M> 11.0, selecting

∼ 10, 000 objects. Using a Bayesian hierarchical inference method, we first investigate the distribution of Sérsic index and size as a function of stellar mass. Then, making use of shear measurements from HSC, we measure the distribution of halo mass as a function of stellar mass, size and Sérsic index.

Results.Our data reveals a steep stellar mass-size relation Re ∝ M βR

∗ , with βRlarger than unity, and a positive correlation between

Sérsic index and stellar mass: n ∝ M0.46

∗ . Halo mass scales approximately with the 1.7 power of the stellar mass. We do not find

evidence for an additional dependence of halo mass on size or Sérsic index at fixed stellar mass.

Conclusions.Our results disfavour galaxy evolution models that predict significant differences in the size growth efficiency of galaxies

living in low and high mass halos.

Key words. Galaxies: elliptical and lenticular, cD – Gravitational lensing: weak – Galaxies: fundamental parameters

1. Introduction

The size evolution of massive quiescent galaxies is one of the open problems in cosmology. On the one hand, observations show that, at fixed stellar mass, the average size of quiescent galaxies has increased by a factor of a few between z ≈ 2 and the present (e.g. Daddi et al. 2005; Trujillo et al. 2006; van Dokkum et al. 2008). On the other hand, theoretical models have struggled to reproduce this trend (see e.g. Hopkins et al. 2010a). Qualita-tively, the size evolution of massive quiescent galaxies is under-stood to be the result of mergers. Minor mergers, in particular, are known to be an efficient mechanism to increase the size of a galaxy (Naab et al. 2009). However, it is not clear whether the observed merger rates are sufficient to explain the size evolution signal, especially in the redshift range 1 < z < 2 (Newman et al. 2012). Moreover, it is very difficult to preserve the tightness of the stellar mass-size relation observed at z ∼ 0 with evolution-ary models based exclusively on dissipationless (dry) mergers (Nipoti et al. 2012).

Clues on the physical mechanisms that are relevant for the size evolution of massive galaxies can be obtained by study-ing the dependence of size on environment (see Shankar et al. 2014b). Although the correlation between environment and galaxy size at fixed stellar mass has been investigated by a rel-atively large number of studies, no clear consensus has been reached, with some works showing evidence for a positive cor-relation Cooper et al. (2012); Papovich et al. (2012); Lani et al. (2013); Yoon et al. (2017); Huang et al. (2018a), and others

? _{Marie Skłodowska-Curie Fellow}

finding results consistent with no dependence Huertas-Company et al. (2013b); Newman et al. (2014); Allen et al. (2015); Dam-janov et al. (2015); Saracco et al. (2017). The signal, if present, is in any case small: reported differences between the average size of galaxies found in clusters and in less massive associations are on the order of 20%.

A possible way forward is looking at correlations between the observed properties of galaxies and those of their host dark matter halos. The evolution of the mass of a dark matter halo traces directly its accretion and merger history. This, in turn, should relate directly to the size evolution of its central galaxy, since mergers are thought to be the main driver of the growth in size, at least at the massive end of the galaxy distribution.

Halo masses can be measured using weak gravitational lens-ing. However, the low strength of the weak lensing signal around typical massive galaxies requires the statistical combination of measurements over hundreds or thousands of lenses, making it challenging to obtain an accurate description of the distribution of halo mass as a function of galaxy properties other than stellar mass.

Recently, Charlton et al. (2017) reported the detection of a positive correlation between halo mass and size at fixed stellar mass, using stacked weak lensing measurements from the CFHTLenS survey (Heymans et al. 2012). Here we use a Bayesian hierarchical inference method to fit for the distribution of halo mass as a function of stellar mass, size and Sérsic in-dex, using photometric and weak lensing data from the Hyper Suprime-Cam (HSC) survey (Aihara et al. 2018a).

We draw the lens galaxies used in our study from the CMASS sample of the Baryon Oscillation Spectroscopic Survey (BOSS Schlegel et al. 2009; Dawson et al. 2013). The CMASS sample has been used for a variety of studies in galaxy evolution and cosmology (Anderson et al. 2014; Beutler et al. 2014; More et al. 2015; Montero-Dorta et al. 2016; Leauthaud et al. 2017; Tinker et al. 2017; Favole et al. 2018, e.g.). Thanks to the im-age quality (0.600typical i−band seeing) and depth of HSC data, our work will allow us to obtain a more accurate description of the light distribution in CMASS galaxies, compared to previ-ous studies. A by-product of our analysis is then a measurement of the stellar mass-size relation of CMASS galaxies, which we generalise to a stellar mass-size-Sérsic index relation, based on unprecedented deep and sharp data.

The structure of this work is the following. In Section 2 we describe the data used for our study. In Section 3 we carry out an investigation on the distribution of Sérsic index and half-light radius as a function of stellar mass. In Section 4 we perform a weak lensing analysis to infer the distribution of dark matter halo mass as a function of various galaxy properties. We discuss our results in Section 5 and conclude in Section 6. We assume a flat ΛCDM cosmology with ΩM= 0.3 and H0= 70 km s−1Mpc−1.

2. Data

2.1. The CMASS sample

The set of galaxies subject of our study is drawn from the CMASS sample of BOSS, which is part of the Sloan Digital Sky Survey III (SDSS-III Eisenstein et al. 2011). The CMASS (’constant mass’) sample has been constructed with the primary goal of selecting luminous red galaxies in the redshift range 0.43 < z < 0.7 for the measurement of the baryon acoustic oscil-lation signal. CMASS galaxies are selected by applying a series of cuts on colours and magnitudes, as measured from SDSS pho-tometry. We refer to Reid et al. (2016) for details on the selection criteria of the CMASS sample.

For each galaxy in the CMASS sample, we have a spectro-scopic redshift measured with the BOSS spectrograph.

2.2. HSC photometry

We select CMASS galaxies that lie in the region covered by the first-year shear catalogue of the HSC survey (Mandelbaum et al. 2018). This corresponds to 136.9 deg2_{imaged in grizy bands at}

the full survey depth (i ∼ 26.4 mag 5σ limit for a point source, Aihara et al. 2018b), resulting in a sample size of approximately 18,000 CMASS galaxies. For each galaxy in the sample, we ob-tain small cutouts of coadded images in each band from the S17a internal data release. Data reduction, including the generation of coadded images, model point spread function (PSF) kernels and sky subtraction is performed using hscPipe (Bosch et al. 2018), a version of the Large Synoptic Survey Telescope software stack (Ivezic et al. 2008; Axelrod et al. 2010; Juri´c et al. 2015). We refer to Bosch et al. (2018) for details on the data reduction pro-cess.

In Figure 1, we show colour-composite images of a set of 20 galaxies drawn randomly from the final sample used for our study (the criteria used to define the final sample will be listed in the rest of this Section).

2.3. Sérsic profile fits

We fit an elliptical Sérsic surface brightness profile (Sersic 1968) to the photometric data in each band:

I(x, y)= I0exp        −b(n) R Re !1/n       , (1)

where x and y are Cartesian coordinates aligned with the major and minor axis of the elliptical isophotes and origin in the centre, Ris the circularised radius,

R2≡ qx2+y

2

q, (2)

nis the Sérsic index, and b(n) is a numerical constant that en-sures that the light enclosed within the isophote with R= Reis

half of the total light (see Ciotti & Bertin 1999).

For each galaxy, we fix the parameters of the Sérsic profile to be the same in all five bands, with the exception of the amplitude I0. In other words, we assume galaxies to have spatially uniform

colours. We run a Markov Chain Monte Carlo (MCMC) to sam-ple the parameter space defined by the galaxy centroid, half-light radius, axis ratio, position angle, Sérsic index, and g−i, r−i, i−z, y − icolours. For each set of values of these structural parame-ters, we find the i−band amplitude that minimises the following χ2_, χ2₌X λ X j         I(mod)_{λ, j} − I_{λ, j}(obs) σ2 j         2 , (3)

where I_{λ, j}(mod)is the PSF-convolved Sérsic model in band λ evalu-ated at pixel j, Iλ, jis the observed surface brightness at the cor-responding pixel and σλ, jis the observational uncertainty. Since the model surface brightness is linear in the i−band amplitude, the above χ2_{can be minimised analytically.}

A circular region of 50 pixel (8.4”) radius centred on each galaxy is used for the fit. This corresponds to a physical radius of 54 kpc at z = 0.55, the median redshift of our sample. The choice of this aperture is a compromise between the need for capturing as much light as possible from each galaxy and at the same time avoiding to select too big of a region, in order not to be affected by systematics in the sky subtraction. We verified that our results do not change if we use a fitting region a factor of two larger or smaller.

The model parameters are initialised to the values of the cmodel_dev fit produced by hscPipe, with n = 4. We run SEx-tractor (Bertin & Arnouts 1996), with a 2σ detection threshold, on the i and g band image to mask objects not associated with the main galaxy. We run a preliminary MCMC on the centroid and colour parameters while keeping the Sérsic index, position angle, axis ratio and half-light radius fixed, run SExtractor once again on the model-subtracted image, then run another MCMC over the full set of parameters. We use the package emcee (Foreman-Mackey et al. 2013) to perform the MCMC. Typical statistical uncertainties are 2% on Reand n, 1% on the axis ratio q, ∼ 1 deg

on the position angle and 0.01 mag on the magnitudes.

Fig. 1. Colour-composite images in HSC-irg bands of a set of 20 galaxies drawn randomly from the sample used for our study. We used the algorithm of Marshall et al. (2016) to create the images.

out photometric measurements in an automated way, such as blended objects, images with significant contamination from a nearby bright star, or strong gravitational lenses. After this vi-sual inspection step, we are left with ∼ 13, 000 elliptical galaxies with clean photometry.

2.4. Stellar masses

We estimate stellar masses by fitting synthetic stellar population models to the observed grizy magnitudes, following the proce-dure described by Auger et al. (2009), with minor modifications. We take simple stellar population (SSPs, or instantaneous-burst) models, produced by Bruzual & Charlot (2003) using semi-empirical stellar spectra from the BaSeL 3.1 library (Westera et al. 2002) and the Padova 1994 evolutionary tracks, assuming a Chabrier initial mass function (IMF, Chabrier 2003). We then use the bc03 code (Bruzual & Charlot 2003) to create composite stellar population (CSP) model spectra, assuming an exponen-tially decaying star formation history, over a four-dimensional grid of age (i.e. time since the initial burst), star formation rate decay time τ, metallicity Z and dust attenuation τV. We use these

synthetic spectra to evaluate the model broadband flux in each HSC filter at each point of the grid, adding redshift as a fifth di-mension. Finally, we fit for stellar mass, age, τ, metallicity and dust attenuation by running an MCMC: for each set of values of the model parameters, we obtain the predicted grizy magnitudes by interpolating over the grid, and compare them to the observed values.

-1.0 0.0 log τV 2 4 6 8 A ge (G yr ) 1 2 τ (G yr ) -2.0 -1.5 lo g Z 11.0 11.2 11.4 log M∗ -1.0 0.0 lo g τV 2 4 6 8 Age (Gyr) 1 2 τ(Gyr) -2.0 -1.5 log Z

Fig. 2. Posterior probability distribution of the stellar population param-eters of an example galaxy. The three different contour levels mark the 68%, 95% and 99.7% enclosed probability regions.

with stellar mass, as indicated by the observational study of Gal-lazzi et al. (2005). In Figure 2, we show the posterior probability distribution of the stellar population parameters for an example galaxy.

2.5. Weak lensing data and cuts

We take weak lensing measurements from the first-year shape catalogue of the HSC survey (Mandelbaum et al. 2018). The source number density of this catalogue is 24.6 arcmin−2

(un-weighted). Shapes are measured on coadded i−band images us-ing the re-Gaussianization PSF correction method of Hirata & Seljak (2003). We refer to Mandelbaum et al. (2018) for details about the shape measurement process and the properties of the shape catalogue.

We use photometric redshifts (photo-zs) obtained with the photo-z code Mizuki (Tanaka 2015) on the data release 1 of the HSC survey (Tanaka et al. 2018).

The main weak lensing analysis is carried out using the Bayesian hierarchical inference method of Sonnenfeld & Leau-thaud (2018, SL18 from now on). This method is based on the isolated lens assumption: each lens galaxy in the sample is as-sumed to be at the centre of its dark matter halo, with no other mass component affecting the shape of the background sources around it. This assumption breaks down for satellite galaxies. It is therefore important to select a sample of lenses with the lowest possible satellite fraction. For this purpose, we first apply a cut in stellar mass, selecting only galaxies with observed stellar mass M∗(obs), defined as the median of the posterior probability

distri-bution in M∗, larger than 1011M. 90% of the CMASS galaxies

left in our sample after the visual inspection step satisfy this re-quirement. We then match our sample with the HSC cluster cat-alogue of (Oguri et al. 2018), and remove objects with a cluster membership probability larger than 50% that are not identified as the brightest cluster galaxy. This step removes ∼ 5% of the objects in the sample.

Following the SL18 method, we model the weak lensing sig-nal produced by each lens on a set of background sources within a cone of a given radius, which we take to be 300 physical kpc at the redshift of the lens, where the isolated lens assumption is more realistic (i.e. where the so-called 1−halo term dominates).

101 ₁₀2 ₁₀3 R (kpc) 100 101 102 103 ∆ Σ (M ⊙ p c − 2) 11.0≥ log M∗< 11.3 11.3≥ log M∗< 11.6 11.6≥ log M∗< 11.9 log M∗≥ 11.9

Fig. 3. Excess surface mass density profile in different bins of observed stellar mass, obtained by stacking the weak lensing signal in circular annuli.

One of the foundations of the SL18 method is the ability of treat-ing the likelihood of the weak lenstreat-ing data around the lenses as independent from each other, which simplifies the problem greatly. However, in case of lenses in close proximity along the line of sight, the lensing cones of different lenses can overlap. In principle, one should simultaneously model the contribution of each lens on all the sources within the overlapping cones. In practice, this is computationally difficult. We avoid this prob-lem by eliminating overlapping pairs, using the same procedure adopted by SL18, as follows. We rank the lenses in decreasing order of observed stellar mass, loop through them, and remove from the sample any lens the lensing cone of which overlaps with that of a more massive galaxy. The idea is that, in case of lenses in close line of sight proximity, the main lensing signal should come from the most massive galaxy. Although in reality the problem is more complex (for instance, the strength of the lensing signal depends not only on the lens mass but also on its redshift), the fraction of objects removed with this procedure is only 10%, and we do not expect this step to introduce any sig-nificant bias to our weak lensing analysis (as already shown by SL18 on mock observations). We then reach our final sample size of 10, 403 CMASS galaxies.

To get a sense of the quality of the weak lensing data avail-able for our study, we carry out a stacked analysis. We make four bins in stellar mass and measure the stacked excess surface mass density ∆Σ in different radial bins, using the HSC weak lens-ing pipeline. We plot the∆Σ profiles in Figure 3. Differences in the lensing signal in different bins are detected with high confi-dence: this suggests that the data should allow us to determine not only an average halo mass of the sample, but also how halo mass scales with stellar mass and, possibly, with size. In any case, Figure 3 is only shown for illustration purposes: all our weak lensing analysis is carried out by forward-modelling the weak lensing signal of individual halos.

2.6. The final sample

11.0 11.2 11.4 11.6 11.8 12.0 12.2 12.4 log M_∗ 0.0 0.5 1.0 1.5 2.0 log Re (kp c) 2000 4000 N 0.0 0.2 0.4 0.6 0.8 1.0 1.2 log n 0.4 0.6 0.8 redshift 500 1000 1500 N S´ersic SerExp de Vauc.

Fig. 4. Distribution in half-light radius, Sérsic index and redshift as a function of stellar mass of the sample of 10, 403 CMASS galaxies sub-ject of our study. Contour levels mark the 68%, 95% and 99.7% en-closed probability regions. Blue contours and histograms refer to the fiducial model, consisting of a single Sérsic surface brightness profile for each galaxy. The distribution in stellar mass and size obtained with the SerExp model, described in subsection 2.7, is plotted in green, while that obtained with the de Vaucouleurs model is marked by solid lines.

also a trend of increasing Sérsic index with stellar mass. In the next Section we will quantify the strength of these correlations.

Figure 5 shows the distribution of half-light radii in angular units. For most objects, the half-light radius is well within the 8.4” region used for the fit, with only 1% of the galaxies in the sample exceeding this limit.

It is also interesting to check how far out from the centre of each galaxy the HSC data allows us to probe. For this pur-pose, we compute, for each object, the radius at which the best-fit i−band Sérsic surface brightness falls below the level of the root mean square fluctuation from the sky background, Rsky. We

plot the distribution of Rsky/Rein Figure 6. The median of this

distribution is 3.6, with only 4% of the objects having a value of Rsky/Resmaller than unity. For the typical galaxy in our

sam-ple, then, HSC i−band data allows us to detect flux out to 3.6

101 102 103 N 0 2 4 6 8 10 12 Re(′′) 0.0 0.2 0.4 0.6 0.8 1.0 N (< Re ( ′′)) /N to t

Fig. 5. Top: Distribution in half-light radius, in angular units, of the galaxies in our sample. Bottom: cumulative distribution in half-light ra-dius.

times the half-light radius. For a de Vaucouleurs profile (n = 4 de Vaucouleurs 1948), the fraction of the total mass enclosed within this aperture is around 80%. This means that ∼ 20% (or 0.08 dex) of the total flux of estimated for a typical galaxy is not directly observed, but rather the result of an extrapolation.

Around 5% of the lenses are identified as brightest central galaxies (BCGs) in the HSC cluster catalogue of Oguri et al. (2018). For these objects, we expect intra-cluster light to con-tribute with a non-negligible fraction to the observed surface brightness. We do not attempt to separate the light of the cen-tral galaxy from the intra-cluster light, as the distinction be-tween these two components is not very well defined. For cluster BCGs, then, our measurements of the surface brightness profile, and the stellar mass and size derived from it, must be interpreted as that of the sum of the central galaxy and the intra-cluster light.

2.7. Alternative models

100 101 102 103 N 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 Rsky/Re 0.0 0.2 0.4 0.6 0.8 1.0 N (< Rsky /R e )/ Nto t

Fig. 6. Top: distribution in the ratio between the radius at which the i−band surface brightness of a galaxy falls below the sky fluctuation level, Rsky, and the half light radius Re. Bottom: cumulative distribution.

degrees of freedom. Bernardi et al. (2014) showed that a SerExp model provides a better description of the surface brightness pro-file of SDSS galaxies, compared to models with a single Sérsic component.

In Figure 7, we show the circularised i−band surface bright-ness profiles as well as the enclosed flux as a function of aper-ture for 20 example galaxies (the same objects plotted in Fig-ure 1), obtained for the Sérsic, the de Vaucouleurs and the Ser-Exp model. In the region probed by HSC data, roughly where the galaxy is brighter than the sky fluctuation level (marked by the horizontal dashed line in each subplot), the differences between the Sérsic (blue curves) and the SerExp model (green curves) are minimal. This is reassuring, as it means that the Sérsic model is sufficiently flexible to capture the complexity in the surface brightness profiles of our sample.

A similar case can be made, for many objects, for the de Vaucouleurs model, although for some systems the data clearly prefers Sérsic indices different from n = 4. These are the objects for which the best-fit de Vaucouleurs profile, shown as a dotted black curve, starts to deviate from the best-fit Sérsic profile well within the region above the background noise level. Typically, the residuals between the best-fit de Vaucouleurs profile and the data show under-subtraction at large radii, for the systems with n> 4, and over-subtraction for the systems with n < 4.

The difference between the large R behaviour of the three models results in different total fluxes and, consequently, differ-ent half-light radii. For most of the objects in our sample, our data does not allow us to decide which of the three models is more accurate, especially between the Sérsic and the SerExp: the differences arise at large radii, where the surface brightness is too faint to be detected by HSC data. We will then carry out the analysis with all models, and check whether the results are stable with respect to the particular assumption on the large R behaviour of the surface brightness profile of our galaxies.

Although the de Vaucouleurs model produces bad fits for a non-negligible fraction of the CMASS galaxies in our sample, we still show results obtained with this particular surface bright-ness profile choice, since many studies in the literature are based on the same model.

3. Generalised mass-size relation

3.1. Basic model

We fit for the mass-size relation of our sample of massive (log M(obs)∗ > 11) elliptical galaxies using a Bayesian

hierarchi-cal approach: we model the distribution of true stellar mass and size of the sample with a functional form, described by a set of hyper-parameters η, which we then infer from the data. We choose the following form for this distribution:

P(M∗, Re|η) = S(M∗)R(Re|M∗). (4)

Here S(M∗) is a skew Gaussian in log M∗:

S(M∗)= 1 √ 2πσ∗ exp ( −(log M∗−µ∗) 2 2σ2 ∗ ) Φ(log M∗), (5) with Φ(log M∗)= 1 + erf s∗ log M∗−µ∗ √ 2σ∗ ! , (6)

and R(Re|M∗) is a Gaussian in log Re

R(Re|M∗)= 1 √ 2πσR exp      −(log Re−µR(M∗)) 2 2σ2 R      , (7)

with a mean that scales with stellar mass as follows

µR(M∗)= µR,0+ βR(log M∗− 11.4). (8)

The hyper-parameters describing this 2-dimensional distribution are:

η ≡ {µ∗, σ∗, s∗, µR,0, σR, βR}. (9)

µ∗, σ∗and s∗describe the stellar mass distribution. If the

skew-ness parameter s∗ is set to 0, this reduces to a Gaussian with

mean µ∗and σ∗. Positive values of s∗correspond to a

distribu-tion with a long tail above µ∗and a sharper cutoff at low masses.

The parameter µR,0is the average value of log Reat the pivot

stel-lar mass log M∗= 11.4, σ∗is the scatter in log Reat fixed stellar

mass, and finally βRis a power-law dependence of Reon M∗.

Our goal is to infer the posterior probability distribution of the set of hyper-parameters η given the data d. Using Bayes’ theorem,

P(η|d) ∝ P(η)P(d|η), (10)

where P(η) is the prior on the hyper-parameters and P(d|η) is the likelihood of observing the data given the value of the hyper-parameters. Assuming measurements on separate galaxies are independent from each other, the latter is the following product over individual objects:

P(d|η)=Y

i

Z

dM∗,idRe,iP(di|M∗,i, Re,i)P(M∗,i, Re,i|η). (11)

28 26 24 22 I (R ) [m ag ar cs ec − 2] n = 1.79 20.5 20.0 19.5 19.0 i( < R ) [m ag ] n = 6.84 n = 5.12 n = 5.33 n = 4.91 28 26 24 22 I (R ) [m ag ar cs ec − 2] n = 8.62 20.5 20.0 19.5 19.0 i( < R ) [m ag ] n = 8.10 n = 2.05 n = 7.74 n = 4.83 28 26 24 22 I (R ) [m ag ar cs ec − 2] n = 4.80 20.5 20.0 19.5 19.0 i( < R ) [m ag ] n = 2.96 n = 4.88 n = 5.78 n = 4.97 28 26 24 22 I (R ) [m ag ar cs ec − 2] n = 4.27 10−1 ₁₀0 ₁₀1 R (arcsec) 20.5 20.0 19.5 19.0 i( < R ) [m ag ] n = 4.93 10−1 ₁₀0 ₁₀1 R (arcsec) n = 7.59 10−1 ₁₀0 ₁₀1 R (arcsec) n = 9.01 10−1 ₁₀0 ₁₀1 R (arcsec) n = 2.88 10−1 ₁₀0 ₁₀1 R (arcsec)

Fig. 7. i−band surface brightness profile (top sub-panel) and enclosed flux as a function of aperture (bottom sub-panel) for the best-fit Sérsic (blue), SerExp (green), and de Vaucouleurs (dotted black) profile for the 20 example objects shown in Figure 1. The dashed horizontal line marks the surface brightness corresponding to the rms sky fluctuation. The dashed vertical lines mark the half-light radius.

of each galaxy. These are necessary to calculate the likelihood. We are not interested in constraining the stellar mass and size of individual galaxies, so we marginalise over all possible values, modulated by P(M∗,i, Re,i|η), which is the distribution introduced

in Equation 4.

The data consists of the observed values of the stellar mass and the half-light radius, M∗(obs)and R(obs)e , and related

uncertain-ties. For simplicity, we neglect the uncertainty on size, since it is on the order of 2%. As a result, the covariance between the mea-surement of the half-light radius and the stellar mass is also set to zero. This is a fair approximation, since the uncertainty on the stellar mass is dominated by systematic uncertainties in the stel-lar population synthesis model. The first term in the integrand of Equation 11 then becomes

P(di|M∗,i, Re,i)= P(M_∗,i(obs)|M∗,i)δ(R(obs)e,i − Re,i), (12)

and Equation 11 reduces to P(d|η)=Y

i

Z

dM∗,iP(M(obs)∗,i |M∗,i)P(M∗,i, Re,i|η). (13)

We approximate the likelihood of obtaining a value of the ob-served stellar mass M_∗,i(obs)given the true value M∗,ias a Gaussian in log M∗,i:

P(M_∗,i(obs)|M∗,i)= √Ai 2π∗,i exp       

−(log M∗,i− log M

(obs) ∗,i )2 2∗,i        . (14)

In the equation above, ∗,i is the standard deviation in log M∗,i,

as obtained from the MCMC chain of the stellar population fit, and Aiis a normalisation constant that ensures that the integral

S´ersic SerExp de Vauc. 1.0 1.2 1.4 βR 0.11 0.13 0.15 σR 0.8 0.85 µR,0 1.0 1.2 1.4 βR 0.11 0.13 0.15 σR log Re∼ N (µR,0+ βR(log M∗− 11.4), σR)

Fig. 8. Posterior probability distribution of the parameters describing the distribution of galaxy sizes as a function of stellar mass, introduced in subsection 3.1. Results based on the Sérsic, SerExp and de Vau-couleurs model are shown.

we applied the cut log M∗(obs)> 11.0 to define our sample, this is

equivalent to

Z ∞ 11.0

dlog M(obs)_∗,i √Ai 2π∗,i exp       

−(log M∗,i− log M

(obs) ∗,i )2 2∗,i        = 1. (15) We sample the posterior probability distribution of the hyper-parameters η given the data d, Equation 10, by running an MCMC. We calculate the integrals in Equation 13 with the importance sampling and Monte Carlo integration method de-scribed in SL18. We assume flat priors on all hyper-parameters except s∗, for which we assume a flat prior on its 10-base

loga-rithm. In Figure 8 we plot the posterior probability distribution relative to the parameters describing the distribution of sizes: µR,0, σR and βR. The median and 68% credible interval of all

hyper-parameters is reported in Table 1. We show results based on the Sérsic, SerExp and de Vaucouleurs models.

We find a remarkably steep mass-size relation, although the exact value of the dependence of Reon M∗depends on the choice

of the photometric model: we obtain βR = 1.37 ± 0.01 if the

Sér-sic photometric model is used and βR = 1.22±0.01 in the SerExp

case. The difference between the models appears to be driven by differences in size at the high mass end of the distribution: at fixed stellar mass, the values of Re obtained with the SerExp

model are slightly smaller compared to the Sérsic values, as can be seen from the bottom-left panel of Figure 4. A similar be-haviour was found by Bernardi et al. (2014) on SDSS galaxies. As a result, the uncertainty on our inference on the slope of the mass-size relation is dominated by systematics related to the par-ticular choice of the surface brightness profile used to describe CMASS galaxies.

The difference between the Sérsic and the de Vaucouleurs sizes is much more evident, hence the much lower value of the mass-size relation slope inferred for this particular model: βR= 0.98 ± 0.01. However, we believe this to be a biased

infer-ence, since a de Vaucouleurs model provides a poor description of the surface brightness profile of a non-negligible fraction of the objects in our sample.

3.2. Distribution of size and Sérsic index

We add complexity to our model by considering the distribution in Sérsic index of our galaxies, in addition to stellar mass and size. We model the distribution in M∗, n and Reas follows:

P(M∗, n, Re|η) = S(M∗)I(n|M∗)R(Re|M∗, n). (16)

The term S(M∗) is the same introduced in Equation 5 and

Equa-tion 6. I(n|M∗) is a Gaussian distribution in the base-10

loga-rithm of the Sérsic index, I(n|M∗)= 1 √ 2πσn exp ( −(log n − µn(M∗)) 2 2σ2 n ) , (17)

with a mean that scales with stellar mass as

µn(M∗)= µn,0+ βn(log M∗− 11.4) (18)

and dispersion σn. Finally, we keep the same form as Equation 7

for the term describing the distribution in half-light radius, R, but update Equation 8 by adding a dependence on the Sérsic index to the average size:

µR(M∗)= µR,0+ βR(log M∗− 11.4)+ νRlog n/4. (19)

At fixed stellar mass, Equation 16 is a bi-variate Gaussian in log n and log Re.

As we have done for the sizes, we assume that the Sérsic in-dices are measured exactly, since the uncertainties on n are very small. In Figure 9 we plot the posterior probability distribution of the hyper-parameters describing the distribution in Sérsic in-dex, µn,0, σnand βn, as well as the hyper-parameters describing

the distribution in size, including the new hyper-parameter νR,

which describes the dependence of size on n. The median and 68% credible values of all hyper-parameters are listed in Table 2. We show only result based on the Sérsic model, as n is not well defined in the case of a SerExp profile and is fixed in the case of a de Vaucouleurs model.

The inferred average value of log n at the pivot stellar mass log M∗= 11.4 is µn,0= 0.704 ± 0.002, corresponding to a Sérsic

index n ≈ 5.1. We find a positive correlation between Sérsic index and stellar mass, quantified by the parameter βn= 0.464 ±

0.009. This is closely related to the wellknown Sérsic index -luminosity correlation (Caon et al. 1993), and tells us that, for instance, the average Sérsic index of galaxies with stellar mass log M∗ = 11.8 is as high as n ≈ 7.8. We also find a correlation

between size and Sérsic index at fixed stellar mass: νR = 0.38 ±

0.01.

The value of parameter βR, the correlation between mass and

size, is smaller compared to the value obtained in the analysis of subsection 3.1, based on the simpler model (1.18 vs 1.37). This is because, in the context of this more complex model, βR

quantifies how size scales with stellar mass at fixed Sérsic index. The combination of this scaling with 1) the dependence of size on n (positive value of parameter νR) and 2) the positive

corre-lation between Sérsic index and stellar mass, produces the steep mass-size relation observed.

At the same time, the inferred value for the average size at the pivot point, µR,0, is smaller than the value obtained previously.

This is because µR,0now refers to the average log Reat a stellar

mass log M∗= 11.4, and Sérsic index n = 4, corresponding to a

de Vaucouleurs profile, which is a smaller value of n with respect to the average for galaxies of that mass.

Table 1. Mass-size relation, as modelled in subsection 3.1. Median values and 68% credible interval of the posterior probability distribution of individual hyper-parameters, marginalised over the rest of the hyper-parameters. Results based on the Sérsic, SerExp and de Vaucouleurs model are reported.

Sérsic model SerExp model de Vauc. model Parameter description

µR,0 0.855 ± 0.002 0.854 ± 0.002 0.774 ± 0.002 Average log Reat stellar mass log M∗= 11.4

σR 0.147 ± 0.002 0.129 ± 0.001 0.112 ± 0.001 Dispersion in log Rearound the average

βR 1.366 ± 0.011 1.218 ± 0.011 0.977 ± 0.011 Power-law dependence of size on M∗

µ∗ 11.249 ± 0.004 11.274 ± 0.005 11.252 ± 0.006 Mean-like parameter of Gaussian term in Equation 5

σ∗ 0.285 ± 0.004 0.254 ± 0.004 0.202 ± 0.004 Dispersion-like parameter of Gaussian term in Equation 5

log s∗ 0.44 ± 0.02 0.31 ± 0.02 0.17 ± 0.03 Log of the skewness parameter in Equation 6

0.4 0.5

β

0.12

0.14

σ

1.15

1.2

β

0.3

0.4

ν

0.69

0.71

µ

0.16

0.17

σ

0.810.82

µ

0.4

0.5

β

0.12 0.14

σ

1.15 1.2

β

0.3

0.4

ν

0.69 0.71

µ

0.16 0.17

σ

log R

∼ N (µ

+ β

(log M

− 11.4) + ν

log [n/4], σ

)

log n

_{∼ N (µ}

+ β

(log M

_{− 11.4), σ}

)

Fig. 9. Posterior probability distribution of the parameters describing the distribution of Sérsic index and galaxy size.

galaxies in the region probed by the data, which roughly cor-responds to the radial range 2kpc. R . 30kpc. The lower limit is set by the finite resolution of HSC data. An extrapolation of the best-fit Sérsic profile to the very inner regions produces very steep inner surface brightness profiles. Such cuspy profiles are typically not observed in the centres of nearby massive galax-ies, but cannot be ruled out in our observations, due to the at-mospheric blurring of the images. The impact of a potentially inaccurate description of the inner surface brightness profile on the stellar mass and size measurement used for our study is nev-ertheless minimal: the mass enclosed in the region significantly affected by the atmospheric seeing is only a small fraction of the total. However, we urge caution when using our derived Sérsic

fits to predict quantities that are more sensitive to the inner stellar distribution, such as the dynamical mass.

The outer limit to the region where the Sérsic profile is most accurate is set by the radius at which the surface brightness of CMASS galaxies falls below the noise level of HSC data. The median value of this radius for our sample is 28kpc.

4. Generalised stellar-to-halo mass relation

Table 2. Mass-Sérsic index-size relation, as modelled in subsection 3.2. Median values and 68% credible interval of the posterior probability distribution of individual hyper-parameters, marginalised over the rest of the hyper-parameters.

Parameter description

µn,0 0.704 ± 0.002 Average log n at stellar mass log M∗= 11.4

σn 0.163 ± 0.001 Dispersion in log n around the average

βn 0.464 ± 0.009 Power-law dependence of Sérsic index on M∗

µR,0 0.817 ± 0.002 Average log Reat stellar mass log M∗= 11.4 and Sérsic index n = 4

σR 0.133 ± 0.002 Dispersion in log Rearound the average

βR 1.184 ± 0.012 Power-law dependence of size on M∗

νR 0.383 ± 0.011 Power-law dependence of size on n

µ∗ 11.249 ± 0.004 Mean-like parameter of Gaussian term in Equation 5

σ∗ 0.285 ± 0.004 Dispersion-like parameter of Gaussian term in Equation 5

log s∗ 0.43 ± 0.02 Log of the skewness parameter in Equation 6

the next subsection, we let halo mass to scale with stellar mass and half-light radius. In subsection 4.2, we generalise this model to allow for a possible additional dependence of halo mass on Sérsic index.

Mirroring the approach adopted in Section 3, we will fit the simpler of the two models to data obtained with the three di ffer-ent parameterisations of the surface brightness profiles, in order to explore the impact of potential systematics associated with the photometry fitting on the derived stellar-to-halo mass relation. We will then fit the full model to our fiducial Sérsic profile-based measurements.

4.1. Halo mass as a function of stellar mass and size

We will start by considering the joint distribution in stellar mass M∗, half-light radius Reand halo mass Mh. Similarly to the

ap-proach adopted in the previous Section, we model this distribu-tion as follows:

P(M∗, Re, Mh)= S(M∗)R(Re|M∗)H (Mh|M∗, Re). (20)

The terms S and R are the same one introduced in subsection 3.1: a skew Gaussian in log M∗and a Gaussian in log Rewith a

mean that scales with stellar mass. The term relative to the halo mass, H , is modelled as a Gaussian in log Mh,

H (Mh|M∗, Re)= 1 √ 2πσh exp        −(log Mh−µh(M∗, Re)) 2 2σ2 h        , (21)

with a mean that scales with stellar mass and size as

µh(M∗, Re)= µh,0+βh(log M∗− 11.4)+ξh(log Re−µR(M∗)), (22)

and dispersion σh. The term µR(M∗) in Equation 22 is the

aver-age log Re for galaxies of mass M∗, introduced in Equation 8.

The parameter ξh, then, is a power-law dependence of halo

mass on ”excess size”, defined as the ratio between the half-light radius of a galaxy and the average size for galaxies of the same stellar mass. This excess size is conceptually similar to the “mass-normalised radius” used by Newman et al. (2012) and Huertas-Company et al. (2013a).

The full list of hyper-parameters of the model is then η ≡ {µh,0, σh, βh, ξh, µ∗, σ∗, s∗, µR,0, σR, βR}. (23)

In order to infer the posterior probability distribution of the hyper-parameters given the data, we need to be able to evaluate the likelihood P(d|η) for any value of η. Under the isolated lens assumption, on which the S18 weak lensing method is based, the

likelihood of the data relative to each lens is independent from each other. We can then write, analogously to Equation 11, P(d|η)=Y

i

P(di|η). (24)

For each galaxy, the data consists of the measured stellar mass, effective radius, and shape measurements of background sources within a cone of 300 physical kpc radius at the redshift of the lens (as previously discussed in subsection 2.5). In order to eval-uate the likelihood of the shape measurements, it is necessary to assume a model for the mass distribution of the lens. We first define the halo mass as the dark matter mass enclosed in a shell with average density equal to 200 times the critical density of the Universe, commonly referred to as M200. We then assume a

spherical Navarro Frenk and White (NFW Navarro et al. 1997) shapefor the dark matter halo density profile:

ρ(r) ∝ 1

r(1+ r/rs)2

. (25)

Finally, we assume a mass-concentration relation from Macciò et al. (2008): given r200 (the radius enclosing a mass equal to

M200), we assume that the concentration ch = r200/rs is drawn

from the following Gaussian in its base-10 logarithm: P(ch|Mh)= 1 √ 2πσc exp ( −(log ch−µc(Mh)) 2 2σc ) , (26)

with mean µcthat scales with halo mass as

µc(Mh)= µc,0+ βc(log Mh− 12) (27)

and dispersion σc. We fix µc,0 = 0.830, βc = −0.098 and σc =

0.10.

For each galaxy, the likelihood of observing the data given the value of the hyper-parameters is obtained by marginalising over all possible values of the concentration, as well as stellar and halo mass and half-light radius:

P(di|η) =

Z

dM∗,idMh,idch,idRe,iP(di|M∗,i, Mh,i, ch,i, Re,i)×

P(ch,i|Mh,i)P(M∗,i, Mh,i, Re,i|η). (28)

Here P(ch,i|Mh,i) is the mass-concentration relation introduced in Equation 26, which acts as a prior on the concentration pa-rameter. As done throughout Section 3, we set the uncertainty on the half-light radius to zero, so that the integral over Re

lens, given the values of M∗, Mh and ch, corrected by an

addi-tive and multiplicaaddi-tive bias on the measurements, obtained via simulations by Mandelbaum et al. (2018). We refer to SL18 and Sonnenfeld et al. (2018, subsection 3.4) for more details on the calculation of the likelihood term for weak lensing data in the SL18 formalism.

In Figure 10, we plot the posterior probability distribution of the hyper-parameters belonging to the halo mass term H in Equation 20. The median and 68% credible interval on the full set of hyper-parameters is given in Table 3. We show results based on the Sérsic, SerExp and de Vaucouleurs model.

The average log Mhfor galaxies of log M∗= 11.4 and

aver-age size for their mass is µh,0 = 12.79 ± 0.03. We then infer a

positive correlation between halo mass and stellar mass at fixed size, βh = 1.70 ± 0.10. However, we do not find any strong

evi-dence for an additional correlation between halo mass and size at fixed stellar mass. The median and 1σ limit of the parameter de-scribing the dependence on size is ξh= −0.14 ± 0.17, consistent

with zero.

These values are obtained using sizes and stellar masses based on the Sérsic photometric model, but the SerExp model produces consistent results, as can be seen in Figure 10. The main differences are that the SerExp model produces a steeper halo mass-stellar mass relation (larger value of parameter βh),

and that the inference on the halo mass-size correlation is broader and closer to zero (ξh= −0.03±0.21). The broadening of

the posterior probability distribution on parameter ξhis related to

the smaller intrinsic scatter around the stellar mass-size relation measured for the SerExp model (σR = 0.13 vs. σR = 0.15

ob-tained using Sérsic): a smaller intrinsic scatter means a shorter lever arm in the measurement of the secondary correlation be-tween Mhand Reat fixed M∗, which in turn results in less

preci-sion.

Using the photometric measurements obtained by imposing a de Vaucouleurs profile, we obtain an even steeper stellar-to-halo mass relation, βh = 2.04 ± 0.15, and an average halo mass

at the pivot stellar mass (parameter µh,0) ∼ 0.1 dex higher

com-pared to the Sérsic case. This difference is related to the differ-ent stellar mass distributions obtained with the two models: a de Vaucouleurs profile tends to under-estimate the stellar mass of the most massive CMASS galaxies, for which the data favours values of the Sérsic index n > 4. This flattening in the slope of the stellar-to-halo mass relation β when going from a de Vau-couleurs to a Sérsic model is a well known effect (Shankar et al. 2014a, 2018; Kravtsov et al. 2018, see e.g.).

In our model of the generalised stellar-to-halo mass relation, stellar mass and size are treated as the independent variables, while halo mass is the dependent variable. For the purpose of comparing our results with other studies, it can be useful to in-vert this relation, and obtain the distribution of size as a function of halo mass. This can be done with a posterior predictive proce-dure: we fix the value of the stellar mass to log M∗= 11.4, draw

samples of the hyper-parameters from the posterior probability distribution, then, for each value of the hyper-parameters, draw samples of halo mass and half-light radius for the corresponding model. In Figure 11, we plot the resulting distribution of Re as

a function of Mh. This can be interpreted as our inference of the

size-halo mass distribution of log M∗= 11.4 galaxies, assuming

our sample is complete in size at this value of the stellar mass. We can then summarise this distribution by fitting a power-law halo mass-size relation to it. We do this for each draw of the values of the hyper-parameters, obtaining Re∝ M_h−0.03±0.03,

cor-responding to the cyan band in Figure 11. This is an alternative

visualisation of our main result, consisting of a lack of a strong correlation between halo mass and size at fixed stellar mass.

4.2. Halo mass as a function of stellar mass, Sérsic index and size

We now add a dimension to our model: the Sérsic index. We modify the galaxy probability distribution Equation 20 by adding a term describing the distribution in n, as follows:

P(M∗, Re, Mh)= S(M∗)I(n|M∗)R(Re|M∗, n)H(Mh|M∗, n, Re).

(29) The term I(n|M∗) is the same introduced in Equation 17: a

Gaus-sian in the base-10 logarithm of the Sérsic index, with mean that scales with stellar mass, as parameterised in Equation 18. Then, as done in subsection 3.2, we let the mean of the size distribu-tion scale with Sérsic index, according to Equadistribu-tion 19. Finally, we modify the mean of the Gaussian distribution in log Mh, to

allow for a power-law dependence of halo mass on Sérsic index, as follows:

µh(M∗, n, Re)=µh,0+ βh(log M∗− 11.4)+ νh(log n − µn(M∗))+

ξh(log Re−µR(M∗, n)).

(30) In the above equation, µn(M∗) is the average base-10 logarithm

of the Sérsic index for galaxies of mass M∗, defined in

Equa-tion 18. The new parameter νh, then, describes how halo mass

scales with “excess Sérsic index”: the ratio between the value of nof a galaxy and the typical value of n for galaxies of the same stellar mass.

The motivation for allowing halo mass to vary with Sérsic index is the following: different values of n for galaxies of the same stellar mass point to different evolutionary histories, possi-bly related to the number and type of mergers experienced. If dif-ferences in evolutionary paths reflect in the halo mass, we should be able to detect a correlation between n and Mh.

We run an MCMC to sample the posterior probability dis-tribution of the Sérsic profile-based model, given the data. The inference on the hyper-parameters describing the halo mass dis-tribution is shown in Figure 12, while the median values and 68% credible regions of the inference on all the model hyper-parameters are listed in Table 4.

Results change little with respect to the inference based on the simpler model used in the previous subsection. Most notably, the inference on the new parameter νh, describing the correlation

between halo mass and Sérsic index, is consistent with zero: νh=

−0.11 ± 0.14. Our data, then, allows us to rule out any strong dependence of halo mass on Sérsic index.

5. Discussion

5.1. The role of observational scatter

The main result of our study is the measurement of the correla-tion, or lack thereof, between halo mass and galaxy size at fixed stellar mass. This is a non-trivial measurement to make: galax-ies are distributed along a relatively narrow mass-size relation, with a 0.15 dex intrinsic scatter in Reat fixed M∗, meaning that

any such signal is to be searched across a small dynamic range in size. Additionally, observational errors cause data points to move in the M∗− Replane from their true position. If not modelled,

S´ersic model

SerExp model

de Vauc. model

-0.5 0.0 0.5

ξ

0.2

0.3

0.4

σ

1.0

1.5

2.0

2.5

β

12.6 12.8 13.0

µ

-0.5

0.0

0.5

ξ

0.2

0.3

0.4

σ

1.0

1.5

2.0

2.5

β

log M

∼ N (µ

+ β

(log M

− 11.4) + ξ

(log R

− µ

(M

)), σ

)

Fig. 10. Posterior probability distribution of the parameters describing the distribution of halo mass, for the model introduced in subsection 4.1. In particular, µh,0is the average log Mh at a stellar mass log M∗ = 11.4 and average size for that mass, σhis the dispersion in log Mharound the

average, βhis a power-law scaling of halo mass with stellar mass, and ξhis a power-law scaling of halo mass with excess size (defined as the ratio

between the size of a galaxy and the average size of galaxies of the same stellar mass).

Table 3. Median values and 68% credible interval of the posterior probability distribution of the hyper-parameters describing the distribution in stellar mass, half-light radius and halo mass of the galaxies in our sample, as inferred by fitting the model introduced in subsection 4.1.

Sérsic model SerExp model de Vauc. model Parameter description

µh,0 12.79 ± 0.03 12.83 ± 0.03 12.91 ± 0.03 Average log Mhat stellar mass log M∗= 11.4 and average size

σh 0.35 ± 0.03 0.32 ± 0.03 0.37 ± 0.04 Dispersion in log Mharound the average

βh 1.70 ± 0.10 1.73 ± 0.11 2.04 ± 0.15 Power-law dependence of halo mass on M∗at fixed size

ξh −0.14 ± 0.17 −0.03 ± 0.21 0.00 ± 0.25 Power-law dependence of halo mass on excess size at fixed M∗

µR,0 0.854 ± 0.002 0.852 ± 0.002 0.775 ± 0.002 Average log Reat stellar mass log M∗= 11.4

σR 0.140 ± 0.002 0.124 ± 0.002 0.108 ± 0.001 Dispersion in log Rearound the average

βR 1.390 ± 0.012 1.239 ± 0.011 1.004 ± 0.011 Power-law dependence of size on M∗

µ∗ 11.246 ± 0.004 11.271 ± 0.004 11.249 ± 0.005 Mean-like parameter of Gaussian term in Equation 5

σ∗ 0.286 ± 0.004 0.256 ± 0.004 0.203 ± 0.004 Dispersion-like parameter of Gaussian term in Equation 5

log s∗ 0.47 ± 0.02 0.34 ± 0.02 0.20 ± 0.03 Log of the skewness parameter in Equation 6

stellar mass, larger size galaxies are on average more massive than their smaller size counterparts, making the interpretation of correlations measured directly on point estimates of observed quantities problematic (see subsection 2.2 of SL18 for a more detailed explanation).

The Bayesian hierarchical formalism on which this study is based allows us to forward model the effects of observational

scatter, and thus obtain an unbiased measurement of the true dis-tribution of the model parameters, provided that the estimates of the observational uncertainties are accurate.

We obtained uncertainties on the values of M∗by fitting

Table 4. Median values and 68% credible interval of the posterior probability distribution of the hyper-parameters describing the distribution in stellar mass, half-light radius, Sérsic index and halo mass of the galaxies in our sample, as inferred by fitting the model introduced in subsection 4.2 to the Sérsic profile-based measurements.

Parameter description

µh,0 12.78 ± 0.03 Average log Mhat stellar mass log M∗= 11.4 and average size

σh 0.37 ± 0.03 Dispersion in log Mharound the average

βh 1.69 ± 0.10 Power-law dependence of halo mass on M∗at fixed Sérsic index and size

νh −0.11 ± 0.14 Power-law dependence of halo mass on excess Sérsic index at fixed M∗

ξh −0.10 ± 0.19 Power-law dependence of halo mass on excess size at fixed M∗and n

µR,0 0.817 ± 0.002 Average log Reat stellar mass log M∗= 11.4 and Sérsic index n = 4

σR 0.128 ± 0.002 Dispersion in log Rearound the average

βR 1.214 ± 0.012 Power-law dependence of size on M∗

νR 0.369 ± 0.011 Power-law dependence of size on n

µn 0.704 ± 0.002 Average log n at stellar mass log M∗= 11.4

σn 0.162 ± 0.001 Dispersion in log n around the average

βn 0.467 ± 0.010 Power-law dependence of Sérsic index on M∗

µ∗ 11.246 ± 0.004 Mean-like parameter of Gaussian term in Equation 5

σ∗ 0.286 ± 0.004 Dispersion-like parameter of Gaussian term in Equation 5

log s∗ 0.46 ± 0.02 Log of the skewness parameter in Equation 6

11.5 12.0 12.5 13.0 13.5 14.0 log Mh 0.2 0.4 0.6 0.8 1.0 1.2 1.4 log Re (kp c) Re∝ M_h−0.03±0.03

Posterior predictive plot, log M_∗= 11.4 galaxies

Fig. 11. Posterior predictive distribution in half-light radius as a func-tion of halo mass for galaxies with log M∗= 11.4. This is obtained by

first drawing values of the hyper-parameters from the posterior prob-ability distribution of our Sérsic model-based inference, then drawing values of Reand Mhfor 1,000 galaxies, given the model specified by

the hyper-parameters. Blue contours correspond to 68% and 95% en-closed probability. The cyan band shows the 68% credible region of the size-halo mass relation, obtained by fitting a power-law relation to the mock distribution of Reand Mhat each draw of the hyper-parameters.

so to avoid unrealistically small errors (see subsection 2.4). As a result, the derived uncertainty on M∗is to some extent arbitrary,

as it reflects our choice on the amplitude of the systematic error added to the data. Moreover, the error on M∗is also set by the

priors on the stellar population parameters entering the fit, such as age or dust attenuation.

The median value of the uncertainty on M∗of the sample is

0.10 dex. In order to test the impact of making a more conser-vative choice on the estimate of the stellar mass uncertainty, we repeated the analysis of subsection 4.1 after adding in quadra-ture a further 0.10 dex systematic error on M∗. The inference on

the hyper-parameters describing the distribution in halo mass is plotted in red in Figure 13, on top of the original inference (in blue). For the sake of simplicity, we only show results based on

the Sérsic model. The inference on the halo mass-size correla-tion parameter ξh changes: with the larger observational

uncer-tainties on M∗, large negative values of ξh are allowed by the

data. However, strong positive correlations between Mh and Re

are still ruled out even with inflated error bars on M∗.

For the sake of completeness, we also show the results ob-tained by setting the uncertainties on M∗ to zero (solid lines in

Figure 13). In this case, the data seems to favour a positive cor-relation between halo mass and size. However, as discussed by SL18, this is a biased inference.

5.2. Comparison with other observational studies 5.2.1. The mass-size relation

A striking result of our study is the steep stellar mass-size re-lation: as shown in subsection 3.1, we find Re ∝ M∗1.37 when

a single Sérsic model is used to describe the surface brightness profile of CMASS galaxies (Re ∝ M∗1.22for the SerExp model).

By comparison, most studies of the mass-size relation in quies-cent galaxies find much shallower slopes, with βR ≈ 0.5 − 0.6

(see e.g. van der Wel et al. 2008; Damjanov et al. 2011; New-man et al. 2012; Huertas-Company et al. 2013a). We think the main reason for this discrepancy lies in the stellar mass distribu-tion of our sample. The median stellar mass is log M∗ = 11.45,

higher than that of the samples used in the aforementioned stud-ies. As shown by Bernardi et al. (2011), the mass-size relation is not a strict power-law at all masses, but becomes steeper above M∗ ∼ 2 × 1011, where most of the galaxies in our sample lie.

We are then probing a region of parameter space where the slope of the M∗− Re relation is steeper, compared to the value that

characterises quiescent galaxies at lower masses.

It is more difficult to reconcile our results with the study of Favole et al. (2018): in their analysis of the mass-size relation of CMASS galaxies, they find a slope βR ≈ 0.2, much shallower

S´ersic model

-0.5 0.0 0.5

ξ

0.2

0.3

0.4

σ

1.0

1.5

2.0

2.5

β

-0.5

0.0

0.5

ν

12.6 12.8 13.0

µ

-0.5

0.0

0.5

ξ

0.2 0.3 0.4

σ

1.0 1.5 2.0 2.5

β

-0.5 0.0 0.5

ν

log M

∼ N (µ

+ β

(log M

− 11.4) + ν

(log n

− µ

(M

)) + ξ

(log R

− µ

(M

, n)), σ

)

Fig. 12. Posterior probability distribution of the parameters describing the distribution of halo mass, for the model including a dependence of Mh

on the Sérsic index introduced in subsection 4.2. In particular, µh,0is the average log Mhat a stellar mass log M∗= 11.4, average n for that mass

and average size given M∗and n. σhis the dispersion in log Mharound the average, βhis a power-law scaling of halo mass with stellar mass, νhis

a power-law scaling of halo mass with excess Sérsic index (defined as the ratio between the Sérsic index of a galaxy and the average n of galaxies of the same stellar mass) and ξhis a power-law scaling of halo mass with excess size.

A key difference in the Favole et al. (2018) study is in the stellar mass measurements: while we fitted stellar population synthesis models to HSC photometry, they used stellar masses from the Portsmouth stellar mass catalogue (Maraston et al. 2013). The Portsmouth stellar masses are obtained from SDSS photometry, much noisier compared to HSC data. This could lead to a higher observational scatter that would flatten the ob-served M∗− Rerelation, if not accounted for.

5.2.2. The stellar-to-halo mass relation

The stellar-to-halo mass relation of CMASS galaxies was stud-ied by Tinker et al. (2017), using galaxy clustering and abun-dance matching, and taking advantage of the stellar mass com-pleteness study of Leauthaud et al. (2016). In Figure 14 we plot the average value of log Mhas a function of log M∗, as measured

by Tinker et al. (2017), on top of the same quantity obtained from our inference, calculated using Equation 22. The stellar masses used by Tinker et al. (2017) were obtained under the assumption of a Kroupa IMF, which we converted into Chabrier IMF-based stellar masses by applying a −0.05 dex shift.

The clustering-based measurement of Tinker et al. (2017) is systematically above our Sérsic profile-based results by ∼ 0.2 − 0.3 dex. This discrepancy is similar to the one observed by Leauthaud et al. (2017), in their comparison between the clus-tering and stacked weak lensing signal of the CMASS sample. Leauthaud et al. (2017) discussed various possible ways to solve this tension, including varying the value of the cosmological pa-rameter S8= σ8

√

Ωm/0.3, allowing for baryonic physics effects

S´ersic model

Inflated M

errors

No M

errors

-0.5 0.0 0.5

ξ

0.2

0.4

0.6

σ

1.2

1.4

1.6

1.8

β

12.6

12.8

µ

-0.5

0.0

0.5

ξ

0.2

0.4

0.6

σ

1.2 1.4 1.6 1.8

β

log M

∼ N (µ

+ β

(log M

− 11.4) + ξ

(log R

− µ

(M

)), σ

)

Fig. 13. Posterior probability distribution of the parameters describing the distribution of halo mass, based on the model of subsection 4.1. The blue contours show the inference obtained using the Sérsic model, as obtained in Section 4 and plotted in Figure 10. The distribution shown in red has been obtained by repeating the analysis of the Sérsic model and adding 0.1 dex in quadrature to the uncertainties on the observed stellar masses. The solid lines mark the distribution obtained by setting all uncertainties on the stellar masses to zero.

Given all these potential systematics affecting the compari-son between clustering and lensing, the discrepancy between our measurement and the Tinker et al. (2017) study is not particularly worrisome. We point out how there is very good agreement on the inference on both the slope of the stellar-to-halo mass rela-tion and the intrinsic scatter, between the two measurements.

5.2.3. Correlation between galaxy size and halo mass Many studies have looked into the correlation between galaxy sizes and the environment they live in, at fixed stellar mass, with somewhat conflicting results (see Section 1). Among these stud-ies, the work by Huang et al. (2018a) is particularly relevant for our analysis, since it is based on a sample of massive galaxies with photometric data from the HSC survey, selected for having a stellar mass enclosed within a radius of 100 kpc larger than 1011.6Mand a spectroscopic redshift in the range 0.3 < z < 0.5.

A good fraction of the galaxies in the Huang et al. (2018a) sam-ple also belong to our samsam-ple of CMASS galaxies.

Huang et al. (2018a) showed how, at fixed value of the stel-lar mass within 100 kpc in projection, M∗,100kpc, galaxies that are identified as the central of a massive cluster (log Mh & 14)

have preferentially smaller values of the stellar mass enclosed within 10 kpc, M∗,10kpc. Using M∗,100kpcas the fiducial value of the stellar mass of a galaxy, this result can be interpreted as the evidence for more extended galaxies (smaller value of M∗,10kpc

at fixed M∗,100kpc) living preferentially in more massive halos. A

recent stacked weak lensing analysis on the same set of objects confirmed the result (Huang et al. 2018b).

We wish to test whether this trend can be seen in our data as well. For this purpose, we perform a posterior predic-tive test: we use our model to generate mock data, and check whether this mock data reproduces the trend observed by Huang et al. (2018a). We proceed as follows: we take the maximum-likelihood values of the hyper-parameters describing the distri-bution in stellar mass, Sérsic index and size, as measured in sub-section 3.2, as well as the hyper-parameters describing the distri-bution in halo mass, as inferred in Section 4 for the Sérsic case. For the sake of a more straightforward interpretation of the re-sults, we set the value of the correlation between halo mass and size to zero, ξh = 0, which is consistent with our inference. We

draw a large set of values of M∗, n, Re, Mh, then compute the

11.0 11.2 11.4 11.6 11.8 12.0 log M∗ 12.0 12.5 13.0 13.5 14.0 log M h SerExp model de Vauc. model S´ersic model Tinker et al. 2017

Fig. 14. Average value of log Mh as a function of M∗. The blue band

and the green and dotted black pairs of lines show the 1σ confidence limit obtained from our analysis based on the Sérsic, SerExp and de Vaucouleurs model, respectively. The solid black line is a measurement from Tinker et al. (2017), obtained from galaxy clustering and abun-dance matching. Error bars indicate the intrinsic scatter in log Mh. The

inner (outer) ticks on the error bars relative to our measurements refer to the 84 (16) percentile of our inference on the scatter parameter σh.

The inference on the intrinsic scatter obtained with the de Vaucouleurs model is omitted, to avoid confusion.

Both our predicted M∗,100kpc and M∗ are model-dependent

quantities: we are only using data within 8.400 to constrain the surface brightness profile of our galaxies, corresponding to a projected physical aperture of 54 kpc at a redshift z 0.55. We are then extrapolating the Sérsic profile to obtain values of M∗,100kpc

and M∗.

In the bottom panel of Figure 15, we plot the values of M∗,10kpcas a function of M∗,100kpc. As done in Figure 9 of Huang et al. (2018a), we select objects with log M∗,100kpc > 11.6, split

the sample between galaxies that lie in halos more massive than log Mh = 14 (in red) and less massive (in grey), then fit the

M∗,10kpc− M∗,100kpcdistribution of each subsample with a linear

relation.

Qualitatively, we see a similar trend to that reported by Huang et al. (2018a): at fixed M∗,100kpc, galaxies at the centre of massive halos have preferentially lower values of M∗,10kpc. This

might seem in contradiction with the way our model was built, since we removed any dependence of the halo mass on size at fixed stellar mass. However, in the context of the Sérsic mod-els at the basis of our analysis, a non-negligible fraction of the total stellar mass of a galaxy comes from stars located beyond 100 kpc in projection, especially for very massive galaxies like the ones considered in this experiment. The fraction of stellar mass beyond 100 kpc is larger for galaxies with larger sizes: this implies that, at fixed M∗,100kpc, galaxies with larger sizes, or with

smaller values of M∗,10kpc, are on average more massive (have a larger value of M∗). As a result, these galaxies tend to live in

more massive dark matter halos, on average. This can verified in the middle panel of Figure 15, where the same data points of the plot in the bottom panel are colour-coded by the total stel-lar mass of each object: stelstel-lar mass increases with increasing M∗,100kpc, but also with decreasing M∗,10kpc. Halo mass follows a

similar trend.

In the top panel of Figure 15, we show a version of the bot-tom plot of the same figure, modified by adding a 0.1 dex random

observational scatter on the stellar mass measurements. This ob-servational scatter is meant to simulate errors in the stellar pop-ulation synthesis fitting, therefore they affect the measurements of the stellar mass at different radii in the same way (i.e. they shift the values of log M∗,100kpc and log M∗,10kpc by the exactly

same amount). The difference in the M∗,10kpc− M∗,100kpcrelation between clusters and lower-mass halos increases when looking at observed quantities, and is similar in amplitude to the signal measured by Huang et al. (2018a). This is related to the distor-tion of the mass-size reladistor-tion due to observadistor-tional scatter dis-cussed in subsection 5.1 and can be understood as follows: given a bin in observed stellar mass, smaller sized galaxies are sta-tistically more likely to have scattered into the bin from lower intrinsic stellar masses, and therefore to live in less massive dark matter halos, compared to larger sized ones (see also subsection 2.2 of SL18).

Both of the effects discussed above follow simply from the existence of a positive correlation between stellar mass and size and do not depend on a particular form of the size or mass-Sérsic index relation. Indeed, we have verified that making the same posterior predictive plot based on de Vaucouleurs fits pro-duces very similar results.

After this test, we conclude that our inference is consistent with the Huang et al. (2018a) analysis. In the context of our model, the correlation between halo mass and the observed val-ues of M∗,10kpcat fixed M∗,100kpccan be explained with the com-bination of two causes: massive galaxies having a non-negligible fraction of their stellar mass beyond 100 kpc, and the effect of observational scatter on the stellar mass measurements.

Charlton et al. (2017) measured the correlation between halo mass and galaxy size at fixed stellar mass, using stacked weak lensing measurements around a large set of galaxies. Charlton et al. (2017) made bins in luminosity and size, and looked at the variation in halo mass, as measured from weak lensing, as a function of size, after accounting for the dependence of halo mass on stellar mass. They measured a value ξh = 0.42 ± 0.12,

inconsistent with our inference. SL18 argued that part of their signal could be due to the effects of observational scatter. This would definitely be the case if their stacked weak lensing analy-sis was carried out in bins of observed stellar mass: as Figure 13 shows, ignoring observational errors on M∗can introduce a

sig-nal of similar magnitude to the reported value. However, Charl-ton et al. (2017) carry our their analysis in luminosity bins. Lu-minosity is measured with a much greater accuracy than stellar mass, and this should lead to a more accurate inference on the parameter ξh. Nevertheless, the Charlton et al. (2017) analysis

still relies on noisy stellar mass measurements from stellar pop-ulation synthesis modelling, and that should have some impact on the inference, if not taken into account. For a fair comparison, it would be important to test the effects of observational scatter in a study of the halo mass-size correlation à la Charlton et al. (2017). This, however, is beyond the scope of this paper.

5.3. Comparison with theoretical models