Constraints on dark energy from TDCOSMO&SLACS lenses

Problems with the cosmological constant model of dark energy motivate the investigation of alternative scenarios. I make the first measurement of the dark energy equation of state using the hierarchical strong lensing time delay likelihood provided by TDCOSMO. I find that the combination of seven TDCOSMO lenses and 33 SLACS lenses is only able to provide a weak constraint on the dark energy equation of state, $w<-1.75$ at 68% confidence, which nevertheless implies the presence of a phantom dark energy component. When the strong lensing time delay data is combined with a collection of cosmic microwave background, baryon acoustic oscillation and Type Ia supernova data, I find that the equation of state is $w = -1.025\pm 0.029$.


INTRODUCTION
The better part of three decades has passed since one of the most profound discoveries in modern cosmology was made: that the expansion rate of the Universe is currently accelerating (Riess et al. 1998;Perlmutter et al. 1999).In the standard model of cosmology, this acceleration is attributed to the cosmological constant Λ acting as a dark energy, with an equation of state  = / = −1, where  is the pressure and  the energy density of the dark energy.However, the true nature of dark energy, and in particular whether its energy density is actually constant in time, remains a subject of debate (Escamilla et al. 2023).
Various dark energy models have been proposed as alternatives to the cosmological constant, most of which incorporate a dark energy whose density changes over time (Copeland et al. 2006).A modelagnostic way of investigating this dynamical kind of dark energy is to use a parameterisation of the equation of state which allows for deviations from  = −1.A popular choice is the Chevalier-Polarski-Linder (CPL) parameterisation (Chevallier & Polarski 2001;Linder 2003), a first-order Taylor expansion in the scale factor .The parameters  0 and   can thus be constrained with data; deviations from  0 = −1 and   = 0 would indicate that the energy density of dark energy is evolving in time.
Myriad observational data have been used to place constraints on the equation of state of dark energy.For example, the Planck satellite's measurements of the cosmic microwave background (CMB) temperature anisotropies, E-mode polarisation and CMB lensing, along with baryon acoustic oscillations (BAO) measured by the 6dF Galaxy Survey and SDSS (Beutler et al. 2011;Ross et al. 2015;Alam et al. 2017) and the Pantheon catalogue of Type Ia supernovae (SNIa) (Scolnic et al. 2018), produced a measurement of  = −1.028±0.031(Aghanim et al. 2020b).A stated goal of Stage IV ★ E-mail: natalie.hogg@lupm.in2p3.frsurveys such as Euclid is to measure the dark energy equation of state with percent-level precision (Amendola et al. 2018).
Strong lensing time delays are the measurement of the different arrival times of the multiple images of a gravitationally lensed source object (Refsdal 1964).This time delay depends on the angular diameter distances between the objects involved and thus directly probes the expansion rate of the Universe,  0 .A distance ladder -a calibration of luminosity distances, typically to SNIa, using parallax measurements and the period-luminosity relation of Cepheid variable stars -is not needed to obtain cosmological information from strong lensing time delays, making them, on paper, a powerful probe of low-redshift cosmology.
Whilst there is no reliance on the distance ladder in strong lensing time delay measurements of  0 , there can be a fairly significant dependence on how the mass profile of the lens galaxy is modelled (Sonnenfeld 2018).This problem has been recognised since the studies of the very first strongly lensed quasar for which time delays were measured, Q0957+561 (Walsh et al. 1979;Falco et al. 1985).
As an example of the effect of different lens modelling on the measurement of  0 , we can compare the H0LiCOW measurement of  0 = 73.3+1.7  −1.8 kms −1 Mpc −1 (Wong et al. 2020), which comes from six strong lensing time delays and is a 2.4% precision measurement, with the TDCOSMO measurement,  0 = 74.5 +5.6  −6.1 kms −1 Mpc −1 (Birrer et al. 2020), which comes from the original six H0LiCOW lenses plus one new lens and is a 9% precision measurement.
The reduction in precision is mainly the result of relaxing certain assumptions in the H0LiCOW lens modelling, namely making a choice of parameterisation for the lens mass profiles which is maximally degenerate with the so-called mass-sheet degeneracy, which I explain further below.To combat this increase in uncertainty, the TDCOSMO team added 33 additional lenses from the SLACS catalogue in order to provide more information on the mass profiles of the TDCOSMO lenses.This data and its nuisance parameters were combined with the original seven lens dataset in a hierarchical Bayesian manner, leading to a measurement of  0 = 67.4 which is in statistical agreement with both the H0LiCOW measurement and the measurement from the seven TDCOSMO lenses alone.While Lewis & Ibata (2002) studied the theoretical effect of dynamical dark energy on strong lensing time delay measurements of  0 , and the TDCOSMO  0 value itself has been used to constrain variations of the fine structure constant (Colaço et al. 2021), axion-photon couplings (Buen-Abad et al. 2022) and interacting dark energy models (Wang et al. 2022), the full likelihood of Birrer et al. (2020) has not yet been used to obtain cosmological constraints beyond those on  0 and the matter density parameter Ω m .The H0LiCOW likelihood of Wong et al. (2020) was used to obtain constraints on various extensions to ΛCDM, including dynamical dark energy, with the previously alluded to caveat that the mass-sheet degeneracy was only broken by the assumption of certain analytic mass density profiles for the time delay lenses, a potential pitfall which the TDCOSMO likelihood avoids.
In this work, I present the first constraint on the dark energy equation of state using the full hierarchical TDCOSMO likelihood.Strong lensing time delays are particularly useful for this task as the constraints in the  − Ω m plane are typically orthogonal to those coming from standard rulers i.e. the CMB and BAO (Motta et al. 2021).I begin by reviewing strong lensing time delays and the construction of the hierarchical TDCOSMO likelihood.I explain my analysis method and then present my results and conclusions.

Theory
Gravitational lensing is the phenomenon which arises due to the deflection of light by massive objects.The strong lensing regime may be defined as that in which multiple images of a single source are produced.This typically occurs on super-galactic scales, with lenses being galaxies and sources being distant and bright objects such as quasars.
Light from the distant source is deflected by the lens galaxy.Different light paths have different lengths, leading to measurable delays between the arrival times of images.This strong lensing time delay is given by where  od is the redshift of the deflector;  od ,  os and  ds are the angular diameter distances between the observer and deflector, observer and source and deflector and source;  is the observed image position;  is the unknown source position; and () is the lensing potential, which carries the information about the mass density profile of the deflector.In a spatially flat Universe (Ω  = 0), the time delay is inversely proportional to  0 via the angular diameter distances involved, where  () ≡  ()/ 0 is the dimensionless Hubble rate.
The difficulty faced by all strong lensing time delay inference is that under any arbitrary linear re-scaling of the source position  → , image positions  are preserved.The lens model is also accordingly transformed,  → + (1−), where  is the deflection angle.This is known as the internal mass-sheet transform or degeneracy (Falco et al. 1985;Schneider & Sluse 2013, 2014).The only way that such a degeneracy can be broken and a measurement of  0 made is by obtaining direct knowledge either of the absolute source size or of the lensing potential itself.The former may be possible by measuring the size of quasar accretion disks using microlensing, but this technique comes with its own difficulties (Chan et al. 2021), so in order to reliably constrain  0 a choice must be made in how to model the deflector mass density profile.

The TDCOSMO likelihood
Whilst the H0LiCOW team used analytic mass profiles for the deflectors to break the mass-sheet degeneracy, the TDCOSMO work used stellar kinematics data to provide information about the lensing potential.This led to the initial reduction in precision of the  0 measurement compared to the H0LiCOW result.The precision was increased again by the addition of further stellar kinematics data from a set of 33 SLACS lenses, specifically selected for their similarity to the original seven TDCOSMO lenses.Note that these additional lenses do not have time delay information.Furthermore, the combination of data was made under the assumption that the TDCOSMO and SLACS lenses were drawn from the same parent population.
The complete likelihood describing this dataset was constructed hierarchically, meaning that a set of hyperparameters were defined which allowed all constraints related to the mass-sheet degeneracy to be inferred on a population level, whilst the lens and light model parameters,  mass and  light , could be inferred on a lens-by-lens basis.Thus, all remaining uncertainty about the mass-sheet degeneracy is propagated to the level of the  0 inference.
Following Birrer et al. (2020), the posterior distribution of the cosmological parameters of interest, , given  sets of individual lens data D  and the model parameters , is given by The nuisance parameter  is divided into the mass and light model parameters which are constrained at the level of each individual lens, and the set of mass-sheet degeneracy hyperparameters which are constrained at the population level,  pop .The hierarchical TDCOSMO likelihood is thus given by where  is the set of angular diameter distances { od ,  os ,  ds } from which all cosmological results are obtained.For a complete discussion of the construction of the hierarchical likelihood, including the details of the population hyperparameters, I refer the reader to Section 3 of Birrer et al. (2020).
The TDCOSMO measurement of  0 = 67.4+4.1 −3.2 kms −1 Mpc −1 is thus the most precise measurement possible from strong lensing time delays which does not make any assumptions about the deflector mass profiles in order to artificially break the mass-sheet degeneracy.It is also the first strong lensing time delay measurement of  0 which used information from other datasets to improve the precision of the measurement.I will now discuss how I used this hierarchical likelihood to measure the dark energy equation of state.

METHOD
I wrote an external likelihood package for the cosmological modelling and sampling software Cobaya (Torrado & Lewis 2021), so that the hierarchical TDCOSMO likelihood can be used to obtain constraints on cosmological model parameters in combination with any other cosmological likelihood and with any choice of Boltzmann code.This package is publicly available. 1 Using the Markov chain Monte Carlo sampler provided by Cobaya, which is adapted from CosmoMC (Lewis & Bridle 2002;Lewis 2013) and uses a fast-dragging procedure to increase sampling speed (Neal 2005), I obtained constraints on three cosmological models: ΛCDM, CDM and  0   CDM.The CDM model allows for a dark energy with a constant equation of state that may differ from  = −1; the  0   CDM model allows for a dynamical dark energy.In a CDM cosmology (extendable to  0   CDM by the CPL parameterisation shown in Equation 1), and recalling that I keep Ω  = 0, the dimensionless Hubble rate is given by where Ω m = Ω b + Ω c , the sum of the baryon and cold dark matter densities, and Ω DE is the dimensionless energy density of dark energy today.
For each cosmological model I considered, I sampled the posterior distributions of  0 , Ω b ℎ 2 and Ω c ℎ 2 , along with the relevant dark energy equation of state parameter(s).In each case, I also sampled the posterior distributions of the hyperparameters associated with the likelihood, and marginalised over them to obtain the posterior distributions on the cosmological parameters of interest.The computation of the angular diameter distances required by the TDCOSMO likelihood was done using the Boltzmann code CAMB (Lewis et al. 2000;Howlett et al. 2012), but I emphasise that I designed the Cobaya interface so that any theory code currently available in or added to Cobaya in the future can be used for this task.Lastly, I used the Parameterised Post-Friedmann framework in CAMB to allow  to cross −1 (Fang et al. 2008).
I validated my implementation of the likelihood package in Cobaya by comparing my ΛCDM results with the original TDCOSMO results, finding an excellent agreement in terms of constraints on both the cosmological and the nuisance parameters.This validation test, plus all of the code needed to reproduce the results and figures in this paper is also publicly available. 2 Besides the TDCOSMO and SLACS datasets, I also obtained constraints on the models with a dataset consisting of the Planck 2018 measurements of the CMB temperature, polarisation and lensing (Aghanim et al. 2020a,c); the BAO measurements from the 6dF Galaxy Survey (Beutler et al. 2011), the SDSS Main Galaxy Sample (Ross et al. 2015) and the SDSS DR12 consensus catalogue (Alam et al. 2017); the Pantheon catalogue of Type Ia supernovae (Scolnic et al. 2018); and the TDCOSMO + SLACS lenses, in order to compare with the constraints obtained just using the strong lensing time delay data and with those of the Planck collaboration (Aghanim et al. 2020b).I will refer to this combination of data as the "full combination" from now on.
The priors I used for the cosmological parameters are listed in Table 1.Following Birrer et al. (2020), I also used the Pantheon prior on Ω m = N (0.298, 0.022) when using the TDCOSMO and SLACS data alone, though it is important to note that in my analysis Ω m 1 https://github.com/nataliehogg/tdcosmo_ext. 2 https://github.com/nataliehogg/slide.

Parameter
Prior Table 1.Prior ranges of the parameters sampled in my analysis.was not directly sampled since it is treated as a derived parameter in CAMB, the posterior being obtained from those of Ω b and Ω c .In the dynamical dark energy cases, I ensured that acceleration occurs by setting a prior on the dark energy equation of state such that  < − 1 3 .

RESULTS
In this section, I present the constraints obtained on ΛCDM, CDM and  0   CDM using the hierarchical TDCOSMO likelihood and the full combination of data listed above.Given the smaller uncertainty on  0 that comes from combining the TDCOSMO + SLACS lenses, I will only present the TDCOSMO alone results for ΛCDM, to demonstrate the replication of Birrer et al. (2020).For the rest of the results, I will show constraints obtained using the complete TD-COSMO + SLACS dataset of 40 lenses.I used GetDist to make the figures and to compute the marginalised parameter values quoted in the text (Lewis 2019).The dark shaded regions of the contour plots represent the 68% confidence limit and the light shaded regions the 95% confidence limit.

ΛCDM
In Figure 1 alone, whilst the red contours show the constraints from the full dataset of 40 lenses (seven TDCOSMO + 33 SLACS lenses).As expected, this result replicates the findings of the TDCOSMO paper, with  0 = 76.8 +6.4 −5.6 kms −1 Mpc −1 from the TDCOSMO lenses and  0 = 68.7 +3.4  −3.9 from the TDCOSMO + SLACS lenses.These values are fully consistent at 1 with the TDCOSMO result.

𝑤CDM
In Figure 2, I show the one and two-dimensional marginalised posterior distributions of  0 ,  and Ω m in a CDM cosmology from the seven TDCOSMO + 33 SLACS lenses.The likelihood is able to provide a weak constraint on the equation of state of dark energy,  < −1.75 at 68% confidence, which implies that the dark energy is 'phantom', the term applied when  < −1.
Some consideration must be given here to what is, at first glance, a surprising result.It is known that in some datasets, such as Planck 2018, large negative values of  correlate with large positive values of  0 (Escamilla et al. 2023).This is due to the so-called geometrical degeneracy, where  0 , Ω m and  can take various values which in combination lead to the same value for the angular diameter distance to the surface of last scattering and hence the same angular size of the sound horizon, provided the physical sound horizon size is kept fixed (Efstathiou & Bond 1999).
Since strong lensing time delays also rely on angular diameter distances to probe cosmology, I infer that a similar degeneracy exists here.This conclusion is supported by the clear correlation between  0 and  in Figure 2. Thus, the strongly negative dark energy equation of state is likely driven by the high central value of  0 = 78.4+8.3  −6.3 obtained in this case -which is nevertheless consistent with the ΛCDM value at 1.Furthermore, the 95% confidence limit is  < −0.74, reflecting the broadness of the constraint obtained.As expected, based on the results of Birrer et al. (2020), the marginalised posterior value of the matter density Ω m = 0.299±0.022 is largely informed by the Pantheon prior.This also acts to somewhat ameliorate the aforementioned degeneracy.
Lastly, I note that a deeply phantom equation of state (for this value of Ω m ,  ≲ −1.4) corresponds to a violation of the null energy condition (Colgáin & Sheikh-Jabbari 2021), which may be problematic depending on the specific cosmological model (Rubakov 2014).

𝑤 0 𝑤 𝑎 CDM
In Figure 3, I show the one and two-dimensional marginalised posterior distributions of  0 ,  0 ,   and Ω m in a  0   CDM cosmology from the seven TDCOSMO + 33 SLACS lenses.The likelihood is again only able to provide loose upper limits on the dark energy equation of state parameters,  0 < −1.86 and   < 0.102 at 68% confidence, and  0 < −0.861,   < 1.97 at 95% confidence.Again the value of  0 is larger than but still consistent with the ΛCDM measurement:  0 = 79.6 +7.5 −6.0 .In Figure 4, I show the two-dimensional marginalised posterior distributions of  0 and Ω m in a ΛCDM cosmology (red), a CDM cosmology (orange) and a  0   CDM cosmology (yellow) from the seven TDCOSMO + 33 SLACS lenses.From this plot, it is clear that whilst the  0 values in the extended cosmologies are large, they are still consistent at 1 with the ΛCDM result.
It is important to note that my results are very similar to those found by Wong et al. (2020).This similarity implies that, for the six original H0LiCOW lenses, the true lensing potential of each lens is wellapproximated by the analytic profiles used in that work i.e. the masssheet degeneracy is making little contribution to the uncertainty on the cosmological parameters.Nevertheless, this may not be true for every lens in the Universe, and therefore the hierarchical likelihood procedure developed in Birrer et al. (2020) is the one which should be used for future cosmological inference involving strong lensing time delays.

Full combination of data
In Figure 5, I show the two-dimensional marginalised posterior distributions of  0 and Ω m in a ΛCDM cosmology (dark purple), a CDM cosmology (medium purple) and a  0   CDM cosmology (light purple) from the full combination of Planck 2018 + BAO + Pantheon SNIa + TDCOSMO + SLACS data.In this case, I did not use the Pantheon prior on Ω m , since the inclusion of the Pantheon dataset provides the same information as that prior.From this plot we can see that the values of  0 and Ω m in the extended cosmologies are completely consistent with the ΛCDM values at 1 when using the full combination of data; the CDM and  0   CDM constraints are virtually identical.
Furthermore, this combination of data inevitably provides a much stronger constraint on ,  0 and   than the TDCOSMO + SLACS data alone, and removes any hint of a phantom dark energy component.In a CDM cosmology, the dark energy equation of state is measured to be  = −1.025± 0.029, a marginal increase in precision compared to the Planck 2018 + BAO + SNIa value of  = −1.028± 0.031.Both of these measurements are consistent with a cosmological constant.
In a  0   CDM cosmology,  0 = −0.985+0.071 −0.091 and   = −0.18+0.33  −0.25 , compared to  0 = −0.957± 0.08 and   = −0.29 can be clearly seen in Figure 6, where I show the two-dimensional marginalised posterior distributions of  0 and   from the TD-COSMO + SLACS data (red) and from the full combination of Planck 2018 + BAO + Pantheon SNIa + TDCOSMO + SLACS (purple).The dashed lines in this plot show the values of the dark energy equation of state parameters which correspond to a cosmological constant,  0 = −1 and   = 0. Lastly, I computed the Δ 2 ≡  2 −  2 ΛCDM for each case studied.The computed values are shown in Table 2. Whilst the Δ 2 is negative for both extended cosmologies in the TDCOSMO + SLACS cases, implying that they are favoured over ΛCDM, the significance of this decrease in  2 must be evaluated using a difference table due to the greater number of degrees of freedom in the CDM cosmology with respect to ΛCDM.With one additional degree of freedom, and for a 95% level of significance, a |Δ 2 | > 3.841 is required for the improvement in fit to be considered significant.This is clearly not the case here.From these results it is also evident that ΛCDM is a better fit to the full combination of data than the extended cosmologies, since the Δ 2 is on parity with ΛCDM in the CDM case and still does not exceed 3.841 for the  0   CDM case.

CONCLUSIONS
In this work, I presented the first constraints on the equation of state of dark energy from the seven TDCOSMO lenses plus 33 SLACS lenses using the hierarchical likelihood provided by TDCOSMO.I wrote an external likelihood package for the Cobaya software, making this likelihood readily available for public use in combination with other cosmological likelihoods and theory codes.
I replicated the original TDCOSMO results in a ΛCDM cosmology and then explored two extended cosmologies, finding that the TDCOSMO + SLACS data was not able to place strong constraints on the equation of state of dark energy, obtaining only weak upper limits on ,  0 and   .The constraints I obtained all implied the presence of a phantom dark energy component,  < −1, at 68% confidence.
The use of the TDCOSMO likelihood in combination with the Planck 2018 likelihood, BAO data and the Pantheon SNIa catalogue yielded more precise constraints, with all the dark energy equation of state parameters consistent with a cosmological constant,  = −1.I computed the Δ 2 to evaluate the fit of each model, finding no preference for the extended cosmologies over the ΛCDM case.
In conclusion, while strong lensing time delays are beginning to provide a competitive (albeit lens-modelling-dependent) measurement of  0 , it is clear from the results of this work, which, in line with expectations from previous literature, indicate that other probes are more useful when studying dark energy.To improve the strong lensing time delay constraint, a larger dataset is certainly needed; perhaps on the order of hundreds or thousands of lenses (Shiralilou et al. 2020).Furthermore, an increase in precision on the  0 inference will naturally lead to a reduction in the effect of the geometrical degeneracy which is contributing to the hints of phantom dark energy in the TDCOSMO data.Fortunately a number of current or near-future experiments, such as JWST, Roman, LSST and Euclid are likely to provide such data in abundance (Collett 2015).

Table 2 .
The two-dimensional marginalised posterior distributions of  0 and   from the TDCOSMO + SLACS data (red) and from the full combination of Planck 2018 + BAO + Pantheon SNIa + TDCOSMO + SLACS (purple).The dashed lines show the values of these parameters which correspond to a cosmological constant,  0 = −1,   = 0.The  2 , degrees of freedom, and Δ 2 values for each case studied.
from Planck 2018 + BAO + SNIa.These are again both almost equivalent in precision and consistent with a cosmological constant.This