H0LiCOW XIII. A 2.4% measurement of $H_{0}$ from lensed quasars: $5.3\sigma$ tension between early and late-Universe probes

TL;DR

This paper measures the Hubble constant using gravitational lensing time delays, achieving 2.4% precision, and finds a significant tension with early-Universe CMB measurements, highlighting potential new physics or systematic issues.

Contribution

It provides a highly precise, independent measurement of H0 from gravitational lensing, confirming the tension with CMB-based estimates and exploring extensions to standard cosmology.

Findings

H0 = 73.3 km/s/Mpc with 2.4% uncertainty

5.3σ tension with Planck CMB measurements

Time-delay cosmography enhances other cosmological constraints

Abstract

We present a measurement of the Hubble constant ($H_{0}$) and other cosmological parameters from a joint analysis of six gravitationally lensed quasars with measured time delays. All lenses except the first are analyzed blindly with respect to the cosmological parameters. In a flat $\Lambda$CDM cosmology, we find $H_{0} = 73.3_{-1.8}^{+1.7}$, a 2.4% precision measurement, in agreement with local measurements of $H_{0}$ from type Ia supernovae calibrated by the distance ladder, but in $3.1\sigma$ tension with $Planck$ observations of the cosmic microwave background (CMB). This method is completely independent of both the supernovae and CMB analyses. A combination of time-delay cosmography and the distance ladder results is in $5.3\sigma$ tension with $Planck$ CMB determinations of $H_{0}$ in flat $\Lambda$CDM. We compute Bayes factors to verify that all lenses give statistically consistent results, showing that we are not underestimating our uncertainties and are able to control our systematics. We explore extensions to flat $\Lambda$CDM using constraints from time-delay cosmography alone, as well as combinations with other cosmological probes, including CMB observations from $Planck$, baryon acoustic oscillations, and type Ia supernovae. Time-delay cosmography improves the precision of the other probes, demonstrating the strong complementarity. Allowing for spatial curvature does not resolve the tension with $Planck$. Using the distance constraints from time-delay cosmography to anchor the type Ia supernova distance scale, we reduce the sensitivity of our $H_0$ inference to cosmological model assumptions. For six different cosmological models, our combined inference on $H_{0}$ ranges from $\sim73$-$78~\mathrm{km~s^{-1}~Mpc^{-1}}$, which is consistent with the local distance ladder constraints.

Tables

Table 1. Table 1: Lenses in the H0LiCOW sample used in this paper.

Lens name	$\alpha$ (J2000)	$\delta$ (J2000)	$z_d$	$z_s$	HST / AO data
B1608$+$656a	16:09:13.96	$+$65:32:29.0	0.6304a	1.394b	HST
RXJ1131$-$1231c	11:31:51.6	$-$12:31:57.0	0.295c	0.654d	HST + AO
HE 0435$-$1223e	04:38:14.9	$-$12:17:14.4	0.4546f,g	1.693h	HST + AO
SDSS 1206$+$4332i	12:06:29.65	$+$43:32:17.6	0.745j	1.789i	HST
WFI2033$-$4723k	20:33:41.9	$-$47:23:43.4	0.6575l	1.662h	HST
PG 1115$+$080m	11:18:16.899	$+$7:45:58.502	0.311n	1.722m	HST + AO

Table 2. Table 2: D Δ ​ t subscript 𝐷 Δ 𝑡 D_{\Delta t} and D d subscript 𝐷 d D_{\rm d} constraints for H0LiCOW lenses.

Lens name	$D_{\Delta t}$ (Mpc)	$D_d$ (Mpc)	Blind analysis	References
B1608$+$656	${5156}_{-236}^{+296}$	${1228}_{-151}^{+177}$	no	Suyu et al. (2010); Jee et al. (2019)
RXJ1131$-$1231	${2096}_{-83}^{+98}$	${804}_{-112}^{+141}$	yesa	Suyu et al. (2014); Chen et al. (2019)
HE 0435$-$1223	${2707}_{-168}^{+183}$	—	yes	Wong et al. (2017); Chen et al. (2019)
SDSS 1206$+$4332	${5769}_{-471}^{+589}$	${1805}_{-398}^{+555}$	yes	Birrer et al. (2019)b
WFI2033$-$4723	${4784}_{-248}^{+399}$	—	yes	Rusu et al. (2019)
PG 1115$+$080	${1470}_{-127}^{+137}$	${697}_{-144}^{+186}$	yes	Chen et al. (2019)

Table 3. Table 3: Approximate D Δ ​ t subscript 𝐷 Δ 𝑡 D_{\Delta t} error budget for H0LiCOW lenses.

Source of uncertainty	B1608$+$656	RXJ1131$-$1231	HE 0435$-$1223	SDSS 1206$+$4332	WFI2033$-$4723	PG 1115$+$080
Time delays	1.7%	1.6%	5.3%	2.3%	2.9%	6.4%
Line-of-sight contribution	6.4%	3.3%	2.8%	2.9%	5.7%	2.7%
Lens model and other sources	3.0%	2.2%	2.5%	8.4%	2.2%	5.7%
Total	5.1%	4.3%	6.5%	9.1%	6.7%	9.0%

Table 4. Table 4: Bayes factor for all pairs of lensed systems (top) and of every individual system relative to the five remaining systems (bottom).

Pairwise Bayes factors
	B1608$+$656	RXJ1131$-$1231	HE 0435$-$1223	SDSS 1206$+$4332	WFI2033$-$4723	PG 1115$+$080
B1608$+$656	—	3.4	10.6	9.1	11.2	3.3
RXJ1131$-$1231		—	5.3	3.1	4.9	6.9
HE 0435$-$1223			—	7.8	9.3	3.9
SDSS 1206$+$4332				—	7.9	2.7
WFI2033$-$4723					—	3.7
PG 1115$+$080						—
Individual lens Bayes factors vs. rest of sample
	B1608$+$656	RXJ1131$-$1231	HE 0435$-$1223	SDSS 1206$+$4332	WFI2033$-$4723	PG 1115$+$080
Bayes Factor	11.1	3.8	11.0	7.5	12.6	4.5

Table 5. Table 5: Description of the cosmological models considered in this work.

Model name	Description	Priors
Time-delay cosmography only
U$\Lambda$CDM	Flat $\Lambda$CDM	$H_0$ uniform in [0, 150] km s-1 Mpc-1
		${\Omega }_m=1-{\Omega }_{\Lambda }$
		${\Omega }_m$ uniform in [0.05, 0.5]
Uo$\Lambda$CDM	Open $\Lambda$CDM	$H_0$ uniform in [0, 150] km s-1 Mpc-1
		${\Omega }_m$ uniform in [0.05, 5]
		${\Omega }_k$ uniform in [$-$0.5, 0.5]
		${\Omega }_{\Lambda }=1-{\Omega }_m-{\Omega }_k>0$
U$w$CDM	Flat $w$CDM	$H_0$ uniform in [0,150] km s-1 Mpc-1
		${\Omega }_m$ uniform in [0.05, 5]
		$w$ uniform in [$-$2.5, 0.5]
		${\Omega }_{\mathrm{DE}}=1-{\Omega }_m$
U$w_0{w}_a$CDM	Flat $w_0{w}_a$CDM	$H_0$ uniform in [0,150] km s-1 Mpc-1
		${\Omega }_m$ uniform in [0.05, 5]
		$w_0$ uniform in [$-$2.5, 0.5]
		$w_a$ uniform in [$-$2, 2]
		${\Omega }_{\mathrm{DE}}=1-{\Omega }_m$
Time-delay cosmography combined with other probes
$\Lambda$CDM	Flat $\Lambda$CDM	${\Omega }_m=1-{\Omega }_{\Lambda }$
		JLA/Pantheon for {$H_0$, ${\Omega }_{\Lambda }$}
o$\Lambda$CDM	Open $\Lambda$CDM	${\Omega }_m=1-{\Omega }_{\Lambda }-{\Omega }_k>0$
		Planck (Section 5.3.1) or JLA/Pantheon (Section 5.4) for {$H_0$, ${\Omega }_{\Lambda }$, ${\Omega }_m$}
$w$CDM	Flat $w$CDM	${\Omega }_m=1-{\Omega }_{\mathrm{DE}}$
		Planck (Section 5.3.2) or JLA/Pantheon (Section 5.4) for {$H_0$, ${\Omega }_{\mathrm{DE}}$, $w$}
$N_{\mathrm{eff}}\Lambda$CDM	Flat $\Lambda$CDM	Planck for {$H_0$, ${\Omega }_{\Lambda }$, $N_{\mathrm{eff}}$}
	Variable $N_{\mathrm{eff}}$
$m_{\nu }\Lambda$CDM	Flat $\Lambda$CDM	Planck for {$H_0$, ${\Omega }_{\Lambda }$, $\sum m_{\nu }$}
	Variable $\sum m_{\nu }$
$N_{\mathrm{eff}}{m}_{\nu }\Lambda$CDM	Flat $\Lambda$CDM	Planck for {$H_0$, ${\Omega }_{\Lambda }$, $N_{\mathrm{eff}}$, $\sum m_{\nu }$ }
	Variable $N_{\mathrm{eff}}$ and $\sum m_{\nu }$
$w_0{w}_a$CDM	Flat $w_0$$w_a$CDM	Planck (Section 5.3.4) or JLA/Pantheon (Section 5.4) for {$H_0$, $w_0$, $w_a$}
o$w$CDM	Open $w$CDM	${\Omega }_m=1-{\Omega }_{\mathrm{DE}}-{\Omega }_k>0$
		JLA/Pantheon for {$H_0$, ${\Omega }_k$, ${\Omega }_{\mathrm{DE}}$, $w$}
o$w_0{w}_a$CDM	Open $w_0$$w_a$CDM	${\Omega }_m=1-{\Omega }_{\mathrm{DE}}-{\Omega }_k>0$
		JLA/Pantheon for {$H_0$, ${\Omega }_k$, ${\Omega }_{\mathrm{DE}}$, $w_0$, $w_a$}

Table 6. Table 6: Cosmological parameters for various cosmologies from time-delay cosmography only.

Model	$H_0$ (km s-1 Mpc-1)	${\Omega }_m$	${\Omega }_{\Lambda }$ or ${\Omega }_{\mathrm{DE}}$	${\Omega }_k$	$w$ or $w_0$	$w_a$
U$\Lambda$CDM	${73.3}_{-1.8}^{+1.7}$	${0.30}_{-0.13}^{+0.13}$	${0.70}_{-0.13}^{+0.13}$	$\equiv$ 0	$\equiv$ $-1$	$\equiv$ 0
Uo$\Lambda$CDM	${74.4}_{-2.3}^{+2.1}$	${0.24}_{-0.13}^{+0.16}$	${0.51}_{-0.18}^{+0.21}$	${0.26}_{-0.25}^{+0.17}$	$\equiv$ $-1$	$\equiv$ 0
U$w$CDM	${81.6}_{-5.3}^{+4.9}$	${0.31}_{-0.10}^{+0.11}$	${0.69}_{-0.11}^{+0.10}$	$\equiv$ 0	$-{1.90}_{-0.41}^{+0.56}$	$\equiv$ 0
U$w_0{w}_a$CDM	${81.3}_{-5.4}^{+5.1}$	${0.31}_{-0.11}^{+0.11}$	${0.69}_{-0.11}^{+0.11}$	$\equiv$ 0	$-{1.86}_{-0.45}^{+0.63}$	$-{0.05}_{-1.37}^{+1.45}$

Table 7. Table 7: Cosmological parameters for various cosmologies from time-delay cosmography combined with other probes. The o Λ Λ \Lambda CDM and w 𝑤 w CDM constraints can be compared with those from time-delay cosmography only (Table  6 ).

Model	H0LiCOW +	$H_0$ (km s-1 Mpc-1)	${\Omega }_m$	${\Omega }_{\Lambda }$ or ${\Omega }_{\mathrm{DE}}$	${\Omega }_k$	$w$ or $w_0$	$w_a$	$N_{\mathrm{eff}}$	$\sum m_{\nu }$ (eV)
o$\Lambda$CDM	Planck	${67.0}_{-1.4}^{+1.3}$	${0.32}_{-0.01}^{+0.01}$	${0.68}_{-0.01}^{+0.01}$	$-{0.002}_{-0.003}^{+0.003}$	$\equiv$ $-1$	$\equiv$ 0	$\equiv$ 3.046	$\equiv$ 0.06
o$\Lambda$CDM	Planck + CMBL	${69.2}_{-1.8}^{+1.7}$	${0.30}_{-0.01}^{+0.02}$	${0.70}_{-0.01}^{+0.01}$	${0.003}_{-0.004}^{+0.003}$	$\equiv$ $-1$	$\equiv$ 0	$\equiv$ 3.046	$\equiv$ 0.06
o$\Lambda$CDM	Planck + BAO	${68.4}_{-0.7}^{+0.7}$	${0.30}_{-0.01}^{+0.01}$	${0.69}_{-0.01}^{+0.01}$	${0.002}_{-0.002}^{+0.002}$	$\equiv$ $-1$	$\equiv$ 0	$\equiv$ 3.046	$\equiv$ 0.06
$w$CDM	Planck	${78.7}_{-2.7}^{+2.7}$	${0.23}_{-0.01}^{+0.02}$	${0.77}_{-0.02}^{+0.01}$	$\equiv$ 0	$-{1.36}_{-0.09}^{+0.09}$	$\equiv$ 0	$\equiv$ 3.046	$\equiv$ 0.06
$w$CDM	Planck + BAO	${70.3}_{-1.5}^{+1.6}$	${0.29}_{-0.01}^{+0.01}$	${0.71}_{-0.01}^{+0.01}$	$\equiv$ 0	$-{1.10}_{-0.06}^{+0.06}$	$\equiv$ 0	$\equiv$ 3.046	$\equiv$ 0.06
$N_{\mathrm{eff}}\Lambda$CDM	Planck	${68.8}_{-1.3}^{+1.2}$	${0.31}_{-0.01}^{+0.01}$	${0.69}_{-0.01}^{+0.01}$	$\equiv$ 0	$\equiv$ $-1$	$\equiv$ 0	${3.20}_{-0.16}^{+0.15}$	$\equiv$ 0.06
$N_{\mathrm{eff}}\Lambda$CDM	Planck + BAO	${69.0}_{-1.1}^{+1.1}$	${0.30}_{-0.01}^{+0.01}$	${0.70}_{-0.01}^{+0.01}$	$\equiv$ 0	$\equiv$ $-1$	$\equiv$ 0	${3.23}_{-0.17}^{+0.17}$	$\equiv$ 0.06
$m_{\nu }\Lambda$CDM	Planck	${68.0}_{-0.7}^{+0.7}$	${0.31}_{-0.01}^{+0.01}$	${0.69}_{-0.01}^{+0.01}$	$\equiv$ 0	$\equiv$ $-1$	$\equiv$ 0	$\equiv$ 3.046	${0.03}_{-0.02}^{+0.05}$
$m_{\nu }\Lambda$CDM	Planck + CMBL	${68.0}_{-0.7}^{+0.6}$	${0.31}_{-0.01}^{+0.01}$	${0.69}_{-0.01}^{+0.01}$	$\equiv$ 0	$\equiv$ $-1$	$\equiv$ 0	$\equiv$ 3.046	${0.01}_{-0.02}^{+0.05}$
$m_{\nu }\Lambda$CDM	Planck + BAO	${68.0}_{-0.5}^{+0.5}$	${0.31}_{-0.01}^{+0.01}$	${0.69}_{-0.01}^{+0.01}$	$\equiv$ 0	$\equiv$ $-1$	$\equiv$ 0	$\equiv$ 3.046	${0.03}_{-0.02}^{+0.04}$
$N_{\mathrm{eff}}{m}_{\nu }\Lambda$CDM	Planck	${69.1}_{-1.3}^{+1.3}$	${0.30}_{-0.01}^{+0.01}$	${0.70}_{-0.01}^{+0.01}$	$\equiv$ 0	$\equiv$ $-1$	$\equiv$ 0	${3.21}_{-0.17}^{+0.20}$	${0.03}_{-0.03}^{+0.05}$
$N_{\mathrm{eff}}{m}_{\nu }\Lambda$CDM	Planck + BAO	${69.0}_{-1.1}^{+1.0}$	${0.30}_{-0.01}^{+0.01}$	${0.70}_{-0.01}^{+0.01}$	$\equiv$ 0	$\equiv$ $-1$	$\equiv$ 0	${3.21}_{-0.17}^{+0.16}$	${0.03}_{-0.02}^{+0.05}$
$w_0{w}_a$CDM	Planck + BAO	${67.1}_{-2.5}^{+3.0}$	${0.32}_{-0.03}^{+0.03}$	${0.68}_{-0.03}^{+0.03}$	$\equiv$ 0	$-{0.73}_{-0.30}^{+0.29}$	$-{1.02}_{-0.88}^{+0.80}$	$\equiv$ 3.046	$\equiv$ 0.06

Table 8. Table 8: Cosmological parameters for various cosmologies when anchoring type Ia supernovae distances with distances from time-delay cosmography.

JLA sample (Betoule et al., 2014) and H0LiCOW
Model	$H_0$ (km s-1 Mpc-1)	${\Omega }_m$	${\Omega }_{\Lambda }$ or ${\Omega }_{\mathrm{DE}}$	${\Omega }_k$	$w$ or $w_0$	$w_a$
$\Lambda$CDM	${73.6}_{-1.8}^{+1.6}$	${0.30}_{-0.03}^{+0.03}$	${0.70}_{-0.03}^{+0.03}$	$\equiv$ 0	$\equiv$ $-1$	$\equiv$ 0
$w$CDM	${73.9}_{-2.7}^{+2.6}$	${0.31}_{-0.10}^{+0.08}$	${0.69}_{-0.08}^{+0.10}$	$\equiv$ 0	$-{1.05}_{-0.28}^{+0.23}$	$\equiv$ 0
$w_0{w}_a$CDM	${74.2}_{-2.6}^{+2.5}$	${0.34}_{-0.10}^{+0.08}$	${0.66}_{-0.08}^{+0.10}$	$\equiv$ 0	$-{1.04}_{-0.28}^{+0.22}$	$-{0.43}_{-1.07}^{+1.10}$
o$\Lambda$CDM	${75.2}_{-2.1}^{+1.8}$	${0.18}_{-0.06}^{+0.08}$	${0.53}_{-0.09}^{+0.12}$	${0.29}_{-0.19}^{+0.14}$	$\equiv$ $-1$	$\equiv$ 0
o$w$CDM	${77.5}_{-2.9}^{+2.5}$	${0.22}_{-0.07}^{+0.07}$	${0.41}_{-0.08}^{+0.14}$	${0.37}_{-0.17}^{+0.09}$	$-{1.46}_{-0.53}^{+0.42}$	$\equiv$ 0
o$w_0{w}_a$CDM	${77.7}_{-2.8}^{+2.4}$	${0.23}_{-0.06}^{+0.07}$	${0.40}_{-0.07}^{+0.13}$	${0.38}_{-0.16}^{+0.09}$	$-{1.49}_{-0.52}^{+0.44}$	$-{0.16}_{-1.24}^{+1.24}$
Pantheon sample (Scolnic et al., 2018) and H0LiCOW
Model	$H_0$ (km s-1 Mpc-1)	${\Omega }_m$	${\Omega }_{\Lambda }$ or ${\Omega }_{\mathrm{DE}}$	${\Omega }_k$	$w$ or $w_0$	$w_a$
$\Lambda$CDM	${73.6}_{-1.8}^{+1.6}$	${0.30}_{-0.02}^{+0.02}$	${0.70}_{-0.02}^{+0.02}$	$\equiv$ 0	$\equiv$ $-1$	$\equiv$ 0
$w$CDM	${74.9}_{-2.4}^{+2.2}$	${0.35}_{-0.06}^{+0.05}$	${0.65}_{-0.05}^{+0.06}$	$\equiv$ 0	$-{1.17}_{-0.22}^{+0.19}$	$\equiv$ 0
$w_0{w}_a$CDM	${75.0}_{-2.3}^{+2.2}$	${0.37}_{-0.07}^{+0.05}$	${0.63}_{-0.05}^{+0.07}$	$\equiv$ 0	$-{1.15}_{-0.23}^{+0.19}$	$-{0.55}_{-1.01}^{+1.22}$
o$\Lambda$CDM	${73.8}_{-2.1}^{+1.9}$	${0.28}_{-0.07}^{+0.07}$	${0.67}_{-0.11}^{+0.11}$	${0.04}_{-0.17}^{+0.18}$	$\equiv$ $-1$	$\equiv$ 0
o$w$CDM	${77.4}_{-3.0}^{+2.5}$	${0.28}_{-0.06}^{+0.07}$	${0.45}_{-0.10}^{+0.17}$	${0.26}_{-0.21}^{+0.16}$	$-{1.52}_{-0.52}^{+0.39}$	$\equiv$ 0
o$w_0{w}_a$CDM	${77.4}_{-3.0}^{+2.5}$	${0.28}_{-0.07}^{+0.08}$	${0.45}_{-0.09}^{+0.17}$	${0.26}_{-0.21}^{+0.15}$	$-{1.50}_{-0.47}^{+0.40}$	$-{0.31}_{-1.17}^{+1.31}$

Equations

$$t(\mathbf{\theta },\mathbf{\beta })=\frac{D_{\Delta t}}{c}\left[\frac{(\mathbf{\theta }-\mathbf{\beta }{)}^2}{2}-\psi (\mathbf{\theta })\right],$$

$$D_{\Delta t}\equiv (1+z_{\mathrm{d}})\frac{D_{\mathrm{d}}{D}_{\mathrm{s}}}{D_{\mathrm{ds}}},$$

$$\Delta t_{ij}=\frac{D_{\Delta t}}{c}\left[\frac{({\mathbf{\theta }}_i-\mathbf{\beta }{)}^2}{2}-\psi ({\mathbf{\theta }}_i)-\frac{({\mathbf{\theta }}_j-\mathbf{\beta }{)}^2}{2}+\psi ({\mathbf{\theta }}_j)\right],$$

$$D_{\Delta t}=\frac{D_{\Delta t}^{\mathrm{model}}}{1-{\kappa }_{\mathrm{ext}}}.$$

$$P(\mathbf{\xi }|{\mathbf{d}}_{\mathrm{img}},\mathbf{\Delta }\mathbf{t},\mathbf{\sigma },{\mathbf{d}}_{\mathrm{LOS}},\mathbf{A})$$

$$P({\mathbf{d}}_{\mathrm{img}},\mathbf{\Delta }\mathbf{t},\mathbf{\sigma },{\mathbf{d}}_{\mathrm{LOS}}|\mathbf{\xi },\mathbf{A})$$

$$P(\mathbf{\pi }|{\mathbf{d}}_{\mathrm{img}},\mathbf{\Delta }\mathbf{t},\mathbf{\sigma },{\mathbf{d}}_{\mathrm{LOS}})$$

$$F_{ij}$$

$$=\frac{⟨L_i{L}_j⟩}{⟨L_i⟩⟨L_j⟩},$$

Full text

††pubyear: 2019††pagerange: H0LiCOW XIII. A $2.4\%$ measurement of $H_{0}$ from lensed quasars: $5.3\sigma$ tension between early and late-Universe probes–LABEL:lastpage

1 Introduction

The flat $\Lambda$ cold dark matter ($\Lambda$CDM) cosmological model has proven to be remarkably successful at describing the Universe as measured by a wide range of experiments, particularly observations of the cosmic microwave background (CMB). The final results from the Planck mission (Planck Collaboration et al., 2018a) provide the most precise constraints on cosmological parameters to date from CMB observations (Planck Collaboration et al., 2018b). However, relaxing the flat $\Lambda$CDM assumption by introducing additional complexities, such as non-zero curvature, an equation of state parameter $w\neq-1$, or a time-varying $w$, leads to much weaker constraints and large degeneracies among the various cosmological parameters. In particular, many parameters become degenerate with the Hubble constant, $H_{0}$, which sets the present-day expansion rate of the universe. $H_{0}$ cannot be constrained directly from CMB observations, but must be inferred by first assuming a cosmological model. In this context, measuring $H_{0}$ independent of CMB observations is one of the most important complementary probes for understanding the nature of the Universe (Weinberg et al., 2013).

The most well-established method for measuring $H_{0}$ is through observations of type Ia supernovae (SNe). Type Ia SNe are “standardizable candles" in that their luminosities, and thus their absolute distances, can be determined by the evolution of their light curves, and therefore can be used to infer $H_{0}$ from the slope of their distance-redshift relation. Type Ia SNe luminosities are typically calibrated via the “distance ladder" (e.g., Sandage et al., 2006; Freedman et al., 2012; Riess et al., 2016, 2018; Riess et al., 2019), in which parallax measurements are used to determine distances to nearby Cepheid variable stars (which have a known period-luminosity relation) and are in turn used to determine distances to type Ia SNe in the Hubble flow.

Recent determinations of $H_{0}$ from the Supernovae, $H_{0}$, for the Equation of State of Dark Energy (SH0ES; Riess et al., 2016) collaboration using this method are in tension with the Planck CMB measurements under the flat $\Lambda$CDM model (e.g., Bernal et al., 2016; Freedman, 2017). The latest SH0ES result finds $H_{0}=74.03\pm 1.42~{}\mathrm{km~{}s^{-1}~{}Mpc^{-1}}$ (Riess et al., 2019), which differs from the Planck flat $\Lambda$CDM result111Baseline $\Lambda$CDM chains with baseline likelihoods (based on plikHM_TTTEEE_lowl_lowE) of $H_{0}=67.4\pm 0.5~{}\mathrm{km~{}s^{-1}~{}Mpc^{-1}}$ (Planck Collaboration et al., 2018b) by $4.4\sigma$. Possible systematic errors in one or both methods may resolve this tension (e.g., Rigault et al., 2015, 2018), but investigations thus far have yet to conclusively identify any such systematic (e.g., Addison et al., 2018; Jones et al., 2018; Roman et al., 2018; Camarena & Marra, 2019; Rose et al., 2019). Furthermore, independent determinations of $H_{0}$ using the “inverse distance ladder" method (e.g., Aubourg et al., 2015; Cuesta et al., 2015; Macaulay et al., 2019) are in agreement with the Planck value, although this depends on assumptions of the physical scale of the sound horizon (e.g., Aylor et al., 2019; Macaulay et al., 2019; Arendse et al., 2019a, b). Other methods, such as CMB polarization measurements (e.g., Henning et al., 2018), galaxy clustering (e.g., Abbott et al., 2018a), water masers (e.g., Herrnstein et al., 1999; Humphreys et al., 2013; Braatz et al., 2018), X-ray observations of SZ galaxy clusters (e.g., Silk & White, 1978; Reese et al., 2002; Bonamente et al., 2006; Kozmanyan et al., 2019), the Balmer line $L-\sigma$ relation of HII galaxies (e.g., Melnick et al., 2000; Chávez et al., 2012; González-Morán et al., 2019), extragalactic background light attenuation (e.g., Salamon et al., 1994; Domínguez & Prada, 2013; Domínguez et al., 2019), type IIP supernova expanding photospheres (e.g., Schmidt et al., 1994; Gall et al., 2016), and gravitational waves (e.g., Abbott et al., 2017; Feeney et al., 2019; Soares-Santos et al., 2019), have yet to resolve the $H_{0}$ discrepancy, as their precision is not yet comparable to Planck or SH0ES, or they require additional assumptions. If unresolved, this tension may force the rejection of the flat $\Lambda$CDM model and indicate new physics that must be incorporated into our understanding of cosmology.

After the submission of this paper, an alternate calibration of the distance ladder using the “tip of the red giant branch" (TRGB) method by the Carnegie-Chicago Hubble Program (CCHP; Beaton et al., 2016) found an intermediate value of $H_{0}=69.8\pm 1.9~{}\mathrm{km~{}s^{-1}~{}Mpc^{-1}}$ (Freedman et al., 2019). However, this measurement is not fully independent of SH0ES since they share some calibrating sources (i.e., galaxies hosting SNe that are close enough for Cepheid and/or TRGB distance measurements), and there is an ongoing debate about the results from this method (e.g., Yuan et al., 2019), further highlighting the need for additional independent probes (see Verde et al., 2019 for a recent review of the field).

Gravitational lensing offers an independent method of determining $H_{0}$. When a background object (the “source") is gravitationally lensed into multiple images by an intervening mass (the “lens"), light rays emitted from the source will take different paths through space-time at the different image positions. Because these paths have different lengths and pass through different gravitational potentials, light rays emitted from the source at the same time will arrive at the observer at different times depending on which image it arrives at. If the source is variable, this “time delay" between multiple images can be measured by monitoring the lens and looking for flux variations corresponding to the same source event. The time delay is related to a quantity referred to as the “time-delay distance", $D_{\Delta t}$, and depends on the mass distribution in the lensing object, the mass distribution along the line of sight (LOS), and cosmological parameters. $D_{\Delta t}$ is primarily sensitive to $H_{0}$, although there is a weak dependence on other parameters (e.g., Coe & Moustakas, 2009; Linder, 2011; Treu & Marshall, 2016). This one-step method is completely independent of and complementary to the CMB and the distance ladder. The distances probed by time-delay cosmography are also larger than those from the distance ladder, making this method immune to a monopole in the bulk velocity field of the local Universe (i.e., a “Hubble bubble").

This method of using gravitational lens time delays to measure $H_{0}$ was first proposed by Refsdal (1964), who suggested using lensed SNe for this purpose. In practice, finding lensed SNe with resolved images is extremely rare, with only two such lenses having been discovered to date (Kelly et al., 2015; Goobar et al., 2017). While the prospect of discovering more lensed SNe in future imaging surveys and measuring their time delays is promising (e.g., Oguri & Marshall, 2010; Goldstein & Nugent, 2017; Goldstein et al., 2018; Huber et al., 2019; Wojtak et al., 2019), lensed quasars have generally been used to constrain $H_{0}$ in this manner (e.g., Vanderriest et al., 1989; Keeton & Kochanek, 1997; Schechter et al., 1997; Kochanek, 2003; Koopmans et al., 2003; Saha et al., 2006; Oguri, 2007; Vuissoz et al., 2008; Fadely et al., 2010; Suyu et al., 2010, 2013; Sereno & Paraficz, 2014; Rathna Kumar et al., 2015; Birrer et al., 2016, 2019; Chen et al., 2016; Wong et al., 2017; Bonvin et al., 2017) due to their brightness and variable nature.

Measuring $H_{0}$ from lensed quasars through this method requires a variety of observational data. Long-term dedicated photometric monitoring of the lens is needed to obtain accurate time delays (e.g., Bonvin et al., 2017, 2018). Several years of monitoring are generally required to overcome microlensing variability, although Courbin et al. (2018) recently demonstrated that delays could be measured from just one year of monitoring owing to high photometric accuracy (milli-mag) and observing cadence (daily). In addition, deep high-resolution imaging of the lens is required to observe the extended images of the quasar host galaxy, which is needed to break degeneracies in the lens modeling between the mass profile and the underlying cosmology (e.g., Kochanek, 2002; Koopmans et al., 2003; Dye & Warren, 2005). Furthermore, to mitigate the effects of the mass-sheet degeneracy (e.g., Falco et al., 1985; Gorenstein et al., 1988; Saha, 2000; Schneider & Sluse, 2013; Xu et al., 2016), it is important to obtain a measurement of the lens galaxy’s velocity dispersion (e.g. Treu & Koopmans, 2002; Koopmans et al., 2003; Koopmans, 2004; Sonnenfeld, 2018). Finally, observational data to constrain the mass along the LOS to the lens are needed to estimate the external convergence, $\kappa_{\rm ext}$, which can bias the inferred $D_{\Delta t}$ if unaccounted for (e.g., Collett et al., 2013; Greene et al., 2013; McCully et al., 2014, 2017; Sluse et al., 2017; Rusu et al., 2017; Tihhonova et al., 2018).

The $H_{0}$ Lenses in COSMOGRAIL’s Wellspring (H0LiCOW) collaboration (Suyu et al., 2017, hereafter H0LiCOW I) has provided the strongest constraints on $H_{0}$ to date from time-delay cosmography. Our most recent measurements had constrained $H_{0}$ to $3.0\%$ precision for a flat $\Lambda$CDM cosmology from a combination of four lensed quasars (Birrer et al., 2019, hereafter H0LiCOW IX). We attain this precision by taking advantage of our substantial dataset, which includes accurate time-delay measurements from the COSmological MOnitoring of GRAvItational Lenses (COSMOGRAIL; Courbin et al., 2005; Eigenbrod et al., 2005; Bonvin et al., 2018) project and radio-wavelength monitoring (Fassnacht et al., 2002), deep Hubble Space Telescope (HST) and/or ground-based adaptive optics (AO) imaging (H0LiCOW I, IX, Chen et al., 2016, 2019), spectroscopy of the lens galaxy to measure its velocity dispersion (e.g., Sluse et al., 2019, hereafter H0LiCOW X), and deep wide-field spectroscopy and imaging to characterize the LOS in these systems (e.g., Sluse et al., 2017; Rusu et al., 2017, hereafter H0LiCOW II and H0LiCOW III, respectively).

In this milestone paper, we present the latest constraints on $H_{0}$ from H0LiCOW from a combined sample of six lensed quasars. Two of the four lenses analyzed previously, HE 0435$-$1223 (Wong et al., 2017, hereafter H0LiCOW IV) and RXJ1131$-$1231 (Suyu et al., 2014), are reanalyzed using new AO data (Chen et al., 2019). We add constraints from two newly-analyzed systems – PG 1115$+$080 (Chen et al., 2019) and WFI2033$-$4723 (Rusu et al., 2019, hereafter H0LiCOW XII) – to provide the tightest $H_{0}$ constraints to date from time-delay cosmography.

This paper is organized as follows. In Section 2, we summarize the theory behind using time-delay cosmography to infer the time-delay distance, which is inversely proportional to $H_{0}$. In Section 3, we present our lens sample and describe how our data and analysis methods allow us to constrain $H_{0}$ to a precision greater than what has previously been possible from time-delay strong lensing. In Section 4, we verify that our lenses are consistent with each other so that we can combine them for our cosmological inference. We present our main results for flat $\Lambda$CDM and more complex cosmologies in Section 5. In Section 6, we discuss the tension between early-Universe and late-Universe probes of $H_{0}$. We summarize our findings in Section 7. Throughout this paper, all magnitudes given are on the AB system. All parameter constraints given are medians and 16th and 84th percentiles unless otherwise stated.

2 Time-Delay Cosmography

In this section, we summarize the theoretical background of time-delay cosmography and how to infer $H_{0}$. We refer readers to recent reviews (e.g., Treu & Marshall, 2016; Suyu et al., 2018) for more details.

When light rays from a background source are deflected by an intervening lensing mass, the light travel time from the source to the observer depends on both their path length and the gravitational potential they traverse. For a single lens plane, the excess time delay of an image at an angular position $\bm{\theta}=(\theta_{1},\theta_{2})$ with corresponding source position $\bm{\beta}=(\beta_{1},\beta_{2})$ relative to the case of no lensing is

[TABLE]

where $D_{\Delta t}$ is the time-delay distance and $\psi(\bm{\theta})$ is the lens potential. The time-delay distance (Refsdal, 1964; Schneider et al., 1992; Suyu et al., 2010) is defined as

[TABLE]

where $z_{\rm d}$ is the lens redshift, $D_{\rm d}$ is the angular diameter distance to the lens, $D_{\rm s}$ is the angular diameter distance to the source, and $D_{\rm ds}$ is the angular diameter distance between the lens and the source. $D_{\Delta t}$ has units of distance and is inversely proportional to $H_{0}$, with weak dependence on other cosmological parameters.

If the alignment between the background source and the foreground lens is close enough, multiple images of the same background source are formed. Light rays reaching the observer will have different excess time delays depending on which image they are observed at. The time delay between two images of such a lens, $\Delta t_{ij}$, is the difference of their excess time delays,

[TABLE]

where $\bm{\theta}_{i}$ and $\bm{\theta}_{j}$ are the positions of images $i$ and $j$, respectively, in the image plane. If the source is variable on short timescales, it is possible to monitor the fluxes of the images and measure the time delay, $\Delta t_{ij}$, between them (e.g., Vanderriest et al., 1989; Schechter et al., 1997; Fassnacht et al., 1999; Fassnacht et al., 2002; Kochanek et al., 2006; Courbin et al., 2011). The lens potentials at the image positions, $\psi(\bm{\theta}_{i})$ and $\psi(\bm{\theta}_{j})$, as well as at the source position, $\bm{\beta}$, can be determined from a mass model of the system. With a measurement of $\Delta t_{ij}$ and an accurate lens model to determine $\psi(\bm{\theta})$, it is possible to determine $D_{\Delta t}$. By further assuming a cosmological model, $D_{\Delta t}$ can be converted into an inference on $H_{0}$.

If there are multiple lenses at different redshifts between the source and the observer, the observed time delays depend on combinations of the angular diameter distances among the observer, the multiple lens planes, and the source. In this case, the image positions are described by the multiplane lens equation (e.g., Blandford & Narayan, 1986; Kovner, 1987; Schneider et al., 1992; Petters et al., 2001; Collett & Auger, 2014; McCully et al., 2014), and the observed time delays are no longer proportional to a single unique time-delay distance. However, it is often the case that the mass in a single lens plane dominates the lensing effect, and the observed time delays are primarily sensitive to the time-delay distance (Equation 2) with the deflector redshift as that of the primary lens plane. This is the case for all of the lenses in the H0LiCOW sample (Section 3.1). The results for any individual system can thus be interpreted as a constraint on $D_{\Delta t}(z_{\mathrm{d}},z_{\mathrm{s}})$, which we refer to as the “effective time-delay distance". Hereafter, $D_{\Delta t}$ refers to the effective time-delay distance (for applicable systems) unless otherwise indicated.

In addition to mass that is explicitly included in the lens model, all other mass along the LOS between the observer and the source contributes to the lens potential that the light rays traverse. This causes additional focusing and defocusing of the rays and can affect the observed time delays (e.g., Seljak, 1994). If unaccounted for, this can lead to biased inferences of $D_{\Delta t}$. If the effects of the perturbing LOS masses are small, they can be approximated by an external convergence in the lens plane $\kappa_{\rm ext}$ (e.g., Keeton, 2003; McCully et al., 2014). The true $D_{\Delta t}$ is related to the time-delay distance inferred from the lens model and measured time delays, $D_{\Delta t}^{\mathrm{model}}$, by the relation

[TABLE]

$\kappa_{\rm ext}$ is defined such that its average value across the sky is zero. In principle, if lenses are randomly distributed, the effect of $\kappa_{\rm ext}$ should average out over a sufficiently large sample. However, the cross section for strong lensing scales as $\sigma^{4}$, where $\sigma$ is the velocity dispersion of the lens galaxy. As a result, lenses are biased toward the most massive galaxies, which are known to cluster (e.g., Dressler, 1980). Indeed, lens galaxies generally lie in overdense environments and lines of sight relative to typical fields (e.g., Treu et al., 2009; Fassnacht et al., 2011; Wong et al., 2018), meaning that $\kappa_{\rm ext}$ will lead to a bias on $D_{\Delta t}$ and needs to be corrected for. $\kappa_{\rm ext}$ cannot, in general, be constrained from the lens model due to the mass-sheet degeneracy (Falco et al., 1985; Gorenstein et al., 1988; Saha, 2000), in which the addition of a uniform mass sheet associated with a rescaling of the mass normalization of the strong lens galaxy and the coordinates in the source plane can modify the product of the time delays and $H_{0}$ but leave other observables unchanged. $\kappa_{\rm ext}$ must instead be estimated through other methods, such as studies of the lens environment or the use of lens stellar kinematics (as noted in the previous section).

With kinematic information on the lens galaxy, it is possible to determine the angular diameter distance to the lens, $D_{\rm d}$, independent of $\kappa_{\rm ext}$ (Paraficz & Hjorth, 2009; Jee et al., 2015). Although the constraints from $D_{\rm d}$ are generally weaker than those from $D_{\Delta t}$, it can break degeneracies among cosmological parameters, particularly for models beyond flat $\Lambda$CDM. In particular, it can break the degeneracy between curvature ($\Omega_{\rm k}$) and the time-varying equation of state parameter of dark energy ($w$) (Jee et al., 2016). The combination of lensing, time delays, and lens kinematic data thus provides a joint constraint on $D_{\Delta t}$ and $D_{\rm d}$ in cases of single strong-lensing planes (see e.g., Birrer et al., 2016; Chen et al., 2019, H0LiCOW IX for more details). These constraints on lensing distances, together with the redshifts of the lenses and sources, then allow us to infer cosmological parameter values for a given cosmological model.

3 Overview of the H0LiCOW Analysis

In this section, we provide a brief summary of the H0LiCOW analysis, including the sample of lenses analyzed to date (Section 3.1), as well as the various components that go into determining cosmological constraints from each lens.

3.1 Lens Sample

Our sample of strongly-lensed quasars comprises six systems analyzed to date by H0LiCOW and collaborators. The six lenses are listed in Table 1, and we show multicolor high-resolution images of them in Figure 1. Each of these systems have been modeled using constraints from high-resolution HST and/or ground-based AO imaging data, time-delay measurements from COSMOGRAIL and Fassnacht et al. (2002), and kinematics from ground-based spectroscopy. In addition, we have constrained $\kappa_{\rm ext}$ in these systems from a wide-field imaging and spectroscopic campaign, as detailed in Section 3.4.

The original H0LiCOW sample (see H0LiCOW I) consists of five lenses, and it was later decided to expand the sample to include four additional systems with HST imaging (PID:14254, PI: T. Treu), placing an emphasis on double-image lens systems (doubles), which yield fewer constraints per system but are more abundant on the sky. Of the five doubles (one from the original sample plus the additional four), SDSS 1206$+$4332 was analyzed first (H0LiCOW IX) because part of the quasar’s host galaxy is quadruply-imaged, providing additional constraints for lens modeling. PG 1115$+$080 was observed with AO imaging from Keck/NIRC2 as part of the Strong lensing at High Angular Resolution Program (SHARP; Fassnacht et al. in preparation) and was incorporated into the H0LICOW sample later (Chen et al., 2019).

The current sample of six systems used in this work include the four quadruple-image lenses (quads) from the original sample, plus SDSS 1206$+$4332 and PG 1115$+$080. These six lenses span a range of lens and source redshifts, as well as a range of image configurations (e.g., double, cross, fold, cusp). Having a sample that spanned a range in these parameters was a consideration in the selection of which lenses to analyze first, as there may be systematics that depend on such factors, and we want to account for them in our analysis (see Ding et al., 2018, who attempt to address these issues based on simulated data).

3.2 Time Delay Measurement

Out of the six lenses of the H0LICOW sample, all except for B1608$+$656 have been monitored in optical by the COSMOGRAIL collaboration from several facilities with 1m and 2m-size telescopes. Several seasons of monitoring are needed in order to disentangle the variations due to microlensing in which brightening or dimming of the quasar images by stars in the lens galaxy can mimic intrinsic features in the light curves.

From the monitoring data, COSMOGRAIL measures time delays using numerical curve-shifting techniques, which fit a function to the light curve of each quasar image and find the time shifts that minimize the differences among them (Tewes et al., 2013a; Bonvin et al., 2019). These techniques are made publicly available as a python package named PyCS 222Available at http://www.cosmograil.org, which also provides tools to estimate the time delays uncertainties in the presence of microlensing. The package was tested on simulated light curves reproducing the COSMOGRAIL data with similar sampling and photometric noise in a blind time delay challenge (Liao et al., 2015). Bonvin et al. (2016) demonstrated the robustness of the PyCS curve-shifting techniques by recovering the time delays at a precision of $\sim$ 3% on average with negligible systematic bias.

Tewes et al. (2013b) applied these techniques to RXJ1131$-$1231  and measured the longest time delay to 1.5% precision (1$\sigma$). The time delay of SDSS 1206$+$4332 was also measured with PyCS; Eulaers et al. (2013) obtained a time delay between the two multiple images of $\Delta t_{AB}=111.3\pm 3$ days, with image A leading image B. Birrer et al. (2019) re-analyzed the same monitoring data with updated and independent curve-shifting techniques and confirmed this result. For HE 0435$-$1223, the latest time delay measurement was obtained with the 13 year-long light curves of the COSMOGRAIL program at 6.5% precision on the longest time delay (Bonvin et al., 2017).

Recently, Courbin et al. (2018) demonstrated that a high-cadence and high signal-to-noise (S/N) monitoring campaign can also disentangle the microlensing variability from the intrinsic variability signal by catching small variations of the quasar that happen on timescales much shorter than the typical microlensing variability. It is therefore possible to disentangle the intrinsic signal of the quasar from the microlensing signal in a single season. High-cadence data were used for WFI2033$-$4723 and PG 1115$+$080 to measure time delays at a few percent precision in one season. These results are in agreement with the time delays measured from decade-long COSMOGRAIL light curves and are combined in the final estimate (Bonvin et al., 2018, 2019).

The remaining lens of the sample, B1608$+$656  was monitored by Fassnacht et al. (1999); Fassnacht et al. (2002) with radio observations from the Very Large Array over three seasons. All three independent time delays between the multiple images were measured to a precision of a few percent.

A complicating factor in converting the observed time delays to a cosmological constraint is the so-called “microlensing time-delay" effect (Tie & Kochanek, 2018). The estimation of this effect is based on the lamp-post model, which predicts delayed emission across the quasar accretion disk from a central driving source. Different regions of the disk can then be magnified by the microlenses differently in each of the multiple images. This reweighting of the delayed emission across the accretion disk could lead to a change in the measured time delay. As the microlensing changes with time, this could lead to a variation in the measured time delays from season to season. There is no evidence of this effect based on our current data, so our main cosmological results do not depend on it. Nonetheless, we quantify this factor for different speculative models (Bonvin et al., 2018, 2019; Chen et al., 2018) in the latest H0LiCOW lens models (Birrer et al., 2019; Chen et al., 2019; Rusu et al., 2019).

3.3 Lens Modeling

The primary lens modeling code used to model the majority of the H0LiCOW lenses is Glee (Suyu & Halkola, 2010; Suyu et al., 2012), although SDSS 1206$+$4332 (H0LiCOW IX) is analyzed using the Lenstronomy code (Birrer et al., 2015; Birrer & Amara, 2018). Both codes model the lens galaxy light as parameterized profiles and fit the lensed quasar image positions and surface brightness distribution of the quasar host galaxy. The primary difference between the codes is that Glee reconstructs the source on a pixelized grid with regularization, whereas Lenstronomy describes the source as a parameterized profile with additional shapelet functions (Refregier, 2003; Birrer et al., 2015). The use of two independent codes is meant to provide a check on lens modeling code systematics. This would ideally require both codes to be tested on the same system, as will be done for future lens analyses (Shajib et al. in preparation; Yıldırım et al. in preparation).

We use two main parameterizations of the lens galaxy in our models: a singular elliptical power-law model, and a composite model consisting of a baryonic component linked to the stellar light distribution plus an elliptical NFW (Navarro et al., 1996) halo representing the dark matter component. For B1608$+$656, which shows two interacting lens galaxies, rather than using the two main parameterizations, we started from the power-law model and performed a pixelated lens potential reconstruction to allow for flexibility, finding small ($\sim 2\%$) potential corrections and thus validating the use of power-law model. Galaxies along the LOS that are deemed to be significant perturbers are included in the model (Section 3.4) through the full multiplane lens equation, and we also include an external shear in the main lens plane. When available, we use the measured velocity dispersions of the lens galaxy and significant perturbers as additional constraints.

To account for systematic effects arising from modeling choices in areas such as the lens parameterization, the source reconstruction, the weighting of the pixels in the image plane, etc., we run multiple models where we vary these choices and combine them in our final inference. In our initial analyses of the first three H0LiCOW lenses (Suyu et al., 2010, 2013, 2014, H0LiCOW IV), we marginalize over these (discrete) modeling choices in deriving the posterior probability density function (PDF) of $D_{\Delta t}$. For models that fit equally well to the data within their modeling uncertainties, we conservatively add their posterior distributions of $D_{\Delta t}$ with equal weight, given uniform prior on these modeling choices. For subsequent analyses, including the reanalysis of RXJ1131$-$1231 and HE 0435$-$1223 (Chen et al., 2019), we weight by the Bayesian Information Criterion (BIC), following the procedure described in H0LiCOW IX.

We have measurements of the velocity dispersions of the lens galaxies in our sample from high-resolution spectroscopy, which are used to mitigate degeneracies in the mass modeling. The velocity dispersion can be combined with the lensing constraints to estimate the angular diameter distancesto the lens (see Section 2), in the systems (B1608$+$656, RXJ1131$-$1231, SDSS 1206$+$4332 and PG 1115$+$080) that could be well-modeled by single lens plane (i.e., without multi lens plane modeling).

3.4 LOS Structure and External Convergence

In accounting for the effects of LOS structure, there are two primary types of perturbations that need to be dealt with. The first is the effect of structures that affect the lens potential significantly enough that they cannot be approximated by their tidal perturbations, but must instead be included explicitly in the lens model (Section 3.3). The second is the combined effect of all other LOS structures, which can be approximated by a $\kappa_{\rm ext}$ term. Quantifying and accounting for both types of perturbers requires spectroscopic and photometric data on LOS objects projected nearby the lens.

The effect of a LOS perturber on the lens potential can be quantified by the “flexion shift" (McCully et al., 2014). In general, objects have a larger flexion shift if they are more massive, projected more closely to the lens, and are either at the redshift of the lens or at a lower redshift, as opposed to a higher redshift (McCully et al., 2017). As a result, we focus our spectroscopic campaign on bright galaxies projected close to the lens and include them in our lens model if their calculated flexion shift is large. When their velocity dispersions are measurable from our data, we use them to set a prior on their Einstein radii during the modeling procedure (e.g., H0LiCOW X). Alternatively, a prior on the Einstein radius is set from scaling relations with luminosity (e.g., H0LiCOW IX). An overview of our spectroscopic campaigns is provided in H0LiCOW I for the majority of our sample of lenses, and in H0LiCOW IX for SDSS 1206$+$4332. For PG 1115$+$080, where we have not conducted spectroscopic follow-up, we use the compilation of redshifts presented by Wilson et al. (2016), which was also useful in providing additional redshifts for some of the other lens fields. Our spectroscopic data provide accurate redshifts for hundreds of bright galaxies as far as $\sim 10\arcmin$ away from the lens systems, allowing us to further quantify the properties of larger structures such as galaxy groups and clusters. Since the analysis of HE 0435$-$1223 (H0LiCOW II), we use an adaptive group-finding algorithm that uses spectroscopic data to identify peaks in redshift space and refines group membership based on proximity between potential group galaxies (e.g., LOS velocity dispersion and centroid; see H0LiCOW X for details). In cases where there are potentially significant effects from such structures, we run systematics tests (Section 3.3) where we include these structures as spherical NFW halos in our models.

To correct for the statistical effect of $\kappa_{\rm ext}$ due to LOS structure, we use a (weighted) galaxy number counts technique (e.g., Greene et al., 2013; Suyu et al., 2013, H0LiCOW III, IX). We count galaxies projected within a fixed distance of the lens and above some flux threshold, weighted by various quantities such as external shear, projected distance, stellar mass, redshift, etc. (see H0LiCOW III, XII for details). We then compare these number counts to those measured analogously along random lines of sight in a control survey to determine the relative over-/under-density of the lens field. Finally, we use the Millennium simulation (Springel et al., 2005) to identify lines of sight that have a similar relative number count densities and build a PDF of $\kappa_{\rm ext}$ determined from ray-tracing (Hilbert et al., 2009), which we apply in post-processing. Alongside our spectroscopic campaign, we have conducted our own multi-band, wide-field imaging campaign to gather data of sufficient depth and spatial coverage to enable this analysis. Our data typically consist of multi-band ultraviolet to infrared observations in good seeing conditions, which we use to perform a galaxy-star classification and to measure physical quantities such as redshifts and stellar masses for the galaxies projected within $<120\arcsec$ of the lenses, and down to $i\sim 23-24$ mag. Exceptions are B1608$+$656 and RXJ1131$-$1231, where we used single-band HST data within $<45\arcsec$ (Fassnacht et al., 2011), and PG 1115$+$080, where we coadded multiple exposures of the data used to measure the time delays. With the exception of the first two lenses we analyzed (B1608$+$656 and RXJ1131$-$1231) where we used archival HST data as control fields, we employed the larger-scale CFHTLenS (Heymans et al., 2012).

Our technique has evolved over the years such that for B1608$+$656 we only employed unweighted number counts to constrain $\kappa_{\rm ext}$, whereas for the remaining lenses we also used constraints from the inferred external shear values of the lens models.333An exception is SDSS 1206$+$4332, where the use of external shear has a negligible effect on $\kappa_{\rm ext}$. Since HE 0435$-$1223 (H0LiCOW III), we have further used combinations of weighted counts to tighten the $\kappa_{\rm ext}$ PDF, which for the three latest lenses have included weighted number counts measured in multiple apertures.

In future work, we plan to return to previous lenses in order to enforce consistency of technique, as well as to further refine our technique by better accounting for lens-lens coupling between the primary lens and LOS structures in the convergence maps from Hilbert et al. (2009), employing other cosmological simulations with different assumed cosmology, and also using new techniques which move away from the statistical approach and have the potential to significantly tighten the $\kappa_{\rm ext}$ PDF (e.g., McCully et al., 2017). In addition, in some of the H0LiCOW systems, we have independently constrained the external convergence using weak lensing (Tihhonova et al. 2018, hereafter H0LiCOW VIII, and Tihhonova et al. submitted, hereafter H0LiCOW XI). The external convergence determined through weak lensing is consistent with our weighted number counts calculation.

3.5 Joint Inference

For our analysis, we make use of multiple datasets, denoted by $\bm{d_{\mathrm{img}}}$ for the HST and (if available) AO imaging data, $\bm{\Delta t}$ for the time delays, $\bm{\sigma}$ for the velocity dispersion of the lens galaxy, and $\bm{d_{\mathrm{LOS}}}$ for the properties of the LOS mass distribution determined from our photometric and spectroscopic data. We want to obtain the posterior PDF of the model parameters $\bm{\xi}$ given the data, $P(\bm{\xi}|\bm{d_{\rm img},\Delta t,\sigma,d_{\mathrm{LOS}},A})$. The vector $\bm{\xi}$ includes the lens model parameters $\bm{\nu}$, the cosmological parameters $\bm{\pi}$, and nuisance parameters representing the external convergence ($\kappa_{\rm ext}$) and anisotropy radius for the lens stellar velocity ellipsoid ($r_{\rm ani}$). ${\bm{A}}$ denotes a discrete set of assumptions about the form of the model, including the choices we have to make about the data modeling region, the set-up of the source reconstruction grid, the treatment of the different deflector mass distributions, etc. In general, ${\bm{A}}$ cannot be fully captured by continuous parameters. From Bayes’ theorem, we have that

[TABLE]

where ${P(\bm{d_{\rm img},\Delta t,\sigma,d_{\mathrm{LOS}}}}|\bm{\xi,A})$ is the joint likelihood function and $P(\bm{\xi}|\bm{A})$ is the prior PDF for the parameters given our assumptions. Since the data sets are independent, the likelihood can be separated,

[TABLE]

We can calculate the individual likelihoods separately and combine them as in Equation (6) to get the final posterior PDF for a given set of assumptions.

For any given lens model, we can vary the content of $\bm{A}$ and repeat the inference of ${\bm{\xi}}$. This can be important for checking the impact of various modeling choices and assumptions, but leaves us with the question of how to combine the results. Depending on the lens, we either combine the models with equal weight, or we can use the Bayesian Information Criterion (BIC) to weight the different models in our final inference (e.g., H0LiCOW IX, XII, Chen et al., 2019). This effectively combines our various assumptions $\bm{A}$ using the BIC so that we obtain $P(\bm{\xi}|\bm{d_{\rm img},\Delta t,\sigma,d_{\mathrm{LOS}}})$. We can further marginalize over the non-cosmological parameters ($\bm{\nu}$, $\kappa_{\rm ext}$, $r_{\rm ani}$) and obtain the posterior probability distribution of the cosmological parameters $\bm{\pi}$:

[TABLE]

In the lens modeling of systems with a single strong lens plane, the parameters associated with cosmology that enter directly into the lens modeling are the two lensing distances $D_{\Delta t}$ and $D_{\rm d}$. In the lens modeling of systems with multiple strong lens planes, we actually vary $H_{0}$, keeping other parameters fixed at $w=-1$, $\Omega_{\rm m}=0.3$, and $\Omega_{\Lambda}=0.7$. This assumes a fixed curvature of the expansion history of the Universe, but not the absolute scale (represented by $H_{0}$ or $D_{\Delta t}$). This is done because there is not a unique $D_{\Delta t}$ when accounting for multiple lens planes, but we convert this to an “effective" $D_{\Delta t}$ that is insensitive to assumptions of the cosmological model. Specifically, given the lens/quasar redshifts and $\bm{\pi}$ (i.e., $H_{0}$ and the other fixed cosmological parameters), we can compute the effective time-delay distance $D_{\Delta t}(\bm{\pi},z_{\rm d},z_{\rm s})$ to obtain the posterior probability distribution of $D_{\Delta t}$, $P(D_{\Delta t}|\bm{d_{\rm img},\Delta t,\sigma,d_{\mathrm{LOS}}})$. In summary, the single-lens plane models yield a joint constraints on $D_{\Delta t}$ and $D_{\rm d}$, whereas multi-lens plane models yield a constraint on the effective $D_{\Delta t}$.

3.6 Blind Analysis

After the development of the lens modeling and analysis methods that were first applied to B1608$+$656, the subsequent five lenses in H0LiCOW are analyzed blindly with respect to the cosmological quantities of interest (i.e., $D_{\Delta t}$, $D_{\rm d}$, $H_{0}$). Throughout the analyses, these values are blinded by subtracting the median of their PDF from the distribution. This allows us to view the shape of the distribution, their relative shifts, as well as covariances with other model parameters without ever seeing the absolute value. This is done to prevent confirmation bias and to remove the tendency for experimenters to stop analyzing systematic errors once they have achieved a result that agrees with a prior “expected" value. When the analysis of a particular H0LiCOW lens is finished and all team members have agreed to show the results, we unblind the relevant parameters and publish the result with no further changes to the analysis.

3.7 Distance Constraints and Error Budget for the Sample

We list the $D_{\Delta t}$ and (when available) $D_{\rm d}$ constraints for each individual lens in Table 2, along with corresponding references. All distances listed here are used in our cosmological inference. Specifically, for B1608$+$656, we use the analytic fit of $P(D_{\Delta t})$ given in Suyu et al. (2010) and of $P(D_{\rm d})$ given in Jee et al. (2019), and multiply these two PDFs since these two distances are uncorrelated for this system. For HE 0435$-$1223, RXJ1131$-$1231 and PG 1115$+$080, we use the resulting Monte Carlo Markov chains (MCMC) of $D_{\Delta t}$ and (if available) $D_{\rm d}$ from the analysis of Chen et al. (2019), which includes the previous HST constraints from Suyu et al. (2014) and H0LiCOW IV, as well as new AO data. For SDSS 1206$+$4332, we use the resulting MCMC chain of $D_{\Delta t}$ and $D_{\rm d}$ from H0LiCOW IX. For WFI2033$-$4723, we use the resulting MCMC chain of $D_{\Delta t}$ from H0LiCOW XII. We use a kernel density estimator to compute $P(D_{\Delta t},D_{\rm d})$ or $P(D_{\Delta t})$ from the chains, allowing us to account for any correlations between $D_{\Delta t}$ and $D_{\rm d}$ in $P(D_{\Delta t},D_{\rm d})$.

We estimate the approximate $D_{\Delta t}$ error budget for each of the lenses in Table 3. The contributions from the time delay measurement and LOS calculation are based on a Gaussian approximation. The remainder of the uncertainty is attributed to the lens model and other sources, which are difficult to disentangle. This breakdown shows that there is no single source of error that dominates the uncertainty from time-delay cosmography in general. Rather, it depends on characteristics of each particular lens that can be effectively random (modulo certain biases that make the lens more likely to be discovered, although such biases affect distance measurements at $\lesssim 1\%$ level; Collett & Cunnington, 2016), such as the image configuration (affects time delays and modeling), the mass/ellipticity of the lens galaxy (affects image separation/time delays, the lens model, and can be linked to the overdensity of the local environment), LOS structure (which is mostly uncorrelated outside of the local lens environment), and other factors. Moving forward, efforts will have to be made to tackle all of these sources of error to improve constraints from the overall sample rather than focusing on a single factor. Alternatively, with a large enough sample, one can pick out a small number of “golden lenses" that have characteristics that make them likely to have small uncertainties from each of the contributing sources of error, although one would have to be careful about potential biases in culling such a subsample.

4 Checking Consistency Among Lenses

We check that all our lenses can be combined without any loss of consistency by comparing their $D_{\Delta t}$ posteriors in the full cosmological parameter space and measuring the degree to which they overlap. We quantify the consistency by using the Bayes factor (or evidence ratio) $F$ in favor of a simultaneous fit of the lenses using a common set of cosmological parameters (e.g., Marshall et al., 2006; Suyu et al., 2013; Bonvin et al., 2017). When comparing data sets $\bf{d_{1}},...,\bf{d_{n}}$, we can either assume the hypothesis $\rm{H^{global}}$ that they can be represented using a common global set of cosmological parameters, or the hypothesis $\rm{H^{ind}}$ that at least one data set is better represented using a different set of cosmological parameters. We emphasize that the latter model would make sense if there is a systematic error present that leads to a vector offset in the inferred cosmological parameters. Parameterizing this offset vector with no additional information would take as many nuisance parameters as there are dimensions in the cosmological parameter space; assigning uninformative uniform priors to each of the offset components is equivalent to using a complete set of independent cosmological parameters for the outlier dataset.

We can compute the Bayes factor between any two lenses,

[TABLE]

where $L_{i}$ and $L_{j}$ are the likelihoods of data sets $\mathrm{{\bf d_{i}}}$ and $\mathrm{{\bf d_{j}}}$, respectively. If the six lenses have Bayes factors $F>1$ for every possible pairwise combination, it means that the lenses are statistically consistent with each other and we can proceed to combine their constraints.

In Table 4, we show that none of the 15 possible pairwise combinations of the six lens systems have a Bayes factor $F<1$. The minimal Bayes factor is obtained for the pair PG 1115$+$080 - SDSS 1206$+$4332 with $F=2.7$, still favoring the $\rm{H^{global}}$ hypothesis. We also test the hypothesis that one out of six systems is better represented in a different set of cosmological parameters than the five remaining lenses. The minimal Bayes factor is obtained for RXJ1131$-$1231 with $F=3.8$, again in favor of the $\rm{H^{global}}$ hypothesis, meaning that all lenses are a consistent realization of the same underlying set of cosmological parameters. We conclude that none of the six lenses is in disagreement with the cosmological parameters inferred from the five other systems. This is an important check of the consistency of our results. If our uncertainties were underestimated, we would not necessarily expect all of our lenses to give statistically consistent results.

5 Results of Cosmographic Analysis

We list the cosmological models considered in our analysis in Table 5. We distinguish between models where we use constraints from strong lenses alone (Sections 5.1 and 5.2) from those in which we combine our constraints with other probes via importance sampling (Section 5.3) or MCMC sampling (Section 5.4), even if the underlying cosmological model is the same. For the analysis using strong lenses only, we adopt uniform priors on the cosmological parameters with ranges indicated in Table 5.

5.1 Flat $\Lambda$CDM

Our baseline model is the flat $\Lambda$CDM cosmology with a uniform prior on $H_{0}$ in the range $[0,150]$ km s*-1* Mpc*-1* and a uniform prior on $\Omega_{\rm m}$ in the range $[0.05,0.5]$. In Figure 2, we show the marginalized constraints on $H_{0}$ from each of the individual H0LiCOW lenses along with the combined constraint from all six systems. We find $H_{0}=73.3_{-1.8}^{+1.7}~{}\mathrm{km~{}s^{-1}~{}Mpc^{-1}}$, a $2.4\%$ precision measurement. We show the median and 68% quantiles of the cosmological parameter distributions in Table 6. This estimate is higher than the Planck Collaboration et al. (2018b) CMB value ($H_{0}=67.4\pm 0.5~{}\mathrm{km~{}s^{-1}~{}Mpc^{-1}}$) by $3.1\sigma$ and in agreement with the latest SH0ES result ($H_{0}=74.03\pm 1.42~{}\mathrm{km~{}s^{-1}~{}Mpc^{-1}}$) from Riess et al. (2019).

Bonvin et al. (2017) noted that the first three H0LiCOW systems showed a trend of lower lens redshift systems having a larger inferred value of $H_{0}$, but could not conclude anything due to the small sample size. With a sample of six lenses, we see that this general trend still remains, as well as a trend of decreasing $H_{0}$ with increasing $D_{\Delta t}$. Even with six lenses, these correlations are not significant enough to conclude whether this is a real effect arising from some unknown systematic, a real physical effect related to cosmology, or just a statistical fluke (see Appendix A). To verify that the low lens redshift systems (RXJ1131$-$1231 and PG 1115$+$080) can safely be combined with the other four, we compute the Bayes factor between these two groups to be $F=1.9$, indicating that there is no statistical evidence that a different set of cosmological parameters is better representing the low redshift lenses. Nevertheless, the persistence of these trends is something to continue to examine as the sample of time-delay lenses increases in the future.

5.2 Extensions to Flat $\Lambda$CDM, Constraints from Time-Delay Cosmography Only

Given the current tension between determinations of $H_{0}$ from CMB observations and local probes, a possibility is that the underlying cosmology that describes our Universe is more complex than the standard flat $\Lambda$CDM model. Here, we present constraints from time-delay cosmography alone in some common single-parameter or two-parameter extensions to flat $\Lambda$CDM. The parameter constraints for the models we test here are given in Table 6.

5.2.1 Open $\Lambda$CDM

A simple modification to the flat $\Lambda$CDM cosmology is an open $\Lambda$CDM cosmology that allows for spatial curvature, $\Omega_{\rm k}\neq 0$. In this model, we have $\Omega_{\rm m}=1-\Omega_{\Lambda}-\Omega_{\rm k}$. We adopt uniform prior on $\Omega_{\rm k}$ in the range $[-0.5,0.5]$, $\Omega_{\Lambda}$ in the range [0,1], and require that $\Omega_{\rm m}>0$ . We still maintain the uniform prior on $H_{0}$ in the range $[0,150]$ km s*-1* Mpc*-1*.

The parameter constraints are given in Table 6. Figure 3 shows the marginalized constraint on $H_{0}$ in an open $\Lambda$CDM cosmology, which we find to be $H_{0}=74.4_{-2.3}^{+2.1}~{}\mathrm{km~{}s^{-1}~{}Mpc^{-1}}$. This is consistent with our flat $\Lambda$CDM constraint, although with a larger uncertainty. This constraint is still inconsistent with the Planck value, indicating that allowing for spatial curvature does not resolve the tension.

In Figure 4, we show a contour plot of the joint constraint on $H_{0}$ and $\Omega_{\rm k}$. The black contour is the constraint from strong lensing alone. We see that $\Omega_{\rm k}$ is very poorly constrained ($\Omega_{\rm k}=0.26_{-0.25}^{+0.17}$). This is not surprising, as the time-delay distance is only weakly sensitive to $\Omega_{\rm m}$ and $\Omega_{\Lambda}$, so we would expect a similar insensitivity to $\Omega_{\rm k}$. However, the fact that time-delay cosmography constrains $H_{0}$ very tightly indirectly imposes a tight constraint on curvature when combined with other probes.

5.2.2 Flat $w$CDM

We consider a flat $w$CDM cosmology in which the dark energy density is not a cosmological constant, but instead is time-dependent with an equation-of-state parameter $w$. We denote the dark energy density parameter as $\Omega_{\rm DE}=1-\Omega_{\rm m}$. The $w=-1$ case corresponds to flat $\Lambda$CDM with $\Omega_{\rm DE}=\Omega_{\Lambda}$. We adopt a uniform prior on $w$ in the range $[-2.5,0.5]$, keeping the same uniform priors on $H_{0}$ and $\Omega_{\rm m}$ as in the flat $\Lambda$CDM model.

We show the parameter constraints in Table 6. In Figure 5, we show the marginalized constraint on $H_{0}$ in this cosmology, which is $H_{0}=81.6_{-5.3}^{+4.9}~{}\mathrm{km~{}s^{-1}~{}Mpc^{-1}}$. The combined constraint on $H_{0}$ appears to be shifted to a higher value than most of the individual lenses, but this is due to the degeneracy between $H_{0}$ and $w$ and the resulting asymmetry in the PDF when projected from the higher-dimensional cosmological parameter space. Indeed, there is a region in the two-dimensional H_{0}$$-$$w plane in which all the individual distributions converge. That region has a higher probability density when performing the joint inference, which in turn drives the marginalized $H_{0}$ value higher.

In Figure 6, we show the joint distribution of $H_{0}$ and $w$. Lensing alone does not constrain $w$ particularly well ($w=-1.90_{-0.41}^{+0.56}$), although there is a degeneracy between $w$ and $H_{0}$, suggesting that combining our constraint with other probes may produce useful constraints.

5.2.3 Flat $w_{0}w_{a}$CDM

The flat $w$CDM cosmology has a time-varying dark energy component with an equation of state parameter $w$. In principle, $w$ itself could be changing with time. We consider a flat $w_{0}w_{a}$CDM cosmology in which the dark energy equation-of-state parameter $w$ is time-varying and parameterized as $w(z)=w_{0}+w_{a}z/(1+z)$ (Chevallier & Polarski, 2001; Linder, 2003). We adopt a uniform prior on $w_{0}$ in the range $[-2.5,0.5]$ and on $w_{a}$ in the range $[-2,2]$, keeping the same uniform priors on $H_{0}$ and $\Omega_{\rm m}$ as in the flat $\Lambda$CDM model.

We show the parameter constraints in Table 6. In Figure 7, we show the joint constraints on $H_{0}$, $w_{0}$, and $w_{a}$ from the six lenses in open black contours. Unsurprisingly, the lenses provide little constraint on $w_{0}$ and effectively no constraint on $w_{a}$, with the posterior PDF on $w_{a}$ spanning the entire prior range of $[-2,2]$. The resulting $H_{0}$ with a time-varying $w(z)$ is similarly high ($H_{0}=81.3_{-5.4}^{+5.1}~{}\mathrm{km~{}s^{-1}~{}Mpc^{-1}}$), as in the case of a flat $w$CDM cosmology, due to parameter degeneracies between $H_{0}$ and $w_{0}$.

5.3 Extensions to Flat $\Lambda$CDM, Combinations with CMB and BAO

The constraints from time-delay cosmography can be combined with other probes to provide joint constraints in the open $\Lambda$CDM flat $w$CDM, and flat $w_{0}w_{a}$CDM cosmologies, as well as a number of other extensions. Although time-delay cosmography is primarily sensitive to $H_{0}$ and only weakly dependent on other parameters, the degeneracies are such that strong lensing is highly complementary to other cosmological probes. Here, we combine strong lensing with CMB observations from Planck (Planck Collaboration et al., 2018b), sometimes including CMB weak lensing (CMBL) results, using the chains provided by the Planck team444http://pla.esac.esa.int/pla/#cosmology. We can also further combine with baryon acoustic oscillation (BAO) constraints at various redshifts (Beutler et al., 2011; Ross et al., 2015; Alam et al., 2017). All parameter constraints for the models investigated in this section are presented in Table 7.

5.3.1 Open $\Lambda$CDM

We test the open $\Lambda$CDM model presented in Section 5.2.1 with combined constraints from time-delay cosmography and Planck. In Figure 4, the grey contours show the constraints from Planck alone, while the orange contours show the combination of Planck with time-delay cosmography. We note that the H0LiCOW and Planck constraints are discrepant with each other. Nonetheless, we show the joint constraints here for completeness. The Planck data alone give a strong degeneracy between $H_{0}$ and $\Omega_{\rm k}$, but the addition of time-delay cosmography information provides important complementary information that constrains the Universe to be nearly flat ($\Omega_{\rm k}=-0.002_{-0.003}^{+0.003}$). This result is also consistent with the combination of Planck+CMBL, as well as Planck+BAO.

5.3.2 Flat $w$CDM

We test the flat $w$CDM model presented in Section 5.2.2 with combined constraints from time-delay cosmography and Planck. The grey contours in Figure 6 show the constraints from Planck alone, while the orange contours show the combination of Planck and time-delay cosmography. Although both Planck and H0LiCOW show degeneracies between the two parameters, the degeneracy directions are slightly different, allowing the combination of the two probes to constrain $w=1.36_{-0.09}^{+0.09}$. There is mild tension at the $\sim 2\sigma$ level between the Planck+H0LiCOW constraints discussed above, and Planck+BAO constraints.

5.3.3 Flat $\Lambda$CDM with variable neutrino species and/or masses

The standard model has an effective number of primordial neutrino species $\mathrm{N_{eff}}=3.046$, but additional relativistic particles in the early Universe prior to recombination could, in principle, add to this quantity. The sum of neutrino masses has not been precisely measured, but limits have been placed by a variety of experiments. In our analysis thus far, we have set $\sum\mathrm{m}_{\nu}=0.06$, the minimum mass allowed by neutrino oscillation experiments (Patrignani, 2016). We test cosmologies in which $\mathrm{N_{eff}}$ or $\sum\mathrm{m}_{\nu}$ are allowed to vary, as well as one in which both are allowed to vary. Time-delay cosmography alone cannot constrain either quantity, but combining it with other probes can help to break degeneracies.

Figure 8 shows a contour plot of $H_{0}$ and $\mathrm{N_{eff}}$ (when $\mathrm{N_{eff}}$ is allowed to vary). Figure 9 shows a contour plot of $H_{0}$ and $\sum\mathrm{m}_{\nu}$ (when $\sum\mathrm{m}_{\nu}$ is allowed to vary). In a model in which both quantities are allowed to vary, we show the parameter constraints in Figure 10. This model gives the combined constraints of $\mathrm{N_{eff}}=3.21_{-0.17}^{+0.16}$ and $\sum\mathrm{m}_{\nu}=0.03_{-0.03}^{+0.05}$, consistent with the standard model values. Although we see that the H0LiCOW and Planck constraints in these cosmologies are somewhat in tension, we provide the joint constraints here for completeness.

5.3.4 Flat $w_{0}w_{a}$CDM

We consider the same flat w_{0}$$w_{a}CDM cosmology as in Section 5.2.3, now including Planck and BAO as external data sets. The joint constraints on $H_{0}$, $w_{0}$, and $w_{a}$ are shown in Figure 7. The combination of H0LiCOW, Planck, and BAO constraints finds $w_{0}=-0.73_{-0.30}^{+0.29}$ and $w_{a}=-1.02_{-0.88}^{+0.80}$. Time-delay cosmography does not add much information due to the large degeneracies.

5.4 Calibrating Type Ia Supernovae Distances with Time-Delay Cosmography

The distance ladder method uses local distance indicators (e.g., Cepheid variables, TRGB) to calibrate the absolute distances to type Ia SNe. In principle, any absolute distance measurement can be used to anchor the distance scale to SNe. The inverse distance ladder method uses distances measured from baryon acoustic oscillations (BAO) in this way (e.g., Aubourg et al., 2015; Cuesta et al., 2015; Macaulay et al., 2019).

It is also possible to anchor SNe distances using either angular diameter distances ($D_{\rm d}$) to lens galaxies (e.g., Jee et al., 2019; Wojtak & Agnello, 2019) or $D_{\Delta t}$ (e.g., Collett et al., 2019; Liao et al., 2019; Taubenberger et al., 2019). This can be used as a complementary probe of $H_{0}$ or the sound horizon, $r_{s}$ (Arendse et al., 2019a, b), that is nearly insensitive to the assumed cosmology.

We follow the methodology of Taubenberger et al. (2019), using the combined $D_{\Delta t}$ and $D_{\rm d}$ measurements from the six lenses analyzed by H0LiCOW (Section 3.7) to anchor measurements of type Ia SNe from the “joint light-curve analysis" (JLA) sample of Betoule et al. (2014). For comparison, we also consider the Pantheon type Ia SNe sample (Scolnic et al., 2018). We use the MontePython v3.1 MCMC sampling package (Audren et al., 2013; Brinckmann & Lesgourgues, 2018), its associated CLASS code (Lesgourgues, 2011), and the JLA and Pantheon SNe samples implemented within MontePython, to combine with our lensing distance measurements and sample cosmological parameters.

We consider six cosmological models that are listed in part of Table 5: $\Lambda$CDM, o$\Lambda$CDM, $w$CDM, $w_{0}w_{a}$CDM, o$w$CDM and o$w_{0}w_{a}$CDM. For the first four cosmological models, we adopt the same uniform prior ranges for the cosmological parameters as in the top part of Table 5 (i.e., same priors as U$\Lambda$CDM, Uo$\Lambda$CDM, U$w$CDM and U$w_{0}w_{a}$CDM, respectively). For o$w$CDM and o$w_{0}w_{a}$CDM, we adopt the same priors for $w$ and $\{w_{0},w_{a}\}$ as those in U$w$CDM and U$w_{0}w_{a}$CDM, respectively, with the remaining parameters ($H_{0}$, $\Omega_{\rm m}$, $\Omega_{\rm k}$ and $\Omega_{\rm DE}$) having the same priors as those in o$\Lambda$CDM.

Table 8 shows the results for the six cosmologies tested, for the JLA sample (top) and the Pantheon sample (bottom). The median values of $H_{0}$ from the JLA sample are within $1.5\,\rm km\,s^{-1}\,Mpc^{-1}$ of those of the Pantheon sample. In comparison to the $H_{0}$ constraints in Section 5.2, particularly for the $w$CDM and $w_{0}w_{a}$CDM models where $H_{0}$ is highly degenerate with $w$, $w_{0}$ and $w_{a}$, the $H_{0}$ from the lenses and SNe in Table 8 are less sensitive to cosmological models, as shown in Figure 11. For the six cosmological models probed in Table 8, the median $H_{0}$ values range from $\sim$73$-$78 $\mathrm{km~{}s^{-1}~{}Mpc^{-1}}$, irrespective of the SNe sample. The tension with Planck in flat $\Lambda$CDM is still $>3\sigma$, similar to our result from time-delay cosmography alone.

We also repeat the JLA analysis using only the marginalized $P(D_{\Delta t})$ from the six lenses as constraints, i.e., omitting the information from $D_{\rm d}$. This allows us to assess the information content on $H_{0}$ from the lensing distances and the added value of measuring $D_{\rm d}$ in addition to $D_{\Delta t}$. Using only $P(D_{\Delta t})$, we find that the $H_{0}$ values are very similar to values in Table 8 (within $0.5\,\rm km\,s^{-1}\,Mpc^{-1}$), and the uncertainties are also only slightly larger (by at most $0.5\,\rm km\,s^{-1}\,Mpc^{-1}$). Therefore, most of the cosmological information, particularly $H_{0}$, is encapsulated in our $D_{\Delta t}$ measurements that are substantially more precise than $D_{\rm d}$ measurements. Nonetheless, future spatially resolved kinematics of the lens galaxy could help tighten the constraints on $D_{\rm d}$, providing more leverage on cosmological parameters (Shajib et al., 2018; Yıldırım et al., 2019).

6 Tension between early-Universe and late-Universe probes of $H_{0}$

As the tension between different probes of $H_{0}$ began to emerge in recent years, a natural direction to look toward in order to resolve this apparent discrepancy has been to examine potential sources of systematic error in the various methods. In addition to exploring possible systematics in the Planck analysis and those based on type Ia SNe calibrated by the distance ladder, having multiple independent probes has proven to be a crucial step in checking these results.

The latest H0LiCOW results presented here, analyzed blindly with respect to cosmological parameters, are the most precise constraints on $H_{0}$ from time-delay cosmography to date, and are independent of both CMB probes (i.e., Planck) and other late-Universe probes such as SH0ES. Our results for a flat $\Lambda$CDM cosmology are in $3.1\sigma$ tension with Planck. In combination with the latest SH0ES result (Riess et al., 2019), we find a $5.3\sigma$ tension between late-Universe determinations of $H_{0}$ and Planck (Figure 12). Other independent methods anchored in the early Universe, such as the analysis of Abbott et al. (2018b) using a combination of clustering and weak lensing, BAO, and big bang nucleosynthesis (BBN), give similar results to Planck. Although the Abbott et al. (2018b) method, as well as the inverse distance ladder, combine measurements from both the early and late Universe, the inferred $H_{0}$ is ultimately set by the sound horizon at recombination, which comes from early-Universe physics.

Given the tension shown here, it is becoming difficult to reconcile the $H_{0}$ discrepancy by appealing to systematic effects. While systematics, especially “unknown unknowns", still cannot be entirely ruled out and should continue to be explored, recent work has only heightened the tension. There also appears to be a growing dichotomy when the different $H_{0}$ probes are split into those anchored by the early-Universe (i.e., CMB), which favor a lower $H_{0}$, and those based on late-Universe probes, which favor a higher $H_{0}$ (e.g., Verde et al., 2019).

As this tension between early-Universe and late-Universe probes continues to grow, we must examine potential alternatives to the standard flat $\Lambda$CDM model. This would be a major paradigm shift in modern cosmology, requiring new physics to consistently explain all of the observational data. We have explored some possible extensions to flat $\Lambda$CDM in Section 5, including spatial curvature, time-varying dark energy (e.g., Di Valentino et al., 2018), and modified neutrino physics such as sterile neutrinos (e.g., Wyman et al., 2014; Gelmini et al., 2019) or self-interacting neutrinos at early times (e.g., Kreisch et al., 2019). Other possible new physics to resolve the discrepancy include an early dark energy component to the Universe that later decays (e.g., Agrawal et al., 2019; Alexander & McDonough, 2019; Aylor et al., 2019; Lin et al., 2019; Poulin et al., 2019), primordial non-Gaussianity (e.g., Adhikari & Huterer, 2019), decaying dark matter (e.g., Pandey et al., 2019; Vattis et al., 2019), and fifth forces (e.g., Desmond et al., 2019).

7 Summary

We have combined time-delay distances and angular diameter distances from six lensed quasars in the H0LiCOW sample to achieve the highest-precision probe of $H_{0}$ to date from strong lensing time delays. Five of the six lenses are analyzed blindly with respect to the cosmological parameters of interest. Our main results are as follows:

• •

We find $H_{0}=73.3_{-1.8}^{+1.7}~{}\mathrm{km~{}s^{-1}~{}Mpc^{-1}}$ for a flat $\Lambda$CDM cosmology, which is a measurement to a precision of $2.4\%$. This result is in agreement with the latest results from measurements of type Ia SNe calibrated by the distance ladder (Riess et al., 2019) and in $3.1\sigma$ tension with Planck CMB measurements (Planck Collaboration et al., 2018b).

• •

Our constraint on $H_{0}$ in flat $\Lambda$CDM is completely independent of and complementary to the latest results from the SH0ES collaboration, so these two measurements can be combined into a late-Universe constraint on $H_{0}$. Together, these are in tension with the best early-Universe (i.e., CMB) determination of $H_{0}$ from Planck at a significance of $5.3\sigma$.

• •

We check that the lenses in our sample are statistically consistent with one another by computing Bayes factors between their $H_{0}$ PDFs. We find that all six lenses are pairwise consistent (i.e., $F>1$), indicating that we are not underestimating our uncertainties and are able to control systematic effects in our analysis.

• •

We compute parameter constraints for cosmologies beyond flat $\Lambda$CDM. In an open $\Lambda$CDM cosmology, we find $\Omega_{\rm k}=0.26_{-0.25}^{+0.17}$ and $H_{0}=74.4_{-2.3}^{+2.1}~{}\mathrm{km~{}s^{-1}~{}Mpc^{-1}}$, which is still in tension with Planck, suggesting that allowing for spatial curvature cannot resolve the discrepancy. In a flat $w$CDM cosmology, we find $H_{0}=81.6_{-5.3}^{+4.9}~{}\mathrm{km~{}s^{-1}~{}Mpc^{-1}}$ and $w=-1.90_{-0.41}^{+0.56}$. In a flat $w_{0}w_{a}$CDM cosmology, we find $H_{0}=81.3_{-5.4}^{+5.1}~{}\mathrm{km~{}s^{-1}~{}Mpc^{-1}}$, but are unable to place meaningful constraints on $w_{0}$ and $w_{a}$.

• •

We combine our constraints with Planck, including CMB weak lensing and BAO constraints. Although time-delay cosmography is primarily sensitive to $H_{0}$, with only a weak dependence on other cosmological parameters, the constraints are highly complementary to other probes such as Planck, CMB weak lensing, and BAO. We test the open $\Lambda$CDM and $w$CDM cosmologies, as well as cosmologies with variable effective neutrino species and/or sum of neutrino masses, and a $w$CDM cosmology with a time-varying $w$. The full parameter constraints for these models when combining H0LiCOW and Planck are given in Table 7.

• •

We use the distance measurements from time-delay cosmography to calibrate the distance scale of type Ia SNe from the JLA and Pantheon samples. This provides a probe of $H_{0}$ that is less dependent on the assumed cosmological model, in comparison to the constraints from lenses alone. We find median $H_{0}$ values ranging from $\sim$73$-$78 $\mathrm{km~{}s^{-1}~{}Mpc^{-1}}$ for a range of cosmologies. The tension with Planck for a flat $\Lambda$CDM cosmology is still $>3\sigma$, similar to the result from time-delay cosmography alone.

Despite efforts to explore and reduce systematic errors in the various methods, the growing tension between early and late-Universe probes of $H_{0}$ has only continued to heighten. If unresolved, this tension may force the rejection of the flat $\Lambda$CDM model in favor of new physics, which would dramatically alter our understanding of the Universe.

While considering the possibility of new physics, we are also continuing to improve the constraints from time-delay cosmography. The current sample of six H0LiCOW systems is already the best-studied sample to date, and a number of additional lenses are being observed with high-resolution imaging (e.g., Shajib et al., 2019) and monitored by COSMOGRAIL. Moving into the future, many new lensed quasars are being discovered in large imaging surveys (e.g., Agnello et al., 2015, 2018a, 2018b; Anguita et al., 2018; Lemon et al., 2018; Lemon et al., 2019; Treu et al., 2018). A sample of $\sim 40$ lenses is needed to constrain $H_{0}$ to the $\sim 1\%$ level (Jee et al., 2016; Shajib et al., 2019), which will be attainable in the near future.

Acknowledgements

We thank the referee, whose suggestions were helpful in improving the clarity of this paper. We thank Chiara Spiniello, Malte Tewes, and Akın Yıldırım for their contributions to the H0LiCOW project. We thank Lodovico Coccato and Johan Richard for their help with the MUSE data used as a part of this project. We thank Aleksi Halkola for support with the Glee lens modeling code. We thank all the observers at the Euler, SMARTS, Mercator and Maidanak telescopes who participated in the queue-mode observations. H0LiCOW and COSMOGRAIL are made possible thanks to the continuous work of all observers and technical staff obtaining the monitoring observations, in particular at the Swiss Euler telescope at La Silla Observatory. The Euler telescope is supported by the Swiss National Science Foundation. This work was supported by World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan. K.C.W. is supported in part by an EACOA Fellowship awarded by the East Asia Core Observatories Association, which consists of the Academia Sinica Institute of Astronomy and Astrophysics, the National Astronomical Observatory of Japan, the National Astronomical Observatories of the Chinese Academy of Sciences, and the Korea Astronomy and Space Science Institute. S.H.S. thanks the Max Planck Society for support through the Max Planck Research Group. S.H.S. and S.T. are supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 771776). 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 No. AST-1715611. M.M, D.S., V.B., F.C., and O.T. are 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 programme (COSMICLENS: grant agreement No 787886). T.T. acknowledges support by the Packard Foundation through a Packard Research fellowship and by the National Science Foundation through NSF grant AST-1714953. A.A. was supported by a grant from VILLUM FONDEN (project number 16599). This project is partially funded by the Danish council for independent research under the project “Fundamentals of Dark Matter Structures”, DFF–6108-00470. S.H. acknowledges support by the DFG cluster of excellence ‘Origin and Structure of the Universe’ (www.universe-cluster.de). A.J.S. acknowledges support by NASA through Space Telescope Science Institute grant HST-GO-15320. L.V.E.K. is partly supported through an NWO-VICI grant (project number 639.043.308). Based on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with programs HST-GO-9375, HST-GO-9744, HST-GO-10158, HST-GO-12889, and HST-14254. Support for programs HST-GO-10158 HST-GO-12889 HST-14254 was provided to members of our team by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. This research made use of Astropy,555http://www.astropy.org a community-developed core Python package for Astronomy (Astropy Collaboration et al., 2013, 2018). This research made use of Matplotlib, a 2D graphics package used for Python (Hunter, 2007). This research made use of emcee, a Python implementation of an affine invariant MCMC ensemble sampler (Foreman-Mackey et al., 2013).

Appendix A Quantifying the trend of $H_{0}$ with $z_{\rm d}$ and $D_{\Delta t}$

There is an apparent trend of decreasing $H_{0}$ inferred from the individual lenses as a function of increasing lens redshift, which is shown in the left panel of Figure 13. The significance of this correlation can be assessed by testing against the null hypothesis, in which the measured $H_{0}$ values are uncorrelated with the lens redshift. To do so, we draw sets of six mock $H_{0}$ values, using each lens’ own uncertainty probability distribution centered around the median joint inference obtained in flat $\Lambda$CDM ($H_{0}=73.3_{-1.8}^{+1.7}~{}\mathrm{km~{}s^{-1}~{}Mpc^{-1}}$). We then fit a linear regression through each mock set. Associating a weight to the mock value from each lens is done according to the following scheme: we first rescale the uncertainties’ probability distributions so that their maximum value equal one. Next, we compute the area under each rescaled distribution, then rescale the areas by their median. Finally, we take the inverse square of the rescaled areas as weight for each mock measurement. The slope of the regression is taken as our measurement of the correlation. We create $10^{5}$ sets of mock values, for which the distribution of the measured slopes is centered around zero, as expected for a null hypothesis. We fit the same kind of linear regression through the data, for which we find a negative slope. We find that the slope of the data falls 1.9$\sigma$ away from the mock slope distribution. We also observe a correlation between $D_{\Delta t}$ and $H_{0}$ (Figure 13, right panel) that deviates from the null hypothesis at a similar significance level of $1.8\sigma$.