© The Authors 2025

1. Introduction

Allan Sandage famously described his life’s work as the “search for two numbers”: the Hubble constant, H₀, and the deceleration parameter, q₀. Our view of cosmology has changed significantly since then, chiefly with the discovery that the Universe is accelerating (Riess et al. 1998; Perlmutter et al. 1999), the explosion in the number of precise cosmological probes, and the emergence of the concordance Λ cold dark matter (ΛCDM) model (e.g., Frieman et al. 2008, and references therein).

The Hubble constant has become even more central in our quest for understanding the composition and expansion history of the Universe. As the uncertainties in the inferred H₀ have shrunk from a factor of two to a percent level, a tension has emerged between direct measurements based on the present-day Universe, and those based on early-Universe probes, extrapolated to the present day under the assumption of a standard flat ΛCDM cosmology. For example, whereas analysis of the cosmic microwave background (CMB) data from the Planck satellite yields 67.4 ± 0.5 km s⁻¹ Mpc⁻¹ in flat ΛCDM (Planck Collaboration VI 2020), the most recent measurement based on the traditional local distance ladder by the SH0ES team, with Cepheids and Type Ia supernova (SN Ia), yields 73.04 ± 1.04 km s⁻¹ Mpc⁻¹ (Riess et al. 2022). Considerable effort has gone into building alternative local distance ladders based on alternatives to Cepheids (e.g. Freedman et al. 2025; Lee et al. 2025; Li et al. 2024a,b); developing independent methods such as cosmic chronometers (Tomasetti et al. 2023); exploiting early-Universe probes independent of the CMB, including Big Bang nucleosynthesis in combination with baryon acoustic oscillations (BAO, e.g., DESI-Collaboration 2025). Overall, this 8% difference, known as the “Hubble tension”, between early-Universe and late-Universe probes, has crossed the traditional gold-standard 5σ threshold in terms of statistical significance (see, e.g., Di Valentino et al. 2025, for a recent review).

If the Hubble tension is real – and not due to unknown systematic uncertainties in multiple measurements – the implications are profound. From a theoretical standpoint, there is no obvious leading contender to reconcile the measurements in tension. Proposed solutions include for example a modification of the early expansion history of the Universe owing to extra relativistic particles beyond the standard model, or an early dark energy phase, but this is by no means the only possibility (e.g., Knox & Millea 2020; Di Valentino et al. 2021; Gelmini et al. 2021; Abdalla et al. 2022; Schöneberg et al. 2022; Vagnozzi 2023; Di Valentino et al. 2025).

2. Overview of current and past time-delay cosmography analyses

2.1. An independent path to H₀

Time-delay cosmography (Refsdal 1964; Treu & Marshall 2016; Treu et al. 2022, 2024; Birrer et al. 2024; Suyu et al. 2024) provides the opportunity to address the tension in a manner completely independent of all other probes. The delays between the arrival times of multiple images of strongly lensed sources measure absolute distances, and thus yield a measurement of H₀, without the need to rely on the local distance ladder for calibration. Furthermore, typical strong lens systems are composed of a deflector at z ∼ 0.5 and a source at z ∼ 1 − 3. They can therefore be considered a late-time probe, unaffected by pre-recombination physics, and yet measure the Universe on the gigaparsec scale, large enough to be insensitive to local deviations from homogeneity and isotropy (Dalang et al. 2023).

The potential of time-delay cosmography was recognized (Refsdal 1964) even before the first strong lenses were discovered (Walsh et al. 1979). As for many other cosmological probes, however, it took decades of development in methodology and improvement in data quality to reach sufficient precision and accuracy to contribute to the measurement of H₀. Breakthroughs were achieved in the measurement of time delays with high-quality light curves (Fassnacht et al. 2002; Courbin et al. 2005; Millon et al. 2020a; Dux et al. 2025), in the modeling of the potential of the main deflector from the quasar host galaxy extended images (Koopmans et al. 2003; Suyu et al. 2009), and in the reconstruction of the effects of the inhomogeneity of the mass distribution along the line of sight (LoS, see, e.g., Suyu et al. 2010; Greene et al. 2013; Treu & Marshall 2016, for a historical perspective). A turning point was the publication of the analysis of six quadruply imaged quasars (hereafter, quads) with state-of-the-art quality data by Wong et al. (2020). By combining the distances obtained from the systems, Wong et al. (2020) measured $H_0=73.3_{-1.8}^{+1.7}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$ with 2.4% precision. Importantly, apart from the first system (Suyu et al. 2010), the analyses of five out of the six quads were performed blindly to prevent experimenter bias. A seventh quad analyzed in a consistent fashion (Shajib et al. 2020) further improved the precision to 2%, i.e. H₀ = 74.2 ± 1.6 km s⁻¹ Mpc⁻¹ (Millon et al. 2020b, hereafter TDCOSMO-1). This 2% measurement was in excellent agreement with the local distance ladder based on Cepheids and SNe Ia from the SH0ES team (Riess et al. 2022), thus strengthening the statistical tension (Wong et al. 2020; Millon et al. 2020b).

2.2. TDCOSMO and the challenge of mass profile systematics

After the Wong et al. (2020) result, many scientists working on time-delay cosmography joined their effort in a newly formed collaboration called the Time-Delay COSMOgraphy (TDCOSMO) collaboration, with the goal of uncovering and mitigating unknown systematic effects (e.g., Millon et al. 2020b; Gilman et al. 2020; Van de Vyvere et al. 2022; Gomer et al. 2022; Shajib et al. 2022).

The main potential source of systematic uncertainty was found to be the description of the mass density profile of the main deflector, stemming from the lensing mass-sheet degeneracy (MSD; Falco et al. 1985). Given a mass model that reproduces the multiple images, one can obtain a family of models that reproduce the images equally well via simple mathematical transformations (Schneider & Sluse 2013, 2014), which modify the inferred source’s intrinsic position, size, and luminosity, as well as the inferred Hubble constant. The MSD can be broken with non-lensing data or by invoking theoretical assumptions. For example, in its simple form, the MSD corresponds to the existence of an infinite sheet of mass in the plane of the lens. Therefore, asserting that the mass density profile of the deflector goes to zero at infinity following a power law breaks the MSD. Wong et al. (2020), Shajib et al. (2020), and Millon et al. (2020b) broke the MSD by assuming two commonly accepted standard forms for the mass density profile of the deflectors, either a power law (Cappellari 2016) or a composite model with stellar mass following the light and a Navarro et al. (1997, hereafter NFW) profile for the dark matter, and anchoring the mass profiles with stellar velocity dispersion.

Given the high stakes, the TDCOSMO collaboration decided to rebuild the inference of H₀ under much weaker and thus conservative assumptions. Birrer et al. (2020, hereafter, TDCOSMO-4) showed that by parametrizing the mass profile in a way that effectively maximizes the MSD and thus the uncertainty on H₀, the precision drops from 2% to ∼9% for the same seven lenses and datasets of Millon et al. (2020b), finding $H_0=74.5_{-6.1}^{+5.6}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$. TDCOSMO-4 also proposed that, by combining time-delay lenses with more numerous non-time-delay lenses (hereafter external lens samples), one can improve the precision, assuming that the lens galaxies are drawn from the same population. As an example, they combined the seven time-delay lenses with 33 non-time delay lenses from the Sloan Lens ACS Survey (SLACS; Bolton et al. 2006, 2008) sample. The combination improved the precision to 5%. While the central value of the posterior shifted to $H_0=67.4_{-3.2}^{+4.1}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$, it remained consistent with the TDCOSMO-only result, within the errors. The shift could be due to statistical noise, differences between the time-delay and non-time-delay deflectors, or it could indicate that some of the previous assumptions or measurements were incorrect. The information available at that time was insufficient to distinguish the three hypotheses. Higher quality stellar kinematics¹, preferably spatially resolved, is needed to break the mass-anisotropy degeneracy (Courteau et al. 2014), and make progress (Treu & Koopmans 2002; Shajib et al. 2018; Birrer & Treu 2021; Yıldırım et al. 2020, 2023). For example, we now know from comparison with Keck spectra with SNR ∼ 160 (Knabel et al. 2025a) that the measurements of stellar velocity dispersion based on relatively low SNR (∼9 Å⁻¹) SDSS spectra of the SLACS lenses used in TDCOSMO-4 suffer from systematic errors at the level of 3.3% and covariance at the level of 2.3%, insufficient for precision cosmology. Therefore, the TDCOSMO+SLACS measurement in TDCOSMO-4 should be discarded and replaced with one based on better data.

2.3. Advancements in this milestone analysis

In this “milestone” paper, we present a major update to the TDCOSMO collaboration’s inference of H₀, following the methodology introduced in TDCOSMO-4 with hierarchical Bayesian inference. We combine multiple improvements in data and analysis as follows.

1. We include the constraints from an eighth time-delay lens, WGD 2038−4008, that was recently analyzed by our team (Wong et al. 2024), in addition to the seven used in our prior analyses (see Table 1). For convenience, this sample of eight time delay lenses is labeled the TDCOSMO-2025 sample.

2. We use improved methods and template libraries to infer stellar velocity dispersions and estimate residual systematic errors and covariance between galaxies and within spatial bins inside galaxies with spatially resolved kinematics (Knabel et al. 2025b). We apply these new methods to derive more precise and accurate stellar kinematics for the time-delay lenses and external lens samples, either by reanalyzing old data or using new data when available.

3. We use new spectroscopy of time-delay lenses to improve their kinematic measurements. For the lens RX J1131−1231, we use JWST-NIRSpec integral field unit (IFU) spectroscopy (Shajib et al. 2025a) and Keck-KCWI ground-based spectroscopy (Shajib et al. 2023) to derive spatially resolved kinematics. For five systems, we use JWST-NIRSpec spectra to infer integrated stellar velocity dispersions. For the two remaining systems, we use new spectra obtained with the MUSE instrument on the VLT.

4. We adopt constant anisotropy as our default parametrization instead of the Osipkov–Merritt (Osipkov 1979; Merritt 1985, hereafter OM) form used in TDCOSMO-4. This is motivated by studies of local elliptical galaxies (Cappellari 2026, Section 3.3), which show that the anisotropy is nearly constant and never fully radial, as in the outer radii of the OM parametrization.

5. We use two external lens samples with improved lens models. In addition to a reanalysis of the SLACS sample (Tan et al. 2024), we consider the Strong Lenses in the Legacy Survey (SL2S) sample (Ruff et al. 2011; Gavazzi et al. 2012; Sonnenfeld et al. 2013a,b; Sheu et al. 2025), which has a larger redshift overlap with the time-delay deflectors. For the SLACS sample, we use spatially resolved stellar kinematics obtained with KCWI for 14 non-time delay lenses (Knabel et al. 2025a). Contrary to TDCOSMO-4, we do not use the SDSS-based stellar velocity dispersion, since Knabel et al. (2025a) showed that the SDSS spectra of the SLACS galaxies suffer from systematic errors on the stellar velocity dispersion at the level of 3.3% and covariance at the level of 2.3%. For the SL2S sample, we use the re-measured stellar velocity dispersion based on the updated methodology and stellar libraries (Mozumdar et al. 2025). These authors show that the SLACS, SL2S, and time-delay lenses follow the same mass-fundamental plane (Bolton et al. 2008; Auger et al. 2009). We are also using new deeper HST observations of the SL2S lenses presented and modeled by Sheu et al. (2025), including a treatment of lens model systematics as a function of the signal-to-noise ratio (SNR) of the observation. We apply rigorous cuts to ensure that the SLACS and SL2S are matched to TDCOSMO-2025 and to enforce high-quality data, ending up utilizing 11 lenses for SLACS and four for SL2S.

6. We account for departures from spherical symmetry in the deflector population and correct for projection effects and the interpretation of the kinematic observables (Huang et al. 2025).

7. We improve and homogenize our treatment of the LoS effects, extending it to the external lens samples (Wells et al. 2024), including a hierarchical inference of the distribution of external convergence values for our external lenses datasets.

8. We improve our description of the surface brightness profiles of the deflector galaxies, both in the lensing analysis and in the dynamical models. We have redone the fits and replaced the Hernquist models with multi-gaussian expansion decompositions based on double-Sérsic fits. For the SL2S lenses, we use Canada–France–Hawaii Telescope Legacy Survey (CFHTLS; Gwyn 2012) imaging in combination with HST-based models of the lensed features to derive the deflector light profile (Sheu et al. 2025). The CFHTLS images are necessary for cases where the deflector is so red that the HST F475X filter does not capture the deflector light.

As in our previous works, our analysis and cosmographic inference are carried out “blinded” to prevent experimenter bias. This is achieved at the software level by removing the mean of all the cosmographically relevant parameters from posterior probability distribution plots. This was done independently for each dataset, so the collaboration was not able to check for mutual consistency between the TDCOSMO-2025, SLACS, and SL2S datasets prior to unblinding. The posterior was unblinded on May 27, 2025, after the analysis was frozen and approved by internal review, and then the results were published without modification².

The paper is structured as follows: In Section 3, we provide a brief overview of the time-delay cosmography methodology. In Section 4, we describe the sample of lenses that we use for the inference. In Section 5 we describe our new stellar kinematic measurements. In Section 6, we describe the selection criteria that we apply to the SLACS and SL2S samples to ensure high-quality measurements and that they match the properties of the time-delay lenses. In Section 7 we describe our joint hierarchical inference of the cosmological parameters and population parameters. In Section 8 we describe our results in the standard flat Λ cold dark matter (CDM) model. In Section 9 we describe our results in alternative cosmological models, including non-flat ΛCDM, flat wCDM, models with dynamical dark energy, such as w_ϕCDM and w₀w_aCDM, as suggested by recent measurements (DESI-Collaboration 2025). In Section 10, we compare our results to our previous work and to other measurements based on time-delay cosmography. In Section 11 we discuss residual systematic errors and selection effects. In Section 12 we provide a brief summary.

The paper was, for the most part, written prior to unblinding. The work was reviewed by the collaboration prior to unblinding, including four designated reviewers (M.C., J.F., D.S., S.H.S.) and five code reviewers (S.B., A.G., X-Y.H., M.M., A.J.S.). After unblinding, numerical values in plots and tables are automatically updated. Sections 8, 9, and the Appendixes were written after unblinding. The abstract, Section 10, Section 11, and Section 12 were mostly written prior to unblinding and completed after unblinding. The other Sections were only copy-edited for clarity, uniformity, and style during the final readings of the complete paper. The final version of the paper was further reviewed by the collaboration and approved by the final reader (T.T.). The likelihood, data, and scripts are public on the TDCOSMO public repository³.

3. Time-delay cosmography

When the luminosity of a strongly lensed background source, such as an active galactic nucleus (AGN) or a SN, varies over time, the variability pattern manifests in each of the multiple images and is delayed in time due to the different light paths of the different images. The arrival-time difference Δt_AB between two images at observed image-plane positions θ_A and θ_B that originated from the source at source-plane position β, is

$$\begin{array}{c}\hfill \mathrm{\Delta }{t}_{\mathrm{AB}}=\frac{D_{\mathrm{\Delta }t}}{c}[\tau ({\mathit{\theta }}_{\mathrm{A}},\mathit{\beta })-\tau ({\mathit{\theta }}_{\mathrm{B}},\mathit{\beta })],\end{array}$$(1)

where

$$\begin{array}{c}\hfill \tau (\mathit{\theta },\mathit{\beta })\equiv [\frac{{(\mathit{\theta }-\mathit{\beta })}^2}{2}-\varphi (\mathit{\theta })]\end{array}$$(2)

is the Fermat potential (Schneider 1985; Blandford & Narayan 1986) with ϕ(θ) being the lensing potential, and

$$\begin{array}{c}\hfill D_{\mathrm{\Delta }t}\equiv (1+z_{\mathrm{d}})\frac{D_{\mathrm{d}}{D}_{\mathrm{s}}}{D_{\mathrm{ds}}}\end{array}$$(3)

is the time-delay distance (Refsdal 1964; Schneider et al. 1992; Suyu et al. 2010). In this expression, z_d is the redshift of the deflector, while D_d, D_s, and D_ds are the angular diameter distances from the observer to the deflector, from the observer to the source, and from the deflector to the source, respectively.

Knowledge of the lensing potential (and hence Fermat Potential, Eq. (2)) paired with a measured time delay between two arriving images can be turned into a measurement of D_Δt. The Hubble constant is inversely proportional to the absolute distances of objects in the Universe for which redshifts are measured and thus scales with D_Δt as

$$\begin{array}{c}\hfill H_0\propto D_{\mathrm{\Delta }t}^{-1}.\end{array}$$(4)

The inverse proportionality between H₀ and the quantity D_Δt makes H₀ the primary cosmological parameter constrained by the time-delay cosmography.

For details on time-delay cosmography, we refer the reader to the reviews, for example, by Treu & Marshall (2016), Treu et al. (2022, 2024), Birrer et al. (2024).

3.1. Mass-sheet degeneracy

A significant challenge in accurately determining the Fermat potential, and thus D_Δt and H₀, is the mass sheet degeneracy (MSD). The MSD arises from a transformation of the lens system that leaves the observed image positions unchanged while rescaling the source position and altering the inferred D_Δt. Specifically, if a lens model with a (normalized) surface mass density or convergence κ(θ) accurately reproduces the observed lensing configuration, then a family of models κ_λ(θ) described by the following transform will also fit the data:

$$\begin{array}{c}\hfill {\kappa }_{\lambda }(\mathit{\theta })=\lambda \kappa (\mathit{\theta })+(1-\lambda ).\end{array}$$(5)

Eq. (5) above is known as the mass-sheet transform (MST) and is a mathematical transformation where the term (1 − λ)≡κ_MST is formally equivalent to a uniform surface density sheet of convergence (or mass) that extends to infinite angular scales. Hence, the name mass-sheet transform and the name of the associated degeneracy, the mass-sheet degeneracy. The term κ_MST can be positive or negative, since it is defined relative to the average positive density of the Universe. The MST, by means of preserving image positions and being linear, also preserves any higher-order relative differentials of the lens equation. Absolute lensing quantities, such as absolute magnification, source size, and luminosity, however, are not preserved by the MST. Only observables related to either the unlensed apparent source size (β vs. λβ), such as the unlensed apparent brightness, or the lensing potential, are capable of breaking the MSD.

The Fermat potential difference between a pair of images A and B (Eq. 2) scales with λ as

$$\begin{array}{c}\hfill \mathrm{\Delta }{\tau }_{\mathrm{AB},\lambda }=\lambda \mathrm{\Delta }{\tau }_{\mathrm{AB}},\end{array}$$(6)

and so does the relative predicted time delay as

$$\begin{array}{c}\hfill \mathrm{\Delta }{t}_{\mathrm{AB},\lambda }=\lambda \mathrm{\Delta }{t}_{\mathrm{AB}}.\end{array}$$(7)

An MSD effect relative to a specified deflector model might be associated with the mass distribution of the main deflector, referred to as the internal MSD scaled with λ_int, or with inhomogeneities along the LoS of the strong lens system, referred to as the external MSD scaled with the external mass-sheet κ_ext. The quantity κ_ext can be directly inferred for a given LoS using photometric and spectroscopic observables of the LoS environment (e.g., Greene et al. 2013; Wells et al. 2024). The effect of all perturbers along a given LoS that are not included explicitly in the lens model itself is quantified with the external convergenceκ_ext. Conceptually, κ_ext can be thought of as the convergence of a mass sheet, which would produce the same cumulative effect as all LoS perturbers if it were placed coplanar with the lens. For a generic lens, κ_ext is usually of order 10⁻² with positive values corresponding to LoSs that are slightly more dense than the background density of the Universe.

The total MST – the relevant transform to constrain for an accurate Fermat potential determination and H₀ measurement – is the product of the internal and external MST (e.g., Schneider & Sluse 2013; Birrer et al. 2016, 2020)

$$\begin{array}{c}\hfill \lambda =(1-{\kappa }_{\mathrm{ext}})\times {\lambda }_{\mathrm{int}}.\end{array}$$(8)

To summarize, the prediction of the time delay (Eq. 1) can be generalized to

$$\begin{array}{c}\hfill \mathrm{\Delta }{t}_{\mathrm{AB},\lambda }=(1-{\kappa }_{\mathrm{ext}}){\lambda }_{\mathrm{int}}\frac{D_{\mathrm{\Delta }t}}{c}\mathrm{\Delta }{\tau }_{\mathrm{AB}},\end{array}$$(9)

where the Fermat potential Δτ_AB is the Fermat potential inferred without the inclusion of internal and external mass sheet.

When transforming a lens model with an MST, the inference of the time-delay distance (Eq. 3) from a measured time delay and previously inferred Fermat potential transforms as

$$\begin{array}{c}\hfill D_{\mathrm{\Delta }t,\lambda }={\lambda }^{-1}{D}_{\mathrm{\Delta }t}.\end{array}$$(10)

In turn, the Hubble constant, when inferred from the time-delay distance D_Δt, transforms as (from Eq. 4)

$$\begin{array}{c}\hfill H_{0,\lambda }=\lambda H_0.\end{array}$$(11)

3.2. Combining lensing and kinematics

Stellar kinematics is the most prominent and commonly used observable to break the total MSD. The collective motion of stars along the LoS is a direct tracer of the 3D gravitational potential and hence provides an independent mass estimate. Joint lensing+dynamics constraints have been used to provide measurements of galaxy mass profiles (e.g., Grogin & Narayan 1996; Romanowsky & Kochanek 1999; Treu & Koopmans 2002; Koopmans et al. 2003; Koopmans 2004; Barnabè et al. 2011, 2012).

The prediction of the stellar velocity dispersion projected along the line of sight σ_los from any model, regardless of the approach, can be decomposed into a cosmology-dependent and cosmology-independent part, as (see e.g., Birrer et al. 2016, 2019)

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{los}}^2=\frac{1-{\kappa }_{\mathrm{s}}}{1-{\kappa }_{\mathrm{ds}}}\frac{D_{\mathrm{s}}}{D_{\mathrm{ds}}}{c}^{2}J({\mathit{\xi }}_{\mathrm{lens}},{\mathit{\beta }}_{\mathrm{ani}},{\lambda }_{\mathrm{int}}),\end{array}$$(12)

where J is a dimensionless quantity dependent on the deflector model parameters (ξ_lens) and the observational conditions (seeing and spectroscopic aperture) of the velocity dispersion measurement (e.g., Binney & Mamon 1982; Treu & Koopmans 2004; Suyu et al. 2010); c is the speed of light; β_ani characterizes the anisotropy profile of the stellar orbits; and κ_s and κ_ds are the LoS external convergences from the observer to the source and the deflector to the source, respectively.

The internal component λ_int describing a physical mass component that influences the motion of the stars is incorporated into the kinematics modeling term J (Teodori et al. 2022). In the case of a sheet-like perturbation much more extended than the effective radius of the deflector, we can approximate the effect of λ_int as

$$\begin{array}{c}\hfill J({\mathit{\xi }}_{\mathrm{lens}},{\mathit{\beta }}_{\mathrm{ani}},{\lambda }_{\mathrm{int}})\approx {\lambda }_{\mathrm{int}}J({\mathit{\xi }}_{\mathrm{lens}},{\mathit{\beta }}_{\mathrm{ani}}).\end{array}$$(13)

Joint lensing and dynamics constraints are sensitive to the combination of terms present in Eq. (12).

Studies of the LoS environment of the lenses and spatially resolved kinematics are needed to break degeneracies between terms within J. For example, by constraining the angular diameter distance ratio D_s/D_ds from relative distance probes, and the LoS contributions κ_s and κ_ds from number counts compared to simulations, it is possible to infer λ_int. In case of an internal component that is less in extent, detailed radially-binned or 2D kinematics maps can further constrain λ_int.

When combining time delays with lensing and dynamics, one needs to use models that match the lensing configuration and constrain them simultaneously with the observed time delays and stellar kinematics, using Eq. (9) and Eq. (12). These two independent equations can be arbitrarily combined algebraically in 2D angular diameter distance constraints (Birrer et al. 2016, 2019; Yıldırım et al. 2020). A convenient transformation of those constraints is the following:

$$\begin{array}{c}\hfill D_{\mathrm{\Delta }t}=\frac{1}{(1-{\kappa }_{\mathrm{ext}}){\lambda }_{\mathrm{int}}}\frac{c\mathrm{\Delta }{t}_{\mathrm{AB}}}{\mathrm{\Delta }{\tau }_{\mathrm{AB}}}\end{array}$$(14)

and

$$\begin{array}{c}\hfill D_{\mathrm{d}}=\frac{1}{1+z_{\mathrm{d}}}\frac{c\mathrm{\Delta }{t}_{\mathrm{AB}}}{\mathrm{\Delta }{\tau }_{\mathrm{AB}}}\frac{c^{2}J({\mathit{\xi }}_{\mathrm{lens}},{\mathit{\beta }}_{\mathrm{ani}},{\lambda }_{\mathrm{int}})}{{\lambda }_{\mathrm{int}}{\sigma }_{\mathrm{los}}^2}.\end{array}$$(15)

When mapped into the D_Δt–D_d plane as outlined above, the projection on constraints in D_d is invariant under any pure external MSD parameter κ_ext (Jee et al. 2015; Birrer et al. 2019; Yıldırım et al. 2023)⁴. If the approximation of Eq. (13) holds, D_d becomes independent of λ_int.

4. Lens samples

Our sample consists of eight time-delay lenses (hereafter, the TDCOSMO-2025 sample; see Fig. 1 for a gallery montage). Those have measured time delays and cosmography-grade lens models from previous work. In addition, we present new high-quality spatially unresolved stellar velocity dispersions (except RX J1131−1231, which has spatially resolved kinematics). To enhance our constraining power on the internal structure of the deflectors, following Birrer et al. (2020) and Birrer & Treu (2021), we added two external samples of lenses that have kinematics but no time delays. The lens models and kinematic data of non-time delay lenses constrain λ_int, and thus indirectly cosmography when combined with time-delay lenses. These two external samples are the SLACS lenses (Bolton et al. 2008) and the SL2S lenses (Gavazzi et al. 2012). We applied cuts to the two samples of non-time-delay lenses to ensure sufficiently high-quality lensing and kinematic data for cosmography and matched them in properties to the TDCOSMO-2025 sample. After the cuts, the SLACS sample consists of 11 lenses with spatially resolved kinematics from Keck-KCWI (Knabel et al. 2025a) and updated lens models (Tan et al. 2024). After the cuts, the SL2S sample consists of four lenses with updated kinematic measurements (Mozumdar et al. 2025) and lens models (Sheu et al. 2025). The samples are described in this section, and the new spectroscopic measurements are described next in Section 5. Then, Section 6 describes the selection cuts.

Fig. 1. Montage of the eight lensed quasar systems in the TDCOSMO-2025 sample. In each panel, the white bar illustrates the 1″scale. The false-color images are made from two or three bands from the HST, Keck-NIRCam, and Chandra X-ray imaging, given their availability. The image for DES J0408−5354 is adapted from Shajib et al. (2020), HE 0435−1223, B1608+656, and WFI 2033−4723 from Suyu et al. (2017), PG 1115+080 and SDSS J1206+4332 from Wong et al. (2020), RX J1131−1231 from Shajib et al. (2024, credits: NASA/ESA/Chandra), and WGD 2038−4008 from Shajib et al. (2022).

4.1. Sample of time-delay lenses

Our TDCOSMO-2025 sample comprises the six original H0LiCOW⁵ lensed quasars – HE 0435−1223, PG 1115+080, RX J1131−1231, SDSS J1206+4332, B1608+656, and WFI 2033−4723 – as presented by Wong et al. (2020) and the references listed in Table 1. In addition, we included DES J0408−5354, presented by Shajib et al. (2020), and WGD 2038−4008, presented by Wong et al. (2024). The first seven systems are shared with the sample previously published by TDCOSMO-1 and TDCOSMO-4. These eight lenses in the TDCOSMO-2025 sample provide cosmological information through their time-delay distance measurement, D_Δt.

In addition to the new data presented in this paper, the analysis presented here relies on data products published in the literature for the D_Δt measurements. This includes time-delay estimations, lens modeling, lens environment characterization, stellar kinematics, and redshift measurements. Individual papers reporting those measurements for each lens are listed in Table 1.

4.2. Samples of non-time-delay lenses

In this section, we describe the external galaxy-galaxy lens samples with imaging and kinematics data from which we extract information about the mass profiles of the deflector galaxies. We describe the lenses in the SLACS sample in Section 4.2.1 and those in the SL2S sample in Section 4.2.2.

4.2.1. SLACS

The original SLACS sample was selected based on SDSS spectra with multiple redshifts, consistent with a background emission line galaxy being lensed by a foreground massive galaxy, and confirmed by imaging with the Hubble Space Telescope (HST; Bolton et al. 2006, 2008). For our analysis, we selected lenses from the subset of 34 lenses that were uniformly modeled by Tan et al. (2024). These lens models adopted the same power-law plus external shear mass model similar to the TDCOSMO-2025 sample, thus putting the two samples on equal footing in terms of the lens models. The uniform modeling of the SLACS lenses was performed with the automated pipeline DOLPHIN (Shajib et al. 2025b), which uses LENSTRONOMY as the modeling engine (Birrer et al. 2016; Birrer & Amara 2018).

Tan et al. (2024) first selected 50 systems out of the full SLACS sample of 85. This subset was chosen by manually examining the lens images and eliminating systems where: (i) there are nearby satellites or galaxies along the LoS, (ii) the source morphology is extremely complex that would require extensive computational resources to fully capture in the model description, (iii) the deflector is a highly flattened disky galaxy, or (iv) the system lacks archival HST imaging data in the visible band (F555W or F606W). Of these initial 50 systems, the homogeneous and automated modeling performed with DOLPHIN produced successful models for the 34 systems that passed this quality assurance procedure. The effective or half-light radius for each of these systems was measured from large cutouts that covered the full extent of the galaxy light above the background level after subtracting the best-fit model-predicted arcs from it. Further selection cuts to match the TDCOSMO-2025 sample and to ensure cosmography grade quality will be discussed in Section 6.

4.2.2. SL2S

The SL2S sample was selected on the basis of ground-based images showing a red galaxy surrounded by a blue ring (Gavazzi et al. 2012), and confirmed as lenses via space-based HST imaging and ground-based spectroscopy (Ruff et al. 2011; Sonnenfeld et al. 2013a,b).

The original papers did not fit power-law mass density profiles to the SL2S data, which is needed in our analysis. Therefore, we adopted a set of total 34 lenses for which power-law models are available, including 13 from Tan et al. (2024) and 21 from Sheu et al. (2025). Tan et al. (2024) initially modeled 31 systems from the SL2S sample based on the availability of HST imaging, single-aperture velocity dispersion measurement, and redshifts for both the deflector and the source. Out of these 31 attempted systems, the automated modeling with DOLPHIN described in Section 4.2.1 produced models that met the quality criteria for 24 systems. However, some of these HST images have low SNR as a result of being short-exposure images obtained from a Snapshot program. To improve the quality of the lens models, we obtained longer-exposure and higher resolution HST images for 21 systems (HST-GO-17130, PI: T. Treu) in the F475X filter, which Sheu et al. (2025) modeled using LENSTRONOMY. Of these 21 newly modeled systems, 11 of these systems overlap with the Tan et al. (2024) set of 24 lenses. Therefore, we utilized 13 systems modeled by Tan et al. (2024), and 21 systems modeled by Sheu et al. (2025). Further selection cuts to match the TDCOSMO-2025 sample and to ensure cosmography grade quality will be discussed in Section 6.

5. New stellar kinematics measurements

In this section, we briefly describe the kinematic measurements that are included in this work, focusing on those that have changed with respect to TDCOSMO-4. Some of the measurements are presented here for the first time, while the rest were previously introduced in previous studies. For those that are discussed elsewhere, we provide only a brief summary and refer the reader to the original papers.

To achieve precise and accurate stellar velocity dispersions for our sample, we undertook a significant effort to refine our measurement techniques and obtain higher-quality data. From the point of view of measurement techniques, Knabel et al. (2025b) show that for data of sufficient SNR and spectral resolution, sub-percent precision and accuracy can be achieved using appropriately clean libraries of stellar templates and state-of-the-art methods. All of the kinematic measurements used in this analysis are based on this new refined methodology. One important outcome of the Knabel et al. (2025b) methodology is the estimation of systematic errors and covariance between individual objects and within spatial bins of kinematic maps, arising from stellar template libraries. We included those in our analysis for each dataset. We neglected the covariance between spectroscopic observations obtained with independent instrumental setups, as those are shown to be negligible for our purposes (Mozumdar et al. 2025).

We stress that the stellar velocity dispersions measurements presented in this section are vastly superior to previous published values in terms of data quality and analysis methods. The previously published values should be considered superseded, and we recommend that they not be used for cosmography. They are fine for applications that do not require the same level of precision and accuracy, such as galaxy formation and evolution studies (Shajib et al. 2021; Tan et al. 2024; Sheu et al. 2025).

5.1. Stellar kinematics of TDCOSMO-2025 systems

In this section, we describe the kinematic measurements of the TDCOSMO-2025 systems. All of these measurements were performed using the SQUIRREL⁶ pipeline (Shajib et al. 2025a) built on the penalized pixel-fitting (PPXF) software program (Cappellari & Emsellem 2004; Cappellari 2017, 2023). The SQUIRREL pipeline streamlines the process of performing kinematic measurements with multiple stellar template libraries and other model setting combinations to extract the associated systematic uncertainties and covariances, following the methodology of Knabel et al. (2025b).

5.1.1. JWST NIRSpec

JWST NIRSpec (Böker et al. 2022) data are available for six out of the eight TDCOSMO-2025 systems. The data were obtained as part of three programs. RX J1131−1231 was observed in Cycle 1 through program JWST-GO-1794 (PI: S. H. Suyu; Co-PIs: A. Yıldırım, T. Treu). PG 1115+080, HE 0435−1223, and WFI 2033−4723 were observed through program JWST-GTO-1198 (PI: M. Stiavelli). B1608+656 and SDSS J1206+4332 were observed in Cycle 2 through program JWST-GO-2974 (PI: T. Treu; Co-PI: A. J. Shajib). The IFU on NIRSpec was used in all cases centered on the deflector galaxy. The grating/filter setup G140M/F100LP was chosen to include the near-infrared Calcium triplet (hereafter, CaT) at the redshift of the deflector.

Achieving our target accuracy and precision required substantial effort to improve data reduction and calibration with respect to the standard pipeline’s output. We used the custom data reduction pipeline REGALJUMPER⁷ developed by Shajib et al. (2025a) to reduce all the JWST data. REGALJUMPER extends the official JWST data reduction pipeline by including additional cleaning steps for artifacts, cosmic rays (using LA-COSMIC; van Dokkum 2001), and 1/f noise (using NSCLEAN; Rauscher 2024). The version of the official JWST pipeline that we used is 1.17.1, with STCAL version 1.11.1 and JWST Calibration References Data System context 1322. Further steps were performed on the reduced datacube of RX J1131−1231 to correct for the oscillatory pattern (also referred to as “wiggles”) introduced by undersampling (e.g., Law et al. 2023) using the software package RACCOON⁸ (Shajib 2025). The wavelength-dependent line spread function (LSF) was obtained from robust fits to the emission lines of a planetary nebula (Shajib et al. 2025c). The effective PSF for the datacube’s spatial directions at the position of the CaT wavelength was directly estimated by Shajib et al. (2025a) through lens modeling of a 2D image obtained by summing the datacube across a narrow wavelength range encompassing the CaT, while simultaneously reconstructing the PSF. We took the full width at half maximum (FWHM) of a circular Gaussian profile fitted to this reconstructed pixelated PSF, finding a FWHM of 0$\stackrel{″}{.}$148. All these steps are described in detail by Shajib et al. (2025a) in the case of RX J1131−1231.

To obtain the 1D resolved kinematics of RX J1131−1231 used in this analysis, we took six annular bins centered on the lens galaxy with the outermost annulus extending to 1$\stackrel{″}{.}$3 (Fig. 2). As noticeable for the outer annuli, the Hα, [N II], and [S II] emission lines from the quasar host galaxy coincidentally fall close to the CaT lines of the lens galaxy. For that reason, following Shajib et al. (2025a), we simultaneously modeled these emission lines from the quasar host galaxy in our kinematic fitting with PPXF. We followed Shajib et al. (2025a) in modeling spurious spike-like artifacts, which could have astrophysical origins, either from faint sources with undetected continuum at other redshifts along the LoS or the zodiac, or they may be artifacts in the data. However, our tests show that not directly modeling these spurious spikes would change the average velocity dispersion only within a sub-percent level. After an initial fit was performed for the spectra of a given annulus with the model settings described above, we boosted the noise levels to achieve χ_red² = 1 before fitting once again, leading to the uncertainty on the measured kinematics being inflated accordingly. The median SNR of the noise-boosted spectra in the annuli is ≳90 Å⁻¹. We produced the covariance matrix of the measured values that includes modeling systematics by marginalizing over the stellar template library choice between Indo-US (Valdes et al. 2004) and X-shooter Spectral Library (XSL, DR3; Verro et al. 2022) following the methodology presented by Knabel et al. (2025b)⁹, detection thresholds for spurious spikes, combination of additive and multiplicative polynomials, and the fitted wavelength range. We used the “cleaned” versions of the libraries provided by Knabel et al. (2025b). The marginalization across model choice combinations was done using Bayesian Information Criterion (BIC) weighting and Bessel correction. The average magnitude of the off-diagonal terms in the covariance matrix is 0.66%. For the spherical Jeans modeling performed in this analysis, we fit to the v_rms profile shown in Fig. 2 (lower left panel), where v_rms² ≡ v_los² + σ_los² with v_los being the LoS velocity after subtracting off the systemic velocity.

Fig. 2. JWST-NIRSpec spectra and kinematic fits for RX J1131−1231. Top left: The six annuli (black contours), from which summed spectra are extracted, are illustrated on top of the NIRSpec white-light image. The regions around the satellite galaxy and the closest quasar image are excluded. The white bar represents 1″ scale, and the North and East directions are pointed with emerald and yellow arrows, respectively. Bottom left: The measured vrms in the six annuli. The horizontal error bars show the annulus widths, and the vertical error bars show the combined statistical and systematic uncertainty for each measurement. The measured values have 0.66% covariance on average. Second and third columns: The six panels show the integrated spectra in each annulus (gray boxes) and the kinematic fit (red line). The height of the gray box represents the nominal uncertainty levels estimated by the JWST pipeline, and the size of the vertical error bars represents the total boosted uncertainty levels to achieve χred2 = 1 for each fit. The vertical orange lines mark the wavelengths of the Calcium triplet lines from the lens galaxy, and the vertical blue lines mark the wavelengths of the Hα, [N II], and [S II] lines from the quasar host galaxy.

It is important to match the light profile used for the luminosity weighting done in the Jeans modeling to that of the kinematic tracer. For the light profile of the kinematic tracer in RX J1131−1231’s NIRSpec data, we adopted a “stitched” surface brightness profile. Within $R\le 1\stackrel{″}{.}3$ in this stitched profile, we took the surface brightness distribution of the lens galaxy from the double Sérsic fit that was obtained from the lens modeling by Shajib et al. (2025a) that was performed with a 2D image obtained from the NIRSpec datacube by summing within 8700–8800 Å in the lens rest-frame. Outside of this radius (i.e., at $R>1\stackrel{″}{.}3$), we adopted the light profile obtained from the double Sérsic fit performed by Shajib et al. (2023) on a large cutout of the HST F814W imaging (with the pivot wavelength corresponding to 6208 Å in the lens rest-frame) after subtracting off the quasar images and the lens-model-predicted arcs. We perform the stitching as the limited field of view of the NIRSpec datacube (roughly 5″ × 5″) does not provide a robust handle on the light profile that is much further out, which is still relevant in the integration of the 3D light profile performed along the LoS for the kinematic modeling. We choose $R=1\stackrel{″}{.}3$ as the stitching radius, as that is the outermost extent of the measured kinematics from the NIRSpec data.

For the remainder of the lenses with JWST spectroscopy, we present, for the first time, spatially integrated kinematics in this paper. Obtaining spatially resolved maps requires additional work to correct for resampling noise or wiggles (Shajib et al. 2025a), which is beyond the scope of this paper. The maps will be presented in future work and will be used in the next TDCOSMO milestone paper.

We obtained the PSF FWHM corresponding to the CaT’s observed wavelength by extrapolating from the estimated FWHM value for RX J1131−1231 using a wavelength-dependent scaling extracted from simulated PSFs with STPSF (Perrin et al. 2014). Similar to the case of RX J1131−1231, the LSF’s FWHM was obtained at the observed wavelength of the CaT using the values provided by Shajib et al. (2025c). The aperture-integrated spectra and fits are shown in Fig. 3. Table 2 lists these integrated measurements, along with the statistical and systematic uncertainties after marginalizing over choices of stellar template libraries and polynomial degrees following the methodology of Knabel et al. (2025b), with BIC weighting and Bessel correction applied. To estimate the noise for each integrated spectrum, we take the nominal uncertainty levels estimated by the JWST pipeline and perform an initial fit for each of 30 combinations of stellar template libraries and polynomial choices. For each of these fits, we take a noise boosting factor required to achieve χ_red² = 1. We average over these factors for each object’s spectrum and rescale the noise according to that value. The median SNR for the noise-boosted spectra of each of the five lens galaxies ranges within ∼34–59 Å⁻¹. This inflated noise is shown in Fig. 3 and used for all subsequent spectral fittings. All spectra were fitted with PPXF using “cleaned” Indo-US and XSL stellar template libraries in the wavelength ranges 8300–8800 Å. Systematic uncertainties due to stellar template library selection are on average 0.77% when isolated from the effect of the continuum polynomial degrees and wavelength range, as shown by Knabel et al. (2025b). We tested additive (multiplicative) continuum polynomial degrees of 2, 3, 4, 5, and 6 (0, 1, and 2), and for each individual object, we selected the best polynomial degrees as the combinations that minimize the difference between the stellar libraries. Only one object, HE 0435−1223, prefers a multiplicative polynomial degree greater than zero, which is clearly needed by the fit to account for a continuum slope issue likely due to imperfect flux calibration. The effect of varying the additive polynomial degree for a fixed template library, wavelength range, and multiplicative polynomial degree is on average 0.55% for the sample. We also tested the choice of fitted wavelength range by shifting it by 25 Å in either direction from the baseline range (8300–8800 Å), resulting in an average 0.38% systematic effect. Since the change in wavelength range involves a change in the noise realization, and the effect is well below the statistical uncertainties, we did not include the wavelength range as an additional systematic term. To estimate the systematic uncertainties associated with the choice of template library and polynomial degrees following Knabel et al. (2025b), we adopted the options of using BIC weights and Bessel correction. The combined systematic uncertainty from the template libraries and polynomial degrees for the sample as a whole is therefore at a level of 0.95% after adding them in quadrature. The total systematic uncertainty was then added in quadrature from these contributions, which is the value reported in Table 2 for each lens. The covariance of stellar velocity dispersion between galaxies is at the level of 0.4%, which we considered negligible given the levels of statistical and systematic variance.

Fig. 3. JWST-NIRSpec integrated spectra (gray bars) and kinematic fits (red lines) for five TDCOSMO-2025 lenses. The height of the gray boxes represents the nominal uncertainty levels estimated by the JWST pipeline, the size of the vertical error bars represents the total boosted uncertainty levels to achieve χred2 = 1, and the width of the boxes represents the wavelength-pixel size. The x-axis shows the lens rest-frame wavelength. The gray-shaded vertical bands were masked during the fits due to contamination by the lensed quasar features.

5.1.2. Reanalysis of Keck-KCWI data for RX J1131−1231

The Keck-KCWI IFU data for RX J1131−1231 were presented by Shajib et al. (2023). Here, we updated the measurement using the new methodology of Knabel et al. (2025b). Shajib et al. (2023) used only one template library, XSL DR2, in their measurement of kinematics. In this paper, we used XSL DR3 and Indo-US libraries to account for the systematic uncertainty stemming from the template library choice (Knabel et al. 2025b). We also accounted for systematics from other model setting choices adopted by Shajib et al. (2023), namely the wavelength range fitted and the combination of polynomial orders, to produce a covariance matrix that encapsulates all the above sources of systematics. Furthermore, we identified a small systematic bias in Shajib et al. (2023), where the XSL DR2 template library was read out in a manner that assumed uniform wavelength sampling across the library, whereas in reality, the wavelength sampling varies slightly across the templates. Thus, ignoring the slight variations in the wavelength sampling led to a slight underestimation of the stellar velocity dispersion. Overall, the reanalysis of the KCWI performed here resulted in a 1.3% increase in velocity dispersion on average (weighted by the inverse variance of the radial annuli; individual changes ranged from 0.74% to 2.5% from the inner annulus to the outer). In Fig. 4, we illustrate the updated v_rms measurements for RX J1131−1231.

Fig. 4. Measured values of vrms for RX J1131−1231 in radial annuli from the JWST NIRSpec (left-hand panel, the same ones shown in Fig. 2) and Keck KCWI (middle panel). Both of these measurements marginalize over the same choices of template libraries, namely the Indo-US and the XSL DR3, in addition to separate choice combinations for polynomial orders and fitted wavelength range. The arrows in these panels illustrate the size of the PSF FWHM for each case ($0\stackrel{″}{.}15$ for JWST-NIRSpec and $0\stackrel{″}{.}96$ for Keck-KCWI): the higher resolution of the JWST-NIRSpec data enabled the identification of the sharp rise of the velocity dispersion profile in the center. The right-hand panel illustrates the 1D surface brightness profile in optical (blue line) and in the IR (red line). The optical light profile was extracted from double Sérsic fitting by Shajib et al. (2023) from a large cutout that encapsulates the full extent of the galaxy in the HST F814W imaging (which corresponds to a pivot wavelength of 6208 Å in the lens rest-frame). The IR light profile was extracted from the double Sérsic light profile fitted as part of the lens modeling done with a 2D image obtained from the NIRSpec datacube by summing within the wavelength range 8700–8800 Å in the lens rest-frame. The surface brightness profiles are normalized by the corresponding amplitudes at $R=1\stackrel{″}{.}3$ (vertical dashed line). For kinematic modeling of the NIRSpec kinematics, we use a stitched light profile that transitions from the IR light profile at $R\le 1\stackrel{″}{.}3$ to the optical light profile shape at $R>1\stackrel{″}{.}3$.

5.1.3. VLT MUSE

Two TDCOSMO-2025 lenses, namely DES J0408−5354 and WGD 2038−4008, were observed with the Multi Unit Spectroscopic Explorer (MUSE) on the European Southern Observatory’s (ESO) Very Large Telescope (VLT). The observations were carried out as part of two ESO programs: 0102.A-0600(E) (PI: A. Agnello) and 105.205V.002 (PI: C. Ducourant). They consist of 7 × 900s exposures for WGD 2038−4008 and 20 × 675s exposures for DES J0408−5354. Observations were performed in wide-field mode with adaptive-optics correction, covering the wavelength ranges 4700–5803 Å and 5966–9350 Å. Details of the data reduction can be found in Buckley-Geer et al. (2020) for DES J0408−5354. The reduction of WGD 2038−4008 follows Sluse et al. (2019), with improvement on sky subtraction and frame combination as prescribed by Bacon et al. (2023).

Similarly to Section 5.1.1, we focused on obtaining aperture-integrated spectra of the deflector, while the analysis of resolved kinematic maps is left for future work. We summed the flux within circular apertures of 1$\stackrel{″}{.}$5 and 1$\stackrel{″}{.}$0 diameters, centered on the lensing galaxy, for WGD 2038−4008 and DES J0408−5354, respectively. The median SNR of the integrated spectra is ∼38 Å⁻¹ for WGD 2038−4008 and ∼33 Å⁻¹ for DES J0408−5354, sufficient to ensure accurate velocity dispersion measurements (Knabel et al. 2025b). For each object, we removed the quasar components before measuring the velocity dispersion. For WGD 2038−4008, the 2D model of the system was constituted of 4 Moffat components for the quasar images and of a PSF-convolved de Vaucouleurs model for the lensing galaxy (see Sluse et al. 2019, for details). For DES J0408−5354, we used LENSTRONOMY iteratively on each slice of the data cube to model the lensed image. During this process, we froze the mass model to the best model used for the cosmographic analysis in Shajib et al. (2020), while allowing the light profiles of the lens, source, and point-like images to vary for each wavelength.

We measured the aperture velocity dispersions with PPXF, following the methodology introduced by Knabel et al. (2025b) for estimating statistical and systematic uncertainties. We marginalized over the choice of template libraries by combining measurements using Indo-US, XSL, and MILES (Valdes et al. 2004; Sánchez-Blázquez et al. 2006; Falcón-Barroso et al. 2011; Verro et al. 2022) stellar template libraries. The fitting was restricted to the region around the Ca H & K absorption lines – the spectral range used spans 3820–4200 Å in the rest frame of the lens galaxy. The LSF of the MUSE instrument over this spectral range was estimated from the empirical relation described by Bacon et al. (2023). The results of these measurements, including their statistical and systematic uncertainties, are summarized in Table 2. Our measurement for DES J0408−5354 represents a methodological improvement over the previous analysis by Buckley-Geer et al. (2020) and is performed within a slightly different aperture.

Fig. 5. VLT-MUSE spectra and kinematics fits with PPXF for DES J0408−5354 (left-hand panel) and WGD 2038−4008 (right-hand panel). The gray rectangles illustrate the data with the height representing the 1σ uncertainty and the width representing the wavelength-pixel size. The red line illustrates the best-fit model. The principal spectral features probing the kinematics are the Ca H & K absorption lines marked with vertical orange lines. The gray-shaded region on the left-hand panel represents the wavelength range excluded from the fit.

5.2. Stellar kinematics of non-time-delay lenses

5.2.1. SLACS

All of the 34 SLACS systems have single-aperture velocity dispersions measured from the SDSS spectra. The fiber diameter for these measurements was 3″ and the average seeing was 1$\stackrel{″}{.}$4, which were used in TDCOSMO-4. Recent analysis by Knabel et al. (2025a) shows that the SDSS spectra of SLACS lenses have SNR (∼9/Å) insufficient for cosmographic analysis, suffering from systematic errors at the level of 3.3% (see also Bolton et al. 2008) and covariance at the level of 2.3%. Therefore, we did not use the SDSS spectra to obtain the stellar velocity dispersions in this analysis.

New kinematic data were obtained for 14 lenses using optical Keck Cosmic Web Imager (KCWI) integral-field spectroscopy on the W. M. Keck Observatory (Morrissey et al. 2012, 2018). The new data are vastly superior to any previous data obtained for this sample, both in terms of SNR (the integrated SNR of the spectra is 160 Å⁻¹ on average over the sample) and spectral resolution. The data were fully described by Knabel et al. (2025a). The effective parent sample for the SLACS external lens sample was thus reduced to these 14 systems. For this analysis, the 2D kinematic maps were rebinned in radial annuli. The radial annuli were determined in terms of fractions of the lens effective radius, as measured from Sérsic lens light models from Sheu et al. (2025). To avoid mismatched centering of the datacube with respect to the HST photometry, we summed the datacube over the wavelength range in which we fit the kinematics (rest frame 3600–4500 Å) and fit the resultant 2D image with a Sérsic profile, accounting for the average value of $0\stackrel{″}{.}9$ arcsecond seeing (in FWHM) over the observations. The background source arcs are blended with the foreground deflector light due to the limited spatial resolution and PSF blurring. Therefore, for this Sérsic profile fit, we focused on a 5 × 5 spaxel grid centered around the brightest one, where the spaxel size is 0$\stackrel{″}{.}$1457.

To rebin the 2D maps from the Voronoi bins (Cappellari & Copin 2003) to the circular annuli, we took a luminosity-weighted average contribution over the Voronoi bins’ spaxels that lie within the circular annuli. Given Voronoi bins j that have some portion lying within annulus k, the weighted combined v_rms, k is

$$\begin{array}{c}\hfill v_{\mathrm{rms},k}^2=\frac{{\sum }_j{I}_{\mathit{jk}}{v}_{\mathrm{rms},j}^2}{{\sum }_j{I}_{\mathit{jk}}}\end{array}$$(16)

where I_jk is the surface brightness of the portion of Voronoi bin j that contributes to annulus k.

In general, each Voronoi bin contributes to multiple annuli. Any annulus that does not contain a luminosity-weighted Voronoi bin center is folded into the adjacent annulus, which in practice only affects the innermost bin with 0$\stackrel{″}{.}$1 radius. Since the Voronoi bins were constructed from an asymmetric subset of spaxels from the datacube (all those spaxels whose SNR per Å > 1), the outer edges of the kinematic maps are not circular. Therefore, to calculate v_rms of the outermost annulus, we averaged over all of the spaxels beyond the inner edge of the annulus, as we did with the other annuli. Then, we assigned the outer edge of this annulus as the radius at which the area enclosed is equivalent to the total area covered by pixels that contribute to that annulus. From Knabel et al. (2025b), we included systematic uncertainties of 0.67% from template libraries and 0.60% from continuum polynomial degrees to be added in quadrature to the statistical uncertainties, as well as a 0.47% covariance term associated with the template library selection. The profiles are shown in Fig. C.1.

5.2.2. SL2S

Mozumdar et al. (2025) reanalyzed all the available Keck and VLT spectra for SL2S galaxies using the methodology introduced by Knabel et al. (2025b). They estimate the statistical and systematic errors and the covariance between sample galaxies. We refer the reader to Mozumdar et al. (2025) for a description of the data and measurements.

6. External lens sample selection

Before we could include non-time-delay lenses (SL2S and SLACS) in the hierarchical inference to improve our constraints on the population-level mass profile of deflector galaxies – and thus on cosmological parameters – we needed to apply some selection criteria. First, they need to have all the relevant data. Thus, we applied a pre-selection based on the existence of kinematic measurements (Knabel et al. 2025a; Mozumdar et al. 2025), lens models (Tan et al. 2024; Sheu et al. 2025), and LoS measurements (Birrer et al. 2020; Wells et al. 2024). We then have 12 lenses left in the SLACS sample and 21 lenses in the SL2S sample. Second, we needed to ensure that the data quality of the non-time delay lenses is sufficient for cosmological applications. The original SL2S and SLACS data were more than adequate for galaxy formation and evolution studies, and cosmography was not an intended goal, in contrast to the data for the TDCOSMO-2025 lenses that were obtained with the specific goal of cosmography. As described in the previous section, some cosmography-grade data for the non-time-delay lenses were subsequently obtained; however, it is important to apply stringent quality cuts. Third, we needed to match the observable properties of the deflectors in the time-delay and non-time-delay lenses in order to minimize residual differences between thepopulations.

6.1. Lens morphology selection

We applied a light-ellipticity cut of q_light > 0.7, where q_light is the apparent axis ratio of the deflector’s isophotes. This cut is applied to match the range of the q_light values of the TDCOSMO-2025 lenses, among which WFI 2033−4723 has a minimum q_light = 0.73. This cut also removes the more flattened galaxies, which have a higher possibility of being spiral galaxies, whose kinematics are not fully modeled by the spherical Jeans equations. This cut removes 1 out of 12 SLACS lenses and 11 out of 21 SL2S lenses.

6.2. Lensing information selection

We also applied a quality cut based on the lensing information quantity ℐ defined in Eqs. (8–9) of Tan et al. (2024), which is a weighted SNR of the lensed arcs from the extended source, with the lens light subtracted. The purpose of this cut was to ensure that the data contained enough information to constrain the power-law slope γ. There are two free dimensionless parameters a and b in the weights, and Sheu et al. (2025) adopted values a = 2 and b = 0.1 that minimize the correlation between logℐ and log σ_γ (σ_γ is the statistical uncertainty of the power-law slope γ). In this work, we also used these values for all SLACS and SL2S lenses. We selected lenses with higher lensing information. For lenses with single-aperture spectroscopy, we applied a lower-end cut ℐ > 150. For lenses with KCWI spatially resolved spectroscopy, which directly constrains the mass density profile slope to better precision than lensing, we did not use the lensing-derived slope. Thus, we did not apply the lensing information cut. As a result, this cut removes 5 out of 21 SL2S lenses, with 2 already removed by the lens morphology selection.

6.3. Kinematics data quality selection

To avoid potential biases in the stellar velocity dispersion and mitigate systematics, we applied a quality cut on the stellar kinematics of the SL2S lenses. We kept only those based on spectra with SNR > 13 Å⁻¹, for which systematic errors are below 2% (Mozumdar et al. 2025). This cut removes 5 out of 21 SL2S lenses, with 2 already removed by the lens morphology selection and 1 already removed by the lensing information selection. In practice, the four final selected systems have spectra with SNR greater than 18 Å⁻¹.

6.4. Velocity dispersion selection

Velocity dispersion is a quantity tightly correlated with mass and all other observables in elliptical galaxies (Auger et al. 2010; Cappellari 2016). For galaxy-scale lenses, measured stellar velocity dispersion σ_los and the lensing velocity dispersion σ_SIS derived from the Einstein radius and cosmological parameter¹⁰ are the same within 7% intrinsic scatter (Treu et al. 2006; Treu 2010). To match the SLACS and SL2S lenses with the TDCOSMO-2025 lenses, we selected subsamples based on σ_SIS ∈ [150, 350] km/s, the range spanned by our time-delay lens sample. This cut does not remove any of the SLACS lenses and removes 1 out of 21 SL2S lenses.

6.5. Environment selection

We applied selection cuts to the local environment of the deflector galaxies, as we require their lensing effects to be described by the composition of a power-law lens mass model and an external, uncorrelated LoS contribution. Therefore, we specifically excluded lenses in an overdense region or those with a perturber nearby.

The local over(under)-density is characterized by the relative over-density quantity ζ_1/R. For the SLACS sample, Birrer et al. (2020) used the r-band data in the DESI Legacy Imaging Surveys (Dey et al. 2019) to count objects with 18 < r < 23 near the lens galaxies within 3″ to 120″ annuli. The perturbative effect of these nearby objects is quantified by the number N_1/R defined as (Greene et al. 2013)

$$\begin{array}{c}\hfill N_{1/R}={\displaystyle \sum _{i,R<2\prime }}\frac{1}{R_i},\end{array}$$(17)

where R_i is the inverse projected distance from the nearby i-th object to the lens galaxy, and the summation was performed within the 3″ to 120″ annulus. To compare with random fields, N_1/R was also calculated for 10,000 random points in the Legacy Imaging Surveys footprint. The relative over-density ζ_1/R is

$$\begin{array}{c}\hfill {\zeta }_{1/R}=\frac{N_{1/R}}{{⟨N_{1/R}⟩}_{\mathrm{rand}}},\end{array}$$(18)

where the denominator is the median of the N_1/R of the random LoSs used for calibration. We refer the reader to Section 6.3.3 of Birrer et al. (2020) for more details on the LoS characteristics of the SLACS sample.

For the SL2S sample, we utilized the LoS measurement of the LoS summary statistics from the CFHTLS, using a magnitude cut of i < 24 in the i-band (Wells et al. 2024). The inverse distances 1/R of nearby objects are calculated within an aperture of 5″ to 120″ for both the SL2S lenses and the reference fields. The over-density ζ_1/R is calculated as

$$\begin{array}{c}\hfill {\zeta }_{1/R}=\mathrm{median}\left[\frac{N_{\mathrm{lens}}^{R<2\prime }\times \mathrm{median}[1/R_{\mathrm{lens}}]}{N_{\mathrm{ref}}^{R<2\prime }\times \mathrm{median}[1/R_{\mathrm{ref}}]}\right],\end{array}$$(19)

where $N_{\mathrm{lens}}^{R<2\prime }$ and $N_{\mathrm{ref}}^{R<2\prime }$ denote the total number of objects within the 120″ aperture under the magnitude cut for the lens and reference field points, and the median in the numerator and the denominator is taken over the aperture. Note that the definition of ζ_1/R is different for the SLACS and SL2S sample, as by visual inspection of ζ_1/R versus κ_ext we find that using the expression in Eq. (18) for the SLACS sample leads to positive correlation between both variables as expected, while for the SL2S sample we do not see the correlation except when using the expression for ζ_1/R in Eq. (19).

For the SLACS sample, we applied an over-density cut ζ_1/R < 2.10, which is the same cut as in Birrer et al. (2020). However, no lens was rejected by this selection cut. For the SL2S sample, we did not apply selection cut on ζ_1/R for two reasons: (i) all the values of the relative over-density ζ_1/R (defined as in Eq. (19)) for our SL2S sample (Wells et al. 2024) do not exceed 2.10, and (ii) the selected SL2S sample lenses are not in overdense regions, as shown in Fig. 6. Fig. 6 also shows the κ_ext distributions for the selected SLACS sample. The population distribution of κ_ext is consistent between our selected SLACS and SL2S samples, although this is not required since we marginalized over κ_ext for individual lenses in the inference.

Fig. 6. Probability distribution function of external convergence κext for the quality samples of the SLACS and the SL2S lenses. Dashed lines are κext of individual lenses, and solid lines are their combination. The median of κext for both samples is shown in the figure legend.

6.6. Selected SLACS and SL2S lenses

We present a summary of the selected SLACS and SL2S lenses used in our joint inference on cosmology and mass density profiles, along with our sample of TDCOSMO-2025 time-delay lenses. After the cuts, we end up with 11 lenses from the SLACS sample and four lenses from the SL2S sample. We reiterate that the SLACS sample is significantly smaller than the one used by Birrer et al. (2020), but it is vastly superior in terms of quality and information content, having replaced the SDSS stellar velocity dispersion with Keck-KCWI spatially resolved measurements, and having updated the lens models (Tan et al. 2024). Their properties are listed in Table A.1.

The key physical quantities of the deflectors together with the source redshifts are also shown in Fig. 7. We compare these quantities of the SLACS and SL2S lenses with those of the TDCOSMO-2025 lenses. Even though the samples are small and shot noise prevents strong statistical conclusions, apparent differences from sample to sample include: (i) the distribution of lens and source redshifts between the SLACS and the SL2S sample, as a result of their selection processes (Bolton et al. 2006; Gavazzi et al. 2012) – the SL2S being better matched to the TDCOSMO-2025; (ii) the half-light radii R_eff are generally smaller for SL2S and TDCOSMO-2025 than for SLACS, consistent with the evolution of the size–mass relation (van Dokkum et al. 2010). Except for these differences, the other properties of the SLACS and the SL2S lenses agree well with the TDCOSMO-2025. We note also that the SL2S, SLACS, and TDCOSMO-2025 lenses follow the same lens-mass fundamental plane, as recently shown by Mozumdar et al. (2025), lending further support to the assumption that they are drawn from the same parent population.

Fig. 7. Properties of selected lens samples that enter the cosmology inference. The sample selection is only based on (i) lens and source existing spectroscopic redshift measurements; (ii) lensing analysis based on imaging data; (iii) LoS analysis. From left to right, the quantities are: lens redshift zlens, background source redshift zsource, stellar velocity dispersion of the deflector σlosap, half-light radius of the deflector Reff, Einstein radius θE, Reff/θE, lensing power-law slope γpl, lensing stellar velocity dispersion σSIS, relative over-density ζ1/R, apparent axis ratio of light profile of the deflector qlight, and lensing information ℐ calculated under the choice of scaling parameters from Sheu et al. (2025). We demonstrate for the selected lenses, there is no significant difference from sample to sample in the parameter spaces of this plot, except for the distribution in redshifts and the difference in effective radius expected from the evolution of the size-mass relation (van Dokkum et al. 2010), see Section 6.6 for more information.

Next in Section 7, we describe how we combined the quality-assured sample of non-time-delay lenses with the TDCOSMO-2025 lenses to obtain better constraints on the mass density slope and ultimately H₀.

7. Hierarchical inference

7.1. General concept

For the hierarchical analysis, we follow the method presented by Birrer et al. (2020). We provide a brief overview of the method here, highlighting the detailed parameterizations, assumptions, and priors that inform our inference. The reader is referred to Birrer et al. (2020) for more details.

In Bayesian language, we want to calculate the probability of the cosmological parameters π, given the strong lensing data set p(π|{𝒟_i}_N), where 𝒟_i is the dataset of an individual lens (including imaging data, time-delay measurements, kinematic observations and the properties of the LoS galaxies) and N is the total number of lenses in the sample.

In addition to π, we denote all the model parameters – including both those at the single lens level and those at the population level – with ξ. Using Bayes’ rule and considering that the data of each individual lens 𝒟_i are independent, we can write:

$$\begin{array}{cc}\hfill p(\mathit{\pi }\mid {\{{\mathcal{D}}_i\}}_N)\propto \mathcal{L}({\{{\mathcal{D}}_i\}}_N|\mathit{\pi })& p(\mathit{\pi })={\displaystyle \int \mathcal{L}}({\{{\mathcal{D}}_i\}}_N\mid \mathit{\pi },\mathit{\xi })p(\mathit{\pi },\mathit{\xi })d\mathit{\xi }\hfill \\ \hfill & ={\displaystyle \int \prod _i^N}\mathcal{L}({\mathcal{D}}_i|\mathit{\pi },\mathit{\xi })p(\mathit{\pi },\mathit{\xi })d\mathit{\xi }.\hfill \end{array}$$(20)

In the following, we divide the nuisance parameter, ξ, into a subset of parameters ξ_i that we constrained independently for each lens, and a set of parameters ξ_pop that require to be sampled across the lens sample population globally. The parameters of each individual lens ξ_i include the lens model, source, and lens light surface brightness and any other relevant parameter of the model needed to predict the imaging data. Hence, we can express the hierarchical inference (Eq. 20) as

$$\begin{array}{cc}\hfill p(\mathit{\pi }\mid {\{{\mathcal{D}}_i\}}_N)\propto {\displaystyle \int \prod _i}& \left[\mathcal{L}({\mathcal{D}}_i\mid D_{\mathrm{d},\mathrm{s},\mathrm{ds}}(\mathit{\pi }),{\mathit{\xi }}_i,{\mathit{\xi }}_{\mathrm{pop}})p({\mathit{\xi }}_i)\right]\hfill \\ \hfill & \times \frac{p(\mathit{\pi },{\{{\mathit{\xi }}_i\}}_N,{\mathit{\xi }}_{\mathrm{pop}})}{{\displaystyle {\prod }_i}p({\mathit{\xi }}_i)}d{\mathit{\xi }}_{\{i\}}d{\mathit{\xi }}_{\mathrm{pop}},\hfill \end{array}$$(21)

where {ξ_i}_N = {ξ₁, ξ₂, ..., ξ_N} is the set of parameters describing the individual lenses and p(ξ_i) are the interim priors on the model parameters in the inference of an individual lens. The cosmological parameters π are fully encompassed in the set of angular diameter distances {D_d, D_s, D_ds}≡D_d, s, ds, and thus, instead of stating π in Eq. 21, we now state D_d, s, ds(π). Up to this point, no approximation was applied to the full hierarchical expression (Eq. 20).

To infer H₀ from a set of lenses, some priors describing the population hyperparameters and their correlations need to be imposed. Our choices and parameterization are presented in Section 7.2.

7.2. Deflector population model

In the following, we describe the model components and their hierarchical treatment for the deflector population. We describe the deflector mass density profiles in Section 7.2.1, the stellar anisotropy and kinematics description in Section 7.2.2, and the LoS convergence in Section 7.2.3.

7.2.1. Mass density profiles

For the mass density profile assumptions, we followed closely the approach of Birrer et al. (2020). The deflectors in the TDCOSMO-2025 lenses are massive elliptical galaxies. These galaxies, observationally, follow a tight relation in a parameter space of luminosity, size, and velocity dispersion (e.g., Faber & Jackson 1976; Auger et al. 2010; Cappellari 2016; Bernardi et al. 2020; Mozumdar et al. 2025), exhibiting a high degree of self-similarity among the population.

For the deflector mass profile, to maximize the uncertainties allowed by the MSD, we chose not to break the MSD based on imaging data and the assumption of a mass profile as in Wong et al. (2020). Instead, we used stellar kinematics to break the MSD and constrain the mass profile. To do so, we chose as a baseline model a power-law elliptical mass distribution (PEMD, Kormann et al. 1994; Barkana 1998)

$$\begin{array}{c}\hfill \kappa ({\theta }_1,{\theta }_2)=\frac{3-{\gamma }_{\mathrm{pl}}}{2}{\left[\frac{{\theta }_{\mathrm{E}}}{\sqrt{q_{\mathrm{m}}{\theta }_1^2+{\theta }_2^2/q_{\mathrm{m}}}}\right]}^{{\gamma }_{\mathrm{pl}}-1},\end{array}$$(22)

where γ_pl is the logarithmic slope of the profile, q_m is the axis ratio of the minor and the major axes of the elliptical profile, and θ_E is the Einstein radius. In the models, we added an additional external shear component. This model is constrained on the lens-by-lens case based on high-resolution imaging data. The PEMD lens profile inherently breaks the MSD, and the parameters of the PEMD profile can be precisely constrained (within a few per cent) by exquisite imaging data. In this work, we avoided describing the PEMD parameters at the population level, such as redshift, mass, or galaxy environment, and utilized the individual lens inference posterior products derived from the data based on uniform priors. To allow for the full flexibility in the radial density profiles that can be produced by an MST, we added a global internal MST distribution specified at the population level. We note that, physically, λ_int does not imply that galaxies contain an infinite sheet of mass. The MST works as long as the “sheet” tapers off sufficiently far from the images (Blum et al. 2020; Birrer et al. 2020). This parameter should be viewed as a parametrization of the deviations from a power-law profile, which is particularly powerful because it is directly degenerate with H₀ and therefore allows us to capture the full uncertainty.

We note that the power-law slope γ_pl of the PEMD profile inferred from imaging data is a local quantity at the Einstein radius of the deflector. The Einstein radius is a geometrical quantity that depends on the mass of the deflector and the redshifts of the lens and the source. Thus, the physical location of the measured γ_pl from imaging data depends on the redshift configuration of the lens system. In a scenario where the mass profiles of massive elliptical galaxies deviate from an MST of a PEMD, resulting in a gradient in the measured slope γ_pl as a function of physical projected distance, a global joint MST correction on top of the individually inferred PEMD profiles may lead to inaccuracies.

To address such potential inaccuracies, we allowed for a radial trend in the applied MST relative to the local quantities inferred from imaging. This was achieved by parameterizing the global MST population with a linear relation in R_eff/θ_E as

$$\begin{array}{c}\hfill {\lambda }_{\mathrm{int}}(R_{\mathrm{eff}}/{\theta }_{\mathrm{E}})={\lambda }_{\mathrm{int},0}+{\alpha }_{\lambda }(\frac{R_{\mathrm{eff}}}{{\theta }_{\mathrm{E}}}-1),\end{array}$$(23)

where λ_int, 0 is the global MST parameter when the Einstein radius is at the half-light radius of the deflector, i.e., R_eff/θ_E = 1, and α_λ is the linear slope in the expected MST parameter as a function of R_eff/θ_E (Birrer et al. 2020). In this form, we assumed self-similarity in the lenses with regard to their half-light radii. In addition to the global MST normalization and trend parameterization, we added a Gaussian distributed intrinsic scatter with standard deviation σ(λ_int) at fixed R_eff/θ_E, to allow for a degree of non-uniformity between the deflectors.

For RX J1131−1231 and for the SLACS lenses, for which we have IFU kinematics data, we also sample individually the power-law slope γ_pl for each of these lenses. The IFU data is constraining enough to supersede the imaging-inferred γ_pl values for the SLACS lenses. For RX J1131−1231 we include the full lensing information and the covariance between inferred D_Δt and γ_pl. Details about the procedure and analysis of RX J1131−1231 can be found in Appendix B.

7.2.2. Stellar anisotropy and kinematics

Strong lensing achieves more constraining power when combined with dynamical mass measurements from kinematics (Treu & Koopmans 2002; Auger et al. 2010; Shajib et al. 2018; Yıldırım et al. 2020, 2023). However, LoS velocity dispersions from single-slit spectroscopy cannot constrain the anisotropy of stellar orbits; spatially-resolved kinematics are required (Cappellari 2008; Barnabè et al. 2009; Shajib et al. 2018; Birrer et al. 2020; Yıldırım et al. 2020; Birrer & Treu 2021; Shajib et al. 2023) to break the degeneracy between the dynamical mass profile and the anisotropy, the mass–anisotropy degeneracy (MAD) (Binney & Mamon 1982; Gerhard 1993; Courteau et al. 2014).

After the early pioneering works (e.g., TIGER, Bacon et al. 1995), at the beginning of the century, spatially resolved kinematics studies of elliptical galaxies started being carried out for statistically significant samples of nearby galaxies (e.g., the SAURON survey: Bacon et al. 2001; de Zeeuw et al. 2002; Emsellem et al. 2004, 2007; Cappellari et al. 2006, 2007). Following the success of these early experiments, integral-field spectroscopy (IFS) became the norm in most major observatories (e.g., MUSE, Bacon et al. 2010). The SAURON and following IFU surveys (ATLAS^3D, Cappellari et al. 2011; CALIFA, Sánchez et al. 2012; SAMI, Bryant et al. 2015; SDSS MaNGA, Bundy et al. 2015) have developed consistent pictures of the kinematic properties of elliptical galaxies in the local Universe (Cappellari 2026). At the same time, methods and systematic uncertainties in the extraction of stellar kinematics at the spaxel level using high-resolution stellar templates have become standardized and well understood (Knabel et al. 2025b), using, for example, PPXF.

Dynamical models have been constructed with the spatially resolved kinematics of the aforementioned IFU surveys by solving 2-integral (Emsellem et al. 1994) and 3-integral Jeans equations under assumptions of spherical or axial symmetry using the Jeans Anisotropic Modeling (JAM) method (Cappellari 2008, 2020) and more generally with Schwarzschild models of orbital superposition (Schwarzschild 1979; Richstone & Tremaine 1988; Rix et al. 1997; van der Marel et al. 1998; Cappellari et al. 2006; Thomas et al. 2004; Lipka & Thomas 2021; Neureiter et al. 2021; Pilawa et al. 2024). Despite some intrinsic model degeneracies, these studies have shown that one can reasonably well constrain the orbital anisotropy of nearby elliptical galaxies using the 2D projection of stellar LoS velocity distributions (LoSVD) and the projected tracer population, i.e., the surface brightness profile (see Cappellari 2026, for a full review).

The velocity anisotropy at a 3D location in the galaxy can be parameterized by β ≡ 1 − σ_t²/σ_r² that describes the ratio of tangential to radial velocity dispersions, where β = 0 corresponds to the isotropic case. Under the assumption of axisymmetry, the velocity ellipsoid (in velocity space) can be assumed to align with the symmetry axis and described in cylindrical coordinates so that β_z = 1 − σ_z²/σ_R². This was found to describe well the stellar kinematics of regularly rotating elliptical galaxies (e.g., Cappellari 2016, Section 3.4). The other extreme assumption for the alignment of the velocity ellipsoid is radial alignment from the center of the galaxy. In this case, the anisotropy can be described in spherical coordinates as β_r = 1 − σ_θ²/σ_r². In the spherical limit, these models converge to spherical models. This is appropriate for the mildly triaxial family of “pressure-supported” non-regular rotators, which are always close to spherical. Spatially resolved kinematics of nearby elliptical galaxies have been shown to be well-fit by a single spatially-uniform nearly-isotropic anisotropy (slightly radial, i.e., β > 0) with the assumption of an oblate axisymmetric geometry. These models are generally preferred over models with radially variant anisotropy profiles (e.g., the OM profile). A general agreement has been demonstrated between axisymmetric and spherical models (Cappellari 2020; Zhu et al. 2023), indicating the robustness of the methods.

For our joint lensing and dynamical analysis, we used the posteriors of the lens mass model as priors and obtained the posteriors of the mass model by requiring the kinematic data to be well described by the updated model. As the dynamical model requires a full 3D description, we accounted for the projection effect introduced by the intrinsic shape of the deflector galaxies. Using a Singular Isothermal Ellipsoid (SIE) mass model and a Jaffe (Jaffe 1983) stellar light tracer profile, Huang et al. (2025) show that the aperture integrated stellar velocity dispersion within R_eff modeled with spherical JAM σ_sph^ap can be biased at the level of around −4% to +1% (Fig. 8 in Huang et al. 2025), depending on the shape and apparent projected axis ratio of the galaxy. We corrected for this bias statistically by a probabilistic deprojection based on the observedprojected axis ratio and a realistic proposed intrinsic shape distribution. Residual random uncertainties are 1% for typical ellipticity values (Huang et al. 2025). We adopted the “axisymmetric correction factor” σ_axi^ap/σ_sph^ap proposed and computed in a forward-modeling fashion by Huang et al. (2025). There, the authors used an axisymmetric intrinsic shape distribution from Li et al. (2018) and assumed random distribution of the inclination angle for the model predictions of σ_axi^ap within R_eff. By circularizing the projected mass and stellar tracer profile, they calculated the model prediction of σ_sph^ap, and constructed a distribution of σ_axi^ap/σ_sph^ap as a function of the projected axis ratio. In this work, we further accounted for the radial dependence of σ_axi^ap/σ_sph^ap for both aperture kinematics and IFU-type stellar velocity dispersion, which are integrated with shell bins. The final σ_axi^ap/σ_sph^ap distribution we applied to the stellar velocity prediction σ_sph^ap as calculated from Eq. (12) is matched to the aperture size or annulus radii for each lens. Furthermore, we included an empirical Gaussian PSF with FWHM = $0\stackrel{″}{.}8$ in the prediction of σ_axi^ap and σ_sph^ap, which reduces numerical stochasticity in the correction factor and better matches the seeing of the data, but does not rescale σ_axi^ap/σ_sph^ap significantly enough so as to introduce biases.

A further improvement can be obtained by modeling them using the state-of-the-art axisymmetric JAM method directly (e.g., Shajib et al. 2023; Yıldırım et al. 2023; Wang et al. 2025) on the 2D velocity dispersion map, gaining more access to the inclination angle and anisotropy, as well as information about the deflector ellipticity. However, applying a full parameter inference with axisymmetric JAM to the complete lensing+kinematics parameter space is extremely computationally demanding due to the difficulty of simultaneously fitting the lensing and dynamical variables in high-dimensional parameter spaces, and it is beyond the scope of this paper. Methods for accelerating the modeling process have been explored (Gomer et al. 2023; Wang et al. 2025), and those methods will be employed in future analyses by our collaboration.

7.2.3. LoS convergence

The external convergence κ_ext is inferred by using weighted number count statistics to estimate the relative density of a given lens field. We used data from the Millennium Simulation (Springel et al. 2005) to map these statistics to similar LoSs in Millennium, which has the value of κ measured at high angular resolution (see Hilbert et al. 2009). These techniques have been used extensively by the TDCOSMO collaboration and its predecessors. For complete overviews, see Rusu et al. (2017) and Wells et al. (2024).

We used the same set of lenses from the SL2S survey, discussed in Section 4.2.2, to build a hierarchical model of the LoS convergence across the sample. This model is discussed in detail by Wells et al. (2024). In brief, the population-level distribution of κ_ext follows a generalized extreme-value distribution with a mean and scatter in their respective parameters. These parameters are then constrained on the population level, which allows the population-level prior to deviate from random LoSs of the Millennium Simulation. The population-level constraints are then applied to the individual LoSs of the lenses, providing a posterior p(κ_ext) that includes as a prior the population-level inference results. Our model demonstrates with high confidence that the LoSs in our sample are drawn from a population that is overdense compared to the population of all LoSs in the Universe, as expected for massive early-type galaxies (Auger 2008; Treu et al. 2009). This demonstrates that correcting for LoS effects on a population level is a necessary step when doing time-delay cosmography with large samples of lenses. For our hierarchical inference, we used the individual posterior distributions of κ_ext (see e.g., Fig. 6).

8. Cosmological results in flat ΛCDM

In this section, we adopted the flat ΛCDM model as our baseline cosmological model. The cosmological constraints in alternative cosmological models are presented in Section 9. We first considered the cosmological constraints obtained only from the lensing information in Section 8.1 before combining them with other cosmological probes in Section 8.2.

8.1. Lensing-only constraints

For the lensing-only analysis, we took uniform priors on H₀ and Ω_m in the range [0, 150] km s⁻¹ Mpc⁻¹ and [0.05, 0.5], respectively. Other priors on the lens population model parameters are summarized in Table 3.

We considered separately the addition of the SLACS and SL2S datasets of galaxy-galaxy lenses to the sample of eight TDCOSMO-2025 time-delay lenses. We list the combinations of external lens samples and external cosmological probes considered in this work in Table 4.

In our baseline ΛCDM analysis, we find $H_0=73.7_{-4.4}^{+4.7}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$ from only the eight TDCOSMO-2025 lenses, at 6.2% precision. The TDCOSMO-2025, SLACS, and SL2S samples yielded consistent constraints in terms of λ_int and can thus be combined. Including the SLACS lenses to further constrain the mass profile of our deflector population, we obtain $H_0=76.8_{-4.6}^{+3.9}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$, at 5.5% precision. The inclusion of the SL2S lenses yields $H_0=74.2_{-4.5}^{+4.6}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$. Combining all lens samples together (TDCOSMO-2025 + SLACS + SL2S) gives similar constraints as TDCOSMO-2025 + SLACS, $H_0=77.8_{-4.7}^{+3.7}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$, at 5.4% precision. The median, 16th and 84th percentiles of the posterior distributions for the cosmological parameters and deflector population model parameters are reported in Table 6.

8.2. Combination with independent cosmological probes

Time-delay cosmography primarily measures H₀ but offers limited information on Ω_m. In this section, to mitigate the degeneracy between H₀ and Ω_m and enhance the precision of ourmeasurements, we included external probes that provide stronger constraints on Ω_m than time-delay cosmography alone. All resulting constraints are summarized in Table 6.

8.2.1. Pantheon+

We used the Pantheon+ SN dataset (Brout et al. 2022; Scolnic et al. 2022), which includes 1701 light curves of 1550 distinct SNe Ia spanning a redshift range from z = 0.001 to 2.26. We incorporated the data and covariance matrix (including systematic uncertainties) provided by Brout et al. (2022) but excluded the SH0ES absolute SN Ia magnitude calibration. Pantheon+ effectively provided a prior on Ω_m (i.e., Ω_m = 0.334 ± 0.018), while the absolute SNe Ia distances were anchored to our sample of time-delay lenses.

In flat ΛCDM with the Pantheon+ SNe, we obtain $H_0=71.6_{-3.3}^{+3.9}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$ using the TDCOSMO-2025 sample alone. When combined with the SLACS and SL2S lenses, the constraint on H₀ improves to $H_0=74.3_{-3.7}^{+3.1}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$, at 4.6% precision. The posterior distribution for the cosmological parameters and deflector population model parameters are shown in Fig. 8. The only parameter that differs between the three lens samples is the intrinsic scatter on λ_int that is found to be consistent with zero for TDCOSMO-2025 but not for SLACS. As discussed in Section 10.1 and Appendix C, the intrinsic scatter somewhat limits the improvement in precision obtained by adding SLACS to TDCOSMO-2025. However, α_λ, i.e., the parameter describing the slope of λ_int with R_eff/θ_E of the deflector, is consistent across all three lens samples: we find ${\alpha }_{\lambda }=0.{00}_{-0.05}^{+0.05}$ for TDCOSMO-2025 alone, ${\alpha }_{\lambda }=0.{03}_{-0.06}^{+0.05}$ for TDCOSMO-2025 + SLACS, and ${\alpha }_{\lambda }=-0.{01}_{-0.05}^{+0.05}$ for TDCOSMO-2025 + SL2S. These values are consistent with zero, providing no evidence that the inferred λ_int varies as a function of the distance from the center of the lens galaxy.

Fig. 8. Posterior distributions for a subset of the parameters described in Table 3. Blue: constraints from the eight TDCOSMO-2025 lenses, including JWST kinematic data for RX J1131−1231. Orange: constraints from the TDCOSMO-2025 sample + 11 SLACS lenses with KCWI IFU kinematics data. Pink: constraints from TDCOSMO-2025 + SLACS + SL2S sample. All contours include Ωm constraint coming from the Pantheon+ sample (Scolnic et al. 2022); Brout:2022. Priors are described in Table 3. The posteriors distribution of H0 and λint, 0 were blinded during the analysis. Results in combination with the DES-SN5YR likelihood are very similar (see Table 6) and are not shown for conciseness.

8.2.2. DES-SN5YR

Similarly to Section 8.2.1, we included the SNe Ia data from the DES (DES Collaboration 2024), without anchors. We adopted the Year-5 data release (i.e., DES-SN5YR¹¹). In flat ΛCDM, the Ω_m prior provided by this dataset is equivalent to Ω_m = 0.352 ± 0.017. For the TDCOSMO-25 sample, we obtain $H_0=71.2_{-3.6}^{+3.7}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$, at 5.1% precision. In combination with the SLACS and SL2S lenses, we obtain a more precise estimate of $H_0=73.9_{-3.0}^{+3.4}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$, at 4.3 % precision.

8.2.3. DESI BAO DR2

The DESI collaboration recently released their baryon acoustic oscillations (BAO) measurement in DESI-Collaboration (2025). We used the 13 distance measurements, including covariances, provided by DESI¹², in seven redshift bins covering the range z = 0.295–2.33. The inferred distances by DESI are relative to the sound horizon r_d through the ratios D_M(z)/r_d and D_H(z)/r_d, where D_M(z) is the transverse comoving distance and the Hubble distance D_H(z) = c/H(z). To break the degeneracy between r_d and H₀, the DESI collaboration adopted a prior either from Big Bang nucleosynthesis (BBN) or from the Planck CMB data to constrain the sound horizon scale.

In our analysis, we used the BAO uncalibrated distance ruler. We adopted a flat uninformative prior on r_d ∼ 𝒰, (0, 300) Mpc. The TDCOSMO-2025 dataset provided the absolute distance measurement necessary to anchor the BAO distances and break the H₀ – r_d degeneracy.

In flat ΛCDM, our measurement, combining DESI BAO with the TDCOSMO-2025, SLACS, and SL2S lensing datasets, yields $H_0=74.8_{-3.4}^{+3.5}$ km s⁻¹ Mpc⁻¹ at 4.6% precision. We also obtain a measurement of the sound horizon scale, $135.6_{-6.1}^{+6.7}$ Mpc, independent of assumptions about early-Universe physics. This measurement differs from that of Planck by ∼1.8σ and is another aspect of the Hubble tension as discussed in, e.g., Arendse et al. (2020).

8.2.4. Planck

As our previous time-delay cosmography analyses (TDCOSMO-1 and TDCOSMO-4) showed significant tension with Planck’s determination of H₀, we refrained from combining our results with Planck in flat ΛCDM. However, in Section 9, we explored different extensions of the ΛCDM model where this tension is less significant, allowing us to combine the TDCOSMO-2025 constraint with the CMB observations from Planck. In this case, we adopted the plikHM_TTTEEE_lowl_lowE chains from Planck Collaboration VI (2020) in various extensions of the ΛCDM model¹³.

9. Alternative cosmological models

Time-delay cosmography is primarily sensitive to H₀, and only weakly dependent, given the current sample size, on parameters used in extensions of flat ΛCDM, such as Ω_k, w₀, and w_a. The sensitivity to H₀ and the different direction of residual parameter degeneracies make it highly complementary to other cosmological probes (Linder 2011). In this section, we report the cosmological constraints provided by TDCOSMO-2025, in combination with the CMB data from Planck, BAO measurements from DESI, and SN Ia data from Pantheon+, in open ΛCDM, flat wCDM, flat w₀w_aCDM, and flat w_ϕCDM cosmologies. All priors used in these alternative cosmological models are listed in Table 4, and the parameter constraints are given in Table 7.

9.1. Open ΛCDM

A simple modification of the flat ΛCDM cosmology is the open ΛCDM cosmology, which allows for spatial curvature, Ω_k ≠0. We adopted a flat prior on Ω_k ∼𝒰 (−0.5, 0.5) and imposed the additional condition Ω_Λ = 1 − Ω_m − Ω_k > 0. The energy density in photons and neutrinos can be neglected with the current level of precision.

Using only the TDCOSMO-2025, SLACS, and SL2S lensing datasets, we obtain $H_0=77.1_{-4.1}^{+4.9}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$. Ω_k is poorly constrained by time-delays alone, with ${\mathrm{\Omega }}_k=-0.{117}_{-0.245}^{+0.305}$. However, H₀ remains well constrained, which provides an indirect constraint on the curvature when combining time-delay cosmography with the CMB data. However, our measurement is in > 3σ tension with that of Planck, which might lead to underestimated uncertainties when combining these two datasets. Still, this combination yields $H_0=65.2_{-2.3}^{+1.8}$ km s⁻¹ Mpc⁻¹ and ${\mathrm{\Omega }}_k=-0.{006}_{-0.006}^{+0.005}$ in open ΛCDM cosmology. We show the posterior distribution of H₀ and Ω_k in Fig. 9 for TDCOSMO, Planck, and their combination.

Fig. 9. Posterior distribution of H0 and Ωk in open ΛCDM cosmology. The black contours show the constraints from TDCOSMO alone, including the TDCOSMO-2025, SLACS, and SL2S lensing datasets. The gray contours show the constraints from Planck. The orange contours correspond to the combination of Planck and TDCOSMO. The two datasets exhibit a tension exceeding 3σ, indicating that this joint constraint should be interpreted with caution. The contour levels represent the 1σ and 2σ constraints.

9.2. Flat wCDM

We considered the flat wCDM cosmological model in which dark energy differs from a cosmological constant and is described by an equation-of-state parameter w. When w equals −1, this model reduces to the flat ΛCDM case. For the parameter w, we imposed a uniform prior over the interval [ − 1.5, 0.5], while maintaining the same uniform priors on H₀ and Ω_m as used in the flat ΛCDM model.

We show the parameters constraints in this cosmology in Fig. 10. The marginalized H₀ posterior for the TDCOSMO-2025+SLACS+SL2S lensing dataset is $H_0=78.0_{-7.3}^{+6.6}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$. This is higher and more uncertain than in the flat ΛCDM cosmology due to the degeneracy between H₀ and w.

Fig. 10. Posterior distribution of H0 and w in flat wCDM cosmology. The black contours show the constraints from TDCOSMO alone, including the TDCOSMO-2025, SLACS, and SL2S lensing datasets. The gray contours show the constraints from Planck. The combination of TDCOSMO with the Planck, DESI BAO, and Pantheon+ is shown with colored contours, as described in the figure legend. The contour levels represent the 1σ and 2σ constraints.

In combination with Planck, our constraints on the dark energy equation-of-state is $\mathit{w}=-1.{38}_{-0.08}^{+0.13}$ and $H_0=79.7_{-4.3}^{+3.2}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$. In this combination, the dark energy equation-of-state deviates from a cosmological constant at ∼4σ. However, this constraint is mainly driven by the tension regarding H₀ between the two datasets, which is already observed in the ΛCDM cosmology.

When combining with the DESI BAO data or Pantheon+, our measurement is in agreement with the ΛCDM model, finding no deviation from a cosmological constant. The marginalized posterior distribution on w is $\mathit{w}=-0.{93}_{-0.09}^{+0.08}$ for TDCOSMO+DESI BAO and $\mathit{w}=-0.{86}_{-0.18}^{+0.14}$ for TDCOSMO + Pantheon+.

9.3. Flat w₀w_aCDM

The flat w₀w_aCDM cosmology includes a dark energy component characterized by an equation-of-state parameter w(z) that can evolve with time. In this model, the dark energy equation-of-state is parameterized as

$$\begin{array}{c}\hfill w(z)=w_0+w_a\frac{z}{1+z},\end{array}$$(24)

following Chevallier & Polarski (2001) and Linder (2003). In this cosmology we adopt flat priors on w₀ ∼𝒰 (−1.5, 0.5) and w_a ∼𝒰 (−10, 10).

Our lensing-only dataset weakly constrains these two parameters, whose posterior distributions are shown in Fig. 11. However, H₀ is constrained to be $H_0=80.9_{-8.0}^{+11.2}$ km s⁻¹ Mpc⁻¹.

Fig. 11. Posterior distribution of H0, w0 and wa in flat w0waCDM cosmology. The combination of TDCOSMO (including the TDCOSMO-2025, SLACS, and SL2S lensing datasets) with the DESI-DR2 BAO and Pantheon+, and the combination of both external datasets, is shown with colored contours described in the figure legend. The contour levels represent the 1σ and 2σ constraints.

The combinations of TDCOSMO individually with Planck, DESI-BAO, or Pantheon+ are consistent with a cosmological constant (i.e., w₀ = −1 and w_a = 0) within 2σ. This is expected given the weak constraints on the dark energy parameters derived only from the TDCOSMO dataset. The combination of TDCOSMO with both Pantheon+ and DESI-BAO does not improve the precision on the dark energy parameters compared to the results presented in DESI-Collaboration (2025), but H₀ is well constrained to be $H_0=73.2_{-3.0}^{+4.0}$ km s⁻¹ Mpc⁻¹ in flat w₀w_aCDM cosmology. The strong preference for dynamical dark energy becomes evident only when multiple other probes are combined, such as BAO+CMB or BAO+SN+CMB (DESI-Collaboration 2025). Larger time-delay lens samples in the future will be able to provide precision on the dark energy parameters comparable to other leading probes (Shajib et al. 2025d).

9.4. Flat w_ϕCDM

The w_ϕCDM model describes the dark energy as originating from a scalar field, such as the pseudo-Nambu-Goldstone bosons (Frieman et al. 1995; Shajib & Frieman 2025). This field is allowed to evolve over time, thus giving rise to a dynamic nature in the dark energy. As a result, this is a physics-inspired model of dynamical dark energy, unlike the w₀w_aCDM model, which is phenomenological and physics-agnostic, thereby allowing for violation of standard physics, such as the null energy condition, within its parameter space. In the w_ϕCDM model, the dynamical behavior of w(z) is well-approximated by the formula (Shajib & Frieman 2025)

$$\begin{array}{c}\hfill w(z)=-1+(1+w_0)exp(-\alpha z),\end{array}$$(25)

where w₀ is the value of w at the current epoch (i.e., z = 0). The α parameter depends on the form of the potential of the scalar field, and it falls within the narrow range of 1.35–1.55, which we adopted as a uniform prior range. However, all current and near-future datasets would not be sufficiently sensitive to constrain α beyond this prior range (Shajib & Frieman 2025). As a result, this w_ϕCDM model is a pseudo-one-parameter extension to the ΛCDM.

In this cosmological model, we obtain $H_0=75.9_{-7.8}^{+6.8}$ km s⁻¹ Mpc⁻¹ and ${\mathit{w}}_0=-0.{95}_{-0.39}^{+0.59}$, using only the TDCOSMO-2025+SLACS+SL2S lensing dataset (Fig. 12). In combination with DESI BAO or Pantheon+, we obtain ${\mathit{w}}_0=-0.{86}_{-0.13}^{+0.11}$ and ${\mathit{w}}_0=-0.{92}_{-0.13}^{+0.11}$, respectively. The combination of TDCOSMO-2025 with DESI BAO or Pantheon+ alone does not exclude the cosmological constant (i.e., the case with w₀ = −1), similar to what is found also for the w₀w_aCDM model above.

Fig. 12. Posterior distribution of H0, Ωm and w0 in flat wϕCDM cosmology. The combinations of TDCOSMO (including the TDCOSMO-2025, SLACS, and SL2S lensing datasets) with the DESI-DR2 BAO and Pantheon+ datasets are shown with colored contours, as described in the figure legend. The contour levels represent the 1σ and 2σ constraints.

10. Comparison with the literature

We now briefly discuss our measurement in comparison with previous work by our collaboration (Section 10.1), by other teams using time-delay cosmography(Section 10.2), and within the broader context of the Hubble tension (Section 10.3).

10.1. Comparison with our previous work

In this section, we compare our results with those from our previous work, as summarized by TDCOSMO-1 and TDCOSMO-4 (see Fig. 13). In addition to performing the comparison, we aim to provide some intuition for the factors that contributed to the change in the overall precision and those that changed the median value of the TDCOSMO+SLACS analysis. The TDCOSMO-1 and TDCOSMO-4 papers were based on seven out of the eight time-delay lenses studied in this work, with the same identical time delays and lens models. The new additional lens (WGD 2038−4008) has considerably larger uncertainties on the time delays (Wong et al. 2024), so we expected it not to have a substantial effect on precision and point estimate. Its central estimate of H₀ was lower than the average for the other seven. Therefore, we expected it to lower H₀ by a modest amount.

Fig. 13. Comparison of the H0 measurements by the TDCOSMO collaboration (and its predecessor H0LiCOW) in chronological order. The measurements at the top by Wong et al. (2020) and Millon et al. (2020b) are “assertive” in the definition of Treu et al. (2022). They are based on power-law and composite models to describe the mass density profile of the deflector galaxies, thus implicitly breaking the MSD and obtaining ∼2% precision. The TDCOSMO-4 results from Birrer et al. (2020) are “conservative” (Treu et al. 2022). They are based on the same data but obtain larger uncertainties by introducing the internal MSD parameter λint, and constraining it with unresolved stellar kinematics. The Birrer et al. (2020) measurement in combination with SLACS is shown as a dashed line for historical purposes but should not be used anymore, as it was later discovered that the stellar velocity dispersions based on low SNR (∼9 Å−1) spectra from SDSS suffer from systematic errors and covariance (Knabel et al. 2025a). The new measurements presented in this work are shown in the bottom panel. They are “conservative” in terms of uncertainties and constrain the MSD using new stellar kinematics based on high SNR JWST-NIRSpec, VLT-MUSE, and Keck-KCWI spectra, as well as an improved methodology (Knabel et al. 2025b).

10.1.1. TDCOSMO-only

Our inferred Hubble constant ($H_0=71.6_{-3.3}^{+3.9}$  km s⁻¹ Mpc⁻¹) is in excellent agreement with the results of both TDCOSMO-1 (H₀ = 74.2 ± 1.6 km s⁻¹ Mpc⁻¹) and TDCOSMO-4 ($H_0=74.5_{-6.1}^{+5.6}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$) in the flat ΛCDM cosmology. The exact comparison depends somewhat on the priors on Ω_m, owing to the mild degeneracy with H₀. TDCOSMO-4 imposed Ω_m = 0.298 ± 0.022 based on the Pantheon SNe, while TDCOSMO-1 used the uniform prior Ω_m ∈ [0.05, 0.5] that is the same one used by Wong et al. (2020). Using the same prior as in TDCOSMO-1 instead of Pantheon+, we obtain $H_0=73.7_{-4.4}^{+4.7}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$.

The TDCOSMO-1 analysis was based on stricter assumptions on the mass profile than those assumed here and by TDCOSMO-4. Effectively, the TDCOSMO-1 power-law result is equivalent to imposing λ_int = 1 in the analysis of TDCOSMO-4, since the data were almost exactly the same. It is thus not surprising that the precision of TDCOSMO-1 was better than that of TDCOSMO-4, and that in the present analysis, which is based on much improved data. TDCOSMO-4 found ${\lambda }_{\mathrm{int}}=1.{02}_{-0.09}^{+0.08}$, resulting in the effectively the same central value of the TDCOSMO-only estimate, but with larger uncertainties (H₀ ∝ λ_int, all other things being equal).

In this work, we made several changes that affect the parameter λ_int, the anisotropy parameter, the power-law slope γ_pl, and ultimately H₀. From a data standpoint, we replaced all stellar velocity dispersions with new measurements based on significantly improved data and methods. Qualitatively, we expected this would change the inference on H₀ in the following manner.

First, the new stellar velocity dispersions for the time-delay lenses have an average total uncertainty of 4.0%, compared with 7.7% for the values used by TDCOSMO-4. The average total uncertainty is defined as < δσ_i/σ_i> _i, averaged over the seven lenses in common with TDCOSMO-4, with δσ_i being the sum in quadrature of average statistical and systematic errors. All other things being equal, this should have considerably improved the precision on H₀, provided that no intrinsic scatter is detected. TDCOSMO-4 did not find any intrinsic scatter in the time-delay lenses. Remarkably, the intrinsic scatter remains consistent with zero even with our improved kinematic measurements. We note that the stellar velocity dispersion is not the only source of uncertainty; the time delays, lens model, and LoS convergence also contribute to the final uncertainty on H₀ (Wong et al. 2020).

Additional information was obtained from the spatially resolved kinematics of RX J1131−1231, which constrain the mass density slope and the internal MSD. As shown by Shajib et al. (2023), the use of the KCWI data alone makes this one system almost as precise ($H_0=77.1_{-7.1}^{+7.3}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$) as the seven lenses combined in TDCOSMO-4 ($H_0=74.5_{-6.1}^{+5.6}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$). The addition of the JWST data was expected to further improve the precision for this system (further gains can be achieved using the full 2D kinematics; Yıldırım et al. 2023; Wang et al. 2025). Therefore, based on these factors alone, we were expecting an improved precision for the TDCOSMO-only measurement, again in the absence of internal scatter on λ_int, depending exactly on the overlap between the individual datasets. As shown in Appendix B, the KCWI and JWST data are in mutual statistical agreement, improving the overall precision stemming from RX J1131−1231. However, for this one system, the uncertainty on κ_ext provides an effective noise floor.

In addition, the average ratio of the newly measured stellar velocity dispersions with respect to previously measured values provides us with a posteriori insight into the agreement between the measurement presented here and that presented in TDCOSMO-4. While we did not blind the estimates of the velocity dispersions, in order to minimize confirmation bias, three different team members¹⁴ measured the velocity dispersions independently and without comparing the results. In some cases, the new results differ significantly from the old ones; four are within 1σ and three are between 1–2σ of the old uncertainties (see Fig. A.1). Moreover, several other factors besides the velocity dispersions contribute to determining H₀, such as spectroscopic apertures and seeing (see the example in Fig. 4). Furthermore, H₀ depends not only on the measured stellar velocity dispersion, but also on the inferred slope, anisotropy parameter, and external convergence (see example in Appendix B) in non-trivial ways. The stellar velocity dispersion measurements were frozen before unblinding and before performing the aggregate comparison described in this section. In a few cases, TDCOSMO-4 combined multiple velocity dispersion measurements for the same system. For the simple comparison we make here, we took an inverse variance weighted average of the four dispersion measurements for DES J0408−5354 used by TDCOSMO-4, adding the covariance in quadrature. For PG 1115+080, we took the central aperture, which is the closest to the JWST one adopted here. For RX J1131−1231, we adopted the stellar velocity dispersion within half of the half-light radius as our point of comparison. In sum, for the seven lenses in common with TDCOSMO-4, the average ratio is σ_new/σ_old = 0.986 ± 0.049. In reality, of course, the weight of each lens is different, and those with the best measured time delays and lens models will carry more weight. This simple check provides some insight as to why the overall mean of the λ_int posterior remained the same within the errors, although only the complete analysis can account for all the changes. In fact, H₀ did not scale exactly as λ_int between TDCOSMO-4 and this work. This is due to two effects. First, the addition of WGD 2038−4008 that reduces H₀ overall, although this is a small contribution as described above, owing to its large uncertainty. Second, and most importantly, the spatially resolved kinematics of RX J1131−1231 prefers a higher slope γ_pl than the lens model, and prefers lower λ_int than the TDCOSMO-2025 sample average.

We have also improved the anisotropy model, now accounting for projection effects that affect the model-predicted stellar velocity dispersion when compared with the observed one. This correction shifts the inferred H₀ by at most 1-2% down (Huang et al. 2025).

In conclusion, the improvement in precision with respect to TDCOSMO-4 can be attributed to new measurements and an enhanced methodology. The good agreement between this blind analysis and TDCOSMO-4 and TDCOSMO-1 provides a powerful validation of the method. Further improvements are possible with improved analysis of 2D kinematics, time delays, LoS characterization, new lens models, and expanded datasets.

10.1.2. TDCOSMO + external lens samples

Importantly, H₀ from the TDCOSMO-2025+SLACS is in excellent agreement with that from the TDCOSMO-2025 alone, given the uncertainties. Compared to TDCOSMO-4, the TDCOSMO+SLACS inference of H₀ moved from $H_0=67.4_{-3.2}^{+4.1}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$ to $76.8_{-4.6}^{+3.9}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$ in our UΛCDM cosmology (H₀=$74.8_{-3.4}^{+3.5}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$ in combination with DESI-DR2 that yields the more similar Ω_m to the one adopted in TDCOSMO-4), to be compared with the TDCOSMO-2025-only result $H_0=73.7_{-4.4}^{+4.7}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$ in the same cosmology.

The improved precision and change in the median value can be understood as follows. We went from using 33 SLACS lenses with SDSS spectra in TDCOSMO-4 to 11 lenses with Keck-KCWI spectra here. We now know that the SDSS-SLACS velocity dispersion measurements is affected by systematic uncertainty at the level of 3.3% and covariance at the level of 2.3% (Knabel et al. 2025a), making them not suited for cosmography and therefore the old TDCOSMO+SLACS result from TDCOSMO-4 suffered from significant underestimation of uncertainty (random and systematic) and should not be considered as valid anymore. However – setting aside systematic errors in the SDSS stellar velocity dispersion – the SNR of the Keck-KCWI spectra is ∼18 times higher per lens than that of SDSS. The covariance between the SLACS lenses’ stellar velocity dispersion profiles, due to stellar template choices, sets a floor of about 1.5-2% to the precision on H₀ attainable with this sample. Overall, we expected the new KCWI-SLACS dataset to reduce the errors on H₀ significantly more than in TDCOSMO-4, even with just 11 lenses, if no intrinsic scatter in λ_int was detected. However, our analysis showed that the SLACS lenses require significant scatter in λ_int, as shown in Appendix C. The improvement in precision provided by the SLACS sample with respect to TDCOSMO-4 is thus limited not by the quality of the data, but by the intrinsic scatter of the sample. We note that the intrinsic scatter for the TDCOSMO-2025 sample is consistent with zero. However, the posterior distribution of σ(λ_int) exhibits a long tail, making it compatible at 95% confidence level, with the intrinsic scatter detected in the SLACS-KCWI sample. It remains to be investigated whether the intrinsic scatter of the SLACS-KCWI can be reduced with more sophisticated models (e.g., 2D dynamical models and more advanced lens models) or whether it is an irreducible property of the sample.

The SL2S dataset is new and consists of only four lenses with integrated stellar velocity dispersion. Thus, we did not expect it to significantly improve the precision of the measurement. Nevertheless, we found it to be in good agreement with the other two samples, which provides an important check of systematics. Further follow-up of the SL2S sample would increase the sample size and potentially contribute more to improving the precision in combination with TDCOSMO, providing a more stringent validation test.

10.2. Comparison with other work

We provide a brief comparison here with the broader literature on time-delay cosmography, referring the reader to recent review articles for more details (e.g., Treu et al. 2022, 2024; Birrer et al. 2024; Suyu et al. 2024).

Measurements based on lensed quasars have been around for decades and have improved considerably since the early days, owing to improvements in data quality and modeling methods. The work presented here is the culmination of years of efforts to perfect the so-called “simply-parametrized” approach, in which the mass distribution or gravitational potential of the main deflector is described by a physically motivated functional form depending on a relatively small number of free parameters (∼10–20 typically), constrained by high resolution images and stellar kinematics (see, e.g., Oguri 2007, for an early result).

The alternative approach is the so-called “free-form” method, in which the mass distribution or gravitational potential of the main deflector is described by pixels or a basis set (Coles 2008; Paraficz & Hjorth 2010). In this approach, the number of free parameters vastly exceeds the number of constraints, requiring regularization to avoid overfitting. These methods have also improved considerably since the early days, with more realistic priors and regularization, and have benefited from better data. With reasonable priors and regularization, Denzel et al. (2021) finds $H_0=71.8_{-3.3}^{+3.9}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$, although the result and precision are heavily dependent on those choices.

A recent and exciting development has been the discovery of multiply imaged SNe (Kelly et al. 2015), the astronomical source originally proposed by Refsdal (1964) for time-delay cosmography. The first measurements of H₀ based on lensed SNe have been published. Kelly et al. (2023) find $H_0=66.6_{-3.3}^{+4.1}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$ based on SN “Refsdal” (see also Grillo et al. 2024) and Pascale et al. (2025) find $H_0=71.8_{-8.1}^{+9.2}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$ based on SN “H0pe”. Although the uncertainties based on just these two systems are too large to settle the Hubble tension, large numbers of lensed SNe are expected to be discovered in the near future (e.g., Arendse et al. 2024) and will further boost the precision of time-delay cosmography (Suyu et al. 2020; Birrer et al. 2022; Shajib et al. 2025d).

The results presented here are overall in excellent agreement with recent work on time-delay cosmography. SN “Refsdal” is the most apparently discrepant result, and yet it is consistent with our measurement within 2σ. Considering the widely different sources, deflectors, and methods involved in these independent studies, the overall agreement can be viewed as a validation of time-delay cosmography.

10.3. Comparison with non-lensing measurements

In Fig. 14, we show a selection of recent measurements of the Hubble constant in flat ΛCDM (for a complete list, see the recent review by Di Valentino et al. 2025). While the selection is by no means unique, we aimed to represent a range of independent methods. As discussed in the Introduction, the early-Universe probes clearly prefer H₀ around 66–68  km s⁻¹ Mpc⁻¹, while late-Universe probes prefer values above 70 km s⁻¹ Mpc⁻¹. Our measurements are in excellent agreement with other late-Universe probes. Importantly, this is true also for our measurement $H_0=77.8_{-4.7}^{+3.7}$  km s⁻¹ Mpc⁻¹ from TDCOSMO-2025+SLACS+SL2S without utilizing any SN information, which is thus completely independent of all other methods. However, given the uncertainties, the results from the early-Universe probes are also consistent within 3σ (see discussion in Section 9 for a more general comparison with the CMB results). Further improvements in the precision of the TDCOSMO measurements are planned and will sharpen the contribution of time-delay cosmography to the “Hubble tension” debate.

Fig. 14. Illustration of the current state of the “Hubble tension” in flat ΛCDM. Only a selection of recent measurements is shown for the early- and late-Universe probes, in order to avoid overcrowding the plot (for an exhaustive list, see Di Valentino et al. 2025). Among the early-Universe probes, we show measurements from the CMB (Planck and ACT-DR6; Planck Collaboration VI 2020; Louis et al. 2025), and baryonic acoustic oscillations (BAO) plus Big Bang nucleosynthesis (BBN) (DESI-Collaboration 2025). Among the late-Universe probes, we show the Cepheid+SN measurement from the SH0ES team (Breuval et al. 2024), the tip of the red giant branch (TRGB) plus SN measurement from the CCHP team (Freedman et al. 2025), the TRGB plus surface brightness fluctuation (SBF) measurement by Jensen et al. (2025), along with the TDCOSMO-2025 and TDCOSMO-2025+SLACS+SL2S time-delay cosmography measurements in combination with the Pantheon+ likelihood (Scolnic et al. 2022); Brout:2022.

11. Systematics

In this section, we discuss residual systematic errors and our efforts over the past few years to quantify and mitigate them (Section 11.1). We also discuss and quantify potential selection effects (Section 11.2).

11.1. Model assumptions

In addition to those discussed in this paper, several other sources of systematic uncertainties should be considered when interpreting our results. A brief summary is provided here, referring the reader to the original references for more details. The important point is that all of these uncertainties are much smaller than those associated with our measurement of H₀, and are random and uncorrelated. Therefore, they will not change our result by more than a fraction of our error bar. Below, we list the sources of uncertainty that are not yet included in our baseline analysis, unless explicitly stated otherwise:

• Parametrization of the mass model: TDCOSMO-1 showed that using a so-called “composite” model, i.e., a constant stellar mass-to-light ratio plus an NFW dark matter halo, changes H₀ by 0.2% with respect to adopting a power-law model. This is important because the composite model has the flexibility to represent mass density profiles that differ from a power law. In practice, however, the data reveal that the two components combine to form an almost perfect power-law profile in the radial range around the Einstein radius, a phenomenon known as the “bulge-halo conspiracy” (Treu & Koopmans 2004; Dutton & Treu 2014).

• Uncertainties on the lens models: Time-delay cosmography requires accurate mass models to infer the difference in the Fermat potential, under the assumption of an elliptical power-law model. Random errors are at the 1–5% level (Wong et al. 2020) for data of our quality. We performed two tests of potential systematics. Shajib et al. (2022) presented a blind analysis of WGD 2038−4008 based on the two independent softwares used by our collaboration (GLEE and LENSTRONOMY, respectively; Suyu & Halkola 2010; Suyu 2012; Birrer & Amara 2018), independent modeling choices, and PSF reconstruction. They found the results to agree within the errors. Williams et al. (2025) modeled new JWST-NIRcam images of WFI 2033−4723 and found the results to be consistent within the error to those obtained by Rusu et al. (2020) using HST images.

• Dark matter substructure: In our mass models, we assumed the mass distribution of the deflector could be described by a smooth function. In reality, cold dark matter is clumpy on subgalactic scales. However, this approximation introduces only ∼1% random scatter in the observed time-delay distance (Gilman et al. 2020), which is negligible given the other sources of uncertainty.

• Peculiar velocities with respect to the CMB reference frame affect the time delay distances by a small amount. Accounting for both the motion of the Milky Way with respect to the CMB and peculiar velocities, Dalang et al. (2023) estimated the effect to be ∼0.25 % for the TDCOSMO-1 sample.

• A central supermassive black hole could affect the interpretation of the stellar velocity dispersion at the JWST-NIRSpec angular resolution (Wang et al. 2025). This is not accounted for in this work and will be included in future works, when we will analyze 2D JWST kinematic maps in detail.

• Departures for spherical symmetry in dynamical models: As discussed by Huang et al. (2025), after we apply our correction for projection effects, residual random errors are at the level of 1–2% in λ_int, which we account for in this work. Impacts on the deprojection from a potential triaxiality are below the 1% level.

• Our time-delay measurement method has been tested in a blind data challenge (Liao et al. 2015). The COSMOGRAIL group’s technique, namely PYCS3 (Tewes et al. 2013; Millon et al. 2020c), has consistently produced accurate and precise time delays, with residual systematics and biases of less than 1% (Bonvin et al. 2016). Systematic errors due to the microlensing time-delay effect (Tie & Kochanek 2018) are on the order of a few hours and have only a small impact on our measurements, given that the typical time delays in our sample of lensed quasars span weeks to months (Chen et al. 2018).

• Testing hierarchical inference:Gomer et al. (2022) tested the ability of the hierarchical framework to infer an unbiased H₀ by using mock lenses with a composite (dark matter + baryons) mass distribution that they modeled with a power-law mass profile. The framework successfully recovered H₀ within 1.5σ of the fiducial value when using time-delay lenses alone. The precision was improved to ∼2% and accuracy to 0.7% median offset when combining time-delay lenses with an external sample of non-time-delay systems at lower z, similar to the SLACS sample. This confirmed that non-time-delay lenses enhance constraints on H₀ when assuming a shared population, despite systematic differences in radii probed by the lensed images (i.e. difference in R_eff/θ_E).

• Departure from elliptical symmetry in the lens models:Van de Vyvere et al. (2022) showed that for realistic values of azimuthal structure higher than order 2 (boxyness and diskyness), the effect on H₀ is less than 1%.

• Nearby perturbers in the SLACS and SL2S lens models: Whereas for the TDCOSMO-2025 lenses we modeled explicitly all the nearby perturbers that contribute more than 5 × 10⁻⁴ in the flexion shift (McCully et al. 2017), this was not done for the SLACS and SL2S lenses. Visual inspection of the selected SLACS and SL2S indicates no major perturber, indicating that the effect is likely to be small. However, this step of the analysis will likely need to be implemented in the external lens datasets as precision improves, to further homogenize the treatment of the two samples.

11.2. Selection effects

Selection effects could potentially bias our results if they are large enough and not properly accounted for (see, e.g. Sonnenfeld et al. 2023, for a comprehensive discussion of strong lensing selection functions and selection effects). There are two kinds of selection effects that could play a role: those affecting the time-delay lenses in an absolute sense and differential ones affecting the combination of time-delay and non-time-delay lenses.

11.2.1. Selection effects in time-delay lenses

The TDCOSMO-2025 sample is composed almost exclusively of quadruply imaged quasars (seven out of eight). This is for good reasons: quads have significantly more information content per system than doubles. Therefore, they have been prioritized for follow-up and analysis. However, quads could, in principle, be a biased subset of galaxies; for example, they preferentially select systems with more ellipticity or external shear that subtend a larger area within the four-image caustic. This would have implications only insofar as priors or other information from non-quads is used. The bias, excluding a caustic area prior in the individual lens models, is estimated to be sub-percent, and even lower for well-constrained lens models based on our high-resolution imaging data (Baldwin & Schechter 2024). In our analysis, this is a small effect, and it is mostly present in the correction for projection along the LoS in the modeling of stellar velocity dispersion (Huang et al. 2025). Since this is a percent-level correction, we expect the effect to be sub-percent on H₀. To quantify the effect more precisely, one needs to know the parent population and the selection function, both of which are difficult to determine. Based on a simple model for the deflector and source population, Collett & Cunnington (2016) estimate the potential bias to be at the sub-percent level. Sonnenfeld et al. (2023) reaches a similar conclusion that the residual bias is probably sub-percent for current samples. Li et al. (2025) claim a larger effect, but their result is not relevant because their model fails to reproduce the actual observations. Their model predicts ${\lambda }_{\mathrm{int}}=0.{87}_{-0.06}^{+0.07}$ for the TDCOSMO-4 quads, in stark contrast to the measured value there, and that obtained in this paper. Their inference of a lower H₀ than that presented in TDCOSMO-4 is almost entirely due to the underestimated λ_int, and insufficiently describes the selection effects that may affect real lenses.

In conclusion, based on the work done so far, it seems that any residual bias due to unmodeled selection effects is smaller than our uncertainties. However, it is clear that more work remains to be done on this front, as the precision improves. On the one hand, even more realistic descriptions of the parent population and selection function can be implemented in theoretical models. On the other hand, samples of doubles with similar precision to those of quads would provide an empirical test of residual biases.

11.2.2. Differential selection effects between time-delay and non-time-delay lenses

A potential source of bias when combining with external lens samples is due to any differences in the selection function of the time-delay and non-time-delay lenses. In this work, we have assumed that the TDCOSMO-2025, SL2S, and SLACS samples have been drawn from the same population, after applying cuts to match observable properties as closely as possible, such as projected ellipticity and lensing velocity dispersion. Potential differences between the populations are described by allowing λ_int to depend on the ratio of effective and Einstein radius (Eq. 23), which is different for the three samples.

However, the selection functions are inherently different for the three samples. For example, Sonnenfeld (2025) shows that the SLACS follow-up procedure – prioritizing systems with higher stellar velocity dispersion from SDSS – introduces a bias in some of the properties of the SLACS lenses compared to the parent population. The most important for our purposes is that the SLACS lenses’ stellar velocity dispersion could be biased high by up to 5% when matching the mass profile. Part of that bias stems from noisy SDSS spectra used for the pre-selection, in the sense that velocity dispersions that were affected by upward noise fluctuations were preferentially selected. That component is eliminated by our analysis as we replaced in this work the SDSS spectra with much higher SNR KCWI ones. A residual bias, due to the intrinsic scatter in the velocity dispersion of elliptical galaxies, could still be present, though a no-bias scenario is also consistent with the data (ee Fig. 4 of Sonnenfeld 2025). It is for this reason that we present separately the results based on the TDCOSMO-2025 lenses only and those in combination with external lens samples.

Our mutually blind analysis shows excellent agreement between the TDCOSMO-2025, SLACS, and SL2S inference of λ_int within the uncertainties. However, as the precision of the measurement improves, more work will be needed to quantify, mitigate, and correct absolute and differential biases due to the lenses’ selection function in order to ensure that they remain below the statistical errors.

12. Summary and conclusions

We presented cosmological constraints based on the analysis of eight multiply imaged quasars with measured time delays and two external samples of non-time delay lenses. This is the culmination of 5 years of work by the TDCOSMO collaboration, including new data, new methods, and many significant improvements with respect to our earlier works based on 7 time delay lenses (Millon et al. 2020b; Birrer et al. 2020). Following Birrer et al. (2020), we parametrized the mass density profile of the main deflector in terms of the mass-sheet-degeneracy parameter λ_int. Using lensing-only observables, this parametrization is completely degenerate with the Hubble constant and thus maximizes the uncertainty, providing a conservative estimate of the error budget.

Within this framework, the degeneracy is broken by stellar velocity dispersion. Therefore, a major effort went into obtaining higher-quality data and improving the methodology to measure stellar velocity dispersion with respect to our previous work Knabel et al. (2025b). We obtained new JWST-NIRspec spectra for six out of eight time-delay lenses and VLT-MUSE spectra for the other two. For the time delay lens RX J1131−1231, spatially resolved velocity dispersion measurements from both JWST-NIRSpec (Shajib et al. 2025a) and Keck-KCWI spectroscopy (Shajib et al. 2023) were used. New spectra and spectral measurements were obtained for the external lens samples. For the SLACS lenses, SDSS spectroscopy was replaced by much higher SNR (∼160 Å⁻¹) spatially resolved kinematics obtained with KCWI at Keck (Knabel et al. 2025a). For the SL2S lenses, a re-analysis of Keck and VLT data was performed (Mozumdar et al. 2025). New lens models were obtained for both SLACS and SL2S datasets, with improved methods and data (Tan et al. 2024; Sheu et al. 2025). Stringent quality cuts were imposed to guarantee the high quality of all measurements, and to match the observable properties of the TDCOSMO-2025 SLACS and SL2S samples. We also accounted for projection effects in the interpretation of the stellar velocity dispersions (Huang et al. 2025) and improved our description of the velocity dispersion anisotropy and surface brightness profiles.

The main results of this study can be summarized as follows:

• In a flat ΛCDM model, using the TDCOSMO-2025 sample in combination with constraints on Ω_m from the Pantheon+ likelihood, we find $H_0=71.6_{-3.3}^{+3.9}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$. A consistent result $H_0=71.2_{-3.6}^{+3.7}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$ is found by combining with the DES-SN5YR sample instead of Pantheon+.

• The TDCOSMO-2025, SLACS, and SL2S samples yield consistent constraints and can thus be combined to improve the precision. In a flat ΛCDM model, using the TDCOSMO-2025+SLACS+SL2S sample in combination with constraints on Ω_m from the Pantheon+ sample, we find $H_0=74.3_{-3.7}^{+3.1}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$.

• The constraints on H₀ from time-delay cosmography alone are robust to the choice of cosmological model. For the TDCOSMO-2025+SLACS+SL2S dataset we find H₀=$77.8_{-4.7}^{+3.7}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$ in flat ΛCDM, $H_0=77.1_{-4.1}^{+4.9}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$ in open ΛCDM, $78.0_{-7.3}^{+6.6}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$ in wCDM, $80.9_{-8.0}^{+11.2}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$ in w₀w_aCDM, and $75.9_{-7.8}^{+6.8}$ in w_ϕCDM.

• The TDCOSMO data are consistent with the DESI-DR2 BAO data, the Pantheon+ SN, and DES-SN5YR data in all the cosmologies considered here.

• In combination with the Planck likelihood, the TDCOSMO data push the CMB posterior towards flatness in open ΛCDM ${\mathrm{\Omega }}_k=-0.{006}_{-0.006}^{+0.005}$, and towards w = −1 in wCDM $\mathit{w}=-1.{38}_{-0.08}^{+0.13}$.

These results illustrate the power of time-delay lenses to constrain H₀ independently of the local distance ladder, and of other cosmological parameters, thereby contributing to the current scientific debate on the “Hubble tension”. The overall precision has almost doubled with respect to TDCOSMO-4, owing mainly to the use of better spectroscopic data and methods for stellar kinematics. Importantly, replacing the SDSS-based spectra with much higher SNR spectra from the Keck-KCWI brought the SLACS sample into excellent agreement with the TDCOSMO-2025 dataset.

Looking ahead, the TDCOSMO collaboration is working on multiple fronts to further improve our precision accuracy towards our ultimate goal of 1%. First, we are expanding the number of lensed quasars with measured time delays (Dux et al. 2025), cosmography grade models, and ancillary data. Second, using recent improvements in modeling speed (Wang et al. 2025), we plan to take advantage of spatially resolved kinematics from space and from the ground to implement axisymmetric Jeans modeling (see Shajib et al. 2023, for a first illustration). Third, we are expanding the sample of non-time delay lenses with data of sufficient quality to be used for cosmography. These improvements will be the subject of future papers.

Acknowledgments

We are grateful to the many friends and colleagues who contributed to the TDCOSMO collaboration over the years, including A. Agnello, T. Anguita, M.W. Auger, R.D. Blandford, V. Bonvin, C.-F. Chen, L. Christensen, T. Collett, S. Ertl, M. Gomer, S. Hilbert, L.V.E. Koopmans, C. Lemon. P.J. Marshall, S. Mukherjee, C.E. Rusu, T. Schmidt, O. Tihhonova, and A. Yıldırım. This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These observations are associated with programs #1198, 1794, and 2974. The specific observations analyzed can be accessed via https://dx.doi.org/10.17909/m1m8-rk48, https://dx.doi.org/10.17909/wjva-0750, and https://dx.doi.org/10.17909/gj7b-p622. Some of the data presented herein were obtained at Keck Observatory, which is a private 501(c)3 non-profit organization operated as a scientific partnership among the California Institute of Technology, the University of California, and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the Native Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain. Based in part on observations collected at the European Southern Observatory data obtained from the ESO Science Archive Facility. This research made use of ASTROPY, a community-developed core PYTHON package for Astronomy (Astropy Collaboration 2013, 2018), the 2D graphics environment MATPLOTLIB (Hunter 2007), EMCEE, a PYTHON implementation of an affine invariant MCMC ensemble sampler (Foreman-Mackey et al. 2013), HIERARC version 1.2.0, a package to hierarchically sample cosmological and lens population parameters (Birrer et al. 2020), LENSTRONOMY version 1.13.1 (Birrer & Amara 2018; Birrer et al. 2021), and the SQUIRREL pipeline (Shajib et al. 2025a) to streamline kinematic measurements using PPXF (Cappellari & Emsellem 2004). S.B. and X.-Y.H. acknowledge support by the Department of Physics & Astronomy, Stony Brook University. F.C. is supported in part by the Swiss National Science Foundation. F.C., D.S., and L.V. acknowledge the support of the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (COSMICLENS: grant agreement No 787886). C.D.F. acknowledges support for this work from the National Science Foundation under Grant No. AST-2407278. A.G. acknowledges support by the SNSF (Swiss National Science Foundation) through mobility grant P500PT_211034. D.L. was supported by research grants (VIL16599, VIL54489) from VILLUM FONDEN. M.M. acknowledges support by the SNSF (Swiss National Science Foundation) through mobility grant P500PT_203114 and return CH grant P5R5PT_225598. T.M. received support from NASA through the STScI grant JWST-GO-3990. E.P. was supported by JSPS KAKENHI Grant Number JP24H00221. V.M. acknowledges support from ANID FONDECYT Regular grant number 1231418, Millennium Science Initiative, AIM23-0001, and Centro de Astrofísica de Valparaíso CIDI N21. Support for this work was provided by NASA through the NASA Hubble Fellowship grant HST-HF2-51492 awarded to A.J.S. by the STScI, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. A.J.S. also received support from NASA through STScI grants HST-GO-16773 and JWST-GO-2974. D.S. acknowledges the support of the Fonds de la Recherche Scientifique-FNRS, Belgium, under grant No. 4.4503.1. M.S. acknowledges partial support from NASA grant 80NSSC22K1294. S.H.S. thanks the Max Planck Society for support through the Max Planck Fellowship. This work is supported in part by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 – 390783311. C.Y.T. was supported by NASA through the Space Telescope Science Institute (STScI) grant HST-AR-16149. T.T. acknowledges support by NSF through grants NSF-AST-1906976, NSF-AST-1836016, and NSF-AST-2407277, and from the Moore Foundation through grant 8548. K.C.W. is supported by JSPS KAKENHI Grant Numbers JP20K14511, JP24K07089, JP24H00221.

References

Appendix A: Main properties of the lenses

In Table A.1 in this appendix, we summarize the main properties of the TDCOSMO time delay lenses and of the selected SLACS and SL2S lenses.

In Fig. A.1, produced after unblinding, we compare the newly measured stellar velocity dispersion to those used by TDCOSMO-4 and Wong et al. (2024).

Fig. A.1. Comparison of the stellar velocity dispersions measurement used in this paper and those in TDCOSMO-4 and Wong et al. (2024, for WGD 2038−4008). Note that the measurement apertures and PSF resolutions differ between those presented in this paper and those in previous ones. Thus, the measurements are not expected to match exactly. For RX J1131−1231, we illustrate here the value reported in Table A.1 that corresponds to the stellar velocity dispersion measured within half the effective radius.

Appendix B: Modeling of RX J1131−1231

In this Appendix, we provide a more detailed description of the analysis of the time delay lens RX J1131−1231, the one for which spatially resolved kinematics are available from JWST and KCWI spectroscopy. One important factor is that the time-delay distance D_Δt based on the lens model by Suyu et al. (2013) depends on the inferred power-law slope γ_pl. This correlation was taken into account when combining the lens model likelihood with the likelihoods from the dynamical analysis, that constrain γ_pl but not D_Δt.

In Fig. B.1, we show the results of the fit to the RX J1131−1231 data alone, without combining with the other lenses in our hierarchical analysis. This analysis is similar to that presented by Shajib et al. (2023), except that one-dimensional velocity dispersion profiles are considered instead of the 2D maps. The figure shows that the KCWI and JWST data yield consistent results, which can be described by a single mass model and combined. The likelihoods can also be combined with the inference of γ_pl from the lens model, yielding the tighter constraints in black. We note that the contours in black are consistent with an isothermal power law (λ_int = 1, γ_pl = 2). The constraints on H₀ from this system alone are consistent with those obtained by Shajib et al. (2023).

Fig. B.1. Posterior distribution function of key parameters describing the time-delay lens RX J1131−1231, based on velocity dispersion profiles presented in this paper from the Keck-KCWI and JWST-NIRSpec data, and in combination with the lens model inferred power-law density slope by Suyu et al. (2013). The individual analysis of this system is shown to demonstrate that the three datasets are mutually consistent and can therefore be combined, and to illustrate the key degeneracies between the parameters. The inferred H0 from this system is consistent with that obtained by Shajib et al. (2023) based on the analysis of 2D kinematic maps from the KCWI. In the hierarchical analysis presented in the main body of the paper, H0 and λint are population-level parameters, jointly constrained by all the lenses in the sample.

Appendix C: Kinematic fits to the KCWI stellar velocity dispersion profiles of the SLACS lenses

In Fig. C.1 we show the velocity dispersion profiles measured for the SLACS lenses based on KCWI data and compare them with the results of the hierarchical fit.

Fig. C.1. For each SLACS lens used in our analysis, we show the observed stellar velocity dispersion profile and the best fit found in the hierarchical analysis. The data for individual lenses is described very well by our model when applied to a single lens. However, in the hierarchical analysis, they were all required to share a single population-level value of λint (Section 7) with the time delay lenses and among each other. The offsets between the data and the best-fit model are due to mismatches between the population-level value of λint and the individual λint that would be required to fit each lens separately. The offsets are reflected in the non-zero scatter in λint found for the SLACS population.

Appendix D: Author contributions and acknowledgments

Author contributions. The authors of this paper are listed in alphabetical order. Here, we specify each author’s contribution to the paper. S.B. developed the hierarchical framework, led the development of HIERARC, contributed to the analysis and writing, and oversaw the Stony Brook group. E.J.B.-G. performed the LoS analysis for DES J0408−5354 and WGD 2038−4008. M.C. contributed to the stellar kinematic measurements, dynamical modeling, and paper writing and served as the internal reviewer of the paper. F.C. established and led the time-delay campaigns, supervised the EPFL group. F.D. measured the time delay for WGD 2038−4008. C.D.F. provided suggestions for data analysis, contributed to the development of environmental analysis code, and supervised the development of one of the velocity dispersion codes as the PhD supervisor of P.M. and P.W., as well as the acquisition of imaging data used for environmental analysis of several lenses. J.F. served as the internal reviewer of the paper and led the DES and U Chicago groups. A.G. served as the internal code reviewer, participated in the development of SQUIRREL, and provided comments on the analysis and the paper. D.G. quantified uncertainties in the time delays from dark-matter substructure. X.-Y.H. managed the hierarchical analysis pipeline, led the axisymmetric correction in kinematic modeling, and participated in writing and code review. S.K. acquired, reduced, and measured KCWI spatially-resolved kinematics of the SLACS sample, reduced and measured aperture-integrated kinematics of JWST-NIRSpec data, contributed to developing the stellar kinematics analysis method used in this work, dynamical modeling, and paper writing. D.L. performed preliminary velocity dispersion measurements for DES J0408−5354 and WGD 2038−4008 from the MUSE data, and provided general indirect contributions to kinematic measurements. H.L. performed data acquisition and analysis for WGD 2038−4008 and DES J0408−5354. M.M. performed the MUSE kinematic measurements, implemented the probe combination analysis, contributed to the time-delay measurement campaigns, and contributed to the analysis and writing. T.M. planned the observations of JWST-GTO-1198 and did initial reduction of the NIRCam imaging and NIRSpec IFU data. V.M. worked on time-delay campaigns and data acquisition. P.M. measured the stellar kinematics of the SL2S sample and contributed to the development of the stellar kinematics analysis method used. E.P. worked on time delay campaigns, data acquisition, and reduction. A.J.S. performed resolved kinematic measurements of RX J1131−1231, wrote part of the paper and made a subset of the figures, developed the SQUIRREL pipeline for kinematics extraction that was used for the TDCOSMO-2025 sample, developed the REGALJUMPER pipeline for JWST data reduction, performed lens modeling of DES J0408−5354 and WGD 2038−4008, oversaw Project Dinos that provided external lens models for the SLACS and SL2S systems, and acquired a subset of the JWST data used as co-PI of JWST-GO-2974. W.S. contributed the model parameters for the systems in the SL2S sample. D.S. performed data reduction and extraction of the MUSE data, designed and contributed to some systematic tests, and conducted the internal review of this paper. A.S. contributed to the velocity dispersion measurements of the SL2S lenses. C.S. contributed to the development of the stellar kinematic measurement method used in this paper. MS carried out the observation and initial reduction of the IFU data for three of the targets and participated in discussions on how to improve the pipeline products for IFU data. S.H.S. led JWST data acquisition of RX J1131−1231, contributed to lens mass modeling, developed mass modeling software GLEE, served as internal reviewer of the paper, and oversaw the TUM/MPA group. C.Y.T. provided lens model parameters for the SLACS and SL2S lens systems. T.T. led the JWST, HST, and Keck data acquisition, contributed to the analysis and stellar kinematic measurements and methods, contributed to the hierarchical analysis and quality control, contributed to the modeling of the TDCOSMO-2025, SLACS, and SL2S lenses, contributed to the writing of the paper, served as final reader of the paper, and oversaw the UCLA group. L.V. worked on lensing degeneracies and extraction of the MUSE data of DES J0408−5354. H.W. worked on testing the preliminary NIRSpec kinematic map of RX J1131−1231 and the degeneracy between anisotropy and black hole mass. P.W. performed the LoS analysis of the SL2S Sample. D.W. verified modeling systematics by modeling WFI 2033−4723’s JWST NIRCam data, and contributed to the determination of external convergence. K.C.W. contributed to lens mass modeling and provided parameter chains for WGD 2038−4008.

All Tables

All Figures

Fig. 1. Montage of the eight lensed quasar systems in the TDCOSMO-2025 sample. In each panel, the white bar illustrates the 1″scale. The false-color images are made from two or three bands from the HST, Keck-NIRCam, and Chandra X-ray imaging, given their availability. The image for DES J0408−5354 is adapted from Shajib et al. (2020), HE 0435−1223, B1608+656, and WFI 2033−4723 from Suyu et al. (2017), PG 1115+080 and SDSS J1206+4332 from Wong et al. (2020), RX J1131−1231 from Shajib et al. (2024, credits: NASA/ESA/Chandra), and WGD 2038−4008 from Shajib et al. (2022).

In the text

Fig. 2. JWST-NIRSpec spectra and kinematic fits for RX J1131−1231. Top left: The six annuli (black contours), from which summed spectra are extracted, are illustrated on top of the NIRSpec white-light image. The regions around the satellite galaxy and the closest quasar image are excluded. The white bar represents 1″ scale, and the North and East directions are pointed with emerald and yellow arrows, respectively. Bottom left: The measured vrms in the six annuli. The horizontal error bars show the annulus widths, and the vertical error bars show the combined statistical and systematic uncertainty for each measurement. The measured values have 0.66% covariance on average. Second and third columns: The six panels show the integrated spectra in each annulus (gray boxes) and the kinematic fit (red line). The height of the gray box represents the nominal uncertainty levels estimated by the JWST pipeline, and the size of the vertical error bars represents the total boosted uncertainty levels to achieve χred2 = 1 for each fit. The vertical orange lines mark the wavelengths of the Calcium triplet lines from the lens galaxy, and the vertical blue lines mark the wavelengths of the Hα, [N II], and [S II] lines from the quasar host galaxy.

In the text

Fig. 3. JWST-NIRSpec integrated spectra (gray bars) and kinematic fits (red lines) for five TDCOSMO-2025 lenses. The height of the gray boxes represents the nominal uncertainty levels estimated by the JWST pipeline, the size of the vertical error bars represents the total boosted uncertainty levels to achieve χred2 = 1, and the width of the boxes represents the wavelength-pixel size. The x-axis shows the lens rest-frame wavelength. The gray-shaded vertical bands were masked during the fits due to contamination by the lensed quasar features.

In the text

Fig. 4. Measured values of vrms for RX J1131−1231 in radial annuli from the JWST NIRSpec (left-hand panel, the same ones shown in Fig. 2) and Keck KCWI (middle panel). Both of these measurements marginalize over the same choices of template libraries, namely the Indo-US and the XSL DR3, in addition to separate choice combinations for polynomial orders and fitted wavelength range. The arrows in these panels illustrate the size of the PSF FWHM for each case ($0\stackrel{″}{.}15$ for JWST-NIRSpec and $0\stackrel{″}{.}96$ for Keck-KCWI): the higher resolution of the JWST-NIRSpec data enabled the identification of the sharp rise of the velocity dispersion profile in the center. The right-hand panel illustrates the 1D surface brightness profile in optical (blue line) and in the IR (red line). The optical light profile was extracted from double Sérsic fitting by Shajib et al. (2023) from a large cutout that encapsulates the full extent of the galaxy in the HST F814W imaging (which corresponds to a pivot wavelength of 6208 Å in the lens rest-frame). The IR light profile was extracted from the double Sérsic light profile fitted as part of the lens modeling done with a 2D image obtained from the NIRSpec datacube by summing within the wavelength range 8700–8800 Å in the lens rest-frame. The surface brightness profiles are normalized by the corresponding amplitudes at $R=1\stackrel{″}{.}3$ (vertical dashed line). For kinematic modeling of the NIRSpec kinematics, we use a stitched light profile that transitions from the IR light profile at $R\le 1\stackrel{″}{.}3$ to the optical light profile shape at $R>1\stackrel{″}{.}3$.

In the text

Fig. 5. VLT-MUSE spectra and kinematics fits with PPXF for DES J0408−5354 (left-hand panel) and WGD 2038−4008 (right-hand panel). The gray rectangles illustrate the data with the height representing the 1σ uncertainty and the width representing the wavelength-pixel size. The red line illustrates the best-fit model. The principal spectral features probing the kinematics are the Ca H & K absorption lines marked with vertical orange lines. The gray-shaded region on the left-hand panel represents the wavelength range excluded from the fit.

In the text

	Fig. 6. Probability distribution function of external convergence κext for the quality samples of the SLACS and the SL2S lenses. Dashed lines are κext of individual lenses, and solid lines are their combination. The median of κext for both samples is shown in the figure legend.
In the text

Fig. 7. Properties of selected lens samples that enter the cosmology inference. The sample selection is only based on (i) lens and source existing spectroscopic redshift measurements; (ii) lensing analysis based on imaging data; (iii) LoS analysis. From left to right, the quantities are: lens redshift zlens, background source redshift zsource, stellar velocity dispersion of the deflector σlosap, half-light radius of the deflector Reff, Einstein radius θE, Reff/θE, lensing power-law slope γpl, lensing stellar velocity dispersion σSIS, relative over-density ζ1/R, apparent axis ratio of light profile of the deflector qlight, and lensing information ℐ calculated under the choice of scaling parameters from Sheu et al. (2025). We demonstrate for the selected lenses, there is no significant difference from sample to sample in the parameter spaces of this plot, except for the distribution in redshifts and the difference in effective radius expected from the evolution of the size-mass relation (van Dokkum et al. 2010), see Section 6.6 for more information.

In the text

Fig. 8. Posterior distributions for a subset of the parameters described in Table 3. Blue: constraints from the eight TDCOSMO-2025 lenses, including JWST kinematic data for RX J1131−1231. Orange: constraints from the TDCOSMO-2025 sample + 11 SLACS lenses with KCWI IFU kinematics data. Pink: constraints from TDCOSMO-2025 + SLACS + SL2S sample. All contours include Ωm constraint coming from the Pantheon+ sample (Scolnic et al. 2022); Brout:2022. Priors are described in Table 3. The posteriors distribution of H0 and λint, 0 were blinded during the analysis. Results in combination with the DES-SN5YR likelihood are very similar (see Table 6) and are not shown for conciseness.

In the text

Fig. 9. Posterior distribution of H0 and Ωk in open ΛCDM cosmology. The black contours show the constraints from TDCOSMO alone, including the TDCOSMO-2025, SLACS, and SL2S lensing datasets. The gray contours show the constraints from Planck. The orange contours correspond to the combination of Planck and TDCOSMO. The two datasets exhibit a tension exceeding 3σ, indicating that this joint constraint should be interpreted with caution. The contour levels represent the 1σ and 2σ constraints.

In the text

Fig. 10. Posterior distribution of H0 and w in flat wCDM cosmology. The black contours show the constraints from TDCOSMO alone, including the TDCOSMO-2025, SLACS, and SL2S lensing datasets. The gray contours show the constraints from Planck. The combination of TDCOSMO with the Planck, DESI BAO, and Pantheon+ is shown with colored contours, as described in the figure legend. The contour levels represent the 1σ and 2σ constraints.

In the text

	Fig. 11. Posterior distribution of H0, w0 and wa in flat w0waCDM cosmology. The combination of TDCOSMO (including the TDCOSMO-2025, SLACS, and SL2S lensing datasets) with the DESI-DR2 BAO and Pantheon+, and the combination of both external datasets, is shown with colored contours described in the figure legend. The contour levels represent the 1σ and 2σ constraints.
In the text

	Fig. 12. Posterior distribution of H0, Ωm and w0 in flat wϕCDM cosmology. The combinations of TDCOSMO (including the TDCOSMO-2025, SLACS, and SL2S lensing datasets) with the DESI-DR2 BAO and Pantheon+ datasets are shown with colored contours, as described in the figure legend. The contour levels represent the 1σ and 2σ constraints.
In the text

Fig. 13. Comparison of the H0 measurements by the TDCOSMO collaboration (and its predecessor H0LiCOW) in chronological order. The measurements at the top by Wong et al. (2020) and Millon et al. (2020b) are “assertive” in the definition of Treu et al. (2022). They are based on power-law and composite models to describe the mass density profile of the deflector galaxies, thus implicitly breaking the MSD and obtaining ∼2% precision. The TDCOSMO-4 results from Birrer et al. (2020) are “conservative” (Treu et al. 2022). They are based on the same data but obtain larger uncertainties by introducing the internal MSD parameter λint, and constraining it with unresolved stellar kinematics. The Birrer et al. (2020) measurement in combination with SLACS is shown as a dashed line for historical purposes but should not be used anymore, as it was later discovered that the stellar velocity dispersions based on low SNR (∼9 Å−1) spectra from SDSS suffer from systematic errors and covariance (Knabel et al. 2025a). The new measurements presented in this work are shown in the bottom panel. They are “conservative” in terms of uncertainties and constrain the MSD using new stellar kinematics based on high SNR JWST-NIRSpec, VLT-MUSE, and Keck-KCWI spectra, as well as an improved methodology (Knabel et al. 2025b).

In the text

Fig. 14. Illustration of the current state of the “Hubble tension” in flat ΛCDM. Only a selection of recent measurements is shown for the early- and late-Universe probes, in order to avoid overcrowding the plot (for an exhaustive list, see Di Valentino et al. 2025). Among the early-Universe probes, we show measurements from the CMB (Planck and ACT-DR6; Planck Collaboration VI 2020; Louis et al. 2025), and baryonic acoustic oscillations (BAO) plus Big Bang nucleosynthesis (BBN) (DESI-Collaboration 2025). Among the late-Universe probes, we show the Cepheid+SN measurement from the SH0ES team (Breuval et al. 2024), the tip of the red giant branch (TRGB) plus SN measurement from the CCHP team (Freedman et al. 2025), the TRGB plus surface brightness fluctuation (SBF) measurement by Jensen et al. (2025), along with the TDCOSMO-2025 and TDCOSMO-2025+SLACS+SL2S time-delay cosmography measurements in combination with the Pantheon+ likelihood (Scolnic et al. 2022); Brout:2022.

In the text

Fig. A.1. Comparison of the stellar velocity dispersions measurement used in this paper and those in TDCOSMO-4 and Wong et al. (2024, for WGD 2038−4008). Note that the measurement apertures and PSF resolutions differ between those presented in this paper and those in previous ones. Thus, the measurements are not expected to match exactly. For RX J1131−1231, we illustrate here the value reported in Table A.1 that corresponds to the stellar velocity dispersion measured within half the effective radius.

In the text

Fig. B.1. Posterior distribution function of key parameters describing the time-delay lens RX J1131−1231, based on velocity dispersion profiles presented in this paper from the Keck-KCWI and JWST-NIRSpec data, and in combination with the lens model inferred power-law density slope by Suyu et al. (2013). The individual analysis of this system is shown to demonstrate that the three datasets are mutually consistent and can therefore be combined, and to illustrate the key degeneracies between the parameters. The inferred H0 from this system is consistent with that obtained by Shajib et al. (2023) based on the analysis of 2D kinematic maps from the KCWI. In the hierarchical analysis presented in the main body of the paper, H0 and λint are population-level parameters, jointly constrained by all the lenses in the sample.

In the text

Fig. C.1. For each SLACS lens used in our analysis, we show the observed stellar velocity dispersion profile and the best fit found in the hierarchical analysis. The data for individual lenses is described very well by our model when applied to a single lens. However, in the hierarchical analysis, they were all required to share a single population-level value of λint (Section 7) with the time delay lenses and among each other. The offsets between the data and the best-fit model are due to mismatches between the population-level value of λint and the individual λint that would be required to fit each lens separately. The offsets are reflected in the non-zero scatter in λint found for the SLACS population.

In the text