Planck satellite, and BAO surveys provide strong support to the standard flat CDM cosmological model (e.g. Komatsu et al. 2011;Hinshaw et al. 2013;Planck Collaboration VI 2018). Under a few strong assumptions, such as flatness and constant dark energy density, these data give sub-per cent precision on the parameters of the standard cosmological model (e.g. Anderson et al. 2014;Kazin et al. 2014;Ross et al. 2015).

Intriguingly, distance measurements from Type Ia supernova (SN) that have been calibrated by the local distance ladder are smaller than the predictions from the CMB data given the flat CDM model (see the illustration in fig. 4 in Cuesta et al. 2015), leading to a ∼4.4σ tension in H 0 between the value predicted by the CMB and the local value (=74.03 ± 1.42 km s -1 Mpc -1 ; Riess et al. 2019). The SN data can also be calibrated by the inverse distance ladder method to yield a model-dependent value of H 0 . Under the assumption of the standard pre-recombination physics, combining BAO and SN with the CMB-calibrated physical scale of the sound horizon gives H 0 = 67.3 ± 1.1 km s -1 Mpc -1 (Aubourg et al. 2015). A recent blind analysis with additional SN data from Dark Energy Survey gives H 0 = 67.77 ± 1.30 km s -1 Mpc -1 (Macaulay et al. 2019). Both results are in excellent agreement with the Planck value (H 0 = 67.27 ± 0.6 km s -1 Mpc -1 ) under the assumption of flat CDM model (Planck Collaboration VI 2018). Furthermore, even without using the CMB anisotropy, the combination of BAO data with light element abundances produces Planck-like H 0 values (Addison et al. 2018). This indicates that systematic errors, especially in the Planck data analysis, most likely are not the main driver of the H 0 discrepancies. Similarly, the local distance ladder analyses have also passed a range of systematic checks (e.g. Efstathiou 2014;Cardona, Kunz & Pettorino 2017;Zhang et al. 2017;Dhawan, Jha & Leibundgut 2018;Feeney, Mortlock & Dalmasso 2018;Follin & Knox 2018;Riess et al. 2019), and rule out the local void scenario (Keenan, Barger & Cowie 2013;Fleury, Clarkson & Maartens 2017;Kenworthy, Scolnic & Riess 2019;Shanks, Hogarth & Metcalfe 2019) There have been several attempts to address this ∼4.4σ tension by extending the standard cosmological model, either by changing the size of the sound horizon in the early Universe (e.g. Heavens, Jimenez & Verde 2014; Wyman et al. 2014;Cuesta et al. 2015;Alam et al. 2017;Agrawal et al. 2019;Kreisch, Cyr-Racine & Doré 2019;Poulin et al. 2019) or by altering the expansion history (Efstathiou 2003;Linder 2004;Moresco et al. 2016;Alam et al. 2017).

Recent studies (e.g. Bernal, Verde & Riess 2016;Joudaki et al. 2018;Aylor et al. 2019;Lemos et al. 2019) have also tried to directly reconstruct H(z) in order to investigate the H 0 tension in the context of a possibly poor understanding of the evolution of dark energy density. From an empirical point of view, the current SN and BAO data sets only support w(z) = -1 within the redshift range where data are available (Cuesta et al. 2016). A very recent and dramatic decrease in w or the presence of strong dark energy at 3 < z < 1000 may escape detection and still generate a high value of H 0 (Riess et al. 2016). Standard sirens could possibly explore the z > 3 range in the future and provide a high-precision H 0 measurement (Chen, Fishbach & Holz 2018b). Nevertheless, it is also important to note that some H 0 -value tension remains even if we do not consider the distance ladder constraints. For example, the high-CMB power spectrum prefers an even lower H 0 value than that from the low-CMB power spectrum (Addison et al. 2016;Planck Collaboration VI 2018).

Given the various tensions across different data sets, any convincing resolution to the H 0 tensions, either due to unknown systematics or due to new physics, needs to simultaneously resolve multiple disagreements. Therefore, comparing the distance measurements among independent and robust methodologies to cross-examine the H 0 tension is probably the only way to shed light on the true answer.

Time-delay cosmography not only is a completely independent technique of the distance ladder methods, but also has the advantage of being a one-step measurement of combined cosmological distances.

In addition, time-delay cosmography is a complementary and costeffective alternative compared to Type Ia SN or BAO (Suyu et al. 2013;Tewes et al. 2013b). In a time-delay gravitational lens, the combined cosmological distance we can measure is called the timedelay distance (Suyu et al. 2010), which is a ratio of the angular diameter distances in the system

where D is the distance to the lensing galaxy, D s is the distance to the background source, and D s is the distance between the lens and the source. Furthermore, we can make a separate determination of D by measuring the velocity dispersion of the lensing galaxy (Jee, Komatsu & Suyu 2015;Birrer, Amara & Refregier 2016;Jee et al. 2016;Birrer et al. 2019). First proposed by Refsdal (1964), the H 0 measurement requires modelling the mass in the lensing galaxy and along the line of sight, and measuring the time delays between multiple images via a monitoring program. The advantage of this method is that D t is primarily sensitive to H 0 and insensitive to the neutrino physics and spatial curvature, but still sensitive to the properties of dark energy (Bonvin et al. 2017;Birrer et al. 2019). The H0LiCOW collaboration 1 is using strong gravitational lens systems to measure cosmological parameters (Suyu et al. 2017). The most recent measurement of H 0 from the collaboration used the doubly lensed quasar system, SDSS J1206+4332, to derive H 0 = 68.8 +5.4 -5.1 km s -1 Mpc -1 for that lens system alone, as well as combining the new system with previous H0LiCOW lenses (Suyu et al. 2009(Suyu et al. , 2010(Suyu et al. , 2013(Suyu et al. , 2014;;Bonvin et al. 2017;Rusu et al. 2017;Sluse et al. 2017;Wong et al. 2017) to obtain a joint inference on H 0 with 3 per cent precision: H 0 = 72.5 +2.1 -2.3 km s -1 Mpc -1 (Birrer et al. 2019). This result agrees with the Riess et al. (2019) value within the 1σ uncertainties. More recently, the collaboration completed its analysis of WFI2033-4723 using the time delays from COSMOGRAIL 2 (Bonvin et al. 2019;Rusu et al. 2019a;Sluse et al. 2019).

Achieving the goal of obtaining a 1 per cent or better measurement of H 0 with time-delay cosmography requires a significantly larger sample of lensed quasars with high-quality data than has been analysed to date. Many new lensed quasars have been discovered (e.g. Lin et al. 2017;Schechter et al. 2017;Agnello et al. 2018; et al. 2014;Agnello 2017;Ostrovski et al. 2017;Petrillo et al. 2017;Lanusse et al. 2018;Spiniello et al. 2018;Treu et al. 2018;Avestruz et al. 2019), and more are expected to be found with the Large Synoptic Survey Telescope (Oguri & Marshall 2010). Hence, a 1 per cent H 0 measurement from time-delay cosmography is a realistic expectation in the near future (e.g. Jee et al. 2015Jee et al. , 2016Jee et al. , 2019;;de Grijs et al. 2017;Shajib, Treu & Agnello 2018;Suyu et al. 2018) if we can control the systematic effects to a sub-per cent level (Dobler et al. 2015;Liao et al. 2015;Ding et al. 2018).

A critical input to achieving precise modelling of the mass distribution in the lensing galaxy is sensitive high-resolution imaging in which extended emission from the background source is detected. For this reason, many models of the lensing potential are based on imaging from Hubble Space Telescope (HST) (e.g. Suyu et al. 2009Suyu et al. , 2010;;Birrer, Amara & Refregier 2015;Birrer et al. 2016Birrer et al. , 2019;;Wong et al. 2017). Adaptive optics (AO) observations with groundbased telescopes can provide imaging with angular resolution that is comparable to or better than HST data, especially for systems that are faint in the optical and bright at near-infrared wavelengths, thus providing an attractive alternative for modelling the lens mass distribution (e.g. Lagattuta et al. 2012;Chen et al. 2016). However, the challenge of using AO data is the unstable point spread function (PSF). Chen et al. (2016) showed that with a new iterative PSF reconstruction method applied to RXJ 1131-1231, the reconstructed PSF allows one to model the AO imaging down to the noise level, as well as providing tighter constraints on the lens model than were obtained from HST imaging of the system. In this new analysis, we apply the PSF reconstruction method to three lens systems, HE 0435-1223, PG 1115+080, and RXJ 1131-1231, which have not only high-resolution AO imaging and HST imaging, but also measured time delays, stellar velocity dispersions, and studies of their environments. We then infer their time-delay distances via detailed lens mass modelling. The AO data for RXJ 1131-1231 have already been analysed by Chen et al. (2016), but the analyses based on the AO imaging of HE 0435-1223 and PG 1115+080 are presented here for the first time. We show the galaxies/group that we explicitly model here. As PG 1115+080 is embedded in a nearby group that consists of 13 galaxies labelled with solid circles, we model not only the main lens but also the group explicitly. The dotted circle represents 1σ uncertainty of the priors of the group position ( RA = 23. 4 and Dec. = 15. 84) measured in Wilson et al. (2016). Since G1 and G2 have the first two largest values of 3 x, we model either G1 or both G1 and G2 explicitly in addition to the main lens. We label G1 and G2 with the dashed circles.

The Keck AO imaging data are part of the Strong-lensing High Angular Resolution Programme (SHARP; Fassnacht et al., in preparation), which aims to study the nature of dark matter using high-resolution AO imaging (e.g. Lagattuta, Auger & Fassnacht 2010;Lagattuta et al. 2012;Vegetti et al. 2012;Hsueh et al. 2016Hsueh et al. , 2017Hsueh et al. , 2018;;Spingola et al. 2018). The time-delay measurements are provided by the COSMOGRAIL group (e.g. Courbin et al. 2005;Vuissoz et al. 2007Vuissoz et al. , 2008;;Courbin et al. 2011;Rathna Kumar et al. 2013;Tewes, Courbin & Meylan 2013a;Tewes et al. 2013b;Bonvin et al. 2017Bonvin et al. , 2018)), which aims to provide the highestprecision measurements of time delays. The lens environment of RXJ 1131-1231 and HE 0435-1223 studies is provided by the H0LiCOW team (Suyu et al. 2014;Rusu et al. 2017;Sluse et al. 2017) The outline of the paper is as follows. In Section 2, we briefly recap the basics of obtaining inferences on cosmography with time-delay lenses and the statistical tools we use for this process. We describe the observations of HE 0435-1223, PG 1115+080, and RXJ 1131-1231 with the AO imaging system at the Keck Observatory in Section 3. In Section 4, we describe the models that we use to analyse the data. In Section 5, we elaborate the detailed lens modelling and the properties of the reconstructed PSF for each system. In Section 6, we present the joint cosmological inference from the AO imaging only as well as from combined AO plus HST. We summarize in Section 7.

We briefly introduce the relation between cosmology and gravitational lensing in Section 2.1 and the joint inference of all information in Section 2.2.

When a compact variable background source, such as an active galactic nucleus (AGN) or SN, sitting inside its host galaxy is strongly lensed by a foreground object, the distorted host galaxy shape combined with the time delay between the multiple images allows one to precisely determine a particular size of the system. One can express the excess time delays as

where θ, β, and ψ(θ) are the image location, the source location, and the projected 2D lensing gravitational potential, respectively (Refsdal 1964;Shapiro 1964). The λ parameter represents the lack of perfect knowledge of the full mass distribution, as discussed later.

The advantage of this formulation is the separability of the cosmographic information, contained in the D t parameter, and the lens modelling. This allows one to infer cosmographic information without the need for cosmological priors in the lens modelling. Although for the lens systems that contain multiple lensing galaxies at different redshifts a particular cosmological model needs to be applied to the analysis, the resulting effective D t is robust against cosmological model assumptions, in that the posterior distribution of the effective D t shows nearly identical distribution (within 1 per cent) for different background cosmological models (Wong et al. 2017;Rusu et al. 2019a).

Because of the mass-sheet transformation (MST; Falco, Gorenstein & Shapiro 1985;Schneider & Sluse 2013, 2014), the determination of D t is subject to an understanding of λ (see also the discussion in Birrer et al. 2019). Thus, additional priors and information from simulations, environmental data, and the stellar velocity dispersion of the lensing galaxy are required to constrain the degeneracy between different mass profiles and the degeneracy between the mass profile and the mass sheet contributed by the environment (Suyu et al. 2010;Fassnacht, Koopmans & Wong 2011;Rusu et al. 2017;Tihhonova et al. 2018). In addition, a prior on the source size (Birrer et al. 2016) or having a background source of known brightness (e.g. a lensed SN; Grillo et al. 2018) can also put constraints on λ. We refer interested readers to Treu & Marshall (2016) and Suyu et al. (2018) for more details.

In this work, we will present our joint inference on D t . We use d i to denote the imaging data, where i = HE, PG, and RXJ represent HE 0435-1223, PG 1115+080, and RXJ 1131-1231, respectively; we use t i for time delays, d ENV i to characterize the lens environments, and σ i for the stellar velocity dispersions of the lensing galaxies. Here, η i are the parameters we want to infer from the data, and A denotes the discrete assumptions that we made in the models (e.g. whether the lensing galaxy is modelled using a power-law or NFW + stellar mass distribution). The posterior of η i can be expressed as

where P (d i , t i , σ i , d ENV i |η i , A i ) is the joint likelihood for each lens. Since we assume that the environment can be decoupled from the lens, and that the data sets are independent,

Note that we do not use the flux ratios of the lensed quasars as constraints on the model.

In order to explore the unmodelled systematic uncertainties that may arise from modelling choices, we vary the content of A for each lens. The marginalized integral can be expressed as

where A i,k and d i,tot represent the different model choices and all data sets for the lens system i, respectively. For ranking the models, we follow Birrer et al. (2019) to estimate the evidence, P ( A i,k ), by using the Bayesian information criterion (BIC), which is defined as

where n is the number of data points including the lens imaging, eight AGN positions, three time delays, and one velocity dispersion, k is the number of free parameters in the lens model that are given uniform priors, plus two source-position parameters, plus one anisotropy radius to predict the velocity dispersion, and L is the maximum likelihood of the model, which is the product of the AGN position likelihood, the time-delay likelihood, the pixelated image plane likelihood, and the kinematic likelihood. The image plane likelihood is the Bayesian evidence of the pixelated-source intensity reconstruction using the arcmask imaging data (see Suyu & Halkola 2010) times the likelihood of the lens model parameters within the image plane region that excludes the arcmask. We follow Birrer et al. (2019) and calculate the relative BIC and weighting for the SPEMD and composite models separately to avoid biases due to our choice of lens model parameterization.

The analysis in this paper is based on new Keck AO and archival HST observations of three gravitational lens systems. In this section, we describe the lens sample and the data acquisition and analysis.

The sample consists of three well-known lensed quasar systems. Images of these lens systems are shown in Fig. 1.

(i) HE 0435-1223: The HE 0435-1223 system (J2000: 4 h 38 m 14. s 9, 12 • 17 14. 4) is a quadruply lensed quasar discovered by Wisotzki et al. (2002). The main lensing galaxy is at a redshift of z = 0.4546 (Morgan et al. 2005), and the source redshift is z s = 1.693 (Sluse et al. 2012). The lens resides inside a galaxy group that contains at least 12 galaxies, with a velocity dispersion σ = 471 ± 100 km s -1 (e.g. Momcheva et al. 2006;Wong et al. 2011;Wilson et al. 2016;Sluse et al. 2017). Wong et al. (2017) measured the stellar velocity dispersion of the lensing galaxy to be σ = 222 ± 15 km s -1 . The time delays of this system were measured by Bonvin et al. (2017) with ∼6.5 per cent uncertainties.

(ii) PG 1115+080: This four-image system was the second strong gravitational lens system to be discovered (Weymann et al. 1980). The system is located at 11 h 18 m 16. s 899, +7 • 45 58. 502 (J2000). The background quasar with a redshift of z s = 1.722 is lensed by a galaxy with z = 0.3098 (Henry & Heasley 1986;Christian, Crabtree & Waddell 1987;Tonry 1998). The lensed images are in a classic 'fold' configuration, with an image pair A1 and A2 near the critical curve. The lens resides inside a galaxy group with 13 known members that has a velocity dispersion of σ = 390 ± 60 km s -1 (Wilson et al. 2016). Tonry (1998) measured the stellar velocity dispersion of the lensing galaxy to be σ = 281 ± 25 km s -1 . The time delays of this system were measured by Bonvin et al. (2018) with ∼6.4 per cent uncertainties.

(iii) RXJ 1131-1231: The RXJ 1131-1231 system (J2000: 11 h 31 m 52 s , -12 • 31 59 ) is a quadruply lensed quasar discovered by Sluse et al. (2003). The spectroscopic redshifts of the lensing galaxy and the background source are at z = 0.295 (Suyu et al. 2013) and z s = 0.657 (Sluse et al. 2007), respectively. Suyu et al. (2013) measured the stellar velocity dispersion of the lensing galaxy to be σ = 323 ± 20 km s -1 . The time delays were measured by Tewes et al. (2013b) with ∼1.5 per cent uncertainties.

All three lens systems were observed at K -band with the Near-Infrared Camera 2 (NIRC2), sitting behind the AO bench on the Keck II Telescope, as part of the SHARP AO effort (Fassnacht et al., in preparation). The targets were observed with either the narrow camera set-up, which provides a roughly 10 arcsec × 10 arcsec field of view and a pixel scale of 9.942 mas, or the wide camera that gives a roughly 40 arcsec × 40 arcsec field of view and a pixel scale of 39.686 mas. Details of the observations are provided in Table 1.

The NIRC2 data were reduced using the SHARP PYTHON-based pipeline, which performs a flat-field correction, sky subtraction, correction of the optical distortion in the images, and a co-addition of the exposures. During the distortion correction step, the images are resampled to produce final pixel scales of 10 mas pixel -1 for the narrow camera and 40 mas pixel -1 for the wide camera. The narrow camera pixels oversample the PSF, which has typical full width at half-maximum (FWHM) values of 60-90 mas. Therefore, to improve the modelling efficiency for the narrow camera data, we perform a 2 × 2 binning of the images produced by the pipeline to obtain images that have a 20 mas pixel -1 scale. Further details on RXJ 1131-1231, which was observed with the NIRC2 wide camera, can be found in Chen et al. (2016).

All three lens systems have been observed by HST (GO-9375, PI: Kochanek; GO-9744, PI: Kochanek; GO-12889, PI: Suyu). The HST imaging of both RXJ 1131-1231 (Suyu et al. 2013) and HE 0435-1223 (Wong et al. 2017) was analysed in the previous work. Therefore, the inferences from these previous models are combined with those from the new AO models in Section 5. In contrast, the HST data for PG 1115+080 have not been modelled using the latest pixelated techniques, although Treu & Koopmans (2002) combined lensing geometry and velocity dispersion to study the content of the luminous matter and dark matter profiles. Therefore, we perform a joint modelling procedure on the AO and HST data for PG 1115+080.

The contribution of the line-of-sight (LOS) mass distribution to the lensing requires deep wide-field imaging of the region surrounding the lens system. For HE 0435-1223 and RXJ 1131-1231, we have used HST/ACS or Subaru SuprimeCam imaging, which have been analysed as part of our previous work on these systems (Suyu et al. 2014;Rusu et al. 2017). To achieve the requisite combination of depth and area for PG 1115+080, we co-added 95 images of the field taken with the Wide Field Imager (Baade et al. 1999) mounted at the Cassegrain focus of the MPIA 2.2 m telescope. The camera provides a pixel scale of 0.238 arcsec. The data were obtained through ESO BB#R c /162 filter. These are the same data used by COSMOGRAIL to measure the time delay for this system (see section 2.1 of Bonvin et al. 2018, for details); each image has been exposed for 330 s and covers a field of view of ∼8 arcmin × 16 arcmin. The data reduction process follows the standard procedure, including master bias subtraction, master flat fielding, sky subtraction, fringe pattern removal, and finally exposure-to-exposure normalization using SEXTRACTOR on field stars prior to co-adding the exposures. The final co-added image has an effective seeing of 0.86 arcsec.

In this section, we describe the models that we use for fitting the high-resolution imaging data, including the lens mass models in Section 4.1, lens light models in Section 4.2, the models for constraining the MST in Section 4.3, and the time-delay prediction models in Section 4.4. We use GLEE, a strong lens modelling code developed by S. H. Suyu and A. Halkola to model the lens arc, lens light, and lens AGNs simultaneously (Suyu & Halkola 2010;Suyu et al. 2012). The source galaxy is reconstructed on a regular grid, and we choose curvature as our regularization form (Suyu et al. 2006). The AO PSF is reconstructed by the method developed in Chen et al. (2016), which uses the four lensed quasars as prior information to iteratively reconstruct the AO PSF from the imaging data. The tests in Chen et al. (2016) showed that the PSF reconstruction method can recover the structure of the host galaxy and PSF in mock data sets. C1. The distributions with G1 or G2 marked are calculated by removing those galaxies from the weighted count constraints, since these galaxies are explicitly included in the lens models. The size of the histogram bin is κ ext = 0.00 055. As the original distributions are noisy, we plot their convolution with a large smoothing window of size 30 × κ ext . In the legend, 'pow' refers to the power-law model, 'com' refers to the composite model, κ refers to the median of the distribution, and σ κ refers to the semi-difference of the 84 and 16 percentiles of the distribution.

The following three analytical functions are used for modelling the main lens, nearby groups, and nearby galaxies:

(i) SPEMD: Many studies have shown that a power-law model provides a good first-order description of the lensing galaxies for galaxy-galaxy lensing (e.g. Koopmans et al. 2006Koopmans et al. , 2009;;Suyu et al. 2009;Auger et al. 2010;Barnabè et al. 2011;Sonnenfeld et al. 2013). Thus, for every lens, we model the mass distribution of the lensing galaxy with a singular power-law elliptical mass distribution (Barkana 1998). The main parameters include radial slope (γ ), Einstein radius (θ E ), and the axial ratio of the elliptical isodensity contour (q).

(ii) Composite: We follow Suyu et al. (2014) and test a composite (baryonic + dark matter) model. The baryonic component is modelled by multiplying the lens surface brightness distribution by a constant M/L ratio parameter (see Section 4.2). For the dark matter component, we adopt the standard NFW profile (Navarro, Frenk & White 1996;Golse & Kneib 2002) with the following parameters: halo normalization (κ s ), halo scale radius (r s ), and halo minor-to-major axial ratio (q), as well as associated position angle (θ q ). Note that the ellipticity is implemented in the potential for the dark matter. The ellipticities of our three lenses are all inside the range where elliptical potential models are a good description of elliptical mass distributions and thus can be reliably applied to observational data (Golse & Kneib 2002).

(iii) SIS: Singular isothermal sphere (SIS) models are used to describe the nearby group and the individual galaxies inside the group of PG 1115+080, the nearby galaxies of HE 0435-1223, and the satellite of RXJ 1131-1231.

The following two analytical functions are used to model the lens light distribution:

(i) 2Sérsic: We model the light distribution of the lens galaxy with two concentric elliptical Sérsic profiles. For all three lens systems, we found that a single Sérsic profile was insufficient for modelling the light distributions.

(ii) 2Chameleon: The chameleon profile is the difference of two isothermal profiles. It mimics a Sérsic profile and enables computationally efficient lens modelling (Dutton et al. 2011). The parametrized Chameleon profile can be found in Suyu et al. (2014). We convert the Chameleon light profile to mass with an additional constant M/L ratio parameter when modelling the composite model described in Section 4.1.

The MST is a known degeneracy in lens modelling, in which one can transform a projected mass distribution, κ(θ), into infinite sets of κ λ (θ) via

without degrading the fit to the imaging. The corresponding timedelay distance changes via

This degeneracy can be produced by both the lens environment and an incorrect description of the mass distribution in the lensing galaxy. We discuss the approach for estimating the contribution from the environment in Section 4.3.1, and for ranking the mass models with kinematic information in Section 4.3.2.

If we perfectly know the true κ(θ), then the role of λ in equation ( 7) can be understood by looking at the behaviour far from the lensing galaxy, i.e. as θ → ∞. In this regime, the mass distribution of the lensing galaxy, κ(θ), approaches 0 and, thus, λ can be interpreted as a constant-density mass sheet contributed by the lens environment.

It is exactly the first-order form produced by LOS structure when its effect is small. Therefore, in our models we identify λ with κ ext , the physical convergence associated with LOS structures that do not affect the kinematics of the strong lens galaxy (e.g. Suyu et al. 2010;Wong et al. 2017;Birrer et al. 2019) In addition, the second-order distortion from the LOS also produces a tidal stretching on the lens images.

We use the following models to capture these two effects.

(i) Shear: Shear distorts the lensed image shapes and thus it can be detected in the modelling process. We express the lens potential in polar coordinates (θ, ϕ) to model the external shear on the imaging plane

where γ ext is the shear strength and ϕ ext is the shear angle. The shear position angle of ϕ ext = 0 • corresponds to a shearing along θ 1 , whereas ϕ ext = 90 • corresponds to shearing along θ 2 .foot_1 (ii) Millennium Simulation: Due to the MST, the lens images do not provide direct information on κ ext . We thus use the results of ray tracing by Hilbert et al. ( 2009) through the Millennium Simulation (Springel et al. 2005) to statistically estimate the mass contribution along the line of sight to our lenses. This technique was first employed by Suyu et al. (2010), who took the ratio of observed galaxy number counts in an aperture around a lens to those in a control survey (from Fassnacht et al. 2011), in order to measure the local over/underdensity of galaxies in the lens fields. They then selected lines of sight of similar over/underdensity from the Millennium Simulation, with their corresponding values of the convergence, thus producing a probability distribution for κ ext : P (κ ext |d ENV ). Suyu et al. (2013) later used, in addition to the number counts, the shear value inferred from lens modelling as an additional constraint on κ ext . Greene et al. (2013) showed that further constraints, in the form of weighted number counts, can be derived by incorporating physical quantities relevant to lensing such as the distance of each galaxy to the lens, redshifts, luminosities, and stellar masses. Birrer et al. (2019) expressed the technique as an application of approximate Bayesian computing and further combined weighted number count constraints from multiple aperture radii. In Sections 5.1.5, 5.2.3, and 5.3.2, we provide an implementation of this technique for each of our three lenses, customized to the nature of the available environment data. For the main numerical and mathematical details of the implementation, we refer the reader to Rusu et al. (2017).

If the mass along the line of sight is large enough such that we cannot ignore the higher-order terms (flexion and beyond), we need to model that mass explicitly. McCully et al. (2014McCully et al. ( , 2017) ) give a quantitative term, flexion shift ( 3 x), which estimates the deviations in lensed image positions due to third-order (flexion) terms, and suggest that if 3 x is higher than 10 -4 arcsec (for the typical galaxy-scale lenses that we are studying here), one should model the perturbers explicitly to avoid biasing H 0 at the (sub-) per cent level. We follow this convention to choose which galaxies are included in our models. Note that the 3 x threshold is based on mock data that do not include extended arcs. Therefore, this criterion should already be conservative since real lens imaging with extended arcs provides more constraining power than the point-source imaging that was used to set the threshold.

Conversely, even if we perfectly know the mass along the line of sight, the value of λ remains uncertain because we do not know the true κ(θ) distribution. As pointed out by Schneider & Sluse (2013), when assuming a specific mass profile, one artificially breaks the internal MST (Koopmans 2004), a special case of the source-position transformation (Schneider & Sluse 2014;Wertz & Orthen 2018;Wertz, Orthen & Schneider 2018). However, recent work based on the Illustris simulation has indicated that this effect may be of less concern for massive galaxies such as those in the H0LiCOW sample (Xu et al. 2016). The MST allows many different mass distributions within the Einstein radius of the lens, as long as the integrated κ within the Einstein ring is preserved. However, the stellar velocity dispersion is sensitive to the integrated value of κ within the effective radius of the lensing galaxy, which is often different from the Einstein radius. Thus, we can use the observed stellar velocity dispersion of the lensing galaxy to rank different mass models (Treu & Koopmans 2002, 2004).

In the lens modelling of RXJ 1131-1231, Suyu et al. (2014) have shown that by including the velocity dispersion, one can obtain a robust D t when considering both power-law and composite models. Sonnenfeld (2018) also shows that velocity dispersion is the key to obtaining an unbiased H 0 measurement. Hence, we follow Suyu et al. (2014) in adopting the composite model and also incorporating the velocity dispersion into the modelling to mitigate this internal mass-profile degeneracy.

There are three different components needed to predict the velocity dispersion.

(i) A 3D mass distribution: Following Suyu et al. (2010), one can obtain the 3D lens mass from the lens modelling by assuming spherical symmetry. 4 In general, the spherically symmetric 3D mass density of the lens can be expressed as

where ρ 0 r n 0 and F n (r) are the normalization and the mass density distribution. By integrating ρ local within a cylinder with radius given by the Einstein radius, R Ein , one obtains where M 2D (R Ein ) is the projected mass within R Ein . The mass contained in M local is

where M ext represents the mass contribution from κ ext and cr = c 2 4πG

is the critical surface mass density. Combining equation ( 12) with equation ( 13), the normalization in equation ( 11) can be expressed as

Substituting this in equation ( 11), we obtain

Although there is (1κ ext ) in equation ( 16), the normalization of the local mass density distribution remains invariant (Yıldırım et al. 2019) as cr can be re-expressed as cr = c 2 4πG

where (1κ ext ) term cancels out in equation ( 16).

(ii) An anisotropy component: We assume the anisotropy component in the form of an anisotropy radius, r ani , in the Osipkov-Merritt (OM) formulation 5 (Osipkov 1979;Merritt 1985)

where r ani = 0 is pure radial orbits and r ani → ∞ is isotropic with equal radial and tangential velocity dispersions.

(iii) A stellar component: We assume a Hernquist profile 6 (Hernquist 1990)

for the power-law model, where I 0 is the normalization term and the scale radius can be related to the effective radius by a = 0.551r eff . 5 We further tested an additional anisotropy model, namely, a two-parameter extension of OM, and found that the uncertainty on the inferred D due to anisotropy models is comparable to that from the choice of different mass models (i.e. power-law versus composite models). Thus, we should not be underestimating the uncertainty on D (Jee et al. 2019). In addition, different anisotropy models have negligible impact on the D t measurement. 6 Suyu et al. (2010) showed that different types of stellar distribution functions, namely, Hernquist versus Jaffe, produced nearly identical PDFs for the cosmological parameters.

For the composite model, the stellar component is represented by the light profile multiplied by a constant mass-to-light ratio.

With the aforementioned three components, we follow Sonnenfeld et al. ( 2012) and calculate the 3D radial velocity dispersion by numerically integrating the solutions of the spherical Jeans equation (Binney & Tremaine 1987)

given the κ ext from Section 4.3.1. Note that since the LOS velocity dispersion has a degeneracy between its anisotropy and the mass profile (Dejonghe 1987), we marginalize the sample of r ani over a uniform distribution [0.5, 5]r eff . To compare with the data, we can get the seeing-convolved luminosity-weighted LOS velocity dispersion

where R is the projected radius, I(R) is the light distribution, P is the PSF convolution kernel (Mamon & Łokas 2005), and A is the aperture. The luminosity-weighted LOS velocity dispersion is given C1. For the powerlaw + G1 and composite + G1 models, we used an inner mask of 5 arcsec around the lens, and for the powerlaw + 5 perturbers model, we used a mask of 12 arcsec radius. The five perturbers are indicated in fig. 3 from Wong et al. (2017). Three additional galaxies enter the 12 arcsec radius inner mask, and we slightly boosted their distance from the lens, in order to avoid masking them. See caption of Fig. 4 for additional details.

by

For a system with significant perturbers at a different redshift from the main lens (e.g. HE 0435-1223), we assume a flat CDM cosmology with H 0 uniform in [0, 150] km s -1 Mpc -1 , m = 0.3, and m = 1to calculate the critical density and rank the models by the predicted velocity dispersion. Note that this assumption does not affect the generality of the conclusion. For our single lens plane systems (i.e. RXJ 1131-1231 and PG 1115+080), we use the measured velocity dispersion to constrain D s /D s and then combine with the measurement of D t to infer the value of D without assuming any cosmological model (Birrer et al. 2016(Birrer et al. , 2019)). The further advantage of this method is that D is not affected by κ ext (Jee et al. 2015).

In Section 2, we showed that the time delays between multiple images are due to the geometry and the gravitational potential that the light passes through. Tie & Kochanek (2018) introduce a possible new microlensing effect on the time delays that can shift the light curves depending on the structure of the accretion disc in the lensed quasar and the density of the stars in the lensing galaxy. They estimated this effect under the assumption of a lamppost model for the accretion disc, where a large part of the disc lights up concurrently on light-travel scales that are a significant fraction of the time delays between the lensed images. The observed time delay could thus be affected by differential magnification of the individual images, resulting from microlensing by stars in the lensing galaxy. However, the lamp-post model is only one choice for how to represent the accretion disc; other accretion disc models for which variability is different from the lamp-post model are possible (e.g. Dexter & Agol 2011). We follow Chen et al. (2018a) and present the D t measurements both with and without the lamppost assumption. However, we only consider the case without the microlensing effect in our final H 0 determination since it is not clear at this point which is the proper disc model to use. Note that it also was not applied in previous H0LiCOW work to infer the final H 0 measurement.

A more detailed description of this effect and how to estimate the probability distribution of the microlensing time-delay effect (MTDE) can be found in Bonvin et al. (2018). We briefly summarize the technique here. We generate magnification maps using GPU-D (Vernardos & Fluke 2014), which incorporates a graphics processing unit implementation of the inverse ray-shooting technique (Kayser, Refsdal & Stabell 1986). All magnification maps have dimension of 8192 × 8192 pixels over a scale of 20 R E , where

We choose the Salpeter initial mass function with mean mass M = 0.3M and the ratio between the upper and lower masses M upper /M lower = 100 (Kochanek 2004). We consider a standard thin disc model (Shakura & Sunyaev 1973). Given the disc size of lensed quasar, the average microlensing time delay at each position on a magnification map can be derived using equation (10) of Tie & Kochanek (2018). The parameters that are used to estimate the probability distribution of the microlensing time delay for each system are listed in Tables 2 and3.

To fold this effect into time-delay modelling, we use

where the first term on the right-hand side is the same as in equation ( 2) and t it j is the extra delay caused by the MTDE between images i and j (see details in Chen et al. 2018a).

Both HE 0435-1223 and RXJ 1131-1231 have been extensively modelled using the extended lensed emission seen in highresolution HST imaging of the systems (Suyu et al. 2014;Wong et al. 2017), but our modelling techniques have not yet been applied to PG 1115+080. Therefore, in this section we begin with a description of the modelling of PG 1115+080 in Section 5.1, and then describe HE 0435-1223 in Section 5.2 and RXJ 1131-1231 in Section 5.3.

For PG 1115+080, we model the HST and AO imaging simultaneously. However, for HE 0435-1223 and RXJ 1131-1231, we only model the AO imaging since the HST imaging has already been modelled (Suyu et al. 2014;Wong et al. 2017), and then combine the two modelling outputs to obtain a joint inference on H 0 (see Section 6.2).

Our analyses of PG 1115+080 and HE 0435-1223 are blind, as in Suyu et al. (2013) and Rusu et al. (2019a), in order to avoid confirmation bias. That is, the values of D t , D (if computed), and H 0 were kept blind until all co-authors came to a consensus to reveal the values during a collaboration telecon on June 5. The analysis was frozen after we unblinded the results and no changes were made to any of the numerical results. The time between unblinding and submission was used to polish the text and figures of the manuscript, and carrying out the detailed comparison of the AO-and HST-based analysis.

In contrast, the RXJ 1131-1231 analysis was not done blindly as the AO data for this system were used to develop the PSF reconstruction technique. On top of the power-law model we have done in Chen et al. (2016), we further test the composite model and use both models to infer D t and D .

To better control the systematics due to the choice of lens modelling technique, we run Markov chain Monte Carlo (MCMC) sampling with different source resolutions in each model. This approach was used because Suyu et al. (2013) have shown that the effects of the pixelated-source grid resolution dominate the uncertainty on the lens modelling when using a modelling code, such as GLEE, that implements the pixelated-source reconstruction technique.

PG 1115+080 is a single-plane lens system embedded in a nearby group that consists of 13 known galaxies (see the solid circles in Fig. 2). If we model the lens without including the group, the mass profile shows a very steep slope (γ ∼ 2.35; note that γ = 2 corresponds to the isothermal profile), which has also been found in previous studies of this system (e.g. Keeton & Kochanek 1997;Treu & Koopmans 2002), and a strong shear (γ ext ∼ 0.15), which comes from the nearby group. Wilson et al. (2016) showed that, compared with the other 11 groups along the light of sight, the nearby group contributes the largest convergence at the lens position. Furthermore, McCully et al. (2017) indicate that the group produces a significant flexion shift. Thus, it is crucial to model not only the main lens but also the group explicitly if we want to obtain an unbiased H 0 measurement.

For the HST imaging, we use TINYTIM (Krist & Hook 1997) to generate the PSFs with different spectral index, α, of a power law from -0.4 to -2.5 and different focusesfoot_3 from 0 to 10. We find that the best fit is the PSF with focus equal to 0 and spectral index equal to -1.6. We use this TINYTIM PSF as the initial guess and then apply the PSF correction method while modelling the HST imaging. For the AO imaging, we follow the criteria described in section 4.4.3 in Chen et al. (2016) and perform eight iterative steps to create the final PSF and make sure the size of the PSF for convolution is large enough so that the results are stable. The FWHM of the AO PSF is 0.07 arcsec, while the FWHM of the HST PSF is 0.15 arcsec. We show the reconstructed AO PSF in Fig. A1.

We follow two approaches to modelling the mass distribution in the main lensing galaxy.

(i) SPEMD + 2Sérsic + shear: We first choose the SPEMD density profile to model the extended arc and reconstruct the source structure on a pixelated grid (Suyu et al. 2006). We found that a single Sérsic profile is not sufficient to describe the light distribution, so we model it with two concentric elliptical Sérsic profiles with free relative position angles and ellipticities. By comparing the mass and light components, we found a similar result to Yoo et al. (2005), namely that the position of the centre of mass is very close to the centre of light, with | r| ≈ 0.015 arcsec. This implies that the offset between the projected centre of dark matter and baryonic matter is small.

(ii) Composite + 2Chameleon + shear: We also model the main lens with composite model. Because the SPEMD + 2Sérsic + shear model indicated that the dark matter and baryonic centroids were consistent, we link the centroid of NFW profile to the centroid of two concentric chameleon profiles. We also follow Wong et al. (2017) and Rusu et al. (2019a) to iteratively update the relative amplitudes of the associated mass components to match those of the light components, as the relative amplitudes of light components can vary when we run the MCMC chains, while the relative amplitudes are fixed in the mass profiles.

Based on the velocity dispersion of the nearby group, σ group = 390 ± 60 km s -1 , the inferred group mass is around 10 13 -10 14 h -1 M (Wilson et al. 2016). Oguri (2006) has shown that in this mass regime the mass profile is too complicated to be described by either a simple NFW profile or SIS profile, as it is a transition between the two. Thus, we use an NFW profile as the fiducial model, but also model the group with an SIS profile as a systematic check. In the following, we show how we determine reasonable priors on the NFW and SIS profiles. 2008) to translate the observed velocity dispersion to scale radius and the normalization of the NFW profile (we compared the priors by assuming WMAP1 and WMAP3 and found that the difference is negligible, so our results are robust to variations in assumed M vir -c vir relation). Here, we briefly recap the process. We use the measured velocity dispersion and its uncertainties, assuming it is a Gaussian distribution, to get a probability distribution for the group virial mass, M vir . Then, we can obtain the concentration, c vir , from the M vir -c vir relationship assuming a reasonable scatter of 0.14 in log c vir (Bullock et al. 2001;Wong et al. 2011). With the critical density and the characteristic overdensity at the lens redshift (Eke, Navarro & Frenk 1998;Eke, Navarro & Steinmetz 2001), we can obtain r vir and a prior probability on the scale radius, r s , via c vir = r vir /r s . The prior probability distribution of the normalization can be calculated by combining r s , the central density of the halo (ρ 0 ), and the critical surface density for lensing ( cr ). (ii) Group (SIS): We convert the velocity dispersion to an Einstein radius via

to get a prior on θ E of 1.4 ± 0.2 arcsec for G1 and 0.4 ± 0.4 arcsec for G2.

For these two models, we also put a prior on the position of the group (see Fig. 2) based on Wilson et al. (2016).

As some of the galaxies inside the group are close to the main lens, these perturbers could individually affect the main lens beyond the second-order distortion terms. We calculate 3 x of the nearby galaxies using the notation and definition in McCully et al. (2017). As 3 x is expressed in terms of the Einstein radius of these perturbers, we convert the measured velocity dispersions of G1 and G2 (250 ± 20 and 130 ± 60 km s -1 , respectively, from Tonry 1998) into corresponding Einstein radii using equation ( 25). For the other galaxies lacking a measurement of the velocity dispersion, we assume that they are located at the group redshift (this assumption maximizes the value of 3 x), and use their relative luminosities compared to either G1 or G2 (depending on the morphology, since G1 is a spiral), to infer a velocity dispersion from the Faber-Jackson relation (Faber & Jackson 1976). We find log 3 x(G1) = -3.68 +0.13 -0.14 (in units of log(arcsec)) and log 3 x(G2) = -4.01 +0.75 -1.07 , whereas the remaining galaxies have log 3 x < -4, and we therefore neglect them. Thus, to test for systematic effects, we model the group either as a single group profile or as a group halo plus either G1, or both G1 and G2, where the galaxies are modelled as SIS mass distributions.

The technique of inferring P (κ ext |d ENV ), based on the Millennium Simulation and observed weighted galaxy number counts, was briefly described in Section 4.3.1. As implemented by Rusu et al. (2017), it requires wide-field, broad-band images to compute photometric redshifts and other physical properties of the galaxies surrounding the lens. The deepest multiband images currently available are provided by the Sloan Digital Sky Survey (SDSS; Adelman-McCarthy et al. 2008) and the Panoramic Survey Telescope and Rapid Response System (Chambers et al. 2016), but these are still relatively shallow, which may result in a biased κ ext (Collett et al. 2013). For our κ ext analysis for PG 1115+080, we therefore use the deep co-added data set from the MPIA 2.2 m telescope described in Section 3.4.

The co-added image has a limiting magnitude of 25.36 ± 0.08, deeper than the control survey, CFHTLenS (r = 24.88 ± 0.16). We perform source detection in this image using SEXTRACTOR (Bertin & Arnouts 1996). For a fair comparison with the control survey, we need to convert our R c magnitudes to r-band magnitudes. However, as we only have a single band, we cannot compute colour terms. Fortunately, a cross-match of the detections in our field with those in SDSS 8 shows that, after correcting for the zero-point offset, the scatter is small, with an rms of ∼0.10 mag, which we add to the photometric error budget. However, we choose a brighter magnitude 8 We ignore the negligible differences between the SDSS and CFHTLenS r-band filters. limit r ≤ 23 mag, in order to be able to use the purely morphological galaxy-star classification of CFHTLenS (Hildebrandt et al. 2012), where objects down to i < 23 mag are classified based solely on the FLUX RADIUS parameter measured by SEXTRACTOR, and because, roughly, ri ∼ 0.5 for our cross-matches with SDSS. Due to the good seeing of our data, we find a clear stellar locus that allows us to determine a good classification threshold for FLUX RADIUS, using the methodology in Coupon et al. (2009). We show the 240 arcsec × 240 arcsec cut-out of the field of view (FOV) in Fig. 3, marking the sources detected down to our magnitude limit.

We compute relative weighted galaxy number counts in terms of simple counts (a weight of unity) as well as using as weight the inverse of the distance between each galaxy in the field and the lens (weighting by 1/r). We do this inside both the 45 and 120 arcsec radius apertures, using the technique from Rusu et al. (2017) and the galaxy catalogue produced earlier, with the exception that we use r-band magnitudes for both the lens field and the CFHTLenS fields, whereas Rusu et al. (2017) used i-band magnitudes. When doing this, we ignore the galaxies confirmed as part of the galaxy group, as the group is explicitly incorporated in our lensing models, and we need to compute κ ext without its contribution. In addition, we account for the galaxies that are expected to be part of the group, but are missed due to the spectroscopic incompleteness, as described in Appendix D. We report our results in Table C1. These numbers are mostly consistent with the unit value, indicating that, after removing the contribution of the galaxy group, the field around the lens is of average density. Finally, we compute P (κ ext |d ENV , γ ) following the technique presented in Birrer et al. (2019), which combines the constraints from both apertures. The combination of apertures results in a tighter distribution, as shown by Rusu et al. (2019a). We show the resulting distributions, corresponding to the various tests of systematics from Section 5.1.6, in Fig. 4 and the summary table in Appendix C.

We summarize the choices that we explore for the mass modelling, including the nearby group/galaxies. For each of the models, we set the weights for the regions containing the AGN images to zero and fix the mass centroid for G1 and G2 at the centre of its light distribution. When modelling the galaxy group as an NFW profile, we use the M vir -c vir relation from Macciò et al. (2008), based on a WMAP5 cosmology (we found that the impact of using different cosmology is negligible). In sum, we explore 160 modelling choices in total, with all different combination of choices among two kinds of main lens models, various mass models for the group, five different resolutions of the reconstructed source, and three different priors of the accretion disc sizes (or we ignore the MTDE). We show the AO imaging reconstruction in Fig. 5 and HST imaging reconstruction in Fig. 6. In

For the HE 0435-1223 system, we model only the AO data, and then combine the results with the HST-based modelling of Wong et al. (2017). This lens system presents a bit of complexity because there are significant contributions to the lensing signal from galaxies at multiple redshifts. In particular, there are five important perturbers (G1-G5) that are close in projection to HE 0435-1223 (see fig. 3 in Wong et al. 2017). Based on the 3 x criterion of McCully et al. (2014,2017), we should include the most massive nearby perturber, G1, explicitly in the model. However, Sluse et al. (2017) show that although the 3 x values of the other four galaxies are not above the threshold when considered individually, when considered together they do show a significant effect. Since G1-G5 are located at different redshifts, we follow Wong et al. (2017) and model this system through the multiplane lens equation (e.g. Blandford & Narayan 1986;Schneider, Ehlers & Falco 1992;Collett & Auger 2014;McCully et al. 2014;Wong et al. 2017). In this case, there is no single time-delay distance, and therefore a particular cosmological model needs to be applied to the analysis. However, if the lens system is dominated by a single primary lens, as is the case for HE 0435-1223, then we can define an effective timedelay distance, D eff t (z , z s ), which is fairly robust to changes in the assumed cosmology.

We follow the same criteria described in Section 5.1.1 and perform 13 iterative correction steps to obtain the final AO PSF of HE 0435-1223. The FWHM of the reconstructed HE 0435-1223 AO PSF is 0.07 arcsec (see Fig. A1).

As we did for PG 1115+080, we model the main lens with either an SPEMD or a composite model. For the composite model, we follow Wong et al. (2017) and set the Gaussian prior for the scale radius to 4.3 ± 2.0 arcsec based on scaling relations derived from the SLACS sample (Gavazzi et al. 2007). The most massive perturber, G1, is modelled as an SIS profile. When modelling G1-G5 simultaneously as SIS distributions, we fix the ratios of their Einstein radii by estimating their stellar masses (Rusu et al. 2017) and then using Bernardi et al. (2011) to convert these to velocity dispersions and then to Einstein radii. We follow Wong et al. (2017) and fix the ratio of Einstein radii, but the global scaling is allowed to vary.

For this system, we have gathered wide-field imaging in a variety of filters, as well as conducting targeted spectroscopy (Rusu et al. 2017;Sluse et al. 2017). Our results on P (κ ext |d ENV , γ ) are presented in Rusu et al. (2017). In this work, we use the shear values determined from our lens modelling of the AO data to update the weighted number counts for the system. Otherwise, we follow the analysis of Rusu et al. (2017) in order to conduct a direct comparison of H 0 from the HST and AO data sets.

We list the systematic tests we have done here. For each of the models, we set the weights for the regions containing the AGN images to zero and fix the mass centroid for G1 at the centre of its light distribution.

(i) A power-law model plus G1 as an SIS. (ii) A composite model plus G1 as an SIS. (iii) A composite model plus the five perturbers (G1-G5). (iv) For all of these models, we test five different source resolutions. See details in Appendix G.

To assess the MST, for each model we perform the importance sampling given the measured velocity dispersion, σ = 222 ± 15 km s -1 inside a 0.54 arcsec × 0.7 arcsec aperture with a seeing of 0.8 arcsec (Wong et al. 2017). We show the AO imaging reconstruction in Fig. 11. The external convergence in Fig. 12. The comparisons between AO results and HST results in two mass models are shown in Fig. 13. and Fig. 14. We show the posteriors of D t in different model choices in Fig. 15 and the posteriors of joint D t in Fig. 16.

Note that for the baryonic component in the composite model, the mass distribution is based on the light distribution in the HST imaging. This is so because an insufficient knowledge about the structure in the wings of the AO PSF introduces a degeneracy between the reconstructed PSF structure and lens galaxy light (see the discussion in Appendix B).

A detailed discussion of the RXJ 1131-1231 lens modelling of the AO imaging can be found in the paper of Chen et al. (2016). To summarize, we used the power-law mass distribution to model the lens potential and used two concentric Sérsic profiles to model the lens light. The satellite galaxy of the main deflector was modelled as an SIS profile. We modelled only the lensing galaxy plus satellite, and did not consider κ ext . The modelling marginalized over five different source resolutions in order to better control the systematics. In this paper, we further explore a different mass model and turn the previous work and the new results in this paper into cosmology.

To add to the previous power-law model, we test a composite model with different source resolutions in this paper. We follow Suyu et al. (2014) to set a Gaussian prior on the NFW scale radius of 18.6 ± 2.6 arcsec, based on the weak lensing analysis of SLACS lenses (Gavazzi et al. 2007) that have similar velocity dispersions to RXJ 1131-1231. For the other parameters, we set uniform priors. We model the satellite light distribution with a circular Sérsic profile, and the satellite mass as an SIS distribution whose centroid is linked to the light centroid.

Note that due to the degeneracy between the reconstructed PSF structure and lens galaxy light, the baryonic mass distribution in the composite model is also based on the light distribution in the HST imaging.

As in Suyu et al. (2013), we use a combination of observations and simulations to estimate the contribution of the LOS mass distribution for RXJ 1131-1231, i.e. P(κ ext ), P(κ ext |γ ), and P (κ ext |d ENV , γ ). Here, γ is the external shear required by the mass models of the main lensing galaxy, while d ENV is the relative overdensity of galaxies within a 45 arcsec aperture that is centred on the lens. This overdensity, ζ 45 arcsec 1 = 1.4 ± 0.05 (following the notation in Birrer et al. 2019), is calculated from galaxies with apparent HST/ACS F814W magnitudes 18.5 ≤ m ≤ 24.5 in both the lens and control samples (Fassnacht et al. 2011). The overdensity and shear values are combined with the simulated lensing data based on the Millennium Simulation (Springel et al. 2005;Hilbert et al. 2009) together with the semi-analytic galaxy model of Henriques et al. (2015), to get the probability distributions for κ ext . We show the results in Fig. 17.

We list the systematic tests we have done including those done in our previous work. In all of the models, the regions near the AGN images are given zero weight.

(i) SPEMD + 2Sérsic lens model. We re-run the model since we did not link the satellite mass position to its light position in Chen et al. (2016).

(ii) A composite model. (iii) For these models, we test five different source resolutions. See details in Appendix H. We use the observed velocity dispersion, 323 ± 20 km s -1 (Suyu et al. 2013), given the κ ext in Section 5.3.2 to sample D t and D without assuming cosmology. We plot the posteriors of D t and D in Fig. 18 and the joint results in Fig. 19.

We present the cosmological inferences based on the distance measurements (see Table 4) of the three gravitational lenses that have AO imaging data, HST imaging data, velocity dispersion measurements, LOS studies, and time-delay measurements. In particular, we present the cosmological inferences based on only the AO imaging data in Section 6.1, while we present the cosmological inferences based on a combination of both the AO and HST imaging in Section 6.2. -4.5 km s -1 Mpc -1 , and RXJ 1131-1231 AO imaging yields H 0 = 77.0 +4.0 -4.6 km s -1 Mpc -1 . The joint analysis of three AO lenses yields H 0 = 75.6 +3.2 -3.3 km s -1 Mpc -1 .

As the AO imaging of HE 0435-1223 and RXJ 1131-1231 is modelled separately from the HST imaging, we have developed a Bayesian approach to properly combine the HST and AO results for these two lenses (see details in Appendix E). In short, since the HST and AO images are independent data sets, we can get the joint probability distribution by multiplying their probability distributions, as long as the prior on the joint parameters is used only once. We can express the joint posterior as 21 presents the marginalized posterior PDF for H 0 assuming flat CDM model. We found that the joint AO + HST of HE 0435-1223 implies a value of the Hubble constant of H 0 = 71.6 +4.7 -4.6 km s -1 Mpc -1 , the joint AO + HST of PG 1115+080 implies H 0 = 81.1 +7.9 -7.1 km s -1 Mpc -1 , and the joint AO + HST of RXJ 1131-1231 implies H 0 = 78.3 +3.4 -3.3 km s -1 Mpc -1 . The combination of three AO + HST lenses yields H 0 = 76.8 +2.6 -2.6 km s -1 Mpc -1 . We found that after combining AO and HST imaging, the dominant sources of uncertainty of both PG 1115+080 and HE 0435-1223 are time-delay measurements, while the dominant source of uncertainty of RXJ 1131-1231 is the LOS mass distribution.

We did the blind analysis on both PG 1115+080 and HE 0435-1223 as well as an extension of our previous analysis of RXJ 1131-1231.

For each system, we combined the AO imaging, HST imaging, the measurements of the lens galaxy's velocity dispersion, the LOS studies from deep wide-area spectroscopic as well as photometric data, and the time-delay measurements from state-of-the-art lightcurve fitting algorithm to infer the value of H 0 . We find that the high S/N AO and HST imaging data yield consistent results, providing an important validation of the AO PSF reconstruction techniques for high-precision lensing work.

This paper demonstrates the ability of using AO imaging to constrain the mass model as well as the value of H 0 . Furthermore, we show that combining AO imaging with HST imaging can further tighten the uncertainties from the lens mass model and thus improve the precision of the determination of H 0 .

In this paper, we infer the value of H 0 under the assumption of a flat CDM model that has a uniform prior on H 0 in the range [0, 150] km s -1 Mpc -1 and a uniform prior on m in the range of [0.05, 0.5]. After unblinding, PG 1115+080 AO imaging yields H 0 = 82.8 +9.4 -8.3 km s -1 Mpc -1 , HE 0435-1223 AO imaging yields H 0 = 70.1 +5.3 -4.5 km s -1 Mpc -1 , and RXJ 1131-1231 AO imaging yields H 0 = 77.0 +4.0 -4.6 km s -1 Mpc -1 . The joint analysis of three AO lenses yields H 0 = 75.6 +3.2 -3.3 km s -1 Mpc -1 . The joint AO + HST of PG 1115+080 yields H 0 = 81.1 +7.9 -7.1 km s -1 Mpc -1 , the joint AO + HST of HE 0435-1223 yields H 0 = 71.6 +4.7 -4.6 km s -1 Mpc -1 , and the joint AO + HST of RXJ 1131-1231 yields H 0 = 78.3 +3.4 -3.3 km s -1 Mpc -1 . The combination of three AO + HST lenses yields H 0 = 76.8 +2.6 -2.6 km s -1 Mpc -1 . We refer the reader to the paper by Wong et al., where the results presented here will be combined with a self-consistent analysis of three previously published systems to carry out a full cosmological investigation.

Here, we show the κ ext distribution of HE 0435-1223 in Fig. 12 and present the summary table of the constraints used to estimate κ ext for all three systems in Table C1. The spectroscopic coverage of the FOV around PG 1115+080 is incomplete. Down to R c ≤ 22.5 and within 120 arcsec radius around the lens, there are 63 galaxies, out of which 33 have spectroscopy (Wilson et al. 2016), 11 of which are part of the galaxy group at z = 0.31, including the lensing galaxy. This means that there may be other galaxies within this magnitude range and radius from the lens that are also part of the galaxy group associated with the lensing galaxy, but that are missed due to spectroscopic incompleteness. 9 In Section 5.1.3, we have specifically computed the lensing properties of this group, based on its physical properties derived by Momcheva et al. (2015). As a result, when we compute κ ext at the location of the lens, using the weighted number counts approach, we must remove the galaxies that are part of this group, as the convergence from the group has already been included in the lensing models, and must not be double counted. While the galaxies that are known to be part of the group can easily be removed, we must also account 9 For RXJ 1131-1231 and HE 0435-1223, we have not incorporated any galaxy groups in the lens models; therefore, we do not need to remove the contribution of such groups from our estimate of the external convergence based on number counts. Hence, the spectroscopic incompleteness is not relevant to our analysis.

for the galaxies expected to be missed due to our spectroscopic incompleteness.

Following the technique presented in Rusu et al. (2019a), we use two different approaches to estimate the number of missing galaxies that are expected to be part of the group, but that are not identified as group members because of our spectroscopic incompleteness. In the first approach, we use the knowledge provided by the number of known group members, the number of galaxies with spectroscopy, and the total number of detected galaxies (within the given magnitude and aperture radius), and we apply Poisson statistics to estimate a number of 10 ± 5 missing galaxies (median, 16th, and 84th percentiles). In the second approach, we use the group velocity dispersion and virial radius from Wilson et al. (2016), and we estimate the expected number of galaxies inside the virial radius using the empirical relation from Andreon & Hurn (2010). Using the measured offset from the group centroid to the lens, and propagating all uncertainties, we measure the expected number of missing galaxies at the intersection of the sphere of virial radius and the 120 arcsec radius cylinder centred on the lens to be 1 +3 -1 . We plot the distributions of these numbers in Fig. D1. The first approach predicts a significantly larger number of missing galaxies than the second. In fact, due to the small value of the velocity dispersion, the second method would only predict a total number of 8 +6 -4 galaxies, therefore less than the number of confirmed group members, unless we enforce this constraint. This discrepancy may be due to the shallow absolute magnitude limit of M V = -20 used by Andreon & Hurn (2010), corresponding to r ∼ 21 at z ∼ 0.3, therefore significantly brighter than our limiting magnitude. We note, however, that the two techniques produced results that were in agreement for a different lens, described in Rusu et al. (2019c).

In view of the above, and also due to the fact that it avoids any physical assumption, we consider the first method to be more reliable. Finally, when computing weighted galaxy counts, we do this by randomly sampling 10 times from the distribution of missing galaxy numbers, and then randomly excluding that number Figure D1. Estimated number of missing galaxy group members inside the ≤120 arcsec radius from the lens system for PG 1115+080, computed with two methods, with or without imposing the prior that the group consists of at least the number of galaxies spectroscopically confirmed to be members. of galaxies from our catalogue of galaxies inside the 120 arcsec apertures.

We present the BIC of the power-law models in Table G1 and composite models in Table G2.