MNRAS 000, 1–?? (2019)

Preprint 12 July 2019

Compiled using MNRAS LATEX style file v3.0

arXiv:1907.04869v1 [astro-ph.CO] 10 Jul 2019

H0LiCOW XIII. A 2.4% measurement of H0 from lensed
quasars: 5.3σ tension between early and late-Universe
probesKenneth C. Wong,1,2? Sherry H. Suyu,3,4,5 Geoff C.-F. Chen,6 Cristian E. Rusu,2,7,6
Martin Millon,8 Dominique Sluse,9 Vivien Bonvin,8 Christopher D. Fassnacht,6Stefan Taubenberger,3 Matthew W. Auger,10 Simon Birrer,11 James H. H. Chan,8Frederic Courbin,8 Stefan Hilbert,12,13 Olga Tihhonova,8 Tommaso Treu,11 Adriano Agnello,14 Xuheng Ding,11 Inh Jee,3 Eiichiro Komatsu,3,1 Anowar J. Shajib,11Alessandro Sonnenfeld,15 Roger D. Blandford,16 Léon V. E. Koopmans,17 Philip J.
Marshall,16 and Georges Meylan8
1 Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan
2 National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan
3 Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, 85748 Garching, Germany
4 Physik-Department, Technische Universität München, James-Franck-Straße 1, 85748 Garching, Germany
5 Academia Sinica Institute of Astronomy and Astrophysics (ASIAA), 11F of ASMAB, No.1, Section 4, Roosevelt Road, Taipei 10617,

Taiwan6 Department of Physics, University of California, Davis, CA 95616, USA
7 Subaru Telescope, National Astronomical Observatory of Japan, 650 N Aohoku Pl, Hilo, HI 96720, USA
8 Institute of Physics, Laboratory of Astrophysics, Ecole Polytechnique Fédérale de Lausanne (EPFL), Observatoire de Sauverny, 1290
Versoix, Switzerland
9 STAR Institute, Quartier Agora - Allée du six Août, 19c B-4000 Liège, Belgium
10 Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK
11 Department of Physics and Astronomy, University of California, Los Angeles, CA 90095, USA
12 Exzellenzcluster Universe, Boltzmannstr. 2, D-85748 Garching, Germany
13 Ludwig-Maximilians-Universität, Universitäts-Sternwarte, Scheinerstr. 1, D-81679 München, Germany
14 DARK, Niels-Bohr Institute, Lyngbyvej 2, 2100 Copenhagen, Denmark
15 Leiden Observatory, Leiden University, Niels Bohrweg 2, 2333 CA Leiden, the Netherlands
16 Kavli Institute for Particle Astrophysics and Cosmology, Stanford University, 452 Lomita Mall, Stanford, CA 94035, USA
17 Kapteyn Astronomical Institute, University of Groningen, PO Box 800, NL-9700 AV Groningen, The Netherlands

Accepted XXX. Received YYY; in original form ZZZABSTRACT© 2019 The Authors

We present a measurement of the Hubble constant (H0 ) and other cosmological parameters from a joint analysis of six gravitationally lensed quasars with measured timedelays. All lenses except the first are analyzed blindly with respect to the cosmologi−1
cal parameters. In a flat ΛCDM cosmology, we find H0 = 73.3+1.7
Mpc−1 , a−1.8 km s
2.4% precision measurement, in agreement with local measurements of H0 from typeIa supernovae calibrated by the distance ladder, but in 3.1σ 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σ tension with Planck CMB
determinations of H0 in flat ΛCDM. We compute Bayes factors to verify that alllenses give statistically consistent results, showing that we are not underestimating
our uncertainties and are able to control our systematics. We explore extensions toflat Λ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.
Using the distance constraints from time-delay cosmography to anchor the type Ia
supernova distance scale, we reduce the sensitivity of our H0 inference to cosmologicalmodel assumptions. For six different cosmological models, our combined inference on
H0 ranges from ∼ 73–78 km s−1 Mpc−1 , which is consistent with the local distanceladder constraints.
Key words: cosmology: observations − cosmology: cosmological parameters − distance scale − gravitational lensing: strong

2
1

K. C. Wong et al.
INTRODUCTION

The flat Λ cold dark matter (Λ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 ΛCDM assumption by introducing additional complexities, such as non-zero curvature, an equation of state parameter w 6= −1, or a time-varying w, leads
to much weaker constraints and large degeneracies among
the various parameters. In particular, many parameters become degenerate with the Hubble constant, H0 , which sets
the present-day expansion rate of the universe. H0 cannot
be constrained directly from CMB observations, but must
be inferred by first assuming a cosmological model. In thiscontext, measuring H0 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 H0
is through observations of type Ia supernovae (SNe). TypeIa SNe are “standardizable candles” in that their luminosities, and thus their absolute distances, can be determinedby the evolution of their light curves, and therefore can beused to infer H0 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, 2019), in which
parallax measurements are used to determine distances tonearby Cepheid variable stars (which have a known periodluminosity relation) and are in turn used to determine distances to type Ia SNe in the Hubble flow.
Recent determinations of H0 from the Supernovae, H0 ,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 ΛCDM
model (e.g., Bernal et al. 2016; Freedman 2017). The latest
SH0ES result finds H0 = 74.03 ± 1.42 km s−1 Mpc−1 (Riess
et al. 2019), which differs from the Planck flat ΛCDM result1 of H0 = 67.4 ± 0.5 km s−1 Mpc−1 (Planck Collaboration et al. 2018b) by 4.4σ. 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; Riess et al. 2019). Furthermore, independent determinations of H0 using the “inverse distanceladder” 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 physicalscale of the sound horizon (e.g., Aylor et al. 2019; Macaulay
et al. 2019; Arendse et al. 2019). 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;
1

Baseline ΛCDM chains with baseline likelihoods (based on
plikHM TTTEEE lowl lowE)

Kozmanyan et al. 2019), the Balmer line L − σ relation of
HII galaxies (e.g., Melnick et al. 2000; Chávez et al. 2012;
González-Morán et al. 2019), extragalactic background lightattenuation (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 H0 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
ΛCDM model and indicate new physics that must be incorporated into our understanding of cosmology.
Gravitational lensing offers an independent method ofdetermining H0 . 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 takedifferent paths through space-time at the different image positions. Because these paths have different lengths and passthrough different gravitational potentials, light rays emitted
from the source at the same time will arrive at the observerat 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∆t , and depends on the mass distribution in the lensing object, the mass distribution along theline of sight (LOS), and cosmological parameters. D∆t is primarily sensitive to H0 , 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 andthe 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 in the local Universe (i.e., a “Hubble bubble”).
This method of using gravitational lens time delays tomeasure H0 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 H0 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 H0 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 fromMNRAS 000, 1–?? (2019)

H0LiCOW XIII: A 2.4% measurement of H0 from lensed quasars
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, κext , which
can bias the inferred D∆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 H0 Lenses in COSMOGRAIL’s Wellspring(H0LiCOW) collaboration (Suyu et al. 2017, hereafter
H0LiCOW I) has provided the strongest constraints on H0
to date from time-delay cosmography. Our most recent measurements had constrained H0 to 3.0% precision for a flat
Λ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 ofthe lens galaxy to measure its velocity dispersion (e.g., Sluseet al. 2019, hereafter H0LiCOW XI), and deep wide-field
spectroscopy and imaging to characterize the LOS in thesesystems (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 H0 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 (Chenet al. 2019) and WFI2033−4723 (Rusu et al. 2019, hereafter
H0LiCOW XII) – to provide the tightest H0 constraints todate from time-delay cosmography.
This paper is organized as follows. In Section 2, we summarize the theory behind using time-delay cosmography toinfer the time-delay distance, which is inversely proportional
to H0 . In Section 3, we present our lens sample and describe
how our data and analysis methods allow us to constrain H0
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 ΛCDM and more complex cosmologies
in Section 5. In Section 6, we discuss the tension between
early-Universe and late-Universe probes of H0 . We summarize our findings in Section 7. Throughout this paper, allmagnitudes given are on the AB system. All parameter conMNRAS 000, 1–?? (2019)3

straints given are medians and 16th and 84th percentiles
unless otherwise stated.

2

TIME-DELAY COSMOGRAPHY

In this section, we summarize the theoretical backgroundof time-delay cosmography and how to infer H0 . We referreaders to recent reviews (e.g., Treu & Marshall 2016; Suyu
et al. 2018) for more details.
When light rays from a background source are deflectedby an intervening lensing mass, the light travel time from thesource to the observer depends on both their path lengthand the gravitational potential they traverse. For a single
lens plane, the excess time delay of an image at an angular
position θ = (θ1 , θ2 ) with corresponding source position β =(β1 , β2 ) relative to the case of no lensing isD∆t (θ − β)2t(θ, β) =− ψ(θ) ,(1)c2
where D∆t is the time-delay distance and ψ(θ) is the lenspotential. The time-delay distance (Refsdal 1964; Schneider
et al. 1992; Suyu et al. 2010) is defined asD∆t ≡ (1 + zd )

Dd Ds
,
Dds

(2)

where zd is the lens redshift, Dd is the angular diameter
distance to the lens, Ds is the angular diameter distance tothe source, and Dds is the angular diameter distance between
the lens and the source. D∆t has units of distance and isinversely proportional to H0 , with weak dependence on othercosmological parameters.
If the alignment between the background source andthe foreground lens is close enough, multiple images of thesame background source are formed. Light rays reaching theobserver will have different excess time delays depending on
which image they are observed at. The time delay between
two images of such a lens, ∆tij , is the difference of theirexcess time delays,(θj − β)2D∆t (θi − β)2∆tij =− ψ(θi ) −+ ψ(θj ) ,
c
2
2
(3)
where θi and θj are the positions of images i and j, respectively, in the image plane. If the source is variable on shorttimescales, it is possible to monitor the fluxes of the images
and measure the time delay, ∆tij , between them (e.g., Vanderriest et al. 1989; Schechter et al. 1997; Fassnacht et al.
1999, 2002; Kochanek et al. 2006; Courbin et al. 2011). The
lens potentials at the image positions, ψ(θi ) and ψ(θj ), as
well as at the source position, β, can be determined from amass model of the system. With a measurement of ∆tij andan accurate lens model to determine ψ(θ), it is possible todetermine D∆t . By further assuming a cosmological model,D∆t can be converted into an inference on H0 .
If there are multiple lenses at different redshifts between
the source and the observer, the observed time delays depend on various combinations of the angular diameter distances among the observer, the multiple lens planes, and thesource. In this case, the image positions are described by themultiplane lens equation (e.g., Blandford & Narayan 1986;

4

K. C. Wong et al.

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 uniquetime-delay distance. However, it is often the case that themass in a single lens plane dominates the lensing effect, andthe observed time delays are primarily sensitive to the timedelay distance (Equation 2) with the deflector redshift asthat of the primary lens plane. This is the case for all of thelenses in the H0LiCOW sample (Section 3.1). The results forany individual system can thus be interpreted as a constraint
on D∆t (zd , zs ), which we refer to as the “effective time-delay
distance”. Hereafter, D∆t refers to the effective time-delay
distance (for applicable systems) unless otherwise indicated.
In addition to mass that is explicitly included in thelens model, all other mass along the LOS between the observer and the source contributes to the lens potential thatthe light rays traverse. This causes additional focusing anddefocusing of the rays and can affect the observed time delays (e.g., Seljak 1994). If unaccounted for, this can lead tobiased inferences of D∆t . If the effects of the perturbing LOS
masses are small, they can be approximated by an externalconvergence in the lens plane κext (e.g., Keeton 2003; McCully et al. 2014). The true D∆t is related to the time-delay
distance inferred from the lens model and measured timemodel
, by the relation
delays, D∆t

D∆t =

model
D∆t
.
1 − κext(4)κext is defined such that its average value across the sky
is zero. In principle, if lenses are randomly distributed, theeffect of κext should average out over a sufficiently large sample. However, the cross section for strong lensing scales as
σ 4 , where σ 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, lensgalaxies generally lie in overdense environments and lines ofsight relative to typical fields (e.g., Treu et al. 2009; Fassnacht et al. 2011; Wong et al. 2018), meaning that κext willlead to a bias on D∆t and needs to be corrected for. κext
cannot, in general, be constrained from the lens model dueto 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 thesource plane can modify the product of the time delays andH0 but leave other observables unchanged. κext must insteadbe estimated through other methods, such as studies of thelens environment or the use of lens stellar kinematics (asnoted in the previous section).
With kinematic information on the lens galaxy, it ispossible to determine the angular diameter distance to thelens, Dd , independent of κext (Paraficz & Hjorth 2009; Jeeet al. 2015). Although the constraints from Dd are generally weaker than those from D∆t , it can break degeneracies
among cosmological parameters, particularly for models beyond flat ΛCDM. In particular, it can break the degeneracy between curvature (Ωk ) and the time-varying equationof state parameter of dark energy (w) (Jee et al. 2016). The
combination of lensing, time delays, and lens kinematic datathus provides a joint constraint on D∆t and Dd in casesof 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 cosmologicalparameter 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 highresolution 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 κext in these systems from awide-field imaging and spectroscopic campaign, as detailedin Section 3.4.
The original H0LiCOW sample (see H0LiCOW I) consists of five lenses, and it was later decided to expand thesample to include four additional systems with HST imaging(PID:14254, PI: T. Treu), placing an emphasis on doubleimage lens systems (doubles), which yield fewer constraints
per system but are more abundant on the sky. Of thefive 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 quadruplyimaged, providing additional constraints for lens modeling. PG 1115+080 was observed with AO imaging fromKeck/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 theoriginal sample, plus SDSS 1206+4332 and PG 1115+080.
These six lenses span a range of lens and source redshifts,
as well as a variety of configurations (e.g., double, cross,fold, cusp). Having a sample that spanned a range in theseparameters was a consideration in the selection of which
lenses to analyze first, as there may be various systematics
that depend on such factors, and we want to account forthem in our analysis (see Ding et al. 2018, who attempt toaddress these issues based on simulated data).

3.2Time 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 and2m-size telescopes. Several seasons of monitoring are neededin order to disentangle the variations due to microlensing inwhich brightening or dimming of the quasar images by stars
MNRAS 000, 1–?? (2019)

H0LiCOW XIII: A 2.4% measurement of H0 from lensed quasars

5

Table 1. Lenses in the H0LiCOW sample used in this paper.
Lens nameα (J2000)δ (J2000)zd

zs

HST / AO data

B1608+656a

16:09:13.96+65:32:29.0

0.6304a1.394b

HST−12:31:57.0

0.295c0.654dHST + AO

−12:17:14.4

0.4546f,g

1.693h

HST + AO

1.789i

HST

RXJ1131−1231c

11:31:51.6

HE 0435−1223e

12:06:29.65+43:32:17.6

0.745jWFI2033−4723k

20:33:41.9−47:23:43.4

0.6575l

1.662hHST

PG 1115+080m

11:18:16.899+7:45:58.502

0.311n1.722m

HST + AOa Myers et al. (1995);
f Morgan et al. (2005);
j Agnello

04:38:14.9

SDSS 1206+4332i

et al. (2016);

1"

b Fassnacht et al. (1996);

c Sluse et al. (2003);

g Eigenbrod et al. (2006);
k Morgan

et al. (2004);

l Sluse1"(a) B1608+656(b) RXJ1131−1231(c) HE 0435−1223(d) SDSS 1206+4332

1"1"(e) WFI2033−4723(f) PG 1115+080

Figure 1. Multicolor images of the six lensed quasars used in
our analysis. The images are created using two or three imaging bands in the optical and near-infrared from HST and/or
ground-based AO data. North is up and east is to the left.
Images for B1608+656, RXJ1131−1231, HE 0435−1223, andWFI2033−4723 are from H0LiCOW I.

et al. (2019);

e Wisotzki et al. (2002);

i Oguri et al. (2005);
m Weymann

et al. (1980);

n Tonry(1998)

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 findthe 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
PyCS2 , 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 ∼ 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σ). The time delay of SDSS 1206+4332 wasalso measured with PyCS; Eulaers et al. (2013) obtained atime delay between the two multiple images of ∆tAB =111.3 ± 3, 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 timedelay (Bonvin et al. 2017).
Recently, Courbin et al. (2018) demonstrated that ahigh-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 fromthe microlensing signal in a single season. High-cadence data
were used for WFI2033−4723 and PG 1115+080 to measuretime delays at a few percent precision in one season. These
results are in agreement with the time delays measured fromdecade-long COSMOGRAIL light curves and are combinedin the final estimate (Bonvin et al. 2018, 2019).

2

MNRAS 000, 1–?? (2019)d Sluse et al. (2007);

h Sluse et al. (2012);

Available at http://www.cosmograil.org

6

K. C. Wong et al.

The remaining lens of the sample, B1608+656 was monitored by Fassnacht et al. (1999, 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 timedelays 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 thedisk can then be magnified by the microlenses differently ineach of the multiple images. This reweighting of the delayedemission across the accretion disk could lead to a change inthe measured time delay. As the microlensing changes withtime, this could lead to a variation in the measured time delays from season to season. This complicating factor is now
quantified (Bonvin et al. 2018, 2019; Chen et al. 2018) inthe latest H0LiCOW lens models (Birrer et al. 2019; Chenet 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 lensgalaxy light as parameterized profiles and fit the lensed
quasar image positions and surface brightness distributionof 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 onlens modeling code systematics.
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 tothe 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, westarted from the power-law model and performed a pixelated lens potential reconstruction to allow for flexibility,
finding small (∼ 2%) potential corrections and thus validating the use of power-law model. Galaxies along the LOS
that are deemed to be significant perturbers are includedin the model (Section 3.4) through the full multiplane lensequation, and we also include an external shear in the mainlens plane. When available, we use the measured velocity
dispersions of the lens galaxy and significant perturbers asadditional constraints.
To account for various systematic effects arising frommodeling choices in areas such as the lens parameterization,
the source reconstruction, the weighting of the pixels in theimage plane, etc., we run a variety of models where we varythese choices and combine them in our final inference. In
our initial analyses of the first three H0LiCOW lenses (Suyuet al. 2010, 2013, 2014, H0LiCOW IV), we marginalize overthese (discrete) modeling choices in deriving the posterior

probability density function (PDF) of D∆t . For models that
fit equally well to the data within their modeling uncertainties, we conservatively add their posterior distributions ofD∆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), weweight 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 break degeneracies in the massmodeling. They are also used to estimate the angular diameter distances to 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 lensplane (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 approximatedby 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 beapproximated by a κ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 canbe 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 projectedclose to the lens and include them in our lens model if theircalculated 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 XI). 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 sampleof 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 brightgalaxies as far as ∼ 100 away from the lens systems, allowing us to further quantify the properties of larger structures such as galaxy groups and clusters (e.g., H0LiCOW
XI). Since the analysis of HE 0435−1223 (H0LiCOW II), a
group-finding algorithm is used to identify and characterize potential group properties (e.g., LOS velocity dispersion
and centroid). 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 κext due to LOS
structure, we use a (weighted) galaxy number counts techMNRAS 000, 1–?? (2019)

H0LiCOW XIII: A 2.4% measurement of H0 from lensed quasars
nique (e.g., Greene et al. 2013; Suyu et al. 2013, H0LiCOW
III, IX). We count galaxies projected within a fixed distanceof the lens and above some flux threshold, weighted by various quantities. We then compare these number counts to
those measured analogously along random lines of sight in acontrol survey, to determine the relative over-/under-densityof the lens field. Finally, we use the Millennium simulation(Springel et al. 2005) to identify lines of sight that have asimilar relative number count densities and build a PDF ofκ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, widefield imaging campaign to gather data of sufficient depth andspatial coverage to enable this analysis. Our data typicallyconsist of multi-band ultraviolet to infrared observations ingood seeing conditions, which we use to perform a galaxystar classification and to measure physical quantities such as
redshifts and stellar masses for the galaxies projected within< 12000 of the lenses, and down to i ∼ 23 − 24 mag. Exceptions are B1608+656 and RXJ1131−1231, where we usedsingle-band HST data within < 4500 (Fassnacht et al. 2011),
and PG 1115+080, where we coadded multiple exposuresof the data used to measure the time delays. With the exception of the first two lenses we analyzed (B1608+656 andRXJ1131−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 κext , whereas for the remaining lenses we also
used constraints from the inferred external shear values ofthe lens models.3 Since HE 0435−1223 (H0LiCOW III), we
have further used combinations of weighted counts to tighten
the κ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 inorder to enforce consistency of technique, as well as to further refine our technique by better accounting for lens-lenscoupling 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 awayfrom the statistical approach and have the potential to significantly tighten the κext PDF (e.g., McCully et al. 2017).
In addition, in some of the H0LiCOW systems, we have independently constrained the external convergence using weaklensing (Tihhonova et al. 2018, hereafter H0LiCOW VIII,
and Tihhonova et al. submitted, hereafter H0LiCOW X).
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, denotedby dimg for the HST and (if available) AO imaging data, ∆t
for the time delays, σ for the velocity dispersion of the lensgalaxy, and dLOS for the properties of the LOS mass distribution determined from our photometric and spectroscopic3

An exception is SDSS 1206+4332, where the use of external
shear has a negligible effect on κext .
MNRAS 000, 1–?? (2019)7

data. We want to obtain the posterior PDF of the model parameters ξ given the data, P (ξ|dimg , ∆t, σ, dLOS , A). The
vector ξ includes the lens model parameters ν, the cosmological parameters π, and nuisance parameters representing
the external convergence (κext ) and anisotropy radius for thelens stellar velocity ellipsoid (rani ). A denotes a discrete setof 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 treatmentof the various deflector mass distributions, etc. In general,
A cannot be fully captured by continuous parameters. From
Bayes’ theorem, we have that
P (ξ|dimg , ∆t, σ, dLOS , A)∝

P (dimg , ∆t, σ, dLOS |ξ, A)P (ξ|A),(5)

where P (dimg , ∆t, σ, dLOS |ξ, A) is the joint likelihood
function and P (ξ|A) is the prior PDF for the parameters
given our assumptions. Since the data sets are independent,
the likelihood can be separated,
P (dimg , ∆t, σ, dLOS |ξ, A)=

P (dimg |ξ, A)×P (∆t|ξ, A)
×P (σ|ξ, A)
×P (dLOS |ξ, A). (6)We can calculate the individual likelihoods separately andcombine 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
A and repeat the inference of ξ. This can be important forchecking the impact of various modeling choices and assumptions, but leaves us with the question of how to combine theresults. Depending on the lens, we either combine the models
with equal weight, or we can use the Bayesian Information
Criterion (BIC) to weight the various models in our final
inference (e.g., H0LiCOW IX, XII, Chen et al. 2019). Thiseffectively combines our various assumptions A using the
BIC so that we obtain P (ξ|dimg , ∆t, σ, dLOS ). We can further marginalize over the non-cosmological parameters (ν,κext , rani ) and obtain the posterior probability distributionof the cosmological parameters π:=

P (π|dimg , ∆t, σ, dLOS )Zdν dκext drani P (ξ|dimg , ∆t, σ, dLOS ).

(7)

In the lens modeling of systems with a single strong lensplane, the parameters associated with cosmology that enter
directly into the lens modeling are the two lensing distancesD∆t and Dd . In the lens modeling of systems with multiple
strong lens planes, we actually vary H0 , keeping other parameters fixed at w = −1, Ωm = 0.3, and ΩΛ = 0.7. This
assumes a fixed curvature of the expansion history of the
Universe, but not the absolute scale (represented by H0 orD∆t ). This is done because there is not a unique D∆t when
accounting for multiple lens planes, but we convert this toan “effective” D∆t that is insensitive to assumptions of thecosmological model. Specifically, given the lens/quasar redshifts and π (i.e., H0 and the other fixed cosmological parameters), we can compute the effective time-delay distanceD∆t (π, zd , zs ) to obtain the posterior probability distribution of D∆t , P (D∆t |dimg , ∆t, σ, dLOS ). In summary, thesingle-lens plane models yield a joint constraints on D∆t

8

K. C. Wong et al.

and Dd , whereas multi-lens plane models yield a constraint
on the effective D∆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 withrespect to the cosmological quantities of interest (i.e., D∆t ,
Dd , H0 ). Throughout the analyses, these values are blindedby 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 isdone 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 for the Sample

We list the D∆t and (when available) Dd constraints foreach individual lens in Table 2, along with correspondingreferences. All distances listed here are used in our cosmological inference. Specifically, for B1608+656, we use theanalytic fit of P (D∆t ) given in Suyu et al. (2010) and of
P (Dd ) given in Jee et al. (2019, in press), 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∆t and (if available) Dd from the analysis ofChen 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∆t and Dd from H0LiCOW IX. For
WFI2033−4723, we use the resulting MCMC chain of D∆tfrom H0LiCOW XII. We use a kernel density estimator tocompute P (D∆t , Dd ) or P (D∆t ) from the chains, allowing
us to account for any correlations between D∆t and Dd in
P (D∆t , Dd ).

4

CHECKING CONSISTENCY AMONGLENSES

We check that all our lenses can be combined without any
loss of consistency by comparing their D∆t posteriors in thefull 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 d1 , ..., dn ,we can either assume the hypothesis Hglobal that they canbe represented using a common global set of cosmologicalparameters, or the hypothesis Hind 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 twolenses,P (di , dj |Hglobal )P (di |Hind )P (dj |Hind )hLi Lj i,=
hLi ihLj i

Fij =

(8)
(9)

where Li and Lj are the likelihoods of data sets di and
dj , respectively. If the six lenses have Bayes factors F > 1
for every possible pairwise combination, it means that thelenses are statistically consistent with each other and we canproceed to combine their constraints.
In Table 3, 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 Hglobal hypothesis. We also test the hypothesis that one
out of six systems is better represented in a different set ofcosmological parameters than the five remaining lenses. The
minimal Bayes factor is obtained for RXJ1131−1231 with
F = 3.8, again in favor of the Hglobal hypothesis, meaningthat all lenses are a consistent realization of the same underlying set of cosmological parameters. We conclude thatnone 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 ouruncertainties 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 4. We distinguish between models where we use constraints from strong lenses alone (Sections 5.1 and 5.2) fromthose 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 isthe same. For the analysis using strong lenses only, we adopt
uniform priors on the cosmological parameters with ranges
indicated in Table 4.
5.1

Flat ΛCDM

Our baseline model is the flat ΛCDM cosmology with auniform prior on H0 in the range [0, 150] km s−1 Mpc−1
and a uniform prior on Ωm in the range [0.05, 0.5]. In Figure 2, we show the marginalized constraints on H0 fromeach of the individual H0LiCOW lenses along with the
combined constraint from all six systems. We find H0 =−1
73.3+1.7
Mpc−1 , a 2.4% precision measurement. We−1.8 km s
show the median and 68% quantiles of the cosmological parameter distributions in Table 5. This estimate is higher than
the Planck Collaboration et al. (2018b) CMB value (H0 =67.4 ± 0.5 km s−1 Mpc−1 ) by 3.1σ and in agreement withMNRAS 000, 1–?? (2019)

H0LiCOW XIII: A 2.4% measurement of H0 from lensed quasars

9

Table 2. D∆t and Dd constraints for H0LiCOW lenses.
Lens nameD∆t (Mpc)

Dd (Mpc)

Blind analysis

References

B1608+656

5156+296−236

1228+177−151

noSuyu et al. (2010); Jee et al. (2019, in press)

RXJ1131−1231

2096+98−83
2707+183−168
5769+589−471
4784+399−248
1470+137−127

804+141−112

yesa

Suyu et al. (2014); Chen et al. (2019)—

yes

Wong et al. (2017); Chen et al. (2019)

1805+555−398

yes

Birrer et al. (2019)b—

yes

Rusu et al. (2019)

697+186−144

yes

Chen et al. (2019)

HE 0435−1223
SDSS 1206+4332
WFI2033−4723
PG 1115+080

Reported values are medians, with errors corresponding to the 16th and 84th percentiles. Values for B1608+656 are calculated from askewed lognormal function fit to the posterior distributions. Values for RXJ1131−1231, HE 0435−1223, SDSS 1206+4332,
WFI2033−4723, and PG 1115+080 are calculated from the MCMC chains of the posterior distributions of the distances. For
RXJ1131−1231, SDSS 1206+4332 and PG 1115+080, we use the joint P (D∆t , Dd ) distributions for cosmographic inferences (seeSection 3.7 for details) to account for correlations between D∆t and Dd .
a The initial HST analysis (Suyu et al. 2013) was performed blindly, but the AO analysis (Chen et al. 2019) was not.
b The values given here are updated values from those presented in Birrer et al. (2019), which had a minor error in the calculation of the 16th and 84th percentiles. The median values are unchanged, while the uncertainties have changed by ∼ 3%.
Table 3. 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

—

3.4

10.6

9.1

11.2

3.3

—

5.3

3.1

4.9

6.9

—

7.8

9.3

3.9

—

7.9

2.7

—

3.7

B1608+656
RXJ1131−1231
HE 0435−1223
SDSS 1206+4332
WFI2033−4723
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

11.1

3.8

11.0

7.5

12.6

4.5

Bayes Factor

The lowest Bayes factor is obtained for PG 1115+080 and SDSS 1206+4332 with an evidence ratio of 2.7, still favoring the hypothesisthat both distributions are two realizations of the same set of cosmological parameters.

the latest SH0ES result (H0 = 74.03 ± 1.42 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 redshiftsystems having a larger inferred value of H0 , but could notconclude 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 H0 with increasing D∆t . Even
with six lenses, these correlations are not significant enoughto 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 thatMNRAS 000, 1–?? (2019)a different set of cosmological parameters is better representing the low redshift lenses. Nevertheless, the persistenceof these trends is something to continue to examine as the
sample of time-delay lenses increases in the future.
5.2

Extensions to Flat ΛCDM, Constraints fromTime-Delay Cosmography Only

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

10

K. C. Wong et al.

Table 4. Description of the cosmological models considered in this work.
Model name

Description

PriorsTime-delay cosmography onlyUΛCDM

Flat ΛCDM

H0 uniform in [0, 150] km s−1 Mpc−1Ωm = 1 − ΩΛΩm uniform in [0.05, 0.5]

UoΛCDM

Open ΛCDM

H0 uniform in [0, 150] km s−1 Mpc−1Ωm uniform in [0.05, 5]Ωk uniform in [−0.5, 0.5]ΩΛ = 1 − Ωm − Ωk > 0

UwCDM

Flat wCDM

H0 uniform in [0,150] km s−1 Mpc−1Ωm uniform in [0.05, 5]
w uniform in [−2.5, 0.5]
ΩDE = 1 − Ωm

Uw0 wa CDM

Flat w0 wa CDM

H0 uniform in [0,150] km s−1 Mpc−1Ωm uniform in [0.05, 5]
w0 uniform in [−2.5, 0.5]
wa uniform in [−2, 2]ΩDE = 1 − Ωm

Time-delay cosmography combined with other probes
ΛCDM

Flat ΛCDMΩm = 1 − ΩΛ
JLA/Pantheon for {H0 , ΩΛ }oΛCDM

Open ΛCDMΩm = 1 − ΩΛ − Ωk > 0
Planck (Section 5.3.1) or JLA/Pantheon (Section 5.4) for {H0 , ΩΛ , Ωm }wCDM

Flat wCDMΩm = 1 − ΩDE
Planck (Section 5.3.2) or JLA/Pantheon (Section 5.4) for {H0 , ΩDE , w}

Neff ΛCDM

Flat ΛCDM
Variable Neff

Planck for {H0 , ΩΛ , Neff }mν ΛCDM

Flat ΛCDM
P
Variablemν

Planck for {H0 , ΩΛ ,Neff mν ΛCDM

Flat ΛCDM P
Variable Neff andmνPlanck for {H0 , ΩΛ , Neff ,w0 wa CDM

Flat w0 wa CDM

Planck (Section 5.3.4) or JLA/Pantheon (Section 5.4) for {H0 , w0 , wa }owCDM

Open wCDMΩm = 1 − ΩDE − Ωk > 0
JLA/Pantheon for {H0 , Ωk , ΩDE , w}

ow0 wa CDM

Open w0 wa CDMΩm = 1 − ΩDE − Ωk > 0
JLA/Pantheon for {H0 , Ωk , ΩDE , w0 , wa }

P

mν }

P

mν }

Planck refers either to the constraints from Planck 2018 Data Release alone, or combined with CMBL or BAO. JLA refers to the joint
light-curve analysis of Betoule et al. (2014). Pantheon refers to the sample of Scolnic et al. (2018).

5.2.1Open ΛCDM

A simple modification to the flat ΛCDM cosmology is an
open ΛCDM cosmology that allows for spatial curvature,Ωk 6= 0. In this model, we have Ωm = 1 − ΩΛ − Ωk . We
adopt uniform prior on Ωk in the range [−0.5, 0.5], ΩΛ inthe range [0,1], and require that Ωm > 0 . We still maintainthe uniform prior on H0 in the range [0, 150] km s−1 Mpc−1 .
The parameter constraints are given in Table 5. Figure 3 shows the marginalized constraint on H0 in an

open ΛCDM cosmology, which we find to be H0 =−1
74.4+2.1
Mpc−1 . This is consistent with our flat
−2.3 km s
Λ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 H0 and Ωk . The black contour is the constraint
from strong lensing alone. We see that Ωk is very poorly
constrained (Ωk = 0.26+0.17−0.25 ). This is not surprising, as the
MNRAS 000, 1–?? (2019)

probability density

H0LiCOW XIII: A 2.4% measurement of H0 from lensed quasars

H0 : 71.0 −+23..93
H0 : 78.2 −+33..44
H0 : 71.7 −+44..85
H0 : 68.9 −+55..41
H0 : 71.6 −+34..89
H0 : 81.1 −+87..01

11

H0 ∈ [0, 150] Ω m ∈ [0.05, 0.5]

All
B1608 (Suyu+2010, Jee+2019)
RXJ1131 (Suyu+2014, Chen+2019)
HE0435 (Wong+2017, Chen+2019)
J1206 (Birrer+2019)
WFI2033 (Rusu+2019)
PG1115 (Chen+2019)

H0 : 73.3 −+11..78

50

60

70
80
H0 [km s −1 Mpc −1 ]90

Figure 2. Marginalized H0 for a flat ΛCDM cosmology with uniform priors. Shown are the H0 posterior PDFs for the individual lenssystems (shaded curves), as well as the combined constraint from all six systems (black line). The median and 16th and 84th percentiles
are shown in the figure legend.

Table 5. Cosmological parameters for various cosmologies from time-delay cosmography only.
Model

H0 (km s−1 Mpc−1 )ΩmΩΛ or ΩDEΩkw or w0

waUΛCDM

73.3+1.7−1.8

0.30+0.13−0.13

0.70+0.13
−0.13

≡0

≡ −1

≡0

UoΛCDM

74.4+2.1−2.3
81.6+4.9−5.3+5.1
81.3−5.4

0.24+0.16−0.13
0.31+0.11
−0.10
+0.11
0.31−0.11

0.51+0.21−0.18
0.69+0.10
−0.11
+0.11
0.69−0.11

0.26+0.17
−0.25

≡ −1

≡0

≡0

−1.90+0.56−0.41−1.86+0.63
−0.45

≡0

UwCDM
Uw0 wa CDM≡0−0.05+1.45−1.37

Reported values are medians, with errors corresponding to the 16th and 84th percentiles.

time-delay distance is only weakly sensitive to Ωm and ΩΛ ,so we would expect a similar insensitivity to Ωk . However,the fact that time-delay cosmography constrains H0 very
tightly indirectly imposes a tight constraint on curvaturewhen combined with other probes.

5.2.2

Flat wCDM

We consider a flat wCDM cosmology in which the dark energy density is not a cosmological constant, but instead is
MNRAS 000, 1–?? (2019)

time-dependent with an equation-of-state parameter w. We
denote the dark energy density parameter as ΩDE = 1 − Ωm .
The w = −1 case corresponds to flat ΛCDM with ΩDE =ΩΛ . We adopt a uniform prior on w in the range [−2.5, 0.5],keeping the same uniform priors on H0 and Ωm as in the flat
ΛCDM model.
We show the parameter constraints in Table 5. In Figure 5, we show the marginalized constraint on H0 in this cos−1
mology, which is H0 = 81.6+4.9
Mpc−1 . The com−5.3 km s
bined constraint on H0 appears to be shifted to a higher

12

K. C. Wong et al.
All
B1608 (Suyu+2010, Jee+2019)
RXJ1131 (Suyu+2014, Chen+2019)
HE0435 (Wong+2017, Chen+2019)
J1206 (Birrer+2019)
WFI2033 (Rusu+2019)
PG1115 (Chen+2019).3
H0 : 71.3 +3−3.5.5
H0 : 78.2 +3−3.5+5
H0 : 71.0 −5..22.7
H0 : 68.1 +5−5.3.6
H0 : 71.0 +4−6.1+8
H0 : 81.2 −7..03

H0 : 74.4 −+22..13

50

60

70
80
−1
H0 [km s Mpc −1 ]90

.9
H0 : 81.6 +4
−5.3

0.8

0.00

70
80
90
−1
−1
H0 [km s Mpc ]

100

wCDM

H0LiCOW
PlanckPlanck + H0LiCOW
Planck + BAO

1.0

0.02

1.2

0.04

w

Ωk

60

Figure 5. Same as Figure 2 for a flat wCDM cosmology.

oΛCDM

1.4
1.6

0.06

1.8

H0LiCOW
PlanckPlanck + H0LiCOW
Planck + CMBL
Planck + BAO

0.08
0.10
45

All
B1608 (Suyu+2010, Jee+2019)
RXJ1131 (Suyu+2014, Chen+2019)
HE0435 (Wong+2017, Chen+2019)
J1206 (Birrer+2019)
WFI2033 (Rusu+2019)
PG1115 (Chen+2019).1
H0 : 70.0 +11−5.2.4
H0 : 79.1 +5−4.2.3
H0 : 73.6 +10−6.9.7
H0 : 71.9 +14−10.1.8H0 : 73.0 +12−7.2+10
H0 : 82.5 −8.3.2

50

Figure 3. Same as Figure 2 for an open ΛCDM cosmology.

0.02

H0 ∈ [0, 150] Ω m ∈ [0.05, 0.5] w ∈ [−2.5, 0.5]

probability density

probability density

H0 ∈ [0, 150] Ω m ∈ [0.05, 0.5] Ω k ∈ [−0.5, 0.5]

50

55

60

−1

65

−1

H0 [km s Mpc ]2.0
70

60 65 70 75 80 85 90 95 100 105

75

Figure 4. H0 -Ωk constraint for an open ΛCDM cosmology.
The black contours show the constraints from H0LiCOW alone,while the grey contours show the constraints from Planck alone.
The colored contours show constraints from Planck combinedwith other probes, as shown in the figure legend. Although the
H0LiCOW and Planck constraints are not consistent with eachother, we show the combination here for completeness. The contour levels represent the 1σ and 2σ constraints.

value than most of the individual lenses, but this is dueto the degeneracy between H0 and w and the resulting
asymmetry in the PDF when projected from the higherdimensional cosmological parameter space. Indeed, there isa region in the two-dimensional H0 −w plane in which all
the individual distributions converge. That region has ahigher probability density when performing the joint inference, which in turn drives the marginalized H0 value higher.
In Figure 6, we show the joint distribution of H0 andH0 [km s −1 Mpc −1 ]Figure 6. H0 -w constraint for a flat wCDM cosmology. The black
contours show the constraints from H0LiCOW alone, while thegrey contours show the constraints from Planck alone. The colored contours show constraints from Planck combined with otherprobes, as shown in the figure legend. The contour levels represent
the 1σ and 2σ constraints.

w. Lensing alone does not constrain w particularly well
(w = −1.90+0.56−0.41 ), although there is a degeneracy between
w and H0 , suggesting that combining our constraint withother probes may produce useful constraints.

5.2.3

Flat w0 wa CDM

The flat wCDM 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
w0 wa CDM cosmology in which the dark energy equationMNRAS 000, 1–?? (2019)

w0

H0LiCOW XIII: A 2.4% measurement of H0 from lensed quasars

0.5
0.0
0.5
1.0
1.5
2.0

w0 wa CDM
H0LiCOW
Planck + BAO
Planck + BAO + H0LiCOW

wa

1
2
3

60

70

−1

80

90

−1

H0 [km s Mpc ]

2

1

w0

0

Figure 7. Constraints on H0 , w0 , and wa for a w0 wa CDM cosmology. The colored contours show constraints from Planck combined with other probes, as shown in the figure legend (no chainsfor Planck alone are available for this cosmology). The contour
levels represent the 1σ and 2σ constraints. H0LiCOW alone places
effectively no constraints on wa with the resulting posterior (openblack contours) spanning the prior range of [−2, 2].

of-state parameter w is time-varying and parameterized asw(z) = w0 + wa z/(1 + z) (Linder 2003). We adopt a uniform
prior on w0 in the range [−2.5, 0.5] and on wa in the range[−2, 2], keeping the same uniform priors on H0 and Ωm asin the flat ΛCDM model.
We show the parameter constraints in Table 5. In
Figure 7, we show the joint constraints on H0 , w0 , andwa from the six lenses in open black contours. Unsurprisingly, the lenses provide little constraint on w0 and effectively no constraint on wa , with the posterior PDF on
wa spanning the entire prior range of [−2, 2]. The resulting H0 with a time-varying w(z) is similarly high (H0 =−1
81.3+5.1
Mpc−1 ), as in the case of a flat wCDM cos−5.4 km smology, due to parameter degeneracies between H0 and w0 .

5.3

results, using the chains provided by the Planck team4 . 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 inTable 6.

5.3.1

0

Extensions to Flat Λ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 ΛCDM flat wCDM, and flat w0 wa CDM cosmologies,
as well as a number of other extensions. Although timedelay cosmography is primarily sensitive to H0 and only
weakly dependent on other parameters, the degeneracies are
such that strong lensing is highly complementary to othercosmological probes. Here, we combine strong lensing with
CMB observations from Planck (Planck Collaboration et al.
2018b), sometimes including CMB weak lensing (CMBL)MNRAS 000, 1–?? (2019)13

Open ΛCDM

We test the open Λ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 givea strong degeneracy between H0 and Ωk , but the additionof time-delay cosmography information provides importantcomplementary information that constrains the Universe to
be nearly flat (Ω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 wCDM

We test the flat wCDM 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 areslightly different, allowing the combination of the two probes
to constrain w = 1.36+0.09−0.09 . There is mild tension at the ∼ 2σ
level between the Planck+H0LiCOW constraints discussed
above, and Planck+BAO constraints.

5.3.3

Flat ΛCDM with variable neutrino species and/or
masses

The standard model has an effective number of primordial
neutrino species Neff = 3.046, but additional relativistic particles in the early Universe prior to recombination could, inprinciple, add to this quantity. The sum of neutrino masses
has not been precisely measured, but limits have been placed
by aPvariety of experiments. In our analysis thus far, we havesetmν = 0.06, the minimum mass allowed by neutrino oscillation experiments
P (Patrignani 2016). We test cosmologies
in which Neff ormν 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 H0 and Neff (when Neffisallowedto vary).
Figure 9 shows a contour plot of H0 andPPmν (whenmν 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 combined4

http://pla.esac.esa.int/pla/#cosmology

14

K. C. Wong et al.

Table 6. Cosmological parameters for various cosmologies from time-delay cosmography combined with other probes. The oΛCDM andwCDM constraints can be compared with those from time-delay cosmography only (Table 5).
H0 (km s−1 Mpc−1 )ΩmΩΛ or ΩDEΩkw or w0

wa

Neff

Planck

67.0+1.3−1.4

0.32+0.01−0.01

0.68+0.01−0.01−0.002+0.003−0.003≡ −1≡0≡ 3.046≡ 0.06oΛCDM

Planck + CMBL

0.003+0.003−0.0040.002+0.002−0.002≡0≡ 3.046≡ 0.06

Planck + BAO

0.70+0.01−0.01
0.69+0.01−0.01
0.77+0.01−0.02
0.71+0.01
−0.01

≡ −1

oΛCDM

0.30+0.02−0.01
0.30+0.01−0.01
0.23+0.02−0.01
0.29+0.01
−0.01

≡ −1

≡0

≡ 3.046≡ 0.06≡0≡0≡ 3.046≡ 0.06≡0−1.36+0.09−0.09−1.10+0.06
−0.06

≡0

≡ 3.046≡ 0.06

≡0

≡ −1

≡0

3.20+0.15−0.16≡ 0.06≡ 0.06Model

H0LiCOW +

oΛCDM

wCDM

PlanckwCDM

Planck + BAO

69.2+1.7−1.8
68.4+0.7−0.7
78.7+2.7−2.7
70.3+1.6−1.5

Neff ΛCDM

Planck

68.8+1.2−1.3

0.31+0.01−0.01

0.69+0.01−0.01

Neff ΛCDM

Planck + BAO

mν ΛCDM

Planckmν ΛCDM

Planck + CMBLmν ΛCDM

Planck + BAO

69.0+1.1−1.1
68.0+0.7−0.7
68.0+0.6−0.7
68.0+0.5−0.5
69.1+1.3−1.3
69.0+1.0−1.1
67.1+3.0−2.5

0.30+0.01−0.01
0.31+0.01−0.01
0.31+0.01−0.01
0.31+0.01−0.01
0.30+0.01−0.01
0.30+0.01−0.01
0.32+0.03−0.03

0.70+0.01−0.01
0.69+0.01−0.01
0.69+0.01−0.01
0.69+0.01−0.01
0.70+0.01−0.01
0.70+0.01−0.01
0.68+0.03−0.03

Neff mν ΛCDM Planck
Neff mν ΛCDM Planck + BAO
w0 wa CDM

Planck + BAOPmν (eV)

≡0

≡ −1

≡0

3.23+3.40
−3.23

≡0

≡ −1

≡0

≡ 3.046

0.03+0.08
−0.03

≡0

≡ −1

≡0

≡ 3.046

0.01+0.05
−0.02

≡0

≡ −1

≡0

≡ 3.046

0.03+0.04
−0.02

≡0

≡ −1

≡0

3.21+0.20−0.17

0.03+0.05
−0.03

≡0

≡ −1

≡0

3.21+0.16−0.17

0.03+0.05
−0.02

≡0

−0.73+0.29−0.30−1.02+0.80−0.88≡ 3.046≡ 0.06

Reported values are medians, with errors corresponding to the 16th and 84th percentiles.

3.6

H0LiCOW
PlanckPlanck + H0LiCOW
Planck + BAO

N eff ΛCDM

0.40
0.35

3.4

H0LiCOW
PlanckPlanck + H0LiCOW
Planck + CMBL
Planck + BAO

m ν ΛCDM

0.30

Σm ν [eV]

3.2

N eff

0.25

3.0

0.20
0.15

2.8

0.10

2.6
2.4

0.45

0.05
64

66

68

70

−1

72

−1

H0 [km s Mpc ]74

76

Figure 8. H0 -Neff constraint for a Neff ΛCDM cosmology. The
black contours show the constraints from H0LiCOW alone, while
the grey contours show the constraints from Planck alone. The
colored contours show constraints from Planck combined withother probes, as shown in the figure legend. The contour levels
represent the 1σ and 2σ constraints.

P
constraints of Neff = 3.21+0.16mν = 0.03+0.05−0.17 and−0.03 , consistent with the standard model values. Although we see thatthe H0LiCOW and Planck constraints in these cosmologies
are somewhat in tension, we provide the joint constraints
here for completeness.

0.00

64

66

68

70

−1

72

−1

H0 [km s Mpc ]

74

76

P
Figure 9. H0 - mν constraint for a mν ΛCDM cosmology. The
black contours show the constraints from H0LiCOW alone, while
the grey contours show the constraints from Planck alone. The
colored contours show constraints from Planck combined withother probes, as shown in the figure legend. The contour levels
represent the 1σ and 2σ constraints.

5.3.4

Flat w0 wa CDM

We consider the same flat w0 wa CDM cosmology as in Section 5.2.3, now including Planck and BAO as external data
sets. The joint constraints on H0 , w0 , and wa are shown in
MNRAS 000, 1–?? (2019)

H0LiCOW XIII: A 2.4% measurement of H0 from lensed quasars

15

Table 7. Cosmological parameters for various cosmologies when anchoring type Ia supernovae distances with distances from time-delay
cosmography.
JLA sample (Betoule et al. 2014) and H0LiCOWModel

H0 (km s−1 Mpc−1 )ΩmΩΛ or ΩDEΩkw or w0wa

ΛCDM

73.6+1.6−1.8

0.30+0.03−0.03

0.70+0.03
−0.03

≡0

≡ −1

≡0

wCDM

73.9+2.6−2.7
74.2+2.5−2.6
75.2+1.8−2.1
77.5+2.5−2.9
77.7+2.4−2.8

0.31+0.08−0.10
0.34+0.08−0.10
0.18+0.08−0.06
0.22+0.07−0.07
0.23+0.07−0.06

0.69+0.10−0.08
0.66+0.10−0.08
0.53+0.12−0.09
0.41+0.14−0.08
0.40+0.13
−0.07

≡0

−1.05+0.23−0.28−1.04+0.22
−0.28

≡0
−0.43+1.10
−1.07

≡ −1

≡0

−1.46+0.42−0.53−1.49+0.44
−0.52

≡0
−0.16+1.24−1.24

w0 wa CDM
oΛCDM
owCDM
ow0 wa CDM≡00.29+0.14−0.19
0.37+0.09−0.17
0.38+0.09−0.16

Pantheon sample (Scolnic et al. 2018) and H0LiCOWModel

H0 (km s−1 Mpc−1 )ΩmΩΛ or ΩDEΩkw or w0wa

ΛCDM

73.6+1.6−1.8

0.30+0.02−0.02

0.70+0.02
−0.02

≡0

≡ −1

≡0

wCDM

74.9+2.2−2.4
75.0+2.2−2.3
73.8+1.9−2.1
77.4+2.5−3.0
77.4+2.5−3.0

0.35+0.05−0.06
0.37+0.05−0.07
0.28+0.07−0.07
0.28+0.07−0.06
0.28+0.08−0.07

0.65+0.06−0.05
0.63+0.07−0.05
0.67+0.11−0.11
0.45+0.17−0.10
0.45+0.17
−0.09

≡0

−1.17+0.19−0.22−1.15+0.19
−0.23

≡0
−0.55+1.22−1.01

0.04+0.18
−0.17

≡ −1

≡0

0.26+0.16−0.21−1.52+0.39
−0.52

≡0

0.26+0.15−0.21−1.50+0.40−0.47−0.31+1.31−1.17

w0 wa CDM
oΛCDM
owCDM
ow0 wa CDM≡0

Reported values are medians, with errors corresponding to the 16th and 84th percentiles.

N eff m ν ΛCDM

3.5

N eff

H0LiCOW
PlanckPlanck + H0LiCOW
Planck + BAO

3.0
2.5

Σm ν [eV]

0.4

0.2

0.0

65

70

75

H0 [km s −1 Mpc −1 ]

2.5

3.0

N eff

3.5

P
Figure 10. Constraints on H0 , Neff , andmν for a mν ΛCDM
cosmology. The black contours show the constraints from
H0LiCOW alone, while the grey contours show the constraints
from Planck alone. The colored contours show constraints from
Planck combined with other probes, as shown in the figure legend.
The contour levels represent the 1σ and 2σ constraints.

MNRAS 000, 1–?? (2019)

Figure 7. The combination of H0LiCOW, Planck, and BAO+0.80
constraints finds w0 = −0.73+0.29−0.30 and wa = −1.02−0.88 .
Time-delay cosmography does not add much information
due to the large degeneracies.
5.4

Calibrating Type Ia Supernovae Distanceswith Time-Delay Cosmography

The distance ladder method uses local distance indicators(particularly Cepheid variables) to calibrate the absolute
distances to type Ia SNe. In principle, any absolute distancemeasurement can be used to anchor the distance scale toSNe. 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 (Dd ) to lens galaxies (e.g., Wojtak & Agnello 2019, Jee et al. 2019, in press) or D∆t (e.g.,
Taubenberger et al. 2019). This can be used as a complementary probe of H0 or the sound horizon, rs (Arendse et al.
2019), that is nearly insensitive to the assumed cosmology.
We follow the methodology of Taubenberger et al.
(2019), using the combined D∆t and Dd measurements fromthe 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

16

K. C. Wong et al.

2011), and the JLA and Pantheon SNe samples implemented
within MontePython, to combine with our lensing distancemeasurements and sample cosmological parameters.
We consider six cosmological models that are listed
in part of Table 4: ΛCDM, oΛCDM, wCDM, w0 wa CDM,
owCDM and ow0 wa CDM. For the first four cosmologicalmodels, we adopt the same uniform prior ranges for the cosmological parameters as in the top part of Table 4 (i.e., samepriors as UΛCDM, UoΛCDM, UwCDM and Uw0 wa CDM,
respectively). For owCDM and ow0 wa CDM, we adopt thesame priors for w and {w0 , wa } as those in UwCDM and
Uw0 wa CDM, respectively, with the remaining parameters(H0 , Ωm , Ωk and ΩDE ) having the same priors as those inoΛCDM.
Table 7 shows the results for the six cosmologies tested,
for the JLA sample (top) and the Pantheon sample (bottom). The median values of H0 from the JLA sample arewithin 1.5 km s−1 Mpc−1 of those of the Pantheon sample.
In comparison to the H0 constraints in Section 5.2, particularly for the wCDM and w0 wa CDM models where H0
is highly degenerate with w, w0 and wa , the H0 from thelenses and SNe in Table 7 are less sensitive to cosmological models, as shown in Figure 11. For the six cosmologicalmodels probed in Table 7, the median H0 values range from∼73−78 km s−1 Mpc−1 , irrespective of the SN sample. The
tension with Planck in flat ΛCDM is still > 3σ, similar to
our result from time-delay cosmography alone.
We also repeat the JLA analysis using only themarginalized P (D∆t ) from the six lenses as constraints, i.e.,omitting the information from Dd . This allows us to assess
the information content on H0 from the lensing distances
and the added value of measuring Dd in addition to D∆t .
Using only P (D∆t ), we find that the H0 values are very
similar to values in Table 7 (within 0.5 km s−1 Mpc−1 ), andthe uncertainties are also only slightly larger (by at most
0.5 km s−1 Mpc−1 ). Therefore, most of the cosmological information, particularly H0 , is encapsulated in our D∆t measurements that are substantially more precise than Dd measurements. Nonetheless, future spatially resolved kinematicsof the lens galaxy could help tighten the constraints on Dd ,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 H0

As the tension between different probes of H0 began toemerge in recent years, a natural direction to look towardin 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 themost precise constraints on H0 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 ΛCDM cosmology are in 3.1σ tension with
Planck. In combination with the latest SH0ES result (Riess

H0LiCOWH0LiCOW + JLA
H0LiCOW + Pantheon

ΛCDM
wCDM
w0 wa CDM
oΛCDM
owCDM
ow0 wa CDM

70

75
80
−1
H0 [km s Mpc−1 ]85Figure 11. Constraints on H0 in various cosmological models(described in Table 4) from the six H0LiCOW lenses (blue triangles, Table 5) and from the combination of lenses and type Ia
SNe (Table 7) using the JLA (red circles) and Pantheon (red diamonds) samples. The points represent the median, with the error
bar showing the 16th and 84th percentiles. In the owCDM and
ow0 wa CDM cosmologies, the cosmological parameter samplingsof the lenses only (without SNe) do not converge due to the multiple ill-constrained parameters, and the H0 values are thus notreported. By anchoring the type Ia SNe distance scale with thelensing distances, the inferred H0 is less sensitive to cosmologicalmodel assumptions, in comparison to the constraints from lenses
alone.

et al. 2019), we find a 5.3σ tension between late-Universe
determinations of H0 and Planck (Figure 12). Other independent methods anchored in the early Universe, such as
the analysis of Abbott et al. (2018b) using a combination ofclustering and weak lensing, BAO, and big bang nucleosynthesis (BBN), give similar results to Planck.
Given the tension shown here, it is becoming difficult
to reconcile the H0 discrepancy by appealing to systematiceffects. While systematics, especially “unknown unknowns”,
still cannot be entirely ruled out and should continue to be
explored, recent work has only heightened the tension. Therealso appears to be a growing dichotomy when the various H0probes are split into those anchored by the early-Universe(i.e., CMB), which favor a lower H0 , and those based onlate-Universe probes, which favor a higher H0 .
As this tension between early-Universe and lateUniverse probes continues to grow, we must examine potential alternatives to the standard flat ΛCDM model. ThisMNRAS 000, 1–?? (2019)

H0LiCOW XIII: A 2.4% measurement of H0 from lensed quasars

17

Figure 12. Comparison of H0 constraints for early-Universe and late-Universe probes in a flat ΛCDM cosmology. The early-Universe
probes shown here are from Planck (orange; Planck Collaboration et al. 2018b) and a combination of clustering and weak lensing data,
BAO, and big bang nucleosynthesis (grey; Abbott et al. 2018b). The late-Universe probes shown are the latest results from SH0ES (blue;
Riess et al. 2019) and H0LiCOW (red; this work). When combining the late-Universe probes (purple), we find a 5.3σ tension with Planck.

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 toflat ΛCDM in Section 5, including spatial curvature, timevarying dark energy (e.g., Di Valentino et al. 2018), andmodified 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; Aylor et al. 2019; Lin et al.
2019; Poulin et al. 2019), primordial non-Gaussianity (e.g.,
Adhikari & Huterer 2019), and decaying dark matter (e.g.,
Pandey et al. 2019; Vattis 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 H0 to date fromstrong lensing time delays. Five of the six lenses are analyzed
blindly with respect to the cosmological parameters of interest. Our main results are:−1
Mpc−1 for a flat ΛCDM
• We find H0 = 73.3+1.7−1.8 km s

MNRAS 000, 1–?? (2019)

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σ tension with Planck CMB
measurements (Planck Collaboration et al. 2018b).
• Our constraint on H0 in flat ΛCDM is completely independent of and complementary to the latest results fromthe SH0ES collaboration, so these two measurements can be
combined into a late-Universe constraint on H0 . Together,
these are in tension with the best early-Universe (i.e., CMB)
determination of H0 from Planck at a significance of 5.3σ.
• We check that the lenses in our sample are statisticallyconsistent with one another by computing Bayes factors between their H0 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 systematiceffects in our analysis.
• We compute parameter constraints for cosmologies beyond flat ΛCDM. In an open ΛCDM cosmology, we find
+2.1−1
Mpc−1 , whichΩk = 0.26+0.17−0.25 and H0 = 74.4−2.3 km s
is still in tension with Planck, suggesting that allowing for
spatial curvature cannot resolve the discrepancy. In a flat
−1
wCDM cosmology, we find H0 = 81.6+4.9
Mpc−1−5.3 km s+0.56
and w = −1.90−0.41 . In a flat w0 wa CDM cosmology, we18

K. C. Wong et al.

−1find H0 = 81.3+5.1
Mpc−1 , but are unable to place−5.4 km s
meaningful constraints on w0 and wa .
• We combine our constraints with Planck, including
CMB weak lensing and BAO constraints. Although timedelay cosmography is primarily sensitive to H0 , with only
a weak dependence on other cosmological parameters, theconstraints are highly complementary to other probes suchas Planck, CMB weak lensing, and BAO. We test the open
ΛCDM and wCDM cosmologies, as well as cosmologies withvariable effective neutrino species and/or sum of neutrino
masses, and a wCDM cosmology with a time-varying w. The
full parameter constraints for these models when combining
H0LiCOW and Planck are given in Table 6.
• We use the distance measurements from time-delay cosmography to calibrate the distance scale of type Ia SNe fromthe JLA and Pantheon samples. This provides a probe of H0
that is less dependent on the assumed cosmological model,
in comparison to the constraints from lenses alone. We findmedian H0 values ranging from ∼73−78 km s−1 Mpc−1 fora range of cosmologies. The tension with Planck for a flat
ΛCDM cosmology is still > 3σ, similar to the result fromtime-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 H0 has only continued toheighten. If unresolved, this tension may force the rejectionof the flat ΛCDM model in favor of new physics, which would
dramatically alter our understanding of the Universe.
While considering the possibility of new physics, we arealso 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 numberof 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,b; Anguita et al. 2018; Lemon et al. 2018,2019; Treu et al. 2018). A sample of ∼ 40 lenses is needed
to constrain H0 to the ∼ 1% level (Jee et al. 2016; Shajib
et al. 2019), which will be attainable in the near future.

ACKNOWLEDGEMENTS
We thank Chiara Spiniello, Malte Tewes, and Akın Yıldırım
for their contributions to the H0LiCOW project. We thankLodovico 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 PremierInternational 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 Instituteof 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 fromVILLUM FONDEN (project number 16599). This projectis 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, HSTGO-12889, and HST-14254. Support for programs HST-GO10158 HST-GO-12889 HST-14254 was provided to membersof 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., underNASA contract NAS 5-26555. This research made use of
Astropy,5 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).

REFERENCES
Abbott B. P., et al., 2017, Nature, 551, 85
Abbott T. M. C., et al., 2018a, Phys. Rev. D, 98, 043526
Abbott T. M. C., et al., 2018b, MNRAS, 480, 3879
Addison G. E., Watts D. J., Bennett C. L., Halpern M., Hinshaw
G., Weiland J. L., 2018, ApJ, 853, 119
Adhikari S., Huterer D., 2019, arXiv e-prints, p. arXiv:1905.02278
Agnello A., et al., 2015, MNRAS, 454, 1260
Agnello A., Sonnenfeld A., Suyu S. H., Treu T., Fassnacht C. D.,
Mason C., Bradač M., Auger M. W., 2016, MNRAS, 458, 3830

5

http://www.astropy.org
MNRAS 000, 1–?? (2019)

H0LiCOW XIII: A 2.4% measurement of H0 from lensed quasars
Agnello A., et al., 2018a, MNRAS, 475, 2086
Agnello A., et al., 2018b, MNRAS, 479, 4345
Agrawal P., Cyr-Racine F.-Y., Pinner D., Randall L., 2019, arXiv
e-prints, p. arXiv:1904.01016
Alam S., et al., 2017, MNRAS, 470, 2617
Anguita T., et al., 2018, MNRAS, 480, 5017
Arendse N., Agnello A., Wojtak R., 2019, arXiv e-prints, p.
arXiv:1905.12000
Astropy Collaboration et al., 2013, A&A, 558, A33
Astropy Collaboration et al., 2018, AJ, 156, 123
Aubourg É., et al., 2015, Phys. Rev. D, 92, 123516
Audren B., Lesgourgues J., Benabed K., Prunet S., 2013, J. Cosmology Astropart. Phys., 2013, 001
Aylor K., Joy M., Knox L., Millea M., Raghunathan S., Kimmy
Wu W. L., 2019, ApJ, 874, 4
Bernal J. L., Verde L., Riess A. G., 2016, J. Cosmology Astropart.
Phys., 2016, 019
Betoule M., et al., 2014, A&A, 568, A22
Beutler F., et al., 2011, MNRAS, 416, 3017
Birrer S., Amara A., 2018, Physics of the Dark Universe, 22, 189
Birrer S., Amara A., Refregier A., 2015, ApJ, 813, 102
Birrer S., Amara A., Refregier A., 2016, J. Cosmology Astropart.
Phys., 8, 020
Birrer S., et al., 2019, MNRAS, 484, 4726
Blandford R., Narayan R., 1986, ApJ, 310, 568
Bonamente M., Joy M. K., LaRoque S. J., Carlstrom J. E., Reese
E. D., Dawson K. S., 2006, ApJ, 647, 25
Bonvin V., Tewes M., Courbin F., Kuntzer T., Sluse D., MeylanG., 2016, A&A, 585, A88
Bonvin V., et al., 2017, MNRAS, 465, 4914
Bonvin V., et al., 2018, A&A, 616, A183
Bonvin V., et al., 2019, arXiv e-prints, p. arXiv:1905.08260
Braatz J., et al., 2018, in Tarchi A., Reid M. J., Castangia P., eds, IAU Symposium Vol. 336, Astrophysical
Masers: Unlocking the Mysteries of the Universe. pp 86–91,
doi:10.1017/S1743921317010249
Brinckmann T., Lesgourgues J., 2018, arXiv e-prints(arXiv:1804.07261), p. arXiv:1804.07261
Chávez R., Terlevich E., Terlevich R., Plionis M., Bresolin F.,
Basilakos S., Melnick J., 2012, MNRAS, 425, L56
Chen G. C.-F., et al., 2016, MNRAS, 462, 3457
Chen G. C.-F., et al., 2018, MNRAS, 481, 1115
Chen G. C. F., et al., 2019, arXiv e-prints, p. arXiv:1907.02533
Coe D., Moustakas L. A., 2009, ApJ, 706, 45
Collett T. E., Auger M. W., 2014, MNRAS, 443, 969
Collett T. E., et al., 2013, MNRAS, 432, 679
Courbin F., Eigenbrod A., Vuissoz C., Meylan G., Magain P.,
2005, in Mellier Y., Meylan G., eds, IAU Symposium Vol.
225, Gravitational Lensing Impact on Cosmology. pp 297–303,
doi:10.1017/S1743921305002097
Courbin F., et al., 2011, A&A, 536, A53
Courbin F., et al., 2018, A&A, 609, A71
Cuesta A. J., Verde L., Riess A., Jimenez R., 2015, MNRAS, 448,
3463
Di Valentino E., Linder E. V., Melchiorri A. r., 2018, Phys.
Rev. D, 97, 043528
Ding X., et al., 2018, arXiv e-prints, p. arXiv:1801.01506
Domı́nguez A., Prada F., 2013, ApJ, 771, L34
Domı́nguez A., et al., 2019, arXiv e-prints, p. arXiv:1903.12097
Dressler A., 1980, ApJ, 236, 351
Dye S., Warren S. J., 2005, ApJ, 623, 31
Eigenbrod A., Courbin F., Vuissoz C., Meylan G., Saha P., Dye
S., 2005, A&A, 436, 25
Eigenbrod A., Courbin F., Meylan G., Vuissoz C., Magain P.,
2006, A&A, 451, 759
Eulaers E., et al., 2013, A&A, 553, A121
Fadely R., Keeton C. R., Nakajima R., Bernstein G. M., 2010,
ApJ, 711, 246
MNRAS 000, 1–?? (2019)19

Falco E. E., Gorenstein M. V., Shapiro I. I., 1985, ApJ, 289, L1
Fassnacht C. D., Womble D. S., Neugebauer G., Browne I. W. A.,
Readhead A. C. S., Matthews K., Pearson T. J., 1996, ApJ,
460, L103
Fassnacht C. D., Pearson T. J., Readhead A. C. S., Browne
I. W. A., Koopmans L. V. E., Myers S. T., Wilkinson P. N.,
1999, ApJ, 527, 498
Fassnacht C. D., Xanthopoulos E., Koopmans L. V. E., Rusin D.,
2002, ApJ, 581, 823
Fassnacht C. D., Koopmans L. V. E., Wong K. C., 2011, MNRAS,
410, 2167
Feeney S. M., Peiris H. V., Williamson A. R., Nissanke S. M.,
Mortlock D. J., Alsing J., Scolnic D., 2019, Physical Review
Letters, 122, 061105
Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013,
PASP, 125, 306
Freedman W. L., 2017, Nature Astronomy, 1, 0121
Freedman W. L., Madore B. F., Scowcroft V., Burns C., Monson
A., Persson S. E., Seibert M., Rigby J., 2012, ApJ, 758, 24
Gall E. E. E., Kotak R., Leibundgut B., Taubenberger S., Hillebrandt W., Kromer M., 2016, A&A, 592, A129
Gelmini G. B., Kusenko A., Takhistov V., 2019, arXiv e-prints,
p. arXiv:1906.10136
Goldstein D. A., Nugent P. E., 2017, ApJ, 834, L5
Goldstein D. A., Nugent P. E., Kasen D. N., Collett T. E., 2018,
ApJ, 855, 22
González-Morán A. L., et al., 2019, MNRAS, p. 1512
Goobar A., et al., 2017, Science, 356, 291
Gorenstein M. V., Shapiro I. I., Falco E. E., 1988, ApJ, 327, 693
Greene Z. S., et al., 2013, ApJ, 768, 39
Henning J. W., et al., 2018, ApJ, 852, 97
Herrnstein J. R., et al., 1999, Nature, 400, 539
Heymans C., et al., 2012, MNRAS, 427, 146
Hilbert S., Hartlap J., White S. D. M., Schneider P., 2009, A&A,
499, 31
Huber S., et al., 2019, arXiv e-prints, p. arXiv:1903.00510
Humphreys E. M. L., Reid M. J., Moran J. M., Greenhill L. J.,
Argon A. L., 2013, ApJ, 775, 13
Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
Jee I., Komatsu E., Suyu S. H., 2015, Journal of Cosmology andAstro-Particle Physics, 2015, 033
Jee I., Komatsu E., Suyu S. H., Huterer D., 2016, Journal of
Cosmology and Astro-Particle Physics, 2016, 031
Jones D. O., et al., 2018, ApJ, 867, 108
Keeton C. R., 2003, ApJ, 584, 664
Keeton C. R., Kochanek C. S., 1997, ApJ, 487, 42
Kelly P. L., et al., 2015, Science, 347, 1123
Kochanek C. S., 2002, ApJ, 578, 25
Kochanek C. S., 2003, ApJ, 583, 49
Kochanek C. S., Morgan N. D., Falco E. E., McLeod B. A., Winn
J. N., Dembicky J., Ketzeback B., 2006, ApJ, 640, 47
Koopmans L. V. E., 2004, preprint, (arXiv:astro-ph/0412596)
Koopmans L. V. E., Treu T., Fassnacht C. D., Blandford R. D.,
Surpi G., 2003, ApJ, 599, 70
Kovner I., 1987, ApJ, 316, 52
Kozmanyan A., Bourdin H., Mazzotta P., Rasia E., Sereno M.,
2019, A&A, 621, A34
Kreisch C. D., Cyr-Racine F.-Y., Doré O., 2019, arXiv e-prints,
p. arXiv:1902.00534
Lemon C. A., Auger M. W., McMahon R. G., Ostrovski F., 2018,
MNRAS, 479, 5060
Lemon C. A., Auger M. W., McMahon R. G., 2019, MNRAS, 483,
4242
Lesgourgues J., 2011, arXiv e-prints, p. arXiv:1104.2932
Liao K., et al., 2015, ApJ, 800, 11
Lin M.-X., Benevento G., Hu W., Raveri M., 2019, arXiv e-prints,
p. arXiv:1905.12618
Linder E. V., 2003, Physical Review Letters, 90, 091301

20

K. C. Wong et al.

Linder E. V., 2011, Phys. Rev. D, 84, 123529
Macaulay E., et al., 2019, MNRAS, 486, 2184
Marshall P., Rajguru N., Slosar A., 2006, Phys. Rev. D, 73, 067302
McCully C., Keeton C. R., Wong K. C., Zabludoff A. I., 2014,
MNRAS, 443, 3631
McCully C., Keeton C. R., Wong K. C., Zabludoff A. I., 2017,
ApJ, 836, 141
Melnick J., Terlevich R., Terlevich E., 2000, MNRAS, 311, 629
Morgan N. D., Caldwell J. A. R., Schechter P. L., Dressler A.,
Egami E., Rix H.-W., 2004, AJ, 127, 2617
Morgan N. D., Kochanek C. S., Pevunova O., Schechter P. L.,
2005, AJ, 129, 2531
Myers S. T., et al., 1995, ApJ, 447, L5
Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563
Oguri M., 2007, ApJ, 660, 1
Oguri M., Marshall P. J., 2010, MNRAS, 405, 2579
Oguri M., et al., 2005, ApJ, 622, 106
Pandey K. L., Karwal T., Das S., 2019, arXiv e-prints, p.
arXiv:1902.10636
Paraficz D., Hjorth J., 2009, A&A, 507, L49
Patrignani C., 2016, Chinese Physics C, 40, 100001
Petters A. O., Levine H., Wambsganss J., 2001, Singularity theory
and gravitational lensing
Planck Collaboration et al., 2018a, preprint, (arXiv:1807.06205)
Planck Collaboration et al., 2018b, preprint, (arXiv:1807.06209)
Poulin V., Smith T. L., Karwal T., Kamionkowski M., 2019, Phys.
Rev. Lett., 122, 221301
Rathna Kumar S., Stalin C. S., Prabhu T. P., 2015, A&A, 580,
A38
Reese E. D., Carlstrom J. E., Joy M., Mohr J. J., Grego L.,
Holzapfel W. L., 2002, ApJ, 581, 53
Refregier A., 2003, MNRAS, 338, 35
Refsdal S., 1964, MNRAS, 128, 307
Riess A. G., et al., 2016, ApJ, 826, 56
Riess A. G., et al., 2018, ApJ, 861, 126
Riess A. G., Casertano S., Yuan W., Macri L. M., Scolnic D.,
2019, ApJ, 876, 85
Rigault M., et al., 2015, ApJ, 802, 20
Rigault M., et al., 2018, arXiv e-prints, p. arXiv:1806.03849
Ross A. J., Samushia L., Howlett C., Percival W. J., Burden A.,
Manera M., 2015, MNRAS, 449, 835
Rusu C. E., et al., 2017, MNRAS, 467, 4220
Rusu C. E., et al., 2019, arXiv e-prints, p. arXiv:1905.09338
Saha P., 2000, AJ, 120, 1654
Saha P., Coles J., Macciò A. V., Williams L. L. R., 2006, ApJ,
650, L17
Salamon M. H., Stecker F. W., de Jager O. C., 1994, ApJ, 423,
L1
Sandage A., Tammann G. A., Saha A., Reindl B., Macchetto
F. D., Panagia N., 2006, ApJ, 653, 843
Schechter P. L., et al., 1997, ApJ, 475, L85
Schmidt B. P., Kirshner R. P., Eastman R. G., Phillips M. M.,
Suntzeff N. B., Hamuy M., Maza J., Aviles R., 1994, ApJ,
432, 42
Schneider P., Sluse D., 2013, A&A, 559, A37
Schneider P., Ehlers J., Falco E. E., 1992, Gravitational Lenses.
Springer, doi:10.1007/978-3-662-03758-4
Scolnic D. M., et al., 2018, ApJ, 859, 101
Seljak U., 1994, ApJ, 436, 509
Sereno M., Paraficz D., 2014, MNRAS, 437, 600
Shajib A. J., Treu T., Agnello A., 2018, MNRAS, 473, 210
Shajib A. J., et al., 2019, MNRAS, 483, 5649
Silk J., White S. D. M., 1978, ApJ, 226, L103
Sluse D., et al., 2003, A&A, 406, L43
Sluse D., Claeskens J. F., Hutsemékers D., Surdej J., 2007, A&A,
468, 885
Sluse D., Hutsemékers D., Courbin F., Meylan G., Wambsganss
J., 2012, A&A, 544, A62

Sluse D., et al., 2017, MNRAS, 470, 4838
Sluse D., et al., 2019, arXiv e-prints, p. arXiv:1905.08800
Soares-Santos M., et al., 2019, ApJ, 876, L7
Sonnenfeld A., 2018, MNRAS, 474, 4648
Springel V., et al., 2005, Nature, 435, 629
Suyu S. H., Halkola A., 2010, A&A, 524, A94
Suyu S. H., Marshall P. J., Auger M. W., Hilbert S., Blandford
R. D., Koopmans L. V. E., Fassnacht C. D., Treu T., 2010,
ApJ, 711, 201
Suyu S. H., et al., 2012, ApJ, 750, 10
Suyu S. H., et al., 2013, ApJ, 766, 70
Suyu S. H., et al., 2014, ApJ, 788, L35
Suyu S. H., et al., 2017, MNRAS, 468, 2590
Suyu S. H., Chang T.-C., Courbin F., Okumura T., 2018, Space
Sci. Rev., 214, 91
Taubenberger S., et al., 2019, arXiv e-prints, p. arXiv:1905.12496
Tewes M., Courbin F., Meylan G., 2013a, A&A, 553, A120
Tewes M., et al., 2013b, A&A, 556, A22
Tie S. S., Kochanek C. S., 2018, MNRAS, 473, 80
Tihhonova O., et al., 2018, MNRAS, 477, 5657
Tonry J. L., 1998, AJ, 115, 1
Treu T., Koopmans L. V. E., 2002, MNRAS, 337, L6
Treu T., Marshall P. J., 2016, A&ARv, 24, 11
Treu T., Gavazzi R., Gorecki A., Marshall P. J., Koopmans
L. V. E., Bolton A. S., Moustakas L. A., Burles S., 2009,
ApJ, 690, 670
Treu T., et al., 2018, MNRAS, 481, 1041
Vanderriest C., Schneider J., Herpe G., Chevreton M., Moles M.,
Wlerick G., 1989, A&A, 215, 1
Vattis K., Koushiappas S. M., Loeb A., 2019, Phys. Rev. D, 99,
121302
Vuissoz C., et al., 2008, A&A, 488, 481
Weinberg D. H., Mortonson M. J., Eisenstein D. J., Hirata C.,
Riess A. G., Rozo E., 2013, Phys. Rep., 530, 87
Weymann R. J., Latham D., Angel J. R. P., Green R. F., Liebert
J. W., Turnshek D. A., Turnshek D. E., Tyson J. A., 1980,
Nature, 285, 641
Wilson M. L., Zabludoff A. I., Ammons S. M., Momcheva I. G.,
Williams K. A., Keeton C. R., 2016, ApJ, 833, 194
Wisotzki L., Schechter P. L., Bradt H. V., Heinmüller J., ReimersD., 2002, A&A, 395, 17
Wojtak R., Agnello A., 2019, MNRAS, 486, 5046
Wojtak R., Hjorth J., Gall C., 2019, MNRAS, 487, 3342
Wong K. C., et al., 2017, MNRAS, 465, 4895
Wong K. C., et al., 2018, ApJ, 867, 107
Wyman M., Rudd D. H., Vanderveld R. A., Hu W., 2014, Phys.
Rev. Lett., 112, 051302
Xu D., Sluse D., Schneider P., Springel V., Vogelsberger M., Nelson D., Hernquist L., 2016, MNRAS, 456, 739
Yıldırım A., Suyu S. H., Halkola A., 2019, arXiv e-prints, p.
arXiv:1904.07237

APPENDIX A: QUANTIFYING THE TRENDOF H0 WITH zd AND D∆t
There is an apparent trend of decreasing H0 inferred fromthe individual lenses as a function of increasing lens redshift,which is shown in the left panel of Figure A1. The significance of this correlation can be assessed by testing against
the null hypothesis, in which the measured H0 values are uncorrelated with the lens redshift. To do so, we draw sets of six
mock H0 values, using each lens’ own uncertainty probability distribution centered around the median joint inference−1
obtained in flat ΛCDM (H0 = 73.3+1.7
Mpc−1 ). We−1.8 km s
then fit a linear regression through each mock set. AssociatMNRAS 000, 1–?? (2019)

H0LiCOW XIII: A 2.4% measurement of H0 from lensed quasars
ing 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 istaken as our measurement of the correlation. We create 105
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 throughthe data, for which we find a negative slope. We find thatthe slope of the data falls 1.9σ away from the mock slope
distribution. We also observe a correlation between D∆t andH0 (Figure A1, right panel) that deviates from the null hypothesis at a similar significance level of 1.8σ.

MNRAS 000, 1–?? (2019)

21

22

K. C. Wong et al.

H0 [km s −1 Mpc −1 ]

90

B1608
RXJ1131
HE0435
J1206
WFI2033
PG1115

85
80
75
70
65
0.3

0.4

0.5

0.6

lens redshift

0.7

2000 3000 4000 5000 6000
D∆t [Mpc]

Figure A1. H0 constraints for the individual H0LiCOW lenses as a function of lens redshift (left) and time-delay distance (right). The
trend of smaller H0 value with increasing lens redshift and with increasing D∆t has significance levels of 1.9σ and 1.8σ, respectively.

MNRAS 000, 1–?? (2019)