Determining Model-Independent $H_{0}$ and Consistency Tests

The flat cosmological constant plus cold dark matter ($\Lambda$CDM) model is currently taken as the concordance cosmological scenario and explains a wide range of observations. However, there is significant tension in the $H_{0}$ value inferred within $\Lambda$CDM from the cosmic microwave background (CMB) observations $H_{0}=67.4\pm 0.5\rm{\ km/s/Mpc}$ (Planck, 2018) and that measured through the Cepheid distance ladder, $H_{0}=74.03\pm 1.42\rm{\ km/s/Mpc}$ (Riess et al., 2019). Other cosmological probes such as baryon acoustic oscillations together with primordial nucleosynthesis constraints agree with the CMB value of $H_{0}$ (Addison et al., 2018; Macaulay et al., 2019; Cuceu et al., 2019; Philcox et al., 2020), while other distance ladder techniques can lie in between (Freedman et al., 2020).

The discrepancy could arise due to unaccounted for systematic errors in observations or reveal new physics significantly different from $\Lambda$CDM.

Independent cosmological probes could provide new perspectives. Strong gravitational lensing by galaxies offers an independent method of determining $H_{0}$ through time-delay lens systems. A typical lensing system consists of a distant active galactic nucleus (AGN) lensed by a foreground elliptical galaxy, forming multiple images along with the arcs of the host galaxy. The images will be magnified and the light will arrive at the Earth delayed by various times. Since AGNs are variable, we can measure the time delays between any two images from their light curves.

From the time delay plus the lens potential measured by the high-resolution imaging and line-of-sight environment we can measure the “time-delay distance” $D_{\Delta t}$. This is a combination of three angular diameter distances, primarily depending on $H_{0}$ though also more weakly depending on the cosmological model, e.g. the matter density, dark energy properties, etc. In addition, with kinematic information on the lens galaxy (for example the velocity dispersion) the angular diameter distance to the lens $D_{\rm{d}}$ can be also determined independent of the external convergence from line of sight perturbers (Paraficz & Hjorth, 2009; Jee et al., 2015, 2016).

Either $D_{\Delta t}$ or $D_{\rm{d}}$, or both jointly, provides a one-step method of determining $H_{0}$, which is substantially independent of and complementary to the CMB, large scale structure, and distance ladder methods. Thus they give a much needed crosscheck. However, like the other cosmological probes, one has to assign a cosmological model when computing the lensing distances. The results may therefore differ for different models.

Rather than computing the distances within a model, one can instead measure the distance-redshift relation. To do so, one needs a probe that is both accurate and samples distance much more densely. Type Ia supernovae (SNe Ia) are superb mappers of the distance-redshift relation. However, they only provide relative distances because the SNe absolute magnitude (and Hubble constant) is unknown. The strengths of the two probes can combine together to remove each weakness: absolute distance measurements from time-delay lensing and relative distances from SNe Ia can anchor each other (for anchoring one type of distance with another, see for example Aubourg et al. (2015); Cuesta et al. (2015); Collett et al. (2019); Pandey, Raveri, Jain (2019)). When combining the two, the results under different cosmological models seem to be stable and consistent (Taubenberger et al., 2019).

One can make the cosmology-model independence (i.e. no form assumed for the expansion history $H(z)$) more explicit than simply checking under different models. In Liao et al. (2019) we applied Gaussian process (GP) regression to SNe Ia data to get a model-independent relative distance-redshift relation, i.e. the shape of distance-redshift function. Anchoring this together with $D_{\Delta t}$ from 4 H0LiCOW lenses resulted in $H_{0}=72.2\pm 2.1\rm{\ km/s/Mpc}$ in a flat universe and $H_{0}=73.0_{-3.0}^{+2.8}\rm{\ km/s/Mpc}$ considering the cosmic curvature density in the range $\Omega_{\rm{k}}=[-0.2,0.2]$.

Currently, the H0LiCOW program has reached its first milestone (Wong et al., 2019). The full dataset consists of 6 lenses, five of which were analyzed blindly, and four of which have both $D_{\Delta t}$ and $D_{\rm{d}}$ measurements. In a flat $\Lambda$CDM cosmology, they found $H_{0}=73.3_{-1.8}^{+1.7}\rm{\ km/s/Mpc}$, consistent with the local Cepheid distance ladder measurement but in $3.1\sigma$ tension with $Planck$ and other CMB and large scale structure measurements. Note the H0LiCOW results vary with the assumed models (Wong et al., 2019), although $H_{0}$ tends to increase with the usual generalized cosmologies. Furthermore, the previously identified trend (Wong et al., 2019; Liao et al., 2019) in derived $H_{0}$ value, or time delay distance excess, with redshift persists. Therefore, it is timely and useful to update our model-independent results.

We briefly introduce the H0LiCOW program and the lensing dataset in Section 2, and present our methodology and updated results in Section 3. In Section 4 we carry out some consistency checks between lensing data, and between lensing and supernova data. We explore the previously identified trend of $H_{0}$ with time-delay distance or lens redshift in Section 5, and present a possible physical explanation based on modified gravity. We summarize and discuss the results and next steps in Section 6.

Strong lensing (SL) by elliptical galaxies is a powerful tool to study both astrophysics and cosmology (Treu, 2010). Lenses with time-delay measurements were proposed to measure $H_{0}$ (Refsdal, 1964; Treu & Marshall, 2016). Specifically, the time delay between any two images is determined by

where the time-delay distance

is a combination of three angular diameter distances $D_{\rm{d}},D_{\rm{s}},D_{\rm{ds}}$, where the subscripts d and s denote the deflector (lens) and the source, respectively. We use units where $c=G=1$ throughout the article. The Fermat potential difference $\Delta\phi$ between the two images is a function of lens mass profile parameters $\boldsymbol{\xi}_{\rm{lens}}$, determined by high-resolution imaging of the host arcs. Note that all other mass along the line of sight could also contribute to the lens potential, causing additional (de)focusing of the light rays and affecting the observed time delays. Considering the effects of the perturber masses are small, they can be approximated by an external convergence $\kappa_{\rm{ext}}$. Then the inferred $D_{\Delta t}$ will be scaled by 1-$\kappa_{\rm{ext}}$. Independent observations such as galaxy counts could break the degenaracy (Rusu et al., 2017). We take the time delay distances as given by H0LiCOW, and refer readers to the systematics treatment in Millon et al. (2019).

Additional information on the lens galaxy such as the light profile $\boldsymbol{\xi}_{\rm{light}}$, the projected stellar velocity dispersion $\sigma^{\rm{P}}$, and the anisotropy distribution of the stellar orbits, parametrized by $\beta_{\rm{ani}}$, can yield the angular diameter distance to the lens (Birrer et al., 2016, 2019):

which correlates with $D_{\Delta t}$. The function $J$ captures all the model components computed from angles measured on the sky (from imaging) and the stellar orbital anisotropy distribution (from spectroscopy). We refer to Section 4.6 of Birrer et al. (2019) for the detailed modeling related to $J$.

The state-of-the-art lensing collaboration H0LiCOW aims at measuring $H_{0}$ with $1\%$ precision using time-delay lenses (Suyu et al., 2017). They take advantage of substantial data consisting of time-delay measurements from the COSMOGRAIL program¹¹1http://www.cosmograil.org and radio-wavelength monitoring, deep HST and ground-based adaptive optics imaging, spectroscopy of the lens galaxy, and deep wide-field spectroscopy and imaging.

In the recent milestone paper (Wong et al., 2019), they gave the latest constraints on $H_{0}$ under different cosmological models with a combined sample of six lenses that span a range of lens and source redshifts, as well as various image configurations (double, cross, fold, and cusp). All lenses except the earliest, B1608+656, were analyzed blindly with respect to the cosmological parameters. Four lenses (RXJ1131-1231, PG 1115+080, B1608+656, SDSS 1206+4332) have both $D_{\Delta t}$ and $D_{\rm{d}}$ measurements. Note that $D_{\Delta t}$ and $D_{\rm{d}}$ measurements for a system are correlated except for B1608+656 whose distance measurements are independent. The distance posterior distributions are released in the form of MCMC chains or skewed log-normal function fits on the H0LiCOW website²²2http://www.h0licow.org. For the case of MCMC chains, we will get the likelihood functions by smoothing the discrete points. We summarize the redshifts and measured distances in Table 1 ordered by the lens redshift. For more detailed information on these lenses, see Wong et al. (2019) and the references therein.

To combine the SNe data with lensing data, we generate samples of unanchored luminosity distance $H_{0}D^{\rm{L}}$ from the posterior of the Pantheon compilation from  Scolnic et al. (2018), calculated with a GP (See, e.g., Koo et al. (2020) for a test of cosmology model independence.) This paper follows the analysis in Liao et al. (2019), which is based on the gphist (Kirkby & Keeley, 2017) code first presented in Joudaki et al. (2018).

Regression using a GP works by generating a set of functions from an infinite dimensional function space characterized by a covariance function. This covariance function is parametrized by a squared-exponential kernel

where $s_{i}=\log(1+z_{i})/\log(1+z_{\rm max})$ and $z_{\rm max}=2.26$ is the maximum redshift of the supernova sample. This has two hyper-parameters that are marginalized over, $\sigma_{f}$ and $\ell$, with $\sigma_{f}$ determining the amplitude of the random fluctuations and $\ell$ determining the coherence length of the fluctuation, equivalently $1/\ell$ is proportional to the number of fluctuations in the range. The priors on these hyper-parameters are scale-invariant and we directly integrate over this space since the dimensionality is small.

The GP function $\gamma(z)=\ln(H^{\rm fid}(z)/H(z))$ involves the expansion history $H(z)$ (which can be integrated to give distances), and $H^{\rm fid}(z)$ is the best-fit $\Lambda$CDM expansion history from the Pantheon dataset and works as the mean function of the GP regression. The GP prior functions are then trained on the Pantheon likelihood, which constrains only the shape of the expansion history, not the absolute scale, so unanchored luminosity distances $H_{0}D^{\rm{L}}(z)$ are the quantities most directly constrained by the Pantheon SNe dataset.

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

1. 1.

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

2. 2.

   Evaluate the values of each of the 1000 $H_{0}D^{\rm{A}}$ curves at the lens and source redshifts of the six (SL) systems to calculate 1000 values of $H_{0}D_{\Delta t}$ for each system using

   $$H_{0}D_{\Delta t}=(1+z_{\mathrm{d}})(H_{0}D_{\mathrm{d}})(H_{0}D_{\mathrm{s}})/(H_{0}D_{\mathrm{ds}})\,;$$

   (5)

3. 3.

   Compute the likelihood, for each of the 1000 realizations, from the H0LiCOW’s $D_{\Delta t}$ combined with $D_{\rm{d}}$ data (if the $D_{\rm{d}}$ measurements are available) for each lens system for many values of $H_{0}$;

4. 4.

   Multiply the six 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 $D_{\rm{ds}}$ between the lens and the source from $H_{0}D^{\rm{A}}$ in step 2, we use the standard distance relation (Weinberg, 1972)

where $\Omega_{k}$ is the dimensionless curvature density. For a spatially flat universe, one simply has $D_{\rm{ds}}=D_{\rm{s}}-[(1+z_{\rm{d}})/(1+z_{\rm{s}})]\,D_{\rm{d}}$.

In step 3, we also evaluate a case with the SL $D_{\rm{d}}$ data alone. In such a case, the anchoring is direct between absolute and relative distances and one does not need to consider the assumption of cosmic curvature.

Figure 1 shows the joint posteriors and the numerical results are summarized in Table 2. First of all, we consider the full data set consisting of 6 lenses, 4 of which have both $D_{\Delta t}$ and $D_{d}$ measurements. The constraint is $H_{0}=72.8_{-1.7}^{+1.6}\rm{\ km/s/Mpc}$ assuming the universe is flat. It is more stringent than our previous result based on 4 lenses with only $D_{\Delta t}$ data, $H_{0}=72.2\pm 2.1\rm{\ km/s/Mpc}$. This model-independent result has comparable value and error bars as that under the flat $\Lambda$CDM model by H0LiCOW, $H_{0}=73.3_{-1.8}^{+1.7}\rm{\ km/s/Mpc}$, supporting the tension with the BAO and CMB results. (Note that the looser constraints expected from removing the $\Lambda$CDM assumption are offset by including supernova distances.)

In addition, we test the contribution of the $D_{\rm{d}}$ data by taking $D_{\Delta t}$ and $D_{\rm{d}}$ separately. The constraint from 4 $D_{\rm{d}}$ data alone is relatively weaker, $H_{0}=81.0_{-6.9}^{+7.1}\rm{\ km/s/Mpc}$, and the 6 $D_{\Delta t}$ data show almost the same constraint power, $H_{0}=72.8_{-1.8}^{+1.7}\rm{\ km/s/Mpc}$, as the full set. Nevertheless, the result from $D_{\rm{d}}$ data is free from assumptions concerning spatial curvature (although note that other lens parameters, such as stellar anisotropy, enter.)

We also consider the non-flat case for the full dataset for completeness. We set the uniform prior of the curvature parameter to be $\Omega_{\rm{k}}=[-0.5,0.5]$ as in Wong et al. (2019). The constraints are $H_{0}=77.3_{-3.0}^{+2.2}\rm{\ km/s/Mpc}$ and $\Omega_{\rm{k}}=0.33_{-0.19}^{+0.12}$. Like the $\Lambda$CDM-assumed case where $H_{0}=74.4_{-2.3}^{+2.1}\rm{\ km/s/Mpc}$, $\Omega_{\rm{k}}=0.26_{-0.25}^{+0.17}$ (Wong et al., 2019), our non-flat result shows a larger $H_{0}$ (possibly the $D_{d}$ data, which give a higher $H_{0}$, have more influence in this case) and slightly favors an open universe. However, the change of $H_{0}$ relative to the flat case in this work is more distinct. It is worth mentioning that the constraint on $\Omega_{\rm{k}}$ in our method is model-independent as well.

Furthermore, to understand the relative contributions, we constrain $H_{0}$ with each lens in the model-independent manner. We use both $D_{\Delta t}$ and $D_{\rm{d}}$ data (if it has). Figure 2 shows the individual posteriors. The numerical results are shown in Table 3, in order of increasing lens redshift. As one can see, a trend may exist: the $H_{0}$ values roughly decrease with the lens redshift or the distance. Our results further confirm the trend noticed in Wong et al. (2019); Liao et al. (2019) and shown in Fig. 5 in Millon et al. (2019) based on $\Lambda$CDM. We explore this further in Section 5.

In this section, we check to what extent the SL data are internally consistent and to what extent they are consistent with the SNe data. This is important to confirm that the results from the combination of lens systems, and from combination of SL and SNe, are robust. As the available SL distances become more numerous and more precise, one can better assess the consistency of these distances with the SNe distances. We begin by comparing the distance posteriors of the SL data, with the values of those quantities predicted by the SNe data.

In Fig. 3, we plot the measured posteriors of the individual systems’ time-delay distances at the 68% and 95% confidence level, as well as the SNe posterior predictions of those systems’ time-delay distances using a GP regression on the Pantheon data. (The uncertainties on the angular distances are currently too large to add useful information to this comparison.) Since the SNe cannot constrain an absolute scale, the SNe posterior predictions for the SL distances are anchored in this figure by our combined $H_{0}$ measurement $H_{0}=72.8^{+1.6}_{-1.7}\rm{\ km/s/Mpc}$, giving the extended green contours, reflecting both the uncertainty in $H_{0}$ and the GP-based unanchored distances ($H_{0}D_{A}$). On the whole, the measurements of the SL distances are consistent with the SNe posterior predictions, in both the 1D and 2D joint posteriors, although since we have ordered the systems by lens redshift one can notice a tendency for the SL posterior to drift rightward (to higher distance and hence lower $H_{0}$) across the SNe posterior with increasing redshift. This will be explored further in Sec. 5.

By forming ratios of time-delay distances, we can cancel out the dependence on $H_{0}$, forming relative distances that the SNe are particularly suited to. In Fig. 4, we plot the measured SL posteriors of the ratios of time-delay distances for certain combinations of the systems, as well as the posterior predictions for those ratios from the SNe data. Indeed, the SNe posterior predictions are nearly pointlike. Rather than showing hundreds of 2D joint combinations of the 15 possible time-delay ratios, we select two cases: the left panel shows the ratios for the two pairs of systems nearly at the same lens redshift (see Table 1), and the right panel for two pairs of systems with the most extreme differences in lens redshift. The pairs nearly at the same redshift (one pair with both $z_{\rm{d}}<0.4$ and one pair with both $z_{\rm{d}}>0.4$) are highly consistent with the very precise SNe predictions. Those pairs at very different redshifts (thus one system has $z_{\rm{d}}<0.4$ and one system has $z_{\rm{d}}>0.4$), are shifted to the boundary of the 68% confidence region, in both 1D and 2D posteriors. This is not statistically significant, but does seem to continue a trend.

To investigate this further, in Fig. 5, we plot the measured SL 1D posteriors of every combination of the ratios of the time-delay distances of sources with lens redshifts both at $z_{\rm{d}}>0.4$ or both at $z_{\rm{d}}<0.4$ (left panel), or one at $z_{\rm{d}}>0.4$ and one at $z_{\rm{d}}<0.4$ (right panel). These 1D posterior ratios are given relative to the predictions of those same ratios from the SNe data, i.e. $(D_{\Delta t,i}/D_{\Delta t,j})/(D_{\Delta t,i}/D_{\Delta t,j})_{\rm SN}$, so consistency gives the value of 1. This can also be viewed equivalently as $(D_{\Delta t,i}/(D_{\Delta t,i})_{\rm SN})/(D_{\Delta t,j}/(D_{\Delta t,j})_{\rm SN}))$, testing that SL cosmology is consistent with SNe cosmology, irrespective of $H_{0}$.

Interestingly, when both lens systems in the ratio have $z_{\rm{d}}>0.4$, or both have $z_{\rm{d}}<0.4$, then the posteriors of the 7 possible combinations all peak very close to unity, showing consistency (see left panel). However, when the lenses in the ratio lie on opposite sides of $z_{\rm{d}}=0.4$, then although the posteriors are still consistent with the value 1, the peaks tend to be $\sim 10\%$ higher. Any one posterior could fluctuate above the consistency value of 1, but we note that all 8 possible combinations all lie high. This is not to say that one should multiply 8 $1\sigma$ deviations (Wong et al. (2019) find this trend to be somewhat under $2\sigma$); we merely consider it odd enough to make it worthwhile looking for physics (remember, this is independent of the value of $H_{0}$) that could account for this. This is what we explore in the next section. Alternately, it could be due to some unidentified observational systematic, and one should see if this putative trend persists with new lens systems and new data.

In going from the four lens systems of Birrer et al. (2019) to the seven lens systems of Wong et al. (2019) and Shajib et al. (2019) (the posterior chains for the last, DES J0408-5354, have not been publicly released yet but we include in this section their quoted $H_{0}$ value from it), the apparent trend of the time-delay distance with respect to the distance relation predicted by the supernova luminosity distance (Liao et al., 2019), or derived Hubble constant (Wong et al., 2019), as a function of time-delay distance or lens redshift has remained. In particular, see Fig. 5 of Millon et al. (2019). The statistical sample is still small, so perhaps this can be just an odd fluctuation despite the extra systems not reducing the trend.

In this section we briefly explore the conjecture that the trend is physical, and that it could be related not to the background expansion history (distances, so SN and BAO are unaffected) but rather the behavior of gravity on light deflection (lensing) evolving with redshift. This is commonly called $G_{\rm light}(z)$ and arises in many modified gravity theories.

The light deflection depends on the sum of the time-time and space-space metric potentials, $\Phi+\Psi$, and is related to the density contrast $\delta\rho/\rho$ by

Could the trend be reflecting $G_{\rm light}(z)$? Note that $G_{\rm light}$ does not affect supernova distances, so those would reflect the actual background expansion history.

Here we simply give a rough analysis, ignoring some subtleties we mention later. The measured lensing time delays depend on

where $\Delta\phi$ is the difference in Fermat potentials. The Fermat potential difference

where $\theta_{i}-\beta$ is the angular difference between the image location and unlensed source location, and $\psi$ is a projected potential. The angular deflection with modified gravity light propagation becomes

and the projected potential is related to the convergence from mass along the line of sight by

Thus the gravitational effects on light propagation give

Since the Hubble constant estimated from a lens system comes from $1/D_{\Delta t}$ then for given measured lens system characteristics we have

That is, if $G_{\rm light}$ is increasing with $a$ then the derived value of $H_{0}$ will increase for lower redshift lens systems.

We can assume this is the cause of the observed trend in derived $H_{0}$ from the lens systems and derive what function $G_{\rm light}(a)$ is needed. At high redshift we expect gravity to restore to general relativity (e.g. to preserve the successes of the cosmic microwave background and primordial nucleosynthesis) so we take values of $H_{0}$ derived at high redshift to be the true values. Figure 6 shows the $H_{0}$ values given in Millon et al. (2019), within $\Lambda$CDM, for the seven lens systems and an illustrative power law in scale factor,

(Note that for many modified gravity models $G_{\rm light}$ actually levels out to a constant de Sitter value not far into the future.)

The rise in measured $H_{0}$ values comes in quite steeply with scale factor, around $a\gtrsim 0.7$ (or $z\lesssim 0.4$). A change so rapid to distances themselves is quite difficult to achieve, since distances are integrals over the expansion rate (and double integrals over dark energy equation of state). And certainly no such dramatic change in seen in supernova distances. However, by changing the gravitational strength affecting light deflection, $G_{\rm light}(a)$, this is less difficult. Indeed, such a rapid change has been shown to occur for some actual modified gravity theories – see for example Fig. 5 (right panel, thin curves) of Linder (2017). Note the numerical solution there shows that indeed the effect on $G_{\rm light}(a)$ can first become significant at low redshifts. While that particular theory (uncoupled Galileon) is ruled out, it does give a proof of principle that modified gravity can act in such a manner.

One might also speculate that the effect on light deflection could show up in weak lensing measurements. For the convergence or shear power spectrum, the relativistic Poisson equation (7) shows that $G_{\rm light}\kappa\sim(\Phi+\Psi)$ and so the measured shear power will be proportional to $G_{\rm light}^{-2}$. There is an extra element in that the growth of structure also depends on $G_{\rm matter}$, but we focus here on $G_{\rm light}$. The value of the mass fluctuation amplitude $\sigma_{8}$, or $S_{8}=\sigma_{8}(\Omega_{m}/0.3)^{1/2}$, derived from the shear power spectrum will thus be proportional to $1/G_{\rm light}(a)$. This was discussed in Daniel & Linder (2010) – “For higher [$G_{\rm light}$], lower values of $\sigma_{8}$ will produce comparable lensing potentials. Larger [$G_{\rm light}$] does not cause $\sigma_{8}$ to decrease per se, rather it brings lower values of $\sigma_{8}$ into better agreement with the data”. Since we take $G_{\rm light}$ to be strengthening at lower redshift for the SL case, this means that the value of $\sigma_{8}$ or $S_{8}$ derived from low redshift surveys should be less than from high redshift surveys (or Planck CMB). This trend does seem consistent with weak lensing survey data (but again, evolution of $G_{\rm matter}$ can overturn this).

Regarding the subtleties we mentioned at the beginning of this section, note that the light deflection occurs all along the path and not just at the lens, but as with general relativistic light deflection one can treat the deflection as occurring at the lens, in a single screen approximation. Right at the lens we might expect the modified gravity to be screened, but further out from the lens the screening vanishes and the dominant part of the path integral is roughly at the lens redshift. Thus we take $\langle G_{\rm light}(a)\Delta\phi\rangle\approx G_{\rm light}(a_{\rm lens})\Delta\phi$.

Finally, it is not clear how $G_{\rm light}(a)$ could help with the Cepheid measurement of higher $H_{0}$, though the value from the tip of the red giant branch technique is more consistent. This whole section is simply speculation, but if the trend in $H_{0}$ measurements with distance persists, we might consider it a question of gravity.

Based on the method we previously proposed, we give a cosmology model-independent determination of $H_{0}$ with the updated H0LiCOW dataset consisting of six lenses. The absolute lensing distances ($D_{\rm{{\Delta t}}}$ and $D_{\rm{d}}$) are used to anchor the Pantheon SNe samples that give the shape of the distance-redshift relation through GP regression. The results are $H_{0}=72.8_{-1.7}^{+1.6}\rm{\ km/s/Mpc}$ for a flat universe and $H_{0}=77.3_{-3.0}^{+2.2}\rm{\ km/s/Mpc}$ for a non-flat universe. These values are consistent with the results assuming a $\Lambda$CDM model, and have comparable uncertainties, though they have the advantage of being cosmology model independent (and include SN). With current data, $D_{\rm d}$ measurements do not play a significant role, though they have the property of being relatable to SN independent of spatial curvature.

We perform several consistency tests of the data, and between the different probes. All show consistency, though an odd systematic trend persists in the value of the derived $H_{0}$ with lens redshift. We illustrate this trend through several methods. In particular, one could interpret it as a transition at $z\approx 0.4$. Irrespective of the value of $H_{0}$, the distances from SL systems lying all below or all above $z\approx 0.4$ are highly consistent with the SNe cosmology (this holds for all 7 such combinations of systems), but comparison of lensing systems on either side of $z\approx 0.4$ all show an offset (admittedly individually statistically minor) from the SN cosmology – this holds, in the same direction, for all 8 such combinations of systems.

This could be a statistical fluke (though it has persisted since the first analysis with fewer systems) or some observational systematic. We speculate about one possible explanation based on physics beyond the standard model, showing how a modification in the effect of gravity on light propagation, $G_{\rm light}(a)$, could account for this. Moreover, models in the literature show that the magnitude and redshift dependence of such an effect is possible, and this could also affect the perceived value of $\sigma_{8}$, possibly bearing on that tension as well. Given the low statistical significance with current data, we merely suggest keeping an eye on whether further, or improved, data continue to support such a physics explanation.

Fortunately time-delay strong lensing is a burgeoning field with the onset of cosmic surveys, and more well-measured and well-analyzed lenses are not far off, with imaging surveys such as DES, ZTF, LSST, and Euclid. Monitoring campaigns, adaptive optics, and spectroscopic followup (including multiobject instruments such as DESI) all play important roles as well. More lens systems at all redshifts – near $z\approx 0.4$, below, and well above will test whether a standard cosmology matches both strong lenses and supernovae. Continued detailed systematics studies of all distance indicators – and all light deflection probes – will be essential for confirming any result.

KL was supported by the National Natural Science Foundation of China (NSFC) No. 11973034. AS would like to acknowledge the support of the Korea Institute for Advanced Study (KIAS) grant funded by the Korea government and Kobayashi-Maskawa Institute (KMI) and Nagoya University for their hospitality in the final stages of this project. 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.