© ESO 2022

1. Introduction

Massive elliptical galaxies, with their smooth stellar distribution and limited presence of cold gas and dust, are relatively simple objects compared to their star-forming counterparts. Qualitatively speaking, the history of these objects consists of two phases: an initial period of intense star formation, followed by a later stage in which most of their growth is due to mergers with other galaxies (van Dokkum et al. 2010; Oser et al. 2010). However, various aspects of their evolution are still unknown at a quantitative level.

One of the biggest open questions in cosmology is what stopped massive galaxies from forming stars and what prevents them from resuming it. Feedback from the central supermassive black hole is believed to play an important role in galaxy quenching (i.e. the shutting off of star formation: see Fabian 2012, and references therein); however, despite recent progress (Cheung et al. 2016; Baron et al. 2018), demonstrating a causal connection between the two phenomena has proven to be difficult (see for example Shangguan et al. 2020).

Their post-quenching evolution is also a subject of debate: at fixed stellar mass, elliptical galaxies at z = 0 are a factor of a few larger, in terms of half-light radius, than their z = 2 analogues (Daddi et al. 2005; Trujillo et al. 2007; Buitrago et al. 2008). Such a size growth is generally attributed to dissipationless mergers (Naab et al. 2009), but both the observed and theoretically predicted merger rates appear too low to account for this phenomenon in its entirety (Newman et al. 2012; Nipoti et al. 2012). The addition of new systems to the quiescent galaxy population due to the quenching of previously star-forming galaxies can alleviate this problem if these recently quenched galaxies have larger sizes than the pre-existing ones (see e.g. Scarlata et al. 2007; Damjanov et al. 2019). Moreover, the presence of radial gradients in stellar population properties, which were not taken into account in several of the studies mentioned above, act in ways that lead to an overestimate of the size growth (Suess et al. 2019a,b). Nevertheless, we are still lacking a physical model able to reproduce the evolution of the population of massive galaxies in all its complexity.

A fundamental ingredient in the evolution of galaxies is their dark matter halo. The dark matter distribution determines the gas accretion rate and temperature in the early stages of star formation, as well as the merger rate with other dark matter halos and their respective galaxies. In particular, the mass of a dark matter halo and its evolution is the most important quantity that determines the properties of the galaxy at its centre (Moster et al. 2018; Behroozi et al. 2019). From a galaxy evolution perspective, having access to measurements of halo mass would then allow one to link the properties of galaxies with their assembly history. Direct constraints on the dark matter distribution can be used to build more accurate models of the relation between galaxies and halos, which in turn can reduce systematic uncertainties in many cosmological probes that rely on galaxies as tracers of the gravitational potential. We refer to Wechsler & Tinker (2018) for a more thorough discussion regarding this point.

One consequence of the connection between the distributions of dark and baryonic matter is the fact that galaxies with higher stellar mass are preferentially found at the centres of the more massive dark matter halos. This is because the stellar mass of a galaxy is related to the mass of the initial gas reservoir, which is in turn proportional to the dark matter mass. The exact mapping between stellar and dark matter halo mass, however, depends on the efficiency of the process of conversion of gas into stars. If this quantity varies as a function of the assembly history of a galaxy, additional non-trivial correlations between the properties of galaxies and halo mass could be introduced. For example, if the gas-to-star conversion efficiency varies with time, the ratio between the stellar and dark matter mass of a galaxy is different among galaxies that formed their stars at different epochs (see for example Moster et al. 2020).

Observationally, searches for correlations between halo mass and galaxy properties have been carried out on various fronts. There have been several efforts aimed at determining how the average halo mass of star-forming galaxies differs from that of quiescent ones at a fixed stellar mass. These include studies based on weak gravitational lensing (Hoekstra et al. 2005; Velander et al. 2014; Hudson et al. 2015; Mandelbaum et al. 2016; Bilicki et al. 2021; Wang et al. 2021), satellite kinematics (More et al. 2011; Lange et al. 2019), gas and globular cluster kinematics (Posti & Fall 2021), galaxy clustering (Zehavi et al. 2011; Cowley et al. 2019), and combinations of multiple such probes (Tinker et al. 2013). While the amplitude and significance of the signal varies among different studies, some of them clearly indicate that passive galaxies tend to be found at the centres of more massive dark matter halos, compared to star-forming galaxies of the same stellar mass. This suggests that halo mass, or more generally halo assembly history, plays a role in determining the quenching time of a galaxy (see Zu & Mandelbaum 2016).

Charlton et al. (2017), Sonnenfeld et al. (2019a), and Huang et al. (2020) have attempted to measure the dependence of halo mass on the stellar profile of massive quiescent galaxies at fixed stellar mass, using weak gravitational lensing observations. They found somewhat conflicting results: Charlton et al. (2017) claim a positive correlation between the halo mass and galaxy half-light radius at a fixed stellar mass, but their measurement is in tension with a study by Sonnenfeld et al. (2019a), which ruled out correlations as strong as those reported by Charlton et al. (2017). Using data, in large part overlapping with the Sonnenfeld et al. (2019a) study, Huang et al. (2020) found a correlation between the halo mass and the ratio between the stellar mass enclosed within two different apertures. They interpreted their measurement as a positive trend between galaxy size and halo mass, qualitatively in agreement with the Charlton et al. (2017) work and in disagreement with the Sonnenfeld et al. (2019a) study. The origin of this discrepancy is unclear: more robust (model-independent) measurements are needed to settle this issue.

Finally, Taylor et al. (2020) used weak lensing data to measure the average dark matter halo mass of galaxies in a narrow stellar mass range around log M_*/M_⊙ ∼ 10.5 as a function of a set of properties, including the stellar density profile and the star formation rate. They found that the halo mass correlates most strongly with the Sérsic index (positive correlation) and with the half-light radius (negative correlation). Their sample, however, included both star-forming galaxies and quiescent ones, so it cannot be compared to the studies cited above directly.

All of these measurements can, in principle, be used to test predictions from theoretical models and hydrodynamical simulations (Shankar et al. 2014; Correa & Schaye 2020; Cui et al. 2021). Moreover, detecting correlations between the inferred halo mass and galaxy properties at a fixed stellar mass would provide further evidence for the existence of dark matter: if stars were the only mass component, it would be very difficult to explain a variation in the strength of the gravitational lensing signal among galaxies of the same baryonic mass (Brouwer et al. 2021). In practice, however, systematic uncertainties in the observations limit our ability to make precise quantitative statements.

A very important source of systematic uncertainty is the measurement of the stellar masses of galaxies. The standard way to obtain stellar masses is by fitting stellar population synthesis models to the observed photometric or spectroscopic data. Choices in various aspects of these models, such as the duration of each stellar evolutionary stage, the stellar initial mass function (IMF), the dust attenuation, the metallicity and its evolution (or lack thereof), and the star formation history, can introduce biases as large as 0.3 dex (Conroy 2013).

This systematic uncertainty limits not only the accuracy of our knowledge of the stellar-to-halo mass ratio, but it can also introduce spurious correlations with secondary parameters or, vice versa, erase the signature of correlations. One illustrative example is the study by Correa & Schaye (2020) of the stellar-to-halo mass relation of Sloan Digital Sky Survey (SDSS; York et al. 2000) galaxies as a function of galaxy morphology. Using group richness as a proxy for halo mass and taking stellar masses from the catalogue of Kauffmann et al. (2003), they found that, at a fixed halo mass, disk galaxies have a higher stellar mass compared to ellipticals. However, the signal disappeared when using the stellar mass catalogue of Chang et al. (2015). Since there is no strong reason to favour one catalogue over the other, it is difficult to draw firm conclusions on the halo mass-galaxy morphology relation with those data.

The goal of this work is to measure the distribution in the dark matter halo mass of central massive elliptical galaxies as a function of observationally robust properties, using weak lensing data from the Kilo Degree Survey (KiDS; de Jong et al. 2013). We choose to work with massive galaxies because they are the objects for which the lensing signal is strongest, and we limit the analysis to ellipticals because the regularity of their surface brightness profile allows for more accurate measurements of their structural properties, compared to spiral galaxies. Finally, we use central galaxies in order to have a cleaner interpretation of the lensing signal. This also implies that our measurements are mostly probing the effects of physical processes internal to the halos, as opposed to the environmental processes to which satellite galaxies are subject (such as stripping or galaxy harassment).

We first focus on the distribution of halo mass as a function of the following quantities: redshift, luminosity, colour, and surface brightness profile. While the use of luminosity instead of stellar mass can be seen as a step back compared to recent studies on the relation between galaxies and halos, galaxy luminosities do not suffer from the numerous sources of systematic uncertainties mentioned above and can be determined much more accurately than stellar masses. On the one hand, this approach shifts the burden of predicting the stellar population properties to the theory side, when comparing our measurements to theoretical models. On the other hand, it makes these comparisons more self-consistent because many of the currently available models of galaxy evolution already make predictions on key stellar population properties, such as the star formation history (and hydrodynamical simulations can predict full spectral energy distributions, see for example Trayford et al. 2017).

In a second stage, we tackle an additional source of systematic errors: the parameterisation of the stellar surface brightness profile. The total luminosity of a galaxy at a cosmological distance is not directly observable because its surface brightness usually extends out to large radii (Tal & van Dokkum 2011) and falls below the detection limit of current photometric survey data. The total flux is then usually obtained by fitting a model surface brightness profile to the observed data and extrapolating the model to infinity or to a very large radius. This procedure can introduce biases on the derived luminosity and half-light radius, in case the true surface brightness deviates from the model one. Sonnenfeld (2020) found that, for a typical massive elliptical galaxy in the SDSS spectroscopic sample imaged at a depth of ∼26 mag, 20% or more of the total flux obtained by fitting a Sérsic profile is due to extrapolation.

We eliminate this source of systematic errors by employing a more robust description of the surface brightness profile, introduced by Sonnenfeld (2020): instead of the total luminosity and half-light radius, we use the flux within a circularised aperture of radius 10 kpc and the light-weighted projected stellar density slope in the same aperture. We then measure the distribution in halo mass as a function of these two quantities, as well as redshift and colour. These measurements will provide a robust benchmark against which predictions from theoretical models can be tested.

The structure of this paper is as follows. In Sect. 2 we describe the data used for our study, including sample selection criteria. In Sect. 3 we determine the distribution in luminosity, half-light radius, and colour of our sample. In Sect. 4 we measure the distribution in dark matter halo mass as a function of redshift, luminosity, half-light radius, and colour. In Sect. 5 we repeat the analysis by replacing the total luminosity and half-light radius with the aperture luminosity within 10 kpc and the light-weighted surface brightness slope. We discuss our results in Sect. 6 and provide our conclusions in Sect. 7.

We assume a flat Λ cold dark matter cosmology with Ω_M = 0.3 and H₀ = 70 km s⁻¹ Mpc⁻¹. Magnitudes are in AB units, and masses and luminosities are in solar units. Quoted uncertainties correspond to the 68% credible region.

2. Data

2.1. Sample selection

Our goal is to select a sample consisting of quiescent galaxies with an elliptical morphology and available spectroscopic data, with which to obtain precise redshift, luminosity, rest-frame colour, and size measurements. We also wish to carry out a weak lensing analysis of the same sample. For this purpose, we require our galaxies to lie in regions of the sky covered by the Data Release 4 (DR4) of KiDS (Kuijken et al. 2019). For the sake of a clean interpretation of our results, particularly from the weak lensing part, we wish to limit our study to central galaxies. We achieved these goals in a series of steps, as described below.

We selected galaxies from the spectroscopic sample of the SDSS Data Release 12 (DR12; Alam et al. 2015). In particular, we considered objects with spectra taken with the SDSS spectrograph (Smee et al. 2013) belonging to the Luminous Red Galaxy sample (Eisenstein et al. 2001) of the SDSS Legacy Survey. We first applied a minimum redshift condition, selecting only objects at z > 0.15. This cut is motivated by the fact that the galaxy stellar mass function of the SDSS Luminous Red Galaxy sample is more or less independent of redshift above z = 0.15, by construction, while it is a strong function of redshift below this limit. This makes it easier to separate trends with redshift from trends with luminosity.

To select quiescent galaxies, we applied a cut to the Hα equivalent width EW_Hα, selecting only objects with EW_Hα > −3 Å (the negative equivalent width corresponds to emission). This step, together with our minimum redshift condition, resulted in a sample of ∼3700 luminous red galaxies at z > 0.15 with KiDS DR4 photometry. We removed obvious disk-like objects by eliminating galaxies with a projected r-band minor-to-major axis ratio smaller than 0.5, as measured from KiDS photometric data (more details on these measurements are provided in Sect. 2.3): they account for 8% of the sample.

We then removed ∼200 objects that are identified as cluster members with a probability being a central galaxy smaller than 0.5 in the red-sequence Matched-filter Probabilistic Percolation (redMaPPER) catalogue (Rykoff et al. 2014). Using KiDS colour-composite images in gri bands, we visually inspected the remaining objects, and we eliminated galaxies with spiral arms, with disk or ring features, or those in very close proximity to objects with a comparable or higher brightness. This visual inspection step removed a further 15% of the galaxies in the sample. Finally, we applied a galaxy isolation cut, which was needed in order to render the weak lensing analysis computationally tractable. For each galaxy, we drew a cone around it with a 500 kpc physical radius at its redshift. From the sample, we then removed any galaxy where the cone overlapped with that of another more luminous galaxy. In Sect. 4 we give more details on the rationale for this choice. Our final sample consists of 2409 galaxies spanning the redshift range 0.15 < z < 0.4. The tenth percentile, median, and 90th percentile of the distribution in log M_* of the sample (measured as described in Sect. 2.4) are 11.47, 11.68, and 11.93, respectively.

2.2. KiDS photometry

Our measurements of the surface brightness profiles of our galaxies and the weak lensing data are based on photometric observations from KiDS DR4. KiDS is a photometric survey carried out from the ESO VLT Survey Telescope (VST; Capaccioli & Schipani 2011; Capaccioli et al. 2012), with the main goal of obtaining cosmological constraints from cosmic shear measurements. KiDS DR4 covers an area of approximately 1000 square degrees, imaged in u, g, r, and i bands with the instrument OmegaCAM (Kuijken 2011). In addition to optical data, KiDS DR4 includes matched-photometry measurements in the near-infrared bands Z, Y, J, H, and K_s obtained from the VISTA telescope as part of the VISTA Kilo-degree INfrared Galaxy survey (VIKING; Edge et al. 2013). The r-band data is the one with the best image quality, with a typical seeing FWHM < 0.8″, while the 5σ mean limiting magnitude within a 2″ aperture is ≈25. For this reason, we used the r-band to measure the structural parameters of our galaxies, as described in the next section. We refer readers to Kuijken et al. (2019) for details on the KiDS DR4 data reduction.

2.3. Surface brightness profile measurements

We describe the surface brightness distribution of each galaxy with an elliptical Sérsic profile,

$$\begin{array}{c}\hfill I(x,y)=I_{0}exp\{-b(n){\left(\frac{R}{R_{\mathrm{e}}}\right)}^{1/n}\},\end{array}$$(1)

where x and y are Cartesian coordinates with an origin at the galaxy centre and aligned with the galaxy major and minor axes, respectively, and R is the circularised radius

$$\begin{array}{c}\hfill R^2\equiv q{x}^2+\frac{y^2}{q},\end{array}$$(2)

where q is the axis ratio, n is the Sérsic index, and b(n) is a numerical constant defined in such a way that the isophote of radius R = R_e encloses half of the total light (see Ciotti & Bertin 1999). Throughout this paper, we refer to R_e as the half-light radius and as the size of a galaxy, interchangeably. The best-fit Sérsic parameters for each galaxy were obtained by fitting seeing-convolved Sérsic models to sky-subtracted and coadded r-band images using the software 2DPHOT (La Barbera et al. 2008), largely following the procedure adopted by Roy et al. (2018) for the analysis of data from previous KiDS data releases.

We obtained magnitudes in all nine bands by considering Galactic extinction-corrected colours, measured with the Gaussian Aperture and point spread function (PSF) technique (GAaP Kuijken et al. 2015), and using them to rescale the Sérsic fit-based r-band flux accordingly. An implicit assumption in this step is that colours do not vary across the spatial extent of a galaxy. We discuss the possible impact of this assumption in Sect. 6.2.

2.4. Rest-frame colour and luminosity measurements

We base our analysis on measurements of the rest-frame r-band luminosity and u − g colour. We obtained these quantities by fitting spectral energy distribution templates, obtained with the software GALAXEV (Bruzual & Charlot 2003), to the observed KiDS multi-band data. In particular, we generated composite stellar population synthesis models assuming a Chabrier stellar IMF (Chabrier 2003) and an exponentially declining star formation history, a range of values of metallicity, and dust attenuation. We then optimised for the parameters of the model given the observed 9-band magnitudes to obtain a best-fit spectral template, with which we derived the rest-frame r-band luminosity and u − g colour.

The model is not able to fit the observed magnitudes down to the noise level, which is very low (typically a few millimags in all bands, except u). The derived u − g colour, then, is sensitive to the spectral energy distribution of the galaxy over the full wavelength range probed by KiDS photometry, and it should therefore be considered as an average measure of galaxy redness over a broad spectral range. To make sure that the results are independent of our way of measuring u − g, we repeated the analysis using a more direct estimator of rest-frame colour, the D_n4000 spectral index, which we obtained from the SDSS spectroscopic catalogue. The D_n4000 index is defined as the ratio between the continuum level redward and blueward of the 4000 Å break (see Balogh et al. 1999, for details). Our measurements of the distribution of halo mass as a function of D_n4000 are indistinguishable from those obtained using u − g. Therefore, we only show the latter in this paper.

One by-product of this spectral energy distribution fitting process is an estimate of the stellar mass of each galaxy. We chose not to use this information for the core of our analysis since the stellar mass parameter is highly degenerate with the stellar IMF, the metallicity, and the star formation history parameters, and because the way that these aspects of our composite stellar population models are set is somewhat arbitrary. The only aspect of our study that relies on these stellar mass measurements is the estimate of the (small) baryonic contribution to the weak lensing signal, in Sect. 4.

2.5. Aperture luminosity and light-weighted slope

We also consider an alternative parameterisation of the luminosity profile of our galaxies, based on quantities that are measured more robustly than the total luminosity and the half-light radius. We use the rest-frame r-band luminosity enclosed within a circularised aperture of physical size 10 kpc, L_r, 10, and the light-weighted surface brightness slope within the same aperture, Γ₁₀. Following Sonnenfeld (2020), the latter is defined as

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{10}\equiv -{\displaystyle \frac{2\pi {\int }_0^{10}R{\displaystyle \frac{\mathrm{d}log{I}_r}{\mathrm{d}logR}}{I}_r(R)\mathrm{d}R}{2\pi {\int }_0^{10}R{I}_r(R)\mathrm{d}R}}=2-{\displaystyle \frac{2\pi {(10)}^2{I}_r(10)}{L_{r,10}}},\end{array}$$(3)

where I_r(R) is the circularised rest-frame r-band surface brightness profile of the galaxy.

To calculate L_r, 10 and Γ₁₀, we assume that the true surface brightness profile of each galaxy is given by its best-fit Sérsic model. This is a good approximation as the typical deviation in the surface brightness profile of massive elliptical galaxies from a Sérsic model, in the region constrained by the data, is a few percent (see Fig. 4 of Sonnenfeld 2020). Although based on the same Sérsic model used to estimate L_r and R_e, the advantage of using L_r, 10 and Γ₁₀ is that they do not rely on any extrapolation to large radii where the surface brightness profile is unconstrained, and they are therefore determined more robustly.

2.6. Background source shape and redshift measurements

Our weak lensing analysis is based on shape measurements of galaxies in the background of the objects in our sample, obtained by Giblin et al. (2021). Shapes are described by two ellipticity components, (ϵ₁, ϵ₂), which were measured from individual r-band exposures with the LENSFIT algorithm (Miller et al. 2007, 2013; Fenech Conti et al. 2017). LENSFIT also assigns a weight w_s to each source, which takes into account the uncertainty in the ellipticity measurements and the intrinsic scatter in the ellipticity distribution of the source sample. The latter dominates the error budget for most sources (see Miller et al. 2013, for details). As we explain in Sect. 4, shape measurements are weighted by their w_s in the weak lensing analysis. We estimated multiplicative shear biases using image simulations from Kannawadi et al. (2019). The uncertainties on these are of the order of 10⁻³.

Background sources were selected on the basis of their photometric redshifts (photo-zs). These were obtained using the Bayesian Photometric Redshift code (BPZ; Benítez 2000; Coe et al. 2006) on KiDS u, g, r, and i and VIKING Z, Y, J, H, and K_s PSF-matched photometry, using a source magnitude prior from Raichoor et al. (2014). We refer readers to Wright et al. (2019) for details on the photo-z measurements. For each lens at redshift z_l, we selected only galaxies with a maximum a posteriori redshift z_B larger than z_l + 0.2. Following Hildebrandt et al. (2021), we then determined the weighted redshift distribution of the resulting sample, n_w(z_s), incorporating the LENSFIT weight. This was done with a set of objects with spectroscopic redshift measurements which was matched to the source sample with a self-organising map approach (Wright et al. 2020). The redshift distribution of the spectroscopic sample was re-weighted to match the relative abundance of source galaxies with different photometric measurements, multiplied by their LENSFIT weight.

3. Distribution in luminosity, half-light radius, and colour

In Fig. 1 we show the distribution in redshift, luminosity, half-light radius, Sérsic index, rest-frame u − g colour, L_r, 10, and Γ₁₀ of our galaxy sample. We can clearly see that R_e correlates positively with L_r and n, while Γ₁₀ anti-correlates with L_r, 10 (as measured by Sonnenfeld 2020, on a similar sample of elliptical galaxies). The distribution in the Sérsic index, however, shows bimodality, with a secondary peak at around n = 10. It is not clear whether this bimodality is real, or if it is a spurious result of the fit with 2DPHOT. Visual inspection did not reveal any obvious difference between galaxies with n > 9 and the rest of the sample, and we found no correlation between n and the χ² of the photometric fit.

Fig. 1. Distribution in redshift, rest-frame r-band luminosity, half-light radius, Sérsic index, rest-frame u − g colour, luminosity within 10 kpc, and light-weighted surface brightness slope within 10 kpc of the galaxy sample. Contour levels enclose 68% and 95% of the sample, respectively.

Galaxies with a very large Sérsic index can be problematic because a large fraction of their model flux is due to an extrapolation of the surface brightness profile to regions that are below the detection limit of the survey (see Fig. 2 of Sonnenfeld 2020). However, as we explained in Sect. 1, the parameterisation in terms of L_r, 10 and Γ₁₀ is robust with respect to this extrapolation. For these reasons, we decided not to apply any further cut, and we carry out the analysis with all galaxies, regardless of their Sérsic index.

As a first step in our analysis, we are interested in compressing these data into a few numbers describing how the size and the colour of our galaxies correlate with each other, how they scale with redshift and luminosity, and what the scatter is around their mean relations. For this purpose, we fitted a bivariate Gaussian in log R_e and rest-frame u − g colour to the data. We chose this model because it is the simplest one allowing us to infer the above quantities. We quantify the goodness of the Gaussian assumption later in this section. Redshift and luminosity are treated as independent variables and their distribution is not modelled.

We write our model as follows (for the sake of keeping the notation compact, we indicate the rest-frame u − g colour as C):

$$\begin{array}{c}\hfill \mathrm{P}(log{R}_{\mathrm{e}},C)=\mathcal{R}(R_{\mathrm{e}}|z,L_r)\mathcal{K}(C|z,L_r,R_{\mathrm{e}}).\end{array}$$(4)

In the above equation, the factor ℛ describes the distribution in half-light radius as a function of redshift and luminosity. This is a Gaussian in log R_e

$$\begin{array}{c}\hfill \mathcal{R}(log{R}_{\mathrm{e}}|z,L_r)=\frac{1}{\sqrt{2\pi }{\sigma }_{\mathrm{R}}}exp\{-\frac{{(log{R}_{\mathrm{e}}-{\mu }_{\mathrm{R}}(z,L_r))}^2}{2{\sigma }_{\mathrm{R}}^2}\},\end{array}$$(5)

with mean

$$\begin{array}{c}\hfill {\mu }_{\mathrm{R}}(z,L_r)={\mu }_{\mathrm{R},0}+{\zeta }_{\mathrm{R}}(z-0.3)+{\beta }_{\mathrm{R}}(log{L}_r-11.3)\end{array}$$(6)

and dispersion σ_R. In Eq. (6), the parameter μ_R, 0 indicates the average value of log R_e at the pivot point z = 0.3 and log L_r = 11.3, while coefficients ζ_R and β_R describe a linear scaling of log R_e with z and log L_r, respectively.

The factor 𝒦 describes the distribution in colour as a function of redshift, luminosity, and half-light radius. This is also a Gaussian distribution

$$\begin{array}{c}\hfill \mathcal{K}(C|z,L_r,log{R}_{\mathrm{e}})=\frac{1}{\sqrt{2\pi }{\sigma }_{\mathrm{C}}}exp\{-\frac{{(C-{\mu }_{\mathrm{C}}(z,L_r,R_{\mathrm{e}}))}^2}{2{\sigma }_{\mathrm{C}}^2}\},\end{array}$$(7)

with the mean that scales with redshift, luminosity, and half-light radius as

$$\begin{array}{cc}\hfill {\mu }_{\mathrm{C}}(z,L_r,R_{\mathrm{e}})=& {\mu }_{\mathrm{C},0}+{\zeta }_{\mathrm{C}}(z-0.3)+{\beta }_{\mathrm{C}}(log{L}_r-11.3)\hfill \\ \hfill & +{\xi }_{\mathrm{C}}(log{R}_{\mathrm{e}}-{\mu }_{\mathrm{R}}(z,L_r))\hfill \end{array}$$(8)

and dispersion σ_C. We would like to point out that in Eq. (8), the pivot point in log R_e space is not fixed, but it is set to the average half-light radius μ_R(z, L_r), defined in Eq. (6), which scales with redshift and luminosity. This means that ξ_C indicates how colour scales with relative size, compared to the average value of log R_e of galaxies with the same redshift and luminosity.

We fitted the parameters of the model with a Markov chain Monte Carlo (MCMC), using the affine-invariant sampling method of Goodman & Weare (2010), implemented in Python with the package EMCEE (Foreman-Mackey et al. 2013). The statistical uncertainties on all of the observed quantities, redshift, luminosity, size, and colour are very small (the typical photometric error is 0.01 mag); therefore, we neglected them when carrying out the fit. After performing the fit, we removed the outliers, which we identified with a sigma-clipping procedure, and we repeated the fit on the new sample, iterating until convergence. To be precise, we labelled as outliers any objects with values of log R_e that are more than 4σ_R away from the average size μ_R(z, L_r) or with values of C that are more than 4σ_C away from the average colour μ_C(z, L_r, R_e). We removed 28 galaxies from the sample with this procedure. We visually inspected these galaxies and found that, for the most part, they are objects with a bright object in their proximity: we suspect that either the GAaP or Sérsic model fits, on which our photometry is based, suffered from a catastrophic failure due to blending.

The median value and 68% credible limits on the marginal posterior probability distribution of each model parameter are reported in Table 1. The results can be summarised as follows. Size scales with luminosity with a steeper-than-linear relation, β_R = 1.381 ± 0.018, and it decreases with increasing redshift at a fixed luminosity: ζ_R = −0.79 ± 0.05. This is a steeper evolution than the change in size at a fixed stellar mass, as reported in the literature: for example, Newman et al. (2012) observed an average size evolution rate at fixed M_* of d⟨log R_e⟩/dz ≈ −0.12 over the redshift range 0 < z < 1. The difference is due to the luminosity evolution of galaxies at a fixed stellar mass: at a fixed luminosity, passive galaxies at a higher redshift are intrinsically less massive than their lower redshift counterparts, and therefore their sizes are smaller on average. Our measurement of the intrinsic scatter in log R_e at a fixed redshift and luminosity is σ_R = 0.183 ± 0.003.

The rest-frame u − g colour decreases with an increasing redshift at a fixed luminosity, ζ_C = −0.412 ± 0.015, which is qualitatively in agreement with a passive evolution scenario, and almost independent of luminosity at a fixed redshift, β_C = −0.009 ± 0.005. Moreover, at a fixed redshift and luminosity, the colour is independent of excess size: ξ_C = −0.007 ± 0.006.

This finding is somewhat surprising, given previous similar studies from the literature. For instance, Kim et al. (2018) found a positive, albeit weak, correlation between the stellar surface mass density of a sample of spheroidal galaxies and the D_n4000 spectral index. Zahid & Geller (2017) also report a correlation between D_n4000 and the size at a fixed stellar mass in a sample of quiescent galaxies. In order to compare our finding with these studies, we repeated our analysis by using D_n4000 measurements in place of the colour. We found that D_n4000 has a weak positive correlation with luminosity, thus still in apparent discrepancy with the studies mentioned above. As we show in Sect. 6.1, however, using stellar mass in place of luminosity leads to an anti-correlation between size and colour. Moreover, the apparent discrepancy with respect to the Zahid & Geller (2017) and Kim et al. (2018) studies could be due to differences in the sample selection criteria: for example, the Zahid & Geller (2017) sample is dominated by galaxies with a stellar mass below 10¹¹ M_⊙, which are excluded by our luminosity cut.

The measurements shown in this section were obtained under the assumption of a bi-variate Gaussian distribution in log R_e and C at a fixed redshift and luminosity. In order to verify the accuracy of this model, in Fig. 2 we show the observed size and colour distributions as a function of the predicted average value of these quantities, obtained via Eqs. (6) and (8). The model is able to correctly reproduce the average size and colour of the majority of the sample. However, there is a subset of galaxies at low luminosity, redder colours, and low redshift that appear to deviate from the 1 : 1 line. We repeated the analysis by removing objects with z < 0.2. While this additional selection caused the inferred distribution in size and colour to change, the final measurements of the halo mass distribution did not vary. Therefore, we present the results obtained with the full sample.

Fig. 2. Comparison between observed and model-predicted sizes and colours. Left: observed half-light radius as a function of the average Re of galaxies with the same redshift and luminosity, as predicted by Eq. (6). Right: observed rest-frame u − g colour as a function of the average colour of galaxies with the same redshift, luminosity, and size, as predicted by Eq. (8). Uncertainties on the model parameters are neglected, as they are very small with respect to the observational scatter. The dashed line marks the 1 : 1 relation.

To check the robustness of the Gaussian scatter assumption, we measured the fraction of galaxies with values of R_e that are within 1σ_R of μ_R(z, L_r). This is 71%, which is slightly larger than the expected value of 68% of a pure Gaussian. Similarly, the fraction of galaxies with colour within 1σ_C of μ_C(z, L_r, R_e) is 71%. These small deviations from a Gaussian distribution are not concerning for the goals of this work.

4. Dark matter halo mass distribution

4.1. Model

We wanted to measure the distribution in dark matter halo mass as a function of redshift, luminosity, size, and colour. We did so with a Bayesian hierarchical inference approach, following the method of Sonnenfeld & Leauthaud (2018). In order to fit the weak lensing data, we needed to model the projected surface mass density distribution of the lens sample. We only used data covering relatively small scales, 50 − 500 kpc, where the lensing signal from the lenses themselves dominates over the contribution from the large-scale structure. Therefore, we modelled each lens simply as the sum of a dark matter halo and a stellar component. We describe the former with a spherical Navarro, Frenk, & White profile (NFW, Navarro et al. 1997), and the latter with a circularised Sérsic profile. We ignored any contribution from gas. We discuss the impact of these assumptions in Sect. 6.2.

The 3D density distribution of a spherical NFW profile, as a function of radius r, is the following:

$$\begin{array}{c}\hfill {\rho }_{\mathrm{DM}}(r)\propto \frac{1}{r{(1+r/r_{\mathrm{s}})}^2}·\end{array}$$(9)

It has two degrees of freedom, related to the normalisation and the scale radius r_s. We describe the normalisation in terms of the mass M_h enclosed within the spherical shell of radius r_h, within which the average density is 200 times the critical density of the Universe at the redshift of the galaxy, ρ_c. We refer to M_h simply as the halo mass. As a second parameter, we used the concentration c_h, defined as the ratio between r_h and the scale radius.

We then modelled the distribution in M_h and c_h of our lens population with the following functional form:

$$\begin{array}{cc}\hfill P(log{M}_{\mathrm{h}},log{c}_{\mathrm{h}}|z,L_r,R_{\mathrm{e}},C)=& \mathcal{H}(log{M}_{\mathrm{h}}|z,L_r,R_{\mathrm{e}},C)\hfill \\ \hfill & \times \mathcal{C}(log{c}_{\mathrm{h}}|M_{\mathrm{h}}).\hfill \end{array}$$(10)

We modelled the first term as a Gaussian in log M_h,

$$\begin{array}{c}\hfill \mathcal{H}(log{M}_{\mathrm{h}}|z,L_r,R_{\mathrm{e}},C)=\frac{1}{\sqrt{2\pi }{\sigma }_{\mathrm{h}}}exp\{-\frac{{(log{M}_{\mathrm{h}}-{\mu }_{\mathrm{h}})}^2}{2{\sigma }_{\mathrm{h}}^2}\},\end{array}$$(11)

with mean

$$\begin{array}{cc}\hfill {\mu }_{\mathrm{h}}(z,L_r,R_{\mathrm{e}},C)=& {\mu }_{\mathrm{h},0}+{\zeta }_{\mathrm{h}}(z-0.3)+{\beta }_{\mathrm{h}}(log{L}_r-11.3)\hfill \\ \hfill & +{\xi }_{\mathrm{h}}(log{R}_{\mathrm{e}}-{\mu }_{\mathrm{R}}(z,L_r))\hfill \\ \hfill & +{\nu }_{\mathrm{h}}(C-{\mu }_{\mathrm{C}}(z,L_r,R_{\mathrm{e}}))\hfill \end{array}$$(12)

and dispersion σ_h. Similarly to our parameterisation of the colour distribution in Sect. 3, in Eq. (12) we allowed for the average log M_h to scale with redshift, luminosity, and excess size, via the parameters ζ_h, β_h, and ξ_h. In addition, we allowed for a correlation between the halo mass and excess colour C − μ_C(z, L_r, R_e), which is the difference between the colour of a galaxy and the average colour of galaxies of the same redshift, luminosity, and size, via parameter ν_h. Sonnenfeld & Leauthaud (2018) show how a log-Gaussian distribution in M_h of the same kind as that found in Eq. (11) is able to accurately capture the average halo mass of a sample of massive galaxies and the correlations between M_h and galaxy properties, when fitted to weak lensing data similar to those of the KiDS survey.

Then, we described the term 𝒞 in Eq. (10) with a fixed mass-concentration relation for the distribution in c_h:

$$\begin{array}{c}\hfill \mathcal{C}(c_{\mathrm{h}}|M_{\mathrm{h}})={\displaystyle \frac{1}{\sqrt{2\pi }{s}_{\mathrm{c}}}}exp\{-{\displaystyle \frac{{(log{c}_{\mathrm{h}}-{\mu }_{\mathrm{c}}(M_{\mathrm{h}}))}^2}{2s_{\mathrm{c}}^2}}\},\end{array}$$(13)

with mean

$$\begin{array}{c}\hfill {\mu }_{\mathrm{c}}(M_{\mathrm{h}})={\mu }_{\mathrm{c},0}+{\beta }_{\mathrm{c}}(log{M}_{\mathrm{h}}-12)\end{array}$$(14)

and dispersion s_c. Following the N-body simulation predictions of Macciò et al. (2008) based on WMAP5 cosmology, we fixed μ_c, 0 = 0.83, β_c = −0.097, and s_c = 0.1. We discuss the impact of this choice in Sect. 6.2.

Our model distribution in dark matter halo mass is then summarised by the following set of free parameters:

$$\begin{array}{c}\hfill \mathit{\eta }\equiv \{{\mu }_{\mathrm{h},0},{\zeta }_{\mathrm{h}},{\beta }_{\mathrm{h}},{\xi }_{\mathrm{h}},{\nu }_{\mathrm{h}},{\sigma }_{\mathrm{h}}\}.\end{array}$$(15)

In principle, μ_R and μ_C are also free parameters that enter the halo mass distribution through Eq. (12). In practice, however, we fixed their values to the median of their marginal posterior distribution obtained from the analysis of Sect. 3, for the purpose of reducing the dimensionality of the problem. This is a fair approximation since the uncertainties on μ_R and μ_C are very small (see Table 1).

We wanted to determine the posterior probability distribution of the model parameters η given the data d. Using Bayes’ theorem, this is given by

$$\begin{array}{c}\hfill P(\mathit{\eta }|\mathit{d})\propto P(\mathit{\eta })P(\mathit{d}|\mathit{\eta }),\end{array}$$(16)

where P(η) is the prior probability on the parameters and P(d|η) is the likelihood of observing the data given the model.

Weak lensing data consist of an ensemble of shape measurements on galaxies located behind the lens. The main observable is the complex ellipticity, ϵ = (ϵ₁, ϵ₂). The expectation value of the tangential component (with respect to the lens galaxy) of the ellipticity ϵ_t is the tangential reduced shear at the image position θ and source redshift z_s, g_t(θ, z_s), corrected for the shear multiplicative bias (1 + m):

$$\begin{array}{c}\hfill E\left[{ϵ}_{\mathrm{t}}\right]=(1+m)g_{\mathrm{t}}(\mathit{\theta },z_{\mathrm{s}}).\end{array}$$(17)

To predict the value of the reduced shear at each source position given our model, we worked under the isolated lens assumption of Sonnenfeld & Leauthaud (2018): we assumed that only the lens galaxy and relative dark matter halo contribute to the lensing signal on sources around it. This allows us to predict the value of the reduced shear g_t given the halo mass, the halo concentration, the redshift, and the stellar mass of a lens galaxy. We then approximated the likelihood of a single shape measurement ${ϵ}_{\text{t}}^{(\text{obs})}$ as a Gaussian centred on Eq. (17) and with a dispersion equal to the inverse square root of the LENSFIT weight, ${\mathit{w}}_{\mathrm{s}}^{-1/2}$. For a given source redshift z_s, the shape measurement likelihood is then

$$\begin{array}{c}\hfill P({ϵ}_{\mathrm{t}}^{(\mathrm{obs})}|M_{\mathrm{h}},c_{\mathrm{h}},z_{\mathrm{s}})=\mathcal{N}((1+m)g_{\mathrm{t}}(\mathit{\theta },z_{\mathrm{s}},M_{\mathrm{h}},c_{\mathrm{h}}),w_{\mathrm{s}}^{-1}),\end{array}$$(18)

where 𝒩(μ, σ²) indicates a Gaussian distribution with mean μ and variance σ². The exact individual redshifts of the source galaxies, however, are not known. We thus marginalise Eq. (18) over the LENSFIT weight-weighted source redshift distribution, n_w(z_s):

$$\begin{array}{c}\hfill P({ϵ}_{\mathrm{t}}^{(\mathrm{obs})}|M_{\mathrm{h}},c_{\mathrm{h}})={\displaystyle \int }\mathrm{d}{z}_{\mathrm{s}}{n}_w(z_{\mathrm{s}})\mathcal{N}((1+m)g_{\mathrm{t}}(\mathit{\theta },z_{\mathrm{s}},M_{\mathrm{h}},c_{\mathrm{h}}),w_{\mathrm{s}}^{-1}),\end{array}$$(19)

where the model reduced shear term in the integrand is zero for values of z_s that are smaller than the lens redshift.

The function n_w(z_s) was determined using a spectroscopic sample of galaxies matched to the entire KiDS DR4 shape measurements catalogue (Hildebrandt et al. 2021). This is not necessarily the same as the redshift distribution of the sources around our lenses as our source sample could be contaminated by galaxies physically associated with the local environment of the lenses. To correct for this difference, we first made five radial bins, then compared the number density of sources around our lens sample with that around a large set of random points covering the same footprint as the lenses. We found that the former is on average 10% larger than the latter, with a mild dependence on the radius. Using these measurements, we corrected the n_w(z_s) of the sources in each radial bin by boosting the number density of sources (which are not lensed) at the same redshift of the lens.

The full likelihood for an ensemble of shape measurements around a single lens is the product over the sources of terms of the form of Eq. (18). We only considered sources in an annulus of physical radius at the lens redshift between 50 and 500 kpc. The isolation criterion applied in the construction of the sample, which is described in Sect. 2.1, ensures that we never use the same shape measurement more than once in different lenses. Additionally, it allows us to write the likelihood as the product of N independent factors, with N being the number of lenses in our sample:

$$\begin{array}{c}\hfill P(\mathit{d}|\mathit{\eta })={\displaystyle \prod _i^{N}P({\mathit{d}}_i|\mathit{\eta }),}\end{array}$$(20)

where d_i are the shape measurements around the ith lens. As explained above, the shear data depend on the values of the halo mass and concentration. To evaluate each term in the product of Eq. (20), it is then necessary to average over all possible values taken by these quantities, as predicted by the following model:

$$\begin{array}{c}\hfill P({\mathit{d}}_i|\mathit{\eta })={\displaystyle \int }\mathrm{d}log{M}_{\mathrm{h}}\mathrm{d}log{c}_{\mathrm{h}}P({\mathit{d}}_i|M_{\mathrm{h}},c_{\mathrm{h}})P(M_{\mathrm{h}},c_{\mathrm{h}}|\mathit{\eta }),\end{array}$$(21)

where P(log M_h, log c_h|η) is the model probability distribution of Eq. (10).

We assumed flat priors on each parameter, with bounds reported in Table 2. We sampled the posterior probability distribution P(η|d) with an MCMC. At each draw of the set of parameters η, we evaluated the integrals of Eq. (21) with Monte Carlo integration and importance sampling, following Sonnenfeld et al. (2019b, see their Appendix A).

4.2. Results

The posterior probability distribution is shown in Fig. 3. The median and 68% credible region of the marginal posterior probability of each parameter is reported in Table 2.

Fig. 3. Posterior probability distribution of the model parameters describing the distribution in dark matter halo mass of the sample of lenses defined in Eq. (10). Contour levels correspond to a 68% and 95% enclosed probability.

We detect a positive correlation between halo mass and luminosity, with logarithmic slope β_h = 0.99 ± 0.29. We do not detect correlations between halo mass and either redshift, size, or colour: parameters ζ_h, ξ_h, and ν_h are all consistent with 0 within 2σ. We interpret these results as follows.

The correlation of halo mass with luminosity is closely related to that with stellar mass: this is confirmed by the fact that the value of β_h that we measured is consistent with the slope of the halo mass-stellar mass correlation measured by Sonnenfeld et al. (2018) on a similarly selected sample of galaxies. The lack of a strong correlation between halo mass and redshift at a fixed luminosity is the result of a combination of a few effects. First of all, there is an intrinsic evolution in the mapping between galaxy and halo masses, as both components grow with time; halo mass grows both via the accretion of matter and as a result of the evolution in the critical density of the Universe, which is an effect referred to as pseudo-evolution. Second, the luminosity evolution of passive galaxies implies that, at fixed luminosity, higher redshift galaxies are less massive than lower redshift ones, and therefore they live in less massive halos. Third, systematic errors in the photometric measurements of the background sources can introduce biases on the redshift evolution parameter ζ_h. For example, if the source redshifts are systematically overestimated as a result of contamination from galaxies physically associated with the lens groups – we only applied a global correction for this effect, which might not be entirely accurate – this biases the inference of the halo mass of galaxies at a different redshift differently due to the non-linear dependence of the lensing critical surface mass density on the lens and source redshift. Our main motivation for introducing a scaling of halo mass with redshift is precisely to separate any evolution signal, whether real or artificial, from the dependence of halo mass on luminosity, size, and colour. Therefore, we do not investigate the nature of the ζ_h measurement further.

The amplitude of the halo mass-size correlation, ξ_h = 0.20 ± 0.29, is consistent with the value found by Sonnenfeld et al. (2019a) on a similarly selected sample of galaxies, ξ_h = −0.14 ± 0.17, and corroborates the evidence against a strong correlation between these two quantities. Besides the data used, the main difference between the two measurements is that, in the present paper, ξ_h describes a correlation at fixed redshift and luminosity; whereas in the case of Sonnenfeld et al. (2019a), ξ_h refers to a correlation at a fixed stellar mass. Our measurement of ξ_h implies that the difference in the average halo mass between galaxies at the 84th and 16th percentile of the size distribution is

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{2{\sigma }_{\mathrm{R}}}⟨log{M}_{\mathrm{h}}⟩\equiv 2{\sigma }_{\mathrm{R}}{\xi }_{\mathrm{h}}=0.07\pm 0.10.\end{array}$$(22)

Similarly, given our constraint on the parameter describing the correlation between halo mass and colour, ν_h, we can quantify the difference in average halo mass between galaxies at the 84th and 16th percentile of the colour distribution at a fixed redshift, luminosity, and size, as

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{2{\sigma }_{\mathrm{C}}}⟨log{M}_{\mathrm{h}}⟩\equiv 2{\sigma }_{\mathrm{C}}{\nu }_{\mathrm{h}}=-0.03\pm 0.11.\end{array}$$(23)

The lack of a strong correlation between halo mass and rest-frame u − g colour at a fixed redshift appears to contrast with the positive trend found by Zu & Mandelbaum (2018) in their weak lensing and clustering study of red SDSS galaxies. We discuss the significance of this result in Sect. 6.3.

Part of the signal captured by ν_h could be due to a variation in the stellar mass-to-light ratio of galaxies with colour: if redder galaxies are intrinsically more massive at a fixed luminosity, their halo mass is also higher in virtue of the steep stellar mass-halo mass correlation. We discuss to what extent the coefficient ν_h is sensitive to this effect in Sect. 6.1.

Finally, we infer a value of the intrinsic scatter in log M_h of σ_h = 0.09 ± 0.08. While the inference is consistent with zero, our measurement does not allow us to put a strong constraint on this quantity: the 95%-ile of the marginal posterior distribution in σ_h is as large as 0.25.

4.3. Posterior predictive check

As a consistency test, we performed a posterior predictive check, which is how the goodness-of-fit is assessed with Bayesian hierarchical models. Given our sample of lens galaxies, we drew realisations of the distribution of halo mass from the posterior probability distribution P(η|d). For each realisation, we predicted the weak lensing signal, which we then compared with the observed data. In order to easily visualise the data, we considered the stacked weak lensing signal in bins for our comparison. We binned our lens sample in three different ways. First, we split it according to the galaxy excess size, dividing the sample into galaxies that lie above or below the average size given their redshift and luminosity, μ_R(z, L_r). Second, we split the sample according to the colour between galaxies that are redder or bluer than the average colour given their properties, μ_C(z, L_r, R_e). Third, we split the sample on the basis of the predicted average halo mass. In particular, we considered the set of model parameters η^(MAP) that maximises the posterior probability distribution and used it to define the maximum a posteriori average halo mass, ${\mu }_{\mathrm{h}}^{(\mathrm{MAP})}$. This is given by

$$\begin{array}{cc}\hfill {\mu }_{\mathrm{h}}^{(\mathrm{MAP})}=& 13.31+1.35(z-0.3)+1.02(log{L}_r-11.3)\hfill \\ \hfill & +0.19(log{R}_{\mathrm{e}}-{\mu }_{\mathrm{R}}(z,L_r))\hfill \\ \hfill & -0.29(C-{\mu }_{\mathrm{C}}(z,L_r,R_{\mathrm{e}})).\hfill \end{array}$$(24)

We then split the sample between galaxies with ${\mu }_{\mathrm{h}}^{(\mathrm{MAP})}$ larger and smaller than the median value, which is 13.35.

We computed the observed signal as follows. We made five radial bins in the range 50 kpc < R < 500 kpc, then we converted the measurement of the tangential ellipticity of each source into an estimate of the excess surface mass density,

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathrm{\Sigma }}^{(\mathrm{obs})}\equiv {ϵ}_{\mathrm{t}}{\mathrm{\Sigma }}_{\mathrm{cr}}\end{array}$$(25)

and related uncertainty, where Σ_cr is the critical surface mass density for gravitational lensing. Since we do not have exact source redshift measurements, we approximated Σ_cr with an average over the source redshift distribution n_z(z_s):

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{cr}}^{-1}\approx \frac{4\pi G}{c^2}D(z_{\mathrm{l}}){\displaystyle \int }\mathrm{d}{z}_{\mathrm{s}}{n}_w(z_{\mathrm{s}})\frac{D(z_{\mathrm{s}},z_{\mathrm{l}})}{D(z_{\mathrm{s}})},\end{array}$$(26)

where D(z) is the angular diameter distance as a function of redshift.

The definition of Eq. (25) is motivated by the fact that, in the case of a single lens, if Σ(R) is the projected surface mass density of the lens and Σ(R)≪Σ_cr, then the expectation value for the right-hand side quantity is the excess surface mass density ΔΣ(R) at the source position. The latter is defined as the difference between the average surface mass density within radius R and Σ(R):

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathrm{\Sigma }(R)=\overline{\mathrm{\Sigma }}(<R)-\mathrm{\Sigma }(R).\end{array}$$(27)

For each radial bin, we obtained a weighted average of ΔΣ^(obs), with an inverse variance weighting scheme, and related uncertainty. The result is shown as error bars in Fig. 4.

Fig. 4. Stacked excess surface mass density profile of the lenses in our sample. The sample is split into two bins according to size (left panel), colour (middle panel), and predicted average halo mass, depending on the maximum a posteriori value of μh defined in Eq. (24) (right panel). The error bars show the excess surface mass density obtained by stacking the observed tangential ellipticity in five radial bins. The bands show the same quantity obtained by posterior prediction from our model, without adding shape noise. For each sub-sample, the figure reports the probability of obtaining a scatter around the posterior predicted ΔΣ profile larger than the observed value, under the assumption that the model describes the truth. Small values for this quantity may indicate that the model is unable to reproduce these observations. Error bars have been offset slightly in the horizontal direction, for clarity.

We subsequently obtained the posterior-predicted expectation value of the stacked excess surface density profile, ΔΣ^(pp), as follows. At each draw of the values of the model parameters from the posterior probability distribution, we randomly drew a value for the halo mass and concentration for each lens in the sample. Then we computed the expectation value of the tangential ellipticity of each source, given its position, redshift, and the mass density profile of its lens. Finally, we obtained the inverse-variance weighted ΔΣ, as was done with the observed data. The 16th and 84th percentile of the posterior-predicted ΔΣ^(pp) profile of the two lens sub-samples for each binning scheme are shown as bands in Fig. 4.

The ΔΣ^(pp) profile reproduces the average observed signal for each sub-sample. In order to quantify the goodness of fit, we considered the standard deviation in the observed ΔΣ around the posterior predicted profile as test statistics: assuming our model to be accurate, we estimated the probability of obtaining a higher standard deviation than the observed value. Very low values for this probability, 5% or lower, might indicate the inability of the model to reproduce these data.

For each posterior draw, we generated mock ΔΣ observations in the five radial bins by adding noise to the posterior predicted profile; computed the standard deviation of the mock data points around ΔΣ^(pp), rms^(mock); and compared this quantity with the standard deviation of the observed data points around ΔΣ^(pp), rms^(obs). In Fig. 4 we report the probability that rms^(mock) > rms^(obs) for each sub-sample. The smallest value of this probability is 8%, which occurs for the sub-sample of higher ${\mu }_{\mathrm{h}}^{(\mathrm{MAP})}$. This means that there is no strong reason to question the ability of our model to reproduce the stacked weak lensing data.

5. Distribution as a function of robust observables

In this section we fit for the distribution of halo mass of the sample as a function of redshift; rest-frame r-band luminosity within an aperture of 10 kpc, L_r, 10; light-weighted projected slope within 10 kpc (defined in Eq. (3)), Γ₁₀; and the rest-frame u − g colour, C. Compared to the analysis of Sect. 4, we replace the total luminosity and the half-light radius with L_r, 10 and Γ₁₀. The motivation for this is to obtain a description of the halo mass distribution that is robust with respect to assumptions on the asymptotic behaviour of the surface brightness profile of the lens galaxies.

Similarly to our approach in Sect. 3, we first fit for the distribution in colour as a function of redshift, L_r, 10 and Γ₁₀. We do this so that we can later model the distribution in halo mass as a function of excess colour, that is the difference between the colour of a galaxy and the average colour of galaxies with the same redshift, L_r, 10 and Γ₁₀. We model the colour distribution with the same functional form as Eq. (7), but now expressed as a function of the new set of parameters:

$$\begin{array}{c}\hfill \mathcal{K}(C|z,L_{r,10},{\mathrm{\Gamma }}_{10})=\frac{1}{\sqrt{2\pi }{\sigma }_{\mathrm{C}}}exp\{-\frac{{(C-{\mu }_{\mathrm{C}}(z,L_{r,10},{\mathrm{\Gamma }}_{10}))}^2}{2{\sigma }_{\mathrm{C}}^2}\}·\end{array}$$(28)

The mean colour μ_C is allowed to scale with z, L_r, 10, and Γ₁₀ as follows:

$$\begin{array}{cc}\hfill {\mu }_{\mathrm{C}}(z,L_{r,10},{\mathrm{\Gamma }}_{10})=& {\mu }_{\mathrm{C},0}+{\zeta }_{\mathrm{C}}(z-0.3)+{\lambda }_{\mathrm{C}}(log{L}_r-10.9)\hfill \\ \hfill & +{\gamma }_{\mathrm{C}}({\mathrm{\Gamma }}_{10}-1.5).\hfill \end{array}$$(29)

We fitted for the parameters of the distribution of Eq. (28): the results are reported in Table 3. Galaxy colour does not depend on L_r, 10 or Γ₁₀ at a fixed redshift. This result matches the finding of Sect. 3: the colour of these galaxies is independent of their surface brightness profile. At a fixed surface brightness profile, colour decreases with increasing redshift with the same trend previously measured in Sect. 3: ζ_C = −0.418 ± 0.017.

We then model the distribution in halo mass and concentration of the sample. We use the same functional form as Eq. (10), expressed in terms of the new set of variables, for this purpose:

$$\begin{array}{cc}\hfill P(log{M}_{\mathrm{h}},log{c}_{\mathrm{h}}|z,L_{r,10},{\mathrm{\Gamma }}_{10},C)=& \mathcal{H}(log{M}_{\mathrm{h}}|z,L_{r,10},{\mathrm{\Gamma }}_{10},C)\hfill \\ \hfill & \times \mathcal{C}(log{c}_{\mathrm{h}}|M_{\mathrm{h}}).\hfill \end{array}$$(30)

Similarly to Eq. (11), the term ℋ is a Gaussian in log M_h, with a mean μ_h that is given by

$$\begin{array}{cc}\hfill {\mu }_{\mathrm{h}}(z,L_{r,10},{\mathrm{\Gamma }}_{10},C)=& {\mu }_{\mathrm{h},0}+{\zeta }_{\mathrm{h}}(z-0.3)+{\lambda }_{\mathrm{h}}(log{L}_r-10.9)\hfill \\ \hfill & +{\gamma }_{\mathrm{h}}({\mathrm{\Gamma }}_{10}-1.5)\hfill \\ \hfill & +{\nu }_{\mathrm{h}}(C-{\mu }_{\mathrm{C}}(z,L_{r,10},{\mathrm{\Gamma }}_{10})),\hfill \end{array}$$(31)

and dispersion σ_h. The quantity μ_C(z, L_r, 10, Γ₁₀) is the average colour given by Eq. (29), with the coefficients ζ_C, λ_C, and γ_C fixed to the median marginal posterior values reported in Table 3. The concentration distribution 𝒞 is the same as Eq. (13).

In Fig. 5 we show the posterior probability distribution of the model parameters (μ_h, 0, ζ_h, λ_h, γ_h, ν_h, σ_h). The median and 68% credible region of the marginal posterior of each parameter is given in Table 4. Our measurement of λ_h does not let us determine whether halo mass increases or decreases with L_r, 10 at a fixed redshift. However, we find that halo mass increases with decreasing Γ₁₀. The values of the parameters describing the redshift and colour dependence of the halo mass at a fixed surface brightness profile, ζ_C and ν_C, are slightly different compared to the inference obtained in Sect. 4, but consistent within the uncertainty. The difference in the average halo mass between galaxies at the 84%-ile and 16%-ile of the colour distribution that we obtain in this case is

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{2{\sigma }_{\mathrm{C}}}⟨log{M}_{\mathrm{h}}⟩\equiv 2{\sigma }_{\mathrm{C}}{\nu }_{\mathrm{h}}=0.00\pm 0.11.\end{array}$$(32)

Fig. 5. Posterior probability distribution of the model parameters describing the distribution in halo mass as a function of redshift, Lr, 10 and Γ10, Eq. (30). Contour levels correspond to a 68% and 95% enclosed probability.

This indicates that our constraint on the correlation between halo mass and colour is somewhat model-dependent, but, given the current uncertainty level, it is still robust with respect to the two parameterisations of the surface brightness profile that we employed. Finally we point out that, although the value of the parameter μ_h, 0 measured here differs from the value found in Sect. 4, this is simply the consequence of adopting a different pivot point for the description of the mean halo mass of Eq. (12), compared to Eq. (31).

6. Discussion

6.1. Expectations from variations in M_*/L

The analysis of Sects. 4 and 5 shows that, at a fixed redshift and surface brightness distribution, there is not a large difference in the average halo mass of redder and bluer elliptical galaxies. If the stellar mass-to-light ratio was identical for all of our galaxies, expressing trends in terms of luminosity or stellar mass would be equivalent. However, colour differences reflect differences in the properties of the stellar population, and therefore we can expect the stellar mass-to-light ratio to vary across our sample. In this section we investigate to what extent our result can change when using stellar mass in place of luminosity.

We proceed as follows. We use composite stellar population synthesis models obtained with GALAXEV (Bruzual & Charlot 2003) to assign stellar masses to galaxies. Rather than fitting for the stellar population parameters to the full set of photometric data, however, we assume a simplified model that results in a one-to-one mapping between the rest-frame u − g colour and stellar mass-to-light ratio. We do this to maximise the correlation between these two properties, so as to obtain an upper limit on the effect of variations in mass-to-light ratio with colour. Then, we repeat the analysis carried out in Sect. 4, using these assigned stellar masses in place of the luminosities.

We assume a model with a fixed metallicity of Z = 0.05 and an exponentially decaying star formation history with decay time τ = 1 Gyr. The time since the initial burst is the only parameter that we leave free to vary. Figure 6 shows how the stellar mass-to-light ratio varies as a function of colour for this model (green curve). Given the value of u − g of each galaxy, we use this curve to assign a value of M_*/L_r and, consequently, a stellar mass to it. For reference, in Fig. 6 we also show curves for two different values of the metallicity covered by the Bruzual & Charlot (2003) stellar population library. Although the value of the metallicity that we adopted is relatively high, we made this choice to ensure that our model can span the whole range of the observed u − g colour.

Fig. 6. Stellar mass-to-light ratio as a function of rest-frame u − g colour, for three composite stellar population models obtained with GALAXEV (Bruzual & Charlot 2003). Each model has a fixed metallicity and an exponentially declining star formation history with a decay time of 1 Gyr. Each curve was obtained by varying the time since the initial burst. The grey lines connect points of fixed time since the initial burst, as indicated by the labels. The model corresponding to a metallicity Z = 0.05 is used in Sect. 6.1 to assign stellar masses to galaxies, given their observed colour.

Following Sect. 3, we first determined the distribution in half-light radius as a function of the redshift and stellar mass of the sample. We found the following expression for the average log R_e:

$$\begin{array}{c}\hfill {\mu }_{\mathrm{R}}(z,M_{\ast })=1.08-0.21(z-0.3)+1.31(log{M}_{\ast }-11.6),\end{array}$$(33)

and an intrinsic scatter in log R_e at a fixed z and M_* of σ_R = 0.20. Then we measured the distribution in colour as a function of redshift, stellar mass, and excess size. This has a best-fit mean given by

$$\begin{array}{cc}\hfill {\mu }_{\mathrm{C}}(z,M_{\ast },R_{\mathrm{e}})=& 1.74-0.40(z-0.3)+0.05(log{M}_{\ast }-11.6)\hfill \\ \hfill & -0.10(log{R}_{\mathrm{e}}-{\mu }_{\mathrm{R}}(z,M_{\ast })),\hfill \end{array}$$(34)

and a scatter of σ_C = 0.05. We point out that, as we anticipated in Sect. 3, colour has a non-negligible anti-correlation with excess size, at fixed stellar mass and redshift.

We fitted for the distribution in halo mass as a function of redshift, stellar mass, excess size, and excess colour. We used the same method and the same model as in Sect. 4, with stellar mass in place of the luminosity. In particular, we assumed a Gaussian distribution in log M_h, with the following mean,

$$\begin{array}{cc}\hfill {\mu }_{\mathrm{h}}(z,M_{\ast },R_{\mathrm{e}},C)=& {\mu }_{\mathrm{h},0}+{\zeta }_{\mathrm{h}}(z-0.3)+{\beta }_{\mathrm{h}}(log{M}_{\ast }-11.6)\hfill \\ \hfill & +{\xi }_{\mathrm{h}}(log{R}_{\mathrm{e}}-{\mu }_{\mathrm{R}}(z,M_{\ast }))\hfill \\ \hfill & +{\nu }_{\mathrm{h}}(C-{\mu }_{\mathrm{C}}(z,M_{\ast },R_{\mathrm{e}})),\hfill \end{array}$$(35)

and dispersion σ_h.

In Table 5 we report the median and 68% credible region of the marginal posterior probability distribution of each model parameter. Figure 7 shows the marginal posterior of the parameter ν_h, which in the context of this analysis quantifies the correlation between halo mass and colour at a fixed stellar mass, redshift, and size (blue histogram). The same figure shows the value of ν_h obtained from the analysis of Sect. 4 (pink histogram).

Fig. 7. Colour dependence of the halo mass. Blue histogram: marginal posterior probability distribution of the parameter νh, describing the dependence of the average halo mass as a function of excess colour at a fixed redshift, stellar mass, and size, as defined in Eq. (35). Pink histogram: marginal posterior probability distribution of the parameter νh, describing the dependence of the average halo mass as a function of excess colour at a fixed luminosity and size, as inferred in Sect. 4.

Compared to the analysis of Sect. 4, the inferred value of ν_h is lower, but still consistent with no correlation: ν_h = −1.2 ± 1.2. This result was obtained under the assumption that the stellar mass-to-light ratio has a perfect positive correlation with colour, and therefore represents an extreme case. The true amplitude of the halo mass-colour correlation at a fixed stellar mass is somewhere between the value obtained using luminosity and that reported in this section. In order to further clarify this issue, however, it is necessary to obtain robust estimates of galaxy stellar mass-to-light ratios. As we explained in Sect. 1, this is currently very challenging.

6.2. Possible sources of systematic errors

Our work consists of measuring the distribution of halo mass as a function of various galaxy properties. We put particular care in ensuring that these properties, including redshift, surface brightness profile, and colour, were measured robustly. Nevertheless, some of these measurements could still be affected by systematic errors.

One possible source of error is the presence of colour gradients in the galaxies in our sample. Massive galaxies are known to have, on average, negative colour gradients (that is, they are redder in the inner parts, see for example Tortora et al. 2010; Szomoru et al. 2013); however, we implicitly assumed spatially uniform colours when deriving sizes, luminosities, and rest-frame colours. This implies that our measurements are sensitive to the choice of the photometric band used for the Sérsic profile fits, the r-band, and the method used to measure colours, which relies on the GAaP method of Kuijken et al. (2015).

A negative colour gradient means that the size of a galaxy is smaller when measured at longer wavelengths. This can introduce a spurious correlation between size and redshift because the observed-frame r-band corresponds to longer wavelengths for galaxies at higher z. Such a correlation is absorbed by the parameter ζ_R of our model, which describes the scaling of size with redshift (see Eq. (6)). However, if gradients are a strong function of luminosity, then that would alter the distribution of galaxies in luminosity-size space, potentially introducing a bias in our estimate of parameters β_h and ξ_h. Similarly, our measurement of the trend of halo mass with colour is robust with respect to gradients, as long as these gradients scale with redshift, luminosity, and size in a way that is captured by the model distribution of Eq. (8). However, if gradients correlate with halo mass in a way that is not captured by our model, then biases are introduced.

The correct approach to dealing with colour gradients would be to adopt a more flexible model for the surface brightness distribution of galaxies, with a definition for galaxy size that is independent of the photometric band used, and with the presence of gradients. The model for the halo mass distribution, then, would need to allow for a scaling of halo mass with the colour gradient of a galaxy. That, however, is beyond the scope of this work.

Another possible source of systematic error is the assumption of an NFW density profile with a mass-concentration relation for the dark matter on which our measurement of the halo mass distribution is based. The choice of the density profile is degenerate with the halo mass: for example, assuming a lower concentration leads to a higher inference on M_h (see Sect. 5.3 of Taylor et al. 2020). In principle, then, a positive trend between halo mass and a galaxy observable could escape our detection if this observable also has an anti-correlation with concentration. It is then important, when comparing our measurements with model predictions, to define the halo masses of the model in as similar a way as possible to ours: for example, by fitting the projected surface mass density profile in the radial range 50 − 500 kpc with an NFW profile and a mass-concentration relation given by Eq. (14).

We also assumed a spherical distribution for the dark matter mass. This is, on average, a good assumption, but individual halos are close to elliptical in projection (see for example Georgiou et al. 2021; Schrabback et al. 2021). To first approximation, ignoring ellipticity worsens the goodness of fit of the lensing data around individual lenses. Our model can compensate for this effect by increasing the intrinsic scatter in the halo mass.

We interpreted the lensing signal under the assumption that all of the galaxies in the sample are at the centre of their host halo. Although we removed probable satellite galaxies with the use of the redMaPPer cluster catalogue, there could still be residual satellites in the sample due to the incompleteness of the catalogue. Nevertheless, Sonnenfeld & Leauthaud (2018) show that the inference method on which our work is based is robust with respect to the presence of satellites: for example, a satellite fraction as large as 16% causes a negligible bias on the inferred average halo mass.

The stellar contribution to the lensing signal was estimated by using stellar masses obtained via stellar population synthesis modelling. As we discussed above, these stellar masses are subject to systematic uncertainties. We repeated the analysis by increasing the stellar masses by 0.25 dex, which corresponds to the difference in stellar mass-to-light ratio between a Salpeter and a Chabrier IMF. The main change that we found is that the average halo mass parameter μ_h, 0 decreases by 0.05, compared to the fiducial analysis.

Finally, we ignored the contribution of gas to the total mass. This means that, implicitly, the gas is described by the dark matter component, in our model. We expect the fraction of gas to total mass to be smaller than 10% and roughly constant with halo mass (Dvorkin & Rephaeli 2015), and therefore it is not significant for the results of our study.

6.3. Comparison with the Zu & Mandelbaum (2018) study

Zu & Mandelbaum (2018), using weak lensing and galaxy clustering, report a positive trend between halo mass and rest-frame g − r colour at a fixed stellar mass for a sample of colour-selected red galaxies from SDSS. They also show how their measurement favours a halo-quenching model, a scenario in which the probability of finding a quenched (that is, red) galaxy in a given halo primarily correlates with the halo mass, as opposed to a model in which the galaxy colour primarily depends on the halo assembly time. Confirming or questioning the trend with our independent measurement can therefore have important implications for galaxy evolution models.

Zu & Mandelbaum (2018) do not explicitly report a distribution of halo mass as a function of colour, and therefore it is difficult to make a direct comparison with our results. Nevertheless, their Fig. 7 shows the weak lensing signal, ΔΣ, measured in three stellar mass bins for two subsets of galaxies split on the basis of their g − r colour. Their data show that, in their highest stellar mass bin, the values of ΔΣ for redder galaxies are ∼50% higher than those of bluer galaxies over the radial range covered by our measurements. At a fixed physical distance from the lens centre and fixed halo concentration, ΔΣ roughly scales with a power 0.6 − 0.7 of the halo mass. Therefore, a 50% difference in ΔΣ corresponds to a halo mass difference of 0.3 dex. Our posterior prediction on the difference in average halo mass between our redder and bluer sub-samples rules out such a large value with more than 99% probability, and the discrepancy becomes even larger when using stellar mass in place of luminosity to split the sample. Therefore, our measurements are in strong contrast with those of Zu & Mandelbaum (2018).

There are, however, several differences between our analysis and that of Zu & Mandelbaum (2018). Firstly, Zu & Mandelbaum (2018) used g − r to distinguish among galaxies of a different colour, while we used u − g. However, for red galaxies, such as those in our sample, the two colours are very strongly correlated, and therefore our results are robust with respect to the particular choice of photometric bands. More important differences are the selection criteria and the stellar mass measurements used: we drew galaxies from the Luminous Red Galaxy sample of SDSS, while Zu & Mandelbaum (2018) used the Main Spectroscopic Sample (Strauss et al. 2002). Although there is overlap between the two samples, the latter is dominated by lower-mass and lower-redshift galaxies. Moreover, we applied a morphological selection and obtained luminosities from Sérsic profile fitting to KiDS data, while Zu & Mandelbaum (2018) used colour alone to define their sample of red galaxies and relied on stellar masses obtained from SDSS.

In light of these differences between the two studies, there can be several possible explanations for the origin of the discrepancy. First, the stellar masses used in one of the two studies could be biased in a colour-dependent way: since stellar mass scales with halo mass, any stellar mass-colour correlation would affect the estimate of ν_h (as Sect. 6.1 shows).

Second, a selection based purely on colour, such as that of Zu & Mandelbaum (2018), results in a sample of red galaxies with a non-negligible fraction of lenticular galaxies, bulge-dominated galaxies with spiral arms, or galaxies with an irregular morphology. If halo mass is a function of morphology, such a dependence could show up as a correlation with colour. Such a correlation could also be observed if the stellar mass measurements are biased in a morphology-dependent way. This could happen, for example, if the presence of spiral arms biases the total flux of a galaxy obtained by fitting a simply parameterised surface brightness profile (Sonnenfeld 2022). We tested this hypothesis by repeating our analysis while relaxing our morphology selection criterion. This resulted in the inclusion of an additional 157 galaxies, representing about 6% of the sample. The median value of the marginal posterior on ν_h changed from ν_h = −0.3 to ν_h = 0.0: the effect of including disk galaxies increases the correlation between halo mass and colour, but not to the level required to explain the discrepancy with the Zu & Mandelbaum (2018) measurement. The fraction of red galaxies with a disky morphology, however, is a strong function of stellar mass. Therefore, the morphological selection could still be important in a sample such as that of Zu & Mandelbaum (2018).

Third, the inference method employed here substantially differs from that of Zu & Mandelbaum (2018): while we fit our model directly to the individual shape measurements around each lens, Zu & Mandelbaum (2018) first binned lenses and then fitted the stacked weak lensing data. In order to understand the impact that binning can have on the measurement, we repeated our analysis by splitting the lens sample into two bins according to colour, using the same procedure as in Sect. 4.3, and then fitting for the average halo mass of each sub-sample separately. Since information on the diversity of galaxies within a bin is usually lost in a binned analysis, we did not account for a dependence of the halo mass on any galaxy property: in other words, we set the model parameters ζ_h, β_h, ξ_h, and ν_h to zero.

We found a value of μ_h, 0 = 13.38 ± 0.08 for the redder sub-sample and μ_h, 0 = 13.26 ± 0.10 for the bluer sub-sample. This result tells us that, if we ignore correlations between halo mass and galaxy properties within a bin, our data favour a model in which redder galaxies live in more massive halos on average; although, solutions in which the two sub-samples share the same halo mass are well within 2σ. The fact that the full model instead favours a slightly negative value of ν_h indicates that the apparent difference in average halo mass found in the binned analysis is more likely the result of correlations between halo mass and galaxy properties other than colour, and most importantly luminosity. Nevertheless, this experiment suggests that a difference in the inferred average halo mass of the order of 0.1 dex with respect to the Zu & Mandelbaum (2018) study could be due to the choice of binning the sample.

Fourth, the colour dependence of halo mass could be a strong function of stellar mass disappearing at the very high masses probed by our sample. Lastly, there could also be biases in the colour or weak lensing measurements in either of the two studies.

In order to rule out any of these hypotheses, it would be necessary to repeat our analysis with the same selection criteria, stellar mass measurements, and colour measurements as those used by Zu & Mandelbaum (2018). That, however, is beyond the scope of this work.

7. Conclusions

We carried out a study of the dark matter halo mass distribution of a sample of massive elliptical galaxies, selected from a combination of spectroscopic and morphological criteria. We first measured the distribution in log M_h as a function of redshift, luminosity, half-light radius, and colour. Then we repeated the analysis by describing the surface brightness distribution of galaxies in terms of the luminosity within a physical aperture of 10 kpc, L_r, 10, and the light-weighted projected slope of the surface brightness profile within the same aperture, Γ₁₀. The first description allows for a more intuitive interpretation of the results. The latter is more robust with respect to assumptions on the surface brightness profile of the galaxies.

Our measurements can be used to test theoretical models for the distribution of dark matter halo masses of massive elliptical galaxies, provided that those models can make predictions on galaxy stellar mass-to-light ratios. For models that can predict the full radial profile of the stellar distribution of galaxies, such as hydrodynamical simulations, we recommend using our description of the halo mass distribution in terms of L_r, 10 and Γ₁₀, presented in Sect. 5, as a reference. This would allow one to avoid introducing biases related to the extrapolation of the Sérsic surface brightness model out to very large radii, where it is not constrained by the data.

Our analysis ruled out strong correlations between halo mass and galaxy size at a fixed luminosity, corroborating previous findings from Sonnenfeld et al. (2019a). In particular, the difference in the average halo mass between galaxies at the 84%-ile and 16%-ile of the distribution in half-light radius at a fixed luminosity and redshift is 0.06 ± 0.11 dex (68% credible region). We also ruled out a strong correlation between halo mass and rest-frame u − g colour. The value of the coefficient describing the halo mass-colour correlation slightly depends on whether we adopted a description in terms of luminosity and size, or L_r, 10 and Γ₁₀. In the latter case, the difference in average halo mass between galaxies at the 84%-ile and 16%-ile of the colour distribution at a fixed redshift and surface brightness distribution is 0.00 ± 0.11 dex.

Using stellar mass instead of luminosity does not shift the inferred amplitude of correlations by more than the uncertainty. We then conclude that our analysis did not find evidence for the existence of a correlation between halo mass and a secondary galaxy property at a fixed stellar mass. Elliptical galaxies with a different size or colour for the same stellar mass have presumably had different assembly histories. Our results indicate that, regardless of these differences, the average efficiency of massive elliptical galaxies in converting gas into stars, which is reflected by the stellar-to-halo mass ratio, is roughly constant as a function of size or colour. The origin of the diversity in the size and colour distribution of these objects must be searched for by examining properties other than halo mass.

Acknowledgments

A.S. acknowledges funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 792916. H. Hoekstra acknowledges support from Vici grant 639.043.512, financed by the Netherlands Organisation for Scientific Research (NWO). M.A. acknowledges support from the European Research Council under grant number 647112. M.B. is supported by the Polish National Science Center through grants no. 2020/38/E/ST9/00395, 2018/30/E/ST9/00698 and 2018/31/G/ST9/03388, and by the Polish Ministry of Science and Higher Education through grant DIR/WK/2018/12. C.H. acknowledges support from the European Research Council under grant number 647112, and support from the Max Planck Society and the Alexander von Humboldt Foundation in the framework of the Max Planck-Humboldt Research Award endowed by the Federal Ministry of Education and Research. H. Hildebrandt is supported by a Heisenberg grant of the Deutsche Forschungsgemeinschaft (Hi 1495/5-1) as well as an ERC Consolidator Grant (No. 770935). K.K. acknowledges support from the Royal Society and Imperial College. N.R.N. acknowledges financial support from the “One hundred top talent program of Sun Yat-sen University” grant N. 71000-18841229. AHW is supported by an European Research Council Consolidator Grant (No. 770935). Based on data products from observations made with ESO Telescopes at the La Silla Paranal Observatory under program IDs 177.A-3016, 177.A-3017 and 177.A-3018, and on data products produced by Target/OmegaCEN, INAFOACN, INAF-OAPD and the KiDS production team, on behalf of the KiDS consortium. OmegaCEN and the KiDS production team acknowledge support by Deutsche Forschungsgemeinschaft, ERC, NOVA and NWO-M grants. Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the US Department of Energy Office of Science. The SDSS-III website is http://www.sdss3.org. SDSS-III is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS-III Collaboration including the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Instituto de Astrofisica de Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University. Author contribution: All authors contributed to the development and writing of this paper. The authorship list is given in two groups: the lead authors (AS, CT, HHo), followed by an alphabetical group which includes those who have either made a significant contribution to the data products, or to the scientific analysis.

References

All Tables

All Figures

	Fig. 1. Distribution in redshift, rest-frame r-band luminosity, half-light radius, Sérsic index, rest-frame u − g colour, luminosity within 10 kpc, and light-weighted surface brightness slope within 10 kpc of the galaxy sample. Contour levels enclose 68% and 95% of the sample, respectively.
In the text

Fig. 2. Comparison between observed and model-predicted sizes and colours. Left: observed half-light radius as a function of the average Re of galaxies with the same redshift and luminosity, as predicted by Eq. (6). Right: observed rest-frame u − g colour as a function of the average colour of galaxies with the same redshift, luminosity, and size, as predicted by Eq. (8). Uncertainties on the model parameters are neglected, as they are very small with respect to the observational scatter. The dashed line marks the 1 : 1 relation.

In the text

	Fig. 3. Posterior probability distribution of the model parameters describing the distribution in dark matter halo mass of the sample of lenses defined in Eq. (10). Contour levels correspond to a 68% and 95% enclosed probability.
In the text

Fig. 4. Stacked excess surface mass density profile of the lenses in our sample. The sample is split into two bins according to size (left panel), colour (middle panel), and predicted average halo mass, depending on the maximum a posteriori value of μh defined in Eq. (24) (right panel). The error bars show the excess surface mass density obtained by stacking the observed tangential ellipticity in five radial bins. The bands show the same quantity obtained by posterior prediction from our model, without adding shape noise. For each sub-sample, the figure reports the probability of obtaining a scatter around the posterior predicted ΔΣ profile larger than the observed value, under the assumption that the model describes the truth. Small values for this quantity may indicate that the model is unable to reproduce these observations. Error bars have been offset slightly in the horizontal direction, for clarity.

In the text

	Fig. 5. Posterior probability distribution of the model parameters describing the distribution in halo mass as a function of redshift, Lr, 10 and Γ10, Eq. (30). Contour levels correspond to a 68% and 95% enclosed probability.
In the text

Fig. 6. Stellar mass-to-light ratio as a function of rest-frame u − g colour, for three composite stellar population models obtained with GALAXEV (Bruzual & Charlot 2003). Each model has a fixed metallicity and an exponentially declining star formation history with a decay time of 1 Gyr. Each curve was obtained by varying the time since the initial burst. The grey lines connect points of fixed time since the initial burst, as indicated by the labels. The model corresponding to a metallicity Z = 0.05 is used in Sect. 6.1 to assign stellar masses to galaxies, given their observed colour.

In the text

Fig. 7. Colour dependence of the halo mass. Blue histogram: marginal posterior probability distribution of the parameter νh, describing the dependence of the average halo mass as a function of excess colour at a fixed redshift, stellar mass, and size, as defined in Eq. (35). Pink histogram: marginal posterior probability distribution of the parameter νh, describing the dependence of the average halo mass as a function of excess colour at a fixed luminosity and size, as inferred in Sect. 4.

In the text