# Evidence for radial variations in the stellar mass-to-light ratio of massive galaxies from weak and strong lensing

Alessandro Sonnenfeld, Alexie Leauthaud, Matthew W. Auger, Raphael Gavazzi, Tommaso Treu, Surhud More, Yutaka Komiyama

► To cite this version:

Alessandro Sonnenfeld, Alexie Leauthaud, Matthew W. Auger, Raphael Gavazzi, Tommaso Treu, et al.. Evidence for radial variations in the stellar mass-to-light ratio of massive galaxies from weak and strong lensing. Monthly Notices of the Royal Astronomical Society, 2018, 481 (1), pp.164-184. 10.1093/mnras/sty2262. hal-01703713

HAL Id: hal-01703713

https://hal.science/hal-01703713v1

Submitted on 10 Jun 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Evidence for radial variations in the stellar mass-to-light ratio of massive galaxies from weak and strong lensing

Alessandro Sonnenfeld ID,<fnref n="1" />★ Alexie Leauthaud,<fnref n="1" />,<fnref n="2" /> Matthew W. Auger,<fnref n="3" />

Raphael Gavazzi,<fnref n="4" /> Tommaso Treu,<fnref n="5" /> Surhud More<fnref n="1" /> and Yutaka Komiyama<fnref n="6" />,<fnref n="7" />

Accepted 2018 August 14. Received 2018 August 09; in original form 2018 January 4

## ABSTRACT

The initial mass function (IMF) for massive galaxies can be constrained by combining stellar dynamics with strong gravitational lensing. However, this method is limited by degeneracies between the density profile of dark matter and the stellar mass-to-light ratio (M/L). In this work, we reduce this degeneracy by combining weak lensing together with strong lensing and stellar kinematics. Our analysis is based on two galaxy samples: 45 strong lenses from the SLACS survey and 1700 massive quiescent galaxies from the Sloan Digital Sky Survey main spectroscopic sample with weak lensing measurements from the Hyper Suprime-Cam survey. We use a Bayesian hierarchical approach to jointly model all three observables. We fit the data with models of varying complexity and show that a model with a radial gradient in the stellar M/L is required to simultaneously describe both galaxy samples. This result is driven by a subset of strong lenses with very steep total density profile that cannot be fitted by models with no gradient. Our measurements are unable to determine whether M*/L gradients are due to variations in stellar population parameters at fixed IMF, or to gradients in the IMF itself. The inclusion of M*/L gradients decreases dramatically the inferred IMF normalization, compared to previous lensing-based studies, with the exact value depending on the assumed dark matter profile. The main effect of strong lensing selection is to shift the stellar mass distribution towards the high-mass end, while the halo mass and stellar IMF distribution at fixed stellar mass are not significantly affected.

**Key words:** galaxies: elliptical and lenticular, cD–gravitational lensing: strong–gravitational lensing: weak.

## 1 INTRODUCTION

The determination of the stellar initial mass function (IMF) of massive early-type galaxies (ETGs) is an outstanding problems in astrophysics. Better constraints on the stellar IMF of massive galaxies will help to improve our understanding of star formation in extreme environments, such as the cores of massive galaxies. Moreover, the stellar IMF is currently the main systematic uncertainty in the determination of galaxy stellar masses. This limits our ability to match the numerical simulations, where the most precise predictions are

made for galaxy mass, and the real Universe, where the most robust observations are those regarding galaxy light. Similarly, our ability to compare theoretical predictions with observations is limited by our ignorance of the dark matter distribution inside and around galaxies.

In recent years, there have been a large number of studies aimed at determining the stellar IMF of massive ETGs (Treu et al. 2010; van Dokkum & Conroy 2010; Cappellari et al. 2012; Dutton et al. 2013; Ferreras et al. 2013; Spiniello et al. 2014; and references therein). One of the observational tools that can be used to measure the IMF of galaxies and their dark matter content at cosmological distances is strong gravitational lensing. In a fundamental study, Auger et al. (2010a) explored the degeneracy between the choice of

<footnote n="1">¹Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan</footnote>

<footnote n="2">²Department of Astronomy and Astrophysics, UCO/Lick Observatory, University of California, 1156 High Street, Santa Cruz, CA 95064, USA</footnote>

<footnote n="3">³Institute of Astronomy, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="4">⁴Institut d'Astrophysique de Paris, UMR 7095 CNRS & Université Pierre et Marie Curie, 98 bis Bd Arago, F-75014 Paris, France</footnote>

<footnote n="5">⁵Department of Physics and Astronomy, 430 Portola Plaza, Los Angeles, CA 90095, USA</footnote>

<footnote n="6">⁶National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan</footnote>

<footnote n="7">⁷SOKENDAI(The Graduate University for Advanced Studies), Mitaka, Tokyo 181-8588, Japan</footnote>

<footnote>* E-mail: alessandro.sonnenfeld@ipmu.jp</footnote>
the IMF and the parametrization of the density profile of dark matter halo using strong lensing, weak lensing, and stellar kinematics data for a sample of 59 galaxies from the Sloan Lens ACS Survey (SLACS Bolton et al. 2006; Auger et al. 2010b). Auger et al. (2010a) found that an IMF normalization slightly heavier than that of a Salpeter IMF and increasing with stellar mass are required to fit the data if a Navarro Frenk & White (NFW; Navarro, Frenk & White 1997) profile is assumed. In this work, we revisit that conclusion.

We use the same data as Auger et al. (2010a) work and complement it with weak lensing and stellar kinematics observations for a set of 1700 massive galaxies from the Sloan Digital Sky Survey (SDSS) spectroscopic sample and shear data from the Hyper Suprime-Cam Subaru Strategic Program (HSC SSP; Aihara et al. 2018, hereafter ‘HSC survey’). Our goal is to find the simplest model for the distribution of dark matter and the IMF of massive galaxies that can simultaneously reproduce all of our observations.

Our weak lensing measurements greatly help to pin down the mean halo mass of the sample, thereby allowing us to partially break the degeneracy between dark matter density profile and stellar IMF. However, it also introduces an additional challenge. The sample of SLACS lenses is too small to obtain high signal-to-noise measurements of weak lensing. We therefore use a large sample of galaxies that are drawn from a similar population as the strong lensing sample, but which are themselves not necessarily strong lensing systems. In order to model both samples simultaneously, we need to model the selection effects of the strong lensing sample.

We adopt the following strategy. We use a Bayesian hierarchical inference method to infer the distribution of halo masses and the normalization of the IMF as a function of stellar mass for SLACS lenses and HSC galaxies independently. We then use the model that best describes the HSC data to generate a mock population of galaxies. We apply a strong lensing selection to this mock sample, and compare its properties with those inferred from the SLACS sample. If the assumed model is accurate, then the prediction from both the HSC sample, and the SLACS sample should match up. We explore three different models, with increasing degrees of complexity, until we reach such an agreement.

This paper is organized as follows. In Section 2, we present the data. In Section 3, we introduce our Bayesian hierarchical inference method. In Section 4, we fit three different models to the data. We discuss our results in Section 5 and conclude in Section 6. We assume a flat Λ cold dark matter (ΛCDM) cosmology with $\Omega_M = 0.3$ and $H_0 = 70 \text{ km s}^{-1} \text{Mpc}^{-1}$. Halo masses, labelled as $M_h$, are defined as the mass enclosed within a spherical shell with average density equal to 200 times the critical density of the Universe. Halo radius is labelled $R_h$ and is also defined assuming a boundary defined with respect to 200 times the critical density. Halo concentration is labelled as $c_h$ and is defined as $c_h = R_h/R_s$, where $R_s$ is the halo scale radius. Stellar and halo masses are expressed in units of solar mass, while sizes and projected distances are expressed in physical units.

# 2 THE DATA

## 2.1 SLACS strong lensing

Our sample of strong lens systems is drawn from the set of 59 SLACS lenses studied by Auger et al. (2010b). The selection of SLACS lenses is described in detail by Bolton et al. (2006) and is briefly summarized here. SLACS lenses were discovered by scanning spectra of quiescent galaxies from the SDSS main spectroscopic sample and looking for multiple emission lines associated

with objects at a higher redshift than the main target galaxy. Strong lens candidates obtained with this method were followed-up with high-resolution imaging from *Hubble Space Telescope* (HST) to confirm their lens nature and to measure their properties.

Auger et al. (2009) present a sample of 84 strong lenses discovered with this procedure. Of these lenses, six are morphologically classified as disc galaxies, and the remaining are classified as ETGs. The Auger et al. (2010b) sample, from which we select the lenses for our study, consists of a subset of 59 lenses drawn from the Auger et al. (2009) sample. These have been selected to have early-type morphology, robust lens models, and reliable stellar kinematics measurements.

For each of these 59 lenses, the following data are available:

(i) Lens and source redshifts and lens velocity dispersions within the SDSS spectroscopic fibre $\sigma_{ap}$ (Auger et al. 2009).

(ii) Rest-frame *I*-band effective radii obtained from de Vaucouleurs profile fit to multiband HST imaging data (Auger et al. 2009).

(iii) Stellar masses obtained from spectral energy distribution (SED) fitting to HST imaging data, assuming a Chabrier IMF, $M_{*}^{(Chab)}$ (Auger et al. 2009).

(iv) Einstein radii $R_{Ein}$ (Auger et al. 2010b).

In addition, for a subset of 33 galaxies, weak lensing shear measurements, obtained with deep data taken with the Advanced Camera for Surveys on HST, are available. We refer to Gavazzi et al. (2007) for a description of these weak lensing measurements.

## 2.2 HSC massive galaxies

The Hyper Suprime-Cam (HSC; Furusawa, Koike & Takata 2018; Kawanomoto et al. 2018; Komiyama, Obuchi & Nakaya 2018; Miyazaki, Komiyama & Kawanomoto 2018) is an optical camera with a 1.5 deg² field of view, installed on the Subaru Telescope. The HSC survey (Aihara et al. 2018) is an ongoing five-band *grizy* photometric survey. The wide layer will cover 1400 deg² to 26.2 mag and with a mean seeing of ~0.6 arcsec in the *i* band.

We wish to build a sample of galaxies selected in a similar way as the parent sample used to build the SLACS sample, in the region of the sky covered by HSC data. We start from the SDSS main spectroscopic sample and select objects with a quiescent spectrum by imposing a maximum H α emission equivalent width of 3 Å. We then eliminate highly elongated objects by selecting objects with a ratio between minor and major axis larger than 0.5, as measured by a fit of a de Vaucouleurs profile (de Vaucouleurs 1948) to the SDSS *r* band. This second step is designed to approximate the selection of objects with an early-type morphology, as done by Auger et al. (2010b). Since we plan to use weak lensing measurements for this sample of galaxies, we also apply a cut in redshift, selecting only objects at $z > 0.05$, so that each galaxy contributes with a non-negligible lensing signal. For the resulting sample of galaxies, we then consider the following data:

(i) Stellar masses from the MPA JHU catalogue (Kauffmann et al. 2003), obtained through stellar population synthesis fits to SDSS broadband photometry and assuming a Chabrier IMF.

(ii) Half-light radii obtained from a de Vaucouleurs profile fit to the SDSS *r* band as provided by SDSS DR12 (Alam et al. 2015).

(iii) Stellar velocity dispersion measurements within the SDSS spectroscopic fibre.
(iv) Weak lensing measurements from HSC (described in the next subsection).

Shape measurements and photo-z measurements for background weak lensing source galaxies are the only data from the HSC survey that we use in this work. Galaxy sizes and stellar masses are obtained from SDSS. Although HSC photometry is deeper and has better image quality than the SDSS, we do not expect it to bring a significant improvement on the measurements of our sample, since these low-redshift massive galaxies are well resolved and detected with high signal-to-noise in SDSS data.

We wish to carry out a weak lensing analysis around the galaxies in this sample, using HSC shear data. We will be using the Bayesian hierarchical inference method developed by Sonnenfeld & Leauthaud (2018). The Sonnenfeld & Leauthaud (2018) method relies on a critical assumption, that of *isolated lenses*: lens galaxies are assumed to be at an infinite projected distance from each other, so that each background source is only lensed by one lens, at most. This approximation was shown by Sonnenfeld & Leauthaud (2018) to be valid in the high-mass regime, around lens galaxies with $log M_* > 11$, but likely breaks down at low masses. We therefore apply a cut in stellar mass to our sample and select galaxies with observed stellar mass larger than $10^{11} M_{\odot}$. We also eliminate galaxies identified as satellites in the redMaPPer cluster catalogue (Rykoff et al. 2014). In particular, we reject objects with a cluster membership probability larger than 50 per cent and a central galaxy probability smaller than 10 per cent. This step eliminates ~5 per cent of the objects. Finally, we remove extreme outliers by fitting a mass-size and mass-velocity dispersion relation to the sample, and removing objects at more than $4\sigma$ away from the mean relation, with an iterative sigma-clipping procedure. These outliers could bias the inference, and most likely correspond to faulty measurements in at least one of the three variables considered: mass, size, or velocity dispersion.

Within the resulting sample, ~2000 galaxies lie in the region of the sky with weak lensing measurements from HSC. We assign background sources to each lens in the sample. For each lens, we consider sources within a projected distance of 300 kpc at the lens redshift. For some lenses, the area within 300 kpc overlaps with that of other objects. In order to avoid using the same source twice, not allowed under the assumption of isolated lenses, we eliminate lenses from the sample until the areas within 300 kpc of the remaining lenses no longer overlap. We do this iteratively, starting from the objects with the smallest observed stellar mass. Around 15 per cent of the lenses are eliminated with this procedure, leaving a final sample of ~1700 galaxies. These are massive ($log M_*^{(Chab)} > 11$) and quiescent galaxies, belonging to the SDSS main spectroscopic sample, with weak lensing measurements from HSC. Throughout this paper, we refer to this sample of galaxies as the 'control sample'.

The control sample has been obtained by applying cuts to the SDSS main spectroscopic sample, motivated by the need for a sample of massive galaxies with a small contamination from satellites, and clear of more massive neighbours within the region used for the weak lensing analysis. For consistency, we apply the same criteria to the SLACS sample. This reduces the number of strong lenses available for the analysis to 45 objects, with respect to the initial 59. Of the 14 rejected lenses, 1 is eliminated exclusively for its low stellar mass ($log M_*^{(Chab)} < 11.0$), while 13 have more massive galaxies in their proximity. Four of these lenses are also identified as satellites by redMaPPer.

Since we will be combining data from this sample with higher quality measurements from the SLACS sample, it is important to make sure that there are no systematic biases between the two

<figure><img src="image_2.png" /><figcaption>**Figure 1.** *Top:* Difference between the SDSS-based stellar mass of SLACS lenses, from MPA-JHU, as a function of the *HST*-based measurements on the same objects by Auger et al. (2009). *Bottom:* Relative difference between the observed *r*-band effective radii from SDSS DR12 as a function of rest-frame I-band effective radii measured with *HST* data by Auger et al. (2009).</figcaption></figure>

data sets. In particular, we want to check if stellar mass and size measurements obtained using ground-based data produce the same values as the measurements obtained from space. This is a relatively straightforward test, since SLACS lenses are in SDSS. In Fig. 1, we show SDSS stellar masses and effective radii of SLACS lenses as a function of the same quantities measured using *HST* data. There is good agreement between the two data sets, particularly in terms of stellar mass. Effective radii derived from SDSS are on average 8 per cent larger than *HST*-based measurements, with a 14 per cent scatter. Since these are strong lenses, though, it is possible that bias and scatter are not fully representative of the typical galaxies in the control sample, due to the presence of a lensed source blended with the main galaxy that could affect the measurement of $R_e$. In any case, we will discuss the implications of such a small discrepancy at the end of Section 4.1.

Finally, we show in Fig. 2 the distribution in observed stellar mass, effective radius, velocity dispersion, redshift, and SDSS rest-frame *g* − *i* colour of the SLACS and the control sample. SLACS galaxies appear to be slightly more compact compared to the control sample, at the low-mass end of the distribution. Their velocity dispersion is correspondingly higher. The distribution of the two samples in the $M_*-z$ plane is similar, with more massive galaxies being preferentially at higher redshifts, a result of the luminosity selection of the SDSS spectroscopic sample, from which both lens sets are drawn. Given these tests, throughout the rest of this paper, we will assume that the SLACS sample and the control sample are drawn from the same population, modulo the strong lensing selection that will be described in Section 3.6.

## 2.3 HSC weak lensing and photo-z

We use weak lensing measurements from the first internal data release (DR1) of the HSC survey shear catalogue (Mandelbaum et al. 2018). The catalogue covers an area of 137 deg$^2$, with a number density of sources of 25 arcmin$^{-2}$ (unweighted). Shape measurements have been obtained on *i*-band images, using the re-Gaussianization (point spread function) correction method of Hirata & Seljak (2003). Calibration for multiplicative and additive bias have been obtained through a set of dedicated simulations, similar
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Observed stellar mass versus effective radius, aperture velocity dispersion, redshift and rest-frame g − i colour, for SLACS lenses and galaxies in the control sample. The dashed lines show the best-fitting relation obtained for the two samples.</figcaption></figure>

to those of the GREAT 3 challenge (Mandelbaum et al. 2014, 2015). We refer to Mandelbaum et al. (2018) for further details.

Photometric redshifts for the background sources are taken from the HSC DR1 (Tanaka et al. 2018). In particular, we use photo-zs obtained through the template fitting-based code MIZUKI (Tanaka 2015). For each lens, we select sources for which the lower 2σ bound on the photo-z is larger than the lens redshift. In other words, we select sources with a 97.5 per cent probability of being in the background of our lenses. For each source, we take the median value of the photo-z probability distribution as the fiducial value for the source redshift. For simplicity, we do not propagate redshift uncertainties to the shear estimate, as this is sub-dominant with respect to the shape noise (see Sonnenfeld & Leauthaud 2018 for a discussion).

## 2.4 Stacked weak lensing

Before introducing the full model used in our analysis we carry out a simple comparison between the stacked weak lensing signal for SLACS lenses and the control sample. The stacked signal for SLACS lenses is obtained by combining, in each radial bin, the measurements of ΔΣ for all the objects in the sample, and converting them to our fiducial cosmology. For the control sample, we combine the signal for all of the ~1700 lenses, but we apply a weighting scheme to match the stellar mass distribution of the HSC galaxies to that of the SLACS sample. In particular, for a given galaxy with

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Stacked weak lensing signal around SLACS lenses and for our control sample of 1700 massive quiescent galaxies with HSC weak lensing. The control sample is matched to the SLACS sample by using a stellar mass-based weighting scheme.</figcaption></figure>

observed stellar mass $M_{*}^{(Chab)}$, we apply the following weight:

$$w_i = \frac{1}{\sqrt{2\pi}\sigma_{*,(\text{SLACS})}} \exp\left\{-\frac{(\log M_{*}^{(\text{Chab})} - \mu_{*,(\text{SLACS})})^2}{2\sigma_{*,(\text{SLACS})}^2}\right\}, \quad (1)$$

where $\mu_{*,(\text{SLACS})}$ and $\sigma_{*,(\text{SLACS})}$ are the mean and standard deviation of the observed log-stellar mass of the SLACS sample. We use the HSC weak lensing pipeline (More et al. in preparation) to calculate ΔΣ in radial bins. The stacked lensing signal for the SLACS and the control sample is shown in Fig. 3. There is good agreement between the two samples. Since the stacked weak lensing signal is mainly sensitive to the halo-mass distribution, this indicates that SLACS lenses have similar halo masses as non-lens galaxies with the same stellar mass. However, because of the small number of SLACS lenses, the SLACS weak lensing measurement is noisy. This prevents us from making more quantitative statements.

## 3 THEORETICAL FRAMEWORK

We wish to describe the distribution of dark matter halos and stellar IMF at the massive end of the galaxy population. We do so with a Bayesian hierarchical inference method. The method takes elements from the analysis of Sonnenfeld et al. (2015), who applied the Bayesian hierarchical inference formalism to strong lensing and stellar kinematic measurements for 80 galaxies, including SLACS lenses, and the more recent formulation of Sonnenfeld & Leauthaud (2018), who introduced a similar method to interpret weak lensing measurements.

One complication in combining the SLACS sample and the control sample is that the two samples are in principle sampling different regions of the distribution in stellar mass, halo mass, and stellar IMF. Due to strong lensing selection (Ruff et al. 2011), SLACS lenses sample the high-mass and high-velocity dispersion end of this distribution, as opposed to our sample of non-lenses, which we defined by applying a lower cut-off in stellar mass. In order to obtain an accurate inference on the halo mass and stellar IMF distribution, it is essential to (1) choose a model that accurately describes how these quantity vary across the region of parameter space spanned by the two samples, and (2) accurately model the strong lensing selection in order to compare the results inferred from the two samples.
We explore different models. In the first model, labelled the ‘vanilla’ model, we assume that dark matter haloes are described by an NFW profile and that halo mass and the normalization of the stellar IMF scale with a power-law dependence on stellar mass. This model is described below. We then generalize the vanilla model in two different ways. The first variation consists in modifying the in-ner density profile of the dark matter halo by allowing for adiabatic contraction or expansion. The second variation instead consists in allowing for a gradient in the stellar mass-to-light ratio (M/L), while keeping the dark matter halo profile fixed to an NFW profile. The adiabatic contraction model and the M*/L gradient model will be described in Sections 4.2 and 4.3, respectively.

## 3.1 Bayesian hierarchical inference

Let us first introduce the concept of Bayesian hierarchical inference in its most basic form. Let us consider a sample of galaxies, each one fully described by the value of a parameter set $\psi$. Let us assume that the individual values of the parameters $\psi$ are drawn from a dis-tribution described by a set of hyper-parameters, $\eta$, the values of which we wish to infer. For a given prior, $P(\eta)$, the posterior prob-ability distribution function (PDF) of the hyper-parameters given a set of measurements $d$ on the galaxy sample is given by Bayes theorem as

$$P(\eta|\mathbf{d}) \propto P(\eta)P(\mathbf{d}|\eta), \quad (2)$$

where $P(\mathbf{d}|\eta)$ is the likelihood of observing the data given the values of the hyper-parameters. If measurements on different galaxies are independent from each other, this can be expanded as

$$P(\mathbf{d}|\eta) = \prod_{i} P(\mathbf{d}_i|\eta). \quad (3)$$

This is not true in general for weak lensing measurements. In principle, each background source is lensed by every foreground galaxy. Equation (3), however, is valid under the approximation of isolated lenses, on which our analysis is based. We treat lens galaxies as if they were at infinite distance from each other, so that the likelihood of a set of shape measurements around a given lens is independent of the model describing other lenses.

Finally, each term in the product can be evaluated by marginal-izing over the individual parameters $\psi$ as

$$P(\mathbf{d}_i|\eta) = \int d\psi_i P(\mathbf{d}_i|\psi_i)P(\psi_i|\eta). \quad (4)$$

Equations (2)–(4) allow us to make an inference on the values of the hyper-parameters $\eta$ by evaluating the posterior PDF $P(\eta|\mathbf{d})$. There is a two-level hierarchy of parameters in this problem: the hyper-parameters $\eta$, at the top level, specify the distribution of lower-level parameters $\psi$, describing individual objects, which are needed to calculate the likelihood.

## 3.2 Individual galaxy parameters

We now apply this formalism to our specific science case. We describe the mass distribution of each galaxy as the sum of two components. A spherically de-projected de Vaucouleurs profile de-scribing the stars, and a spherical NFW profile for the dark matter halo. We then introduce the IMF mismatch parameter, defined as the ratio between the true stellar mass of a galaxy, and the stellar mass an observer would infer assuming a Chabrier IMF and having otherwise perfect knowledge of the stellar population parameters of

the galaxy, $M_*^{(\text{Chab})}$:

$$\alpha_{\text{IMF}} \equiv \frac{M_*}{M_*^{(\text{Chab})}}. \quad (5)$$

We use the following parameters to describe each galaxy:

$$\psi \equiv \{M_*^{(\text{Chab})}, \alpha_{\text{IMF}}, M_h, c_h\}, \quad (6)$$

where $M_h$ is the halo mass and $c_h$ the halo concentration.

## 3.3 Hyper-parameters of the galaxy population distribution

At this point, we need to specify a form for the distribution of galaxy parameters, $\psi$ described by a given set of hyper-parameters. We assume the following form:

$$P(\psi|\eta) = S(M_*^{(\text{Chab})})H(M_h|M_*^{(\text{Chab})}) \times C(c_h|M_h)I(\alpha_{\text{IMF}}|M_*^{(\text{Chab})}). \quad (7)$$

The term $S(M_*^{(\text{Chab})})$ describes the distribution in stellar mass. Following Sonnenfeld & Leauthaud (2018), we use a skew Gaussian distribution, which is shown to be appropriate for the description of a sample obtained by applying a stellar mass cut:

$$S(M_*^{(\text{Chab})}) = \frac{1}{\sqrt{2\pi}\sigma_*} \exp \left\{ -\frac{(\log M_*^{(\text{Chab})} - \mu_*)^2}{2\sigma_*^2} \right\} \times \Phi(\log M_*^{(\text{Chab})}), \quad (8)$$

with

$$\Phi(\log M_*^{(\text{Chab})}) = 1 + \text{erf} \left( s_* \frac{\log M_*^{(\text{Chab})} - \mu_*}{\sqrt{2}\sigma_*} \right). \quad (9)$$

The values of the parameters $\mu_*$, $\sigma_*$, and $s_*$ will naturally be dif-ferent for the two samples of galaxies, the SLACS lenses and the control sample.

The next term in equation 7 is $H(M_h|M_*^{(\text{Chab})})$, describing the dis-tribution of halo masses as a function of stellar mass, also modelled as a Gaussian:

$$H(M_h|M_*^{(\text{Chab})}) = \frac{1}{\sqrt{2\pi}\sigma_h} \exp \left\{ -\frac{(\log M_h - \mu_h)^2}{2\sigma_h^2} \right\}. \quad (10)$$

As anticipated, we allow the mean of this Gaussian to be a function of stellar mass:

$$\mu_h = \mu_{h,0} + \beta_h(\log M_*^{(\text{Chab})} - 11.3). \quad (11)$$

The halo mass–concentration relation is modelled as

$$P(c_h|M_h) = \frac{1}{\sqrt{2\pi}\sigma_c} \exp \left\{ -\frac{(\log c_h - \mu_c(M_h))^2}{2\sigma_c^2} \right\} \quad (12)$$

with

$$\mu_c(M_h) = c_{h,0} + \beta_c(\log M_h - 12). \quad (13)$$

We fix the values of the parameters of the mass–concentration rela-tion to $c_{h,0} = 0.830$, $\beta_c = -0.098$, and $\sigma_c = 0.1$, based on results from numerical simulations by Macciò, Dutton & van den Bosch (2008).

Finally we introduce a distribution for $\alpha_{\text{IMF}}$:

$$I(\alpha_{\text{IMF}}|M_*^{(\text{Chab})}) = \frac{1}{\sqrt{2\pi}\sigma_{\text{IMF}}} \exp \left\{ -\frac{(\log \alpha_{\text{IMF}} - \mu_{\text{IMF}})^2}{2\sigma_{\text{IMF}}^2} \right\}. \quad (14)$$

Similarly to the halo mass distribution, we allow the mean of this Gaussian to scale with stellar mass

$$\mu_{\text{IMF}} = \mu_{\text{IMF},0} + \beta_{\text{IMF}}(\log M_*^{(\text{Chab})} - 11.3). \quad (15)$$
We ignore redshift evolution in the model parameters, since the redshift range spanned by the lens samples is quite small. In summary, the full list of hyper-parameters is

$$ \eta = \{\mu_{*}, \sigma_{*}, s_{*}, \mu_{h,0}, \beta_{h}, \mu_{\text{IMF},0}, \beta_{\text{IMF}}, \sigma_{\text{IMF}}\}. \quad (16) $$

A brief description of each parameter is provided in Table 1.

## 3.4 The likelihood term

In order to evaluate the posterior PDF of the hyper-parameters we need to calculate the likelihood term $P(\mathbf{d}|\psi)$ for each galaxy. The data consist of the following:

(i) Observed stellar mass from stellar population synthesis, $M_{*}^{(\text{obs})}$.

(ii) Line-of-sight velocity dispersion measured within the SDSS spectroscopic aperture, $\sigma_{\text{ap}}$.

(iii) Weak lensing shape measurements, {WL}

(iv) Einstein radius $R_{\text{Ein}}$, defined as the radius of a circular aperture within which the mean surface mass density equals the critical density of the lens-source pair (SLACS lenses only).

The likelihood is then given by

$$ \begin{aligned} P(\mathbf{d}|\psi) = & P(M_{*}^{(\text{obs})}|M_{*}^{(\text{Chab})})P(\sigma_{\text{ap}}|M_{*}^{(\text{Chab})}, \alpha_{\text{IMF}}, M_{\text{h}}, c_{\text{h}}) \\ & \times P(\{\text{WL}\}|M_{*}^{(\text{Chab})}, \alpha_{\text{IMF}}, M_{\text{h}}, c_{\text{h}}) \\ & \times P(R_{\text{Ein}}|M_{*}^{(\text{Chab})}, \alpha_{\text{IMF}}, M_{\text{h}}, c_{\text{h}}). \end{aligned} \quad (17) $$

The stellar mass term $P(M_{*}^{(\text{obs})}|M_{*}^{(\text{Chab})})$ is a Gaussian in $M_{*}^{(\text{Chab})}$. The velocity dispersion term is evaluated using the spherical Jeans equation to construct a model seeing-convolved luminosity-weighted line-of-sight velocity dispersion within the spectroscopic aperture. The strong and weak lensing terms are calculated by producing model mass density profiles composed of the sum of a de Vaucouleurs profile of mass $\alpha_{\text{IMF}}M_{*}^{(\text{Chab})}$ for the stars and a dark matter halo, assuming that the two have the same centre and no other mass component contributes to the lensing signal in the region probed by the data.

The likelihood for HSC weak lensing is given by the following product over individual lensed sources:

$$ P(\{\text{WL}\} | M_{*}^{(\text{Chab})}, \alpha_{\text{IMF}}, M_{\text{h}}, c_{\text{h}}) = \prod_{i} P(e_{t,i} | M_{*}^{(\text{Chab})}, \alpha_{\text{IMF}}, M_{\text{h}}, c_{\text{h}}). \quad (18) $$

Here, $e_{t,i}$ indicates the tangential component of the observed ellipticity of the $i$th source, relative to the lens. Unlike ellipticity measurements based on other methods, the expectation value of $e_{t,i}$ defined with the re-Gaussianization method of Hirata & Seljak (2003) is not equal to the tangential component of the shear, $g_t$. Instead, the following identity holds:

$$ 2Rg_t = \langle e_t \rangle, \quad (19) $$

where R is called the ‘shear responsivity’, $R \approx 1 - e_{\text{rms}}^2$, and $e_{\text{rms}}$ is the intrinsic RMS scatter in source ellipticities. equation (19) is valid for unbiased ellipticity estimates. In practice, the expectation value of the ellipticity is corrected with a multiplicative and additive bias, as follows:

$$ 2R[(1+m)g_t + c_t] = \langle e_t \rangle, \quad (20) $$

where $c_t$ is the tangential component of the additive bias. m and c are calibrated on simulations, as described by Mandelbaum et al. (2018). The likelihood of observing a given value of the tangential

ellipticity given the model is then

$$ P(e_t | M_{*}^{(\text{Chab})}, \alpha_{\text{IMF}}, M_{\text{h}}, c_{\text{h}}) = \frac{1}{\sqrt{2\pi\sigma_{e,\text{tot}}^2}} \exp \left\{ -\frac{(2R[(1+m)g_t + c_t] - e_t)^2}{2\sigma_{e,\text{tot}}^2} \right\}, \quad (21) $$

where $\sigma_{e,\text{tot}}^2$ is the sum in quadrature of the intrinsic shape noise $e_{\text{rms}}$ and the shape measurement error, $\sigma_e$. The reduced shear $g_t$ is a function of the model parameters, as well as the redshifts of lens and sources.

## 3.5 Goodness of fit estimation

We now discuss how the goodness of fit of our model is assessed. In Bayesian hierarchical inference methods, this is usually done through *posterior predictive tests*: we use the posterior probability distribution to generate mock data, which we then compare with the observed data. In principle, the goodness of fit can then be determined by defining a metric quantifying how well the model prediction matches the observations. For simplicity, however, we will limit our evaluation to a qualitative comparison.

In our case, the relevant data consists of weak lensing shape measurements, Einstein radii and central velocity dispersions. For both galaxy samples, we will attempt to reproduce their stacked weak lensing signal. For the strong lens sample, we will also predict Einstein radii and stellar velocity dispersion measurements. We will first compare predicted and observed values of $R_{\text{Ein}}$ and $\sigma_{\text{ap}}$ separately, and then we will combine these to obtain a measurement of the slope of the total density profile, $\gamma'$. The total density slope has been measured for lenses in the SLACS sample by fitting power-law density profiles, $\rho(r) \propto r^{-\gamma'}$, to the Einstein radius and velocity dispersion measurements (Koopmans et al. 2006; Auger et al. 2010b). We will use the same procedure to determine posterior predicted values of $\gamma'$.

The density slope has a more intuitive physical meaning compared to the Einstein radius. Whereas $R_{\text{Ein}}$ indicates the total mass enclosed within an aperture of arbitrary radius, determined by the geometry of the lens-source system, the total density slope $\gamma'$ provides information on the internal structure of a lens. Although the value $\gamma'$ is not a pure observable, and depends somewhat on the assumption of a power-law density profile (see Xu et al. 2017 for a detailed study), its use should be seen as a way of mapping observables from the two-dimensional space defined by $R_{\text{Ein}}$ and $\sigma_{\text{ap}}$ to a one-dimensional variable. At fixed $\sigma_{\text{ap}}$, increasing $R_{\text{Ein}}$ corresponds to a smaller value of $\gamma'$ (a shallower density profile). At fixed $R_{\text{Ein}}$, increasing $\sigma_{\text{ap}}$ results in a larger $\gamma'$ (a steeper density profile).

In practice, we proceed as follows. We first draw a subset of points from the MCMC sample of the posterior probability distribution for a given model. For each point, corresponding to a set of hyper-parameter values, we generate a mock lens population, with values of stellar mass, stellar IMF, halo mass and concentration drawn from the model distribution. Additionally, we draw values for the lens redshift distribution, which we model as a Gaussian centred at 0.18 with 0.07 dispersion, which is a good approximation for both samples. Finally, we draw values for effective radii from a mass-size relation. This is modelled as a Gaussian in $R_e$ with mean given by

$$ \mu_R = \mu_{R,0} + \beta_R (\log M_*^{(\text{Chab})} - 11.3), \quad (22) $$

with $\mu_{R,0} = 0.85$, $\beta_R = 0.65$, and dispersion $\sigma_R = 0.13$. These values are determined by fitting the mass-size relation of the control sample.
<table><caption>**Table 1.** Vanilla model. Median values and 68 per cent confidence interval of the posterior probability distribution of individual hyper-parameters, marginalized over the rest of the hyper-parameters. The third column lists values of the hyper-parameters obtained by applying a strong lensing selection correction to a model population corresponding to the maximum-likelihood of the control sample inference, as described in Section 3.6. In the last row, we list the inferred median stellar mass of the different samples.</caption><thead><tr><th></th><th>SLACS</th><th>Control</th><th>Control, SL pred.</th><th>Parameter description</th></tr></thead><tbody><tr><td>μ<sub>h, 0</sub></td><td>12.31 ± 0.17</td><td>13.02 ± 0.04</td><td>13.05</td><td>Average log M<sub>h</sub> at stellar mass log M<sub>*</sub><sup>(Chab)</sup> = 11.3</td></tr><tr><td>σ<sub>h</sub></td><td>0.52 ± 0.15</td><td>0.30 ± 0.05</td><td>0.31</td><td>Dispersion in log M<sub>h</sub> around the average</td></tr><tr><td>β<sub>h</sub></td><td>0.64 ± 0.91</td><td>1.51 ± 0.15</td><td>1.56</td><td>Power-law dependence of halo mass on M<sub>*</sub><sup>(Chab)</sup></td></tr><tr><td>μ<sub>IMF, 0</sub></td><td>0.32 ± 0.03</td><td>0.27 ± 0.01</td><td>0.30</td><td>Average log α<sub>IMF</sub> at stellar mass log M<sub>*</sub><sup>(Chab)</sup> = 11.3</td></tr><tr><td>σ<sub>IMF</sub></td><td>0.05 ± 0.02</td><td>0.08 ± 0.01</td><td>0.08</td><td>Dispersion in log α<sub>IMF</sub> around the average</td></tr><tr><td>β<sub>IMF</sub></td><td>0.29 ± 0.13</td><td>0.08 ± 0.02</td><td>0.08</td><td>Power-law dependence of IMF normalization on M<sub>*</sub><sup>(Chab)</sup></td></tr><tr><td>μ<sub>*</sub></td><td>11.29 ± 0.09</td><td>11.23 ± 0.03</td><td>11.38</td><td>Average parameter in equation (8)</td></tr><tr><td>σ<sub>*</sub></td><td>0.22 ± 0.05</td><td>0.23 ± 0.02</td><td>0.22</td><td>Dispersion parameter in equation (8)</td></tr><tr><td>log s<sub>*</sub></td><td>-0.54 ± 0.85</td><td>-0.21 ± 0.19</td><td>-0.89</td><td>Log of the skewness parameter in equation (9)</td></tr><tr><td>Median log M<sub>*</sub></td><td>11.37 ± 0.03</td><td>11.33 ± 0.01</td><td>11.41</td><td>Median stellar mass (not a hyper-parameter)</td></tr></tbody></table>

Given a mock lens population, we compute the average $\Delta\Sigma$ at the location of the radial bins used in Fig. 3. For consistency with Section 2.4, we apply the same stellar mass-based weighting scheme to $\Delta\Sigma$ for the mock based on the inference from the control sample.

For the SLACS-based mock, we also generate a population of background sources. We draw source redshifts from a Gaussian distribution centred at 0.60 and with a dispersion of 0.18, which is a good approximation of the redshift distribution of SLACS strongly-lensed sources. We then compute Einstein radii and, by solving the spherical Jeans equation, the seeing-convolved surface brightness-weighted line-of-sight velocity dispersion within a circular aperture of 1.5 arcsec radius. Finally, we fit a power-law density profile to the mock Einstein radius and velocity dispersion data to obtain the density slope $\gamma'$.

We then compare the distribution in $\Delta\Sigma$, $R_{Ein}$, $\sigma_{ap}$, and $\gamma'$ obtained from mocks generated, using 1000 randomly drawn points from the MCMC chain with the observed stacked weak lensing and density slope measurements.

## 3.6 The strong lensing selection function

In order to compare the inferences from both galaxy samples, strong lensing selection effects need to be accounted for. There are two separate selections at play. The first is due to different galaxies having different cross-sections for strong lensing. This is relatively easy to model, as will be shown later. The second is due to some strong lenses being easier to discover than others, depending on their properties and on the efficiency of the lens search algorithm.

In principle, strong lensing selection can be explicitly modelled in our Bayesian formalism. This would be achieved by describing the distribution of SLACS galaxies as the product of the distribution of the control sample times a term proportional to the product between the lensing cross-section and the lens finding efficiency, which skews it towards the position in parameter space occupied by strong lenses. This approach was used in the strong lensing and stellar dynamics analysis of Sonnenfeld et al. (2015). However, the relatively high dimensionality of our model makes it technically challenging to implement this approach here. Instead, we carry out a posterior predictive test. We start from the maximum likelihood model inferred from the control sample, which we assume to be the truth, and use it to draw a mock population. We populate these galaxies with background sources with the same redshift distribution as that of SLACS lensed sources and calculate the strong lensing cross-section and lens finding efficiency of each lens-source pair. We then

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Schematic representation of the strategy adopted to compare the inferences from the SLACS sample and the control sample, while taking into account strong lensing selection effects.</figcaption></figure>

draw a random subsample of galaxies, weighting each galaxy by its strong lensing cross section and lens finding efficiency. This subsample is meant to simulate a sample of strong lenses. Finally, we fit our model to this subsample. The inferred hyper-parameters will represent the distribution of SLACS galaxies as predicted from the control sample. A schematic representation of the process is shown in Fig. 4.

Strictly speaking, the procedure we just described is not a posterior predictive test because we only focus on the maximum-likelihood set of hyper-parameters. In principle, one should perform a posterior predictive measurement at each point of the hyper-parameter space sampled by the control sample posterior probability distribution, and obtain a new distribution of hyper-parameters to be compared with the SLACS inference. In practice, this would require too much computational time. We then only consider how
the lensing selection shifts the peak of the posterior, and assume
that the same shift applies to the whole distribution.

We define the strong lens cross-section as the area in the source plane that gets mapped into sets of at least two multiple images with a minimum magnification of $|\mu_{\min}| = 0.5$. This last criterion is chosen in order to exclude lens configurations with very de-magnified images, although in principle these are still strong lens systems, in practice they cannot be identified as such. We assume circular symmetry for the cross-section computation. We verified that the results are stable with respect to variations in the value of $\mu_{\min}$.

The lens finding efficiency is more difficult to model. SLACS lenses were found by means of a spectroscopic search for emission lines from lensed background sources, using the 1.5 arcsec radius SDSS fibre. Qualitatively, we understand how the lens finding efficiency varies as a function of the Einstein radius: it goes to zero for very small values of $R_{Ein}$, since the flux from the lensed source is correspondingly very small, peaks for values of $R_{Ein}$ close to the fibre radius, and goes down for very large values of $R_{Ein}$, since most of the flux falls outside of the spectroscopic fibre for these lensed sources. However, a quantitative assessment of the lens finding efficiency requires dedicated simulations, including the creation of synthetic spectra, the exploration of different mock lens configurations, and modelling the effects of the atmospheric seeing.

Arneson, Brownstein & Bolton (2012) carried out such a study in the context of power-law lens models. Their main result is that, for a SLACS-like survey, the lens finding efficiency is zero for very small values of $R_{Ein}$, is approximately constant for values of the Einstein radius up to $R_{Ein} \simeq 2.0$ arcsec, then falls with increasing $R_{Ein}$ for lenses with a density profile steeper than isothermal (see the middle left-hand panel of fig. 5 of Arneson et al. 2012). With this result in mind, we approximate the lens finding efficiency as a constant for values of $R_{Ein}$ down to 0.5 arcsec and zero for smaller values. This is a good approximation, since most of the Einstein radii of the population of mock strong lenses have values $R_{Ein} < 2.0$ arcsec. We verified that none of our results change when the lens finding efficiency is set to zero for $R_{Ein} > 2.0$ arcsec.

# 4 RESULTS

## 4.1 NFW halos and spatially constant IMF

We fit the model described in Section 3 to the observed stellar masses, velocity dispersions, background source shapes. For SLACS galaxies, we also use Einstein radii. The SLACS sample and the control sample are fit separately. The posterior PDF is evaluated by running a Markov chain Monte Carlo. We assume uniform priors on all hyper-parameters except for $s_*$, For $s_*$, we assume a uniform prior on its logarithm in the range (−1, 1). Integrals over individual lens parameters in equation (4) are computed with Monte Carlo integration and importance sampling, following Sonnenfeld & Leauthaud (2018; see also Schneider et al. 2015).

Fig. 5 displays the posterior PDF relative to the hyper-parameters describing the halo mass and stellar IMF distribution, inferred from the SLACS and control samples separately. Median values with 68 per cent confidence region of the inference on the full list of hyper-parameters are reported in Table 1. The values of the hyper-parameters resulting from applying the strong lensing selection correction to the maximum-likelihood model for the control sample, as described in Section 3.6, are plotted as black triangles in Fig. 5, and listed in Table 1. The shifts with respect to the maximum-likelihood model due to lensing selection are minimal. The largest effect is

seen for the average IMF normalization parameter, $\mu_{\text{IMF},0}$, which increases by ~0.05 and is brought into better agreement with the SLACS inference. However, lensing selection has a very small effect on halo mass. The reason for this is that the projected dark matter mass within the Einstein radius of galaxies, typically on the scale of 5 kpc, is a very mild function of halo mass.

Let us first consider the inference on the hyper-parameters obtained from the control sample. The derived average in $\log M_h$ at the pivot stellar mass $\log M_*^{(\text{Chab})} = 11.3$ is $\mu_{h,0} = 13.02 \pm 0.04$, while the average of the distribution in $\log \alpha_{\text{IMF}}$ is $\mu_{\text{IMF},0} = 0.27 \pm 0.01$. This result is consistent with the stacked weak lensing measurement of Schulz, Mandelbaum & Padmanabhan (2010), who concluded that a stellar IMF consistent with a Salpeter IMF (which would correspond to $\log \alpha_{\text{IMF}} = 0.25$) is needed to match the dynamical masses (i.e. the velocity dispersion measurements) of a sample of SDSS galaxies selected in a similar way to our control sample.

The inference based on the SLACS sample is for the most part consistent with that based on the control sample, with one important exception: halo mass. SLACS data favours a very low value for the average halo mass: $\mu_{h,0} = 12.31 \pm 0.17$. This value is not only inconsistent with the inference from HSC data, but, together with the inferred value of the IMF normalization $\mu_{\text{IMF},0} = 0.32$, would correspond to a baryonic fraction of 15 per cent of the total mass, very close to the cosmological value of ~16 per cent (Planck Collaboration XIII 2016), and therefore unphysical.

This result might seem in contradiction with the plot of Fig. 3, in which the stacked weak lensing signal measured for SLACS lenses is observed to be consistent with that of the control sample, once matched in stellar mass. However, when fitting our model to the whole set of available data for the SLACS lenses, the likelihood is mostly determined by the constraints on small scales, strong lensing and stellar kinematics, as the precision on the Einstein radius and velocity dispersion measurements is much higher compared to the weak lensing measurements. The low value of the inferred average halo mass, then, indicates that the model is not a good description of the structure of massive galaxies at all the scales probed by the data.

The discrepancy between the inferences for the two data sets is further illustrated in the top panel of Fig. 6, where we plot the average stellar, dark matter, and total density profile inferred for a galaxy of stellar mass $\log M_*^{(\text{Chab})} = 11.3$ and effective radius $R_e = 7$ kpc. The inner density profile is quite similar in the two cases. This is because, for both data sets, the model is fit to the SDSS velocity dispersion, which is mostly sensitive to the total density profile on the scales of a few kpc. On the other hand, the two density profiles diverge drastically at large radii, due to the large difference in the inferred halo mass.

Although in principle, for this comparison, we should be considering the average density profile obtained after applying the strong lensing selection correction to the control sample, the differences with respect to the uncorrected profile are minimal, and we do not plot it in Fig. 6 to minimize confusion.

Let us now consider the goodness of fit. In the top panel of Fig. 7, we plot the posterior predictive stacked weak lensing signal, generated as described in Section 3.5, compared with the observed data. For the control sample, the model is a good description of the data. The same can be said for the model based on the SLACS sample: although the inferred halo mass is very low, the predicted $\Delta\Sigma$ profile goes through the data points, in virtue of the large uncertainties of the SLACS weak lensing measurements.

Fig. 8 focuses on small-scale observables: the Einstein radius, the velocity dispersion and the inner slope of the total density profile, $\gamma'$,
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Vanilla model. Posterior probability distribution of the parameters describing the distribution in halo mass and IMF mismatch parameter. Red region: inference from the control sample. Blue lines: inference from SLACS galaxies. Black triangles: values of the hyper-parameters obtained by applying a strong lensing selection correction to the maximum-likelihood control sample model, as described in Section 3.6.</figcaption></figure>

derived from the combination of the two as described in Section 3.5. As can be seen from the first two plots in the top row, the SLACS-based model produces a reasonable distribution in Einstein radius and velocity dispersion (blue histograms), compared to the observed distributions (black histograms).

The corresponding density slope $\gamma'$ is plotted in the upper right-hand panel as a function of stellar mass density $\Sigma_* = M_{*}^{(\text{Chab})}/(2\pi R_c^2)$, which has been shown to correlate with $\gamma'$ (Sonnenfeld et al. 2013). The posterior predictive distribution of $\gamma'$ based on the SLACS sample seems to reproduce reasonably well the observed distribution, including the trend with $\Sigma_*$. In this respect, we recover the result of Shankar et al. (2017), who showed that the correlation of the density slope with stellar mass and size can be explained with NFW + de Vaucouleurs models and a simple stellar-to-halo mass relation. However, there are a few SLACS lenses that cannot be described by the posterior predictive $\gamma'$ distribution. These are objects with $\log \Sigma_* \approx 9.0$ and large values of the density slope $\gamma' \approx 2.3$. The vanilla model is unable to reproduce such steep density profiles.

This is particularly the case for the values of the hyper-parameters derived using HSC data (the control sample): the predictive distribution generated from the maximum-likelihood point of the HSC-based inference, corrected for strong lensing selection effects and shown as red circles in the upper right-hand panel of Fig. 8, covers only the low-$\gamma'$ end of the observed distribution. This is another

way of understanding why a simple NFW + de Vaucouleurs model is not a good description of the data: if the model is tuned to reproduce the weak lensing signal at large radii, the predicted inner density profile is too shallow. The low values of $\gamma'$ are in turn a result of the predicted Einstein radii being too large compared to the observed distribution, as can be seen by comparing the red and black histograms in the upper left-hand panel of Fig. 8: for a given value of the velocity dispersion, a larger Einstein radius corresponds to a shallower density profile, and vice versa.

The statistical errors on the hyper-parameters are significantly smaller for the inference based on the control sample, compared to the SLACS sample. This is especially the case for the hyper-parameter describing the average stellar IMF normalization, $\mu_{\text{IMF,0}}$, for which we obtained a 0.01 dex uncertainty. The stellar IMF is constrained, for the control sample, mostly by the SDSS velocity dispersion measurements. This is because the velocity dispersion is mostly sensitive to the mass distribution within the region probed by the stars, which is dominated by the baryons. Dark matter parameters, on the other hand, are mostly constrained by weak lensing data, at larger radii. Our ‘vanilla’ model has only one degree of freedom in the stellar component: its mass, given by the product $M_{*}^{(\text{Chab})}\alpha_{\text{IMF}}$. The model is then well constrained by the data, given that observational uncertainties on the velocity dispersion are on the order of 10 km s$^{-1}$ (~5 per cent relative uncertainty on $\sigma_{\text{ap}}$). The combination of 1700 measurements then brings the statistical
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Average density profile of a galaxy with log M<sub>*</sub><sup>(Chab)</sup> = 11.3 and R<sub>e</sub> = 7 kpc inferred by fitting the vanilla (top), adiabatic contraction (middle), and M<sub>*</sub>/L gradient (bottom) model to the SLACS (blue) and control sample (red) data sets. Shaded region: 68 per cent confidence region on the total density profile. Dashed lines: 84 per cent and 16 per cent enclosed probability limits on the dark matter density. Dotted lines: 84 per cent and 16 per cent enclosed probability limits on the stellar mass density. In each panel, the inset shows in greater detail the region of the plot delimited by the rectangle.</figcaption></figure>

error on the mean IMF normalization to a very small value. How-
ever, the true uncertainty on the IMF normalization parameter is
most likely dominated by systematic uncertainties related to our
particular model choice. This will become evident in Section 4.3.

One possible source of systematic uncertainty is the measure-
ment of the effective radius. We have shown in Fig. 1 how the
SDSS-based values of R_e, which we use for our control sample,
appear to be systematically larger by 8 per cent compared to HST-
based measurements on the same galaxies. Assuming this offset to

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Posterior predictive tests based on the stacked weak lensing profile, to assess the goodness of fit of the vanilla (top), adiabatic contraction (middle), and M<sub>*</sub>/L gradient (bottom) model. Observed data points from Fig. 3 are shown as error bars. The shaded bands correspond to regions of 68 per cent enclosed probability in the average ΔΣ. In the case of the control sample, both the posterior predictive and the observed distribution are obtained by applying a stellar mass weighting scheme based on the stellar mass distribution of SLACS lenses. The shaded region indicates the radial range not used to constrain the model.</figcaption></figure>

be real, and not due to contamination from the lensed background
source (we used SLACS lenses for this comparison), we can esti-
mate the impact on the inference on α<sub>IMF</sub> with a simple dynamical
argument. The IMF normalization that we derive from the control
sample is essentially determined by the ratio between a dynamical
mass, obtained from the velocity dispersion measurements, and a
stellar mass obtained from photometry. From the virial theorem,
the dynamical mass scales with the effective radius and the velocity
dispersion as follows:

$$
M_{\text{dyn}} \sim \frac{R_e \sigma^2}{G}. \tag{23}
$$
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Posterior predictive test on strong lensing and stellar kinematics observables, to assess the goodness of fit of the vanilla (top), adiabatic contraction (middle), and *M*<sub>*</sub>/L gradient (bottom) model. **Left-hand panel:** Distribution in Einstein radius. The black histogram is the observed distribution from Auger et al. (2009), the blue one is obtained from the posterior inferred using SLACS data, the red one is obtained from the peak of the posterior obtained from the control sample, corrected by strong lensing selection (i.e. weighted by strong lensing cross-section). Histograms are normalized to the total number of SLACS lenses. **Middle panel:** Distribution in SDSS velocity dispersion. Different colours correspond to the observed distribution of SLACS lenses from Auger et al. (2009) (black), the observed distribution of the control sample (grey), and the posterior predicted distributions based on the SLACS (blue) and control sample (red) inferences. **Right-hand panel:** Total density slope, obtained by fitting a power-law density profile to R<sub>Ein</sub> and σ<sub>ap</sub>, as a function of stellar mass density. Black error bars correspond to values obtained by Auger et al. (2010b) on the SLACS sample. The shaded region is the distribution obtained by drawing samples of galaxies from the posterior given by the SLACS inference. Red circles are samples drawn from the peak of the control sample posterior, corrected by strong lensing selection.</figcaption></figure>

Therefore, errors on R_e propagate linearly to the dynamical mass, and, to first approximation, to the IMF normalization. An 8 per cent systematic error on the effective radius then corresponds to a similar error on α IMF. We therefore assume a systematic uncertainty of ~0.03 on the hyper-parameter μ IMF,0, to be added in quadrature to the statistical uncertainty.

## 4.2 Modified dark matter profile

The vanilla model produces an unphysical solution when fitted to SLACS data. This prompts us to propose a different model. We then relax the assumption of an NFW profile for the dark matter halo and allow for the effect of adiabatic contraction or expansion. We define adiabatically contracted/expanded haloes following Dutton
<figure><img src="image_10.png" /><figcaption>Figure 9. Density profile of adiabatically contracted/expanded halo for different values of the response efficiency parameter $v$. For $v = 0$, the profile is that of an NFW halo. In this example the halo mass is $\log M_h = 13$, the halo concentration is $c_h = 5$, the stellar mass is $\log M_* = 11.5$, and the effective radius is 7 or 4 kpc. The density profile of the baryonic component for the two effective radii is also plotted.</figcaption></figure>

et al. (2007). For a given halo mass, stellar mass and effective radius, a fully contracted halo profile is calculated using the adiabatic contraction prescription of Blumenthal et al. (1986). The Blumenthal et al. (1986) treatment assumes that baryons and dark matter share initially the same density profile, which we assume to be an NFW halo with a standard mass-concentration relation. As the baryons contract to the observed distribution (which we describe with a spherically de-projected de Vaucouleurs profile), the orbits of dark matter particles are modified while conserving the product $rM(r)$, where $r$ is the distance from the centre of the potential and $M(r)$ is the total mass enclosed within the spherical shell of radius $r$. The product $rM(r)$ is an adiabatic invariant in the case of circular orbits.

Starting from the fully contracted dark matter profile we then build a family of profiles by introducing the contraction efficiency parameter $v$, which is allowed to range between 1, corresponding to a fully contracted profile à la Blumenthal et al. (1986) and $-1$, corresponding to an expanded profile.

The resulting dark matter density profile for different values of the parameter $v$ and for two different values of the galaxy effective radius is plotted in Fig. 9. The effects of adiabatic contraction or expansion are visible only in the inner ~50 kpc. We also point out how, for the same halo mass, stellar mass and contraction efficiency parameter $v$, the effect of adiabatic contraction is stronger for more compact galaxies.

We model the distribution of this new parameter, $v$, as a Gaussian, truncated between $-1 < v < 1$, which we multiply to equation (7). The new probability distribution of the individual galaxy parameters given the hyper-parameters becomes

$$P(\psi|\eta) = S(M_*^{(\text{Chab})})H(M_h|M_*^{(\text{Chab})}) \times C(c_h|M_h)\bar{\mathcal{I}}(\alpha_{\text{IMF}}|M_*^{(\text{Chab})})\mathcal{N}(v), \quad (24)$$

with

$$\mathcal{N}(v) = \frac{A_v}{\sqrt{2\pi}\sigma_v} \exp \left\{ -\frac{(v - \mu_v)^2}{2\sigma_v^2} \right\}; \quad -1 < v < 1. \quad (25)$$

The new hyper-parameters $\mu_v$ and $\sigma_v$ are the mean and dispersion of this truncated Gaussian, on which the normalization constant

$A_v$ depends. We assume the following exponential prior on the parameter $\sigma_v$

$$P(\sigma_v) \propto \exp\left(-\frac{\sigma_v}{0.1}\right), \quad (26)$$

in order to penalize solutions with very large values of the intrinsic scatter in the adiabatic contraction efficiency, which would correspond to a flat distribution in $v$. In other words, we assert a model in which the distribution in the adiabatic contraction efficiency is centred around a well-defined mean.

We fit this new model to SLACS data and the control sample separately. The inference on the hyper-parameters describing halo mass, stellar IMF and adiabatic contraction efficiency is plotted in Fig. 10, while median and 68 per cent confidence interval of all hyper-parameters are listed in Table 2.

Allowing for adiabatic contraction or expansion does help bringing the average halo mass inferred from SLACS lenses in a $2\sigma$ agreement with the control sample inference. However, a careful look at the $\mu_v - \mu_h$ degeneracy contour (first column, fourth row in Fig. 10) shows that strong lensing data allow for reasonably large average halo masses only for negative values of the adiabatic contraction efficiency parameter, corresponding to adiabatic expansion. The control sample, instead, favours slight contraction. The inference on the IMF normalization is also correspondingly different: much higher for SLACS galaxies.

The discrepancy between the two inferences can also be seen by looking at the density profile, plotted in the middle panel of Fig. 6: although the total density profile of the two models agree, the inferred central dark matter density differs by at least a factor of a few between the two data sets. Once again, our model is unable to simultaneously describe both data sets with the same values of the hyper-parameters.

The middle panels of Figs 7 and 8 show posterior predictive tests for this model, through which we can check the goodness of fit. Compared to the vanilla model, in this case the model based on SLACS data provides a better description of the stacked weak lensing signal. The adiabatic contraction model, however, is still unable to describe the inner density profile of a subset of SLACS lenses with relatively large values of $\gamma'$, as can be seen from the middle right-hand panel of Fig. 8. These objects, when fitted with a bulge + halo model, prefer solutions with very little dark matter (see e.g. Posacki et al. 2015). In fact, even mass-follows-light models for some of these lenses under-predict their measured stellar velocity dispersion. As a reference, a strong lens consisting only of a de Vaucouleurs profile (no dark matter) corresponds to $\gamma' \approx 2.2$ (see red triangles in fig. 6 of Sonnenfeld et al. 2013), a value exceeded by 16 out of the 45 strong lenses in our sample.

Adiabatically contracted or expanded dark matter profiles do not help improve the match with the data in these cases, because this family of models can only produce dark matter profiles that are shallower than the stellar component, as can be seen from Fig. 9. As a result, the data favours solutions with as little dark matter in the inner regions as possible, either via a low average halo mass or through adiabatic expansion. We have verified this conjecture by repeating the fit over a smaller sample of strong lenses, obtained by eliminating the 16 objects with a power-law density slope $\gamma'$ steeper than 2.2. In this case, we find values of the hyper-parameters that are consistent with those inferred from the control sample, indicating that the problem lies in the objects with a steep density profile.

In principle, we could allow for even greater freedom on the choice of the dark matter profile. For instance, we could also vary the halo concentration, which at the moment is fixed by a relatively
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Adiabatic contraction model. Posterior probability distribution of the parameters describing the distribution in halo mass, IMF mismatch parameter and adiabatic contraction efficiency. Red region: Inference from the control sample. Blue lines: Inference from SLACS lenses. Black triangles: Values of the hyper-parameters obtained by applying a strong lensing selection correction to the maximum-likelihood control sample model, as described in Section 3.6.</figcaption></figure>

<table><caption>**Table 2.** Adiabatic contraction model. Median values and 68 per cent confidence interval of the posterior probability distribution of individual hyper-parameters, marginalized over the rest of the hyper-parameters. The third column lists values of the hyper-parameters obtained by applying a strong lensing selection correction to a model population corresponding to the maximum-likelihood of the control sample inference, as described in Section 3.6.</caption><thead><tr><th></th><th>SLACS</th><th>Control</th><th>Control, SL pred.</th><th>Parameter description</th></tr></thead><tbody><tr><td>μ<sub>h, 0</sub></td><td>12.45 ± 0.25</td><td>13.03 ± 0.03</td><td>13.03</td><td>Average log M<sub>h</sub> at stellar mass log M<sub>*</sub><sup>(Chab)</sup> = 11.3</td></tr><tr><td>σ<sub>h</sub></td><td>0.39 ± 0.16</td><td>0.17 ± 0.02</td><td>0.18</td><td>Dispersion in log M<sub>h</sub> around the average</td></tr><tr><td>β<sub>h</sub></td><td>1.28 ± 0.92</td><td>1.29 ± 0.10</td><td>1.42</td><td>Power-law dependence of halo mass on M<sub>*</sub><sup>(Chab)</sup></td></tr><tr><td>μ<sub>ν</sub></td><td>-0.38 ± 0.44</td><td>0.47 ± 0.05</td><td>0.51</td><td>Average adiabatic contraction efficiency parameter</td></tr><tr><td>σ<sub>ν</sub></td><td>0.31 ± 0.20</td><td>0.15 ± 0.02</td><td>0.17</td><td>Dispersion in adiabatic contraction efficiency parameter</td></tr><tr><td>μ<sub>IMF, 0</sub></td><td>0.33 ± 0.03</td><td>0.15 ± 0.02</td><td>0.18</td><td>Average log α<sub>IMF</sub> at stellar mass log M<sub>*</sub><sup>(Chab)</sup> = 11.3</td></tr><tr><td>σ<sub>IMF</sub></td><td>0.05 ± 0.02</td><td>0.09 ± 0.01</td><td>0.08</td><td>Dispersion in log α<sub>IMF</sub> around the average</td></tr><tr><td>β<sub>IMF</sub></td><td>0.26 ± 0.14</td><td>0.01 ± 0.03</td><td>0.01</td><td>Power-law dependence of IMF normalization on M<sub>*</sub><sup>(Chab)</sup></td></tr><tr><td>μ<sub>*</sub></td><td>11.30 ± 0.09</td><td>11.25 ± 0.03</td><td>11.34</td><td>Average parameter in equation (8)</td></tr><tr><td>σ<sub>*</sub></td><td>0.22 ± 0.05</td><td>0.23 ± 0.01</td><td>0.23</td><td>Dispersion parameter in equation (8)</td></tr><tr><td>log s<sub>*</sub></td><td>-0.58 ± 0.82</td><td>-0.32 ± 0.18</td><td>-0.45</td><td>Log of the skewness parameter in equation (9)</td></tr><tr><td>Median log M<sub>*</sub></td><td>11.38 ± 0.03</td><td>11.33 ± 0.01</td><td>11.40</td><td>Median stellar mass (not a hyper-parameter)</td></tr></tbody></table>

narrow prior. We suspect that models for which the inner slope of the dark matter density profile is steeper than the stellar component would finally allow us to obtain a good fit. This could be achieved by setting the scale radius of the dark matter halo to be comparable to the Einstein radius, or, in other words, by increasing the concentration by a factor of 10 or more. However, such steep density profiles would be difficult to justify theoretically: any baryonic

effect additional to adiabatic contraction would tend to flatten the profile.

## 4.3 Gradients in stellar M/L

In this subsection, we introduce a new model, with more freedom in the stellar mass distribution, and with the dark matter model reverted
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Top: Ratio between light-weighted stellar M/L and stellar M/L at the effective radius, as a function of the M/L gradient slope $\gamma$. Centre: Ratio between half-mass radius and half-light radius as a function of gradient slope $\gamma$. Bottom: Central velocity dispersion for a galaxy with stellar mass log $M_*$ = 11.5, halo mass log $M_h$ = 13, and effective radius $R_e$ = 7 kpc, as a function of $M_*/L$ radial slope. The model velocity dispersion is calculated by solving the spherical Jeans equation under the assumption of isotropic orbits.</figcaption></figure>

to an NFW profile. We introduce a new parameter describing a radial gradient in stellar M/L. While until now we have implicitly assumed that the stellar mass follows exactly the light distribution; here, we relax this assumption and allow for the following M/L profile

$$ \Upsilon_*(R) = \Upsilon_{*,e} \left(\frac{R}{R_e}\right)^{\gamma}, \qquad (27) $$

where $\Upsilon_{*,e}$ is the M/L at the effective radius, and $-1 < \gamma < 0$. With this definition, for negative values of $\gamma$ the M/L in the inner regions of the galaxy is larger than in the outskirts. Although this parametrization of the $M_*/L$ profile is unphysical at very small and large radii, where it tends to infinity and zero, respectively, this is not a problem because the contribution to the total budget of the stellar mass in these regions is very small. In other words, the results would not change if we add an upper and lower limit to the M/L. However, the above parametrization allows for a relatively easy implementation of an M/L gradient in our analysis.

There can be two origins for $M_*/L$ gradients. One is the presence of gradients in age, dust, or metallicity across the galaxy. Another possibility is the presence of gradients in stellar IMF. Gradients in stellar population properties at fixed IMF can in principle be measured with spatially resolved spectroscopic data, or estimated with high-quality multiband photometry (see e.g. San Roman et al. 2018). Our stellar population measurements do not provide us with spatially resolved information. The stellar population synthesis-based stellar masses used so far, $M_*^{(Chab)}$, are obtained fitting synthetic stellar population models to total magnitudes of galaxies. The resulting stellar M/L is therefore a light-weighted average over the extent of the whole galaxy. We then define the light-weighted M/L as follows:

$$ \langle \Upsilon_* \rangle_{\text{LW}} = \frac{\int_{0}^{\infty} I(R) R \Upsilon_*(R) dR}{\int_{0}^{\infty} I(R) R dR}, \qquad (28) $$

where $I(R)$ is the surface brightness profile of the galaxy, in our case assumed to be a de Vaucouleurs profile.

In Fig. 11, we plot the relation between $\Upsilon_{*,e}$ and $\langle \Upsilon_* \rangle_{\text{LW}}$. In addition, we compute the ratio between half-mass radius and half-

light radius, as well as the impact of the gradient on the central velocity dispersion, plotted in Fig. 11 as a function of $\gamma$.

This parametrization allows for a steeper total density profile compared to our previous stellar mass-follows-light model. Strong lensing and stellar kinematics data are particularly sensitive to this new parametrization, since it can modify significantly the total density profile in the inner few kpc. As can be seen in Fig. 11, changing the slope of the $M_*/L$ gradient can change significantly the velocity dispersion of a galaxy (see also Bernardi et al. 2018 for a related discussion).

Similarly to the adiabatic contraction efficiency case, we describe the distribution of the gradient parameter $\gamma$ as a truncated Gaussian of mean $\mu_\gamma$ and dispersion $\sigma_\gamma$. The new probability distribution of the individual galaxy parameters given the hyper-parameters becomes

$$ P(\psi|\eta) = S(M_*^{(\text{Chab})})H(M_h|M_*^{(\text{Chab})}) \times C(c_h|M_h)\bar{\mathcal{I}}(\alpha_{\text{IMF}}|M_*^{(\text{Chab})})\mathcal{G}(\gamma), \qquad (29) $$

with

$$ \mathcal{G}(\gamma) = \frac{A_\gamma}{\sqrt{2\pi}\sigma_\gamma} \exp \left\{ -\frac{(\gamma - \mu_\gamma)^2}{2\sigma_\gamma^2} \right\}. \qquad (30) $$

We assume an exponential prior on the intrinsic scatter in the gradient parameter:

$$ P(\sigma_\gamma) \propto \exp \left( -\frac{\sigma_\gamma}{0.1} \right). \qquad (31) $$

We fit the model to SLACS lenses and the control sample. The inference on the hyper-parameters describing halo mass, stellar IMF and $M_*/L$ gradient is plotted in Fig. 12, while the median and 68 per cent enclosed probability values of all hyper-parameters are listed in Table 3. This time there is good agreement between the inference based on SLACS data and on HSC weak lensing data, with or without strong lensing selection correction. Both data sets find the need for a gradient in $M_*/L$, an average log halo mass of $\mu_{h,0} \approx 13$ and an average IMF intermediate between Chabrier and a Salpeter IMF (which would correspond to $\mu_{\text{IMF},0} = 0.25$). The inference is consistent with no mass dependence of the IMF (parameter $\beta_{\text{IMF}} \sim 0$). The agreement between the two inferences can also be seen by looking at the predicted average density profile, plotted in the bottom panel of Fig. 6.

The goodness of fit is also improved, compared to the previous two models, as can be seen by comparing the bottom rows of Figs 7 and 8 with the middle and top panels. In particular, the range of values of the inner density slope $\gamma'$ predicted by the $M_*/L$ gradient model is now a good match to the observations.

By fitting the posterior predicted $\gamma'$ distribution with a Gaussian with mean dependent on stellar surface mass density,

$$ P(\gamma'|\Sigma_*) = \frac{1}{\sqrt{2\pi}\sigma_{\gamma'}} \exp \left\{ -\frac{(\gamma' - \mu_{\gamma'}(\Sigma_*))^2}{2\sigma_{\gamma'}^2} \right\}, \qquad (32) $$

$$ \mu_{\gamma'} = \gamma_0 + \eta_{\gamma'}(\log \Sigma_* - 9.0), \qquad (33) $$

we find an average $\gamma'$ at the pivot surface mass density $\gamma_0 = 2.13$, a dependence of $\gamma'$ on $\Sigma_*$ of $\eta_{\gamma'} = 0.32$, and an intrinsic scatter $\sigma_{\gamma'} = 0.09$. The values inferred by Sonnenfeld et al. (2013) for the SLACS sample are $\gamma_0 = 2.11 \pm 0.02$ (taking into account the redshift-dependence of $\gamma'$), $\eta_{\gamma'} = 0.38 \pm 0.07$, and $\sigma_{\gamma'} = 0.12 \pm 0.02$, in good agreement with the model prediction.

The predicted distribution in Einstein radius is closer to the observed one, with respect to previous models, although still unable to reproduce the sharp peak around 1.2 arcsec. This suggests that
<figure><img src="image_13.png" /><figcaption>**Figure 12.** *M*<sub>*</sub>/L gradient model. Posterior probability distribution of the parameters describing the distribution in halo mass, IMF mismatch parameter, and IMF gradient. *Red region:* Inference from the control sample. *Blue region:* Inference from SLACS galaxies. *Black triangles:* Values of the hyper-parameters obtained by applying a strong lensing selection correction to the maximum-likelihood control sample model, as described in Section 3.6.</figcaption></figure>

<table><caption>**Table 3.** *M*<sub>*</sub>/L gradient model. Median values and 68 per cent confidence interval of the posterior probability distribution of individual hyper-parameters, marginalized over the rest of the hyper-parameters.</caption><thead><tr><th></th><th>SLACS</th><th>Control</th><th>Control, SL pred.</th><th>Parameter description</th></tr></thead><tbody><tr><td>μ<sub>h,0</sub></td><td>13.12 ± 0.19</td><td>13.04 ± 0.04</td><td>13.11</td><td>Average log M<sub>h</sub> at stellar mass log M<sub>*</sub><sup>(Chab)</sup> = 11.3</td></tr><tr><td>σ<sub>h</sub></td><td>0.37 ± 0.16</td><td>0.31 ± 0.04</td><td>0.23</td><td>Dispersion in log M<sub>h</sub> around the average</td></tr><tr><td>β<sub>h</sub></td><td>0.54 ± 0.69</td><td>1.48 ± 0.15</td><td>1.43</td><td>Power-law dependence of halo mass on M<sub>*</sub><sup>(Chab)</sup></td></tr><tr><td>μ<sub>γ</sub></td><td>-0.24 ± 0.04</td><td>-0.26 ± 0.05</td><td>-0.26</td><td>Average M<sub>*</sub>/L gradient parameter</td></tr><tr><td>σ<sub>γ</sub></td><td>0.09 ± 0.02</td><td>0.06 ± 0.01</td><td>0.07</td><td>Dispersion in the gradient parameter</td></tr><tr><td>μ<sub>IMF,0</sub></td><td>0.11 ± 0.04</td><td>0.02 ± 0.04</td><td>0.04</td><td>Average log α<sub>IMF</sub> at stellar mass log M<sub>*</sub><sup>(Chab)</sup> = 11.3</td></tr><tr><td>σ<sub>IMF</sub></td><td>0.05 ± 0.02</td><td>0.06 ± 0.01</td><td>0.05</td><td>Dispersion in log α<sub>IMF</sub> around the average</td></tr><tr><td>β<sub>IMF</sub></td><td>0.23 ± 0.15</td><td>0.03 ± 0.02</td><td>0.03</td><td>Power-law dependence of IMF normalization on M<sub>*</sub><sup>(Chab)</sup></td></tr><tr><td>μ<sub>*</sub></td><td>11.28 ± 0.09</td><td>11.24 ± 0.03</td><td>11.34</td><td>Average parameter in equation (8)</td></tr><tr><td>σ<sub>*</sub></td><td>0.22 ± 0.05</td><td>0.23 ± 0.01</td><td>0.24</td><td>Dispersion parameter in equation (8)</td></tr><tr><td>log s<sub>*</sub></td><td>-0.44 ± 0.86</td><td>-0.31 ± 0.16</td><td>-0.30</td><td>Log of the skewness parameter in equation (9)</td></tr><tr><td>Median log M<sub>*</sub></td><td>11.37 ± 0.03</td><td>11.33 ± 0.01</td><td>11.43</td><td>Median stellar mass (not a hyper-parameter)</td></tr></tbody></table>

either the model or the strong lensing selection correction are not a perfect description of reality. Nevertheless, the *M*<sub>*</sub>/L gradient model performs significantly better compared to the vanilla and the adiabatic contraction, and from now on we will adopt it as our fiducial model.

There is a strong degeneracy between the parameter describing the average stellar IMF normalization, μ<sub>IMF,0</sub> and the one describing the average *M*<sub>*</sub>/L gradient, μ<sub>γ</sub>, for both the SLACS and the control sample data sets (see panel in row 6 and column 4 in Fig. 12). This is because, as we explained in Section 4.1, the parameters describing
the stellar component are mostly constrained by stellar velocity dis-
persion data, for which we only have single aperture measurements
from SDSS. With only one stellar kinematics data point per galaxy,
it is not possible to constrain both the IMF normalization and the
M*,L gradient on individual objects. This degeneracy can be broken
only partially with the statistical combination of measurements over
many galaxies.

5 DISCUSSION

We fitted models for the structure of massive ETGs to two sets
of massive quiescent galaxies, one consisting of strong lenses and
the other one drawn from the general population of galaxies, using
weak lensing, strong lensing and stellar kinematics data. We con-
sidered three different models. The first two models, consisting of
a de Vaucouleurs stellar bulge and a dark matter halo with varying
degrees of freedom in the density profile, fail to self-consistently
describe strong lenses and the sample of massive galaxies, even after
accounting for lensing selection effects. These results highlight the
challenge of finding an accurate description for the density profile
of massive galaxies over the two decades in radius probed by our
data, 1 kpc < R < 300 kpc.

One of the requirements for a successful model is the ability to reproduce the range of total density slopes measured for the strong lenses. About a third of the SLACS lenses have a total density profile steeper than what can be accounted for with mass-follows-light models ($\gamma' > 2.2$). The only dark matter models that would be able to explain these observations are the ones for which the inner density profile is allowed to be steeper than that of the stars. However, we reject such models based on our physical prior that there are no known mechanism that can produce such steep dark matter profiles.

A model for the mass distribution consisting of NFW haloes
and de Vaucouleurs profiles with a radial gradient in stellar M/L
is instead able to describe SLACS strong lenses and the control
sample of massive galaxies with the same sets of hyper-parameters.
The maximum-likelihood value of the average M*/L gradient slope
found for SLACS galaxies is µγ = -0.24. Is this value of the slope
reasonable? How much of this gradient can be attributed to gradients
in stellar population properties at fixed IMF, and how much could
be due to a gradient in IMF?

The galaxies in both our samples show colour gradients: they
are redder in the centre with respect to the outskirts. This can be
seen by looking at the distribution of effective radii measured in
different photometric bands, plotted in Fig. 13: half-light radii in
redder bands are systematically smaller than in bluer bands.

Colour gradients can be the result of a gradient in age, dust con-
tent, or metallicity. Massive ETGs such as the objects in our samples
have typically very little dust. This leaves age and metallicity as the
most probable origins of the gradients. Varying age at fixed metal-
licity or vice versa results in a higher M*/L in correspondence with
redder colour. Therefore, regardless of the interpretation, the colour
gradients observed in our sample are qualitatively in agreement with
the negative M*/L slopes we inferred dynamically.

In the literature, there are more quantitative estimates of M*/L
gradients in ETGs. Tortora et al. (2011) carried out a stellar popula-
tion synthesis analysis on a sample of ~10 000 galaxies with SDSS
photometry, allowing for gradients in age and metallicity in their
model. They found that massive ETGs show on average positive
M*/L gradients, i.e. a higher M/L in the outskirts with respect to
the central parts, in contrast with our inference, although they find

<figure><img src="image_14.png" /><figcaption>**Figure 13.** Observed frame i-band effective radius as a function of g-band effective radius, from SDSS photometry, for SLACS lenses and the control sample.</figcaption></figure>

negative gradients if they limit the analysis to systems with older
stellar populations.

Szomoru et al. (2013) used an empirical colour-Υ* relation and SDSS photometric data to constrain gradients in M/L on a sample of 220 high mass galaxies with HST data. They quote an average ratio between half-mass radius and half-light radius of R_m/R_e ~ 0.75 for this sample. This value corresponds to an M/L slope of γ ≈ −0.13, as can be seen from Fig. 11.

Newman, Ellis & Treu (2015) estimated M*/L gradients in a sample of ten strong lens ETGs in the centre of massive groups, converting colour gradients into M*/L gradients assuming a fixed age. They found a median M*/L slope of −0.15.

Finally, Poci, Cappellari & McDermid (2017) carried out a study
of the radial profile of M*/L of galaxies in the Atlas 3D sample
(Cappellari et al. 2011), using spatially resolved spectroscopic data.
They found an average M*/L profile very close to flat.

Overall, these studies reveal M/L gradients that are somewhat
shallower than what is required to account for the average slope
we measured from SLACS data, μγ = -0.24 ± 0.04, although our
inference on the slope is strongly degenerate with the IMF mismatch
parameter, as we already pointed out. The measurements listed
above have been obtained with models with a fixed IMF, but steeper
M*/L slopes can be obtained by allowing for a gradient in IMF. If we
allow for the IMF to be not universal, then we should expect IMF
gradients in massive ETGs, given our understanding of their mass
assembly history. The cores of these objects are believed to form
at high redshift, while the outskirts are accreted via mergers with
smaller systems down to present times. Stars in different regions
of massive ETGs are then believed to have formed in different
environments, and therefore could have different IMFs.

Radial gradients in the IMF have been claimed by Martín-Navarro
et al. (2015), La Barbera et al. (2016), and van Dokkum et al. (2017),
based on spatially resolved spectroscopic studies of nearby ETGs.
On M87, a gradient in IMF has been detected both via dynamical
modelling (Oldham & Auger 2018a) and from a spectral analysis
(Sarzi et al. 2018). Collett et al. (2018) measured an IMF gradient
in a nearby lens, through a joint lensing and stellar kinematics
analysis. A similar result was found for a cluster brightest central
galaxy by Smith, Lucey & Edge (2017b), although they allow for
an alternative interpretation consisting of a spatially constant IMF
and a very massive central black hole. At the moment, we cannot
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Top: Enclosed IMF normalization profile, defined as the ratio between the true stellar mass and the stellar population synthesis stellar mass, calculated as a function of projected distance from the centre, for a galaxy with $log M_{*}^{(Chab)} = 11.3$ and $R_{e} = 7$ kpc. The band spans the $1\sigma$ confidence region, as inferred by fitting the control sample. The purple band corresponds to our fiducial model, introduced in Section 4.3, while the green band is inferred by fitting a model with a $M_{*}/L$ gradient and dark matter haloes described by an adiabatically contracted profile with contraction efficiency $v = 0.5$. Bottom: Projected dark matter fraction profile for the two models.</figcaption></figure>

claim nor rule out IMF gradients based on our analysis. In order to determine whether IMF gradients are needed to reproduce our measurements, we need to carry out a spatially resolved stellar population synthesis analysis of the galaxies in our samples. We leave this for a future study.

Regardless of its origin, allowing for a gradient in $M/L$ changes dramatically the average IMF inferred from the combination of lensing and stellar kinematics. While Auger et al. (2010a) and Sonnenfeld et al. (2015) found an average IMF normalization heavier than a Salpeter IMF in their study of SLACS lenses, our current analysis yields a value of $\mu_{IMF,0} = 0.02 \pm 0.04$, consistent with a Chabrier IMF. In parallel, the inferred dark matter masses are revised upwards. Sonnenfeld et al. (2015) measured values of the projected dark matter mass fraction within 5 kpc below 20 per cent, significantly smaller than predictions from numerical simulations (see e.g. fig. 10 of Xu et al. 2017). In the current analysis, the dark matter fraction for the same objects is increased to 30 per cent.

We point out, however, that the value of the IMF normalization we infer is particularly sensitive to two of our model assumptions: the functional form of the $M_{*}/L$ profile, which we assume to be a power law at all radii (equation 27), and the choice of the dark matter density profile, which in our fiducial model is fixed to an NFW model. To gauge the impact of the assumed $M_{*}/L$ gradient, we plot in the top panel of Fig. 14 the inferred enclosed IMF mismatch parameter profile, $\alpha_{IMF}(<R)$, defined as the ratio between the true stellar mass enclosed within projected radius $<R$ and the stellar population synthesis-derived stellar mass enclosed within the same aperture, for an average galaxy with $log M_{*}^{(Chab)} = 11.3$ and $R_{e} = 7$ kpc. For our fiducial model (purple contour), the value of $\alpha_{IMF}(<R)$ converges to the inferred value of $\mu_{IMF,0} = 0.02 \pm 0.04$ at large radii. However, around the scale of the half-light radius, 5 kpc $< R < 10$ kpc, the region our dynamical constraints are most sensitive to, we infer values of $\alpha_{IMF}(<R)$ that are ~0.1 dex higher, roughly in between the values corresponding to a Chabrier and a

Salpeter IMF. Since the stellar mass profile is relatively unconstrained at large radii, we consider the inferred value of $\alpha_{IMF}$ in the inner regions to be more robust, compared to the value of $\mu_{IMF,0}$: an alternative model in which the $\alpha_{IMF}(<R)$ profile remains flat at $R > 10$ kpc would most likely be indistinguishable from our model based on a pure power-law $M_{*}/L$ profile, given our data.

To estimate how much our results depend on the assumed dark matter density profile, we repeat the analysis assuming adiabatically contracted dark matter halos, following the model introduced in Section 4.2, with a value of the contraction efficiency parameter fixed to $v = 0.5$. We pick this value because it is close to the average contraction efficiency we inferred when fitting the control sample with the model with no $M_{*}/L$ gradient. The main effect of assuming a dark matter model with a steeper inner density profile is a change in the value of the average IMF normalization, which for the control sample is inferred to be $\mu_{IMF,0} = -0.06 \pm 0.04$. This change is illustrated in the top panel of Fig. 14, where the enclosed $\alpha_{IMF}$ profile for this new model is plotted in green. We point out how, while the average IMF normalization goes down in a model with adiabatically contracted dark matter halos, the inferred value of the $M_{*}/L$ gradient does not change significantly. The need for a gradient in the model appears to be robust to the particular choice of the dark matter density profile.

In the bottom panel of Fig. 14, we plot the inferred projected dark matter fraction profile for the two models. This quantity is, quite naturally, particularly sensitive to the assumed dark matter density profile.

Relaxing the assumption of a spatially constant $M_{*}/L$ is the key ingredient for reconciling strong lensing observations, which probe the mass distribution at the scales of a few kpc, with weak lensing constraints, probing scales in the range 10–300 kpc. In principle, there are other model assumptions that, if relaxed, could help mitigate the tension between the SLACS and HSC measurements. As is often done in joint strong lensing and stellar dynamic studies, we assumed spherical symmetry for the gravitational potential and isotropic orbits for the stars. If we allow, for example, for radial anisotropy, the same mass distribution would produce a higher line-of-sight velocity dispersion. As pointed out earlier, failure to match the velocity dispersion of some SLACS lenses is the reason why fits with our vanilla and adiabatically contracted/expanded model return very low halo masses. Therefore, in principle, allowing for radial anisotropy could help bring halo masses inferred from strong lensing closer to the values inferred from weak lensing.

Dynamical models that allow for anisotropy have been fitted to SLACS data by Posacki et al. (2015). Even by allowing for anisotropy, the inferred dark matter fractions are very close to zero for a significant number of objects, suggesting that anisotropy alone is unlikely to be the solution to the apparent tension between SLACS and HSC measurements observed in the context of the vanilla model.

Another effect that could increase the predicted velocity dispersions is allowing for elongation in the mass distribution along the line of sight (see e.g. fig. 17 of Sonnenfeld et al. 2012). Deviations from spherical symmetry of SLACS lenses have been explored by Barnabè et al. (2011) using spatially resolved stellar kinematics measurements. They found that the inferred density profiles are consistent with those derived from a simple spherical Jeans equation approach.

The presence of a central supermassive black hole, not included in our fiducial model, would also boost the central velocity dispersion, increasing the mass enclosed in the central regions in a similar way as a radially decreasing IMF normalization does (see e.g. Smith, Lucey & Edge 2017a).
To test the effects of supermassive black holes on our measurement, we consider a ‘vanilla + black hole’ model, consisting of NFW dark matter halos, de Vaucouleurs stellar profiles with spatially constant $M_*/L$ and a central black hole with mass $M_{BH}$. Fitting individual SLACS lenses with this model, assuming uninformative (flat) priors on the model parameters, can lead to very large, unphysical, values of the ratio between black hole mass and stellar mass, especially for the lenses with steep density slopes ($\gamma' > 2.3$). We then put a prior on $M_{BH}$ assuming a power-law relation between black hole mass and the central velocity dispersion, with Gaussian scatter in $\log M_{BH}$ at fixed velocity dispersion, as measured by McConnell & Ma (2013). We also set an upper limit on the ratio between black hole mass and stellar mass: $M_{BH}/M_* < 0.01$.

We fit the SLACS sample with this model, using the same form for the distribution of halo mass, stellar mass and IMF normalization as that used to describe the vanilla model. We obtain an average log-halo mass of $\mu_{h,0} = 12.51 \pm 0.19$. This is slightly larger than the value obtained with the vanilla model, but still inconsistent with the halo mass constraints from weak lensing. We then conclude that accounting for the presence of central black holes does not change significantly our inference.

Finally, some of our results could, in principle, change if we drop the assumption of a de Vaucouleurs profile for the light distribution. For instance, adopting the more general Sérsic model (Sérsic 1963) could lead to overall steeper density profiles that would in part mimic the effects of a $M/L$ gradient. Shankar et al. (2018) showed how an evolving Sérsic index with time can help match simple semi-empirical models with lensing observations of the redshift evolution of the density slope $\gamma'$. In order to robustly assess the impact of the use of a Sérsic model, we need to obtain self-consistent stellar mass, Sérsic indices, and effective radii for both the SLACS lenses and the control sample. We leave this exploration for future work. Nevertheless, the more flexible analysis of Posacki et al. (2015), who modelled the surface brightness distribution of SLACS lenses with a multi-Gaussian expansion, produced similar results to our study, suggesting that the effect of deviations from a pure de Vaucouleurs profile are probably small.

We point out how very low inferred dark matter fractions are not a peculiarity of our study of SLACS lenses, but appear to be a common occurrence in stellar dynamics measurements on ETGs, even when more sophisticated modelling tools compared to our isotropic spherical Jeans analysis are used. This is the case, for example, for the analysis of Atlas 3D galaxies by Cappellari et al. (2013). Cappellari et al. (2013) used Jeans axisymmetric dynamical models with orbital anisotropy to fit stellar kinematics data from integral field spectroscopy. They found a median dark matter fraction within the 3D half-light radius as low as 13 per cent, indicating that low inferred dark matter masses are robust to the inclusion of anisotropy or deviations from spherical symmetry. The Cappellari et al. (2013), however, assumed a spatially constant $M/L$ for the stellar component. As pointed out by Bernardi et al. (2018), allowing for a gradient would change significantly the inferred stellar and dark matter masses of the Atlas 3D sample of galaxies, as is the case for our analysis of SLACS lenses. Recently, Pechetti et al. (2017) re-analysed a subsample of 27 low-mass ($\log M_* < 10.5$) galaxies from the Atlas 3D sample, using high resolution imaging from HST. They found evidence for an increase in the stellar $M/L$ in the inner regions of these galaxies, which they interpret as the combined effect of the presence of over-massive central black holes and IMF gradients.

<figure><img src="image_16.png" /><figcaption>**Figure 15.** Stellar mass versus effective radius and velocity dispersion of mock galaxies. **Grey circles:** Sample of 1000 mock galaxies generated from the maximum-likelihood point of the $M_*/L$ gradient model, fitted to the control sample. **Blue circles:** Subsample of 100 objects selected with a random draw weighted by the strong lensing cross-section of each galaxy. **Dashed lines:** Observed best-fitting $M_*(\text{Chab}) - \sigma_{\text{ap}}$ (top) and $M_*(\text{Chab}) - R_e$ relation (bottom) for the SLACS and the control sample, as plotted in Fig. 2. **Solid lines:** Best-fitting $M_*(\text{Chab}) - \sigma_{\text{ap}}$ (top) and $M_*(\text{Chab}) - R_e$ relation (bottom) inferred from the posterior predictive sample, with (blue) or without (black) the strong lensing selection correction. The $M_*(\text{Chab}) - R_e$ relation of the mock sample is omitted from the bottom panel, since it is set equal to the observed relation by construction.</figcaption></figure>

## 5.1 Differences between strong lenses and non-lenses

A fundamental feature of our analysis is the treatment of strong lensing selection, needed to compare the inference based on the SLACS sample of lenses with the control sample, and described in Section 3.6. Regardless of the model used, a common result of strong lensing selection is to increase the average stellar mass of the sample. Lensing selection has instead a very limited effect on the distribution of halo masses at fixed stellar mass. As discussed in Section 4.1, the reason for this is the low sensitivity of the lensing cross-section to the halo mass. The lensing cross-section is determined by the projected mass within a few kpc, a region dominated by the baryons. As a result, the main effect of strong lensing selection is to modify the distribution in stellar mass.

The change in IMF normalization due to lensing selection is also quite small, in our fiducial model with $M/L$ gradients. This is because the inferred value of the intrinsic scatter in the IMF normalization is small, $\sigma_{\text{IMF}} = 0.06 \pm 0.01$. With a small intrinsic scatter, the effect of lensing selection is also small.

The mass-size and mass-velocity dispersion relation are also relatively unchanged by strong lensing selection. This is shown in Fig. 15, where we plot the $M_*(\text{Chab}) - R_e$ and $M_*(\text{Chab}) - \sigma$ relation of the mock sample, generated from the fiducial model fitted to the control sample, together with the distribution of a sample of lenses, drawn from the same mock sample with a probability proportional to the strong lensing cross-section of each galaxy. Strong lenses are only slightly more compact, ~0.02 dex on average, at fixed stellar mass, with respect to the general population, and have higher ve-
<figure><img src="image_17.png" /><figcaption>**Figure 16.** Shaded regions: Predicted average line-of-sight velocity dispersion profile, calculated for two galaxy with log $M_*^{(Chab)}$ = 11.3, 11.6, $R_e$ = 7, 10 kpc, with average values of the IMF mismatch parameter and $M_*/L$ gradient slope, and halo mass corresponding to the 16, 50, and 84 percentile of the distribution, as inferred from the control sample. The bands delimit 68 per cent enclosed probability regions. The profiles are calculated using the spherical Jeans equation under the assumption of isotropic orbits. Grey lines: velocity dispersion profiles measured by Veale et al. (2018) for a set of ETGs from the MASSIVE survey, with $M_K > -25.6$.</figcaption></figure>

locity dispersion. This is in qualitative agreement with observations of the SLACS sample, although we are unable to reproduce in detail the observed relations, particularly the smaller average sizes of SLACS lenses observed at the low-mass end. However, the difference between our prediction and the observed $M_*^{(Chab)} - R_e$ relation is of the same magnitude as the 8 per cent systematic uncertainty on our measurements of $R_e$, and therefore not very significant.

## 5.2 Model prediction: velocity dispersion profile

The $M_*/L$ gradient model provides a good fit to the data used in this study: Einstein radii, central velocity dispersions and, on larger scales, weak lensing. In principle, we could add more constraints to the analysis, in order to test even more complex models than the current one. One possibility would be to make use of spatially resolved kinematics measurements, which could, for example, allow us to relax the assumption of isotropic orbits. Although there are integral field spectroscopy observations available for 17 of the SLACS lenses (Czoske et al. 2012), their analysis would require a significant amount of work, and we leave it for future study. Nevertheless, we can use our model to predict the velocity dispersion profile of a typical galaxy, and make a qualitative comparison with existing measurements.

In Fig. 16, we plot the line-of-sight velocity dispersion profile of two galaxies with log $M_*^{(Chab)}$ = 11.3, 11.6, $R_e$ = 7, 10 kpc, average IMF normalization and $M_*/L$ gradient slope for their stellar mass, and halo mass corresponding to the 16, 50, and 84 percentile of the inferred distribution. We assume orbital isotropy. The predicted profile shows a relatively sharp drop in the line-of-sight velocity dispersion, followed by either a slowly falling or rising outer profile, depending on the halo mass. The integral field spectroscopy data of SLACS lenses show a qualitatively similar behaviour in the inner regions (see fig. 5 of Czoske et al. 2012), but lack observations at radii larger than $R_e$, where the profile is predicted to flatten.

We can then compare our prediction with observations of $z \sim 0$ ETGs. Veale et al. (2018) studied in detail the velocity dispersion profile of 90 ETGs from the MASSIVE survey (Ma et al. 2014). Although the average mass of the Veale et al. (2018) sample is somewhat larger than that of our samples, there is significant overlap. We over-plot in Fig. 16 the velocity dispersion profile of a subset of MASSIVE galaxies, selected by setting a limit on their K-band magnitude, $M_K > -25.6$. These are the objects plotted in the left-hand panel of fig. 1 of Veale et al. (2018). Most of the MASSIVE galaxies exhibit a velocity dispersion profile that is falling with radius in the inner regions and then rising or staying constant at larger radii, in qualitative agreement with our prediction.

However, the Veale et al. (2018) reveals a much wider variety of velocity dispersion profiles compared to what our model is able to describe: there are objects with flat inner profile, and objects with steeply rising or falling outer profile, which are difficult to reproduce in the context of our model. This seems to indicate that the real structure of massive galaxies is most likely more complex than our description. The recent strong lensing study of Oldham & Auger (2018b) appears to corroborate this idea. Oldham & Auger (2018b) were able to disentangle the luminous and dark matter on individual objects for a set of 12 lenses. Even by allowing for the presence of $M_*/L$ gradients, which they detect in a quarter of their objects, they find large variations in the inner dark matter density slope across the population, ranging from cored density profiles to profiles as steep as isothermal.

It is possible that by relaxing our assumption of a universal dark matter density profile we can produce a broader range of velocity dispersion profiles, in better agreement with the observations from the MASSIVE survey. For a quantitative comparison between our model and spatially resolved kinematics data, however, additional ingredients, such as orbital anisotropy and rotation, must be taken into account.

# 6 CONCLUSIONS

We analysed two samples of massive quiescent galaxies, selected in stellar mass by applying a cut log $M_*^{(Chab)}$ > 11.0. The first sample consists of ~1700 galaxies in the SDSS legacy spectroscopic sample, with weak lensing measurements from HSC. The second sample consists of 45 strong lenses from the SLACS survey. We fitted models for the distribution of stellar and dark matter mass and stellar IMF across the population of galaxies to weak lensing, strong lensing, and stellar kinematics data on the two samples. We then compared the results obtained on the two samples, after correcting for strong lensing selection effects, with the goal of finding a simple model for structure of massive galaxies that can reproduce observations on both data sets. We found the following:

(i) A model in which dark matter halos have an NFW density profile and the stellar mass follows the light distribution does not provide a good description of SLACS lenses: the inferred halo masses are too low compared to galaxies of the same stellar mass from the control sample, a discrepancy that cannot be ascribed to strong lensing selection.

(ii) Modifying the dark matter profile by allowing for adiabatic contraction or expansion does not solve the discrepancy between the two samples: SLACS lenses want strong adiabatic expansion while HSC galaxies prefer mild contraction. This result is driven by lenses with a particularly steep total density profile, for which even mass-follows-light models struggle to fit both the strong lensing and stellar kinematics constraints.
(iii) Allowing for a gradient in stellar M/L, results in a good match between the SLACS sample and HSC galaxies, and shifts downwards the inferred average IMF normalization, the exact value depending on the assumed dark matter density profile.

(iv) In light of the current analysis, such a gradient can be interpreted both in terms of a gradient in stellar population properties at fixed IMF, or in terms of a gradient in IMF itself.

(v) The main effect of strong lensing selection is to shift the median stellar mass towards higher values. Galaxy properties at fixed stellar mass show little differences between strong lens samples and the general galaxy population.

The main limitation of our study is the relatively low number of constraints for individual galaxies, which forces us to adopt rather simple mass models. Another limitation is the necessity to correct for strong lensing selection effects when mapping the HSC sample, used for the bulk of our weak lensing analysis, to the SLACS sample of strong lenses. The need for the use of two different samples of galaxies is dictated by the low signal-to-noise ratio of the weak lensing data available for the strong lens sample. This situation can be avoided by using a much larger sample of strong lenses, so that both the strong and the weak lensing analysis can be carried out on the same set of objects. Current surveys such as HSC, the Kilo Degree Survey (de Jong et al. 2015) and the Dark Energy Survey (Diehl et al. 2014), are enabling the discovery of hundreds of new lenses (Diehl et al. 2017; Petrillo et al. 2017; Sonnenfeld et al. 2018), and in the future, it will be possible to perform a joint strong and weak lensing analysis on a large sample of galaxies.

## ACKNOWLEDGEMENTS

We thank Melanie Veale for kindly providing observed velocity dispersion profiles of MASSIVE galaxies. This work was supported by World Premier International Research Center Initiative (WPI Initiative), Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan. AS is partly supported by KAKENHI Grant Number JP17K14250.

## REFERENCES

* Aihara H. et al., 2018, PASJ, 70, S4

* Alam S., et al., 2015, ApJS, 219, 12

* Arneson R. A., Brownstein J. R., Bolton A. S., 2012, ApJ, 753, 4

* Auger M. W., Treu T., Bolton A. S., Gavazzi R., Koopmans L. V. E., Marshall P. J., Bundy K., Moustakas L. A., 2009, ApJ, 705, 1099

* Auger M. W., Treu T., Bolton A. S., Gavazzi R., Koopmans L. V. E., Marshall P. J., Moustakas L. A., Burles S., 2010b, ApJ, 724, 511

* Auger M. W., Treu T., Gavazzi R., Bolton A. S., Koopmans L. V. E., Marshall P. J., 2010a, ApJ, 721, L163

* Barnabè M., Czoske O., Koopmans L. V. E., Treu T., Bolton A. S., 2011, MNRAS, 415, 2215

* Bernardi M., Sheth R. K., Dominguez-Sanchez H., Fischer J.-L., Chae K.-H., Huertas-Company M., Shankar F., 2018, MNRAS, 477, 2560

* Blumenthal G. R., Faber S. M., Flores R., Primack J. R., 1986, ApJ, 301, 27

* Bolton A. S., Burles S., Koopmans L. V. E., Treu T., Moustakas L. A., 2006, ApJ, 638, 703

* Cappellari M. et al., 2011, MNRAS, 413, 813

* Cappellari M. et al., 2012, Nature, 484, 485

* Cappellari M. et al., 2013, MNRAS, 432, 1709

* Collett T. E. et al., 2018, Science, 360, 1342

* Czoske O., Barnabè M., Koopmans L. V. E., Treu T., Bolton A. S., 2012, MNRAS, 419, 656

* de Jong J. T. A. et al., 2015, A&A, 582, A62

* de Vaucouleurs G., 1948, Ann. Astrophys., 11, 247

* Diehl H. T. et al., 2014, in *Observatory Operations: Strategies, Processes, and Systems V. Conference proceeding of SPIE*, the International society for Optics and Photonics, p. 91490V

* Diehl H. T. et al., 2017, ApJS, 232, 15

* Dutton A. A., Macciò A. V., Mendel J. T., Simard L., 2013, MNRAS, 432, 2496

* Dutton A. A., van den Bosch F. C., Dekel A., Courteau S., 2007, ApJ, 654, 27

* Ferreras I., La Barbera F., de la Rosa I. G., Vazdekis A., de Carvalho R. R., Falcón-Barroso J., Ricciardelli E., 2013, MNRAS, 429, L15

* Furusawa H., Koike M., Takata T., et al., 2018, PASJ, 70, S3

* Gavazzi R., Treu T., Rhodes J. D., Koopmans L. V. E., Bolton A. S., Burles S., Massey R. J., Moustakas L. A., 2007, ApJ, 667, 176

* Hirata C., Seljak U., 2003, MNRAS, 343, 459

* Kauffmann G. et al., 2003, MNRAS, 341, 33

* Kawanomoto S., Miyazaki S., Komiyama Y. Y. X., 2018, PASJ, 70, S66

* Komiyama Y., Obuchi Y., Nakaya H., et al., 2018, PASJ, 70, S2

* Koopmans L. V. E., Treu T., Bolton A. S., Burles S., Moustakas L. A., 2006, ApJ, 649, 599

* La Barbera F., Vazdekis A., Ferreras I., Pasquali A., Cappellari M., Martín-Navarro I., Schönebeck F., Falcón-Barroso J., 2016, MNRAS, 457, 1468

* Ma C.-P., Greene J. E., McConnell N., Janish R., Blakeslee J. P., Thomas J., Murphy J. D., 2014, ApJ, 795, 158

* Macciò A. V., Dutton A. A., van den Bosch F. C., 2008, MNRAS, 391, 1940

* Mandelbaum R. et al., 2014, ApJS, 212, 5

* Mandelbaum R. et al., 2015, MNRAS, 450, 2963

* Mandelbaum R. et al., 2018, PASJ, 70, S25

* Martín-Navarro I., Barbera F. L., Vazdekis A., Falcón-Barroso J., Ferreras I., 2015, MNRAS, 447, 1033

* McConnell N. J., Ma C.-P., 2013, ApJ, 764, 184

* Miyazaki S., Komiyama Y., Kawanomoto S., et al., 2018, PASJ, 70, S1

* Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493

* Newman A. B., Ellis R. S., Treu T., 2015, ApJ, 814, 26

* Oldham L., Auger M., 2018a, MNRAS, 474, 4169

* Oldham L. J., Auger M. W., 2018b, MNRAS, 476, 133

* Pechetti R., Seth A., Cappellari M., McDermid R., den Brok M., Mieske S., Strader J., 2017, ApJ, 850, 15

* Petrillo C. E. et al., 2017, MNRAS, 472, 1129

* Planck Collaboration XIII, 2016, A&A, 594, A13

* Poci A., Cappellari M., McDermid R. M., 2017, MNRAS, 467, 1397

* Posacki S., Cappellari M., Treu T., Pellegrini S., Ciotti L., 2015, MNRAS, 446, 493

* Ruff A. J., Gavazzi R., Marshall P. J., Treu T., Auger M. W., Brault F., 2011, ApJ, 727, 96

* Rykoff E. S. et al., 2014, ApJ, 785, 104

* San Roman I. et al., 2018, A&A, 609, A20

* Sarzi M., Spiniello C., La Barbera F., Krajnović D., van den Bosch R., 2018, MNRAS, 478, 4084

* Schneider M. D., Hogg D. W., Marshall P. J., Dawson W. A., Meyers J., Bard D. J., Lang D., 2015, ApJ, 807, 87

* Schulz A. E., Mandelbaum R., Padmanabhan N., 2010, MNRAS, 408, 1463

* Shankar F. et al., 2017, ApJ, 840, 34

* Shankar F. et al., 2018, MNRAS, 475, 2878

* Smith R. J., Lucey J. R., Edge A. C., 2017a, MNRAS, 467, 836

* Smith R. J., Lucey J. R., Edge A. C., 2017b, MNRAS, 471, 383

* Sonnenfeld A., Leauthaud A., 2018, MNRAS, 477, 5460

* Sonnenfeld A., Treu T., Gavazzi R., Marshall P. J., Auger M. W., Suyu S. H., Koopmans L. V. E., Bolton A. S., 2012, ApJ, 752, 163

* Sonnenfeld A., Treu T., Gavazzi R., Suyu S. H., Marshall P. J., Auger M. W., Nipoti C., 2013, ApJ, 777, 98

* Sonnenfeld A., Treu T., Marshall P. J., Suyu S. H., Gavazzi R., Auger M. W., Nipoti C., 2015, ApJ, 800, 94

* Sonnenfeld A. et al., 2018, PASJ, 70, S29

* Spiniello C., Trager S., Koopmans L. V. E., Conroy C., 2014, MNRAS, 438, 1483

* Szomoru D., Franx M., van Dokkum P. G., Trenti M., Illingworth G. D., Labbé I., Oesch P., 2013, ApJ, 763, 73
Sérsic J. L., 1963, BAAA, 6, 41

Tanaka M., 2015, ApJ, 801, 20

Tanaka M. et al., 2018, PASJ, 70, S9

Tortora C., Napolitano N. R., Romanowsky A. J., Jetzer P., Cardone V. F., Capaccioli M., 2011, MNRAS, 418, 1557

Treu T., Auger M. W., Koopmans L. V. E., Gavazzi R., Marshall P. J., Bolton A. S., 2010, ApJ, 709, 1195

van Dokkum P., Conroy C., Villaume A., Brodie J., Romanowsky A. J., 2017, ApJ, 841, 68

van Dokkum P. G., Conroy C., 2010, Nature, 468, 940

Veale M., Ma C.-P., Greene J. E., Thomas J., Blakeslee J. P., Walsh J. L., Ito J., 2018, MNRAS, 473, 5446

Xu D., Springel V., Sluse D., Schneider P., Sonnenfeld A., Nelson D., Vogelsberger M., Hernquist L., 2017, MNRAS, 469, 1824

This paper has been typeset from a *TEX/LATEX* file prepared by the author.