© ESO 2020

1. Introduction

The time-delay method applied to gravitationally lensed quasars (Refsdal 1964) provides a perhaps unrivalled combination of high sensitivity to the Hubble constant H₀, and minimal dependence on the other cosmological parameters, while relying only on well known physics (i.e., gravity). These qualities make this method particularly important in the present context, where there is growing evidence for tension in H₀ measurements using cosmological probes based on the early Universe and the late Universe (Verde et al. 2019). The power of the method in providing reliable H₀ measurements depends on three main factors: 1- precise time-delay measurements between multiple images of the background source, 2- well constrained models of the dominant primary and nearby lens galaxies, and 3- an estimate of the combined lensing effect of all the mass along the line of sight up to the redshift of the lensed quasar.

Precise and accurate time-delay measurements are available, for example, from the COSMOGRAIL collaboration, using long-term photometric monitoring of selected lensed quasars (e.g., Courbin et al. 2018; Bonvin et al. 2018, 2019). The precision and accuracy of the COSMOGRAIL technique have been verified via a blind time-delay challenge (Dobler et al. 2015; Liao et al. 2015; Bonvin et al. 2016). The time-delays were then used to constrain cosmological parameters with detailed modeling of the potential well of the lens using the constraining power of sharp Hubble Space Telescope (HST) images (e.g., Suyu et al. 2010, 2014; Wong et al. 2017; Birrer et al. 2019; Rusu et al. 2020) or Keck AO imaging (e.g., Chen et al. 2019). The measured stellar kinematics of the lensing galaxy were used to mitigate the impact of well-known lensing degeneracies on the cosmological inference (e.g., Treu & Koopmans 2002). Finally, multi-band wide-field imaging and/or spectroscopy (e.g., Rusu et al. 2017; Sluse et al. 2019) was used to constrain the combined lensing effect of the line-of-sight objects and large-scale structures in a statistical way (Greene et al. 2013; Rusu et al. 2017). Tihhonova et al. (2018) also show that these estimates of the line-of-sight effects are compatible with the ones obtained with weak gravitational lensing.

Adopting these data and methodology, the H0LiCOW collaboration (Suyu et al. 2017) is analyzing a sample of lenses suitable for high-precision H₀ measurements. The latest results based on six systems are summarized by Wong et al. (2020). We stress that the H0LiCOW results are obtained through blind analyses, in the sense that the mean value of all the observed cosmological parameters is hidden to the investigators until the analysis is complete and the papers have been written¹. The goal of this procedure is to avoid conscious or unconscious bias from the experimenters. We note that the six measurements that have been published thus far are statistically consistent with each other, in the sense that the scatter between the measurements is as expected from the estimated uncertainties. This means that if there are any unknown uncorrelated sources of error, those are subdominant with respect to the ones currently considered.

The resulting value of the Hubble constant in a flat Λ cold dark matter (ΛCDM) universe, $H_0=73.3_{-1.8}^{+1.7}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$ (2.4% precision), is 3σ higher than the early-Universe results (Planck Collaboration VI 2020), adopting the same ΛCDM cosmological model, and is in very good agreement with other independent local measurements (e.g., Riess et al. 2019). When combined with completely independent results from other local measurements of H₀, the tension with the early-Universe probes range between 4 and 6σ (Verde et al. 2019), depending on the combination of probes. Very recently, Pandey et al. (2020) also carried out statistical tests independent of any underlying cosmology, showing that the distances measured with strong lensing time delays and with supernovae, which are both local but independent measurements, are fully compatible (see also Wojtak & Agnello 2019). Although they cannot exclude that supernovae and lenses share exactly the same systematics, these systematic biases would also have to be preserved across redshift, which seems unlikely.

The blind analysis of a seventh lens system using very similar methods for the lens modeling, time-delay measurement, external convergence estimation and kinematics modeling to those adopted by H0LiCOW has recently been published by the STRIDES collaboration (Shajib et al. 2020). This work finds 74.2${}_{-3.0}^{+2.7}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$, in agreement with the H0LiCOW result (an independent analysis adopting a different modeling software is currently under way). This most recent system is particularly interesting since it has two sets of multiple images at different redshifts, which help break some of the degeneracies, and results in the most precise individual measurement so far. In order to make further progress in this important area, members of the COSMOGRAIL, H0LiCOW, SHARP and STRIDES collaborations interested in time-delay cosmography of lensed quasars have decided to join forces with other scientists and form a new “umbrella” collaboration named TDCOSMO² (Time-Delay COSMOgraphy).

The high statistical significance of the tension between early and late Universe probes has prompted two lines of investigation. On the one hand, theorists have been trying to find ways to reconcile the measurements by considering models beyond the standard ΛCDM one (e.g., Knox & Millea 2020). On the other hand, observing teams have been focusing on increasing the precision of each method while carrying out tests of potential systematic uncertainties to ensure that the tension is real. After all, “extraordinary claims require extraordinary evidence”.

In this work, the first by the TDCOSMO collaboration, we explore a number of potential systematic uncertainties that may affect the time-delay cosmography method, after reviewing its methodology and implementation by TDCOSMO in Sect. 2 and the inference procedure in Sect. 3. First, in Sect. 4 we explore potential biases introduced by systematic uncertainties in the modeling and measurement of the deflector stellar velocity dispersion. Second, in Sect. 5, we study uncertainties in the modeling of the line-of-sight contribution. Third, in Sect. 6, we address the long standing issue of the mass-sheet degeneracy and the flexibility of lensing models. It is very well known that assumptions must be made on the form of the main deflector mass distribution to break the mass-sheet degeneracy. As many authors have pointed out (Falco et al. 1985; Read et al. 2007; Schneider & Sluse 2013; Xu et al. 2016; Sonnenfeld 2018; Kochanek 2020), if the models adopted are insufficiently flexible, the resulting uncertainties are underestimated and potentially biased. Section 7 offers a summary and conclusions.

We address these three sources of potential systematic uncertainties using a combination of observational tests and simulations. We stress that a full simulation of the observational setup and lens modeling procedure is needed if one wants to obtain quantitative estimates of the uncertainties. Previous works (Schneider & Sluse 2013; Sonnenfeld 2018; Kochanek 2020) were based on idealized, often spherical models. Those are useful to gain intuition of the problem, but by their very nature cannot provide quantitative estimates due to the extreme approximation and the limited information utilized to constrain them, often just the Einstein Radius and an integrated velocity dispersion. The only way to obtain a faithful estimate of the uncertainties is to reproduce the measurement using the same amount of information (thousands of pixels from imaging, multiple time-delays, stellar kinematics) and modeling techniques. The simulated dataset shown in this paper are produced using the pipeline developed by Ding et al. (2017a,b, 2018). In order to isolate and quantify the uncertainties associated with the lens mass modeling procedure, the simulated data consist of high-resolution images of lens systems comparable to the real observations, high-precision time delays (higher precision than those of real lenses so that the time-delay uncertainties are subdominant compared to the modeling uncertainties that we aim to quantify), and do not include the line-of-sight structures. The lens fitting procedure that we use to analyze these simulated data resembles as closely as possible that of the TDCOSMO collaboration.

2. Background

2.1. Time-delay cosmography and the mass-sheet degeneracy

Time delays in gravitationally lensed quasars provide a direct measurement of the so-called “time-delay distance”, which is a combination of angular diameter distances to the source, D_s, to the deflector, D_d, from the deflector to the source, D_ds, and the redshift of the deflector z_d:

$$\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}$$(1)

(Refsdal 1964; Schneider et al. 1992; Suyu et al. 2010).

This quantity is related to the relative time delay between two multiple images A and B, Δt_AB, by:

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

where θ is the image position on the plane of the sky, β is the (unobservable) source position, c is the speed of light and ψ is the lensing potential which is defined such that the deflection angle α(θ) is given by α(θ) ≡ ∇ψ(θ). From Eq. (2), we see that D_Δt depends on the geometry of the lensed system and on the potential well of the lensing galaxy. The mass profile is expressed as a dimensionless surface mass density, κ(θ), called the convergence. It is related to how the light beams from the source are stretched or squeezed, leading to an apparent (de)magnification and can be expressed as half of the Laplacian of the lensing potential:

$$\begin{array}{c}\hfill \kappa (\mathit{\theta })=\frac{1}{2}{\mathrm{\nabla }}^2\psi (\mathit{\theta }).\end{array}$$(3)

We can also define the Fermat potential ϕ (Schneider 1985; Blandford & Narayan 1986) as

$$\begin{array}{c}\hfill \varphi (\mathit{\theta })\equiv \frac{{(\mathit{\theta }-\mathit{\beta })}^2}{2}-\psi (\mathit{\theta }).\end{array}$$(4)

Using this definition, Eq. (2) reduces to

$$\begin{array}{c}\hfill \mathrm{\Delta }{t}_{\mathrm{AB}}=\frac{D_{\mathrm{\Delta }t}}{c}[\varphi ({\mathit{\theta }}_{\mathbf{A}})-\varphi ({\mathit{\theta }}_{\mathbf{B}})]\equiv \frac{D_{\mathrm{\Delta }t}}{c}\mathrm{\Delta }{\varphi }_{\mathrm{AB}},\end{array}$$(5)

where Δϕ_AB is the difference of Fermat potentials at the positions of the multiple images. Based on the multipole decomposition of the gravitational potential, Kochanek (2002) shows that the time-delay distance, D_Δt, depends on the mean surface density ⟨κ⟩ at the Einstein radius θ_E, specifically over the annulus defined by image positions.

An inherent limitation of the lensing models to infer D_Δt is the so called Mass-Sheet Transformation (MST, e.g., Falco et al. 1985) and its generalization (Saha 2000; Saha & Williams 2006; Liesenborgs & De Rijcke 2012; Schneider & Sluse 2014; Wagner 2018; Wertz et al. 2018). The MST transforms the projected mass distribution and the source plane position according to:

$$\begin{array}{cc}& \kappa (\theta )\to {\kappa }_{\lambda }(\theta )=(1-\lambda )\times \kappa (\theta )+\lambda ,\hfill \\ \hfill & \beta \to \beta \prime =\lambda \beta ,\hfill \end{array}$$(6)

where β is the (unknown) source position on the sky prior to lensing. In other words, one can add a mass sheet to any model and apply a scaling factor, λ, without changing the lensing observables except the time delays and therefore the inferred cosmology.

The time-delay distance given by any model is affected by MST as follows :

$$\begin{array}{c}\hfill D_{\mathrm{\Delta }t}^{\mathrm{MST}}=D_{\mathrm{\Delta }t}^{\mathrm{true}}\times (1-\lambda ).\end{array}$$(7)

In the TDCOSMO analyses, this scaling factor λ is identified with the external convergence factor κ_ext which accounts for the contribution of all the mass along the line of sight (LOS). It is estimated independently from the lens modeling by comparing the relative number of galaxies weighted by physically relevant priors such as the distance to the lens, the stellar mass and the redshift in a large aperture around the strong lens system with simulated LOS extracted from numerical simulations with similar statistical properties (Rusu et al. 2017). Alternatively, the external convergence can be estimated from a weak lensing analysis (Tihhonova et al. 2018).

In addition to the MST above due to external mass sheets (i.e., external mass structures that do not affect the stellar dynamics of the foreground lens galaxy), MST can also manifest itself approximately as a change in the radial mass profile of the foreground lens galaxy. We describe this as an “internal” mass sheet. To mitigate the effects of the internal mass sheet, we consider different families of models and further use kinematic measurements of the foreground lens that provide additional constraints on the lens mass models. In particular, the goodness of fit to the kinematic data, especially spatially-resolved lens stellar velocity dispersion, allows us to distinguish between otherwise degenerate lensing mass models (e.g., Yıldırım et al. 2020).

The lens stellar velocity dispersion of the foreground lens galaxy allow the inference of the angular diameter distance, D_d, to the lens, in addition to the time-delay distance (Paraficz & Hjorth 2009; Jee et al. 2015, 2019). The inference of D_d depends on the anisotropy of stellar orbits (Jee et al. 2015), but this additional distance measurement provides more leverage on constraining cosmological models (Jee et al. 2016; Shajib et al. 2018).

2.2. Two-distance inference

In the most recent analysis of SDSS J1206+4332, PG 1115+080, RX J1131−1231, B1608+656 and DES J0408−5354 (Birrer et al. 2016, 2019; Chen et al. 2019; Shajib et al. 2020; Wong et al. 2020), the time-delay distance D_Δt and the angular diameter distance to the lens D_d are jointly inferred. Following the method developed in Birrer et al. (2016), the luminosity weighted LOS velocity dispersion within an aperture 𝒜 of the main deflector σ_v can be expressed as:

$$\begin{array}{c}\hfill {\sigma }_v^2=(1-{\kappa }_{\mathrm{ext}})\frac{D_{\mathrm{s}}}{D_{\mathrm{ds}}}{c}^{2}J({\xi }_{\mathrm{lens}},{\xi }_{\mathrm{light}},{\beta }_{\mathrm{ani}}),\end{array}$$(8)

where ξ_lens is the set of all parameters contained in the lens mass model, ξ_light is the parameter of the light models and J is a function that captures all dependencies on the modeling parameters and the anisotropy profile β_ani. Using Eqs. (1), (5) and (7), we have :

$$\begin{array}{c}\hfill \frac{D_{\mathrm{d}}{D}_{\mathrm{s}}}{D_{\mathrm{ds}}}=\frac{c\mathrm{\Delta }{t}_{\mathrm{AB}}}{(1+z_d)(1-{\kappa }_{\mathrm{ext}})\mathrm{\Delta }{\varphi }_{\mathrm{AB}}({\xi }_{\mathrm{lens}})}·\end{array}$$(9)

Combining Eqs. (8) and (9), we obtain an expression for the angular diameter distance to the lens which is independent of the external convergence:

$$\begin{array}{c}\hfill D_{\mathrm{d}}=\frac{c^3\mathrm{\Delta }{t}_{\mathrm{AB}}J({\xi }_{\mathrm{lens}},{\xi }_{\mathrm{light}},{\beta }_{\mathrm{ani}})}{(1+z_{\mathrm{d}}){\sigma }_v^2\mathrm{\Delta }{\varphi }_{\mathrm{AB}}({\xi }_{\mathrm{lens}})}·\end{array}$$(10)

We immediately see that the angular diameter distance D_d varies as $\frac{1}{{\sigma }_{\mathit{v}}^2}$. The dependence of D_d to a change in the measurement of σ_v can therefore be computed analytically :

$$\begin{array}{c}\hfill \frac{\mathrm{d}{D}_{\mathrm{d}}}{D_{\mathrm{d}}}=-2\frac{\mathrm{d}{\sigma }_v}{{\sigma }_v},\end{array}$$(11)

whereas D_Δt is left unchanged when varying the velocity dispersion. The final H₀ measurement is obtained by combining these two distance measurements. As a consequence, the importance of the velocity dispersion in the final H₀ value depends on the relative precision between the angular diameter distance and the time-delay distance, and on the mapping between the parameters. The D_Δt measurement is typically more constraining of H₀ than D_d given the current observational data. Future observations with spatially resolved kinematics are expected to improve substantially the D_d constraints (Yıldırım et al. 2020).

Two of the lens systems in the TDCOSMO sample, HE 0435−1223 and WFI 2033−4723, have nearby massive perturbing galaxies at a different redshift from the strong lensing galaxy, and thus required multi-lens-plane mass modeling. The single-lens-plane Eqs. (8)–(9) are thus not directly applicable, given the additional angular diameter distances involved in the multiple lens planes. Nonetheless, the mass model of the lens galaxy can still be used to predict the velocity dispersion to compare to the measured value, so the kinematic measurement can be used to further constrain the mass model. It turns out that an effective time-delay distance could be derived for these two lens systems, but the inference of D_d accounting for the multi-lens planes is deferred to future work.

2.3. The current TDCOSMO model families

The collaborations within TDCOSMO currently consider two classes of models (composite and power-law), to reconstruct the mass distribution of the main lens, with the exception of the first system analyzed B1608+656 (Koopmans et al. 2003; Suyu et al. 2010). B1608+656 was modeled only using a power-law, as Suyu et al. (2009) showed that deviations to a smooth potential using pixellated corrections were negligible. The fact that the corrections are so small, even though the deflector in this complex lens is an obvious merger between two galaxies, is a remarkable indication of the degree of smoothness of the overall gravitational potential. This is also supported by the analysis of extended rings used to detect substructures in lenses through their impact on the smoothness of Einstein rings. Aside from specific features arising from well-identified substructures in any given lens, no statistically significant correction to simple parametric lens models is found by Vegetti et al. (2014).

For the above reasons, the TDCOSMO analyses consider purely analytical lens models with sufficient degrees of freedom to catch a broad range of observables given current imaging capabilities with HST or adaptive optics. More specifically, the TDCOSMO analyses considers elliptical power-law and composite models, with the addition of external shear.

2.3.1. Power-law model

Power-law models have a constant projected mass slope over the entire profile. The convergence of the power-law elliptical mass distribution (Barkana 1998) is described by :

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

where γ is the slope of the profile, q_m is the axis ratio of the elliptical profile and θ_E is the Einstein radius. The coordinate system is defined such that the θ₁ and θ₂ coordinates are along the major and minor axis respectively. The cored power-law profile is a natural extension of this model which introduces an additional free parameter, namely the core radius in the center of the profile θ_c and is defined as :

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

This profile has therefore a shallower slope in the center to reproduce the core of galaxies. A complete description of this mass model can be found in Barkana (1998). Although not used by the TDCOSMO collaboration, except in the analysis of RX J1131−1231 by Suyu et al. (2014) who found negligible core size, we tested cored power-law profiles on simulated lenses in Sect. 6.

2.3.2. Composite model

The second family of mass models used by the TDCOSMO collaboration are the so-called composite models, which consist of baryonic matter and dark matter components. For the dark matter, a Navarro–Frenk–White (NFW) profile is used. The spherical NFW density distribution is given by:

$$\begin{array}{c}\hfill {\rho }_{\mathrm{NFW}}(r)=\frac{{\rho }_{\mathrm{s}}}{(r/r_{\mathrm{s}}){(1+r/r_{\mathrm{s}})}^2},\end{array}$$(14)

where r_s is the scale radius and ρ_s is a normalization factor (Navarro et al. 1997). For the baryonic component, the TDCOSMO collaboration adopts the Chameleon profile, which is the difference between two singular isothermal ellipsoids and closely mimics a Sérsic profile. A complete description of this model can be found in Dutton et al. (2011) and Suyu et al. (2014). This family of mass model allows more flexible mass distribution than power-law models since the slope of the projected mass profile is not constant over the whole lens galaxy.

3. Inference procedure and limitations of toy models

The next step required to derive a H₀ measurement from the data is a statistical inference. The collaborations contributing to TDCOSMO adopt a Bayesian framework and compute the posterior probability distribution function of all the cosmological and nuisance parameters given the data.

The imaging and spectroscopic data contain huge amounts of information, well beyond the position of the quasar images. Setting aside the line of sight, which is constrained independently, the main sources of constraints for the main deflector(s) mass models are: the pixels of the high resolution images (of order 10⁴); independent time delays (up to three for a quad); stellar velocity dispersion of the main deflector and nearby perturbers, if present. The inference required to extract all the information from the data is computationally very intensive. Taking into account the need to explore multiple and flexible models to marginalize over modeling choices, the TDCOSMO analysis required up to a million CPU hours per lens.

In the recent past, simplified toy models, that is, models in which either (i) the lens systems are not simulated with sufficient complexity, or (ii) the inference procedure does not exploit the full information content, have been used to investigate systematic uncertainties in time-delay cosmography (Schneider & Sluse 2013; Sonnenfeld 2018; Kochanek 2020). These models are certainly a useful illustration, and it is encouraging that they conclude that a precision within the range 3–10% can be reached with their simplified approach and limited constraints. However, owing to their limitations, those models cannot provide the quantitative answers that are needed to understand whether there are biases at the 2% level, which is the current achievement of time-delay cosmography. Chief among the limitations of previous works is the use of spherical models. Spherical models are inherently inappropriate to model quads (e.g., Kochanek 2006), because they cannot even produce four images and thus are intrinsically less constrained by the data than observed quads.

The bulk of the lensing information comes from the radial extent and surface brightness distribution of the lensed images, which constrains directly the radial dependency of the mass distribution, the key parameter driving the inference of H₀. Toy models neglect this information (e.g., Kochanek 2020), and are mostly spherical and constrained solely by the position of the quasar images spanning just 10% on either side of the Einstein radius. Furthermore, they are constrained only by the positions of the multiple images of the quasars and not using the full information content of the lensed host galaxy, often amounting to thousands of high signal-to-noise ratio pixels (see Sect. 6 and Appendix B for details). These constraints would have no way to detect significant departures from a power law for example, which could instead be detected in real-life cases as variations in the distortion of the images spanning a much larger significant radial range. Indeed, most of the HST data used in time-delay cosmography display prominent Einstein rings, spanning several tenths of arcseconds radially. In other words, the radial width of the ring is significant compared with the Einstein radius itself, hence constraining the potential well radially. This is clearly illustrated with the case of RX J1131−1231 in, for example, Suyu et al. (2014). In addition, toy models typically condense the information in a few parameters and thus cannot realistically explore the degeneracies between true model parameters and how uncertainties in the actual data translate into inference.

Last but not least, it is crucial to have the ability to assess the goodness of the models, in both absolute and relative terms. This is to our knowledge the most powerful way to establish whether the chosen parametrization is an appropriate description of the data. The power of goodness of fit estimates depends on the realism and information content of the models. Models that are based on key summary statistics are able to use goodness of fit only on those statistics. More realistic models that produce and fit, for example, image positions will have a few more observables to assess goodness of it. Models that produce the full surface brightness distribution of the lenses source(s) and other observables will have access to many more observables, and thus have significantly larger power for model exploration and selection.

4. Influence of kinematics data on H₀ measurement

One important ingredient to mitigate the impact of the MST is to use the kinematics of the deflector as an independent mass estimator (Treu & Koopmans 2002; Koopmans 2004), since within a cosmological model D_d and D_Δt are related to each other. So far, the central stellar velocity dispersion integrated within an aperture, σ_v, has been used even though additional and substantial gains can be obtained by including spatially resolved information that helps break the mass-anisotropy degeneracy (Barnabè et al. 2011; Czoske et al. 2012; Shajib et al. 2018; Yıldırım et al. 2020).

The inference of the Hubble constant is driven by a combination of observables, including the extended images used in the lens model, multiple time delays if available, and kinematic information. Thus, the dependency of H₀ on kinematics data defined by

$$\begin{array}{c}\hfill \xi \equiv \frac{\delta H_0/\overline{H_0}}{\delta {\sigma }_v/\overline{{\sigma }_v}}\end{array}$$(15)

cannot be estimated with simple dimensional arguments or toy models, but needs to be computed by repeating the inference while varying the input kinematics data. The result will depend on the details of the analysis as well as on the relative quality and constraining power of the kinematic and nonkinematic data, and on how the D_Δt − D_d plane maps into H₀ as a result of the deflector and source redshifts. Each of these factors varies from lens to lens as we show below and thus cannot be simply derived from a toy model and generalized to every lens.

4.1. The TDCOSMO analysis and its sensitivity to the measured velocity dispersion

Simple models such as the Singular Isothermal Sphere (SIS) models, can have a very strong dependency on the velocity dispersion. This dependency could be on the order of ξ ∼ 1, which means that a 1% change in the velocity dispersion σ_v leads roughly to a 1% change in H₀. The high sensitivity of SIS mass models to a change in the velocity dispersion arises from the fact that they have only one free parameter (the normalization). If galaxies were all SIS, then such a high sensitivity would allow us to better constrain the mass model through more precise and accurate kinematic measurements.

In this section we show that the TDCOSMO measurements, which use models more flexible than SIS and constrain them with a wealth of data, are less sensitive to the kinematics information than SIS. In order to quantify how the error on σ_v propagates into H₀, we recomputed the posterior distributions for D_Δt and D_d after changing arbitrarily the median value for our σ_v distribution. We perform the test for four values of the shift, that is, δσ_v/σ_v = ± 5% and δσ_v/σ_v = ±10%, for each individual lens in the TDCOSMO sample as well as for the joint H₀ inference. Throughout this section, the H₀ inference was performed in flat ΛCDM cosmology with a uniform prior on Ω_m ∈ [0.05, 0.5].

Figure 1 summarizes the results, where we define $\overline{H_0}$ and $\overline{{\sigma }_{\mathit{v}}}$ as the inferred H₀ value of the system and its measured aperture velocity dispersion. The models used in Fig. 1 include both composite and power-law mass models³ combined according to the standard procedure described in previous papers (e.g., Suyu et al. 2014; Chen et al. 2019; Birrer et al. 2019; Rusu et al. 2020). We first discuss in this section the general trend between σ_v and H₀ for the combination of the two model families. Then, we discuss the specifics of each model family separately.

Fig. 1. Sensitivity of the inferred Hubble constant as a function of fractional change in the measured lens velocity dispersion, σv (see Eq. (15)). Each color corresponds to one of the seven strong lens systems of the current TDCOSMO sample. The dotted lines display the best linear fit to the data. The joint inference performed on the seven lenses is shown in black. The error bars correspond to the 16th and 84th percentile of the posterior distributions. The two bottom panels show the sensitivity of H0 to a change in the measured lens velocity dispersion for power-law (left) and composite (right) models independently. The sensitivity of the joint inference, ⟨ξ⟩ is indicated on each panel.

The slope ξ quantifies the sensitivity of the inferred H₀ value to a change in velocity dispersion. It is computed by performing a linear regression to the points (Table 1). We observe large variation of measured slopes from object to object. However, for the full sample, the joint H₀ inference leads to a mean sensitivity of ⟨ξ⟩ = 0.07 ± 0.02. In other words, a systematic increase (decrease) of 10% on the velocity dispersion increases (decreases) H₀ by approximately 0.7%.

PG 1115+080 and DES J0408−5354 differ from the other lenses with a slightly negative slope of ξ = −0.04 ± 0.01 and ξ = −0.01 ± 0.01 respectively. For the other lenses, increasing the velocity dispersion leads to a smaller angular diameter distance D_d and therefore to a higher H₀ (Eq. (11)). This behavior could be explained for DES J0408−5354 as this lens is a complex system with several sources located at two different redshifts. Thus, the reduced dependency on velocity dispersion could be due to the extraordinary azimuthal and radial extent of the lensing information, and the fact that multiple redshift sources might help limit the effects of MST. In this regime, the kinematics information only brings very limited constraints on the mass model. The measurement of H₀ is therefore almost insensitive to the kinematics.

In the case of PG 1115+080, the time-delay distance D_Δt, which does not depend on the kinematics data, has a much larger constraining power on H₀ than the angular diameter distance D_d. As a result, PG 1115+080 is also almost insensitive to the velocity dispersion. The same effect explains, to a lesser extend, the low sensitivity of RX J1131−1231. We note that SDSS J1206+4332 has the largest sensitivity to a change in σ_v, with an increase of 10% in velocity dispersion leading to an increase of H₀ by 4.2%. We interpret this as the effect of D_Δt being less well constrained by the lensing data on their own. The more limited lensing constraints with respect to other systems are probably because this is the only doubly imaged quasar in the sample – all the others are quadruply imaged. Last but not least, we note that SDSS J1206+4332 and PG 1115+080 have the largest relative uncertainty on σ_v among the TDCOSMO sample. Therefore, the zero points on the x-axis of Fig. 1 for these two objects are the most uncertain.

4.2. Sensitivity to kinematics of the different mass model families

We repeat the experiment for power-law models and composite model separately to check the sensitivity to kinematics data of each family of mass models. We do not use B1608+656 when computing the sensitivity to kinematics of the composite models since this system has a pixelated potential correction performed on the power-law model, but no composite model. Bottom panels of Fig. 1 show the result of this test.

We obtain ⟨ξ_composite⟩ = 0.06 ± 0.02 and ⟨ξ_PL⟩ = 0.07 ± 0.01. The value of the joint inference is similar for both the composite and the power-law cases but each lens behaves differently. While WFI 2033−4723 becomes more sensitive to the kinematics when modeled only with a composite model, SDSS J1206+4332 has its sensitivity almost halved. We can explain this behavior as due to the relative precision of the two families of models, which is different from one lens to the other. The time-delay distance of SDSS J1206+4332 is better constrained by composite models ($D_{\mathrm{\Delta }t}={5690}_{-356}^{+449}$ Mpc at 7.1 % precision) than with power-law models ($D_{\mathrm{\Delta }t}={5873}_{-659}^{+659}$ Mpc at 11.2 % precision). The relative weight of the D_Δt compared to the D_d in the final value of H₀ is therefore more important in the composite model case.

WFI 2033−4723 experiences the opposite behavior; it has tighter constrains with power-law models ($D_{\mathrm{\Delta }t}={4701}_{-204}^{+242}$ Mpc at 4.74 % precision) than with composite models ($D_{\mathrm{\Delta }t}={4909}_{-319}^{+485}$ Mpc at 8.2% precision). WFI 2033−4723 is therefore more sensitive to the kinematics in the composite model case.

In summary, there is no evidence that one family of mass models is significantly more sensitive to the kinematics than the other. For individual lenses, we observe differences but they can be explained by the relative precision that each of the models can achieve on the D_Δt measurement with respect to their D_d measurement, based on the relative weight of the lensing and kinematic constraints and on the redshift of deflector and source that determine how the D_Δt − D_d constraint maps into H₀.

5. Search for correlations between H₀ and physically independent observables

The inference of H₀ relies on many independent ingredients and observables, such as the velocity dispersion of the deflector and the relative density of galaxies in the line of sight up to the background quasar. Those quantities do not have any physical reason to be correlated with H₀. Thus, any evidence of a correlation between these observables and the inferred value of H₀ across the TDCOSMO sample, beyond the expected error covariance, would be an indication of underlying systematic errors. In this section, we carry out a number of empirical tests, correlating H₀ with observables and properties of the instrumental setup, and find no evidence for any statistically significant dependency.

5.1. Dependency on the characteristic scale of the lens system and spectroscopic aperture.

Figure 2 shows the inferred Hubble constant for each of the seven TDCOSMO lenses for several combinations of characteristic scales of the lens systems and the aperture used for spectroscopic follow-up. In the left panel, we use the ratio between the Einstein and the effective radii to investigate any departure from the assumed description of the radial mass density profile. The ratio between the effective radius and the Einstein radius is used as a diagnostic of the relative spatial distribution of luminous and total matter. If the TDCOSMO models were insufficiently flexible, one may expect a trend in this ratio because the sum of the dark and luminous component would produce different shape of the total mass profile and a lack of flexibility in the mass model would not be able to reproduce the correct underlying distribution. In the middle panel are shown the ratios between Einstein radius and the spectroscopic aperture, which compare the spatial scales at which the lensing and kinematic information is obtained. Finally, the right panel of Fig. 2, shows the ratio between the effective radius and the radius of the spectroscopic aperture, which could potentially be affected if the stellar kinematics were incorrectly modeled. One expects trends in all the above quantities if, for example, the assumptions about orbital anisotropy were systematically wrong.

Fig. 2. Effective radius θeff, Einstein radius θE and radius of the spectroscopic aperture θaperture of the TDCOSMO lenses. We show the ratios of these three quantities and the corresponding H0 value inferred for each system. We do not observe significant correlations between the characteristic sizes of the lens, the spectroscopic aperture and H0. The horizontal lines indicate the latest H0LiCOW 2019 (dotted orange, Wong et al. 2020) and Planck (dashed blue, Planck Collaboration VI 2020) results along with the 1σ uncertainties.

The Spearman’s rank correlation coefficient between H₀ and θ_E/θ_eff, θ_E/θ_aperture and θ_eff/θ_aperture are respectively −0.11, 0.42 and 0.67. The probability that an uncorrelated data set produces such correlation coefficients (i.e., the p-value) is 0.82, 0.33 and 0.10. Therefore we conclude that in all three cases, there is no statistically significant correlation, even though the dynamical range on the x-axis is a factor of 3–6. While the absence of correlations does not prove that all systematic errors are below the statistical uncertainties, this is an important sanity check for our current models and for future work as the statistical precision improves with growing sample size.

In addition, observational and modeling effects such as the choice of stellar template, the choice of anisotropy model, or the PSF modeling could potentially bias the measured velocity dispersion of the main deflector and thus H₀. The net effect of all these possible sources of systematic errors is difficult to quantify exactly but they typically scale with the effective radius of the lens θ_eff or the aperture radius of the spectroscopic observation θ_aperture. The absence of any trend in Fig. 2 is reassuring in this regard. Moreover, as we showed in Sect. 4, even ∼5% systematic bias on the measured velocity dispersion, or equivalently on the modeled quantities due to incorrect anisotropy assumptions, will only produce an average 0.35% bias on H₀. Furthermore, as shown above, the direction and amplitude of the error would be different for each lens and therefore this systematic uncertainty would also show as a source of scatter or trend across the sample, which are not observed.

5.2. Dependency on intrinsic parameters of the deflector traced by the velocity dispersion

An additional potential concern is whether systematic differences between our assumptions and the internal structure of early-type galaxies could give rise to measurable biases. For example, the so-called “tilt” of the fundamental mass plane is believed to arise primarily from the increase in dark-to-stellar matter ratio, a systematic change in stellar initial mass function with galaxy stellar mass, and possibly a small subdominant contribution from systematic variations in stellar orbits anisotropy (Auger et al. 2010; Cappellari 2016). The stellar initial mass function is not a concern in the TDCOSMO analysis, since the stellar mass to light in the composite models is a free parameter. However, in principle the other two sources of “tilt” could introduce a potential systematic effect in TDCOSMO analysis, where each system is analyzed independently and with the same priors, rather than with priors that depend on the stellar mass.

In Fig. 3 we show the inferred H₀ as a function of stellar velocity dispersion, a redshift independent proxy of position along the fundamental plane. In this case, we found a Spearman’s rank correlation coefficient of 0.07 with a p-value of 0.88. Hence, we conclude that there is no statistically significant trend in these data, indicating that any residual velocity dispersion dependent bias is smaller than the measurement uncertainties, and thus not significant at this stage. As for the plots shown in the previous (and next) section, this sanity test should be repeated as the sample size and individual measurement precision increase.

Fig. 3. Hubble constant as a function of the measured velocity dispersion of the main lens. The horizontal lines indicate the latest H0LiCOW 2019 (dotted orange, Wong et al. 2020) and Planck (dashed blue, Planck Collaboration VI 2020) results along with the 1σ uncertainties.

5.3. Dependency on the external convergence and lens redshift

In the previous sections, the focus is on how the lens velocity dispersion influences H₀ measurements. But there is also an external contribution of all objects along the line of sight to the main lensing potential. This external convergence, κ_ext, is estimated in all TDCOSMO systems from galaxy counts, in combination with spectroscopy for obtaining redshifts for galaxies and quantifying coherent structures (e.g., groups and clusters). Tihhonova et al. (2018) showed that this measurement is compatible with the constraints obtained on κ_ext with weak lensing. κ_ext is directly related to the time-delay distance D_Δt, as shown in Eq. (7). Similarly, the effect of the external convergence on the inferred H₀ can be written as :

$$\begin{array}{c}\hfill H_0^{\mathrm{uncorr}}=\frac{H_0^{\mathrm{corr}}}{(1-{\kappa }_{\mathrm{ext}})},\end{array}$$(16)

where $H_0^{\text{uncorr}}$ ($H_0^{\text{corr}}$) is the value of H₀ before (after) correction from κ_ext. The effect of this external MST can be mitigated by directly inferring κ_ext. To test the presence of residual external Mass-Sheet Degeneracy (MSD) not entirely removed by the measurement of κ_ext, we investigate the presence of correlation between the estimated κ_ext and the inferred H₀ value for the seven lenses of the TDCOSMO sample. The top panel of Fig. 4 shows the relation between the H₀ measurements before correction for the mass along the line of sight $H_0^{\text{uncorr}}$, and the estimated convergence. A trend is visible between these two quantities indicating that the measurement is indeed sensitive to the lens environment. If no correction is applied, the lenses located in over-dense regions (positive κ_ext) tend to have a higher $H_0^{\text{uncorr}}$ than lenses in under-dense regions (negative κ_ext). We fit a linear model to the uncorrected data, and measure a slope of a^uncorr = 88.9 ± 29.1 km s⁻¹ Mpc⁻¹, well compatible with the expected slope of $a^{\text{uncorr}}=H_0^{\text{corr}}=73.7\text{km}{\text{s}}^{-1}{\text{Mpc}}^{-1}$. Both the uncorrected and corrected data are well fitted by our linear model, with a reduced χ² of 0.61 and 0.95 respectively.

Fig. 4. Measured Hubble constant, before (upper panel) and after (lower panel) correction for the mass along the line of sight as a function of the estimated external convergence. $H_0^{\text{uncorr}}$ and $H_0^{\text{corr}}$ are related according to Eq. (16). The dashed black lines show the best linear fit, and the shaded gray envelopes correspond to the 1σ uncertainties. The dotted blue lines represent the relation expected from the theory between $H_0^{\text{uncorr}}$, $H_0^{\text{corr}}$ and κext.

As shown on the bottom panel of Fig. 4, this trend disappears when correcting for the external convergence and there is no evidence for residual correlation between $H_0^{\text{corr}}$ and κ_ext. In fact, the best-fit slope coefficient in this case is a^corr = −5.1 ± 23.7 km s⁻¹ Mpc⁻¹, consistent with no correlation. This is an indication that the external convergence correction makes the trend disappear, which is what would be expected if our correction were accurately accounting for κ_ext. The present data set shows no evidence of residual systematic bias involving the LOS mass density.

As first mentioned by Wong et al. (2020), the H0LiCOW collaboration reported the presence of a possible trend between the lens redshift and the inferred $H_0^{\text{corr}}$ value at low statistical significance level (∼1.9σ). When adding DES J0408−5354 to the six H0LiCOW lenses, the significance of the trend is slightly reduced to ∼1.7σ. We note that, having tested multiple correlations, it might be expected to find one at marginal significance, as a result of the “look elsewhere effect”. This trend is still present before correction for the external convergence as shown on Fig. 5. The data are well-fitted by a linear model both before and after the LOS correction with a reduced χ² of 1.52 and 0.24. The significance level of this correlation before LOS correction is still on the order of ∼2σ. Hence, there is no direct indication that the trend is due to unaccounted systematics in κ_ext.

Fig. 5. H0 constraints for the TDCOSMO lenses as a function of lens redshift before (top) and after (bottom) correction for the external convergence. The best linear fits and their 1σ envelopes are shown in shaded gray. The tentative (1.7σ significance) trend is not introduced by the LOS contribution as it is still visible before correcting for the external convergence.

6. Impact of the choice of families of mass model

In this section we quantify how much the inference on H₀ depends on the choice of the mass density profile adopted for the lens modeling. We first use the six systems for which both power-law and composite mass models have been performed and compare the results. We show that even though the two model families have sufficient flexibility to produce a broad range of profile shapes, in practice when applied to real elliptical massive galaxies, they form mass density profiles close to a simple power law. As we see below, this is likely due to the “bulge-halo” conspiracy (Treu & Koopmans 2004; Dutton & Treu 2014).

Then, we carry out end-to-end simulations in order to quantify the flexibility of our models and how the data actually allows us to constrain them. Meeting this goal requires the simulated properties of lenses to be close enough to those of real galaxies. About 90% of galaxy-scale lenses are early-type galaxies (Auger et al. 2009), which satisfy very tight correlations between their observable properties (Auger et al. 2010). This indicates a high degree of regularity in the relative distribution of dark and luminous matter, often referred as the “bulge-halo conspiracy”. This bulge-halo conspiracy results in the total mass density profile of lenses being very close to a singular isothermal ellipsoid (e.g., Koopmans et al. 2006, 2009; van de Ven et al. 2009; Cappellari 2016), even out to large radii (Gavazzi et al. 2007; Lagattuta et al. 2010).

Importantly, the simulations we use all consider spatially extended lensing information, spanning a large range in radial extension. This radial extent must provide sufficient leverage to inform us about any possible departures from a simple power law within the actual range of observables. A goodness-of-fit criterion is then used to verify that the model adopted is indeed a good description of the data. Models that are exclusively based on the positions of two or four multiple quasar images, rather than the full surface brightness distribution of its spatially-extended host galaxy, cannot provide an accurate account of the uncertainties from surface brightness modeling. Therefore, models based on two or four image positions cannot satisfy the above goodness-of-fit requirement, even if they include time delays and stellar-velocity-dispersion measurements. In the following, we describe our set of simulated lenses in Sect. 6.2, present the results in Sect. 6.3 and discuss our findings in Sect. 6.4.

6.1. H₀ inference per model family

The TDCOSMO collaboration uses both composite and power-law models in their analysis, except for B1608+656 (see Sect. 2.3). Apart from this exception, the published estimates of H₀ correspond to the marginalization over the two model families as a way to account for modeling uncertainties (Wong et al. 2020; Shajib et al. 2020). The sample size of real lenses is now sufficiently large to infer H₀ by model family and to test whether this choice makes a difference at the 2% precision level of the statistical uncertainty. This is illustrated in Fig. 6, where the priors on the cosmological parameters are the same as adopted by Wong et al. (2020): H₀ ∈ [0, 150] km s⁻¹ Mpc⁻¹, Ω_m ∈ [0.05, 0.5] and Ω_m = 1 − Ω_Λ.

Fig. 6. Marginalized H0 posteriors for power-law (left panel) and composite models (right panel). The cosmological inferences are for a flat ΛCDM cosmology with uniform priors. The posterior probability distributions for each individual system are shown with shaded color curves and the combined constraint from the six systems corresponds to the solid black curve. The legend indicates the median, 16th and 84th percentiles of the H0 distributions.

The H₀ values vary with the model family for individual objects, and this testifies to the flexibility of the families of models. However, the choice of model family changes the combined value by much less than the estimated statistical uncertainty. Quantifying these statements, the combined value from the six lenses is $H_0=74.2_{-1.6}^{+1.6}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$ when we use exclusively power-law models and $H_0=74.0_{-1.8}^{+1.7}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$ when we use only composite model. This corresponds only to a 0.2% difference. Individual objects can have larger differences between power-law and composite models than the combined estimate, but the two posterior probability distributions always remain compatible. The largest differences are found for PG 1115+080 (5%) and SDSS J1206+4332 (4%), which still have the two distributions compatible at the ∼0.6σ level.

Last but not least, there is no indication in the current sample of six lenses that one given family of models systematically gives a lower or higher H₀ value. For example, WFI 2033−4723 has a higher H₀ value when modeled with a power law rather than a composite, while the opposite behavior is found for SDSS J1206+4332; and other such examples can be easily found in Fig. 6.

In conclusion, even though our two families of models are flexible enough to produce a broad range of H₀ values, in practice they do not. In the following, we investigate with simulated lens systems the reasons why composite and power-law models provide comparable estimates of H₀ in spite of allowing for flexibility. We also investigate under which circumstances gravitational lenses can be modeled with both composite and power-law models and still yield the same H₀.

6.2. Simulations

We generate six mock lens systems chosen to illustrate the range of possible outcomes, labeled by IDs #1 through #6. We describe the process of the simulations in this section. In addition to the power-law and composite models typically used by TDCOSMO we also include cored power laws to explore the effects of adding extra flexibility to the models.

The simulated HST images are produced using the pipeline described by Ding et al. (2017a, 2018). The image frame size is chosen to be 99 × 99 pixels, with a pixel scale of $0\stackrel{″}{.}08$ to mimic the realistic HST WFC3/F160W drizzled resolution. Mass profile parameters are chosen such that the Einstein radius is roughly at the scale of 1″ as typical for galaxy-scale lenses. The noise in each pixel is composed of the Gaussian background noise and the Poisson noise. For Gaussian background noise, we assume an rms of 0.003, which is directly measured from empty regions in the real data; the Poisson noise is added, based on a total exposure time of 2400 s. For computational speed, the PSF is assumed as a Gaussian kernel with $FWHM=0\stackrel{″}{.}25$.

Three mass models, including power-law (ID #1, #2), cored power-law (ID #3, #4), and composite (ID #5, #6) mass density profiles, are adopted to generate the six mock systems. All of the systems are elliptical in projection in order to allow for quad-like configurations by construction. For each family of mass distribution, we generate two mock lensed systems, one with the source lying close to a fold of a caustic (“fold” configuration) and one with the source lying close to the lens-optical axis (“cross-like” configuration). The “cross” represents a worst case scenario because the radial ranges and the differences in the time delays are limited by symmetry. The simulated lens systems are shown in Fig. 7.

Fig. 7. Sample of simulated lenses: three pairs are generated from power-law, cored power-law, and composite lens models. The color scale is logarithmic and is the same for all images. Identifiers associated to each lens are also indicated. Refer to Sect. 6.2 for a description of these simulations. Model #6, although composite, is chosen so that the total mass profile resembles a power law in the region of the Einstein radius.

For the composite model, the total mass consists of a baryonic elliptical Hernquist profile (Hernquist 1990), and a dark matter elliptical NFW profile (see Eq. (14) and Navarro et al. 1997). The baryonic part is linked to the lens surface brightness through a constant mass-to-light ratio. While we use the same axis ratios for the baryonic and dark matter components, we allow for slight offsets in their position angles; the total projected mass profile is therefore not elliptical. We note that the system (ID #6) is chosen to describe a scenario similar to realistic galaxies, in which luminous and dark matter conspire to produce a total mass model very close to a power-law profile. This is consistent with the findings of the H0LiCOW, SHARP, and STRIDES collaborations so far (Suyu et al. 2014; Wong et al. 2017; Birrer et al. 2019; Chen et al. 2019; Rusu et al. 2020; Shajib et al. 2020). Other cored power-law and composite systems (ID #3 – #5) are designed on purpose to depart significantly from a single power law in order to test the effect on H₀ and investigate whether the information contained in the data can capture this discrepancy. For all the lenses, the deflector surface brightness is simulated as an elliptical Hernquist profile. The ellipticity of the simulated lens galaxy corresponds to an axis ratios of q ∼ 0.9 ± 0.01. We use an elliptical Sérsic profile (Sérsic 1963) to simulate the extended part of the source light, which is sufficient for our purpose. Lensed quasar images are modeled as point spread functions centered on the images of the host galaxy.

The simulated time delays are calculated within a fiducial flat ΛCDM cosmology with Ω_m = 0.27, and Ω_Λ = 0.73, and Hubble constant $H_0^{\text{fiducial}}=70.7\text{km}{\text{s}}^{-1}{\text{Mpc}}^{-1}$, which was chosen randomly. For the time-delay uncertainties, we assume an unbiased random error with rms level set as the largest value between Δt × 1% and 0.25 days. The uncertainties on the time delays are chosen to be smaller than current uncertainties of real data in order to focus mainly on the modeling uncertainties.

Since the tests in this section focus on the mass reconstruction of the main deflector, we do not include in the simulations the effects of the galaxies along the line of sight, which are treated separately in real data. Likewise, we simulate and model the velocity dispersion using spherical Jeans equations following Suyu et al. (2010) and Birrer et al. (2019), and assume an anisotropy radius equal to the lens half-light radius. This is a simplification of the stellar kinematics treatment with respect to the analysis of real systems where TDCOSMO marginalizes over the unknown anisotropy. In this exercise where we aim to illustrate the constraining power of the images while saving computing time, we do not use the LOS velocity dispersion as a direct constraint in the modeling but rather only calculate the modeled values to make the comparison with measured values. The relevant key properties of the six simulated lenses are summarized in Table A.1.

6.3. Results

The six mock lenses are modeled using the public strong lensing modeling package LENSTRONOMY⁴ (Birrer et al. 2015; Birrer & Amara 2018), which was used for the latest analysis of the real systems SDSS J1206+4332 and DES J0408−5354 (Birrer et al. 2019; Shajib et al. 2020). The exact/known input PSF is used as the effect of PSF imperfections is not investigated in this work. The light profile of the lens and of the source are modeled as Hernquist and Sérsic profiles respectively. We fit three types of analytical elliptical mass profiles to the simulated data, namely a power-law, a cored power-law and a composite profile. Specifically for the composite model, we emphasize that no strong prior is applied on the scale radius of the dark matter component. Instead, we use a noninformative uniform prior r_s ∼ 𝒰(5″, 40″), so that the dark matter component effectively has two degrees of freedom in the radial direction. The 99 × 99 pixels contained in the images and three independent time delays are used for the fit. We, however, mask a central region corresponding to three pixels (i.e., $0\stackrel{″}{.}24$) since we do not want to form any central image which could lead to extra constraints on the lens model (see also Tagore et al. 2018; Mukherjee et al. 2018, 2019). The resulting fitted models are used to infer only H₀ (Ω_m is kept fixed to 0.27) from the time-delay distance alone. The lens velocity dispersion is computed only for comparison but is not included in the H₀ inference, to highlight the information content of the images.

We use the Bayesian Information Criterion (BIC) to evaluate the quality of the fit. The BIC is defined by

$$\begin{array}{c}\hfill \mathrm{BIC}=k\times \ln (n)-2\times \ln (\widehat{L}),\end{array}$$(17)

where k is the number of free parameters, $\widehat{L}$ is the maximum likelihood of the model and n is the number of data points. The likelihood used for the fit uses only the imaging and time-delay information so that n corresponds to the number of nonmasked pixels in the image plus the three time delays. Our models have 25 free parameters for the power-laws, 26 for the cored power-laws and 29 for the composite models.

The recovered H₀ value, integrated LOS velocity dispersion within a square aperture of side 1″ and the BIC values are given in Table 2. The corresponding image residuals of the lens modeling are shown in Table 3. As expected, we recover the correct H₀ value within the 1-σ errors of the posterior distribution when fitting the same mass model family as used in the simulation. This case corresponds to the diagonal of Tables 2 and 3.

Interestingly, the core size of the cored power-law profile is well constrained by the data. Indeed, when a cored power-law profile is fitted to data generated with power law with no core, the core size is well constrained and shrinks to zero. If there is a core in the simulation (e.g mock lenses #3 and #4), the core-size is recovered within 2.2% accuracy and within < 3.0% precision with a cored power-law model. This indicates that the lensing data are sensitive to the presence of a sizeable core in galaxies. The sensitivity stems from the robust constraint on the mass enclosed within the Einstein radius that indirectly depends on the core size.

We deliberately choose not to present the results of the composite models fitted to power-law and cored power-law simulations. This is because, by construction, the lens light profile of these simulations does not necessarily correspond to their mass profile. In the power-law and cored power-law profiles, the lens light profile bears no relation to the mass distribution, and is only used as a tracer of the stars when computing the stellar velocity dispersion. As a result, we cannot have a meaningful comparison between power-law and composite models if we assume that the baryonic component of the composite model is traced by the arbitrary lens light in the power-law model. This limitation is inherent to these simulations and we do not expect that this happens for real galaxies, because it is unlikely that the stellar light is not tracing at all the baryonic mass component.

The tests performed on composite simulated lenses #5 & #6 show that the ability of a power law or a cored power law to recover the correct H₀ depends on the characteristics of the composite lenses. In both cases, the power-law models give much poorer fits to the data than the true composite models (ΔBIC = 434 for #5 and ΔBIC = 4455 for #6). Adding one more degree of freedom by using a cored power law instead of a power law improves the fit but it is still significantly poorer than the composite models (ΔBIC = 95 for #5 and ΔBIC = 1049 for #6 in the case of a cored power law). We note that the image residuals in lens #6 are worse than that in lens #5, since #6 is in a fold configuration with higher lensing magnifications and thus produces correspondingly higher amounts of image residuals. The recovered H₀ is compatible with the true value for the lens #6, but in lens #5 it is biased toward lower H₀ by 9.4%. In short, the different behavior arises because of intrinsic differences in the composite mass density profile. While mock #5 is chosen to be different from a power law, mock #6 is chosen to be similar to a power law. When the truth is a composite similar to a power law, the inferred H₀ is the same. When it is not, the two models lead to different inferences. As discussed in Sect. 6.1 the real universe is similar to #6 and dissimilar to #5.

As an additional test, we model the simulated data using only the four lensed image positions, the lensing galaxy position and the time delays to investigate the effect of neglecting the other information. We find that, as mentioned in Sect. 3, this is not sufficient to constrain all the lens model parameters. A reduced χ²< 1 can be obtained for all the mocks using a power law model, except for mock #6 for which the best reduced χ² is ∼1.9. Even when the true mass distribution is a power law (e.g., mock#1 and #2), the maximum likelihood models are associated with power-law indices substantially different from the input one, yielding a bias on H₀, that can reach 90% (see Appendix B for details). This is well understood as the multipole components of the lens potential can compensate for large changes in the monopole structure which are only poorly constrained by the few image positions. This test highlights the necessity of using the full information provided by the high resolution images to better constrain the lens potential. In particular, the multiple images of the lensed host galaxy are critical to pin down the uncertainty on the average mass density at the image positions (Kochanek et al. 2001).

6.4. Discussion

In this section we discuss the results of the simulations with the goal of providing an intuitive physical understanding of the quantities that are relevant for time-delay cosmology and how they are constrained by the data. As noted by Kochanek (2002), the time delay is mainly determined by the mean convergence ⟨κ⟩ in an annulus between the multiple images. Figure 8 shows the radial convergence profiles of the models averaged over the azimuth angle. The shaded gray contour corresponds to the separation between the multiple images. The quality of the fit in this region determines the accuracy on H₀. The Einstein radius is typically very well constrained by any lens model, so the only way to modify the mean ⟨κ⟩ at the positions of the multiple images is to change the slope of the convergence profile while keeping constant the integrated mass within the Einstein radius. This is a well-known problem in time-delay cosmography called the profile slope degeneracy (Witt et al. 2000; Wucknitz 2002; Suyu 2012).

Fig. 8. Azimuthally averaged radial convergence profiles, for the different lens models applied to fit the sample of six mock lenses (Fig. 7). Upper part of each panel: true profiles are shown in dotted lines; power-law, cored power-law and composite models are shown in blue continuous, green dashed and red dot-dashed lines, respectively. The spectroscopic (square) aperture used for computing velocity dispersions is indicated as a vertical dotted line, and the true Einstein radius location is indicated as a vertical dashed line. The gray area encloses lensed quasar image positions. For each model, the inferred H0 values are indicated (in km s−1 Mpc−1), and must be compared to the input value $H_0^{\text{fiducial}}=70.7\text{km}{\text{s}}^{-1}{\text{Mpc}}^{-1}$. Lower part of each panel: relative error computed as (truth − model)/truth. The pixel size in the simulated images is $0\stackrel{″}{.}08$.

As argued by Sonnenfeld (2018), assuming a too rigid model such as a power law model, can lead to a bias up to ∼10% if the true underlying profile contains a change of slope within the Einstein radius. Sonnenfeld (2018) concluded that at least three degrees of freedom are required in the lens model to recover an un-biased result if no kinematics information is used. With the addition of kinematics, uncertainty can be reduced to 1% (in accuracy) even within the simplified constraints considered in that study. Based on a sample of simulated galaxies from the Illustris simulation, Xu et al. (2016) studied how physically motivated numerical density profiles are transformed into an approximate power-law (in the region where lensed images are formed) by means of a mass-sheet like transformation. They reported that a large range of transformation was allowed, which would translate into a large scatter and possible bias on the inferred H₀. They concluded that the amplitude of the bias depends on the (logarithmic) curvature of the mass density profile of the simulated galaxies. This behavior was previously illustrated in Schneider & Sluse (2014) and more recently in Gomer & Williams (2019).

We recover the findings of Sonnenfeld (2018), Xu et al. (2016), Schneider & Sluse (2014) and Gomer & Williams (2019) with our simulated lens #5, where the combination of the Hernquist and NFW profile is designed to produce an inflection point in the radial profile of the convergence within the Einstein radius. For this system, the composite and power-law models are discrepant, thus providing an indication that the power-law model is indeed too rigid. This rigidity results in a significant difference in goodness of fit (ΔBIC = 434), as well as on the inferred H₀.

For the lens system #6, the radial convergence profile does not have inflection points and therefore it is impossible to change the slope of the profile while keeping the Einstein radius identical. In this case, the recovered value of H₀ is compatible with the true value for both the composite and power-law model. The fact that the two families of models are providing compatible H₀ indicates that the convergence profile is well-recovered in the annulus around the Einstein radius.

The TDCOSMO collaboration has systematically tested both model families in their analysis after the first and only nonblind published system B1608+656. The tight agreement between the composite and power-law models in the TDCOSMO analyses supports the hypothesis that, as a result of the bulge-halo conspiracy, the kind of real galaxies that act as strong lenses are similar to our #6 mock. The mass density profile is well approximated by a power law. In this case, the stellar component and the extended NFW halo combine to form a profile very close to a power law near the quasar images. If this had not been the case, we have shown in this work that the composite and power law would not have produced the same mean convergence ⟨κ⟩ at the image position, and thus would have yielded very different H₀ values. If TDCOSMO had found this discrepancy, it would have been accounted for in the error budget of each individual lens since they marginalize over model families. In contrast, both classes of models produce fits with comparable goodness-of-fit and H₀ in the real world, resulting in high precision, including modeling errors in the error budget. Our analysis in Sect. 6.1 shows that in practice the two models agree even at the current sample precision of < 2.4%.

7. Conclusion

As the statistical precision of time-delay cosmology improves with the analysis and publication of multiple gravitational lens systems by the H0LiCOW, SHARP, and STRIDES collaborations, a parallel effort must be undertaken to ensure that systematic uncertainties remain subdominant. In this first paper of the TDCOSMO collaboration (i.e., COSMOGRAIL, H0LiCOW, SHARP, STRIDES members), we investigate and quantify a number of potential systematic uncertainties that could affect the analysis. Before we summarize the main results of this work, it is important to highlight a few general points that are relevant to the estimation of systematic errors in time-delay cosmography:

1. The TDCOSMO analyses are carried out blindly to cosmological parameters, with the exception of the first system B1608+656 (Suyu et al. 2010) in order to avoid implicit/explicit experimenter bias.

2. The TDCOSMO estimates of H₀ are obtained independently for each lens, and they are found to be statistically consistent with each other (Wong et al. 2020). The statistical consistency demonstrates that uncorrelated systematic errors are negligible with respect to statistical errors. So any investigation of systematic errors must focus on correlated errors that would affect many systems in the same way.

3. Toy models based on simplified assumptions and constraints cannot offer any quantitative estimates of systematic errors given the current state-of-the-art data-sets and lens models. The only way to estimate quantitative errors is to carry out an analysis that is very similar to the one performed on real data, using the full extent of the available information, including the high-resolution images, multiple time delays (if available), and stellar kinematics. For example, the dependency on the inferred distances on stellar velocity dispersion is nontrivial, it varies from lens to lens, depending on the precision of the various constraints, the lensing configuration, the source and deflector redshift, and the spectroscopic aperture used for the kinematic measurement. On average over the current TDCOSMO sample, uncertainty in velocity dispersion δσ_v/σ_v translates into approximately δH₀/H₀ ∼ 0.07 × δσ_v/σ_v.

Keeping these general considerations in mind, the main results of this work are as follows:

• No evidence is found for any correlation between the measured value of H₀ and observables related to the internal structure of the lens galaxies (e.g., velocity dispersion, effective radius), or to the size of the spectroscopic aperture. If our assumptions about the kinematic field of the lens galaxies had been significantly wrong, then we would have expected to detect trends in these parameters, since our deflectors and spectroscopic observations span a significant range of configurations. Of course absence of evidence is not evidence of absence and more work remains to be done in this area, even though the weak dependency of the inferred H₀ on kinematic data implies that systematic uncertainties in this area will have a subdominant impact on H₀.

• No evidence is found for any correlation between the measured value of H₀ and the external convergence estimated from galaxy number counts and numerical simulations. In contrast, if no external convergence is applied, H₀ is found to depend on the overdensity of galaxies in the field, a clearly unphysical result.

• Tests based on mock lens systems that have simulated data comparable in quality to real lens systems show that the current approach of considering different mass profiles has sufficient flexibility in the mass model to infer a wide range of H₀ values, should the data require it.

• Mock lens galaxies composed of baryons and dark matter whose total mass distribution is not well approximated by a power law produce discrepant H₀ inferences and significant differences in image residuals when comparing power-law and composite mass models. In contrast, mock lens galaxies whose baryonic and dark matter components conspire to form a power law lead to comparable H₀ inferences between power-law and composite mass models.

• The comparison of power-law and composite mass models allows us to quantify deviations in H₀ due to our mass model assumptions. By using these two families of models and marginalising over them, the resulting H₀ accounts for modeling uncertainties. Future measurements of spatially resolved kinematics of the lens system would provide highly constraining measurements of the lens mass distribution that potentially allow us to distinguish/rank mass models, removing the need to marginalize over degenerate lensing mass models.

• The similarity of H₀ constraints from power-law and composite models of TDCOSMO lenses shows that the total mass profiles of galaxies are close to power laws due to the bulge-halo conspiracy. For the six lenses that have been analyzed with both power-law and composite models we find $H_0=74.2_{-1.6}^{+1.6}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$ and $74.0_{-1.8}^{+1.7}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$ respectively. The difference between the two model families is much smaller than the inferred statistical errors. The similar H₀ from the different families of models thus made the current H₀ measurement with ∼2% uncertainty from TDCOSMO achievable.

Based on a number of tests carried out in this paper, we find no evidence that the error budget reported by the H0LiCOW/SHARP/STRIDES (TDCOSMO) collaborations is significantly underestimated. We emphasize that our tests reproduce very closely the TDCOSMO inference procedure, in contrast to previous work in the literature that does not have the fidelity to investigate this issue.

While investigating potential sources of systematic uncertainties was an important first step, meeting the goal of 1% precision and accuracy with time-delay cosmography (e.g., Suyu 2012; Treu & Marshall 2016), requires additional and continued efforts over the coming years. Aside from expanding sample sizes and improving statistical precision per system, some of the clear steps along the way are: (i) exploring broader model families and the impact of departures from elliptical symmetry and including spatially variable mass-to-light ratio, (ii) explore in more detail the bulge-halo conspiracy based on high resolution data for local early-type galaxies, (iii) explore the effect of allowing for gradients in stellar mass-to-light ratios (e.g., Sonnenfeld et al. 2018) in composite models; (iv) carrying out a full Bayesian hierarchical analysis of existing samples of lenses in order to constrain parameters that cannot be inferred on single lens but require an inference at the population level, (v) accounting for measurement and modeling covariance, and (vi) performing realistic data challenges such as the one proposed by Ding et al. (2018), with increasing level of realism and complexity as data also improve. These steps are nontrivial from a modeling point of view, considering that the analysis of any single system currently requires a year of expert investigator time and on the order of a million CPU hours (e.g., Shajib et al. 2020). Substantial advances in automation and speed are required in order to carry out those next steps, but given their importance for the determination of H₀, they are worth undertaking.

Acknowledgments

This program is supported by the Swiss National Science Foundation (SNSF) and by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (COSMICLENS: grant agreement No 787886). Additional funding is provided by the Packard Foundation through a Packard Research Fellowship to T.T., by the National Science Foundation through grant AST-1906976 to T.T. and by NASA through HST grants HST-GO-10158, HST-GO-12889, HST-GO-14254, HST-GO-15320, HST-GO-15652. S.H.S. thanks the Max Planck Society for support through the Max Planck Research Group. G.C.-F.C. acknowledges support from the Ministry of Education in Taiwan via Government Scholarship to Study Abroad (GSSA). C.D.F. and G.C.-F.C. acknowledge support for this work from the National Science Foundation under Grant Nos. AST-1715611 and AST-1907396. This work was supported by World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan. L.V.E.K. has been supported through an NWO-VICI grant (project number 639.043.308). 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) and emcee, a Python implementation of an affine invariant MCMC ensemble sampler Foreman-Mackey et al. (2013).

References

Appendix A: Properties of simulated lenses

As a complement to Sect. 6, we show in Table A.1 a subset of important properties of the simulated lenses. Characteristic radii are indicated: half-light radius, effective Einstein radius, core radius in the case of cored power-law profiles, and scale radius of the dark matter profile for composite models. The ratio of the lens half-light radius and Einstein radius is also computed. Additionally, the input logarithmic slope of the convergence profile, the lens mass ellipticity, true time delays and LOS velocity dispersion are indicated. The spectroscopic aperture used for simulating and modeling kinematics is a square with side 1″. For composite models, we provide the dark matter fraction within the Einstein radius. Lastly, to ease comparison with previous studies, we add the measure of the curvature of the total mass ξ, as defined in Xu et al. (2016). Given this definition, a concave-upward (convex-downward) radial convergence profile has curvature greater than one (lower than one), and a perfect power-law have curvature equal to one. Xu et al. (2016) conclude that galaxies close to isothermal and those having a curvature parameter close to one provide the smallest bias on H₀. We recover these findings only partially with these simulated lens systems. We find that the curvature criterion is the most important criterion to ensure a low bias on H₀ even if the slope differs significantly from isothermal, as illustrated with our simulated galaxy #6.

Appendix B: Models only based on lensed quasar positions and time delays

We model the simulated dataset as generated in Sect. 6.2, using only the lensed quasar positions, the lensing galaxy position and the relative time delays between the lensed images. We assume an uncertainty $\sigma =0\stackrel{″}{.}004$ on the point-source positions, $\sigma =0\stackrel{″}{.}01$ on the lensing galaxy centroid, and the same uncertainty on the time-delay as in Sect. 6.3. Similarly to extended source modeling performed in Sect. 6.3, we employ the lens modeling package LENSTRONOMY. We adopt both the Singular Isothermal Ellipsoid (SIE) model (i.e., fix slope value γ = 2.0) and the Power-law model in this test. An independent modeling has been carried out with lensmodel (Keeton 2001, 2011). We obtained similar inference with both packages and therefore only report hereafter results obtained with LENSTRONOMY.

We use the true parameters as the input values to start performing the modeling. A careful choice of the likelihood and sampling options has to be carried out to ensure that image position constraints arise from the same source. In practice, we sample the source plane and evaluate the positional likelihood in the image plane, but adding a source plane likelihood term to ensure that each image arises from the same source within $\sigma =0\stackrel{″}{.}001$. A notebook implementing our fitting strategy is available online⁵.

In Fig. B.1, we show the corner plots of the inference based on the mock system #1. When using a SIE model where the mass slope value is fixed to the truth, we could obtain an unbiased H₀ with uncertainty at the ∼10% level. However, when the slope is a free parameter, the inference broadens significantly as the data are not sufficient to constrain that parameter. In particular, the uncertainty on H₀ increases by a factor of ∼3 and the maximum likelihood deviates from the truth by up to 100%. This contrasts with the same model constrained by the point-source and extended images from the source. Those features constrain accurately the position of the centroid of the lens potential and its ellipticity, breaking degeneracies between those quantities and H₀.

Fig. B.1. Corner plot of the inference of modeling mock system #1 using only the lensed quasar position and time delay. The SIE model and Power-law model are adopted on the left and right, separately. The blue lines indicate the true values in the simulation.

Qualitatively similar behavior is observed for the other systems, but we do not report inferred parameters in those cases due to the difficulty to achieve convergence of the MCMC chains for the power-law model. This is due to the degeneracies observed between q, γ and H₀ which implies a sampling of a large region of the parameter space, enforcing exploration of parameter values for which results of LENSTRONOMY modeling has not been fully tested (e.g., γ >  2.5).

All Tables

All Figures

Fig. 1. Sensitivity of the inferred Hubble constant as a function of fractional change in the measured lens velocity dispersion, σv (see Eq. (15)). Each color corresponds to one of the seven strong lens systems of the current TDCOSMO sample. The dotted lines display the best linear fit to the data. The joint inference performed on the seven lenses is shown in black. The error bars correspond to the 16th and 84th percentile of the posterior distributions. The two bottom panels show the sensitivity of H0 to a change in the measured lens velocity dispersion for power-law (left) and composite (right) models independently. The sensitivity of the joint inference, ⟨ξ⟩ is indicated on each panel.

In the text

Fig. 2. Effective radius θeff, Einstein radius θE and radius of the spectroscopic aperture θaperture of the TDCOSMO lenses. We show the ratios of these three quantities and the corresponding H0 value inferred for each system. We do not observe significant correlations between the characteristic sizes of the lens, the spectroscopic aperture and H0. The horizontal lines indicate the latest H0LiCOW 2019 (dotted orange, Wong et al. 2020) and Planck (dashed blue, Planck Collaboration VI 2020) results along with the 1σ uncertainties.

In the text

	Fig. 3. Hubble constant as a function of the measured velocity dispersion of the main lens. The horizontal lines indicate the latest H0LiCOW 2019 (dotted orange, Wong et al. 2020) and Planck (dashed blue, Planck Collaboration VI 2020) results along with the 1σ uncertainties.
In the text

Fig. 4. Measured Hubble constant, before (upper panel) and after (lower panel) correction for the mass along the line of sight as a function of the estimated external convergence. $H_0^{\text{uncorr}}$ and $H_0^{\text{corr}}$ are related according to Eq. (16). The dashed black lines show the best linear fit, and the shaded gray envelopes correspond to the 1σ uncertainties. The dotted blue lines represent the relation expected from the theory between $H_0^{\text{uncorr}}$, $H_0^{\text{corr}}$ and κext.

In the text

	Fig. 5. H0 constraints for the TDCOSMO lenses as a function of lens redshift before (top) and after (bottom) correction for the external convergence. The best linear fits and their 1σ envelopes are shown in shaded gray. The tentative (1.7σ significance) trend is not introduced by the LOS contribution as it is still visible before correcting for the external convergence.
In the text

Fig. 6. Marginalized H0 posteriors for power-law (left panel) and composite models (right panel). The cosmological inferences are for a flat ΛCDM cosmology with uniform priors. The posterior probability distributions for each individual system are shown with shaded color curves and the combined constraint from the six systems corresponds to the solid black curve. The legend indicates the median, 16th and 84th percentiles of the H0 distributions.

In the text

Fig. 7. Sample of simulated lenses: three pairs are generated from power-law, cored power-law, and composite lens models. The color scale is logarithmic and is the same for all images. Identifiers associated to each lens are also indicated. Refer to Sect. 6.2 for a description of these simulations. Model #6, although composite, is chosen so that the total mass profile resembles a power law in the region of the Einstein radius.

In the text

Fig. 8. Azimuthally averaged radial convergence profiles, for the different lens models applied to fit the sample of six mock lenses (Fig. 7). Upper part of each panel: true profiles are shown in dotted lines; power-law, cored power-law and composite models are shown in blue continuous, green dashed and red dot-dashed lines, respectively. The spectroscopic (square) aperture used for computing velocity dispersions is indicated as a vertical dotted line, and the true Einstein radius location is indicated as a vertical dashed line. The gray area encloses lensed quasar image positions. For each model, the inferred H0 values are indicated (in km s−1 Mpc−1), and must be compared to the input value $H_0^{\text{fiducial}}=70.7\text{km}{\text{s}}^{-1}{\text{Mpc}}^{-1}$. Lower part of each panel: relative error computed as (truth − model)/truth. The pixel size in the simulated images is $0\stackrel{″}{.}08$.

In the text

	Fig. B.1. Corner plot of the inference of modeling mock system #1 using only the lensed quasar position and time delay. The SIE model and Power-law model are adopted on the left and right, separately. The blue lines indicate the true values in the simulation.
In the text