© The Authors 2025

1. Introduction

Strong gravitational lensing is one of the most reliable methods for studying the mass distribution of galaxies. It can provide a measurement of the projected mass of a lens galaxy with a precision and accuracy of a few percent, regardless of the dynamical state of the lens. The best-studied sample of lenses so far is the Sloan Lens ACS (SLACS) survey (Bolton et al. 2006), consisting of about one hundred strong lenses, most of which are massive early-type galaxies. The standard strategy for studying the inner structure of the SLACS lenses so far has been to combine information from strong lensing and stellar kinematics. Joint lensing and stellar dynamics analyses have revealed that SLACS lenses have, on average, a total density profile slightly steeper than isothermal (Koopmans et al. 2006; Auger et al. 2010a; Barnabè et al. 2011) and that the steepness of the density profile correlates positively with the stellar mass density (Sonnenfeld et al. 2013). Lensing and dynamics have also been used to determine the relative contribution of stars and dark matter to the total mass budget of massive galaxies (Treu et al. 2010; Auger et al. 2010b; Sonnenfeld et al. 2015; Posacki et al. 2015; Shajib et al. 2021; Sheu et al. 2024). According to a large part of previous investigations, the stars appear to have a relatively large mass-to-light ratio, similar to that predicted by a Salpeter stellar initial mass function (IMF), while the dark matter mass enclosed within the inner regions seems to be smaller than that predicted by numerical simulations (see e.g. Fig. 11 of Shajib et al. 2021).

Interpreting the stellar kinematics data of a galaxy is difficult. It requires modelling its three-dimensional mass structure and the orbits of its stars, which are two properties that are not directly observable. While some of the degeneracies affecting dynamical models can be broken with the help of integral field spectroscopic data, for most of the SLACS lenses, the only stellar kinematics information currently available comes from Sloan Digital Sky Survey (SDSS) single-aperture spectra. This scarcity of spatially resolved spectroscopic data has forced past investigations to rely on simplifying assumptions. Nearly all of the studies mentioned above have assumed spherical symmetry, and some of them have asserted orbital isotropy (Auger et al. 2010a; Sonnenfeld et al. 2015). One notable exception is the work of Barnabè et al. (2011), who fitted models with axial symmetry and orbits described by a two-integral distribution function to spatially resolved kinematics data of 16 lenses. Barnabè et al. (2011) modelled the radial density profile of the lenses as a power law, ρ(r)∝r^−γ, and measured the density slope γ for each of the lenses. Their values of γ are systematically lower, by approximately 0.05, than those obtained by Auger et al. (2010a) with an isotropic spherical Jeans model. Such a discrepancy is larger than the current uncertainty on the average γ of the SLACS lens sample, and should therefore warn against assuming spherical symmetry or isotropy when studying the SLACS lenses.

The problem of analysing SLACS stellar kinematics data is made worse by selection effects. SLACS lenses were selected on the basis of the observed velocity dispersion (Bolton et al. 2006). This caused lenses with a larger observed velocity dispersion to be overrepresented, with respect to a purely mass-selected sample (Sonnenfeld 2024). In order for a dynamical analysis to accurately capture the mass and orbital structure of a SLACS lens, these selection effects must be taken into account. So far, this has only been done by Sonnenfeld (2024, hereafter Paper I), albeit in the simplified context of spherical and isotropic models. Since the velocity dispersion correlates with the three-dimensional structure and orbital anisotropy of a galaxy, a more complex and realistic dynamical model of a SLACS lens would require a non-trivial prior on these quantities. For instance, models with a larger radial anisotropy and elongation along the line of sight would have to be upweighted.

In this paper, we re-analyse the SLACS sample with a more conservative approach. We focus exclusively on projected quantities, which are directly observable, and use only gravitational lensing data as constraints on the mass distribution. We combine strong lensing measurements of SLACS with weak lensing observations obtained for the parent galaxy sample from which the lenses were drawn while explicitly accounting for selection effects. The primary goal of this work is to put constraints on the stellar mass-to-light ratio of massive galaxies and on their inner dark matter density profile. We want to determine to what extent gravitational lensing data, on its own, can support previous claims of a non-universal IMF in early-type galaxies and whether it provides evidence for the contraction of dark matter, a phenomenon ubiquitously observed in hydrodynamical simulations at the mass scales of the SLACS lenses (Gnedin et al. 2004; Duffy et al. 2010; Schaller et al. 2015; Cautun et al. 2020) but for which there is still little experimental evidence.

Our secondary goal is to better understand how the velocity dispersion of the SLACS lenses differs from that of massive galaxies of the same mass, and to corroborate the claim made in Paper I that the values of the observed velocity dispersion are overestimated. This is an important point because the SLACS lenses, with their stellar kinematics data, are being used in support of time delay measurements of the Hubble constant (Birrer et al. 2020), as a source of information on the internal structure of strong lenses. If the SLACS lenses are biased in a unique way with respect to the general population of galaxies, they might also be biased with respect to the time-delay lenses, and therefore bias the measurement of the Hubble constant.

This paper is structured as follows. In Section 2 we describe the methods and data used for the analysis. In Section 3 we show the results. We discuss and conclude in Section 4. Throughout this work, we assume a flat Λ cold dark matter model, with H₀ = 70 km s⁻¹ Mpc⁻¹ and Ω_m = 0.3. All masses are in solar units, lengths are in units of kiloparsecs, and velocities are in kilometres per second.

2. Methods

2.1. Analysis formalism

Our analysis method follows in large part that of Paper I. The starting point is a simplified version of the fundamental equation of statistical strong lensing,

$$\begin{array}{c}\hfill {\mathrm{P}}_{\mathrm{SL}}^{(\mathrm{eff})}({\mathit{\psi }}_{\mathrm{g}},z_{\mathrm{s}})\propto {\mathrm{P}}_{\mathrm{g}}({\mathit{\psi }}_{\mathrm{g}}|\mathit{\eta }){\mathrm{P}}_{\mathrm{s}}^{(\mathrm{eff})}(z_{\mathrm{s}}|\mathit{\eta }){\mathrm{P}}_{\mathrm{sel}}({\mathit{\psi }}_{\mathrm{g}},z_{\mathrm{s}}|\mathit{\eta }),\end{array}$$(1)

which describes the probability distribution of lens-source pairs (left-hand side) as the product of the distribution of foreground galaxies P_g, the empirically derived effective distribution of background sources P_s^(eff), and the lens selection probability P_sel. In the equation above, ψ_g indicates the set of parameters describing a foreground galaxy, z_s is the source redshift, while η is the ensemble of population-level parameters, also referred to as hyper-parameters of the model. Equation (1) was derived in Paper I by means of a key approximation: the probability of a galaxy-source pair to give rise to a detectable strong lens separates cleanly into a factor that depends purely on geometry and another factor that depends purely on the source brightness. That approximation allowed us to eliminate any explicit dependence on the source brightness from Equation (1).

The goal of the analysis is to infer the population-level parameters η given the data d of the SLACS lenses. The posterior probability of η is

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

and the likelihood is the product over all of the SLACS lenses:

$$\begin{array}{c}\hfill \mathrm{P}(\mathbf{d}|\mathit{\eta })={\displaystyle \prod _i}\mathrm{P}({\mathbf{d}}_{\mathbf{i}}|\mathit{\eta }).\end{array}$$(3)

To evaluate each factor, it is necessary to marginalise over all possible values of the individual lens parameters and source redshift, the prior of which is given by the distribution of Equation (1):

$$\begin{array}{c}\hfill \mathrm{P}({\mathbf{d}}_{\mathbf{i}}|\mathit{\eta })={\displaystyle \int d}{\mathit{\psi }}_{\mathrm{g},i}d{z}_{\mathrm{s}}i\mathrm{P}({\mathbf{d}}_{\mathbf{i}}|{\mathit{\psi }}_{\mathrm{g},i},z_{\mathrm{s}}i,\mathit{\eta }){\mathrm{P}}_{\mathrm{SL}}^{(\mathrm{eff})}({\mathit{\psi }}_{\mathrm{g},i},z_{\mathrm{s}}i|\mathit{\eta }).\end{array}$$(4)

2.2. Data

The sample and data used are the same as in Paper I. These are 59 early-type galaxy strong lenses from the SLACS sample, each with the following measurements:

• The lens and source spectroscopic redshift, z_g and z_s.

• The stellar population synthesis-based stellar mass, M_*^(obs), based on a Chabrier IMF and a de Vaucouleurs stellar profile, and its uncertainty.

• The half-light radius R_e.

• The Einstein radius, ${\theta }_{\mathrm{E}}^{(\mathrm{obs})}$, with a conservative 5% uncertainty.

• The line-of-sight stellar velocity dispersion from SDSS, σ_ap^(obs), and its uncertainty.

These data were taken from Auger et al. (2009). Uncertainties on redshifts and half-light radii are small and were ignored in the analysis.

2.3. Mass model

We adopted a two-component mass model, consisting of a stellar bulge and a dark matter halo. The stars follow a de Vaucouleurs profile, described by two parameters: the total stellar mass, M_*, and the half-light radius, R_e. Although the surface brightness profile of the SLACS lenses is better described by more complex models, such as the sum of multiple Sérsic components (Etherington et al. 2022), a simple two-parameter model such as a de Vaucouleurs profile can already provide a good description of the inner stellar distribution of early-type galaxies, with typical errors in the inner 10 kpc on the order of a few percent (Sonnenfeld 2020). In addition to M_*, we considered the stellar population synthesis stellar mass, M_*^(sps), which we defined as the stellar mass that an observer would infer by fitting the stellar population synthesis model of Auger et al. (2009) to perfect photometric data. The true stellar mass M_* is related to M_*^(sps) by means of the stellar population synthesis mismatch parameter, α_sps, defined as

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{sps}}=\frac{M_{\ast }}{M_{\ast }^{(\mathrm{sps})}}·\end{array}$$(5)

This means that, if the Auger et al. (2009) stellar population synthesis model is accurate, then α_sps = 1. The stellar population mismatch parameter is one of the main quantities that we seek to constrain in this work.

To describe the dark matter density profile, we emulated the effect of the gravitational contraction resulting from the infall of baryons. The idea is that the dark matter and the baryonic matter that eventually formed stars followed initially the same distribution. Subsequently, the baryons sank to the centre due to dissipative processes, altering the gravitational potential, and the dark matter particles responded by shrinking their orbits. We modelled this process under the approximation of adiabatic contraction of particles in circular orbits, as originally suggested by Blumenthal et al. (1986). Although neither the contraction of dark matter is truly adiabatic nor the orbits are circular, this is a widely used approach in semi-analytic models of dark matter halos. It produces results that are qualitatively similar to those of hydrodynamical simulations (see e.g. Cautun et al. 2020), and has been used successfully in lensing and dynamical analyses of the dark matter distribution in galaxies (Auger et al. 2010b; Dutton et al. 2011; Oguri et al. 2014; Shajib et al. 2021).

For a collisionless dark matter particle in a circular orbit at radius r, if the total mass M(r) enclosed within the sphere with radius r varies slowly with respect to the orbital time, then the quantity rM(r) is conserved. If r_i and r_f are respectively the initial (pre-contraction) and final radius of a dark matter particle, then the conservation of rM(r) implies that

$$\begin{array}{c}\hfill r_i{M}_{\mathrm{DM},i}(r_i)(1+\frac{M_{\ast }}{M_{\mathrm{h}}})=r_f[M_{\mathrm{DM},f}(r_f)+M_{\ast }(r_f)].\end{array}$$(6)

In the above equation, M_h is the total halo mass, while M_DM, i(r) and M_DM, f indicate the initial and final dark matter density profiles. Shells at different r do not cross each other, so M_DM, i(r_i) = M_DM, f(r_f), and Equation (6) becomes

$$\begin{array}{c}\hfill r_i{M}_{\mathrm{DM},i}(r_i)(1+\frac{M_{\ast }}{M_{\mathrm{h}}})=r_f[M_{\mathrm{DM},i}(r_i)+M_{\ast }(r_f)].\end{array}$$(7)

This equation describes a dark matter halo that responds solely to the baryons that end up forming the central galaxy. In reality, star formation is not 100% efficient, and part of the baryons involved in the star formation process are ejected by winds, supernovae or active galactic nuclei. Any gas that is ejected from the inner regions would cause the dark matter halo to expand, rather than contract. Consequently, we expect the model described above to be an upper limit to the amount of contraction. In order to allow for a range of contraction scenarios, we modified Equation (7) by introducing an adiabatic contraction efficiency parameter, ϵ:

$$\begin{array}{cc}\hfill r_i{M}_{\mathrm{DM},i}(r_i)& (1+\frac{M_{\ast }}{M_{\mathrm{h}}})\hfill \\ \hfill & =r_f[M_{\mathrm{DM},i}(r_i)(1+(1-ϵ)\frac{M_{\ast }}{M_{\mathrm{h}}})+ϵM_{\ast }(r_f)].\hfill \end{array}$$(8)

In words, Equation (8) is equivalent to a model in which dark matter responds to the contraction of only a fraction ϵ of the stellar mass of the galaxy. This parameterisation of the contraction efficiency is an alternative to the widely adopted one of Dutton et al. (2007), which is based on first obtaining the contraction factor r_f/r_i of each shell from Equation (7), and then rescaling it by taking its power (r_f/r_i)^ν. Our model is easier to implement numerically and, we believe, more intuitive.

Given the initial dark matter distribution and the final stellar mass distribution, Equation (6) allows us to determine the final dark matter profile, by solving for r_i at any value of r_f (or vice-versa). We allowed the value of ϵ to vary continuously from ϵ = 1, which corresponds to the maximally contracted case, to ϵ = 0, for which the solution to Equation (8) is simply r_f = r_i. We modelled the initial dark matter distribution with a Navarro, Frenk and White (NFW) profile. We set M₂₀₀ = M_h, where M₂₀₀ is the mass within a sphere with an average density equal to 200 times the critical density of the Universe, and adopted the following mass-concentration relation from Dutton & Macciò (2014) to determine the initial scale radius:

$$\begin{array}{c}\hfill log{c}_{200}=0.905-0.101(log{M}_{200}-12+logh),\end{array}$$(9)

with h = 0.7. We neglected scatter in concentration at fixed halo mass, to keep the dimensionality of the problem low.

The dark matter profile that is obtained by solving Equation (8) is a non-analytic function, which makes it cumbersome to use to predict lensing-related quantities. In order to speed up computations, we approximated the resulting profile with a generalised Navarro, Frenk and White (gNFW) model:

$$\begin{array}{c}\hfill \rho (r)\propto \frac{1}{r^{{\gamma }_{\mathrm{DM}}}{(1+\frac{r}{r_s})}^{3-{\gamma }_{\mathrm{DM}}}}·\end{array}$$(10)

A gNFW profile has three degrees of freedom: the mass normalisation, the scale radius r_s and the inner density slope γ_DM. To determine these quantities, we fixed the virial mass to be the same as that of the uncontracted halo, and set r_s and γ_DM so as to match the density and the logarithmic density slope, dlog ρ/dlog r, at r = R_e. Figure 1 shows examples of dark matter density profiles determined by solving Equation (8), along with their gNFW approximations, for a few values of ϵ, M_*, and R_e, at fixed halo mass. For a given ϵ, dark matter profiles get steeper with increasing stellar mass and with decreasing half-light radius. The gNFW approximation matches the original dark matter profile to better than 10%, at most radii. The gNFW profile can approximate well also the enclosed projected mass, which is the main quantity probed by lensing, as shown in the bottom of Figure 1. Differences between the original dark matter profile and the gNFW model are around a few percent at radii close to the median of the SLACS lenses (4.2 kpc).

Fig. 1. Stellar and dark matter distribution of example model galaxies. Top: Density profile. Bottom: Enclosed projected mass. In each top subpanel, the solid black line shows the stellar mass distribution, which follows a spherically de-projected de Vaucouleurs profile, with half-light radius indicated by the vertical dotted line. Solid coloured lines are dark matter profiles obtained by applying the adiabatic contraction prescription of Equation (8), with different values of the contraction efficiency parameter ϵ. Dashed lines are gNFW profile approximation of the contracted profile. The bottom subpanels show the relative difference between the gNFW and the corresponding adiabatically contracted profile. The legends at the bottom indicate the values of γDM and rs of each gNFW model. All galaxies have a halo mass of M200 = 1013 M⊙. Stellar mass and half-light radius vary as indicated on the top panels.

In summary, we described the dark matter halo by means of two parameters: M₂₀₀ and ϵ. To compute its lensing properties, we used the gNFW profile corresponding to the values of M₂₀₀, ϵ, M_* and R_e of the galaxy. For positive values of ϵ, the central dark matter density correlates positively with the stellar mass and negatively with the half-light radius.

2.4. Population distribution

Our inference method requires us to model the distribution of foreground galaxies and background sources that generate the SLACS lenses: P_g and P_s^(eff) in Equation (1). The former describes the parent population of galaxies among which SLACS lenses were searched for, which are early-type galaxies in the SDSS spectroscopic sample. The foreground galaxy parameters that are relevant for the problem are

$$\begin{array}{c}\hfill {\mathit{\psi }}_{\mathrm{g}}\equiv \{z_{\mathrm{g}},m_{\ast },r_{\mathrm{e}},{\alpha }_{\mathrm{sps}},m_{\mathrm{h}},ϵ,s_{\mathrm{ap}}\},\end{array}$$(11)

where we introduced the following compact notation:

$$\begin{array}{cc}\hfill m_{\ast }& \equiv log\frac{M_{\ast }^{(\mathrm{sps})}}{M_{\odot }},\hfill \end{array}$$(12)

$$\begin{array}{cc}\hfill r_{\mathrm{e}}& \equiv log\frac{R_{\mathrm{e}}}{\mathrm{kpc}},\hfill \end{array}$$(13)

$$\begin{array}{cc}\hfill m_{\mathrm{h}}& \equiv log\frac{M_{200}}{M_{\odot }},\hfill \end{array}$$(14)

$$\begin{array}{cc}\hfill s_{\mathrm{ap}}& \equiv log\frac{{\sigma }_{\mathrm{ap}}}{\mathrm{km}{\mathrm{s}}^{-1}}·\hfill \end{array}$$(15)

We chose the following form for their distribution:

$$\begin{array}{cc}\hfill {\mathrm{P}}_{\mathrm{g}}({\mathit{\psi }}_{\mathrm{g}})& =\mathcal{M}(z,m_{\ast })\mathcal{R}(r_{\mathrm{e}}|m_{\ast })\delta ({\alpha }_{\mathrm{sps}}-{\overline{\alpha }}_{\mathrm{sps}})\mathcal{H}(m_{\mathrm{h}}|m_{\ast })\hfill \\ \hfill & \times \delta (ϵ-\overline{ϵ})\mathcal{S}(s_{\mathrm{ap}}|m_{\ast },r_{\mathrm{e}},m_{\mathrm{h}}).\hfill \end{array}$$(16)

We set the factors ℳ and ℛ, describing the distribution in stellar mass and half-light radius, to be the same as in Paper I. The former is a Schechter function with a redshift-dependent truncation, while the latter is a Gaussian in r_e with the mean that scales quadratically with m_*. We asserted α_sps to be the same across all population. For the distribution in halo mass, we adopted the same model used by Sonnenfeld et al. (2018):

$$\begin{array}{c}\hfill \mathcal{N}(m_{\mathrm{h}}|m_{\ast })={\mathcal{N}}_{m_{\mathrm{h}}}({\mu }_{\mathrm{h}}(m_{\ast }),{\sigma }_{\mathrm{h}}^2),\end{array}$$(17)

where the notation 𝒩_x(μ, σ²) indicates a Gaussian distribution in x with mean μ and variance σ², and

$$\begin{array}{c}\hfill {\mu }_{\mathrm{h}}={\mu }_{\mathrm{h},0}+{\beta }_{\mathrm{h}}(m_{\ast }-11.3).\end{array}$$(18)

Sonnenfeld et al. (2018) measured the distribution in halo mass of the parent population of the SLACS lenses, by fitting the model of Equation (17) to weak lensing observations from the Hyper Suprime-Cam Subaru Strategic Program (HSC SSP Aihara et al. 2018a,b). They constrained the parameters of Equation (17) and Equation (18) to be μ_h, 0 = 13.04 ± 0.04, β_h = 1.48 ± 0.15 and σ_h = 0.31 ± 0.04, with very small covariance between them. We used these measurements as a prior on the three parameters, thus incorporating weak lensing information on the halo mass distribution of the foreground galaxies. We approximated this prior as a trivariate Gaussian with diagonal covariance matrix. One implicit assumption is that the halo mass distribution is independent of the contraction parameter ϵ. This is justified by the fact that Sonnenfeld et al. (2018) fitted both adiabatically contracted and NFW models to the weak lensing data, finding negligible differences in the inferred halo mass distribution parameters between the two models. As with the parameter α_sps, we also asserted that ϵ is the same among all galaxies. For simplicity of notation, we replace $({\overline{\alpha }}_{\mathrm{sps}},\overline{ϵ})$ with (α_sps, ϵ) from here on.

The last factor in Equation (16) is a term describing the distribution in velocity dispersion of parent population galaxies. This term is needed because the SLACS selection function depends explicitly on the observed velocity dispersion. In Paper I we predicted σ_ap by means of dynamical models, assuming a spherical mass distribution and isotropic orbits. Here we take a less assumption-heavy approach by simply proposing an empirical relation in which the velocity dispersion scales with stellar mass, size, and halo mass:

$$\begin{array}{c}\hfill \mathcal{S}(s_{\mathrm{ap}})={\mathcal{N}}_{s_{\mathrm{ap}}}({\mu }_{\sigma }(m_{\ast },r_{\mathrm{e}},m_{\mathrm{h}}),{\sigma }_{\sigma }^2),\end{array}$$(19)

with

$$\begin{array}{cc}\hfill {\mu }_{\sigma }(m_{\ast },r_{\mathrm{e}},m_{\mathrm{h}})& ={\mu }_{\sigma ,0}+{\beta }_{\sigma }(m_{\ast }-11.3)+{\xi }_{\sigma }(r_{\mathrm{e}}-{\mu }_R(m_{\ast }))\hfill \\ \hfill & +{\nu }_{\sigma }(m_{\mathrm{h}}-{\mu }_{\mathrm{h}}(m_{\ast })).\hfill \end{array}$$(20)

In the above equation, μ_R is the average r_e for galaxies of a given stellar mass, which we set to the quadratic mass-size relation of Hyde & Bernardi (2009). Equation (19) can be thought of as an extension of the fundamental plane relation that includes the halo mass, so we refer to it as the fundamental hyper-plane.

The factor P_s^(eff) in Equation (1) describes an effective redshift distribution of background sources that can be lensed into detectable multiple images, averaged over the source brightness. We used the same distribution as in Paper I, which is a Gaussian in z_s:

$$\begin{array}{c}\hfill {\mathrm{P}}_{\mathrm{s}}^{(\mathrm{eff})}(z_{\mathrm{s}})={\mathcal{N}}_{z_{\mathrm{s}}}({\mu }_{z_{\mathrm{s}}},{\sigma }_{z_{\mathrm{s}}}^2).\end{array}$$(21)

2.5. Selection function

The model for the selection probability term, P_sel, is the same as in Paper I. It is the product between the probability of detection of a lens-source pair and the probability of finding a lens in the SLACS data, given that it is detectable:

$$\begin{array}{c}\hfill {\mathrm{P}}_{\mathrm{sel}}({\mathit{\psi }}_{\mathrm{g}},z_{\mathrm{s}})=g({\mathit{\psi }}_{\mathrm{g}},z_{\mathrm{s}}){\mathrm{P}}_{\mathrm{find}}({\mathit{\psi }}_{\mathrm{g}},z_{\mathrm{s}}).\end{array}$$(22)

The factor g, which describes the probability of detection, is proportional to the strong lensing cross-section for a reference background source. Because the detection of a SLACS strong lens involved both spectroscopy and photometry, the choice of reference source requires defining both its surface brightness and spectrum. We chose as reference a point source with intrinsic broadband flux equal to the detection limit of the photometric data used for imaging follow-up, and emission line flux equal to a third of the detection limit for spectroscopy. This is the same definition used in Paper I, where we found that the exact choice of reference source did not affect the results. We then computed g by simulating the production of photometric and spectroscopic data of lensed sources, as a function of the lens parameters. As is stated in Section 2.1 and explained in greater detail in Paper I, we made the assumption that the lensing cross-section for sources different than the reference are simply obtained by rescaling g by a factor that depends on the source spectral energy distribution, and that is implicitly included in the source distribution P_s^(eff).

The lens finding probability P_find describes a selection in the estimated Einstein radius that was applied by Bolton et al. (2006) to prioritise the photometric follow-up of spectroscopically detected lens candidates. Following Paper I, we modelled it as

$$\begin{array}{c}\hfill {\mathrm{P}}_{\mathrm{find}}\left({\theta }_{\mathrm{E}}^{(\mathrm{est})}\right)=\frac{1}{1+exp[-a({\theta }_{\mathrm{E}}^{(\mathrm{est})}-{\theta }_0)]},\end{array}$$(23)

where ${\theta }_{\mathrm{E}}^{(\mathrm{est})}$ is the Einstein radius of a singular isothermal sphere with a velocity dispersion equal to the observed one:

$$\begin{array}{c}\hfill {\theta }_{\mathrm{E}}^{(\mathrm{est})}=4\pi {\left(\frac{{\sigma }_{\mathrm{ap}}^{(\mathrm{obs})}}{c}\right)}^2\frac{D_{\mathrm{ds}}}{D_{\mathrm{s}}}·\end{array}$$(24)

2.6. Practical implementation

The hyper-parameters of the model are

$$\begin{array}{c}\hfill \mathit{\eta }=\{{\alpha }_{\mathrm{sps}},ϵ,{\mu }_{\mathrm{h},0},{\beta }_{\mathrm{h}},{\sigma }_{\mathrm{h}},{\mu }_{\sigma ,0},{\beta }_{\sigma },{\xi }_{\sigma },{\nu }_{\sigma },{\sigma }_{\sigma },{\mu }_{z_{\mathrm{s}}},{\sigma }_{z_{\mathrm{s}}},{\theta }_0,a\}.\end{array}$$(25)

Table 1 provides a brief description of each of them. In order to obtain the posterior probability P(η|d), we need to evaluate integrals of the kind of Equation (4). With out parameterisation, Equation (4) becomes

$$\begin{array}{cc}\hfill \mathrm{P}(\mathbf{d}|\mathit{\eta })& ={\displaystyle \int d}{m}_{\mathrm{h}}d{m}_{\ast }d{s}_{\mathrm{ap}}\mathrm{P}({\theta }_{\mathrm{E}}^{(\mathrm{obs})}|m_{\ast },r_{\mathrm{e}},m_{\mathrm{h}},{\alpha }_{\mathrm{sps}},ϵ)\mathrm{P}\left(m_{\ast }^{(\mathrm{obs})}|m_{\ast }\right)\hfill \\ \hfill & \times \mathrm{P}\left(s_{\mathrm{ap}}^{(\mathrm{obs})}|s_{\mathrm{ap}}\right){\mathrm{P}}_{\mathrm{SL}}^{(\mathrm{eff})}(z_{\mathrm{g}},m_{\ast },r_{\mathrm{e}},m_{\mathrm{h}},s_{\mathrm{ap}},z_{\mathrm{s}}|\mathit{\eta }),\hfill \end{array}$$(26)

where we integrated over the lens and source redshift, the likelihood of which is a Dirac delta function. Then, we approximated the likelihood in s_ap as a Gaussian, which allowed us to compute the following integral over s_ap,

$$\begin{array}{cc}& I_{\sigma }(m_{\ast },r_{\mathrm{e}},m_{\mathrm{h}},\mathit{\eta })={\displaystyle \int d}{s}_{\mathrm{ap}}\mathrm{P}\left(s_{\mathrm{ap}}^{(\mathrm{obs})}|s_{\mathrm{ap}}\right)\mathcal{S}(s_{\mathrm{ap}}|m_{\ast },r_{\mathrm{e}},m_{\mathrm{h}},\mathit{\eta })\hfill \\ \hfill & =\frac{1}{\sqrt{2\pi ({\sigma }_{\sigma }^2+\mathrm{\Delta }{s}_{\mathrm{ap}}^2)}}exp\{-\frac{{({\mu }_{\sigma }(m_{\ast },r_{\mathrm{e}},m_{\mathrm{h}})-s_{\mathrm{ap}}^{(\mathrm{obs})})}^2}{2({\sigma }_{\sigma }^2+\mathrm{\Delta }{s}_{\mathrm{ap}}^2)}\},\hfill \end{array}$$(27)

where Δs_ap is the observational uncertainty associated with $s_{\mathrm{ap}}^{(\mathrm{obs})}$. We were then left with the following integral over the stellar and halo mass:

$$\begin{array}{cc}\hfill \mathrm{P}(\mathbf{d}|\mathit{\eta })& ={\displaystyle \int d}{m}_{\mathrm{h}}d{m}_{\ast }\mathrm{P}\left(m_{\ast }^{(\mathrm{obs})}|m_{\ast }\right)\mathrm{P}({\theta }_{\mathrm{E}}^{(\mathrm{obs})}|m_{\ast },r_{\mathrm{e}},m_{\mathrm{h}},{\alpha }_{\mathrm{sps}},ϵ)\hfill \\ \hfill & \times \mathcal{M}(z_{\mathrm{g}},m_{\ast })\mathcal{R}(r_{\mathrm{e}}|m_{\ast })\mathcal{H}(m_{\mathrm{h}}|m_{\ast },\mathit{\eta })I_{\sigma }(m_{\ast },r_{\mathrm{e}},m_{\mathrm{h}},\mathit{\eta })\hfill \\ \hfill & \times g({\mathit{\psi }}_{\mathrm{g}}^{(\mathrm{E})},z_{\mathrm{s}}){\mathrm{P}}_{\mathrm{find}}(s_{\mathrm{ap}}^{(\mathrm{obs})},z_{\mathrm{g}},z_{\mathrm{s}}).\hfill \end{array}$$(28)

We computed these integrals via Monte Carlo integration and importance sampling, and sampled the posterior probability with a Markov chain Monte Carlo (MCMC).

In order for P_SL^(eff) to represent a proper probability distribution, it needs to be normalised:

$$\begin{array}{c}\hfill {\displaystyle \int d}{\mathit{\psi }}_{\mathrm{g}}d{z}_{\mathrm{s}}{\mathrm{P}}_{\mathrm{SL}}^{(\mathrm{eff})}({\mathit{\psi }}_{\mathrm{g}},z_{\mathrm{s}}|\mathit{\eta })=1.\end{array}$$(29)

For each draw of values of the hyper-parameters, we computed integrals of the kind of Equation (29) via Monte Carlo integration, and rescaled the factor g in Equation (28) accordingly.

2.7. Priors

We assumed a uniform prior on log α_sps, over the range (0.0, 0.3). The lower bound corresponds to a Chabrier IMF, while the upper bound is slightly above the value corresponding to a Salpeter IMF, which is log α_sps = 0.25. For the contraction parameter ϵ, we assumed a uniform prior over the range (0, 1). The lower bound corresponds to NFW dark matter halos, while the upper bound describes maximally contracted halos. We also assumed uniform priors on the parameters describing the fundamental hyper-plane of Equation (19). For the halo mass parameters we assumed Gaussian priors, as explained in Section 2.4. We assumed a flat prior on the lens finding parameters θ₀ and log a. Finally, we fixed the source redshift distribution parameters to μ_zs = 0.48 and σ_zs = 0.215, which are the values measured in Paper I. This allowed us to greatly simplify the calculation of the integrals of Equation (29). In Paper I we found that none of the galaxy mass parameters correlated with them. Hence, we believe that this choice does not introduce any significant bias. These prior choices are summarised in Table 1.

3. Results

3.1. Mass parameters

Figure 2 shows the posterior probability in the two main parameters of the model: α_sps and ϵ. Individually, the two are essentially unconstrained, as they are degenerate with each other. This was expected since increasing the efficiency of dark matter contraction raises the inner dark matter density, which, in order to fit the Einstein radius of a lens, needs to be compensated by a corresponding decrease in stellar mass. In the limit ϵ = 0.8, the data requires log α_sps ≈ 0, which is the value of α_sps corresponding to a Chabrier IMF. At the opposite limit, models with ϵ = 0 and 0.20 ≲ log α_sps ≲ 0.25 (i.e. slightly below a Salpeter IMF) are allowed by the data.

Fig. 2. Posterior probability in the dark matter contraction and stellar population synthesis parameters. The fiducial inference is shown by the thick solid contours. Filled contours show the inference obtained by ignoring selection effects altogether. Contour levels correspond to 68% and 95% enclosed probability.

In order to gauge the importance of selection effects, we also fitted a model with no lensing selection-related terms, essentially dropping P_s^(eff) and P_sel from Equation (1). The resulting inference on (ϵ, α_sps) is shown in Figure 2 as dashed contours. At fixed ϵ, α_sps is overestimated by approximately 0.04 dex, an amount that is larger than the 1σ uncertainty. This shows that selection effects are an important component of strong lensing measurements, already at the current level of precision.

3.2. Velocity dispersion

Figure 3 shows the posterior probability in the parameters describing the distribution in velocity dispersion, or the fundamental hyper-plane. In order to interpret these results, we also predicted the values of these parameters by using the spherical Jeans equation, assuming orbital isotropy, for a few pairs of values of (ϵ, α_sps) that are consistent with the data. The Jeans model prediction can match the inferred values of μ_σ, 0, ξ_σ and ν_σ, which describe, respectively, the average velocity dispersion of galaxies at log M_*^(sps) = 11.3, and the scaling of velocity dispersion with size and halo mass. However, the Jeans model underpredicts the value of β_σ: the data favour a significantly steeper scaling of velocity dispersion with stellar mass. This mismatch could be an indication of a mass-dependent α_sps, which would be consistent with previous findings from lensing and dynamics analyses (Treu et al. 2010; Sonnenfeld et al. 2015; Posacki et al. 2015). But a trend of anisotropy with stellar mass could also explain the discrepancy.

Fig. 3. Posterior probability of the parameters describing the distribution in velocity dispersion. These are defined in Equation (19) and Equation (20). Filled purple contours: Posterior probability of the model. Solid pink contours: Posterior predicted fundamental hyperplane of the SLACS lenses. Dashed black contours: Posterior predicted fundamental hyperplane of the SLACS lenses, based on the observed (noisy) velocity dispersion. The mean and standard deviation of the marginal posterior distribution in each parameter is indicated above the corresponding histogram. Contour levels indicate regions of 68% and 95% enclosed probabilities. The three points are values of the parameters obtained by means of dynamical modelling via the isotropic spherical Jeans equation, for different pairs of values of (αsps, ϵ).

Spherical and isotropic Jeans models predict σ_σ = 0 by construction, as the choice of mass profile determines the velocity dispersion uniquely. A fundamental hyper-plane with zero intrinsic scatter is not completely ruled out by the SLACS data, but is highly disfavoured. We do expect some non-zero scatter in velocity dispersion, for the reason that galaxies are not spherical, their orientation is random, and their anisotropy profile is unlikely to be the same across the whole population.

As we explained in Section 1, one of the goals of this work is to understand whether the velocity dispersion of the SLACS lenses is biased with respect to their parent population, and whether the observed values of σ_ap are biased with respect to the truth. We did so by first obtaining the fundamental hyper-plane parameters of SLACS-like samples of lenses from posterior prediction. We drew sets of values of the model parameters from the posterior probability; then, at each draw, we generated sets of SLACS-like lenses and fitted the distribution in velocity dispersion of Equation (19) to these mock data. The resulting distribution in fundamental hyper-plane parameters is shown in Figure 3 as solid pink contours. The posterior predicted fundamental hyper-plane of the SLACS lenses has a larger value of μ_σ, 0 than that of the parent population, meaning a larger velocity dispersion at fixed stellar mass, size and halo mass. This bias was not present in the model of Paper I, because the mass distribution determined uniquely the velocity dispersion, by means of the stellar dynamics model. In other words, in the context of Paper I, SLACS lenses differed from parent population galaxies in terms of their inner structure, but had the same velocity dispersion at fixed density profile. Here we allowed for scatter between the mass properties of the lens and the velocity dispersion, therefore trading part of the bias on the mass for a bias on the velocity dispersion. We believe this new model to be more realistic, as it does not rely on the assumptions of spherical symmetry and orbital isotropy.

In Paper I we found that the velocity dispersion of the SLACS lenses is systematically overestimated. To verify whether that finding holds with this model as well, we fitted the model of Equation (19) directly to the posterior predicted noisy $s_{\mathrm{ap}}^{(\mathrm{obs})}$, instead of the noiseless s_ap. The resulting distribution is shown in Figure 3 as dashed black contours. As expected from the results of Paper I, the value of μ_σ, 0 obtained in this way is yet larger than that based on the true distribution, meaning that the observed velocity dispersion of SLACS lenses is on average larger than the true value.

In order to better illustrate the difference between these various values of μ_σ, 0 and account for their covariance, we show in Figure 4 the distribution of μ_σ, 0 of the true and observed velocity dispersion of the SLACS lenses, as a function of the corresponding value for the parent population. Averaged over the uncertainties, our model predicts that SLACS lenses have a 3% larger velocity dispersion, compared to that of galaxies of the same stellar mass, size and halo mass, and their observed velocity dispersion is a further 2% higher. Hence, intrinsic and observational bias contribute roughly equally to the bias in velocity dispersion of the SLACS lenses.

Fig. 4. Posterior predicted distribution in the observed and true μσ, 0 of SLACS lenses, as a function of the corresponding value of the parent population. The quantity μσ, 0 describes the average log σap of galaxies with log M*(sps) = 11.3, average size and average halo mass for their stellar mass.

The MCMC of the two inferences are available online¹, along with posterior predicted mock samples of the galaxy and strong lens population. These mocks can be used to verify, further investigate and better quantify any difference between SLACS and the parent population.

3.3. Goodness of fit

We checked the goodness of fit of the model by running posterior predictive tests. This means choosing a series of scalar test quantities derived from the data that we wish the model to reproduce, generating mock observations by sampling from the posterior probability, and measuring the probability of the model to produce more extreme values than the test quantities. A very low or very large probability indicates that the model is unlikely to reproduce those aspects of the data. We defined the test quantities by focusing on the Einstein radius. We chose the average, the standard deviation, the 10%- and the 90%-ile of the observed distribution in Einstein radius. Figure 5 shows the posterior predicted distribution of the four test quantities. In all cases, the observed values are well within the posterior predicted distributions. Hence, the model is able to reproduce these data.

Fig. 5. Posterior predictive test of the goodness-of-fit. Vertical lines indicate the observed values of the four test quantities. In each panel, the percentage to the left and right of the vertical line indicate the fraction of posterior predicted samples for which the mock test quantity is smaller or larger than observed.

3.4. Prediction of other lensing observables

Since the SLACS sample and the data used for our analysis are not sufficient to individually determine α_sps and ϵ, it can be useful to investigate whether current or future data from SLACS or SLACS-like samples can break this degeneracy. We did this, too, with posterior prediction. We first focused on the lensing-only power-law slope, which is the value of the density slope γ_PL that is obtained by fitting a lens model with a power-law radial density profile, ρ(r)∝r^−γ_PL, to imaging data of a strong lens. Shajib et al. (2021), Etherington et al. (2022) and Tan et al. (2024) measured γ_PL for 21, 29, and 28 of the SLACS lenses included in this work, respectively. Some of these measurements could, in principle, be used as additional constraints. But since our model is not a power-law, we must first understand what exact property of the lenses do the γ_PL capture. Sonnenfeld (2018) and Kochanek (2020) argued that, when a power-law model is fitted to well-resolved, high signal-to-noise ratio lensed images around nearly circularly symmetric lenses, γ_PL is constrained by the ratio of radial magnification between the main arc and its counter-image. Gomer et al. (2023) showed by means of simulations that this is the case, even for moderately elliptical lenses, provided that the fitted model has sufficient flexibility in the azimuthal direction. The radial magnification ratio of images produced by a circular lens can be expressed as a function of the following ratio of radial derivatives of the lensing potential ψ, evaluated at the Einstein radius:

$$\begin{array}{c}\hfill {\xi }_{\mathrm{rad}}={\theta }_{\mathrm{E}}\frac{{\psi }_{\mathrm{E}}″\prime }{1-{\psi }_{\mathrm{E}}″}·\end{array}$$(30)

For a power-law model,

$$\begin{array}{c}\hfill {\xi }_{\mathrm{rad}}={\gamma }_{\mathrm{PL}}-2.\end{array}$$(31)

Therefore, we can define an equivalent lensing-only power-law slope for our model, by computing ψ_E″ and ψ_E″′, obtaining ξ_rad from Equation (30), and setting γ_PL = ξ_rad + 2.

We did this for our model, while fixing the value of the dark matter contraction parameters to two extreme values: ϵ = 0.0 and ϵ = 0.8. Figure 6 shows the posterior predicted distribution in γ_PL, as a function of the stellar surface mass density, Σ_*^(sps) = M_*^(sps)/(2πR_e²). The two quantities are positively correlated: this is expected, since Σ_*^(sps) correlates with the total density slope γ obtained from lensing and dynamics (Sonnenfeld et al. 2013), which should be closely related to γ_PL. The ϵ = 0.0 and ϵ = 0.8 model lie on different regions of the Σ_*^(sps) − γ_PL parameter space. Although there is overlap between the two posterior predicted distributions, most of their widths are due to intrinsic scatter within the population of lenses rather than the parameter uncertainty. Therefore, the two models stand apart when comparing the average γ_PL of the lens population. This means that γ_PL observations could already help distinguish between uncontracted or significantly contracted dark matter profiles. Unfortunately, as Figure 6 shows, the distribution of current measurements of γ_PL of the SLACS lenses bears very little resemblance to that predicted with our model. The observed γ_PL appear to be anti-correlated with Σ_*^(sps). This is at odds with lensing and dynamics constrains, as already pointed out by Etherington et al. (2023).

Fig. 6. Posterior predicted distribution in the lensing-only power-law slope, as a function of stellar mass density. Simulated values of γPL have been obtained from Equation (31). The two sets of contours correspond to the posteriors obtained by fitting models with a fixed value of ϵ, as indicated in the legend. Error bars are measurements from Shajib et al. (2021), Etherington et al. (2022), Tan et al. (2024).

We quantified the discrepancy with another posterior predictive test, in which we generated noisy measurements of γ_PL; fitted them with a linear relation of the kind

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{PL}}^{\left(\mathrm{obs}\right)}\left({\mathrm{\Sigma }}_{\ast }^{(\mathrm{obs})}\right)={\mathcal{N}}_{{\gamma }_{\mathrm{PL}}^{(\mathrm{obs})}}({\gamma }_{\mathrm{PL},0}+{\beta }_{{\gamma }_{\mathrm{PL}}}({\mathrm{\Sigma }}_{\ast }^{(\mathrm{obs})}-9),{\sigma }_{{\gamma }_{\mathrm{PL}}}^2);\end{array}$$(32)

and counted the fraction of times in which the parameters have more extreme values than those obtained by fitting Equation (32) to the observed samples of Shajib et al. (2021), Etherington et al. (2022), and Tan et al. (2024). The results are shown in Figure 7. Only the intercept γ_PL, 0 can be reproduced by the model, and only for two of the datasets. All three of the observed datasets show an anti-correlation between γ_PL and Σ_*^(sps), corresponding to a negative value of β_γPL. Out of 10 000 posterior predicted mocks, only a handful have a value of β_γPL that is smaller than the observed ones.

Fig. 7. Posterior predicted distribution of the coefficients describing the γPL − Σ*(sps) relation of SLACS lenses. The coefficients are defined in Equation (32). For each coefficient, three distributions are shown, each to be compared with the value obtained by fitting the measurements of Shajib et al. (2021), Etherington et al. (2022) and Tan et al. (2024). The colour of each distribution corresponds to that of the data points in Figure 6. The distributions are different because the observed samples and the observational uncertainties vary from one work to the other. In each panel, the percentages to the left and right of the vertical lines indicate the fraction of posterior predicted samples for which the mock test quantity is smaller or larger than observed.

This discrepancy could be explained by the model being inaccurate. However, given the fact that many of the measurements of γ_PL from the different works disagree on individual lenses (see Figure 6 of Tan et al. 2024), we believe that there are problems with the observations, or their interpretation. One possible source of error could lie in the choice of the model for the measurement of γ_PL. Shajib et al. (2021), Etherington et al. (2022) and Tan et al. (2024) fitted lens models with elliptical symmetry to the SLACS lenses. If the true azimuthal structure of the lenses deviates from elliptical symmetry, then the value of γ_PL obtained can no longer interpreted as a measurement of the radial profile via Equation (31). As was argued by Kochanek (2021) and shown in simulations by Van de Vyvere et al. (2022) and Gomer et al. (2023), fits of elliptical models to lenses with more complex azimuthal structure lead to biases. Other sources of error, which could explain the differences between the three measurements on the same lenses, are the modelling of the point spread function and of the background source (see Galan et al. 2024, for a discussion). We conclude that, at the present, we cannot rely on existing measurements of γ_PL, until it is shown that they are reproducible and robust to the effects mentioned above. Nevertheless, measurements of the radial magnification ratio on large sample of lenses have the potential to break the degeneracy between the stellar IMF and the dark matter profile, as also shown by Sonnenfeld & Cautun (2021).

Another piece of information that could help break the degeneracy between α_sps and ϵ is the number density of lenses. As has been shown by Zhou et al. (2024), the number of lenses that are found in a given area of the sky can provide complementary constraints on the structural parameters of the galaxy population. In the context of our model, the number density of lenses in a given survey is given by the following product:

$$\begin{array}{c}\hfill n_{\mathrm{lens}}=n_{\mathrm{g}}{n}_{\mathrm{s}}{\displaystyle \int d}{\mathit{\psi }}_{\mathrm{g}}d{z}_{\mathrm{s}}{\mathrm{P}}_{\mathrm{g}}({\mathit{\psi }}_{\mathrm{g}}){\mathrm{P}}_{\mathrm{s}}^{(\mathrm{eff})}(z_{\mathrm{s}})g({\mathit{\psi }}_{\mathrm{g}},z_{\mathrm{s}}){\mathrm{P}}_{\mathrm{find}}({\mathit{\psi }}_{\mathrm{g}},z_{\mathrm{s}}),\end{array}$$(33)

where n_g and n_s are the projected number density of foreground galaxies and background sources, respectively, and g is the geometrical component of the lensing cross-section. We do not have precise estimates of n_g, which depends on the effective area of the SLACS survey, and especially of n_s. Nevertheless, we can study how n_lens varies as a function of the lens population parameters for arbitrary values of n_g and n_s. We computed this quantity by Monte Carlo integration, in the same way as we obtained the normalisation of P_SL^(eff). Figure 8 shows the posterior predicted n_lens as a function of ϵ. The two quantities are positively correlated: when averaged over the uncertainty on the other parameters, the model with ϵ = 0.8 predicts a 20% larger n_lens compared to the ϵ = 0.0 case. Although at fixed ϵ the predicted distribution in n_lens has a relatively large spread, most of this spread is due to the uncertainty on the selection function parameters, in particular θ₀. This means that, if the selection function was known exactly and both n_g and n_s were known precisely, a sample of ∼100 lenses would allow us to distinguish between the two extreme scenarios allowed by the current constraints. Although this might be difficult to achieve for the SLACS sample, future surveys with better-defined selection criteria might be more fortunate in this regard.

Fig. 8. Posterior predicted number density of SLACS lenses, as a function of the contraction efficiency parameter. The dashed line shows the average nlens, obtained by marginalising over the uncertainty on the model parameters.

4. Discussion and conclusions

We used the SLACS sample to revisit the issue of the relative contribution of baryonic and dark matter to the mass budget of early-type galaxies, in view of new insight on its selection function obtained in Paper I. We tried to constrain the stellar mass-to-light ratio, described by the average stellar population synthesis mismatch parameter α_sps, and the amount of dark matter contraction with respect to an NFW profile, described by the contraction efficiency parameter ϵ. Although we were not able to constrain either of them individually, we put a tight bound on their combination that allowed us to rule out large regions in the α_sps − ϵ parameter space. A stellar IMF heavier than Salpeter (i.e. that produces a larger mass-to-light ratio), corresponding to log α_sps > 0.25, is disfavoured by the data. Although the shape of the ϵ − α_sps correlation suggests that larger values of α_sps might be allowed if we let ϵ to extend below zero, which would correspond to expansion of the halo, we have a strong prior against such a solution, because hydrodynamical simulations predict contraction.

At fixed ϵ, our inferred value of α_sps is generally lower than that of most other works in the literature. For instance, Shajib et al. (2021) measured log α_sps = 0.26 on average, for a halo profile with no contraction or expansion, in their analysis of 21 SLACS lenses. The difference with respect to our measurement, 0.04 dex, is identical to the shift in the posterior probability due to the inclusion of the terms related to the selection function. Hence, it can be attributed entirely to it. However, Sonnenfeld et al. (2019), using the same method employed in this work, measured an average log α_sps − 0.04 ± 0.11 on a different sample of lenses, under the assumption of an NFW profile for the dark matter distribution. Their value should be compared with our ϵ = 0 (i.e. NFW) solution, which is significantly larger at 0.20 < log α_sps < 0.24. This discrepancy could be due to intrinsic differences between the two galaxy samples, or, more likely, to systematic errors in the measurement of the stellar population synthesis stellar mass. Independent analyses of the two samples are needed to clarify this issue. In light of this difference, we recommend caution when comparing the value of α_sps from this work to samples of galaxies different from the SLACS lenses.

We also investigated how the velocity dispersion of the SLACS lenses compares to that of parent population galaxies, at fixed mass profile. We found that SLACS lenses have a 3% larger velocity dispersion, and that their observed velocity dispersion is biased upwards by another 2%. These biases pose challenges for the use of SLACS lenses in the aid of time-delay cosmography measurements. Birrer et al. (2020) used the SLACS sample to constrain the mass-sheet transformation parameter, λ_mst, which describes how the density profile of the lenses deviates from a pure power-law model. A positive bias on the velocity dispersion leads to an overestimate of λ_mst, which in turn leads to H₀ being overestimated. The proper way to correct for this bias would be to use priors that take the SLACS selection function into account, although that might be difficult in practice. One might hope to remove at least the observational bias with higher signal-to-noise spectroscopic measurements that supersede the SDSS-based ones on which our analysis is based. However, whether this can be done or not depends on what exact process generated the observational scatter in the SDSS velocity dispersion measurements. In the ideal case of uncorrelated statistical noise, then the uncertainty should go to zero in the limit of high signal-to-noise. However, if systematic effects are at play, such as a mismatch between the stellar templates used for the velocity dispersion fit and the intrinsic spectra of the lenses, correcting for this bias might be difficult, as it might still be present in better quality data.

Knabel et al. (2024) have recently obtained spatially-resolved spectroscopic data on 13 SLACS lenses with the Keck Cosmic Web Imager (KCWI), and compared SDSS-based measurements of σ_ap with the same quantity obtained from these higher quality data. They found that the SDSS σ_ap^(obs) are under-estimated by a few percent. This bias goes in the opposite direction than the one predicted by our model. Assuming that the same bias applies to all SDSS spectra of massive early-type galaxies, this does not change the conclusions of our study, as it would simply shift the overall value of μ_σ, 0, but not the relative difference between the SLACS lenses and the parent population.

Finally, we explored ways to break the degeneracy between ϵ and α_sps with future lensing data. One observable quantity that could help in this regard is the lensing-only density slope, γ_PL, to be obtained by measuring the ratio of radial magnification between the multiple images. Although such measurements already exist for a fraction of the SLACS lenses, they are completely inconsistent with our model. In particular, while the model predicts a positive correlation between stellar mass density and γ_PL, the observations show an anti-correlation. While we have reasons to believe that there are problems with some of these measurements, it would be interesting to also explore whether it is possible to find models that better match the observations. We leave that to future work.

In summary, this analysis constitutes a shift in strategy in the analysis of samples of strong lenses, moving from joint lensing and dynamics to purely lensing-based modelling. This choice allowed us to focus exclusively on projected quantities, avoiding the need to predict difficult properties to measure, such as the stellar anisotropy and the three-dimensional structure of lens galaxies. We modelled the distribution of stellar and dark matter of the SLACS lenses, together with their selection function. This gave us a twofold advantage: it allowed us to correct for selection effects, which we found to have a significant impact when compared to the observational uncertainty, and to incorporate weak lensing information from the parent population of SLACS lens galaxies. The approach used in this work can serve as a solid foundation on which to build our understanding of galaxy structure with upcoming lensing data. Combined with a larger number of lenses, a yet better understanding of selection effects, and number density or reliable radial magnification information, it will make it possible to break the degeneracy between the stellar IMF and the dark matter density profile in early-type galaxies.

Acknowledgments

This work was supported by the National Key R&D Program of China (No. 2023YFA1607800, 2023YFA1607802).

References

All Tables

All Figures

Fig. 1. Stellar and dark matter distribution of example model galaxies. Top: Density profile. Bottom: Enclosed projected mass. In each top subpanel, the solid black line shows the stellar mass distribution, which follows a spherically de-projected de Vaucouleurs profile, with half-light radius indicated by the vertical dotted line. Solid coloured lines are dark matter profiles obtained by applying the adiabatic contraction prescription of Equation (8), with different values of the contraction efficiency parameter ϵ. Dashed lines are gNFW profile approximation of the contracted profile. The bottom subpanels show the relative difference between the gNFW and the corresponding adiabatically contracted profile. The legends at the bottom indicate the values of γDM and rs of each gNFW model. All galaxies have a halo mass of M200 = 1013 M⊙. Stellar mass and half-light radius vary as indicated on the top panels.

In the text

	Fig. 2. Posterior probability in the dark matter contraction and stellar population synthesis parameters. The fiducial inference is shown by the thick solid contours. Filled contours show the inference obtained by ignoring selection effects altogether. Contour levels correspond to 68% and 95% enclosed probability.
In the text

Fig. 3. Posterior probability of the parameters describing the distribution in velocity dispersion. These are defined in Equation (19) and Equation (20). Filled purple contours: Posterior probability of the model. Solid pink contours: Posterior predicted fundamental hyperplane of the SLACS lenses. Dashed black contours: Posterior predicted fundamental hyperplane of the SLACS lenses, based on the observed (noisy) velocity dispersion. The mean and standard deviation of the marginal posterior distribution in each parameter is indicated above the corresponding histogram. Contour levels indicate regions of 68% and 95% enclosed probabilities. The three points are values of the parameters obtained by means of dynamical modelling via the isotropic spherical Jeans equation, for different pairs of values of (αsps, ϵ).

In the text

	Fig. 4. Posterior predicted distribution in the observed and true μσ, 0 of SLACS lenses, as a function of the corresponding value of the parent population. The quantity μσ, 0 describes the average log σap of galaxies with log M*(sps) = 11.3, average size and average halo mass for their stellar mass.
In the text

	Fig. 5. Posterior predictive test of the goodness-of-fit. Vertical lines indicate the observed values of the four test quantities. In each panel, the percentage to the left and right of the vertical line indicate the fraction of posterior predicted samples for which the mock test quantity is smaller or larger than observed.
In the text

Fig. 6. Posterior predicted distribution in the lensing-only power-law slope, as a function of stellar mass density. Simulated values of γPL have been obtained from Equation (31). The two sets of contours correspond to the posteriors obtained by fitting models with a fixed value of ϵ, as indicated in the legend. Error bars are measurements from Shajib et al. (2021), Etherington et al. (2022), Tan et al. (2024).

In the text

Fig. 7. Posterior predicted distribution of the coefficients describing the γPL − Σ*(sps) relation of SLACS lenses. The coefficients are defined in Equation (32). For each coefficient, three distributions are shown, each to be compared with the value obtained by fitting the measurements of Shajib et al. (2021), Etherington et al. (2022) and Tan et al. (2024). The colour of each distribution corresponds to that of the data points in Figure 6. The distributions are different because the observed samples and the observational uncertainties vary from one work to the other. In each panel, the percentages to the left and right of the vertical lines indicate the fraction of posterior predicted samples for which the mock test quantity is smaller or larger than observed.

In the text

	Fig. 8. Posterior predicted number density of SLACS lenses, as a function of the contraction efficiency parameter. The dashed line shows the average nlens, obtained by marginalising over the uncertainty on the model parameters.
In the text