Issue 
A&A
Volume 639, July 2020



Article Number  A101  
Number of page(s)  19  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/201937351  
Published online  16 July 2020 
TDCOSMO
I. An exploration of systematic uncertainties in the inference of H_{0} from timedelay cosmography
^{1}
Institute of Physics, Laboratory of Astrophysics, Ecole Polytechnique Fédérale de Lausanne (EPFL), Observatoire de Sauverny, 1290 Versoix, Switzerland
email: martin.millon@epfl.ch
^{2}
Department of Physics and Astronomy, University of California, Los Angeles, CA 90095, USA
^{3}
MaxPlanckInstitut für Astrophysik, KarlSchwarzschildStr. 1, 85748 Garching, Germany
^{4}
PhysikDepartment, Technische Universität München, JamesFranckStraße 1, 85748 Garching, Germany
^{5}
Academia Sinica Institute of Astronomy and Astrophysics (ASIAA), 11F of ASMAB, No.1, Section 4, Roosevelt Road, Taipei 10617, Taiwan
^{6}
Kavli Institute for Particle Astrophysics and Cosmology and Department of Physics, Stanford University, Stanford, CA 94305, USA
^{7}
Department of Physics, University of California, Davis, CA 95616, USA
^{8}
STAR Institute, Quartier Agora – Allée du six Août, 19c B4000 Liège, Belgium
^{9}
Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 2778583, Japan
^{10}
DARK, Niels Bohr Institute, Lyngbyvej 2, 2100 Copenhagen, Denmark
^{11}
Kapteyn Astronomical Institute, University of Groningen, PO Box 800, 9700 AV Groningen, The Netherlands
^{12}
Instituto de Física y Astronomía, Facultad de Ciencias, Universidad de Valparaíso, Avda. Gran Bretaña 1111, Valparaíso, Chile
^{13}
National Astronomical Observatory of Japan, 2211 Osawa, Mitaka, Tokyo 1818588, Japan
^{14}
Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK
^{15}
Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK
^{16}
Leiden Observatory, Leiden University, Niels Bohrweg 2, 2333 CA, Leiden, The Netherlands
^{17}
University of Portsmouth, Institute of Cosmology and Gravitation, Portsmouth PO1 3FX, UK
^{18}
INAF – Osservatorio Astronomico di Capodimonte, Salita Moiariello, 16, 80131 Napoli, Italy
^{19}
European Southern Observatory, KarlSchwarschildStr. 2, 85748 Garching, Germany
^{20}
Fermi National Accelerator Laboratory, PO Box 500, Batavia, IL 60510, USA
Received:
18
December
2019
Accepted:
20
May
2020
Timedelay cosmography of lensed quasars has achieved 2.4% precision on the measurement of the Hubble constant, H_{0}. As part of an ongoing effort to uncover and control systematic uncertainties, we investigate three potential sources: 1 stellar kinematics, 2 lineofsight effects, and 3 the deflector mass model. To meet this goal in a quantitative way, we reproduced the H0LiCOW/SHARP/STRIDES (hereafter TDCOSMO) procedures on a set of real and simulated data, and we find the following. First, stellar kinematics cannot be a dominant source of error or bias since we find that a systematic change of 10% of measured velocity dispersion leads to only a 0.7% shift on H_{0} from the seven lenses analyzed by TDCOSMO. Second, we find no bias to arise from incorrect estimation of the lineofsight effects. Third, we show that elliptical composite (stars + dark matter halo), powerlaw, and cored powerlaw mass profiles have the flexibility to yield a broad range in H_{0} values. However, the TDCOSMO procedures that model the data with both composite and powerlaw mass profiles are informative. If the models agree, as we observe in real systems owing to the “bulgehalo” conspiracy, H_{0} is recovered precisely and accurately by both models. If the two models disagree, as in the case of some pathological models illustrated here, the TDCOSMO procedure either discriminates between them through the goodness of fit, or it accounts for the discrepancy in the final error bars provided by the analysis. This conclusion is consistent with a reanalysis of six of the TDCOSMO (real) lenses: the composite model yields H_{0} = 74.0_{−1.8}^{+1.7} km s^{−1} Mpc^{−1}, while the powerlaw model yields 74.2_{−1.6}^{+1.6} km s^{−1} Mpc^{−1}. In conclusion, we find no evidence of bias or errors larger than the current statistical uncertainties reported by TDCOSMO.
Key words: gravitational lensing: strong / methods: data analysis
© ESO 2020
1. Introduction
The timedelay method applied to gravitationally lensed quasars (Refsdal 1964) provides a perhaps unrivalled combination of high sensitivity to the Hubble constant H_{0}, 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_{0} 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_{0} measurements depends on three main factors: 1 precise timedelay 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 timedelay measurements are available, for example, from the COSMOGRAIL collaboration, using longterm 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 timedelay challenge (Dobler et al. 2015; Liao et al. 2015; Bonvin et al. 2016). The timedelays 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 wellknown lensing degeneracies on the cosmological inference (e.g., Treu & Koopmans 2002). Finally, multiband widefield imaging and/or spectroscopy (e.g., Rusu et al. 2017; Sluse et al. 2019) was used to constrain the combined lensing effect of the lineofsight objects and largescale structures in a statistical way (Greene et al. 2013; Rusu et al. 2017). Tihhonova et al. (2018) also show that these estimates of the lineofsight 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 highprecision H_{0} 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^{1}. 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, (2.4% precision), is 3σ higher than the earlyUniverse 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_{0}, the tension with the earlyUniverse 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, timedelay 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, 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 timedelay cosmography of lensed quasars have decided to join forces with other scientists and form a new “umbrella” collaboration named TDCOSMO^{2} (TimeDelay 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 timedelay 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 lineofsight contribution. Third, in Sect. 6, we address the long standing issue of the masssheet 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 masssheet 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 timedelays, 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 highresolution images of lens systems comparable to the real observations, highprecision time delays (higher precision than those of real lenses so that the timedelay uncertainties are subdominant compared to the modeling uncertainties that we aim to quantify), and do not include the lineofsight 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. Timedelay cosmography and the masssheet degeneracy
Time delays in gravitationally lensed quasars provide a direct measurement of the socalled “timedelay 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}:
(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:
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:
We can also define the Fermat potential ϕ (Schneider 1985; Blandford & Narayan 1986) as
Using this definition, Eq. (2) reduces to
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 timedelay 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 MassSheet 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:
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 timedelay distance given by any model is affected by MST as follows :
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 spatiallyresolved 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 timedelay 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. Twodistance 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 timedelay 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:
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 :
Combining Eqs. (8) and (9), we obtain an expression for the angular diameter distance to the lens which is independent of the external convergence:
We immediately see that the angular diameter distance D_{d} varies as . The dependence of D_{d} to a change in the measurement of σ_{v} can therefore be computed analytically :
whereas D_{Δt} is left unchanged when varying the velocity dispersion. The final H_{0} measurement is obtained by combining these two distance measurements. As a consequence, the importance of the velocity dispersion in the final H_{0} value depends on the relative precision between the angular diameter distance and the timedelay distance, and on the mapping between the parameters. The D_{Δt} measurement is typically more constraining of H_{0} 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 multilensplane mass modeling. The singlelensplane 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 timedelay distance could be derived for these two lens systems, but the inference of D_{d} accounting for the multilens 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 powerlaw), 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 powerlaw, 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 wellidentified 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 powerlaw and composite models, with the addition of external shear.
2.3.1. Powerlaw model
Powerlaw models have a constant projected mass slope over the entire profile. The convergence of the powerlaw elliptical mass distribution (Barkana 1998) is described by :
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 θ_{1} and θ_{2} coordinates are along the major and minor axis respectively. The cored powerlaw 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 :
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 powerlaw 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 socalled 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:
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 powerlaw 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_{0} 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^{4}); 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 timedelay 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 timedelay 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_{0}. 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 signaltonoise 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 reallife cases as variations in the distortion of the images spanning a much larger significant radial range. Indeed, most of the HST data used in timedelay 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_{0} 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 massanisotropy 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_{0} on kinematics data defined by
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_{0} 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_{0}. 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_{0}, 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_{0} inference. Throughout this section, the H_{0} 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 and as the inferred H_{0} value of the system and its measured aperture velocity dispersion. The models used in Fig. 1 include both composite and powerlaw mass models^{3} 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_{0} 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 H_{0} to a change in the measured lens velocity dispersion for powerlaw (left) and composite (right) models independently. The sensitivity of the joint inference, ⟨ξ⟩ is indicated on each panel. 

Open with DEXTER 
The slope ξ quantifies the sensitivity of the inferred H_{0} 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_{0} 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_{0} by approximately 0.7%.
Summary of the H_{0} values (Col. 2) reported in Wong et al. (2020) and Shajib et al. (2020).
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_{0} (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_{0} is therefore almost insensitive to the kinematics.
In the case of PG 1115+080, the timedelay distance D_{Δt}, which does not depend on the kinematics data, has a much larger constraining power on H_{0} 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_{0} 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 xaxis 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 powerlaw 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 powerlaw 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 powerlaw 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 timedelay distance of SDSS J1206+4332 is better constrained by composite models ( Mpc at 7.1 % precision) than with powerlaw models ( Mpc at 11.2 % precision). The relative weight of the D_{Δt} compared to the D_{d} in the final value of H_{0} is therefore more important in the composite model case.
WFI 2033−4723 experiences the opposite behavior; it has tighter constrains with powerlaw models ( Mpc at 4.74 % precision) than with composite models ( 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_{0}.
5. Search for correlations between H_{0} and physically independent observables
The inference of H_{0} 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_{0}. Thus, any evidence of a correlation between these observables and the inferred value of H_{0} 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_{0} 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 followup. 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 H_{0} value inferred for each system. We do not observe significant correlations between the characteristic sizes of the lens, the spectroscopic aperture and H_{0}. 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. 

Open with DEXTER 
The Spearman’s rank correlation coefficient between H_{0} 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 pvalue) 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 xaxis 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_{0}. 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_{0}. 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 earlytype galaxies could give rise to measurable biases. For example, the socalled “tilt” of the fundamental mass plane is believed to arise primarily from the increase in darktostellar 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_{0} 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 pvalue 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. 

Open with DEXTER 
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_{0} 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 timedelay distance D_{Δt}, as shown in Eq. (7). Similarly, the effect of the external convergence on the inferred H_{0} can be written as :
where () is the value of H_{0} 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 MassSheet Degeneracy (MSD) not entirely removed by the measurement of κ_{ext}, we investigate the presence of correlation between the estimated κ_{ext} and the inferred H_{0} value for the seven lenses of the TDCOSMO sample. The top panel of Fig. 4 shows the relation between the H_{0} measurements before correction for the mass along the line of sight , 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 overdense regions (positive κ_{ext}) tend to have a higher than lenses in underdense 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^{−1} Mpc^{−1}, well compatible with the expected slope of . Both the uncorrected and corrected data are well fitted by our linear model, with a reduced χ^{2} 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. and 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 , and κ_{ext}. 

Open with DEXTER 
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 and κ_{ext}. In fact, the bestfit slope coefficient in this case is a^{corr} = −5.1 ± 23.7 km s^{−1} Mpc^{−1}, 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 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 wellfitted by a linear model both before and after the LOS correction with a reduced χ^{2} 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. H_{0} 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. 

Open with DEXTER 
6. Impact of the choice of families of mass model
In this section we quantify how much the inference on H_{0} depends on the choice of the mass density profile adopted for the lens modeling. We first use the six systems for which both powerlaw 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 “bulgehalo” conspiracy (Treu & Koopmans 2004; Dutton & Treu 2014).
Then, we carry out endtoend 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 galaxyscale lenses are earlytype 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 “bulgehalo conspiracy”. This bulgehalo 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 goodnessoffit 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 spatiallyextended 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 goodnessoffit requirement, even if they include time delays and stellarvelocitydispersion 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_{0} inference per model family
The TDCOSMO collaboration uses both composite and powerlaw models in their analysis, except for B1608+656 (see Sect. 2.3). Apart from this exception, the published estimates of H_{0} 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_{0} 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} ∈ [0, 150] km s^{−1} Mpc^{−1}, Ω_{m} ∈ [0.05, 0.5] and Ω_{m} = 1 − Ω_{Λ}.
Fig. 6. Marginalized H_{0} posteriors for powerlaw (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 H_{0} distributions. 

Open with DEXTER 
The H_{0} 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 when we use exclusively powerlaw models and when we use only composite model. This corresponds only to a 0.2% difference. Individual objects can have larger differences between powerlaw 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_{0} value. For example, WFI 2033−4723 has a higher H_{0} 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_{0} values, in practice they do not. In the following, we investigate with simulated lens systems the reasons why composite and powerlaw models provide comparable estimates of H_{0} in spite of allowing for flexibility. We also investigate under which circumstances gravitational lenses can be modeled with both composite and powerlaw models and still yield the same H_{0}.
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 powerlaw 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 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 galaxyscale 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 .
Three mass models, including powerlaw (ID #1, #2), cored powerlaw (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 quadlike 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 lensoptical axis (“crosslike” 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 powerlaw, cored powerlaw, 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. 

Open with DEXTER 
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 masstolight 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 powerlaw 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 powerlaw 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_{0} 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 , which was chosen randomly. For the timedelay 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 halflight 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^{4} (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 powerlaw, a cored powerlaw 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., ) 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_{0} (Ω_{m} is kept fixed to 0.27) from the timedelay distance alone. The lens velocity dispersion is computed only for comparison but is not included in the H_{0} 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
where k is the number of free parameters, 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 timedelay 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 powerlaws, 26 for the cored powerlaws and 29 for the composite models.
The recovered H_{0} 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_{0} 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.
BIC value, reduced χ^{2} of the image fit, measured H_{0}, tension relative to the true value of and integrated stellar velocity dispersion.
Residual maps of the lens modelling, i.e. normalized χ^{2} per pixel.
Interestingly, the core size of the cored powerlaw profile is well constrained by the data. Indeed, when a cored powerlaw 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 coresize is recovered within 2.2% accuracy and within < 3.0% precision with a cored powerlaw 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 powerlaw and cored powerlaw simulations. This is because, by construction, the lens light profile of these simulations does not necessarily correspond to their mass profile. In the powerlaw and cored powerlaw 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 powerlaw and composite models if we assume that the baryonic component of the composite model is traced by the arbitrary lens light in the powerlaw 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_{0} depends on the characteristics of the composite lenses. In both cases, the powerlaw 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_{0} is compatible with the true value for the lens #6, but in lens #5 it is biased toward lower H_{0} 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_{0} 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 χ^{2}< 1 can be obtained for all the mocks using a power law model, except for mock #6 for which the best reduced χ^{2} 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 powerlaw indices substantially different from the input one, yielding a bias on H_{0}, 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 timedelay 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_{0}. 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 wellknown problem in timedelay 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; powerlaw, cored powerlaw and composite models are shown in blue continuous, green dashed and red dotdashed 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 H_{0} values are indicated (in km s^{−1} Mpc^{−1}), and must be compared to the input value . Lower part of each panel: relative error computed as (truth − model)/truth. The pixel size in the simulated images is . 

Open with DEXTER 
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 unbiased 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 powerlaw (in the region where lensed images are formed) by means of a masssheet 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_{0}. 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 powerlaw models are discrepant, thus providing an indication that the powerlaw 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_{0}.
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_{0} is compatible with the true value for both the composite and powerlaw model. The fact that the two families of models are providing compatible H_{0} indicates that the convergence profile is wellrecovered 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 powerlaw models in the TDCOSMO analyses supports the hypothesis that, as a result of the bulgehalo 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_{0} 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 goodnessoffit and H_{0} 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 timedelay 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 timedelay cosmography:

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.

The TDCOSMO estimates of H_{0} 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.

Toy models based on simplified assumptions and constraints cannot offer any quantitative estimates of systematic errors given the current stateoftheart datasets 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 highresolution 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_{0}/H_{0} ∼ 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_{0} 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_{0} on kinematic data implies that systematic uncertainties in this area will have a subdominant impact on H_{0}.

No evidence is found for any correlation between the measured value of H_{0} and the external convergence estimated from galaxy number counts and numerical simulations. In contrast, if no external convergence is applied, H_{0} 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_{0} 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_{0} inferences and significant differences in image residuals when comparing powerlaw 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_{0} inferences between powerlaw and composite mass models.

The comparison of powerlaw and composite mass models allows us to quantify deviations in H_{0} due to our mass model assumptions. By using these two families of models and marginalising over them, the resulting H_{0} 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_{0} constraints from powerlaw and composite models of TDCOSMO lenses shows that the total mass profiles of galaxies are close to power laws due to the bulgehalo conspiracy. For the six lenses that have been analyzed with both powerlaw and composite models we find and respectively. The difference between the two model families is much smaller than the inferred statistical errors. The similar H_{0} from the different families of models thus made the current H_{0} 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 timedelay 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 masstolight ratio, (ii) explore in more detail the bulgehalo conspiracy based on high resolution data for local earlytype galaxies, (iii) explore the effect of allowing for gradients in stellar masstolight 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_{0}, they are worth undertaking.
The first H0LiCOW lens, namely B1608+656, was modeled with a powerlaw model and pixellated potential corrections, which were found to be small. A composite model was not applied, so we use only the powerlaw model in our analysis (see Suyu et al. 2010, for details).
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 AST1906976 to T.T. and by NASA through HST grants HSTGO10158, HSTGO12889, HSTGO14254, HSTGO15320, HSTGO15652. 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. AST1715611 and AST1907396. This work was supported by World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan. L.V.E.K. has been supported through an NWOVICI grant (project number 639.043.308). This research made use of Astropy, a communitydeveloped 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 ForemanMackey et al. (2013).
References
 Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Astropy Collaboration (PriceWhelan, A. M., et al.) 2018, AJ, 156, 123 [NASA ADS] [CrossRef] [Google Scholar]
 Auger, M. W., Treu, T., Bolton, A. S., et al. 2009, ApJ, 705, 1099 [NASA ADS] [CrossRef] [Google Scholar]
 Auger, M. W., Treu, T., Bolton, A. S., et al. 2010, ApJ, 724, 511 [NASA ADS] [CrossRef] [Google Scholar]
 Barkana, R. 1998, ApJ, 502, 531 [NASA ADS] [CrossRef] [Google Scholar]
 Barnabè, M., Czoske, O., Koopmans, L. V. E., Treu, T., & Bolton, A. S. 2011, MNRAS, 415, 2215 [NASA ADS] [CrossRef] [Google Scholar]
 Birrer, S., & Amara, A. 2018, Phys. Dark Univ., 22, 189 [NASA ADS] [CrossRef] [Google Scholar]
 Birrer, S., Amara, A., & Refregier, A. 2015, ApJ, 813, 102 [NASA ADS] [CrossRef] [Google Scholar]
 Birrer, S., Amara, A., & Refregier, A. 2016, J. Cosmol. Astropart. Phys., 2016, 020 [NASA ADS] [CrossRef] [Google Scholar]
 Birrer, S., Treu, T., & Rusu, C. E. 2019, MNRAS, 206 [Google Scholar]
 Blandford, R., & Narayan, R. 1986, ApJ, 310, 568 [NASA ADS] [CrossRef] [Google Scholar]
 Bonvin, V., Tewes, M., Courbin, F., et al. 2016, A&A, 585, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bonvin, V., Chan, J. H. H., Millon, M., et al. 2018, A&A, 616, A183 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bonvin, V., Millon, M., Chan, J. H. H., et al. 2019, A&A, 629, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cappellari, M. 2016, ARA&A, 54, 597 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, G. C. F., Fassnacht, C. D., Suyu, S. H., et al. 2019, MNRAS, 490, 1743 [NASA ADS] [CrossRef] [Google Scholar]
 Courbin, F., Bonvin, V., BuckleyGeer, E., et al. 2018, A&A, 609, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Czoske, O., Barnabè, M., Koopmans, L. V. E., Treu, T., & Bolton, A. S. 2012, MNRAS, 419, 656 [CrossRef] [Google Scholar]
 Ding, X., Liao, K., Treu, T., et al. 2017a, MNRAS, 465, 4634 [NASA ADS] [CrossRef] [Google Scholar]
 Ding, X., Treu, T., Suyu, S. H., et al. 2017b, MNRAS, 472, 90 [CrossRef] [Google Scholar]
 Ding, X., Treu, T., & Shajib, A. J. 2018, ArXiv eprints [arXiv:1801.01506] [Google Scholar]
 Dobler, G., Fassnacht, C. D., Treu, T., et al. 2015, ApJ, 799, 168 [NASA ADS] [CrossRef] [Google Scholar]
 Dutton, A. A., & Treu, T. 2014, MNRAS, 438, 3594 [CrossRef] [Google Scholar]
 Dutton, A. A., Brewer, B. J., Marshall, P. J., et al. 2011, MNRAS, 417, 1621 [NASA ADS] [CrossRef] [Google Scholar]
 Falco, E. E., Gorenstein, M. V., & Shapiro, I. I. 1985, ApJ, 289, L1 [NASA ADS] [CrossRef] [Google Scholar]
 ForemanMackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [NASA ADS] [CrossRef] [Google Scholar]
 Gavazzi, R., Treu, T., Rhodes, J. D., et al. 2007, ApJ, 667, 176 [NASA ADS] [CrossRef] [Google Scholar]
 Gomer, M. R., & Williams, L. L. R. 2019, ArXiv eprints [arXiv:1907.08638] [Google Scholar]
 Greene, Z. S., Suyu, S. H., Treu, T., et al. 2013, ApJ, 768, 39 [NASA ADS] [CrossRef] [Google Scholar]
 Hernquist, L. 1990, ApJ, 356, 359 [NASA ADS] [CrossRef] [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
 Jee, I., Komatsu, E., & Suyu, S. H. 2015, J. Cosmol. Astropart. Phys., 2015, 033 [NASA ADS] [CrossRef] [Google Scholar]
 Jee, I., Komatsu, E., Suyu, S. H., & Huterer, D. 2016, J. Cosmol. Astropart. Phys., 2016, 031 [NASA ADS] [CrossRef] [Google Scholar]
 Jee, I., Suyu, S. H., Komatsu, E., et al. 2019, Science, 365, 1134 [CrossRef] [Google Scholar]
 Keeton, C. R. 2001, ArXiv eprints [arXiv:astroph/0102340] [Google Scholar]
 Keeton, C. R. 2011, GRAVLENS: Computational Methods for Gravitational Lensing [Google Scholar]
 Knox, L., & Millea, M. 2020, Phys. Rev. D, 101, 043533 [CrossRef] [Google Scholar]
 Kochanek, C. S. 2002, ApJ, 578, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Kochanek, C. S. 2006, in SaasFee Advanced Course 33: Gravitational Lensing:Strong, Weak and Micro, eds. G. Meylan, P. Jetzer, P. North, et al., 91 [NASA ADS] [CrossRef] [Google Scholar]
 Kochanek, C. S. 2020, MNRAS, 493, 1725 [CrossRef] [Google Scholar]
 Kochanek, C. S., Keeton, C. R., & McLeod, B. A. 2001, ApJ, 547, 50 [NASA ADS] [CrossRef] [Google Scholar]
 Koopmans, L. V. E. 2004, ArXiv eprints [arXiv: astroph/0412596] [Google Scholar]
 Koopmans, L. V. E., Treu, T., Fassnacht, C. D., Bland ford, R. D., & Surpi, G. 2003, ApJ, 599, 70 [NASA ADS] [CrossRef] [Google Scholar]
 Koopmans, L. V. E., Treu, T., Bolton, A. S., Burles, S., & Moustakas, L. A. 2006, ApJ, 649, 599 [NASA ADS] [CrossRef] [Google Scholar]
 Koopmans, L. V. E., Bolton, A., Treu, T., et al. 2009, ApJ, 703, L51 [NASA ADS] [CrossRef] [Google Scholar]
 Lagattuta, D. J., Fassnacht, C. D., Auger, M. W., et al. 2010, ApJ, 716, 1579 [NASA ADS] [CrossRef] [Google Scholar]
 Liao, K., Treu, T., Marshall, P., et al. 2015, ApJ, 800, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Liesenborgs, J., & De Rijcke, S. 2012, MNRAS, 425, 1772 [NASA ADS] [CrossRef] [Google Scholar]
 Mukherjee, S., Koopmans, L. V. E., Metcalf, R. B., et al. 2018, MNRAS, 479, 4108 [CrossRef] [Google Scholar]
 Mukherjee, S., Koopmans, L. V. E., Metcalf, R. B., et al. 2019, ArXiv eprints [arXiv:1901.01095] [Google Scholar]
 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493 [NASA ADS] [CrossRef] [Google Scholar]
 Pandey, S., Raveri, M., & Jain, B. 2020, Phys. Rev. D, 102, 023505 [CrossRef] [Google Scholar]
 Paraficz, D., & Hjorth, J. 2009, A&A, 507, L49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VI. 2020, A&A, in press, https://doi.org/10.1051/00046361/201833910 [Google Scholar]
 Read, J. I., Saha, P., & Macciò, A. V. 2007, ApJ, 667, 645 [NASA ADS] [CrossRef] [Google Scholar]
 Refsdal, S. 1964, MNRAS, 128, 307 [NASA ADS] [CrossRef] [Google Scholar]
 Riess, A. G., Casertano, S., Yuan, W., Macri, L. M., & Scolnic, D. 2019, ApJ, 876, 85 [NASA ADS] [CrossRef] [Google Scholar]
 Rusu, C. E., Fassnacht, C. D., Sluse, D., et al. 2017, MNRAS, 467, 4220 [NASA ADS] [CrossRef] [Google Scholar]
 Rusu, C. E., Wong, K. C., Bonvin, V., et al. 2020, MNRAS, in press, [arXiv:1905.09338] [Google Scholar]
 Saha, P. 2000, AJ, 120, 1654 [NASA ADS] [CrossRef] [Google Scholar]
 Saha, P., & Williams, L. L. R. 2006, ApJ, 653, 936 [NASA ADS] [CrossRef] [Google Scholar]
 Schneider, P. 1985, A&A, 143, 413 [Google Scholar]
 Schneider, P., & Sluse, D. 2013, A&A, 559, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schneider, P., & Sluse, D. 2014, A&A, 564, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schneider, P., Ehlers, J., & Falco, E. E. 1992, Gravitational Lenses [Google Scholar]
 Sérsic, J. L. 1963, Boletin de la Asociacion Argentina de Astronomia La Plata Argentina, 6, 41 [Google Scholar]
 Shajib, A. J., Treu, T., & Agnello, A. 2018, MNRAS, 473, 210 [NASA ADS] [CrossRef] [Google Scholar]
 Shajib, A. J., Birrer, S., Treu, T., et al. 2020, MNRAS, 494, 6072 [CrossRef] [Google Scholar]
 Sluse, D., Rusu, C. E., & Fassnacht, C. D. 2019, MNRAS, 2136 [Google Scholar]
 Sonnenfeld, A. 2018, MNRAS, 474, 4648 [NASA ADS] [CrossRef] [Google Scholar]
 Sonnenfeld, A., Leauthaud, A., Auger, M. W., et al. 2018, MNRAS, 481, 164 [NASA ADS] [CrossRef] [Google Scholar]
 Suyu, S. H. 2012, MNRAS, 426, 868 [NASA ADS] [CrossRef] [Google Scholar]
 Suyu, S. H., Marshall, P. J., Blandford, R. D., et al. 2009, ApJ, 691, 277 [NASA ADS] [CrossRef] [Google Scholar]
 Suyu, S. H., Marshall, P. J., Auger, M. W., et al. 2010, ApJ, 711, 201 [NASA ADS] [CrossRef] [Google Scholar]
 Suyu, S. H., Treu, T., Hilbert, S., et al. 2014, ApJ, 788, L35 [NASA ADS] [CrossRef] [Google Scholar]
 Suyu, S. H., Bonvin, V., Courbin, F., et al. 2017, MNRAS, 468, 2590 [NASA ADS] [CrossRef] [Google Scholar]
 Tagore, A. S., Barnes, D. J., Jackson, N., et al. 2018, MNRAS, 474, 3403 [NASA ADS] [CrossRef] [Google Scholar]
 Tihhonova, O., Courbin, F., Harvey, D., et al. 2018, MNRAS, 477, 5657 [NASA ADS] [CrossRef] [Google Scholar]
 Treu, T., & Koopmans, L. V. E. 2002, MNRAS, 337, L6 [NASA ADS] [CrossRef] [Google Scholar]
 Treu, T., & Koopmans, L. V. E. 2004, ApJ, 611, 739 [NASA ADS] [CrossRef] [Google Scholar]
 Treu, T., & Marshall, P. J. 2016, A&ARv, 24, 11 [NASA ADS] [CrossRef] [Google Scholar]
 van de Ven, G., Mandelbaum, R., & Keeton, C. R. 2009, MNRAS, 398, 607 [NASA ADS] [CrossRef] [Google Scholar]
 Vegetti, S., Koopmans, L. V. E., Auger, M. W., Treu, T., & Bolton, A. S. 2014, MNRAS, 442, 2017 [CrossRef] [Google Scholar]
 Verde, L., Treu, T., & Riess, A. G. 2019, Nat. Astron., 3, 891 [NASA ADS] [CrossRef] [Google Scholar]
 Wagner, J. 2018, A&A, 620, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wertz, O., Orthen, B., & Schneider, P. 2018, A&A, 617, A140 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Witt, H. J., Mao, S., & Keeton, C. R. 2000, ApJ, 544, 98 [NASA ADS] [CrossRef] [Google Scholar]
 Wojtak, R., & Agnello, A. 2019, MNRAS, 486, 5046 [NASA ADS] [CrossRef] [Google Scholar]
 Wong, K. C., Suyu, S. H., Auger, M. W., et al. 2017, MNRAS, 465, 4895 [NASA ADS] [CrossRef] [Google Scholar]
 Wong, K. C., Suyu, S. H., Chen, G. C. F., et al. 2020, MNRAS, in press [arXiv:1907.04869] [Google Scholar]
 Wucknitz, O. 2002, MNRAS, 332, 951 [NASA ADS] [CrossRef] [Google Scholar]
 Xu, D., Sluse, D., Schneider, P., et al. 2016, MNRAS, 456, 739 [NASA ADS] [CrossRef] [Google Scholar]
 Yıldırım, A., Suyu, S. H., & Halkola, A. 2020, MNRAS, 493, 4783 [CrossRef] [Google Scholar]
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: halflight radius, effective Einstein radius, core radius in the case of cored powerlaw profiles, and scale radius of the dark matter profile for composite models. The ratio of the lens halflight 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 concaveupward (convexdownward) radial convergence profile has curvature greater than one (lower than one), and a perfect powerlaw 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_{0}. 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_{0} even if the slope differs significantly from isothermal, as illustrated with our simulated galaxy #6.
Key properties of the simulated lenses described in Sect. 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 on the pointsource positions, on the lensing galaxy centroid, and the same uncertainty on the timedelay 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 Powerlaw 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 . A notebook implementing our fitting strategy is available online^{5}.
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_{0} 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_{0} 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 pointsource 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_{0}.
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 Powerlaw model are adopted on the left and right, separately. The blue lines indicate the true values in the simulation. 

Open with DEXTER 
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 powerlaw model. This is due to the degeneracies observed between q, γ and H_{0} 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
Summary of the H_{0} values (Col. 2) reported in Wong et al. (2020) and Shajib et al. (2020).
BIC value, reduced χ^{2} of the image fit, measured H_{0}, tension relative to the true value of and integrated stellar velocity dispersion.
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 H_{0} to a change in the measured lens velocity dispersion for powerlaw (left) and composite (right) models independently. The sensitivity of the joint inference, ⟨ξ⟩ is indicated on each panel. 

Open with DEXTER  
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 H_{0} value inferred for each system. We do not observe significant correlations between the characteristic sizes of the lens, the spectroscopic aperture and H_{0}. 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. 

Open with DEXTER  
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. 

Open with DEXTER  
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. and 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 , and κ_{ext}. 

Open with DEXTER  
In the text 
Fig. 5. H_{0} 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. 

Open with DEXTER  
In the text 
Fig. 6. Marginalized H_{0} posteriors for powerlaw (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 H_{0} distributions. 

Open with DEXTER  
In the text 
Fig. 7. Sample of simulated lenses: three pairs are generated from powerlaw, cored powerlaw, 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. 

Open with DEXTER  
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; powerlaw, cored powerlaw and composite models are shown in blue continuous, green dashed and red dotdashed 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 H_{0} values are indicated (in km s^{−1} Mpc^{−1}), and must be compared to the input value . Lower part of each panel: relative error computed as (truth − model)/truth. The pixel size in the simulated images is . 

Open with DEXTER  
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 Powerlaw model are adopted on the left and right, separately. The blue lines indicate the true values in the simulation. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.