# Cosmology from Large Populations of Galaxy-Galaxy Strong Gravitational
Lenses

Tian Li, Thomas E. Collett, Coleman M. Krawczyk and Wolfgang Enzi
Institute of Cosmology and Gravitation, University of Portsmouth, Burnaby Rd,
Portsmouth PO1 3FX, UK
<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present a forecast analysis on the feasibility of measuring the
cosmological parameters with a large number of galaxy-galaxy scale strong
gravitational lensing systems. Future wide area surveys are expected to
discover and measure the properties of more than 10 000 strong lensing
systems. We develop a hierarchical model that can simultaneously constrain the
lens population and cosmological parameters by combining Einstein radius
measurements with stellar dynamical mass estimates for every lens.
Marginalizing over the lens density profiles and stellar orbital anisotropies,
we find that $w$ can be constrained to a precision of $0.11$ with 10 000
galaxy-galaxy lens systems, which would be better than any existing single-
probe constraint. We test our method on 161 existing lenses, finding
$w=-0.96\pm 0.46$. We also show how to mitigate against the potential
systematic of redshift evolution in the mean lens density profile of the
population.

###### keywords:

(cosmology:) cosmological parameters – cosmology: observations – gravitational
lensing: strong – galaxies: structure

††pubyear: 2023††pagerange: Cosmology from Large Populations of Galaxy-Galaxy
Strong Gravitational Lenses–A

##  1 Introduction

The acceleration of the universe has been discovered through Type-Ia
supernovae (Perlmutter et al., 1999; Riess et al., 1998), and has been
measured by several observational probes, including Cosmic Microwave
Background (CMB) Anisotropies, Baryon Acoustic Oscillations, Weak
Gravitational Lensing, Galaxy Clustering, and Redshift Space Distortion
(Mortonson et al., 2013). These observations concluded that under flat
$\Lambda$CDM model, the dark energy density makes up 70$\%$ of the universe
today, and has an equation of state of $w\approx-1$. However, the exact nature
of dark energy and dark matter still remains unknown. The 5 $\sigma$
discrepancy of Hubble Constant between Planck’s CMB measurement (Planck
Collaboration et al., 2014, 2016, 2020) and local measurements of supernovae
(Riess et al., 2011; Freedman et al., 2012; Riess et al., 2016) also suggest
potential new physics beyond the $\Lambda$CDM model.

In addition to the above methods, galaxy scale strong gravitational lensing
provides an independent probe to constrain the cosmological parameters. Strong
gravitational lensing occurs when two galaxies align perfectly to our line of
sight, such that the light from the background source galaxy will be distorted
and magnified by the foreground lens galaxy, resulting in multiple sources or
an Einstein ring in the observed image (Einstein, 1936; Zwicky, 1937a, b). The
radius of the Einstein ring (Eintein radius), relative positions, flux ratios,
and time delays between multiple images depend both on the gravitational
potential of the lens galaxy and angular diameter distances between the
observer, lens galaxy, and source galaxy. It is the dependence on the angular
diameter distances that makes strong lensing sensitive to cosmological
parameters. Figure 1 illustrates the sensitivity of the Einstein radius to the
equation of state of dark energy. The most well-studied method of constraining
cosmology with strong lensing is time-delay cosmography, which uses the
temporal variation of a gravitationally lensed quasar or supernova to
constrain the Hubble constant (Refsdal, 1964; Birrer et al., 2022; Treu &
Marshall, 2016; Treu et al., 2022). Alternatively, the equation of state of
dark energy, $w$, can be measured from systems with sources at multiple
redshifts (Gavazzi et al., 2008): Collett & Auger (2014) used a single double
source plane lens to infer $w=-0.99_{-0.22}^{+0.19}$. Aside from single lens
analyses, statistical analyses of the ensemble of lens systems can also place
bounds on cosmological parameters. Oguri et al. (2008), Chae et al. (2002),
and Chae (2007) studied lenses from Cosmic Lens All-sky Survey (CLASS, Myers
et al., 2003) and The Sloan Digital Sky Survey Quasar Lens Search (SQLS, Inada
et al., 2012). They measured $w$ through comparing the empirical distribution
of image separations in observed samples of lenses with theoretical models.
Since the number of lens systems is low, the constraint on $w$ using this
method is weak (e.g., $w=-1.1\pm 0.6_{-0.5}^{+0.3}$, Oguri et al. (2008)).

The current cosmology analyses with strong lensing systems are limited by
several factors. The primary issue is that the sample of known galaxy-scale
lenses is only a few hundred systems, discovered in several surveys with
heterogeneous selection functions. Once a sample of time-delay or compound
lenses has been selected only a handful of suitable systems have adequate data
for precision cosmography, e.g the latest TDCOSMO sample has only 6 time-delay
lenses (Birrer et al., 2020). Accurate and efficient lens modelling is also a
challenge. The mass distribution in the lens must be inferred to convert
lensing observables into cosmological constraints. This challenge is
compounded by the mass-sheet degeneracy, where different mass models can
produce identical strong lensing observables but imply different cosmological
parameters (Schneider & Sluse, 2013).

Figure 1: The Einstein radius as a function of the source redshift, for
different values of the equation of state of dark energy. All other parameters
are fixed. We assume the redshift of lens galaxy is 0.2, and has an Einstein
mass of $10^{11}\mathrm{M}\odot$.

Future optical imaging surveys, including Euclid and LSST, are predicted to
discover more than $10^{6}$ galaxy-galaxy strong lensing systems (Collett,
2015). Euclid will provide high-resolution images, such that lens modelling
could plausibly be performed for every lens, without the need for additional
imaging data. However, almost all strong lensing science requires accurate
lens and source redshifts, and the 4MOST Strong Lensing Spectroscopic Legacy
Survey (4SLSLS) is expected to obtain the spectrum of $\approx$ 10 000 strong
lensing systems (Collett et al., 2023).

Combining lensing and stellar dynamics opens up a new statistical method to
constrain cosmology. Gravitational lensing determines the mass inside Einstein
Radius, and stellar kinematics determines the gravitational potential within
which the stars are moving (Koopmans, 2006). Assuming general relativity, the
dynamic mass enclosed with the Einstein radius must equal the gravitational
mass measured using the Einstein radius.

Thus the combination of lensing and dynamical observables are sensitive to the
mass profile of the lens, the orbital properties of its stars, and the
cosmological distances (Futamase & Yoshida, 2001; Grillo et al., 2008). The
challenge inherent to this method is that the mass profile of lenses and the
orbital profiles of their stars are not well known. To produce cosmological
constraints, either assumptions must be made about the lenses, or the
cosmological parameters and lens properties must be inferred simultaneously.
Biesiada et al. (2010) applied the lensing plus dynamics method on 20 lens
systems. By assuming a SIS mass profile for all lens galaxies, they found
$\Omega_{M}$ = 0.27 $\pm$ 0.28, and $w=-0.63\pm 0.45$. Cao et al. (2015a) and
Chen et al. (2019) improved on this by fitting $\sim 100$ lenses with powerlaw
density profiles.

In this paper, we investigate how well the combination of Euclid and 4MOST
data for 10 000 lenses can constrain the cosmological parameters. We employ a
Bayesian hierarchical model and simultaneously fit for the cosmological
parameters and the ensemble properties of the lens galaxies, including the
density profile slope and the stellar orbital anisotropy. It is hard to
perform detailed modelling for 10 000 lens systems, so we assume that we can
only use catalog-level data in our analysis.

The rest of the paper is organized as follows. In Section 2, we present the
mass model of the lens galaxy and the equation that relates the galaxy’s
velocity dispersion to cosmological parameters. We then introduce the
properties of the mock data and discuss potential future surveys for data
acquisition. In Section 3, we describe the hierarchical model used to
simultaneously fit the lens galaxy properties and cosmology. In Section 4, we
present the results obtained under different cosmological models and data
measurement accuracies. The final section summarizes the main conclusions. In
this paper, the fiducial cosmology for the mock data set is as follows:
$\Omega_{\mathrm{M}}=0.3$, $\Omega_{\Lambda}=0.7$, $\Omega_{k}=0$, and $w=-1$.

##  2 Theory

The mass enclosed within the Einstein ring can be related to Einstein radius
through:

$\theta_{\mathrm{E}}=\sqrt{\frac{4GM\left(\theta_{\mathrm{E}}\right)}{c^{2}}\frac{D_{\mathrm{ls}}}{D_{\mathrm{l}}D_{\mathrm{s}}}}$
(1)

where $D_{l}$ is the angular diameter distance of the lens galaxy, $D_{ls}$ is
the angular diameter distance between lens and source, and $D_{s}$ is that
between observer and source. $\theta_{E}$ is the angular Einstein radius.
Mcyan($\theta_{\mathrm{E}}$) is the galaxy mass enclosed within the Einstein
radius, and the angular diameter distance is :

$\displaystyle
D_{\mathrm{ij}}=\frac{c/H_{0}}{\left(1+z_{\mathrm{j}}\right)}\left(\frac{\operatorname{sinn}\left(\sqrt{\left|\Omega_{k}\right|}\int_{z_{\mathrm{i}}}^{z_{\mathrm{i}}}\frac{\mathrm{d}z}{E(z)}\right)}{\sqrt{\left|\Omega_{k}\right|}}\right)$
(2)

where sinn(x) = sin(x), x, or sinh(x) for open ($\Omega_{k}<0$), flat
($\Omega_{k}=0$), or closed ($\Omega_{k}>0$) universes respectively, and E(z)
is the normalised Hubble parameter:

$\displaystyle
E(z)=\sqrt{\Omega_{\mathrm{M}}(1+z)^{3}+\Omega_{k}(1+z)^{2}+\Omega_{\Lambda}(1+z)^{3(1+w)}}$
(3)

When combining lensing and dynamics, the cosmological model is not directly
probed by the measurement of a single distance, but instead through a ratio of
distances $\frac{D_{s}}{D_{ls}}$. However, to make this measurement we need a
model that can connect the dynamical data with the lensing mass. Since lenses
are typically elliptical galaxies (ETGs) with E/S0 morphologies (Oguri &
Marshall, 2010), we use a power profile for both total mass density and
stellar luminosity profile (Koopmans, 2006):

$\displaystyle\begin{array}[]{l}\rho(r)=\rho_{0}\left(\frac{r}{r_{0}}\right)^{-\gamma}\\\
\nu(r)=\nu_{0}\left(\frac{r}{r_{0}}\right)^{-\delta}\\\
\beta(r)=1-\frac{\sigma_{\theta}^{2}}{\sigma_{r}^{2}}\end{array}.$ (7)

where $\rho(r)$ is the mass density (include dark matter) distribution
function. $\nu(r)$ is the luminosity density of stars. $\beta(r)$ is the
anisotropy of the stellar velocity dispersion (stellar orbital anisotropy),
where $\sigma_{\theta}$ and $\sigma_{r}$ are the tangential and radial
velocity dispersion, respectively. $\beta$ ranges from +1 to $-\infty$, where
$\beta=0$ corresponds to the "isotropic" case, $\beta=1$ corresponds to a
galaxy with pure circular stellar movement, $\beta=-\infty$ means that stars
in the galaxy only have radial movement. In the scenario where
$\gamma=\delta=2$ and $\beta=0$, the mass model reduces to the Singular
Isothermal Sphere (SIS) model, which is a commonly used approximation for the
mass profiles of elliptical galaxies (e.g. Auger et al. (2010)).

After solving the spherical Jeans equation and substituting the dynamical mass
into equation (2) we get (Koopmans, 2006):

$\sigma_{\|}^{2}(R_{A})=\frac{c^{2}}{2\sqrt{\pi}}\frac{D_{s}}{D_{ls}}\theta_{E}\times
f\left(\gamma,\delta,\beta\right)\left(\frac{\theta_{A}}{\theta_{E}}\right)^{2-\gamma},$
(8)

where

$\begin{array}[]{r}f\left(\gamma,\delta,\beta\right)=\frac{3-\delta}{(\xi-2\beta)(3-\xi)}\left[\frac{\Gamma[(\xi-1)/2]}{\Gamma(\xi/2)}-\beta\frac{\Gamma[(\xi+1)/2]}{\Gamma[(\xi+2)/2]}\right]\\\
\times\frac{\Gamma(\gamma/2)\Gamma(\delta/2)}{\Gamma[(\gamma-1)/2]\Gamma[(\delta-1)/2]}\end{array}$
(9)

$\sigma_{\|}^{2}(R_{A})$ is the luminosity averaged line-of-sight velocity
dispersion (LOSVD) measured in a circular fibre of radius $\theta_{A}$ 111For
4MOST the fibre diameter is 1.45 arcseconds (de Jong et al., 2012).. $\Gamma$s
are Gamma functions, and $\xi=\gamma+\delta-2$.

Solving Equation 8, we find that for a fixed galaxy mass and surface
brightness profile, a steeper density profile (larger $\gamma$) and a higher
velocity anisotropy (larger $\beta$) lead to a higher stellar velocity
dispersion. Figure 2 shows the value of Equation 9 as a function of $\gamma$
and $\beta$, assuming $\delta=2.173$, which is the typical value of SLACS
lenses.

Figure 2: The value of the function that links Einstein radius and velocity
disperion as a function of the orbital anisotropy and lens density profile
slope. It is given in equation 9 as a function of $\gamma$ and $\beta$. We
have fixed $\delta=2.173$. The white regions are unphysical as the luminosity-
weighted mass diverges here. Figure 3: Schematic of our hierarchical model.
Solid blue boxes are parameters of our model, and solid black boxes are
observables or predictions of our model. The black dashed box includes the
hyperparameters of the model, the parent lens population parameters and the
cosmological parameters. The green dashed box contains the properties of
individual lens systems, including the nuisance parameters $\gamma_{i}$ and
$\beta_{i}$ of individual lens galaxies, which are drawn from the population
parent distribution. Modeled velocity dispersions and Einstein radii can be
computed through combining the parameters of the individual lenses and
cosmology. Comparing these modelled values to the observed values, a
likelihood function can be built up, which we sample with a Hamiltonian Monte
Carlo analysis.

###  2.1 Mock data

To build our mock sample of lenses, we use the simulated strong lensing
population forecast to be observed by Euclid from Collett (2015). We remove
all systems with a source redshift greater than 1.5 since the [OII] doublet is
redshifted out of the 4MOST wavelength range at this redshift. The sample is
built assuming lenses are uniformly distributed in co-moving volume, follow
the observed velocity dispersion function of SDSS (Choi et al., 2007), and
have an SIS density profile. The source properties use the LSST sky
simulations, with redshifts and number counts matched to observations
(Connolly et al., 2010). These assumptions produce a realistic lens
population, but lack the complexities of a non-SIS density profile, or of
stellar anisotropy. Therefore, we take only the lens galaxy redshift, source
galaxy redshift, and the Einstein radius in the data set. We assign $\gamma$,
$\beta$, and $\delta$ values to each lens galaxy to produce new velocity
dispersions.

Observations shows that $\gamma$ has a distribution of $2.078\pm 0.16$ (Auger
et al., 2010), $\delta$ has a distribution of $2.173\pm 0.085$ (Chen et al.,
2019), and $\beta=0.18\pm 0.13$ (Schwab et al., 2010; Bolton et al., 2006). We
generate mock galaxies using these distributions. The true values of velocity
dispersion are generated through equation (6). Then, we add random noise to
these values to simulate measurement error. The measurement error of Einstein
radius is set as 0.01 arcsec, and the error on velocity dispersion of 4MOST is
10km/s. We neglect errors on $\delta$ and the redshifts since they can be
measured with high accuracy. In our fiducial model, we assume that $\beta$ and
$\gamma$ for each lens are unknown. The full setup of mock data is shown in
Table 1.

Property name | Distribution | Measurement error
---|---|---
$\gamma$ | 2 $\pm$ 0.16 | 0.02*
$\delta$ | 2.173 $\pm$ 0.085 | None
$\beta$ | 0.18 $\pm$ 0.13 | $-$
$\sigma_{v}$ | Equation 8 | 10 km/s
$\theta_{\mathrm{E}}$ | Collett (2015) | 0.01”
$z_{l}$, $z_{s}$ | Collett (2015) | None
*$\gamma$ is only treated as an observable in Section 5
Table 1: The parameters of the mock strong lensing systems used in this work.
Einstein radius and redshifts are taken from Collett (2015), with source
redshifts below $z_{s}<1.5$. The velocity dispersions are computed from the
other parameters using equation 8. Among the above properties, $\delta$,
redshift, Einstein radius, and Velocity dispersion are treated as observables.
The measurement error of $\delta$ and redshift are negligible. $\gamma$ can
also be measured through detailed lens modelling, but we treat it as a free
parameter except in the end of Section 5.

##  3 Hierarchical Model

Since we do not know the mass profile or orbital anisotropy of individual
lenses apriori, and they are not easily measured without detailed lens
modelling and integral field unit kinematics, we instead require a
hierarchical model to connect the measured Einstein radius and the velocity
dispersion of every lens with the underlying cosmological parameters of the
Universe. In fact, even the population properties for lens density profiles
and stellar anisotropies are not well known once the cosmology is allowed to
be a free parameter. For our population model, we assume that density profile
slopes and stellar anisotropies follow a Gaussian distribution. Our
hierarchical model must therefore fit for the ensemble mean density profile
slope, $\langle\gamma\rangle$, and scatter, $\sigma_{\gamma}$, the ensemble
mean anisotropy $\langle\beta\rangle$, and scatter, $\sigma_{\beta}$, the
individual slopes, $\gamma_{i}$, and anisotropies, $\beta_{i}$, of each lens,
and the underlying cosmological model parameters. Figure 3 illustrate the
structure of our hierarchical model, with the observed Einstein radii and
velocity dispersions linked to the model parameters through equation 8.

###  3.1 Sampling the parameters of our hierarchical model

Our hierarchical model has a large number of free parameters ( $\sim$ 3 per
observed lens $\times$ 10 000 lenses), traditional MCMC methods are not
sufficient for exploring such a posterior distribution. To that end, we make
use of the computational frameworks initially developed for large machine
learning models, specifically the NumPyro (Phan et al., 2019; Bingham et al.,
2019) probabilistic programming language (PPL).

NumPyro is an extension of the Pyro (Bingham et al., 2019) framework that uses
JAX222http://github.com/google/jaxfor automatic differentiation and adds
various Hamiltonian Monte-Carlo (HMC) (Duane et al., 1987; Brooks et al.,
2011) sampling methods. Automatic differentiation is a technique that
efficiently computes the partial derivative of a function without the need to
write any new code (Baydin et al., 2018), allowing for the computation of
gradients of the likelihood function to be evaluated for all free parameters
in the model.

We estimate the posterior distribution through NumPyro’s No-U-Turn-Sampler
(NUTS). This is an HMC sampler that uses gradient information to find new draw
proposals from the likelihood. Its main advantages over traditional MCMC
methods are it produces chains with small auto-correlation values (e.g. each
draw is independent), requires relatively few warm-up steps, and can draw
samples from very high dimensional distributions without issue.

The priors of the model parameters are as follows:

Cosmology $\displaystyle\text{: }\left\\{\begin{aligned}
\Omega_{\mathrm{M}}&\sim\mathcal{U}(0,1)\\\
\Omega_{k}&\sim\mathcal{U}(-1,1)\\\ w&\sim\mathcal{U}(-2,0)\\\
w_{a}&\sim\mathcal{U}(-3,1)\\\
w_{0}&\sim\mathcal{U}(-2,2)\end{aligned}\right.$ Parent Lens Distribution
$\displaystyle\text{: }\left\\{\begin{aligned} \text{
}\gamma&\sim\mathcal{U}(1.5,2.5)\\\ \text{
}\sigma_{\gamma}&\sim\mathcal{HN}(0.5)\\\ \text{
}\beta&\sim\mathcal{U}(-0.6,1)\\\ \text{
}\sigma_{\beta}&\sim\mathcal{HN}(0.5)\\\ \end{aligned}\right.$ Individual
Lenses $\displaystyle\text{: }\left\\{\begin{aligned} \text{
}\gamma_{i}&\sim\mathcal{U}(\gamma,\sigma_{\gamma})\\\ \text{
}\beta_{i}&\sim\mathcal{U}(\beta,\sigma_{\beta})\\\ \end{aligned}\right.$

Where $\mathcal{U}$ stands for uniform distribution, and $\mathcal{HN}$ is
half normal distribution. Observables like Einstein radius and velocity
dispersion have a Gaussian prior centered on their observed values and have a
standard deviation that comes from the observational uncertainty. For
NumPyro’s NUTS sampler we use "$init\\_to\\_median$" as the initialization
strategy with $target\\_accept\\_prob$ = 0.99. We run 10 chains with 1000
warm-up steps and 1500 sampling steps. Under this set-up, a typical HMC run
takes less than 1 hour on an A100 GPU. The typical $\hat{r}$ for each fit
parameter is smaller than 1.05 with no divergent samples, indicating that all
the chains independently converged on the same values.

##  4 Results

Parameter | $\Omega_{M}$ | $w$ | $\Omega_{k}$ | $\gamma$ | $\sigma_{\gamma}$ | $\beta$ | $\sigma_{\beta}$
---|---|---|---|---|---|---|---
Fiducial | 0.3 | -1 | 0 | 2 | 0.16 | 0.18 | 0.13
$\Lambda$CDM | 0.301 $\pm$ 0.015 | - | - | 2.004 $\pm$ 0.003 | 0.162 $\pm$ 0.003 | 0.176 $\pm$ 0.010 | 0.123 $\pm$ 0.005
$w$CDM | 0.300 $\pm$ 0.015 | -0.98 $\pm$ 0.11 | - | 2.004 $\pm$ 0.003 | 0.161 $\pm$ 0.003 | 0.173 $\pm$ 0.012 | 0.124 $\pm$ 0.005
o$w$CDM | 0.288 $\pm$ 0.094 | -0.98 $\pm$ 0.11 | 0.02 $\pm$ 0.16 | 2.004 $\pm$ 0.004 | 0.161 $\pm$ 0.003 | 0.173 $\pm$ 0.015 | 0.124 $\pm$ 0.005

Figure 4: Posterior distributions and 1D marginalized posterior distributions
of the cosmological parameters and the lens population parameters for 10 000
lens systems. Grey assumes an o$w$CDM universe, red is for $w$CDM, and blue is
for $\Lambda$CDM universe. $\gamma$ and $\sigma_{\gamma}$ are the population
density profile slope mean and standard deviation respectively. $\beta$ and
$\sigma_{\beta}$ are the mean and standard deviation values of the
population’s velocity anisotropy. The contour shows the 68$\%$ and 95$\%$
confidence levels, while the grey dashed line represents the fiducial value
which was used to generate our mock data.

We apply our hierarchical approach to model up to 10 000 lens systems and fit
three different cosmological models. Aside from the $\Lambda$CDM model, we fit
the $w$CDM model, where the equation of state of dark energy is not fixed at
-1. We also fit the o$w$CDM model, which allows the universe to have positive
or negative curvature.

Figure 4 shows the posterior distributions of the $\Lambda$CDM, $w$CDM, and
o$w$CDM models. The marginalized results are also shown in the figure. The
parameters of lens populations are accurately recovered within 68$\%$
confidence level. Especially, the mean of $\gamma$ is constrained with high
precision ($\pm$ 0.03 level). Figure 5 shows the offset between the recovered
and input values of the nuisance parameters in our model, i.e. the
$\gamma_{i}$ and $\beta_{i}$ The standard deviation of the differences are
both 0.12. This number is smaller than the intrinsic scatter of $\gamma$ and
$\beta$ (0.16 and 0.13), and demonstrates that the model is accurately
recovering these nuisance parameters as well as the population and
cosmological parameters that we are primarily interested in.

Figure 5: Histograms of the difference between input and recovered individual
nuisance parameters. $\gamma^{\mathrm{fit}}_{i}$ and
$\beta^{\mathrm{fit}}_{i}$ are the means of the posteriors for each individual
lens. The standard deviation of both fitting errors are comparable than the
intrinsic scatters of the population.

Every model successfully recovers the input cosmological parameters. For
$w$CDM, we find $w=-0.98\pm 0.11$. This is more precise than current
individual cosmological probes: the Pantheon Type Ia supernovae data find
$w=-0.90\pm 0.14$ (Brout et al., 2022a); baryon acoustic oscillations (BAO)
from eBOSS give $w=-0.69\pm 0.15$ (Alam et al., 2021). However, our results
are less constraining on $w$ than the combination of both datasets and cosmic
microwave background data from Planck, which yields $w$ = -1.03 $\pm$ 0.03
(Planck Collaboration et al., 2020).

We find some strong covariances between the cosmological parameters and the
lens population parameters of our hierarchical model. The most noticeable are
covariances between $\langle\beta\rangle$ with $\omega_{M}$, $\Omega_{k}$ and
$w$. This degeneracy arises from the fact that when mass is held constant,
varying $\langle\beta\rangle$ leads to different velocity dispersions and,
consequently, requires different cosmological parameters. Similarly, the mean
density profile slope, $\langle\gamma\rangle$, is covariant with cosmology,
since it governs the mapping from velocity dispersion to Einstein radius.

Figure 6: Top: 1 $\sigma$ uncertainty on the equation of state of dark energy
as a function of numbers of lenses and observational velocity dispersion
measurement uncertainty.
Bottom: At a fixed LOSVD error of 10 km/s, the standard deviation of $w$ as
function of number of lenses. The black dots represents the model result
compared against the solid curve of 10 times the square root of the number of
lenses.

###  4.1 The impact of data quality on the cosmological constraining power of
lens samples

The previous section was based on the assumption that 4MOST delivers 10 000
strong lenses with a velocity dispersion measured to a precision level of 10
km/s. In practice, it is difficult to accurately constrain the velocity
dispersion of the lens galaxy in a spectroscopic survey. First, not all lens
galaxies will be bright enough for their spectra to have sufficient signal-to-
noise ratio (SNR) to measure with 4MOST (Iovino et al., 2023). Second,
measuring velocity dispersions from spectra is intrinsically hard, since we do
not perfectly know the correct choice of stellar templates for the spectral
energy distributions of lens galaxies (Collett et al., 2018; Newman et al.,
2017). Thirdly, continuum contamination from the source galaxy may be
challenging to remove (but see Turner et al 2023, submitted)

In this section, we explore how the constraining power for the $w$CDM
cosmological model varies when we have different sample sizes or velocity
dispersion measurement errors. The results are shown in the top panel of
figure 6. As expected, the constraining power on $w$ improves as we increase
the number of samples or improve the velocity dispersion measurement. The $w$
constraint versus sample size roughly follows $\sigma_{w}=10\times N^{-0.5}$
as shown in the bottom panel of figure 6. For 4MOST observations, there will
be a trade-off between observing more lens systems with relatively low SNR
(resulting in larger LOSVD errors) and obtaining longer exposure time for each
target (resulting in higher SNR but lower sample size). We find that
generally, decreasing the LOSVD error by a factor of two (improving $w$ by a
factor of $\approx$ 1.3-1.7) is comparably helpful for $w$ as increasing the
sample size by a factor of two.

As shown in figure 5, the individual $\beta_{i}$ are constrained comparably
well to direct measurements of the circular velocity curve using HST (0.12 vs
0.05-0.4, Gerhard et al., 2001). The individual $\gamma_{i}$ error (0.12) is
about 5 times larger than the typical result from lens modeling (e.g., 2.11
$\pm$ 0.04 by Dye & Warren (2005), 1.96 $\pm$ 0.02 by Dye et al. (2008), and
2.08 $\pm$ 0.03 by Suyu et al. (2010)). This indicates that adding constraints
on individual $\gamma_{i}$ through lens modeling, instead of leaving them as
nuisance parameters would substantially improve the cosmological constraining
power.

Figure 7: Constraints on a flat cosmological model with evolving $w$ where $w$
= w0 \+ wa(1-a). The grey dashed line is the fiducial value of our mocks. The
blue contour is the posteriors distribution of $w_{a},w_{0}$ with 10 000
strong lensing systems. The green and red lines show current results from Type
Ia supernova and CMB+BAO, respectively (Brout et al., 2022b; Planck
Collaboration et al., 2020).

###  4.2 Evolving Equation of State of Dark Energy

The equation of state of dark energy ($w$) is defined as the ratio of pressure
over energy density. In many theoretical cosmological models, $w$ evolves with
redshift (e.g. Peebles & Ratra (1988); Caldwell (2002); Feng et al. (2005)).
We test our method on an evolving $w$ where $w(z)$ = $w_{0}$ \+
$w_{a}\frac{z}{1+z}$ (Caldwell et al., 1998). Since our mock data has fixed
$w=-1$, the fiducial value this model is $w_{0}=-1$, $w_{a}=0$ (Chevallier &
Polarski, 2001; Linder, 2003).

Our forecast constraints for 10 000 lenses are shown in Figure 7. The 1D
marginalized posteriors give $w_{a}=0.0^{+0.8}_{-1.4}$,
$w_{0}=-1.0^{+0.5}_{-0.4}$. Since these parameters are strongly covariant, we
use the figure of merit (FoM) from equation 6 in Mortonson et al. (2010) to
quantify the overall dark energy constraining power:

$\displaystyle\mathrm{FoM}=\frac{6.17\pi}{A_{95}},$ (10)

Where $A_{95}$ is the Area enclosed within the 95$\%$ confidence contour in
the $w_{0}-w_{a}$ plane. Larger FoM’s therefore imply better dark energy
constraining power. The FoM of our 10 000 lens result is 15. As a comparison,
the 10-year LSST forecast for the combination of weak lensing and Large Scale
Structure which has a FoM of 49 (The LSST Dark Energy Science Collaboration et
al., 2018), for LSST 10-year supernovae the FoM is expected to be 32. We find
that the galaxy-galaxy strong lensing constraints have very similar
$w_{0}-w_{a}$ covariance as those forecast for weak lensing. This is to be
expected since both methods constrain the distance ratio $D_{s}/D_{ls}$.

##  5 The impact of the evolution of Lens Galaxy Population

The results of the previous section were based on mock data where $\gamma$ and
$\beta$ follow Gaussian distributions that do not evolve with redshift.
However, this assumption may not be true of the real Universe. Lens galaxies
might have evolving density slopes due to various reasons. At high redshift,
dissipative processes like gas-rich mergers dominate, leading to steeper total
density slopes (large $\gamma$, see Remus et al. (2013); Sonnenfeld et al.
(2014)). At redshifts lower than $z\sim 2$ (which is all the lenses in our
sample), dissipationless processes like gas-poor mergers can flatten the
density slopes (e.g., Hilz et al. (2012); Hilz et al. (2013)).

Hydrodynamic simulations show that the mass-density slope of early-type
galaxies (ETGs) becomes flatter at lower redshift (Johansson et al., 2012;
Remus et al., 2013; Remus et al., 2017; Wang et al., 2020). However, most
strong lensing studies find no evidence for evolving density slopes with
redshift (Koopmans, 2006; Auger et al., 2010; Ruff et al., 2011; Bolton et
al., 2012; Barnabè et al., 2011; Cui et al., 2017), or a slight decrease in
density slopes with redshift (Sonnenfeld et al., 2013; Cao et al., 2015b; Cao
et al., 2016; Holanda et al., 2017). This might suggest that strong lensing
surveys are biased towards certain lens populations, or that there is no
evolution at all at low redshift. Since the possibility of an evolving mass
profile has not been ruled out, we explore the robustness of our method when
applied to an evolving lens population.

According to the simulations of Remus et al. (2017), the evolving total mass
profile can be approximated by a linear relation: $\gamma_{z}=0.21z+2.03$.
This equation gives unreasonably steep profiles compared to observed lenses,
so for our new mock data we use the evolution of Remus et al. (2017), but fix
the mean and scatter to give a total population that has the same average
slope and scatter as in Section 4. Thus we draw our density slopes from
$\langle\gamma_{z}\rangle=0.21z+1.89$ with a scatter of 0.155.

Figure 8: The impact on the inferred $w$ and $\Omega_{\mathrm{M}}$ of redshift
evolution in the mean lens density profile. The blue contours represent a
model that does not allow for redshift evolution in $\gamma$. The red contours
represent the result for a model which fits for linear evolution of the mean
population $\gamma$. The grey contour results when individual $\gamma_{i}$ are
measured through detailed lens modelling of every lens.

We first try to fit the redshift evolving lens dataset with the non-evolving
hierarchical model used in Section 4. The resulting cosmology posteriors
distributions for 10 000 lens systems are shown in blue in figure 8 (full
result in appendix A). The population parameter of $\langle\gamma\rangle$
population is recovered, but the cosmology parameters and the stellar
anisotropies are systematically incorrect. Thus, a more complex model is
needed to deal with this scenario.

We can generalize our hierarchical model to account for population redshift
evolution by fitting additional parameters that describe the evolution of the
population. In order to keep our parameters as similar as possible to those in
section 4, we parameterise our density profile redshift evolution as follows:

$\gamma_{z}=\langle\gamma\rangle+\Delta_{\gamma}\times(z-0.47)\pm\sigma_{\gamma}.$
(11)

where 0.47 is the mean lens redshift of our population, we set
$\langle\gamma\rangle=2$ as the ensemble mean density profile slope and
$\Delta_{\gamma}=0.21$ as its linear evolution with redshift.

Unlike the model where we ignore the evolution of the population, we find that
fitting for the evolution does reproduce the input lens galaxy population
parameters (Figure 11). We recover the correct redshift density evolution:
$\Delta_{\gamma}=0.195^{+0.015}_{-0.017}$. Using this method, we are able to
constrain $w$ to -1.1 $\pm$ 0.07 level (Figure 8), which is better than the
$\pm$ 0.11 error that we found for $w$CDM without $\gamma$ evolution. For a
flat $w_{0}$$w_{a}$CDM cosmology, the figure of merit improves to 17 as shown
in the top panel of the figure 10. The modest improvements are due to the
slightly reduced intrinsic scatter of our mock population.

Alternatively, the population parameters of describing density profiles can be
constrained on a lens-by-lens basis using detailed lens modelling of the arcs
observed in each system, without the need to know the cosmological model. This
is because the slope of the mass density determines the radial derivative of
the deflection angles, and thus the radial width of the arc is sensitive to
the density profile slope (Dye & Warren, 2005; Suyu et al., 2010; Collett et
al., 2018). With image quality comparable to that of HST and Euclid, one can
constrain $\gamma$ to a precision level of $\sigma_{\gamma}\approx 0.02$ (Meng
et al., 2015), although this neglects the impact of the mass-sheet degeneracy
(Gorenstein et al., 1988; Saha, 2000; Wucknitz, 2002; Liesenborgs & De Rijcke,
2012; Schneider & Sluse, 2013). Whilst lens modelling at scale is currently
challenging, adding a precise prior on $\gamma_{i}$ for each lens should
significantly improve the constraints on cosmological parameters. If we assume
that we are able to pre-determine the value of $\gamma$ for each lens system
with a precision level of 0.02, we can treat $\gamma$ as an observable,
similar to the Einstein radius. The measurement of "$w$" is greatly improved
in this scenario: $w=-1.01\pm 0.06$ for a 10 000 lenses with no evolution of
the population, or $w=-1.08\pm 0.07$ for the population where
$\Delta_{\gamma}=0.21$. Respectively, the figure of merit improves to 28 and
64 for 10 000 lenses assuming a flat $w_{0}$$w_{a}$CDM cosmology (see the
bottom two panels of the figure 10).

##  6 Larger scatter on lens population

In this study, our mock lenses are created, assuming the distribution of
$\sigma_{\gamma}$ and $\sigma_{\beta}$ match those inferred from SLACS data
(Auger et al., 2010) and a set of nearby elliptical galaxies (Gerhard et al.,
2001). This assumption is critical to our forecasts since $\sigma_{\gamma}$
and $\sigma_{\beta}$ quantify how standardizable each lens is. If
$\sigma_{\gamma}$ or $\sigma_{\beta}$ are significantly larger in the real
Universe then the power of this method to constrain cosmography will be
greatly diminished. It should be noted that the lens systems were selected
from ground-based sky surveys, which may introduce biases compared to space-
based surveys like Euclid. Additionally, the nearby elliptical galaxies might
not be a perfect representation of lens elliptical galaxies. For instance, Xu
et al. (2017) analyzed elliptical galaxies in the Illustris simulation and
found that the standard deviation of velocity anisotropy can reach 0.3, while
in our work, we employed a value of 0.13.

We generate mock galaxies with a wider range of $\gamma$ (0.1 - 0.3) and
$\beta$ (0.1 - 0.3) to test the effectiveness of our method. Table 2 presents
the constraints on $w$ ($\sigma_{w}$) with 1,000 lens systems in a flat
universe under different distributions of $\gamma$ and $\beta$. The
constraints deteriorate as the scatter in either $\gamma$ or $\beta$
increases, but even in the most pessimistic scenario the constraints degrade
by a factor of 2.17. On the other hand, our forecasts may be pessimistic if
lenses can be further standardized by understanding the physical properties
that drive the scatter of $\gamma_{i}$ and $\beta_{i}$ away from the
population mean.

| $\mathrm{\sigma}_{\gamma}=0.1$ | $\sigma_{\gamma}=0.2$ | $\sigma_{\gamma}=0.3$
---|---|---|---
$\sigma_{\beta}=0.1$ | 0.90 | 1.03 | 1.33
$\sigma_{\beta}=0.2$ | 1.15 | 1.26 | 1.47
$\sigma_{\beta}=0.3$ | 1.52 | 1.79 | 2.17
Table 2: The relative change of the 68 percent confidence interval of $w$ as a
function of how standardizable lenses are. We assume with 10 000 lens systems,
but change the intrinsic scatter of the lens population $\gamma$ and $\beta$.

##  7 Application on existing data

Parameter | $\Omega_{M}$ | $w$ | $\gamma$ | $\sigma_{\gamma}$ | $\beta$ | $\sigma_{\beta}$
---|---|---|---|---|---|---
$w$CDM | 0.32 $\pm$ 0.25 | -0.96 $\pm$ 0.46 | 1.89 $\pm$ 0.05 | 0.18 $\pm$ 0.03 | -0.17 $\pm$ 0.16 | 0.21 $\pm$ 0.09

Figure 9: Posterior distribution and 1D marginalized distribution of the
cosmological and lens population parameters assuming a $w$CDM universe and 161
real lenses from Cao et al. (2015a) and Chen et al. (2019).

As a simple example of our method, we apply it to existing data to evaluate
how well our power-law mass model can describe real lens galaxies.
Additionally, we can compare our results with previous research that used a
similar method and the $w$CDM model (Cao et al., 2015b; Chen et al., 2019). We
use a dataset of 161 lens systems selected by Chen et al. (2019) (see Table
A1), which were obtained from surveys including LSD, SL2S, SLACS, and S4TM.
They measured the luminosity density slope and calculated the equivalent fiber
radius of the spectrum for each galaxy. As we lack the specific $\delta$ value
for each individual galaxy, we treat $\delta$ as an unknown value and draw
from its measured distribution: $\delta=\mathcal{N}(2.173,0.085)$. The average
velocity dispersion error for this dataset is 22 km/s.

The resulting posterior distribution is shown in figure 9. In our results, we
find that $w$ = -0.90 $\pm$ 0.45, which is an even better measurement compared
to the theoretical result shown in Figure 6. Regarding galaxy population
properties, we obtain $\langle\gamma\rangle=1.89\pm+0.05$,
$\sigma{\gamma}=0.18\pm 0.03$, $\langle\beta\rangle=-0.17\pm 0.16$, and
$\sigma_{\beta}=0.21\pm 0.09$ at a 68$\%$ confidence level. The predicted
$\gamma$ and $\beta$ populations are both smaller than those reported in other
studies. This is likely due to the fact that we lack accurate $\delta$ values.
Additionally, these 161 lens systems were obtained from four different
surveys, and the population (selection function) of lens systems from
different surveys can vary in which case fitting the ensemble with a single
lens population would be incorrect. Understanding the selection function
between different surveys is crucial for accurately measuring the lens
population because the future strong lensing sample is likely to be the
combination of LSST and Euclid discoveries.

##  8 Discussion and conclusion

In this work, we have constructed a hierarchical model of galaxy-galaxy lenses
and the underlying cosmological parameters. We have employed Hamiltonian Monte
Carlo through JAX-based NumPyro modelling to efficiently perform the analysis.
Our findings indicate that we can simultaneously constrain cosmological
parameters and lens galaxy properties under different cosmological models,
including $\Lambda$CDM, $w$CDM, and o$w$CDM. With a sample size of 10 000 lens
systems, our method should achieve a 68$\%$ confidence interval of $\pm$0.11
on $w$. These levels of constraint are comparable to other cosmic probes such
as the Cosmic Microwave Background (CMB), standard candles, weak lensing, and
galaxy clustering. Furthermore, we have shown that the evolution of the lens
population can also be simultaneously constrained. We also tested the ability
of galaxy-galaxy lenses to constrain evolving dark energy. Our forecast Figure
of Merit is 15, which is within a factor of 3 of each individual probe
forecast for the LSST 10-year cosmological constraints (Mortonson et al.,
2010). These results rely only on measurements of the Einstein radius,
redshifts, velocity dispersion, and luminosity profile of each lens.
Additional constraints on the lens density profile from detailed lens
modelling, would improve the cosmological constraining power by a factor of
$\sim$2, however, our method can still work without this huge investment in
detailed modelling.

Additionally, we applied our model to 161 real lenses in Section 7, finding
$w=-0.90\pm 0.45$. As we discuss below, this result is likely systematics
dominated but it illustrates that the method does work on real data.

One of the major simplifications made in this study is that the mass and light
models used may not accurately represent real galaxies. The power law light
profile we employ is a simplification of the more general Sersic profile
(Sérsic, 1963), which can describe the light profile of most elliptical
galaxies. The effect of the PSF is also not considered in our model, which
will bring extra complexity to the velocity dispersion measurement.
Furthermore, galaxies are made of dark and luminous matter, and whilst the
’bulge-halo conspiracy’ gives total density profiles that are close to
powerlaw, the absolute truth is undoubtedly more complex than we have assumed.
Our powerlaw assumptions, make the mathematics of our problem analytic
(Equation 8), but should not hugely change the final constraints, since the
model is only relevant in mapping the aperture dynamical mass onto the mass
within the Einstein radius.

In this study, we made the assumption that all galaxy properties follow
Gaussian distributions, this might not be true in reality. Additionally, due
to a lack of observational evidence, we assumed that there are no correlations
between $\gamma$, $\beta$, and $\delta$. However, in actuality, these
properties are likely to have initial scaling relations. For example, Auger et
al. (2010) found that a steeper mass density profile implies a higher central
surface mass density. Furthermore, Cappellari et al. (2007) discovered that
elliptical galaxies can be categorized as slow and fast rotators, and they
exhibit different $\beta$ populations. There are also some outliers that have
very low $\beta$ found by (Gerhard et al., 2001). In our results, the orbital
velocity anisotropy of lens galaxies exhibits a strong degeneracy with
cosmological parameters. Therefore, understanding these correlations will aid
in better fitting the mass models of individual galaxies. Another potential
observational bias arises from our assumption that elliptical galaxies have a
constant $\beta$, whereas observations suggest that elliptical galaxies have
varying $\beta$ with radius. Specifically, the central region of typical
massive elliptical galaxies tends to be isotropic or mildly radially
anisotropic (Gerhard et al., 2001; Cappellari et al., 2007). In spectroscopic
surveys, the fiber has a fixed radius, which means that the region from which
we measure the velocity dispersion will have a different size depending on the
half-light radius. Consequently, any observed evolution in the properties may
simply be a side effect resulting from the evolution of the lens galaxy’s
angular size. Similarly, we have ignored the potential for redshift evolution
in the population anisotropy, which may be expected from simulations (Xu et
al., 2017).

The simulated lens population in this study was limited to a redshift cutoff
of 1.5. In real observations, the lens population is subject to additional
limitations. For instance, observing the lens galaxy and the source galaxy
simultaneously in a spectroscopic survey might not be feasible for lens
systems with a large Einstein radius. As a result, some low-redshift lens
galaxies may be ruled out from the analysis. Additionally, faint lens galaxies
can lead to large velocity dispersion measurement errors, making them less
likely to be included in the analysis. This can exclude low-mass galaxies and
consequently affect the $\gamma$ population estimation. Conversely, luminous
galaxies beyond the redshift cutoff might be included in the analysis, which
can introduce further complexities. These observational limitations and
selection effects should be taken into account when interpreting the results
and considering the true lens galaxy population.

The method used in this study can also be valuable for the analysis of time-
delay cosmography, such as gravitational lens quasars and supernovae. Similar
to our approach, galaxy density profiles can be simultaneously inferred with
the Hubble constant, as demonstrated in previous work (e.g., Birrer et al.
(2020)). But their mass models are heavily prior dominated, the bias in
$\gamma$ will introduce a significant bias in $\mathrm{H}_{0}$ in LSST era
where the number of lensed quasar/SNe are large (Collett & Cunnington, 2016).
The lens galaxy population parameters we obtained can serve as priors when
modeling gravitational lens quasars or supernovae. However, it is important to
note that normal lens systems may have different population parameters
compared to lensed quasar or supernova systems. Conducting a detailed
investigation into the selection function between these systems will be
necessary to mitigate any biases.

Whilst this work is a simplified investigation of how to simultaneously
constrain the astrophysics of lenses and the underlying cosmological
parameters, it has shown that the potential constraining power is competitive
to more established cosmological probes. These encouraging results motivate
further work on modelling the selection function, fitting more general density
and dynamical profiles, and gathering the required data.

## Acknowledgements

We are grateful to Andy Lundgren, Giovanni Granata, Hannah Turner, Simon
Birrer, Sydney Erickson, Shuo Cao, Shude Mao, Shawn Knabel, Phil Marshall,
Russell Smith, and Leon Koopmans for helpful conversations that have enriched
this work.

Numerical computations were done on the Sciama High Performance Compute (HPC)
cluster which is supported by the ICG, SEPNet, and the University of
Portsmouth.

This work has received funding from the European Research Council (ERC) under
the European Union’s Horizon 2020 research and innovation programme (LensEra:
grant agreement No 945536). TC is funded by the Royal Society through a
University Research Fellowship. For the purpose of open access, the authors
have applied a Creative Commons Attribution (CC BY) license to any Author
Accepted Manuscript version arising.

##  9 Data Availability

Forecast lens population data are available at
https://github.com/tcollett/LensPop. Model posterior chains are available from
the corresponding author on request. The parameters of 161 Lens systems are
available at Chen et al. (2019) in https://doi.org/10.1093/mnras/stz1902.

## References

  * Alam et al. (2021) Alam S., et al., 2021, Phys. Rev. D, 103, 083533
  * Auger et al. (2010) Auger M. W., Treu T., Bolton A. S., Gavazzi R., Koopmans L. V. E., Marshall P. J., Moustakas L. A., Burles S., 2010, ApJ, 724, 511
  * Barnabè et al. (2011) Barnabè M., Czoske O., Koopmans L. V. E., Treu T., Bolton A. S., 2011, MNRAS, 415, 2215
  * Baydin et al. (2018) Baydin A., Pearlmutter B., Radul A., Siskind J., 2018, Journal of Machine Learning Research, 18, 1
  * Biesiada et al. (2010) Biesiada M., Piórkowska A., Malec B., 2010, MNRAS, 406, 1055
  * Bingham et al. (2019) Bingham E., et al., 2019, J. Mach. Learn. Res., 20, 28:1
  * Birrer et al. (2020) Birrer S., et al., 2020, A&A, 643, A165
  * Birrer et al. (2022) Birrer S., Millon M., Sluse D., Shajib A. J., Courbin F., Koopmans L. V. E., Suyu S. H., Treu T., 2022, arXiv e-prints, p. arXiv:2210.10833
  * Bolton et al. (2006) Bolton A. S., Rappaport S., Burles S., 2006, Phys. Rev. D, 74, 061501
  * Bolton et al. (2012) Bolton A. S., et al., 2012, ApJ, 757, 82
  * Brooks et al. (2011) Brooks S., Gelman A., Jones G., Meng X.-L., 2011, Handbook of Markov Chain Monte Carlo. CRC press
  * Brout et al. (2022a) Brout D., et al., 2022a, ApJ, 938, 110
  * Brout et al. (2022b) Brout D., et al., 2022b, ApJ, 938, 110
  * Caldwell (2002) Caldwell R. R., 2002, Physics Letters B, 545, 23
  * Caldwell et al. (1998) Caldwell R. R., Dave R., Steinhardt P. J., 1998, Ap&SS, 261, 303
  * Cao et al. (2015a) Cao S., Biesiada M., Gavazzi R., Pió rkowska A., Zhu Z.-H., 2015a, The Astrophysical Journal, 806, 185
  * Cao et al. (2015b) Cao S., Biesiada M., Gavazzi R., Piórkowska A., Zhu Z.-H., 2015b, ApJ, 806, 185
  * Cao et al. (2016) Cao S., Biesiada M., Yao M., Zhu Z.-H., 2016, MNRAS, 461, 2192
  * Cappellari et al. (2007) Cappellari M., et al., 2007, MNRAS, 379, 418
  * Chae (2007) Chae K.-H., 2007, Astrophys. J. Lett., 658, L71
  * Chae et al. (2002) Chae K.-H., et al., 2002, Phys. Rev. Lett., 89, 151301
  * Chen et al. (2019) Chen Y., Li R., Shu Y., Cao X., 2019, Mon. Not. Roy. Astron. Soc., 488, 3745
  * Chevallier & Polarski (2001) Chevallier M., Polarski D., 2001, International Journal of Modern Physics D, 10, 213
  * Choi et al. (2007) Choi Y.-Y., Park C., Vogeley M. S., 2007, Astrophys. J., 658, 884
  * Collett (2015) Collett T. E., 2015, ApJ, 811, 20
  * Collett & Auger (2014) Collett T. E., Auger M. W., 2014, MNRAS, 443, 969
  * Collett & Cunnington (2016) Collett T. E., Cunnington S. D., 2016, MNRAS, 462, 3255
  * Collett et al. (2018) Collett T. E., et al., 2018, Science, 360, 1342
  * Collett et al. (2023) Collett T. E., et al., 2023, The Messenger, 190, 49
  * Connolly et al. (2010) Connolly A. J., et al., 2010, in Angeli G. Z., Dierickx P., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 7738, Modeling, Systems Engineering, and Project Management for Astronomy IV. p. 77381O, doi:10.1117/12.857819
  * Cui et al. (2017) Cui J.-L., Li H.-L., Zhang X., 2017, Science China Physics, Mechanics, and Astronomy, 60, 80411
  * Duane et al. (1987) Duane S., Kennedy A. D., Pendleton B. J., Roweth D., 1987, Physics Letters B, 195, 216
  * Dye & Warren (2005) Dye S., Warren S. J., 2005, ApJ, 623, 31
  * Dye et al. (2008) Dye S., Evans N. W., Belokurov V., Warren S. J., Hewett P., 2008, MNRAS, 388, 384
  * Einstein (1936) Einstein A., 1936, Science, 84, 506
  * Feng et al. (2005) Feng B., Wang X., Zhang X., 2005, Physics Letters B, 607, 35
  * Freedman et al. (2012) Freedman W. L., Madore B. F., Scowcroft V., Burns C., Monson A., Persson S. E., Seibert M., Rigby J., 2012, ApJ, 758, 24
  * Futamase & Yoshida (2001) Futamase T., Yoshida S., 2001, Prog. Theor. Phys., 105, 887
  * Gavazzi et al. (2008) Gavazzi R., Treu T., Koopmans L. V. E., Bolton A. S., Moustakas L. A., Burles S., Marshall P. J., 2008, ApJ, 677, 1046
  * Gerhard et al. (2001) Gerhard O., Kronawitter A., Saglia R. P., Bender R., 2001, AJ, 121, 1936
  * Gorenstein et al. (1988) Gorenstein M. V., Falco E. E., Shapiro I. I., 1988, ApJ, 327, 693
  * Grillo et al. (2008) Grillo C., Lombardi M., Bertin G., 2008, Astron. Astrophys., 477, 397
  * Hilz et al. (2012) Hilz M., Naab T., Ostriker J. P., Thomas J., Burkert A., Jesseit R., 2012, MNRAS, 425, 3119
  * Hilz et al. (2013) Hilz M., Naab T., Ostriker J. P., 2013, MNRAS, 429, 2924
  * Holanda et al. (2017) Holanda R. F. L., Pereira S. H., Jain D., 2017, MNRAS, 471, 3079
  * Inada et al. (2012) Inada N., et al., 2012, AJ, 143, 119
  * Iovino et al. (2023) Iovino A., et al., 2023, The Messenger, 190, 22
  * Johansson et al. (2012) Johansson P. H., Naab T., Ostriker J. P., 2012, ApJ, 754, 115
  * Koopmans (2006) Koopmans L. V. E., 2006, EAS Publ. Ser., 20, 161
  * Liesenborgs & De Rijcke (2012) Liesenborgs J., De Rijcke S., 2012, MNRAS, 425, 1772
  * Linder (2003) Linder E. V., 2003, Phys. Rev. Lett., 90, 091301
  * Meng et al. (2015) Meng X.-L., Treu T., Agnello A., Auger M. W., Liao K., Marshall P. J., 2015, J. Cosmology Astropart. Phys., 2015, 059
  * Mortonson et al. (2010) Mortonson M. J., Huterer D., Hu W., 2010, Phys. Rev. D, 82, 063004
  * Mortonson et al. (2013) Mortonson M. J., Weinberg D. H., White M., 2013, arXiv e-prints, p. arXiv:1401.0046
  * Myers et al. (2003) Myers S. T., et al., 2003, MNRAS, 341, 1
  * Newman et al. (2017) Newman A. B., Smith R. J., Conroy C., Villaume A., van Dokkum P., 2017, ApJ, 845, 157
  * Oguri & Marshall (2010) Oguri M., Marshall P. J., 2010, MNRAS, 405, 2579
  * Oguri et al. (2008) Oguri M., et al., 2008, Astron. J., 135, 512
  * Peebles & Ratra (1988) Peebles P. J. E., Ratra B., 1988, ApJ, 325, L17
  * Perlmutter et al. (1999) Perlmutter S., et al., 1999, ApJ, 517, 565
  * Phan et al. (2019) Phan D., Pradhan N., Jankowiak M., 2019, arXiv preprint arXiv:1912.11554
  * Planck Collaboration et al. (2014) Planck Collaboration et al., 2014, A&A, 571, A16
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 594, A13
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A6
  * Refsdal (1964) Refsdal S., 1964, Mon. Not. Roy. Astron. Soc., 128, 307
  * Remus et al. (2013) Remus R.-S., Burkert A., Dolag K., Johansson P. H., Naab T., Oser L., Thomas J., 2013, ApJ, 766, 71
  * Remus et al. (2017) Remus R.-S., Dolag K., Naab T., Burkert A., Hirschmann M., Hoffmann T. L., Johansson P. H., 2017, MNRAS, 464, 3742
  * Riess et al. (1998) Riess A. G., et al., 1998, AJ, 116, 1009
  * Riess et al. (2011) Riess A. G., et al., 2011, ApJ, 730, 119
  * Riess et al. (2016) Riess A. G., et al., 2016, ApJ, 826, 56
  * Ruff et al. (2011) Ruff A. J., Gavazzi R., Marshall P. J., Treu T., Auger M. W., Brault F., 2011, ApJ, 727, 96
  * Saha (2000) Saha P., 2000, AJ, 120, 1654
  * Schneider & Sluse (2013) Schneider P., Sluse D., 2013, A&A, 559, A37
  * Schwab et al. (2010) Schwab J., Bolton A. S., Rappaport S. A., 2010, ApJ, 708, 750
  * Sérsic (1963) Sérsic J. L., 1963, Boletin de la Asociacion Argentina de Astronomia La Plata Argentina, 6, 41
  * Sonnenfeld et al. (2013) Sonnenfeld A., Treu T., Gavazzi R., Suyu S. H., Marshall P. J., Auger M. W., Nipoti C., 2013, ApJ, 777, 98
  * Sonnenfeld et al. (2014) Sonnenfeld A., Nipoti C., Treu T., 2014, ApJ, 786, 89
  * Suyu et al. (2010) Suyu S. H., Marshall P. J., Auger M. W., Hilbert S., Blandford R. D., Koopmans L. V. E., Fassnacht C. D., Treu T., 2010, ApJ, 711, 201
  * The LSST Dark Energy Science Collaboration et al. (2018) The LSST Dark Energy Science Collaboration et al., 2018, arXiv e-prints, p. arXiv:1809.01669
  * Treu & Marshall (2016) Treu T., Marshall P. J., 2016, A&ARv, 24, 11
  * Treu et al. (2022) Treu T., Suyu S. H., Marshall P. J., 2022, A&ARv, 30, 8
  * Wang et al. (2020) Wang Y., et al., 2020, MNRAS, 491, 5188
  * Wucknitz (2002) Wucknitz O., 2002, MNRAS, 332, 951
  * Xu et al. (2017) Xu D., Springel V., Sluse D., Schneider P., Sonnenfeld A., Nelson D., Vogelsberger M., Hernquist L., 2017, MNRAS, 469, 1824
  * Zwicky (1937a) Zwicky F., 1937a, Phys. Rev., 51, 290
  * Zwicky (1937b) Zwicky F., 1937b, Phys. Rev., 51, 679
  * de Jong et al. (2012) de Jong R. S., et al., 2012, in McLean I. S., Ramsay S. K., Takami H., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 8446, Ground-based and Airborne Instrumentation for Astronomy IV. p. 84460T (arXiv:1206.6885), doi:10.1117/12.926239

## Appendix A

The appendix contains 2D distributions that are not shown in the main body.

Figure 10: Constraints on $w_{0}$ and $w_{a}$ for a redshift evolving
population mean density profile slope. The top panel results for a model which
fits for linear evolution of the mean population $\gamma$ (FoM=17). The bottom
panel results when individual $\gamma_{i}$ are measured through detailed lens
modelling of every lens (FoM=28). Figure 11: The full 2D posterior
distributions obtained for a redshift evolving population mean density profile
slope. The blue contours result from neglecting redshift dependence of
$\gamma$ whilst the red contours allow $\gamma$ to evolve linearly with
redshift.