A SHARP view of H0LiCOW: $H_{0}$ from three time-delay gravitational lens systems with adaptive optics imaging

TL;DR

This paper measures the Hubble Constant using three gravitational lens systems with adaptive optics and HST imaging, demonstrating consistent results and the viability of AO data for lens-based cosmology.

Contribution

It introduces a method combining AO and HST imaging for precise lens modeling to measure H0, expanding the use of AO in time-delay cosmography.

Findings

Joint AO+HST analysis yields H0=76.8 km/s/Mpc with 3% uncertainty.

AO-only analysis produces consistent H0 estimates across three lenses.

AO data effectively complement HST in high-resolution lens modeling.

Abstract

We present the measurement of the Hubble Constant, $H_0$, with three strong gravitational lens systems. We describe a blind analysis of both PG1115+080 and HE0435-1223 as well as an extension of our previous analysis of RXJ1131-1231. For each lens, we combine new adaptive optics (AO) imaging from the Keck Telescope, obtained as part of the SHARP AO effort, with Hubble Space Telescope (HST) imaging, velocity dispersion measurements, and a description of the line-of-sight mass distribution to build an accurate and precise lens mass model. This mass model is then combined with the COSMOGRAIL measured time delays in these systems to determine $H_{0}$. We do both an AO-only and an AO+HST analysis of the systems and find that AO and HST results are consistent. After unblinding, the AO-only analysis gives $H_{0}=82.8^{+9.4}_{-8.3}~\rm km\,s^{-1}\,Mpc^{-1}$ for PG1115+080, $H_{0}=70.1^{+5.3}_{-4.5}~\rm km\,s^{-1}\,Mpc^{-1}$ for HE0435-1223, and $H_{0}=77.0^{+4.0}_{-4.6}~\rm km\,s^{-1}\,Mpc^{-1}$ for RXJ1131-1231. The joint AO-only result for the three lenses is $H_{0}=75.6^{+3.2}_{-3.3}~\rm km\,s^{-1}\,Mpc^{-1}$. The joint result of the AO+HST analysis for the three lenses is $H_{0}=76.8^{+2.6}_{-2.6}~\rm km\,s^{-1}\,Mpc^{-1}$. All of the above results assume a flat $\Lambda$ cold dark matter cosmology with a uniform prior on $\Omega_{\textrm{m}}$ in [0.05, 0.5] and $H_{0}$ in [0, 150] $\rm km\,s^{-1}\,Mpc^{-1}$. This work is a collaboration of the SHARP and H0LiCOW teams, and shows that AO data can be used as the high-resolution imaging component in lens-based measurements of $H_0$. The full time-delay cosmography results from a total of six strongly lensed systems are presented in a companion paper.

Tables

Table 1. Table 1: Details of the AO observations.

Lens	Instrument	Date	$t_{\exp }$
HE0435$-$1223	Keck/NIRC2-N	2010-12-02	8100
…	Keck/NIRC2-N	2012-01-01	2400
PG1115+080	Keck/NIRC2-N	2012-05-15	900
…	Keck/NIRC2-N	2017-04-11	900
RXJ1131$-$1231	Keck/NIRC2-W	2012-05-16	1800
…	Keck/NIRC2-W	2012-05-18	1800

Table 2. Table 2: Lensing parameters for creating the magnification maps of microlensing time-delay effect. The values of κ 𝜅 \kappa , γ 𝛾 \gamma , and κ ∗ / κ subscript 𝜅 𝜅 \kappa_{*}/\kappa of HE 0435 − - 1223 andPG 1115 + + 080 are from our lens models

Name	img	$\kappa$	$\gamma$	${\kappa }_{\ast }/\kappa$
HE 0435$-$1223	A	0.473	0.358	0.347
(This work)	B	0.630	0.540	0.361
	C	0.494	0.327	0.334
	D	0.686	0.575	0.380
PG 1115$+$080	A1	0.424	0.491	0.259
(This work)	A2	0.451	0.626	0.263
	B	0.502	0.811	0.331
	C	0.356	0.315	0.203
RXJ 1131$-$1231	A	0.526	0.410	0.429
(This work)	B	0.459	0.412	0.434
	C	0.487	0.306	0.414
	D	0.894	0.807	0.581

Table 3. Table 3: Quasar source parameters: the black hole mass ( M bh subscript 𝑀 bh M_{\rm bh} ), the Eddington ratio ( L / L E 𝐿 subscript 𝐿 E L/L_{\rm E} ), the accretion efficiency ( η 𝜂 \eta ), the inclination angle (incl.), and the references.

Name	$M_{\mathrm{bh}}$ [$M_{☉}$]	$L/L_E$	$\eta$	incl.[deg]	reference
HE 0435$-$1223	$5.75\times {10}^8$	0.1	0.1	0	Sluse et al. (2012)
PG 1115$+$080	$1.2\times {10}^9$	0.1	0.1	0	Morgan et al. (2010)
RXJ 1131$-$1231	$1.3\times {10}^8$	0.1	0.1	0	Dai et al. (2010)

Table 4. Table 4: Joined constraints used to estimate κ ext subscript 𝜅 ext \kappa_{\textrm{ext}}

Lens	model name	${\gamma }_{\mathrm{ext}}$	${\zeta }_1^{45″}$	${\zeta }_{1/r}^{45″}$	${\zeta }_1^{120″}$	${\zeta }_{1/r}^{120″}$
RXJ 1131$-$1231	power law G	$0.083\pm 0.003$	${1.4}_{-0.05}^{+0.05}$
RXJ 1131$-$1231	composite G	$0.058\pm 0.005$	${1.4}_{-0.05}^{+0.05}$
HE 0435$-$1223	power law G + G1	$0.041\pm 0.018$
HE 0435$-$1223	composite G + G1	$0.026\pm 0.026$
HE 0435$-$1223	composite G + 5 perturbers	$0.056\pm 0.026$
HE 0435$-$1223	power law G + G1 + LOS	$0.041\pm 0.018$	${1.27}_{-0.05}^{+0.05}$	${1.31}_{-0.05}^{+0.05}$
HE 0435$-$1223	composite G + G1 + LOS	$0.026\pm 0.026$	${1.27}_{-0.05}^{+0.05}$	${1.31}_{-0.05}^{+0.05}$
HE 0435$-$1223	composite G + 5 perturbers + LOS	$0.056\pm 0.026$	${1.21}_{-0.05}^{+0.05}$	${1.17}_{-0.05}^{+0.05}$
PG 1115$+$080	power law G + NFW group	$0.027\pm 0.009$	${1.11}_{-0.11}^{+0.01}$	${1.27}_{-0.23}^{+0.01}$	${1.10}_{-0.09}^{+0.05}$	${0.98}_{-0.10}^{+0.03}$
PG 1115$+$080	power law G + SIS group	$0.061\pm 0.009$	${1.11}_{-0.11}^{+0.01}$	${1.27}_{-0.23}^{+0.01}$	${1.10}_{-0.09}^{+0.05}$	${0.98}_{-0.10}^{+0.03}$
PG 1115$+$080	power law G + NFW group + G1	$0.054\pm 0.009$	${1.00}_{-0.11}^{+0.01}$	${1.02}_{-0.12}^{+0.03}$	${1.09}_{-0.09}^{+0.04}$	${0.96}_{-0.10}^{+0.04}$
PG 1115$+$080	power law G + NFW group + G1 + G2	$0.058\pm 0.012$	${0.89}_{-0.11}^{+0.01}$	${0.84}_{-0.24}^{+0.06}$	${1.07}_{-0.09}^{+0.05}$	${0.94}_{-0.10}^{+0.04}$
PG 1115$+$080	composite G + NFW group	$0.048\pm 0.009$	${1.11}_{-0.11}^{+0.01}$	${1.27}_{-0.23}^{+0.01}$	${1.10}_{-0.09}^{+0.05}$	${0.98}_{-0.10}^{+0.03}$
PG 1115$+$080	composite G + SIS group + G1	$0.072\pm 0.008$	${1.00}_{-0.11}^{+0.01}$	${1.02}_{-0.12}^{+0.03}$	${1.09}_{-0.09}^{+0.04}$	${0.96}_{-0.10}^{+0.04}$
PG 1115$+$080	composite G + NFW group + G1	$0.072\pm 0.008$	${1.00}_{-0.11}^{+0.01}$	${1.02}_{-0.12}^{+0.03}$	${1.09}_{-0.09}^{+0.04}$	${0.96}_{-0.10}^{+0.04}$
PG 1115$+$080	composite G + NFW group + G1 + G2	$0.060\pm 0.008$	${0.89}_{-0.11}^{+0.01}$	${0.84}_{-0.24}^{+0.06}$	${1.07}_{-0.09}^{+0.05}$	${0.94}_{-0.10}^{+0.04}$

Table 5. Table 5: Total 20 power-law models of PG1115+080 ordered in increased Δ Δ \Delta BIC value.

Main lens model	perturbers	$S_r$	$\Delta$ BIC	posterior weight
SPEMD	group (NFW)	$37\times 37$	0	1
SPEMD	group (NFW)+G1+G2	$37\times 37$	11	0.9820
SPEMD	group (NFW)+G1+G2	$39\times 39$	18	0.9469
SPEMD	group (NFW)+G1	$37\times 37$	31	0.8967
SPEMD	group (NFW)	$38\times 38$	35	0.8338
SPEMD	group (NFW)+G1+G2	$33\times 33$	41	0.7614
SPEMD	group (NFW)+G1	$41\times 41$	46	0.6827
SPEMD	group (NFW)+G1+G2	$35\times 35$	46	0.6827
SPEMD	group (NFW)+G1	$39\times 39$	54	0.5198
SPEMD	group (SIS)   +G1	$36\times 36$	72	0.4415
SPEMD	group (SIS)   +G1	$34\times 34$	74	0.3684
SPEMD	group (NFW)+G1	$35\times 35$	88	0.3019
SPEMD	group (NFW)+G1+G2	$41\times 41$	90	0.2433
SPEMD	group (SIS)   +G1	$32\times 32$	91	0.1930
SPEMD	group (NFW)	$32\times 32$	92	0.1510
SPEMD	group (NFW)	$36\times 36$	95	0.1172
SPEMD	group (NFW)	$34\times 34$	114	0.0911
SPEMD	group (NFW)+G1	$43\times 43$	126	0.0721
SPEMD	group (SIS)   +G1	$40\times 40$	127	0.0597
SPEMD	group (SIS)   +G1	$38\times 38$	148	0.0536

Table 6. Table 6: Total 20 composite models of PG1115+080 ordered in increased Δ Δ \Delta BIC value.

Main lens model	perturbers	$S_r$	$\Delta$ BIC	posterior weight
COMPOSITE	group (SIS)   +G1	$37\times 37$	0	1
COMPOSITE	group (NFW)	$30\times 30$	34	0.9799
COMPOSITE	group (NFW)	$32\times 32$	55	0.9409
COMPOSITE	group (NFW)	$36\times 36$	67	0.8853
COMPOSITE	group (NFW)+G1+G2	$39\times 39$	69	0.8162
COMPOSITE	group (NFW)+G1+G2	$37\times 37$	80	0.7374
COMPOSITE	group (NFW)+G1	$33\times 33$	88	0.6529
COMPOSITE	group (NFW)+G1	$35\times 35$	90	0.5664
COMPOSITE	group (NFW)+G1	$32\times 32$	96	0.4815
COMPOSITE	group (NFW)+G1+G2	$33\times 33$	118	0.4010
COMPOSITE	group (NFW)+G1+G2	$35\times 35$	122	0.3273
COMPOSITE	group (NFW)	$34\times 34$	132	0.2620
COMPOSITE	group (SIS)   +G1	$31\times 31$	138	0.2057
COMPOSITE	group (SIS)   +G1	$35\times 35$	145	0.1584
COMPOSITE	group (NFW)+G1	$39\times 39$	146	0.1200
COMPOSITE	group (SIS)   +G1	$33\times 33$	149	0.0897
COMPOSITE	group (SIS)   +G1	$39\times 39$	156	0.0669
COMPOSITE	group (NFW)+G1+G2	$31\times 31$	163	0.0506
COMPOSITE	group (NFW)+G1	$37\times 37$	193	0.0402
COMPOSITE	group (NFW)	$38\times 38$	232	0.0352

Table 7. Table 7: Total 5 powerlaw models of HE 0435 − - 1223 ordered in increased Δ Δ \Delta BIC value.

Main lens model	perturbers	$S_r$	$\Delta$ BIC	posterior weight
SPEMD	G1	$50\times 50$	0	1
SPEMD	G1	$40\times 40$	5	0.9743
SPEMD	G1	$45\times 45$	7	0.9328
SPEMD	G1	$35\times 35$	14	0.8914
SPEMD	G1	$30\times 30$	16	0.8658

Table 8. Table 8: Total 10 composite models of HE 0435 − - 1223 ordered in increased Δ Δ \Delta BIC value.

Main lens model	perturbers	$S_r$	$\Delta$ BIC	posterior weight
COMPOSITE	G1–G5	$43\times 43$	0	1
COMPOSITE	G1	$43\times 43$	10	0.9855
COMPOSITE	G1	$41\times 41$	11	0.9580
COMPOSITE	G1–G5	$41\times 41$	12	0.9201
COMPOSITE	G1	$39\times 39$	16	0.8756
COMPOSITE	G1–G5	$39\times 39$	19	0.8288
COMPOSITE	G1–G5	$37\times 37$	20	0.7844
COMPOSITE	G1	$35\times 35$	21	0.7466
COMPOSITE	G1	$37\times 37$	21	0.7466
COMPOSITE	G1–G5	$35\times 35$	23	0.7046

Table 9. Table 9: Total 5 power-law model of RXJ 1131 − - 1231 ordered in increased Δ Δ \Delta BIC value.

Main lens model	perturbers	$S_r$	$\Delta$ BIC	posterior weight
SPEMD	satellite (SIS)	$79\times 79$	0	1
SPEMD	satellite (SIS)	$77\times 77$	85	0.4618
SPEMD	satellite (SIS)	$75\times 75$	176	0.0912
SPEMD	satellite (SIS)	$73\times 73$	499	0.0071
SPEMD	satellite (SIS)	$71\times 71$	795	0.0002

Table 10. Table 10: Total 5 composite models of RXJ 1131 − - 1231 ordered in increased Δ Δ \Delta BIC value.

Main lens model	perturbers	$S_r$	$\Delta$ BIC	posterior weight
COMPOSITE	satellite (SIS)	$71\times 71$	0	1
COMPOSITE	satellite (SIS)	$70\times 70$	398	0.4050
COMPOSITE	satellite (SIS)	$73\times 73$	915	0.0575
COMPOSITE	satellite (SIS)	$74\times 74$	1641	0.0025
COMPOSITE	satellite (SIS)	$72\times 72$	1798	0.0000

Equations

$$D_{\Delta t}\equiv \left(1+z_{\ell }\right)\frac{D_{\ell }{D}_s}{D_{\ell s}}\propto H_0^{-1},$$

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

$$\begin{array}{rl}& P({\mathbf{\eta }}_i|{\mathbf{d}}_i,{\mathbf{\Delta }\mathbf{t}}_i,{\sigma }_i,{\mathbf{d}}_{{\text{ENV}}_i},\mathbf{A})\\ & \propto P({\mathbf{d}}_i,{\mathbf{\Delta }\mathbf{t}}_i,{\sigma }_i,{\mathbf{d}}_{{\text{ENV}}_i}|{\mathbf{\eta }}_i,\mathbf{A})P({\mathbf{\eta }}_i|\mathbf{A}),\end{array}$$

$$\begin{array}{rl}& P({\mathbf{d}}_i,{\mathbf{\Delta }\mathbf{t}}_i,{\sigma }_i,{\mathbf{d}}_{{\text{ENV}}_i}|{\mathbf{\eta }}_i,{\mathbf{A}}_i)\\ & =P({\mathbf{d}}_i|{\mathbf{\eta }}_i,{\mathbf{A}}_i)P({\mathbf{\Delta }\mathbf{t}}_i|{\mathbf{\eta }}_i,{\mathbf{A}}_i)P({\sigma }_i|{\mathbf{\eta }}_i,\mathbf{A})P({\mathbf{d}}_{{\text{ENV}}_i}|{\mathbf{\eta }}_i,{\mathbf{A}}_i).\end{array}$$

$$\begin{array}{rl}P({\mathbf{\eta }}_i|{\mathbf{d}}_{i,\text{tot}})& =\int P({\mathbf{\eta }}_i|{\mathbf{d}}_{i,\text{tot}},{\mathbf{A}}_i)P({\mathbf{A}}_i)d{\mathbf{A}}_i\\ & \approx \sum _{k}P({\mathbf{\eta }}_i|{\mathbf{d}}_{i,\text{tot}},{\mathbf{A}}_{i,k})P({\mathbf{A}}_{i,k}),\end{array}$$

$$\text{BIC}=\text{ln}(n)k-2\text{ln}(\hat{L}),$$

$${\kappa }_{\lambda }(\theta )=(1-\lambda )\kappa (\theta )+\lambda ,$$

$$D_{\Delta t}=D_{\Delta t,\lambda }^{\mathrm{true}}=\frac{D_{\Delta t}^{\mathrm{model}}}{1-\lambda }.$$

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

$${\psi }_{\text{ext}}(\theta ,\phi )=\frac{1}{2}{\gamma }_{\text{ext}}{\theta }^2\cos 2(\phi -{\phi }_{\text{ext}}),$$

$${\rho }_{\text{local}}(r)={\rho }_0{r}_0^n{F}_n(r),$$

$$\begin{array}{rl}{M}_{\text{local}}& =4\pi {\int }_0^{\infty }dz{\int }_0^{R_{\text{Ein}}}{\rho }_0{r}_0^n{F}_n(\sqrt{s^2+z^2})sds\\ & ={\rho }_0{r}_0^n{M}_{\text{2D}}(R_{\text{Ein}}),\end{array}$$

$$M_{\text{local}}=M_{\text{Ein}}-M_{\text{ext}}=\pi R_{\text{Ein}}^2{\Sigma }_{\text{cr}}(1-{\kappa }_{\text{ext}}),$$

$${\Sigma }_{\text{cr}}=\frac{c^2}{4\pi G}\frac{D_s}{D_{\ell }{D}_{\ell s}}$$

$${\rho }_0{r}_0^n=\frac{\pi R_{\text{Ein}}^2{\Sigma }_{\text{cr}}(1-{\kappa }_{\text{ext}})}{M_{\text{2D}}(R_{\text{Ein}})}.$$

$${\rho }_{\text{local}}=\frac{\pi R_{\text{Ein}}^2{\Sigma }_{\text{cr}}(1-{\kappa }_{\text{ext}})}{M_{\text{2D}}(R_{\text{Ein}})}{F}_n(r).$$

$${\Sigma }_{\text{cr}}=\frac{c^2}{4\pi G}\frac{D_{\Delta t}}{1+z_{\ell }}\frac{1}{D_{\ell }^2}=\frac{c^2}{4\pi G}\frac{D_{\Delta t}^{\text{model}}}{(1+z_{\ell })(1-{\kappa }_{\text{ext}})}\frac{1}{D_{\ell }^2},$$

$${\beta }_{\text{ani}}=\frac{r^2}{r_{\text{ani}}^2+r^2},$$

$${\rho }_{\ast }=\frac{I_{0}a}{2\pi r(r+a{)}^3},$$

$$\frac{1}{{\rho }_{\ast }}\frac{d({\rho }_{\ast }{\sigma }_{\text{r}})}{dr}+2\frac{{\beta }_{\text{ani}}{\sigma }_{\text{r}}}{r}=-\frac{GM(r)}{r^2},$$

$$({\sigma }^{\text{P}}{)}^2=\frac{{\int }_{\mathcal{A}}[I(R){\sigma }_s^2\ast \mathcal{P}]d\mathcal{A}}{{\int }_{\mathcal{A}}[I(R)\ast \mathcal{P}]d\mathcal{A}},$$

$$I(R){\sigma }_s^2=2{\int }_R^{\infty }(1-{\beta }_{\text{ani}}\frac{R^2}{r^2})\frac{{\rho }_{\ast }{\sigma }_{\text{r}}^{2}rdr}{\sqrt{r^2-R^2}}.$$

$$⟨R_{\text{E}}⟩=\sqrt{\frac{D_s{D}_{\ell s}}{D_{\ell }}\frac{4G⟨M_{\star }⟩}{c^2}}$$

$$\Delta t_{ij}=(D_{\Delta t}/c)\Delta {\tau }_{ij}+(t_i-t_j),$$

$${\sigma }^2={\theta }_E\frac{c^2}{4\pi }\frac{D_s}{D_{\ell s}},$$

$$\begin{array}{rl}P(\mathbf{\eta }& |{\mathbf{d}}_{\text{HST}},{\mathbf{d}}_{\text{AO}},\sigma ,{\mathbf{d}}_{\text{ENV}},\text{H},\text{MS})\\ & \propto \frac{P(\mathbf{\eta }|{\mathbf{d}}_{\text{HST}},\sigma ,{\mathbf{d}}_{\text{ENV}},\text{H},\text{MS})P(\mathbf{\eta }|{\mathbf{d}}_{\text{AO}},\text{H},\text{MS})}{P(\mathbf{\eta }|\text{H},\text{MS})},\end{array}$$

$$\begin{array}{rl}P(& {\varphi }_{ij}^{\text{H}},{\kappa }_{\text{ext}},{\gamma }_{\text{ext}},\mathbf{\zeta }|{\mathbf{d}}_{\text{HST}},{\mathbf{d}}_{\text{AO}},\sigma ,{\mathbf{d}}_{\text{ENV}},\text{H},\text{MS})\\ =& P({\mathbf{d}}_{\text{HST}},{\mathbf{d}}_{\text{AO}},\sigma ,{\mathbf{d}}_{\text{ENV}}|{\varphi }_{ij}^{\text{H}},{\kappa }_{\text{ext}},{\gamma }_{\text{ext}},\mathbf{\zeta },\text{H},\text{MS})\\ & \cdot \frac{P({\varphi }_{ij}^{\text{H}},{\kappa }_{\text{ext}},{\gamma }_{\text{ext}},\mathbf{\zeta }|\text{H},\text{MS})}{P({\mathbf{d}}_{\text{HST}},{\mathbf{d}}_{\text{AO}},\sigma ,{\mathbf{d}}_{\text{ENV}}|\text{H},\text{MS})},\end{array}$$

$$\begin{array}{rl}P& ({\varphi }_{ij}^{\text{H}},{\kappa }_{\text{ext}},{\gamma }_{\text{ext}},\mathbf{\zeta }|\text{H},\text{MS})\\ & =P({\varphi }_{ij}^{\text{H}}|\text{H})P({\kappa }_{\text{ext}}|\text{H,MS})P({\gamma }_{\text{ext}}|\text{H})P(\mathbf{\zeta }|\text{H})\end{array}$$

$$\begin{array}{rl}P(& {\mathbf{d}}_{\text{HST}},{\mathbf{d}}_{\text{AO}},\sigma ,{\mathbf{d}}_{\text{ENV}}|{\varphi }_{ij}^{\text{H}},{\kappa }_{\text{ext}},{\gamma }_{\text{ext}},\mathbf{\zeta },\text{H},\text{MS})\\ =& P({\mathbf{d}}_{\text{HST}},\sigma ,{\mathbf{d}}_{\text{ENV}}|{\varphi }_{ij}^{\text{H}},{\kappa }_{\text{ext}},{\gamma }_{\text{ext}},\mathbf{\zeta },\text{H},\text{MS})\\ & \cdot P({\mathbf{d}}_{\text{AO}}|{\varphi }_{ij}^{\text{H}},{\kappa }_{\text{ext}},{\gamma }_{\text{ext}},\mathbf{\zeta },\text{H},\text{MS}).\end{array}$$

$$\begin{array}{rl}P(& {\varphi }_{ij}^{\text{H}},{\kappa }_{\text{ext}},{\gamma }_{\text{ext}},\mathbf{\zeta }|{\mathbf{d}}_{\text{HST}},{\mathbf{d}}_{\text{AO}},\sigma ,{\mathbf{d}}_{\text{ENV}},\text{H},\text{MS})\\ =& P({\mathbf{d}}_{\text{HST}},\sigma ,{\mathbf{d}}_{\text{ENV}}|{\varphi }_{ij}^{\text{H}},{\kappa }_{\text{ext}},{\gamma }_{\text{ext}},\mathbf{\zeta },\text{H},\text{MS})\\ & \cdot \frac{P({\mathbf{d}}_{\text{AO}}|{\varphi }_{ij}^{\text{H}},{\kappa }_{\text{ext}},{\gamma }_{\text{ext}},\mathbf{\zeta },\text{H},\text{MS})P({\varphi }_{ij}^{\text{H}},{\kappa }_{\text{ext}},{\gamma }_{\text{ext}},\mathbf{\zeta }|\text{H},\text{MS})}{P({\mathbf{d}}_{\text{HST}},{\mathbf{d}}_{\text{AO}},\sigma ,{\mathbf{d}}_{\text{ENV}}|\text{H},\text{MS})},\end{array}$$

Full text

††pubyear: 2019††pagerange: A SHARP view of H0LiCOW: $H_{0}$ from three time-delay gravitational lens systems with adaptive optics imaging–10

1 Introduction

1.1 Distance Measurement Discrepancy: A 4.4$\sigma$ Tension on the Value of $H_{0}$

The temperature anisotropies of the Cosmic Microwave Background (CMB) and density correlations of Baryon Acoustic Oscillations (BAO) obtained with the Wilkinson Microwave Anisotropy Probe, Planck satellite, and BAO surveys provide strong support to the standard flat $\Lambda$CDM cosmological model (e.g., Komatsu et al., 2011; Hinshaw et al., 2013; Planck Collaboration et al., 2018). Under a few strong assumptions, such as flatness and constant dark energy density, these data give sub-percent 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 $\Lambda$CDM model (see the illustration in Fig. 4 in Cuesta et al., 2015), leading to a $\sim 4.4\sigma$ tension in $H_{0}$ between the value predicted by the CMB and the local value ($=74.03\pm 1.42~{}\rm 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\pm 1.1~{}\rm 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\pm 1.30~{}\rm km\,s^{-1}\,Mpc^{-1}$ (Macaulay et al., 2019). Both results are in excellent agreement with the Planck value ($H_{0}=67.27\pm 0.6~{}\rm km\,s^{-1}\,Mpc^{-1}$) under the assumption of flat $\Lambda$CDM model (Planck Collaboration et al., 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 also have passed a range of systematic checks (e.g., Efstathiou, 2014; Cardona et al., 2017; Zhang et al., 2017; Follin & Knox, 2018; Feeney et al., 2018; Dhawan et al., 2018; Riess et al., 2019), and rule out the local void scenario (Keenan et al., 2013; Fleury et al., 2017; Shanks et al., 2019; Kenworthy et al., 2019).

There have been several attempts to address this $\sim 4.4\sigma$ tension by extending the standard cosmological model, either by changing the size of the sound horizon in the early Universe (e.g., Heavens et al., 2014; Wyman et al., 2014; Cuesta et al., 2015; Alam et al., 2017; Kreisch et al., 2019; Poulin et al., 2018; Agrawal 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 et al., 2016; Joudaki et al., 2018; Lemos et al., 2019; Aylor et al., 2019) also have 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 et al., 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-$\ell$ CMB power spectrum prefers an even lower $H_{0}$ value than that from the low-$\ell$ CMB power spectrum (Addison et al., 2016; Planck Collaboration et al., 2018).

Given the various tensions across different datasets, any convincing resolution to the $H_{0}$ tensions, either due to unknown systematics or 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.

1.2 Distance Measurement from Time-Delay Cosmography

Time-delay cosmography is not only a completely independent technique of the distance ladder methods, but it also has the advantage of being a one-step measurement of combined cosmological distances. In addition, time-delay cosmography is a complementary and cost-effective 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 time-delay distance (Suyu et al., 2010), which is a ratio of the angular diameter distances in the system:

[TABLE]

where $D_{\ell}$ is the distance to the lensing galaxy, $D_{s}$ is the distance to the background source, and $D_{\ell s}$ is the distance between the lens and the source. Furthermore, we can make a separate determination of $D_{\ell}$ by measuring the velocity dispersion of the lensing galaxy (Jee et al., 2015, 2016; Birrer et al., 2016, 2019). First proposed by Refsdal (1964), the $H_{0}$ measurement requires modeling 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_{\Delta 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 collaboration111$H_{0}$ Lenses in COSMOGRAIL’s Wellspring, www.h0licow.org 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\begin{subarray}{c}+5.4\\ -5.1\end{subarray}~{}\rm 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, 2013, 2014; Sluse et al., 2017; Rusu et al., 2017; Wong et al., 2017; Bonvin et al., 2017) to obtain a joint inference on $H_{0}$ with 3% precision: $H_{0}=72.5\begin{subarray}{c}+2.1\\ -2.3\end{subarray}~{}\rm km\,s^{-1}\,Mpc^{-1}$ (Birrer et al., 2019). This result agrees with the Riess et al. (2019) value within the $1\sigma$ uncertainties. More recently, the collaboration completed its analysis of WFI2033$-$4723 using the time delays from COSMOGRAIL222COSmological MOnitoring of GRAvItational Lenses (Bonvin et al., 2019; Sluse et al., 2019; Rusu et al., 2019a).

Achieving the goal of obtaining a 1% 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 analyzed to date. Many new lensed quasars have been discovered (e.g., Lin et al., 2017; Schechter et al., 2017; Agnello et al., 2018; Ostrovski et al., 2018; Williams et al., 2018; Rusu et al., 2019b; Lemon et al., 2019) using state-of-the-art lens-finding techniques applied to current large sky surveys (e.g. Joseph et al., 2014; Agnello, 2017; Petrillo et al., 2017; Ostrovski 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% $H_{0}$ measurement from time-delay cosmography is a realistic expectation in the near future (e.g., Jee et al., 2015, 2016; de Grijs et al., 2017; Suyu et al., 2018; Shajib et al., 2018, Jee et al. 2019 submitted) if we can control the systematic effects to a sub-percent level (Dobler et al., 2013; Liao et al., 2015; Ding et al., 2018).

1.3 Lens Modeling with Adaptive Optics Data

A critical input to achieving robust modeling 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 HST (e.g., Suyu et al., 2009; Suyu et al., 2010; Birrer et al., 2015; Wong et al., 2017; Birrer et al., 2016, 2019). Adaptive optics (AO) observations with ground-based telescopes can provide imaging with comparable or better angular resolution to the HST data, especially for systems that are faint in the optical and bright at near infrared wavelengths, thus providing an attractive alternative for modeling 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.

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 et al., 2010, 2012; Vegetti et al., 2012; Hsueh et al., 2016, 2017; Hsueh 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., 2007, 2008; Courbin et al., 2011; Tewes et al., 2013a, b; Rathna Kumar et al., 2013; Bonvin et al., 2017, 2018), which aims to provide the highest-precision measurements of time delays. The lens environment of RXJ 1131$-$1231 and HE 0435$-$1223 studies are provided by the H0LiCOW team (Suyu et al., 2014; Sluse et al., 2017; Rusu 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 analyze the data. In Section 5, we elaborate the detailed lens modeling 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.

2 Basic theory

In this section, 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.

2.1 Time-Delay Cosmography

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

[TABLE]

where $\bm{\theta}$, $\bm{\beta}$, and $\psi\left(\bm{\theta}\right)$ are the image location, the source location, and the projected two-dimension lensing gravitational potential, respectively (Shapiro, 1964; Refsdal, 1964). The $\lambda$ parameter represents the lack of perfect knowledge of the full mass distribution, as discussed below.

The advantage of this formulation is the separability of the cosmographic information, contained in the $D_{\Delta t}$ parameter, and the lens modeling. This allows one to infer cosmographic information without the need for cosmological priors in the lens modeling.

However, because of the Mass-Sheet Transformation (MST, Falco et al., 1985; Schneider & Sluse, 2013, 2014), the determination of $D_{\Delta t}$ is subject to an understanding of $\lambda$ (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 et al., 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 $\lambda$. We refer interested readers to Treu & Marshall (2016) and Suyu et al. (2018) for more details.

2.2 Joint Inference

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

[TABLE]

where $P(\bm{d}_{i},\bm{\Delta t}_{i},\sigma_{i},\bm{d}_{\textrm{ENV}_{i}}|\bm{\eta}_{i},\bm{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,

[TABLE]

In order to explore the unmodeled systematic uncertainties that may arise from modeling choices, we vary the content of $\bm{A}$ for each lens. The marginalized integral can be expressed as

[TABLE]

where $\bm{A}_{i,k}$ and $\bm{d}_{i,\textrm{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(\bm{A}_{i,k})$, by using the Bayesian Information Criterion (BIC), which is defined as

[TABLE]

where $n$ is the number of data points including the lens imaging, 8 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 $\hat{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.

3 Data

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.

3.1 The Sample

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

HE0435$-$1223: The HE 0435$-$1223 system (J2000: 04 ${}^{\textrm{h}}$38${}^{\textrm{m}}$14$\aas@@fstack{s}$9, −12°17′14$\aas@@fstack{\prime\prime}$4) is a quadruply-lensed quasar discovered by Wisotzki et al. (2002). The main lensing galaxy is at a redshift of $z_{\ell}=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 $\sigma=471\pm 100\rm 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 $\sigma=222\pm 15~{}\rm km\,s^{-1}$. The time delays of this system were measured by Bonvin et al. (2017) with $\sim 6.5\%$ uncertainties. 2. 2.

PG1115+080: This four-image system was the second strong gravitational lens system to be discovered (Weymann et al., 1980). The system is located at 11h18m16$\aas@@fstack{s}$899 +07°45′58$\aas@@fstack{\prime\prime}$502 (J2000). The background quasar with a redshift of $z_{s}$ =1.722 is lensed by a galaxy with $z_{\ell}$ = 0.3098 (Henry & Heasley, 1986; Christian et al., 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 $\sigma=390\pm 60\rm km\,s^{-1}$ (Wilson et al., 2016). Tonry (1998) measured the stellar velocity dispersion of the lensing galaxy to be $\sigma=281\pm 25~{}\rm km\,s^{-1}$. The time delays of this system were measured by Bonvin et al. (2018) with $\sim 8.5\%$ uncertainties. 3. 3.

RXJ1131$-$1231: The RXJ 1131$-$1231 system (J2000: 11${}^{\textrm{h}}$31${}^{\textrm{m}}$52${}^{\textrm{s}}$, $−$12°31′59″) is a quadruply-lensed quasar discovered by Sluse et al. (2003). The spectroscopic redshifts of the lensing galaxy and the backgound source are at $z_{\ell}=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 $\sigma=323\pm 20\rm km\,s^{-1}$. The time delays were measured by Tewes et al. (2013b) with $\sim 1.5\%$ uncertainties.

3.2 Keck Adaptive Optics Imaging

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 prep). The targets were observed either with the narrow camera setup, which provides a roughly 10$\times$10″ field of view and a pixel scale of 9.942 milliarcsec (mas), or the wide camera that gives a roughly 40$\times$40″ field of view and 39.686 mas pixels. 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 coadditon of the exposures. During the distortion correction step, the images are resampled to produce final pixel scales of 10 mas pix*-1* for the narrow camera and 40 mas pix*-1* for the wide camera. The narrow camera pixels oversample the PSF, which has typical FWHM values of 60–90 mas. Therefore, to improve the modeling efficiency for the narrow camera data, we perform a 2$\times$2 binning of the images produced by the pipeline to obtain images that have a 20 mas pix*-1* scale. Further details on RXJ 1131$-$1231, which was observed with the NIRC2 wide camera, can be found in Chen et al. (2016).

3.3 Hubble Space Telescope Imaging

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) were analyzed in 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 modeled 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 modeling procedure on the AO and HST data for PG 1115$+$080.

3.4 MPIA 2.2 m Imaging

The contribution of the line-of-sight 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 analyzed 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 coadded 95 images of the field taken with the Wide Field Imager (Baade et al., 1999) mounted at the Cassegrain focus of the MPIA 2.2m telescope. The camera provides a pixel scale of $0\aas@@fstack{\prime\prime}238$. 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 Sec. 2.1 of Bonvin et al., 2018, for details); each image has been exposed for 330 seconds and covers a field of view of $\sim 8^{\prime}\times 16^{\prime}$. 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 coadding the exposures. The final coadded image has an effective seeing of $0\aas@@fstack{\prime\prime}86$.

4 Lens models

In this section, we describe the models that we use for fitting the hight 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 modeling 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), and reconstructed the AO PSF by using the PSF correction method developed in Chen et al. (2016).

4.1 Mass models

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

• •

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., 2006; Koopmans 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 ($\gamma^{\prime}$), Einstein radius ($\theta_{\rm E}$), and the axis ratio of the elliptical isodensity contour ($q$).

• •

Composite: we follow Suyu et al. (2014) and test a composite (baryonic + dark matter) model. The baryonic component is modeled 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 et al., 1996; Golse & Kneib, 2002) with the following parameters: halo normalization ($\kappa_{s}$), halo scale radius ($r_{s}$), and halo minor-to-major axis ratio ($q$) as well as associated position angle ($\theta_{q}$). Note that the ellipticity is implemented in the potential for the dark matter.

• •

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.

4.2 Lens light models

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

• •

2S$\acute{\text{e}}$rsic: we model the light distribution of the lens galaxy with two concentric elliptical S$\acute{\text{e}}$rsic profiles. For all three lens systems we found that a single Sérsic profile was insufficient for modeling the light distributions.

• •

2Chameleon: the chameleon profile is the difference of two isothermal profiles. It mimics a S$\acute{\text{e}}$rsic profile and enables computationally efficient lens modeling (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 modeling the composite model described in Section 4.1.

4.3 Strategies for mitigating the Mass-sheet Transformation

The mass-sheet transformation is a known degeneracy in lens modeling, in which one can transform a projected mass distribution, $\kappa(\theta)$, into infinite sets of $\kappa_{\lambda}(\theta)$ via

[TABLE]

without degrading the fit to the imaging. The corresponding time-delay distance changes via

[TABLE]

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.

4.3.1 Mass along the line of sight

If we perfectly know the true $\kappa(\theta)$, then the role of $\lambda$ in Equation (7) can be understood by looking at the behaviour far from the lensing galaxy, i.e., as $\theta\rightarrow\infty$. In this regime, the mass distribution of the lensing galaxy, $\kappa(\theta)$, approaches 0 and, thus, $\lambda$ can be interpreted as a constant-density mass sheet contributed by the lens environment. It is exactly the first-order form produced by line-of-sight (LOS) structure when its effect is small. Therefore, in our models we identify $\lambda$ with $\kappa_{\textrm{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):

[TABLE]

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.

• •

Shear: Shear distorts the lensed image shapes and thus it can be detected in the modeling process. We express the lens potential in polar coordinates ($\theta$, $\varphi$) to model the external shear on the imaging plane:

[TABLE]

where $\gamma_{\text{ext}}$ is the shear strength and $\varphi_{\text{ext}}$ is the shear angle. The shear position angle of $\varphi_{\text{ext}}=0^{\circ}$ corresponds to a shearing along $\theta_{1}$ whereas $\varphi_{\text{ext}}=90^{\circ}$ corresponds to shearing along $\theta_{2}$.333Our (right-handed) coordinate system $(\theta_{1},\theta_{2})$ has $\theta_{1}$ along the East-West direction and $\theta_{2}$ along the North-South direction.

• •

Millennium Simulation: Due to the MST, the lens images do not provide direct information on $\kappa_{\textrm{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/under-density of galaxies in the lens fields. They then selected lines of sight of similar over/under-density from the Millennium Simulation, with their corresponding values of the convergence, thus producing a probability distribution for $\kappa_{\rm ext}$: $P(\kappa_{\textrm{ext}}|\bm{d}_{\mathrm{ENV}})$. Suyu et al. (2013) later used, in addition to the number counts, the shear value inferred from lens modeling as an additional constraint on $\kappa_{\textrm{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. (2014, 2017) give a quantitative term, flexion shift ($\Delta_{3}x$), which estimates the deviations in lensed image positions due to third-order (flexion) terms, and suggest that if $\Delta_{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)-percent level. We follow this convention to choose which galaxies are included in our models. Note that the $\Delta_{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 provide more constraining power than the point-source imaging that was used to set the threshold.

4.3.2 Lens Kinematics

Conversely, even if we perfectly know the mass along the line of sight, the value of $\lambda$ remains uncertain because we do not know the true $\kappa(\theta)$ 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 (SPT;, Schneider & Sluse, 2014; Wertz et al., 2018; Wertz & Orthen, 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 $\kappa$ within the Einstein ring is preserved. However, the stellar velocity dispersion is sensitive to the integrated value of $\kappa$ 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 modeling of RXJ 1131$-$1231, Suyu et al. (2014), have shown that by including the velocity dispersion, one can obtain a robust $D_{\Delta t}$ when considering both power-law and composite model. 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 modeling to mitigate this internal mass-profile degeneracy.

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

• •

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

[TABLE]

where $\rho_{0}r_{0}^{n}$ and $F_{n}(r)$ are the normalization and the mass density distribution. By integrating $\rho_{\textrm{local}}$ within a cylinder with radius given by the Einstein radius, $R_{\textrm{Ein}}$, one obtains

[TABLE]

where $M_{\textrm{2D}}(R_{\textrm{Ein}})$ is the projected mass within $R_{\textrm{Ein}}$. The mass contained in $M_{\textrm{local}}$ is

[TABLE]

where $M_{\textrm{ext}}$ represents the mass contribution from $\kappa_{\textrm{ext}}$ and

[TABLE]

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

[TABLE]

Substituting this in Equation (11), we obtain

[TABLE]

Although there is $(1-\kappa_{\textrm{ext}})$ in Equation (16), the normalization of the local mass density distribution remains invariant (Yıldırım et al., 2019) as $\Sigma_{\textrm{cr}}$ can be re-expressed as

[TABLE]

where $(1-\kappa_{\textrm{ext}})$ term cancels out in Equation (16).

• •

An anisotropy component: we assume the anisotropy component in the form of an anisotropy radius, $r_{\textrm{ani}}$, in the Osipkov-Merritt formulation (Osipkov, 1979; Merritt, 1985),

[TABLE]

where $r_{\textrm{ani}}=0$ is pure radial orbits and $r_{\textrm{ani}}\rightarrow\infty$ is isotropic with equal radial and tangential velocity dispersions.

• •

A stellar component: we assume a Hernquist profile (Hernquist, 1990),

[TABLE]

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_{\textrm{eff}}$. For the composite model, the stellar component is represented by the light profile multiplied by a constant mass-to-light ratio.

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

[TABLE]

given the $\kappa_{\textrm{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_{\textrm{ani}}$ over a uniform distribution $[0.5,5]r_{\textrm{eff}}$. To compare with the data, we can get the seeing-convolved luminosity-weighted line-of-sight velocity dispersion,

[TABLE]

where $R$ is the projected radius, $I(R)$ is the light distribution, $\mathcal{P}$ is the PSF convolution kernel (Mamon & Łokas, 2005), and $\mathcal{A}$ is the aperture. The luminosity-weighted line-of-sight velocity dispersion is given by

[TABLE]

For a system with significant perturbers at a different redshift from the main lens (e.g., HE 0435$-$1223), we assume a flat $\Lambda$CDM cosmology with $H_{0}$ uniform in [0, 150]$~{}\rm km\,s^{-1}\,Mpc^{-1}$, $\Omega_{\textrm{m}}=0.3$, and $\Omega_{\textrm{m}}=1-\Omega_{\Lambda}$ to 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_{\ell s}$ and then combine with the measurement of $D_{\Delta t}$ to infer the value of $D_{\ell}$ without assuming any cosmological model (Birrer et al., 2016, 2019). The further advantage of this method is that $D_{\ell}$ is not affected by $\kappa_{\textrm{ext}}$ (Jee et al., 2015).

4.4 Microlensing time-delay prediction models

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 lamp-post 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_{\Delta t}$ measurements both with and without the lamp-post 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 wasn’t applied in previous H0LiCOW work to infer the final H0 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 (GPU) implementation of the inverse ray–shooting technique (Kayser et al., 1986). All magnification maps have dimension of $8192\times 8192$ pixels over a scale of $20\langle R_{\textrm{E}}\rangle$, where

[TABLE]

We choose the Salpeter initial mass function with mean mass $\langle M_{\star}\rangle=0.3M_{\sun}$ and the ratio between the upper and lower masses $M_{\rm upper}/M_{\rm lower}=100$ (Kochanek, 2004). We consider a standard thin disk model (Shakura & Sunyaev, 1973). Given the disk 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 which are used to estimate the probability distribution of the microlensing time-delay for each system are listed in Table 2 and Table 3.

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

[TABLE]

where the first term on the right hand side is the same as in Equation (2) and $t_{i}-t_{j}$ is the extra delay caused by the MTDE between images $i$ and $j$ (see details in Chen et al., 2018a).

5 lens modeling

Both HE 0435$-$1223 and RXJ 1131$-$1231 have been extensively modeled using the extended lensed emission seen in high-resolution HST imaging of the systems (Wong et al., 2017; Suyu et al., 2014), but our modeling techniques have not yet been applied to PG 1115$+$080. Therefore, in this section we begin with a description of the modeling 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 modeled (Wong et al., 2017; Suyu et al., 2014), and then combine the two modeling 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_{\Delta t}$, $D_{\ell}$ (if computed), and $H_{0}$ were kept blind until all coauthors came to a consensus to reveal the values during a collaboration telecon on June 5th. 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_{\Delta t}$ and $D_{\ell}$.

To better control the systematics due to the choice of lens modeling 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 modeling when using a modeling code, such as GLEE, that implements the pixelated-source reconstruction technique.

5.1 PG1115+080 Modeling

PG 1115$+$080 is a single-plane lens system embedded in a nearby group which consists of 13 known galaxies (see the solid circles in Figure 2). If we model the lens without including the group, the mass profile shows a very steep slope ($\gamma\sim 2.35$; note that $\gamma=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 ($\gamma_{\textrm{ext}}\sim 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) indicates 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.

5.1.1 The PSF of PG1115

For the HST imaging, we use Tinytim (Krist & Hook, 1997) to generate the PSFs with different spectral index, $\alpha$, of a power-law from -0.4 to -2.5 and different focuses444The flux per unit frequency interval is F$\nu=C\nu^{\alpha}$, where $\alpha$ is the power-law index and C is a constant; focus is related to the breathing of the second mirror, which is between $0\sim 10$. 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 modeling the HST imaging. For the AO imaging, we follow the criteria described in Section 4.4.3 in (Chen et al., 2016) and perform 8 iterative steps to create the final PSF and make sure the size of the PSF for convolution is large enough so that the results is stable. The full width half maximum (FWHM) of the AO PSF is 0.07″, while the FWHM of the HST PSF is 0.15″. We show the reconstructed AO PSF in Figure 22.

5.1.2 Main Lens

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

• •

SPEMD+2S$\acute{\text{e}}$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$\acute{\text{e}}$rsic profile is not sufficient to describe the light distribution, so we model it with two concentric elliptical S$\acute{\text{e}}$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 center of mass is very close to the center of light, with $|\Delta r|\approx 0.015^{\prime\prime}$. This implies that the offset between the projected center of dark matter and baryonic matter is small.

• •

Composite+2Chameleon+shear: We also model the main lens with composite model. Because the SPEMD+2S$\acute{\text{e}}$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, whilst the relative amplitudes are fixed in the mass profiles.

5.1.3 Nearby Group

Based on the velocity dispersion of the nearby group, $\sigma_{\textrm{group}}=390\pm 60~{}\rm km\,s^{-1}$, the inferred group mass is around $10^{13}\sim 10^{14}h^{-1}M_{\odot}$ (Wilson et al., 2016). Oguri (2006) has shown that in this mass regime the mass profile is too complicated to be described by a either a simple NFW profile or SIS profile, as it is a transition between the two. Thus, we use a NFW profile as the fiducial model, but also model the group with a SIS profile as a systematic check. In the following, we show how we determine reasonable priors on the NFW and SIS profiles.

• •

group (NFW): we follow Wong et al. (2011) and use $M_{\textrm{vir}}-c_{\textrm{vir}}$ relationships based on WMAP5 results in Macciò et al. (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_{\textrm{vir}}-c_{\textrm{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_{\textrm{vir}}$. Then we can obtain the concentration, $c_{\textrm{vir}}$, from the $M_{\textrm{vir}}$-$c_{\textrm{vir}}$ relationship assuming a reasonable scatter of 0.14 in $\textrm{log}~{}c_{\textrm{vir}}$ (Bullock et al., 2001; Wong et al., 2011). With the critical density and the characteristic overdensity at the lens redshift (Eke et al., 1998, 2001), we can obtain $r_{\textrm{vir}}$ and a prior probability on the scale radius, $r_{s}$, via $c_{\textrm{vir}}=r_{\textrm{vir}}/r_{s}$. The prior probability distribution of the normalization can be calculated by combining $r_{s}$, the central density of the halo ($\rho_{0}$), and the critical surface density for lensing ($\Sigma_{\textrm{cr}}$).

• •

group (SIS): we convert the velocity dispersion to an Einstein radius via

[TABLE]

to get a prior on $\theta_{\rm E}$ of $1.4\arcsec\pm 0.2\arcsec$ for G1 and $0.4\arcsec\pm 0.4\arcsec$ for G2.

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

5.1.4 Nearby Perturbing Galaxies: G1 and G2

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 $\Delta_{3}x$ of the nearby galaxies using the notation and definition in McCully et al. (2017). As $\Delta_{3}x$ is expressed in terms of the Einstein radius of these perturbers, we convert the measured velocity dispersions of G1 and G2 ($250\pm 20~{}\rm km\,s^{-1}$ and $130\pm 60~{}\rm 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 $\Delta_{3}x$), and use their relative luminosities compared to either G1 and 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\Delta_{3}x(\textrm{G1})=-3.68^{+0.13}_{-0.14}$ (in units of log(arcseconds)) and $\log\Delta_{3}x(\textrm{G2})=-4.01^{+0.75}_{-1.07}$, whereas the remaining galaxies have $\log\Delta_{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 modeled as SIS mass distributions.

5.1.5 LOS Analysis and the External Convergence

The technique of inferring $P(\kappa_{\textrm{ext}}|\bm{d}_{\mathrm{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 multi-band 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 $\kappa_{\textrm{ext}}$ (Collett et al., 2013). For our $\kappa_{\rm ext}$ analysis for PG 1115$+$080 we therefore use the deep coadded data set from the MPIA 2.2 m Telescope described in §3.4.

The coadded image has a limiting magnitude of $25.36\pm 0.08$, deeper than the control survey, CFHTLenS ($r=24.88\pm 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 color terms. Fortunately, a cross-match of the detections in our field with the ones in SDSS555We ignore the negligible differences between the SDSS and CFHTLenS $r$-band filters. shows that, after correcting for the zero-point offset, the scatter is small, with an rms of $\sim 0.10$ mag, which we add to the photometric error budget. However, we choose a brighter magnitude limit $r\leq 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, $r-i\sim 0.5$ for our cross-matches with SDSS. Due to the good seeing of our data, we find a clear stellar locus which allows us to determine a good classification threshold for FLUX_RADIUS, using the methodology in Coupon et al. (2009). We show the $240\arcsec\times 240\arcsec$ cutout of the field of view (FOV) in Figure 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\arcsec$- and $120\arcsec$-radius apertures, using the technique from Rusu et al. (2017) and the galaxy catalogue produced above, 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 $\kappa_{\textrm{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 4. 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(\kappa_{\textrm{ext}}|\bm{d}_{\mathrm{ENV}},\gamma)$ 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 Figure 4 and the summary table in Appendix (C).

5.1.6 Systematics Tests and Unblinding Results

We summarize below the choices that we explore for the mass modeling, 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 center of its light distribution. When modeling the galaxy group as a NFW profile, we use the $M_{\rm vir}-c_{\rm vir}$ from Macciò et al. (2008), based on a WMAP5 cosmology (we found that the impact of using different cosmology is negligible.).

$\bullet$

SPEMD + 2S$\acute{\text{e}}$rsic + NFW group.

$\bullet$

SPEMD + 2S$\acute{\text{e}}$rsic + NFW group + G1

$\bullet$

SPEMD + 2S$\acute{\text{e}}$rsic + SIS group + G1

$\bullet$

SPEMD + 2S$\acute{\text{e}}$rsic + NFW group + G1 + G2

$\bullet$

Composite lens + NFW group

$\bullet$

Composite lens + NFW group + G1

$\bullet$

Composite lens + SIS group + G1

$\bullet$

Composite lens + NFW group + G1 + G2

$\bullet$

For all of the above models, we test five different source resolutions.

To address the MST, we use the $P(\kappa_{\textrm{ext}}|\bm{d}_{\mathrm{ENV}},\gamma)$ from the previous section and importance sample the measured velocity dispersion, $\sigma=281\pm 25~{}\rm km\,s^{-1}$ (Tonry, 1998), which was obtained inside a 1.0-arcsec2 aperture with a seeing of $0.8\arcsec$. When sampling $D_{\Delta t}$ with the time-delay measurements from Bonvin et al. (2018), we follow Chen et al. (2018a) to sample the $D_{\Delta t}$ both with and without considering the MTDE. The parameters for estimating the star density and accretion disk model can be found in Chen et al. (2018a). In sum, we explore 160 modeling 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 disk sizes (or we ignore the MTDE). We show the AO imaging reconstruction in Figure 5 and HST imaging reconstruction in Figure 6. In Figure 7 and Figure 8, we present the posteriors of the important parameters of AO-only results and HST+AO results for the composite model and the power-law model, respectively. The $D_{\Delta t}$ distributions from the 160 different model choices are shown in Figure 9. The final $D_{\Delta t}$-$D_{\ell}$ without considering MTDE is shown in Figure 10. We show the $\Delta$BIC value of each model without considering MTDE in Appendix (F). We found that the $D_{\Delta t}$ inferred from various model choices are statistically consistent. The uncertainties of the final marginalized $D_{\Delta t}$ are $\sim 9.6\%$ and $D_{\ell}$ are $\sim 29\%$.

5.2 HE0435-1223 modeling

For the HE 0435$-$1223 system we model only the AO data, and then combine the results with the HST-based modeling 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 $\Delta_{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 $\Delta_{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 multi-plane lens equation (e.g., Blandford & Narayan, 1986; Schneider et al., 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 time delay distance, $D_{\Delta t}^{\textrm{eff}}\left(z_{\ell},z_{s}\right)$, which is fairly robust to changes in the assumed cosmology.

5.2.1 The AO PSF of HE 0435$-$1223

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″(see Figure 22).

5.2.2 Lens Model Choices

As we did for PG 1115$+$080, we model the main lens with either a SPEMD or composite model. For the composite model, we follow Wong et al. (2017) and set the Gaussian prior for the scale radius to $4.3\arcsec\pm 2.0\arcsec$ based on scaling relations derived from the SLACS sample (Gavazzi et al., 2007). The most massive perturber, G1, is modeled as a SIS profile. When modeling 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.

5.2.3 LOS Analysis and the External Convergence

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(\kappa_{\textrm{ext}}|\bm{d}_{\mathrm{ENV}},\gamma)$ are presented in Rusu et al. (2017). In this work, we use the shear values determined from our lens modeling 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 datasets. We show our updated results in Figure 12.

5.2.4 Systematics Tests and Unblinding results

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 center of its light distribution.

$\bullet$

A power-law model plus G1 as a SIS.

$\bullet$

A composite model plus G1 as a SIS.

$\bullet$

A composite model plus the five perturbers (G1 – G5).

$\bullet$

For all of the above models, we test five different source resolutions.

We show an example of the AO imaging reconstruction in Figure 11 and present the posteriors of $D_{\Delta t}$ in Figure 15.

To assess the MST, for each model we perform the importance sampling given the measured velocity dispersion, $\sigma=222\pm 15~{}\rm km\,s^{-1}$ inside a $0.54\arcsec\times 0.7\arcsec$ aperture with a seeing of $0.8\arcsec$ (Wong et al., 2017).

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 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).

5.3 RXJ1131-1231 modeling

A detailed discussion of the RXJ 1131$-$1231 lens modeling 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$\acute{\text{e}}$rsic profiles to model the lens light. The satellite galaxy of the main deflector was modeled as an SIS profile. We modeled only the lensing galaxy plus satellite, and did not consider $\kappa_{\rm ext}$. The modeling 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.

5.3.1 Main Lens and Satellite

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\arcsec\pm 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$\acute{\text{e}}$rsic profile, and the satellite mass as a 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.

5.3.2 LOS Analysis and the External Convergence

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(\kappa_{\textrm{ext}})$, $P(\kappa_{\textrm{ext}}|\gamma)$, and $P(\kappa_{\textrm{ext}}|\bm{d}_{\mathrm{ENV}},\gamma)$. Here $\gamma$ is the external shear required by the mass models of the main lensing galaxy, while $\bm{d}_{\mathrm{ENV}}$ is the relative overdensity of galaxies within a 45″ aperture that is centred on the lens. This overdensity, $\zeta_{1}^{45\arcsec}=1.4\pm 0.05$ (following the notation in Birrer et al. (2019)), is calculated from galaxies with apparent HST/ACS F814W magnitudes $18.5\leq m\leq 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 $\kappa_{\rm ext}$.

5.3.3 Systematics Tests

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.

$\bullet$

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

$\bullet$

A composite model.

$\bullet$

For the above models, we test five different source resolutions.

We use the observed velocity dispersion, $323\pm 20\rm km\,s^{-1}$ (Suyu et al., 2013) given the $\kappa_{\textrm{ext}}$ in Section 5.3.2 to sample $D_{\Delta t}$ and $D_{\ell}$ without assuming cosmology. We plot the posteriors of $D_{\Delta t}$ and $D_{\ell}$ in Figure 18.

6 Cosmological inference

We present the cosmological inferences based on our sample of three gravitational lenses that have AO imaging data, HST imaging data, velocity dispsersion measurements, line-of-sight 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.

6.1 Cosmological inference from AO Strong Lensing

Figure 20 presents the marginalized posterior PDF for $H_{0}$ assuming a flat $\Lambda$CDM model which has a uniform prior on $H_{0}$ in the range [0, 150] $\rm km\,s^{-1}\,Mpc^{-1}$ and a uniform prior on $\Omega_{\textrm{m}}$ in the range of [0.05, 0.5]. After unblinding, we find that PG 1115$+$080 AO imaging yields $H_{0}=82.8\begin{subarray}{c}+9.4\\ -8.3\end{subarray}~{}\rm km\,s^{-1}\,Mpc^{-1}$; HE 0435$-$1223 AO imaging yields $H_{0}=70.1\begin{subarray}{c}+5.3\\ -4.5\end{subarray}~{}\rm km\,s^{-1}\,Mpc^{-1}$. RXJ 1131$-$1231 AO imaging yields $H_{0}=77.0\begin{subarray}{c}+4.0\\ -4.6\end{subarray}~{}\rm km\,s^{-1}\,Mpc^{-1}$. The joint analysis of three AO lenses yields $H_{0}=75.6\begin{subarray}{c}+3.2\\ -3.3\end{subarray}~{}\rm km\,s^{-1}\,Mpc^{-1}$.

6.2 Cosmological Inference from AO and HST Strong Lensing imaging

As the AO imaging of HE 0435$-$1223 and RXJ 1131$-$1231 are modeled 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 are used only once. We can express the joint posterior as

[TABLE]

where $P(\bm{\eta}|\bm{d}_{\textrm{HST}},\sigma,\bm{d_{\textrm{ENV}}},\textrm{H},\textrm{MS})$ reflects the previous modeling results, $P(\bm{\eta}|\bm{d}_{\textrm{AO}},\textrm{H},\textrm{MS})$ is obtained by this work, $P(\bm{\eta}|\textrm{H},\textrm{MS})$ is the prior used in both analysis, H is the lens model, and MS is the Millennium Simulation.

Figure 21 presents the marginalized posterior PDF for $H_{0}$ assuming flat-$\Lambda$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\begin{subarray}{c}+4.7\\ -4.6\end{subarray}~{}\rm km\,s^{-1}\,Mpc^{-1}$; the joint AO+HST of PG 1115$+$080 implies $H_{0}=81.1\begin{subarray}{c}+7.9\\ -7.1\end{subarray}~{}\rm km\,s^{-1}\,Mpc^{-1}$. The joint AO+HST of RXJ 1131$-$1231 implies $H_{0}=78.3\begin{subarray}{c}+3.4\\ -3.3\end{subarray}~{}\rm km\,s^{-1}\,Mpc^{-1}$. The combination of three AO+HST lenses yields $H_{0}=76.8\begin{subarray}{c}+2.6\\ -2.6\end{subarray}~{}\rm 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 sources of uncertainty of RXJ 1131$-$1231 is the line-of-sight mass distribution.

7 Conclusions

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 line-of-sight studies from deep wide-area spectroscopic as well as photometric data, and the time-delay measurements from state-of-the-art light curve 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 and 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 $\Lambda$CDM model which has a uniform prior on $H_{0}$ in the range [0, 150] $\rm km\,s^{-1}\,Mpc^{-1}$ and a uniform prior on $\Omega_{\textrm{m}}$ in the range of [0.05, 0.5]. After unblinding, PG 1115$+$080 AO imaging yields $H_{0}=82.8\begin{subarray}{c}+9.4\\ -8.3\end{subarray}~{}\rm km\,s^{-1}\,Mpc^{-1}$; HE 0435$-$1223 AO imaging yields $H_{0}=70.1\begin{subarray}{c}+5.3\\ -4.5\end{subarray}~{}\rm km\,s^{-1}\,Mpc^{-1}$. RXJ 1131$-$1231 AO imaging yields $H_{0}=77.0\begin{subarray}{c}+4.0\\ -4.6\end{subarray}~{}\rm km\,s^{-1}\,Mpc^{-1}$. The joint analysis of three AO lenses yields $H_{0}=75.6\begin{subarray}{c}+3.2\\ -3.3\end{subarray}~{}\rm km\,s^{-1}\,Mpc^{-1}$.

The joint AO+HST of PG 1115$+$080 yields $H_{0}=81.1\begin{subarray}{c}+7.9\\ -7.1\end{subarray}~{}\rm km\,s^{-1}\,Mpc^{-1}$. The joint AO+HST of HE 0435$-$1223 yields $H_{0}=71.6\begin{subarray}{c}+4.7\\ -4.6\end{subarray}~{}\rm km\,s^{-1}\,Mpc^{-1}$. The joint AO+HST of RXJ 1131$-$1231 yields $H_{0}=78.3\begin{subarray}{c}+3.4\\ -3.3\end{subarray}~{}\rm km\,s^{-1}\,Mpc^{-1}$. The combination of three AO+HST lenses yields $H_{0}=76.8\begin{subarray}{c}+2.6\\ -2.6\end{subarray}~{}\rm 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.

Acknowledgements

G. C.-F. Chen acknowledges support from the Ministry of Education in Taiwan via Government Scholarship to Study Abroad (GSSA). C.D.F. and G. C.-F. Chen acknowledge support for this work from the National Science Foundation under Grant No. AST-1715611. S.H. acknowledges support by the DFG cluster of excellence ‘Origin and Structure of the Universe’ (www.universe-cluster.de). L. V. E. K. acknowledges the support of an NWO-VIDI grant (nr. 639.043.308).

This work was supported by World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan. K.C.W. is supported in part by an EACOA Fellowship awarded by the East Asia Core Observatories Association, which consists of the Academia Sinica Institute of Astronomy and Astrophysics, the National Astronomical Observatory of Japan, the National Astronomical Observatories of the Chinese Academy of Sciences, and the Korea Astronomy and Space Science Institute. SHS thanks the Max Planck Society for support through the Max Planck Research Group. J. Chan, V. Bonvin and F. Courbin acknowledge support from the Swiss National Science Foundation (SNSF). This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (COSMICLENS: grant agreement No 787866). X.D. and T.T. acknowledges support from the Packard Foundation through a Packard Research Fellowship and from the NSF through grant AST-1714953. G. C.-F. Chen thank Jen-Wei Hsueh for the useful discussions and feedback. The data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California, and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.

This work is based on observations obtained with MegaPrime/MegaCam, a joint project of CFHT and CEA/IRFU, at the Canada-France-Hawaii Telescope (CFHT) which is operated by the National Research Council (NRC) of Canada, the Institut National des Sciences de l’Univers of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. This research used the facilities of the Canadian Astronomy Data Centre operated by the National Research Council of Canada with the support of the Canadian Space Agency. CFHTLenS data processing was made possible thanks to significant computing support from the NSERC Research Tools and Instruments grant program.

This work used computational and storage services associated with the Hoffman2 Shared Cluster provided by UCLA Institute for Digital Research and Education’s Research Technology Group. Data analysis was in part carried out on common use data analysis computer systems at the Astronomy Data Center, ADC, of the National Astronomical Observatory of Japan.

Appendix A Comparison of the reconstructed AO PSF from three lenses

We show the comparison of the reconstructed AO PSF in Figure 22. The 2D plots of the three PSF clearly show that the atmosphere disturbance produce various structures of PSF with a core plus a wing, while the radial average intensity indicates that the intensity gradient of the cores are very similar inside 0.1″.

Appendix B Degeneracy between AO PSF wing and lens light

The PSF-reconstructed method developed in Chen et al. (2016) allow us to model the AO imaging down to the noise level and recover $D_{\Delta t}$ when we adopt the power-law model. However, we found that there exist a degeneracy between AO PSF wing and the lens light if we reconstruct the AO PSF from the AO imaging only. While it does not cause a problem for the power-law model since the power-law mass is constrained by the arc, the composite model instead could yield a baised $H_{0}$ because we assume the baryonic matter distribution follow the lens light. Thus, in the case of HE 0435$-$1223 and RXJ 1131$-$1231 where we do not reconstruct the PSF simultaneously from HST and AO imaging, we fix the mass distribution of the baryonic matter to that baryonic mass inferred from HST imaging.

In the case of PG 1115$+$080, we found that the degeneracy can be broken by simultaneously modeling the AO imaging and HST imaging as they share the same mass model. The additional constraint allow us to better characterize the AO PSF wing. Thus, given the same PSF, both power-law model and composite model can be modeled down to the noise level. In the future, 2D kinematic data could also potentially further break the degeneracy between the PSF wing and lens light, as it provides the information to characterize the dark and baryonic matter content (Cappellari et al., 2013).

Appendix C Constraints used to estimate $\kappa_{\textrm{ext}}$

Here we show the $\kappa_{\textrm{ext}}$ distribution of HE 0435$-$1223 in Figure 12 and present the summary table of the constraints used to estimate $\kappa_{\textrm{ext}}$ for all three systems in Table 4.

77footnotetext: The five perturbers are indicated in Figure 3 from Wong et al. (2017). Three additional galaxies enter the 12″-radius inner mask, and we slightly boosted their distance from the lens, in order to avoid masking them.

Appendix D Accounting for the missing galaxy group members in PG1115+080 due to spectroscopic incompleteness

The spectroscopic coverage of the FOV around PG 1115$+$080 is incomplete. Down to $R_{c}\leq 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 which are also part of the galaxy group associated with the lensing galaxy, but which are missed due to spectroscopic incompleteness. 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 $\kappa_{\mathrm{ext}}$ at the location of the lens, using the weighted number counts approach, we must remove the galaxies which 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 which are known to be part of the group can easily be removed, we must also account 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 part of the group. 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\pm 5$ missing galaxies (median and 16th, 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 centered on the lens to be $1^{+3}_{-1}$. We plot the distributions of these numbers in Figure 23. 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\sim 21$ at $z\sim 0.3$, therefore significantly brighter than our limiting magnitude. We note, however, that the two techniques produced results which were in agreement for a different lens, described in Rusu et. al., submitted.

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 of galaxies from our catalogue of galaxies inside the $120\arcsec$ apertures.

Appendix E Combining the independent measurements

Consider that we have following sets of information, AO imaging ($\bm{d}_{\textrm{AO}}$), HST imaging ($\bm{d}_{\textrm{HST}}$), velocity dispersion ($\sigma$) and environment data ($\bm{d_{\textrm{ENV}}}$), which can be used to constrain the difference of the fermat potential. We already have two independent measurements, $P(\bm{\eta}|\bm{d}_{\textrm{AO}},\textrm{H},\textrm{MS})$ and $P(\bm{\eta}|\bm{d}_{\textrm{HST}}$, $\sigma$, $\bm{d_{\textrm{ENV}}}$,$\textrm{H},\textrm{MS})$, The goal is to combine them and obtain the joint constraint, $P(\bm{\eta}|\bm{d}_{\textrm{HST}},\bm{d}_{\textrm{AO}},\sigma$, $\bm{d_{\textrm{ENV}}},\textrm{H},\textrm{MS})$, where $\bm{\eta}=(\phi_{ij}^{\textrm{H}},\kappa_{\textrm{ext}},\gamma_{\textrm{ext}},\bm{\zeta})$, $\phi_{ij}^{\textrm{H}}$ are the difference of the model fermat potentials (not the true fermat potentail, $\phi_{ij}^{\textrm{True}}$) at imaging i and j which we are interested, $\kappa_{\textrm{ext}}$ is the convergence from the line of sight, $\gamma_{\textrm{ext}}$ is the external shear inferred from the imaging data, and $\bm{\zeta}$ are the other parameters which we want to marginalize upon. H is one kind of the lens models (i.g., the power-law model or the composite model) and MS is Millennium Simulation.

We start with the joint constraint and step by step link to the independent measurements. With Bayes’ theorem, the joint constraint can be expressed as

[TABLE]

where

[TABLE]

are the priors on the parameters of the mass model and MS. Note that Millennium Simulation naturally provides more $\kappa_{\textrm{ext}}$ which close to mean density and lens model implicitly assumes $\kappa_{\textrm{ext}}<1$, so $P(\kappa_{\textrm{ext}}|\textrm{H,MS})$ is a non-flat prior with an upper bound ($<1$).

The next step is to separate the datasets into [$\bm{d}_{\textrm{HST}},\sigma,\bm{d_{\textrm{ENV}}}$] and [$\bm{d}_{\textrm{AO}}$]. Since the data are all independent, we can write down

[TABLE]

Although $\bm{d}_{\textrm{AO}}$ do not have the direct constraint power on $\kappa_{\textrm{ext}}$, the shear value inferred from $\bm{d}_{\textrm{AO}}$ implicitly help constrain $\kappa_{\textrm{ext}}$. Furthermore, we leave $\kappa_{\textrm{ext}}$ in the last term in Equation (29) because we want to link to $P(\bm{\eta}|\bm{d}_{\textrm{AO}},\textrm{H},\textrm{MS})$ in the future steps. Equation (27) become

[TABLE]

where

[TABLE]

The last term in Equation (31) tells us that the shear value inferred from the lens imaging can also help us to further constrain the convergence because of the correlation between the $\gamma_{\textrm{ext}}$ and $\kappa_{\textrm{ext}}$ in MS.

Bayes theorem tells us that

[TABLE]

and

[TABLE]

By substituting Equation (32) and Equation (33) into Equation (30), finally we obtain

[TABLE]

The first term in the right hand side in Equation (34) is done by Wong et al. (2017), while the numerator in the second term is from AO data alone. Thus, based on Equation (34), in order to get the joint constraint, we need to multiply this two posteriors and divide it by the non-uniform priors (e.g., $P(\kappa_{\textrm{ext}}|\textrm{H,MS})$) used in both datasets. This is because we need to get rid of the extra constraining power from doubly using the same non-uniform priors. Therefore, the denominator in Equation (34) become $P(\kappa_{\textrm{ext}}|\textrm{H,MS})$. Note that $\phi_{ij}^{\textrm{H}}$ is the model fermat potential, but what we want to obtain is $P(\phi_{ij}^{\textrm{True}})$, which can be expressed as

[TABLE]

where

[TABLE]

Appendix F SUMMARY OF PG1115+080 LENS MODELS WITH RESPECT TO THE BIC VALUE

We present the BIC of the power-law models in Table 5 and composite models in Table 6.

Appendix G SUMMARY OF HE 0435$-$1223 LENS MODELS WITH RESPECT TO THE BIC VALUE

We present the BIC of the power-law models in Table 7 and composite models in Table 8.

Appendix H SUMMARY OF RXJ 1131$-$1231 LENS MODELS WITH RESPECT TO THE BIC VALUE

We present the BIC of the power-law models in Table 9 and composite models in Table 10.