A model-independent determination of the Hubble constant from lensed quasars and supernovae using Gaussian process regression

TL;DR

This paper presents a model-independent approach using Gaussian process regression to measure the Hubble constant from lensed quasars and supernovae, achieving high precision without assuming a cosmological model.

Contribution

It introduces a novel Gaussian process-based method to determine H0 directly from lensing and supernova data, avoiding reliance on the ΛCDM model.

Findings

H0 measured as 72.2 ± 2.1 km/s/Mpc assuming flat universe

H0 measured as 73.0 +2.8/-3.0 km/s/Mpc with curvature prior

Achieves 2.9% precision in H0 measurement

Abstract

Strongly lensed quasar systems with time delay measurements provide "time delay distances", which are a combination of three angular diameter distances and serve as powerful tools to determine the Hubble constant $H_0$. However, current results often rely on the assumption of the $\Lambda$CDM model. Here we use a model-independent method based on Gaussian process to directly constrain the value of $H_0$. By using Gaussian process regression, we can generate posterior samples of unanchored supernova distances independent of any cosmological model and anchor them with strong lens systems. The combination of a supernova sample with large statistics but no sensitivity to $H_0$ with a strong lens sample with small statistics but $H_0$ sensitivity gives a precise $H_0$ measurement without the assumption of any cosmological model. We use four well-analyzed lensing systems from the state-of-art lensing program H0LiCOW and the Pantheon supernova compilation in our analysis. Assuming the Universe is flat, we derive the constraint $H_0=72.2 \pm 2.1\,$km/s/Mpc, a precision of $2.9\%$. Allowing for cosmic curvature with a prior of $\Omega_{k}=[-0.2,0.2]$, the constraint becomes $H_0=73.0_{-3.0}^{+2.8}\,$km/s/Mpc.

Tables

Table 1. Table 1 Lens and source redshifts for the four strong lens systems ordered by distance (see Taubenberger et al. 2019 and the references therein).

Order	Name	$z_L$	$z_S$
1	RXJ1131-1231	0.295	0.654
2	HE 0435-1223	0.4546	1.693
3	B1608+656	0.6304	1.394
4	SDSS 1206+4332	0.745	1.789

Equations

$$\Delta t=D_{\Delta t}\Delta \varphi ({\mathbf{\xi }}_{lens}),$$

$$D_{\Delta t}=(1+z_l)\frac{D_l{D}_s}{D_{ls}},$$

$$P(D_{\Delta t})=\frac{1}{\sqrt{2\pi }(x-{\lambda }_D){\sigma }_D}\exp \left\{-\frac{[\ln (x-{\lambda }_D)-{\mu }_D{]}^2}{2{\sigma }_D^2}\right\},$$

$$⟨\gamma (z_1)\gamma (z_2)⟩={\sigma }_f^2\exp \{-[s(z_1)-s(z_2){]}^2/(2{\ell }^2)\},$$

$$H_0{D}^L(z)=(1+z){\Omega }_k^{-1/2}\sinh \left[{\Omega }_k^{1/2}{\int }_0^{z}d{z}^{\prime }/[H(z^{\prime })/H_0]\right],$$

$$\begin{array}{rl}{D}_{ls}& =D_s\sqrt{1+{\Omega }_k(1+z_l{)}^2(H_0{D}_l{)}^2}\\ & -\frac{1+z_l}{1+z_s}{D}_l\sqrt{1+{\Omega }_k(1+z_s{)}^2(H_0{D}_s{)}^2}.\end{array}$$

Full text

1. Introduction

The standard cosmological model of $\Lambda$CDM has achieved great successes in explaining a wide array of cosmological observations, including the expansion history from distances of Type Ia supernovae (SNe) and baryon acoustic oscillations (BAO), the cosmic growth history from galaxy surveys, and the cosmic microwave background (CMB). However, there are some discrepancies, the clearest of which is the value of the Hubble constant $H_{0}$ as measured locally, from SNe Ia calibrated by Cepheid variable stars (Riess et al., 2019) or the tip of the red giant branch of stars (Freedman et al., 2019), and as derived cosmologically, from CMB (Planck Collaboration et al., 2018) or BAO with or without CMB or indeed any early universe information (Addison et al., 2018; Macaulay et al., 2019; Cuceu et al., 2019).

Note that the cosmological estimation of $H_{0}$ is not a direct measurement but a derivation as one of a number of cosmological model parameters. Usually the model assumed is the $\Lambda$CDM when inferring $H_{0}$. Therefore, the discrepancy reveals either new physics beyond the standard model or unknown systematic errors in the observations. To better understand this tension problem, one may use new cosmological approaches to determine $H_{0}$ model-dependently or even model-independently. For example, gravitational wave events from inspiraling compact objects plus their electromagnetic counterparts are very promising for an alternate local direct measurement (Abbott et al., 2017).

Strong gravitational lensing systems with time delays offer a method for cosmological determination of $H_{0}$ that is partly direct (time delays are proportional to $1/H_{0}$) and partly derived from a cosmological parameter fit, but independent of both the local and the early universe. A typical system consists of a lensed quasar at cosmological distance, lensed by a foreground elliptical galaxy, forming multiple images of the active galactic nucleus (AGN) and the arcs of its host galaxy. With years of monitoring on the AGN light curves, one can measure the time delay between any two images corresponding to different paths, due to the geometric and Shapiro effects following the Fermat principle. The time delay thus depends on both the geometry of the Universe and the gravity field of the lens galaxy (lensing potential).

With ancillary data to measure the lensing potential, for example high-resolution imaging, stellar dynamic measurements, and line-of-sight environment measurements, we can measure the geometry of the Universe in terms of the “time delay distance” $D_{\Delta t}$. This is a ratio of three angular diameter distances, and depends on $H_{0}$ and also other cosmological parameters. In this way strong lensing has been used to determine $H_{0}$. The lensing program H0LiCOW (Suyu et al., 2017) with the Hubble Space Telescope can measure $D_{\Delta t}$ with several percent precision for each system, from percent level measurements of both time delays and lens modeling.

With only four well-measured systems, they constrained $H_{0}$ with $3\%$ precision in a flat $\Lambda$CDM model (Birrer et al., 2019), and looser precision in other cosmological models where the dark energy equation of state is allowed to differ from $-1$. A more ambitious aim (Suyu et al., 2017) is to measure $H_{0}$ within $1\%$ uncertainty such that the current tension may be investigated with high statistical significance. Though the measured $H_{0}$ would be very important for the community, it is worth noting that these measurements rely on assuming certain cosmological models. In Taubenberger et al. (2019), they found the results could be stable with respect to different models when combining strong lensing with SNe Ia as another cosmological probe111One might wonder whether some of this could possibly be due to using JLA SNe data in some tension with lensing data, such that their combination pulled in on $\Lambda$CDM where all the models considered were similar. Wong et al. (2019) did find that using Pantheon data increased the spread in the mean $H_{0}$ value from flat $\Lambda$CDM to flat $w_{0}w_{a}$CDM from 0.6 for JLA to 1.4 for Pantheon, which is still not significant compared to the uncertainties. .

The Inverse distance ladder method (Aubourg et al., 2015; Cuesta et al., 2015) provides a more model-independent way to infer $H_{0}$. The idea is to anchor the relative distances from SNe Ia with the absolute distance measurements from other cosmological approaches. In this work, we propose a model-independent approach to determining $H_{0}$ using a Gaussian process (GP) and by anchoring the SNe Ia with strong lensing. This combines the strengths of each technique: the large sample of SNe constrains the individual distances, freeing the smaller sample of strong lenses to lock down $H_{0}$ from the distance ratio. Having $H_{0}$ determined in this way could be more direct and informative for understanding the Hubble tension problem since the results are not biased by any model or parametric assumption – and are independent of both the local universe and the very high redshift, early universe. This gives a new angle on the problem.

We should note that a recent work by Collett et al. (2019) used SNe Ia and strong lensing to determine $H_{0}$ by implementing fourth order polynomial fitting to the supernovae data. Comparison of our results will be interesting, especially because using parametric approaches such as polynomial fitting to different cosmology data can be prone to issues such as instabilities, e.g. see (Jönsson et al., 2004; Shafieloo et al., 2006; Shafieloo, 2007; Holsclaw et al., 2010a, b, 2011; Shafieloo, 2012; Shafieloo et al., 2012). It is therefore useful to crosscheck, especially as using model-independent or non-parametric reconstruction approaches, as we do in this work, reduces potential bias from the forms of the parametric or model assumptions.

This Letter is organized as follows: In Section 2 we introduce the time delay cosmology and the latest lensing data. Then we combine lensing and SNe Ia to give a model-independent constraint on $H_{0}$ in Section 3. We summarize and make discussions in Section 4. Throughout this Letter, we use the natural units of c=G=1 in all equations. $H_{0}$ in units of km/s/Mpc will be recovered in the results.

2. Time delay distances from lensing

According to the theory of strong lensing, the time delay between two images of the AGN is determined by both the geometry of the Universe and the gravity field of the lens galaxy through:

[TABLE]

where $\Delta\phi=[(\boldsymbol{\theta}_{A}-\boldsymbol{\beta})^{2}/2-\psi(\boldsymbol{\theta}_{A})-(\boldsymbol{\theta}_{B}-\boldsymbol{\beta})^{2}/2+\psi(\boldsymbol{\theta}_{B})]$ is the difference of Fermat potentials at two images, $\boldsymbol{\theta}_{A}$ and $\boldsymbol{\theta}_{B}$ denote the angular positions of the images while $\boldsymbol{\beta}$ denotes the position of the source (supposing it is unlensed). $\psi$ is the two-dimensional lensing potential via the Poisson equation $\nabla^{2}\psi=2\kappa$, where $\kappa$ is the surface mass density of the lens in units of critical density $\Sigma_{\mathrm{crit}}=D_{s}/(4\pi D_{l}D_{ls})$. Assuming a lens model, $\Delta\phi$ is determined by the parameters $\boldsymbol{\xi}_{lens}$ therein. $D_{\Delta t}$ is the “time delay distance” formed from three angular diameter distances:

[TABLE]

where $l,s$ stands for lens and source. Note that $D_{\Delta t}$ is primarily sensitive to $H_{0}$, providing a powerful and independent way to determine it.

The $\Delta t$ can be measured by comparing the light curve shift of the two AGN images. With current techniques and the quality of the light curves, the precision of $\Delta t$ can be up to percent levels (see, e.g., Tewes et al. (2013) among many others). The upcoming LSST will discover thousands of lensed quasars, some of which will have well-measured light curves. The Time Delay Challenge (Liao et al., 2015) showed about 400 systems will have robust time delay measurements with average precision $3\%$, making time delay cosmography very promising. Meanwhile, the Fermat potentials can be measured by high-resolution imaging from space telescopes or ground-based adaptive optics, together with the stellar dynamics and the structure along the line-of-sight. The precision of $\Delta\phi$ is comparable with that of $\Delta t$, resulting in few percent level determination of $D_{\Delta t}$ for each system.

The state-of-art lensing project H0LiCOW aims at measuring $H_{0}$ with precision $\lesssim 1\%$ based on a small sample of well-observed lenses in the near future. Currently, under a flat $\Lambda$CDM model, they get a result with $3\%$ precision including systematics, from only 4 lenses: RXJ1131-1231, HE 0435-1223, B1608+656 and SDSS 1206+4332. Table 1 lists the lens and source redshifts for these systems. The data from another two lenses should also be released soon (Chen et al., 2019; Rusu et al., 2019; Wong et al., 2019). The posteriors of the time delay distances for the four lenses are given in the H0LiCOW papers and website222http://www.h0licow.org. For the first three of them, the posterior probability distributions of $D_{\Delta t}$ are described by the following analytic fit:

[TABLE]

where $x=D_{\Delta t}/(1\,{\rm Mpc})$ and the parameters $(\lambda_{D},\sigma_{D},\mu_{D})$ can be found in Table 3 of Bonvin (2017). For the fourth lens SDSS 1206+4332, the posterior is given in the form of the Markov chain Monte Carlo (MCMC) (Birrer et al., 2019).

However, since the angular distances depend on the cosmological model, the time delay distances also change in different cosmological models. The value of $H_{0}$ can then vary considerably between models. Taubenberger et al. (2019) claim to achieve stable results for $H_{0}$, i.e. insensitive to cosmological model, by using strong lensing to anchor Type Ia supernova data. However, that analysis was still within a set of cosmological models. Furthermore, it used the JLA supernovae data, which mildly prefer a dark energy equation of state $w>-1$ (see Fig. 16 of Betoule et al. (2014)), while strong lenses mildly prefer $w<-1$ (Wong et al. (2019)). The combination therefore may lie close to the cosmological constant $w=-1$ where the models they consider are equivalent. Thus, it is useful to try a model independent analysis for $H_{0}$, with the more recent Pantheon supernovae data set.

3. Methodology and results

SNe Ia observations led to the discovery of cosmic acceleration. They are incisive probes of cosmology through determining the shape of the cosmic distance-redshift relation. That is, they determine distances in a relative sense, but the absolute distance is convolved with a combination of the absolute magnitude of SNe Ia and the Hubble constant. The SNe data can be combined with (“anchored by”) an absolute distance probe (for example Cepheid variable stars, the tip of the red giant branch stars, or gravitational wave sirens in the local universe) to form an absolute distance probe. In addition, absolute distance probes at cosmological distances, e.g. strong lens systems, can anchor SN. This greatly benefits the leverage of the absolute probe in constraining cosmology since SNe tend to be much more numerous as well as very good probes of dark energy properties, mapping a wide range of cosmic expansion history. Thus the combination of strong lensing time delays and SNe can determine both $H_{0}$ and the cosmic expansion history.

3.1. Data and Method

For the SNe data we use the most recent and largest data set, the Pantheon compilation (Scolnic et al., 2018). To combine the Pantheon SNe and the H0LiCOW strong lenses datasets, we generate a posterior sampling of the $H_{0}$-independent quantity $H_{0}D^{L}(z)$ from the Pantheon dataset. To do the posterior sampling in a manner independent of a cosmological model, we use GP regression (Holsclaw et al., 2010a, b, 2011; Shafieloo et al., 2012, 2013). The GP regression used here is based on the GPHist code (Kirkby & Keeley, 2017) first used in Joudaki et al. (2018). GP regression works by generating a family of functions over an infinite dimensional function space as determined by a kernel. We use

[TABLE]

where $\sigma_{f}$ and $\ell$ are hyperparameters, respectively characterizing the amplitude of variations with redshift and their correlation scale. The hyperparameters play important roles for both physical insight and error control, and must be fit for not fixed. The priors on the GP hyperparameters are scale-invariant i.e. flat in the log of the hyperparameters. Since the dimmensionality of the hyperparameters is small, we directly integrate over the hyperparameters.

We then use GP on the SNe data to generate expansion histories $H(z)/H_{0}$ where $\gamma(z)=\ln([H^{\rm fid}(z)/H_{0}]/[H(z))/H_{0}])$. Here $H^{\rm fid}(z)/H_{0}$ is taken to be the best fit $\Lambda$CDM model for the Pantheon data and serves the role of the mean function for GP regression. Such prewhitening is standard practice and extensive tests show the resulting median inference does not depend on the details of the mean function (Shafieloo et al., 2012, 2013; Aghamousa et al., 2017). With $H(z)/H_{0}$ in hand, we can calculate the unanchored SNe luminosity distances

[TABLE]

(where $\sinh$ is a complete function valid for all signs of $\Omega_{k}$) and any corresponding angular diameter distances $H_{0}D^{A}=H_{0}D^{L}/(1+z)^{2}$ we will need later for the strong lensing systems, where $\Omega_{k}=1-\Omega_{\rm total}$ is the curvature energy density in units of the critical density. The likelihood of how well these SNe distances fit the Pantheon data are then used as weights in randomly selecting 1000 samples used for $D^{A}$ in the strong lens analysis.

Example GP curves are shown in Fig. 1. The data is dense and precise enough to provide a well constrained distance-redshift relation. The spread is about 2% at $z=1$, and the GP is not constrained to follow the input mean function – in fact it deviates from it by about 1.5% at $z=1$. Note that the GP covers the full range of the strong lensing system redshifts so there is no extrapolation needed.

To summarize the method for constraining $H_{0}$:

Draw 1000 unanchored luminosity distance curves $H_{0}D^{L}$ from the GP fit to the SNe data, and convert to unanchored angular diameter distances $H_{0}D^{A}$; 2. 2.

Evaluate the values of each of the 1000 $H_{0}D^{A}$ curves at the lens and source redshifts of the four strong lens systems to calculate 1000 values of $H_{0}D_{\Delta t}$ for each system using $H_{0}D_{\Delta t}=(1+z_{l})(H_{0}D_{l})(H_{0}D_{s})/(H_{0}D_{ls})$; 3. 3.

Compute the likelihood, for each of the 1000 realizations, from the H0LiCOW’s $D_{\Delta t}$ data for each lens system for many values of $H_{0}$; 4. 4.

Multiply the four likelihoods to form the full likelihood for each realization, for each value of $H_{0}$; 5. 5.

Marginalize over the realizations to form the posterior distribution of $H_{0}$.

Note that to obtain the angular diameter distance between the lens and the source we use the standard distance relation (Weinberg, 1972)

[TABLE]

Note that for a spatially flat universe one simply has $H_{0}D^{L}(z)=(1+z)\int_{0}^{z}dz^{\prime}/[H(z^{\prime})/H_{0}]$ and $D_{ls}=D_{s}-[(1+z_{l})/(1+z_{s})]\,D_{l}$.

3.2. Results

The final posterior distribution for $H_{0}$ in the flat universe case is shown in Fig. 2. Our model-independent constraint is $H_{0}=72.2\pm 2.1\,$ km/s/Mpc (median value plus the $16^{th}$ and $84^{th}$ percentiles around this). This can be compared to the time delay plus SNe result of Taubenberger et al. (2019) of $73.1^{+2.1}_{-2.2}$ km/s/Mpc within $\Lambda$CDM or $73.1\pm 3.0$ km/s/Mpc within $w_{0}w_{a}$CDM, or of Wong et al. (2019) of $73.6^{+1.6}_{-1.8}$ km/s/Mpc and $75.0^{+2.2}_{-2.3}$ km/s/Mpc respectively. Note that the first set uses JLA SNe rather than Pantheon and the second set uses six lens systems rather than four. In addition, the first set uses $D_{\Delta t}$ only while the second uses the combination of $D_{\Delta t}$ and $D_{l}$. The main point, however, is that the uncertainties for a model independent analysis combining strong lensing time delays and a wide ranging distance probe such as supernovae can be comparable to those assuming a specific model, while reducing possible bias.

For the case where spatial curvature is included, we set the flat prior on $\Omega_{k}=[-0.2,0.2]$ as in Taubenberger et al. (2019). The posterior distribution is shown in Fig. 2 also, and the marginalized distributions give $H_{0}=73.0_{-3.0}^{+2.8}\,$ km/s/Mpc and $\Omega_{k}=0.07_{-0.14}^{+0.09}$. The weaker constraint is due to the covariance between $H_{0}$ and $\Omega_{k}$. While consistent with the flat universe results, note the posterior distribution is more nonGaussian.

When trying to resolve tensions, it is worthwhile checking for the internal consistency of data combinations used. We explore first consistency between the time delay distances and SNe distances, and then consistency within the set of time delay systems, assuming a flat universe. To check the consistency between the distances from the SNe reconstructions and the time delay distances from the strong lenses, we use the GP fit to the expansion history derived from the SNe luminosity distances to compute predicted time delay distances with appropriate $z_{l}$ and $z_{s}$. As the SNe luminosity distances are unanchored, the conversion from $H_{0}D_{\Delta t}$ to $D_{\Delta t}$ stretches out the joint contours along this degeneracy. Figure 3 shows the results of this consistency check. The distances are indeed consistent, and hence the combination of the two probes is justified.

3.3. Consistency

One could also evaluate the consistency between the best fit time delay distances from the lensing data and the SNe reconstruction using the best fit $H_{0}$ from the combination. For only four points, a $\chi^{2}$ has limited significance but we can mention that $\chi^{2}=2.28,0.15,0.80,0.24$ for RXJ1131-1231, HE0435-1223, B1608+656, SDSS1206+4332, respectively, for the data plotted with respect to the GP distance relation. We note that with the systems ordered by increasing time delay distance, one does see a trend where the lensing data monotonically climbs higher above the relation predicted by the SNe sample with distance. For such a small lensing data sample it is difficult to tell whether this reflects a real systematic. (Note that Wong et al. (2019) shows another monotonic trend, in the derived $H_{0}$, with $D_{\Delta t}$, using six strong lens systems). It would be interesting as data sets get larger to study not only the mean $H_{0}$ derived, but whether any trends exist, which could potentially point to systematics with distance (or other physical characteristics) such as line of sight mass corrections or stellar dynamics scale effects.

To avoid assuming any value for $H_{0}$ to employ the SNe data, we can also consider ratios of time delay distances, which are independent of $H_{0}$. Furthermore, to explore the possibility of trends we plot these against each other. Rather than show all 15 plots of pair combinations of the 6 ratios, we show only 2 in Fig. 4: neighbors – $D_{\Delta t,1}/D_{\Delta t,2}$ vs $D_{\Delta t,3}/D_{\Delta t,4}$, and extremes – $D_{\Delta t,1}/D_{\Delta t,4}$ vs $D_{\Delta t,2}/D_{\Delta t,3}$.

We see that the data and reconstructed distances are consistent at the 68.3% confidence level. Most of the apparent trend with distance seems to be due to the nearest lens, RXJ1131-1231. This could be simply a statistical fluctuation but note that Birrer et al. (2016) found that the use of a different lens galaxy kinematic prior for this system could shift the value of $H_{0}$ from it down by 14%. In addition, Chen et al. (2019), using adaptive optics imaging only, reduced its value of $H_{0}$ from $78.2\pm 3.4$ km/s/Mpc to $77.0^{+4.0}_{-4.6}$ km/s/Mpc (though not a significant shift).

To explore the potential impact of an outlier, we repeat our model independent analysis using only three of the time delay distances at a time, and investigate whether the remaining system caused a shift in the final $H_{0}$ constraint. Figure 5 shows the results. Removal of RXJ1131-1231 does indeed have the greatest impact in alleviating tension with the Planck value of $H_{0}$, and it was the lens system most in tension with the SNe reconstruction. Note that as this is the lowest redshift system, one might expect that the SNe, which densely sample those distances, should provide an accurate result. Finally, we can mention that independent measurement of the time delays themselves shows a small, but interesting effect. The COSMOGRAIL team has been extremely open and helpful about releasing the lightcurve monitoring data. Hojjati et al. (2013) used a GP method in 2013 to measure the time delay from RXJ1131-1231 and found a value 1.4% larger than COSMOGRAIL’s (still within $1\sigma$). H0LICOW (Suyu, 2013) kindly put this value into their analysis pipeline and found this decreased the value of $H_{0}$ from this system by 1.4% (The magnitudes don’t need to be the same because changing $\Delta t$ also changes the lens modeling.) So at least at a minor level, the values could potentially come into greater agreement.

4. Conclusion and discussions

We carry out a model independent analysis (with no distance parametrizations or assumptions about any dark energy model) using Gaussian process regression to determine the Hubble constant $H_{0}$ by anchoring the Pantheon Type Ia supernovae with time delay distances from the four publicly released, robust H0LiCOW lenses. The analysis is done for both flat and non-flat universes.

Furthermore, we explore both the internal consistency of the time delay systems estimations of $H_{0}$ and the external consistency of the time delay distances with the GP expansion history derived from SNe. This includes new probability distribution functions for $H_{0}$ for each lens system in the model independent analysis. No statistically significant tensions are found, although there are some $\sim 1\sigma$ trends that could be checked with future data.

Our model independent results are $H_{0}=72.2\pm 2.1\,$km/s/Mpc for a flat Universe and $73.0^{+3.0}_{-2.8}$ km/s/Mpc allowing curvature. These are consistent with the time delay and time delay plus SNe results of Wong et al. (2019); Taubenberger et al. (2019); Collett et al. (2019), made within specific cosmological models or with polynomial fitting of the distance relation. Collett et al. (2019) find $H_{0}=74.2^{+3.0}_{-2.9}$ km/s/Mpc for a flat Universe and $H_{0}=75.7^{+4.5}_{-4.4}$ km/s/Mpc allowing for curvature. Despite not assuming a specific model, the uncertainties in our constraints are comparable to these, while reducing possible bias. They are also consistent with local distance measures of $H_{0}$, lying midway between Freedman et al. (2019) and Riess et al. (2019). The strong lensing time delay plus SNe model independent method looks quite promising as further data on strong lens time delay systems becomes available.

For future studies, current surveys like the Dark Energy Survey (DES) (Treu et al., 2018) and the Hyper SuprimeCam Survey (HSC) (More et al., 2017), and the upcoming surveys like the Large Synoptic Survey Telescope (LSST) (Oguri & Marshall, 2010) and Euclid and WFIRST satellites (Barnacka, 2018; Petrushevska et al., 2018) will bring us thousands of lensed quasars and over one hundred lensed SNe Ia, a part of which will have well-measured time delays (Liao et al., 2015). With high-quality ancillary observations, some dozens of systems will give us time delay distances at percent levels. Moreover, the angular diameter distances will be well-measured as well (Jee et al., 2016; Liao, 2019). In combination of these two kinds of distances, the Hubble constant could be constrained to sub-percent precision. However, at that stage, systematic errors should be important. Dedicated analysis should be applied to individual lenses such that the combinations are robust.

Supernovae data will continue to improve as well, playing an important role as a dense sampler of cosmic expansion history over a wide range of redshifts. In addition to strong lensing, further local distance measurements such as gravitational waves from mergers of binary stars as standard sirens could join this method and provide strong anchoring ability to model-independently determine the Hubble constant.

Acknowledgments

We thank Stefan Taubenberger for providing the results for comparison. KL and AS would like to thank Zong-Hong Zhu and Beijing Normal University for the hospitality where the discussions on this project initiated. KL was supported by the National Natural Science Foundation of China (NSFC) No. 11603015 and the Fundamental Research Funds for the Central Universities (WUT:2018IB012). AS would like to acknowledge the support of the Korea Institute for Advanced Study (KIAS) grant funded by the Korea government. EL is supported in part by the Energetic Cosmos Laboratory and by the U.S. Department of Energy, Office of Science, Office of High Energy Physics, under Award DE-SC-0007867 and contract no. DE-AC02-05CH11231.