H0LiCOW XII. Lens mass model of WFI2033-4723 and blind measurement of its time-delay distance and $H_0$

TL;DR

This paper models the gravitational lens WFI2033-4723 using 14 years of data to measure the time-delay distance and estimate the Hubble constant with high precision, accounting for systematic uncertainties.

Contribution

It provides a detailed lens mass model and a blind cosmographical analysis for WFI2033-4723, improving the measurement of the Hubble constant from gravitational lensing.

Findings

Effective time-delay distance: 4784 Mpc with 6.6% precision

Hubble constant estimate: 71.6 km/s/Mpc with ~5% uncertainty

Systematic effects from environment and line-of-sight included

Abstract

We present the lens mass model of the quadruply-imaged gravitationally lensed quasar WFI2033-4723, and perform a blind cosmographical analysis based on this system. Our analysis combines (1) time-delay measurements from 14 years of data obtained by the COSmological MOnitoring of GRAvItational Lenses (COSMOGRAIL) collaboration, (2) high-resolution $\textit{Hubble Space Telescope}$ imaging, (3) a measurement of the velocity dispersion of the lens galaxy based on ESO-MUSE data, and (4) multi-band, wide-field imaging and spectroscopy characterizing the lens environment. We account for all known sources of systematics, including the influence of nearby perturbers and complex line-of-sight structure, as well as the parametrization of the light and mass profiles of the lensing galaxy. After unblinding, we determine the effective time-delay distance to be $4784_{-248}^{+399}~\mathrm{Mpc}$, an average precision of $6.6\%$. This translates to a Hubble constant $H_{0} = 71.6_{-4.9}^{+3.8}~\mathrm{km~s^{-1}~Mpc^{-1}}$, assuming a flat $\Lambda$CDM cosmology with a uniform prior on $\Omega_\mathrm{m}$ in the range [0.05, 0.5]. This work is part of the $H_0$ Lenses in COSMOGRAIL's Wellspring (H0LiCOW) collaboration, and the full time-delay cosmography results from a total of six strongly lensed systems are presented in a companion paper (H0LiCOW XIII).

Tables

Table 1. Table 1: Lens Model Parameters

Parameter	BIC-weighted Marginalized Constraints
Singular Power Law Ellipsoid Model
${\theta }_E\left(″\right)$777Spherical-equivalent Einstein radius	${0.929}_{-0.016}^{+0.015}$
$q$	${0.79}_{-0.01}^{+0.01}$
${\theta }_q$ (∘)	${33.1}_{-0.9}^{+0.8}$
${\gamma }^{\prime }$	${1.95}_{-0.01}^{+0.02}$
${\gamma }_{\mathrm{ext}}$	${0.112}_{-0.004}^{+0.006}$
${\theta }_{\gamma }$ (∘)	${83.6}_{-2.1}^{+7.1}$
X ${\theta }_E\left(″\right)$	${0.001}_{-0.001}^{+0.001}$
G2 ${\theta }_E\left(″\right)$	${0.932}_{-0.062}^{+0.027}$
G2 $q$	${0.66}_{-0.02}^{+0.04}$
G2 ${\theta }_q$ (∘)	${38.5}_{-4.1}^{+4.4}$
Composite Model
Stellar M/L ($M_{\odot }/L_{\odot }$)888M/L within ${\theta }_E$ for rest-frame $V$ band. The point source component of the lens light is assumed to be from low-level AGN emission as opposed to stellar light and is not included in the calculation. The given uncertainties are a combination of statistical effects and a systematic uncertainty equal to the difference between the calculated M/L with and without the point source contribution. The stellar mass is calculated assuming $H_0=70\mathrm{km}{s}^{-1}{\mathrm{Mpc}}^{-1}$, ${\Omega }_m=0.3$, ${\Omega }_{\Lambda }=0.7$, but changes in the cosmology affect the M/L by a negligible amount.	${2.1}_{-0.2}^{+0.2}$
Chameleon1 $q$	${0.762}_{-0.003}^{+0.003}$
Chameleon1 ${\theta }_q$ (∘)	${23.3}_{-0.2}^{+0.2}$
Chameleon2 $q$	${0.771}_{-0.003}^{+0.002}$
Chameleon2 ${\theta }_q$ (∘)	${26.2}_{-0.5}^{+0.3}$
NFW ${\kappa }_{0,h}$	${0.147}_{-0.012}^{+0.002}$
NFW $r_s\left(″\right)$	${11.15}_{-0.09}^{+0.21}$
NFW $q$	${0.89}_{-0.01}^{+0.01}$
NFW ${\theta }_q$ (∘)	${72.9}_{-1.3}^{+0.8}$
${\gamma }_{\mathrm{ext}}$	${0.133}_{-0.001}^{+0.002}$
${\theta }_{\gamma }$ (∘)	${89.4}_{-0.4}^{+0.4}$
X ${\theta }_E\left(″\right)$	${0.018}_{-0.002}^{+0.002}$
G2 ${\theta }_E\left(″\right)$	${1.008}_{-0.004}^{+0.019}$
G2 $q$	${0.93}_{-0.01}^{+0.01}$
G2 ${\theta }_q$ (∘)	${39.4}_{-10.3}^{+1.6}$

Table 2. Table 2: Weighted galaxy count ratios ζ q ¯ ¯ subscript 𝜁 𝑞 \overline{\zeta_{q}} for WFI2033 − - 4723

	$45″$	$45″$	$45″$	$120″$	$120″$	$120″$
Weight $q$	fiducial	fiducial +	fiducial +	fiducial	fiducial +	fiducial +
		$z=0.49$ group	$z=0.49,0.66$ groups		$z=0.49$ group	$z=0.49,0.66$ groups
$1$	${1.44}_{-0.08}^{+0.06}$	${1.20}_{-0.10}^{+0.03}$	${1.17}_{-0.07}^{+0.04}$	${1.55}_{-0.11}^{+0.07}$	${1.34}_{-0.08}^{+0.05}$	${1.27}_{-0.08}^{+0.04}$
$z$	${1.70}_{-0.13}^{+0.09}$	${1.36}_{-0.07}^{+0.07}$	${1.34}_{-0.08}^{+0.07}$	${1.65}_{-0.12}^{+0.12}$	${1.37}_{-0.10}^{+0.09}$	${1.31}_{-0.09}^{+0.08}$
$M_{\star }$	${1.45}_{-0.28}^{+0.11}$	${0.78}_{-0.04}^{+0.10}$	${0.80}_{-0.03}^{+0.10}$	${2.39}_{-0.29}^{+0.11}$	${1.68}_{-0.09}^{+0.30}$	${1.48}_{-0.09}^{+0.27}$
$M_{\star }^2$	${1.42}_{-0.49}^{+0.26}$	${0.51}_{-0.06}^{+0.12}$	${0.54}_{-0.05}^{+0.13}$	${3.69}_{-0.78}^{+0.43}$	${2.12}_{-0.17}^{+0.76}$	${1.75}_{-0.14}^{+0.64}$
$M_{\star }^3$	${1.40}_{-0.65}^{+0.45}$	${0.33}_{-0.06}^{+0.13}$	${0.37}_{-0.06}^{+0.14}$	${5.75}_{-1.75}^{+0.95}$	${2.73}_{-0.38}^{+1.53}$	${2.12}_{-0.31}^{+1.18}$
$1/r$	${1.33}_{-0.08}^{+0.01}$	${1.07}_{-0.06}^{+0.04}$	${1.04}_{-0.05}^{+0.04}$	${1.55}_{-0.11}^{+0.07}$	${1.35}_{-0.09}^{+0.05}$	${1.25}_{-0.08}^{+0.05}$
$z/r$	${1.49}_{-0.14}^{+0.03}$	${1.17}_{-0.06}^{+0.06}$	${1.14}_{-0.05}^{+0.06}$	${1.59}_{-0.10}^{+0.09}$	${1.33}_{-0.08}^{+0.06}$	${1.24}_{-0.07}^{+0.05}$
$M_{\star }/r$	${1.69}_{-0.33}^{+0.34}$	${0.69}_{-0.03}^{+0.11}$	${0.75}_{-0.05}^{+0.10}$	${2.08}_{-0.13}^{+0.21}$	${1.53}_{-0.09}^{+0.29}$	${1.38}_{-0.08}^{+0.27}$
$M_{\star }^2/r$	${1.97}_{-0.64}^{+0.73}$	${0.49}_{-0.07}^{+0.10}$	${0.56}_{-0.08}^{+0.11}$	${3.24}_{-0.55}^{+0.52}$	${1.89}_{-0.20}^{+0.85}$	${1.52}_{-0.11}^{+0.78}$
$M_{\star }^3/r$	${2.07}_{-1.00}^{+0.93}$	${0.32}_{-0.07}^{+0.12}$	${0.38}_{-0.07}^{+0.14}$	${5.30}_{-1.49}^{+1.00}$	${2.53}_{-0.44}^{+1.78}$	${1.86}_{-0.30}^{+1.42}$
$M_{\star ,\mathrm{rms}}^2$	${1.19}_{-0.23}^{+0.11}$	${0.71}_{-0.04}^{+0.08}$	${0.74}_{-0.04}^{+0.08}$	${1.92}_{-0.21}^{+0.11}$	${1.46}_{-0.06}^{+0.24}$	${1.32}_{-0.05}^{+0.23}$
$M_{\star ,\mathrm{rms}}^3$	${1.12}_{-0.21}^{+0.11}$	${0.69}_{-0.04}^{+0.08}$	${0.72}_{-0.04}^{+0.08}$	${1.79}_{-0.20}^{+0.09}$	${1.40}_{-0.07}^{+0.22}$	${1.28}_{-0.06}^{+0.21}$
$M_{\star }^2/r_{,\mathrm{rms}}$	${1.40}_{-0.25}^{+0.24}$	${0.70}_{-0.06}^{+0.07}$	${0.75}_{-0.06}^{+0.07}$	${1.80}_{-0.16}^{+0.14}$	${1.38}_{-0.08}^{+0.28}$	${1.23}_{-0.04}^{+0.29}$
$M_{\star }^3/r_{,\mathrm{rms}}$	${1.27}_{-0.25}^{+0.17}$	${0.68}_{-0.05}^{+0.08}$	${0.73}_{-0.06}^{+0.07}$	${1.74}_{-0.18}^{+0.11}$	${1.36}_{-0.08}^{+0.27}$	${1.23}_{-0.07}^{+0.26}$
$M_{\star }/r^3$	${1.13}_{-0.27}^{+0.26}$	${0.46}_{-0.04}^{+0.05}$	${0.48}_{-0.03}^{+0.06}$	${1.86}_{-0.23}^{+0.12}$	${1.37}_{-0.16}^{+0.14}$	${1.25}_{-0.15}^{+0.10}$
$M_{\star }/r^2$	${1.41}_{-0.29}^{+0.41}$	${0.54}_{-0.02}^{+0.09}$	${0.58}_{-0.04}^{+0.10}$	${2.04}_{-0.18}^{+0.14}$	${1.47}_{-0.09}^{+0.11}$	${1.31}_{-0.10}^{+0.11}$
$\sqrt{M_{\star }}/r$	${1.44}_{-0.19}^{+0.15}$	${0.81}_{-0.03}^{+0.06}$	${0.83}_{-0.04}^{+0.06}$	${1.78}_{-0.11}^{+0.06}$	${1.40}_{-0.06}^{+0.05}$	${1.29}_{-0.07}^{+0.04}$
$\sqrt{M_h}/r$	${1.57}_{-0.39}^{+0.29}$	${1.08}_{-0.22}^{+0.57}$	${1.07}_{-0.21}^{+0.69}$	${1.81}_{-0.28}^{+0.09}$	${1.39}_{-0.06}^{+0.19}$	${1.27}_{-0.05}^{+0.18}$

Table 3. Table 3: Effective time-delay distance and BIC weighting for different lens models

Model	$D_{\Delta t}$ (Mpc)	$\Delta$BIC	BIC weight
SPEMD fiducial	${4640}_{-195}^{+238}$	17	0.674
SPEMD AGN mask + 1pix	${4631}_{-199}^{+247}$	49	0.235
SPEMD AGN mask weight=0	${4611}_{-195}^{+258}$	1984	0.000
SPEMD arcmask+1pix, 60x60 source	${4894}_{-190}^{+219}$	306	0.000
SPEMD lens group halo	${4732}_{-203}^{+228}$	16	0.686
SPEMD lens + z=0.49 group halo	${4740}_{-190}^{+251}$	0	1.000
SPEMD chameleon light profiles	${4312}_{-185}^{+241}$	1962	0.000
All SPEMD (BIC-weighted)	${4703}_{-203}^{+245}$	$-$	$-$
Composite fiducial	${4731}_{-246}^{+416}$	202	0.000
Composite AGN mask + 1pix	${4836}_{-235}^{+418}$	206	0.000
Composite AGN mask weight=0	${4743}_{-238}^{+417}$	637	0.000
Composite arcmask+1pix, 60x60 source	${4944}_{-258}^{+437}$	559	0.000
Composite lens group halo	${4779}_{-283}^{+475}$	16	0.770
Composite lens + z=0.49 group halo	${5009}_{-305}^{+471}$	0	1.000
All Composite (BIC-weighted)	${4913}_{-320}^{+489}$	$-$	$-$
All	${4784}_{-248}^{+399}$	$-$	$-$

Table 4. Table 4: SPEMD Model Parameters

Parameter	Marginalized Constraints
	Fiducial	AGNmask+1	AGNwht=0	Arc+1,60src	Group	Group + $z$=0.49 group	Chameleon
${\theta }_E\left(″\right)$131313Spherical-equivalent Einstein radius	${0.944}_{-0.004}^{+0.004}$	${0.943}_{-0.004}^{+0.004}$	${0.950}_{-0.005}^{+0.004}$	${0.933}_{-0.004}^{+0.004}$	${0.912}_{-0.003}^{+0.003}$	${0.927}_{-0.006}^{+0.005}$	${0.946}_{-0.003}^{+0.003}$
$q$	${0.80}_{-0.01}^{+0.01}$	${0.80}_{-0.01}^{+0.01}$	${0.81}_{-0.01}^{+0.01}$	${0.79}_{-0.01}^{+0.01}$	${0.79}_{-0.01}^{+0.01}$	${0.78}_{-0.01}^{+0.01}$	${0.78}_{-0.01}^{+0.01}$
${\theta }_q$ (∘)	${32.7}_{-0.9}^{+0.8}$	${32.8}_{-0.9}^{+0.8}$	${34.7}_{-1.0}^{+1.1}$	${31.9}_{-0.7}^{+0.8}$	${32.9}_{-0.7}^{+0.6}$	${33.5}_{-0.7}^{+0.8}$	${31.5}_{-0.7}^{+0.7}$
${\gamma }^{\prime }$	${1.96}_{-0.02}^{+0.02}$	${1.98}_{-0.02}^{+0.02}$	${2.01}_{-0.02}^{+0.01}$	${1.90}_{-0.01}^{+0.02}$	${1.94}_{-0.01}^{+0.01}$	${1.95}_{-0.01}^{+0.01}$	${2.02}_{-0.01}^{+0.01}$
${\gamma }_{\mathrm{ext}}$	${0.117}_{-0.004}^{+0.004}$	${0.120}_{-0.004}^{+0.003}$	${0.125}_{-0.004}^{+0.004}$	${0.109}_{-0.003}^{+0.003}$	${0.110}_{-0.004}^{+0.003}$	${0.110}_{-0.003}^{+0.003}$	${0.126}_{-0.003}^{+0.003}$
${\theta }_{\gamma }$ (∘)	$-{89.6}_{-1.3}^{+0.9}$	$-{89.1}_{-1.3}^{+0.8}$	$-{87.7}_{-0.6}^{+0.5}$	${88.0}_{-0.6}^{+0.7}$	${82.4}_{-1.0}^{+0.9}$	${82.6}_{-1.9}^{+1.8}$	$-{88.4}_{-0.5}^{+0.5}$
X ${\theta }_E\left(″\right)$	${0.001}_{-0.001}^{+0.001}$	${0.001}_{-0.001}^{+0.001}$	${0.001}_{-0.001}^{+0.001}$	${0.001}_{-0.001}^{+0.001}$	${0.001}_{-0.001}^{+0.001}$	${0.001}_{-0.001}^{+0.001}$	${0.001}_{-0.000}^{+0.001}$
G2 ${\theta }_E\left(″\right)$	${0.926}_{-0.027}^{+0.029}$	${0.939}_{-0.025}^{+0.023}$	${0.930}_{-0.036}^{+0.036}$	${0.945}_{-0.027}^{+0.021}$	${0.868}_{-0.008}^{+0.008}$	${0.947}_{-0.014}^{+0.017}$	${0.929}_{-0.014}^{+0.015}$
G2 $q$	${0.66}_{-0.02}^{+0.02}$	${0.67}_{-0.02}^{+0.02}$	${0.67}_{-0.03}^{+0.02}$	${0.68}_{-0.01}^{+0.02}$	${0.70}_{-0.02}^{+0.01}$	${0.65}_{-0.02}^{+0.01}$	${0.69}_{-0.02}^{+0.01}$
G2 ${\theta }_q$ (∘)	${41.5}_{-4.7}^{+4.0}$	${39.8}_{-3.4}^{+4.1}$	${38.2}_{-5.3}^{+6.0}$	${35.2}_{-3.4}^{+2.9}$	${40.4}_{-1.9}^{+2.6}$	${35.4}_{-2.1}^{+2.4}$	${46.2}_{-4.0}^{+3.7}$

Table 5. Table 5: Composite Model Parameters

Parameter	Marginalized Constraints
	Fiducial	AGNmask+1	AGNwht=0	Arc+1,60src	Group	Group + $z$=0.49 group
Stellar M/L ($M_{\odot }/L_{\odot }$)141414M/L within ${\theta }_E$ for rest-frame $V$ band. The point source component of the lens light is assumed to be from low-level AGN emission as opposed to stellar light and is not included in the calculation. The given uncertainties are a combination of statistical effects and a systematic uncertainty equal to the difference between the calculated M/L with and without the point source contribution. The stellar mass is calculated assuming $H_0=70\mathrm{km}{s}^{-1}{\mathrm{Mpc}}^{-1}$, ${\Omega }_m=0.3$, ${\Omega }_{\Lambda }=0.7$, but changes in the cosmology affect the M/L by a negligible amount.	${2.2}_{-0.2}^{+0.2}$	${2.1}_{-0.2}^{+0.2}$	${2.1}_{-0.2}^{+0.2}$	${2.1}_{-0.2}^{+0.2}$	${2.1}_{-0.2}^{+0.2}$	${2.1}_{-0.2}^{+0.2}$
Chameleon1 $q$	${0.759}_{-0.001}^{+0.003}$	${0.761}_{-0.001}^{+0.001}$	${0.761}_{-0.001}^{+0.001}$	${0.758}_{-0.001}^{+0.001}$	${0.763}_{-0.001}^{+0.001}$	${0.765}_{-0.000}^{+0.000}$
Chameleon1 ${\theta }_q$ (∘)	${23.1}_{-0.2}^{+0.2}$	${23.1}_{-0.2}^{+0.2}$	${23.3}_{-0.1}^{+0.2}$	${23.3}_{-0.2}^{+0.1}$	${23.3}_{-0.1}^{+0.1}$	${23.4}_{-0.1}^{+0.1}$
Chameleon2 $q$	${0.771}_{-0.001}^{+0.001}$	${0.772}_{-0.001}^{+0.001}$	${0.774}_{-0.001}^{+0.001}$	${0.770}_{-0.001}^{+0.001}$	${0.770}_{-0.001}^{+0.001}$	${0.767}_{-0.001}^{+0.001}$
Chameleon2 ${\theta }_q$ (∘)	${26.1}_{-0.3}^{+0.3}$	${26.4}_{-0.3}^{+0.3}$	${25.5}_{-0.3}^{+0.3}$	${25.8}_{-0.3}^{+0.3}$	${26.4}_{-0.2}^{+0.2}$	${26.3}_{-0.2}^{+0.2}$
NFW ${\kappa }_{0,h}$	${0.143}_{-0.007}^{+0.003}$	${0.148}_{-0.003}^{+0.003}$	${0.147}_{-0.004}^{+0.003}$	${0.159}_{-0.004}^{+0.005}$	${0.136}_{-0.003}^{+0.002}$	${0.149}_{-0.001}^{+0.001}$
NFW $r_s\left(″\right)$	${10.46}_{-0.14}^{+0.07}$	${10.31}_{-0.08}^{+0.12}$	${10.53}_{-0.08}^{+0.08}$	${10.32}_{-0.08}^{+0.07}$	${11.35}_{-0.10}^{+0.05}$	${11.10}_{-0.05}^{+0.05}$
NFW $q$	${0.88}_{-0.01}^{+0.02}$	${0.90}_{-0.00}^{+0.00}$	${0.92}_{-0.01}^{+0.01}$	${0.88}_{-0.01}^{+0.01}$	${0.88}_{-0.01}^{+0.01}$	${0.89}_{-0.02}^{+0.01}$
NFW ${\theta }_q$ (∘)	${71.6}_{-1.0}^{+1.1}$	${73.5}_{-0.4}^{+0.4}$	${73.6}_{-0.7}^{+0.5}$	${71.3}_{-0.8}^{+0.7}$	${71.7}_{-0.4}^{+0.3}$	${73.4}_{-0.4}^{+0.4}$
${\gamma }_{\mathrm{ext}}$	${0.138}_{-0.002}^{+0.001}$	${0.137}_{-0.001}^{+0.001}$	${0.135}_{-0.001}^{+0.001}$	${0.137}_{-0.001}^{+0.002}$	${0.134}_{-0.002}^{+0.002}$	${0.133}_{-0.001}^{+0.002}$
${\theta }_{\gamma }$ (∘)	$-{89.1}_{-0.3}^{+0.3}$	$-{88.6}_{-0.2}^{+0.1}$	$-{88.1}_{-0.1}^{+0.1}$	${89.9}_{-0.2}^{+0.3}$	${89.5}_{-0.4}^{+0.3}$	${89.2}_{-0.4}^{+0.5}$
X ${\theta }_E\left(″\right)$	${0.018}_{-0.002}^{+0.002}$	${0.015}_{-0.002}^{+0.002}$	${0.010}_{-0.002}^{+0.002}$	${0.021}_{-0.002}^{+0.002}$	${0.018}_{-0.002}^{+0.003}$	${0.018}_{-0.002}^{+0.002}$
G2 ${\theta }_E\left(″\right)$	${1.034}_{-0.005}^{+0.004}$	${1.036}_{-0.003}^{+0.004}$	${1.046}_{-0.005}^{+0.005}$	${1.034}_{-0.004}^{+0.004}$	${1.026}_{-0.002}^{+0.002}$	${1.005}_{-0.003}^{+0.002}$
G2 $q$	${0.94}_{-0.01}^{+0.01}$	${0.91}_{-0.01}^{+0.01}$	${0.90}_{-0.01}^{+0.01}$	${0.95}_{-0.01}^{+0.01}$	${0.93}_{-0.01}^{+0.01}$	${0.93}_{-0.00}^{+0.00}$
G2 ${\theta }_q$ (∘)	${39.7}_{-2.1}^{+1.3}$	${40.5}_{-0.7}^{+0.6}$	${40.3}_{-0.8}^{+0.8}$	${40.5}_{-0.8}^{+1.0}$	${29.4}_{-0.8}^{+1.0}$	${40.5}_{-1.0}^{+0.8}$

Table 6. Table 6: Effective time-delay distance and BIC weighting for different source resolutions

Model	$D_{\Delta t}$ (Mpc)	$\Delta$BIC	Relative BIC weight
SPEMD fiducial, $47\times 47$ source	${4711}_{-184}^{+232}$	92	0.316
SPEMD fiducial, $48\times 48$ source	${4706}_{-180}^{+227}$	76	0.470
SPEMD fiducial, $49\times 49$ source	${4685}_{-210}^{+241}$	139	0.049
SPEMD fiducial, $50\times 50$ source	${4640}_{-195}^{+238}$	86	0.365
SPEMD fiducial, $51\times 51$ source	${4631}_{-184}^{+229}$	62	0.605
SPEMD fiducial, $52\times 52$ source	${4712}_{-193}^{+244}$	94	0.295
SPEMD fiducial, $53\times 53$ source	${4691}_{-202}^{+243}$	59	0.641
SPEMD fiducial, $54\times 54$ source	${4645}_{-172}^{+230}$	32	0.864
SPEMD fiducial, $56\times 56$ source	${4686}_{-186}^{+231}$	0	1.000
SPEMD fiducial, $58\times 58$ source	${4724}_{-188}^{+238}$	14	0.956
SPEMD fiducial, $60\times 60$ source	${4844}_{-183}^{+243}$	12	0.964
Composite fiducial, $47\times 47$ source	${4933}_{-253}^{+401}$	0	1.000
Composite fiducial, $48\times 48$ source	${4753}_{-251}^{+418}$	10	0.507
Composite fiducial, $49\times 49$ source	${5011}_{-259}^{+420}$	114	0.000
Composite fiducial, $50\times 50$ source	${4732}_{-247}^{+417}$	26	0.175
Composite fiducial, $51\times 51$ source	${4810}_{-254}^{+417}$	118	0.000
Composite fiducial, $52\times 52$ source	${4956}_{-261}^{+422}$	60	0.013
Composite fiducial, $53\times 53$ source	${4893}_{-258}^{+431}$	69	0.006
Composite fiducial, $54\times 54$ source	${4819}_{-258}^{+422}$	114	0.000
Composite fiducial, $56\times 56$ source	${4886}_{-251}^{+423}$	112	0.000
Composite fiducial, $58\times 58$ source	${4885}_{-250}^{+424}$	129	0.000
Composite fiducial, $60\times 60$ source	${4843}_{-248}^{+415}$	189	0.000

Equations

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

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

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

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

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

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

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

$$I({\theta }_1,{\theta }_2)=A\exp \left[-k\left({\left(\frac{\sqrt{{\theta }_1^2+{\theta }_2^2/q_{\mathrm{L}}^2}}{r_{\mathrm{eff}}}\right)}^{1/n}-1\right)\right],$$

$$\frac{{\sigma }_{\theta }^2}{{\sigma }_r^2}=1-\frac{r^2}{r_{\mathrm{ani}}^2+r^2},$$

$$P({\sigma }_{\mathrm{LOS}}|\mathbf{\nu },\mathbf{\pi },{\kappa }_{\mathrm{ext}},r_{\mathrm{ani}})$$

$$\mathrm{BIC}=\ln (n)k-2\ln (\hat{L}).$$

$$f_{\mathrm{BIC}}(x)=h(x,{\sigma }_{\mathrm{BIC}})\ast \exp \left(-\frac{x-\mathrm{BI}{\mathrm{C}}_{\min }}{2}\right),$$

$$P({\kappa }_{\mathrm{ext}}|{\mathbf{d}}_{\mathbf{LOS}})\equiv P({\kappa }_{\mathrm{ext}}|{\zeta }_q,...)$$

$$=\int \prod _{q}d{\zeta }_q{P}_{\mathrm{MS}}({\kappa }_{\mathrm{ext}}|{\zeta }_q^{\mathrm{MS}}\equiv {\zeta }_q,...)P({\zeta }_q,...|{\mathbf{d}}_{\mathbf{LOS}})$$

$$D_{\Delta t}^{ij}\equiv (1+z_i)\frac{D_i{D}_j}{D_{ij}},$$

$$t=\sum _{i=1}^{s-1}\frac{D_{\Delta t}^{i,i+1}}{c}\left[\frac{({\mathbf{\theta }}_{i+1}-{\mathbf{\theta }}_i{)}^2}{2}-{\beta }_{i,i+1}{\psi }_i({\mathbf{\theta }}_i)\right],$$

$${\beta }_{ij}=\frac{D_{ij}{D}_s}{D_j{D}_{is}}.$$

Full text

††pubyear: 2019††pagerange: H0LiCOW XII. Lens mass model of WFI2033$-$4723 and blind measurement of its time-delay distance and $H_{0}$–17

1 Introduction

The flat $\Lambda$CDM cosmological model, characterized by spatial flatness, dark energy in the form of a cosmological constant, and cold dark matter, is considered to be the standard cosmological model today. Although this model is known as the concordance model, sources of tension have nonetheless begun to appear as the results of different cosmological experiments have grown in precision. Most notably, the tension between the measurements of the Hubble constant from the analysis of the cosmic microwave background (CMB) by the Planck mission (under the strict assumption of flat $\Lambda$CDM) and of Type Ia supernovae standard candles calibrated using the local distance ladder by the Supernovae, H0, for the Equation of State of Dark Energy collaboration (SH0ES; Riess et al., 2016) has recently increased from $3.4\sigma$ (Planck Collaboration et al., 2015; Riess et al., 2016) to $4.4\sigma$ (Planck Collaboration et al., 2018; Riess et al., 2019). The latest results are $H_{0}=67.4\pm 0.5~{}\mathrm{km~{}s^{-1}~{}Mpc^{-1}}$ from Planck, and $H_{0}=74.03\pm 1.42~{}\mathrm{km~{}s^{-1}~{}Mpc^{-1}}$ from SH0ES.

At present, sources of systematic error in either of these measurements that are significant enough to explain the discrepancy have not been demonstrated to exist. This opens up the intriguing possibility of having to extend the standard cosmological model by allowing for curvature, more general dark energy, or increasing the number of neutrinos (see Figure 13 in Riess et al., 2019), or to consider exotic alternatives, such as a vacuum phase transition (Di Valentino et al., 2018), early dark energy models (Poulin et al., 2018), self-interacting neutrinos (Kreisch et al., 2019) or decaying dark matter (Vattis et al., 2019). The various parameters of such extensions are highly degenerate with the value of $H_{0}$, and therefore a high-precision determination, with a technique independent of, and therefore not subject to the same systematics of either Planck or SH0ES, is in demand (e.g., Hu, 2005; Suyu et al., 2012a; Weinberg et al., 2013). Some proposed independent methods, such as water masers (e.g., Gao et al., 2016; Braatz et al., 2018), extragalactic background light attenuation (e.g., Domínguez et al., 2019), and gravitational waves (e.g., Feeney et al., 2019), etc. have yet to resolve the $H_{0}$ discrepancy, as their precision is not yet comparable to Planck or SH0ES.

The time-delay cosmography technique uses gravitational lens time delays to measure $H_{0}$. This technique rests on the fact that light rays from a multiply-imaged source will take different paths as they propagate through spacetime, with different geometrical lengths and gravitational potential depths. This will introduce an offset in arrival times, which can be measured through monitoring, if the source brightness varies in time. The measured time delays are used to infer the “time delay distance", primarily sensitive to $H_{0}$, which therefore provides a one-step way of measuring $H_{0}$ (e.g., Vanderriest et al., 1989; Keeton & Kochanek, 1997; Oguri, 2007; Suyu et al., 2010). Although proposed more than half a century ago by Refsdal (1964) in the context of lensed supernovae, the original idea has only recently been implemented (Grillo et al., 2018). Far more common is the use of gravitationally lensed quasars, given the sample of 250 such systems known to date (e.g., Lemon et al., 2019).

In practice, an accurate measurement of $H_{0}$ through this method requires extensive observational data for each system, as well as the development of advanced modeling techniques (Suyu & Halkola, 2010; Suyu et al., 2012b; Tewes et al., 2013; Birrer et al., 2015; Bonvin et al., 2016; Birrer & Amara, 2018), and has only become feasible in the current decade. Our collaboration, $H_{0}$ Lenses in COSMOGRAIL’s Wellspring (H0LiCOW; Suyu et al., 2017, hereafter H0LiCOW I) is designed to perform such measurements. We have precise, long-term time-delay measurements from the COSmological MOnitoring of GRAvItational Lenses (COSMOGRAIL; Courbin et al., 2005; Eigenbrod et al., 2005; Bonvin et al., 2018) project. We use deep high-resolution imaging from the Hubble Space Telescope (HST) or adaptive optics that provide constraints on the lens model not only from the point-source positions, but also from the extended arcs of the lensed quasar host galaxy. Finally, we have velocity dispersion measurements of the lens galaxies and characterizations of their environments and line of sight (LOS; e.g., Collett et al., 2013; Greene et al., 2013; McCully et al., 2014, 2017; Tihhonova et al., 2018), in order to reduce the mass-sheet degeneracy (e.g., Falco et al., 1985; Schneider & Sluse, 2013).

With four lenses, we measured $H_{0}=72.5_{-2.3}^{+2.1}~{}\mathrm{km~{}s^{-1}~{}Mpc^{-1}}$ with a precision of 3.0% (Birrer et al., 2019, hereafter H0LiCOW IX) including systematic uncertainties, achieving our previous goal of the program of reaching $<3.5\%$ precision from the five separate lenses in the base H0LiCOW sample (see H0LiCOW I) and finding good agreement with SH0ES. We have thus shown that we are on track to measure $H_{0}$ with a precision of 1% from a future sample of $\sim 40$ lenses with comparable precision per system (e.g., Treu & Marshall, 2016; Shajib et al., 2018), a result which will have significant implications for understanding the current tension with the CMB value. Time-delay cosmography is therefore a very effective technique, in the sense that only a relatively small number of systems is required to achieve a tight precision. The efficiency is similar to that expected for gravitational wave detections with optical counterparts (Chen et al., 2018b). As we work towards the 1% precision goal from a sample of lenses, it is important to keep systematics in the inference of $H_{0}$ from individual systems within the 1% threshold, in order to insure accuracy, and also to test for biases by using multiple codes (Birrer et al., 2019) and data challenges (Liao et al., 2015; Ding et al., 2018).

In this paper, we present the results of a detailed lens modeling analysis of the gravitational lens WFI2033$-$4723 (J2000: 20^{\mathrm{h}}33^{\mathrm{m}}41\mbox{.!!^{\mathrm{s}}}9, $-47^{\circ}23\arcmin 43\aas@@fstack{\prime\prime}4$), a quadruply-lensed quasar discovered by Morgan et al. (2004). The source redshift is $z_{\rm s}=1.662$ (Sluse et al., 2012), and the main deflector is a massive elliptical galaxy at a redshift of $z_{\rm d}=0.6575\pm 0.0002$ (Sluse et al., 2019, hereafter H0LiCOW X), updating the $z_{\rm d}=0.661\pm 0.001$ measurement from Eigenbrod et al. (2006)). Bonvin et al. (2019b) (hereafter COSMOGRAIL XVIII) measure the time delays between the quasar images based on 14 years of monitoring, and H0LiCOW X study the environment and LOS to the lens, based on multi-band imaging and targeted spectroscopy. Our work supersedes the models presented in Vuissoz et al. (2008) (hereafter COSMOGRAIL VII), which are based on monitoring of shorter duration and constrained only by the positions of the quasar images.

This is the fifth H0LiCOW system analyzed in this manner, following B1608$+$656 (Suyu et al., 2010), RXJ1131$-$1231 (Suyu et al., 2013, 2014), HE 0435$-$1223 (Wong et al., 2017, hereafter H0LiCOW IV), and SDSS 1206$+$4332 (H0LiCOW IX), with a sixth lens, PG 1115$+$080, analyzed simultaneously (Chen et al., 2019). A H0LiCOW milestone paper (Wong et al., 2019) presents the results of a conjoined cosmographical analysis of these lenses.

This paper is organized as follows. We give a brief overview of using time-delay lenses for cosmography in Section 2. In Section 3, we describe the observational data used in our analysis. We describe our lens modeling procedure in Section 4. In Section 5 we quantify the effect of the lens environment in terms of an external convergence. The time-delay distance results and their implications for cosmology are presented in Section 6. We summarize our main conclusions in Section 7.

2 Summary of Time-Delay Cosmography

2.1 Time-delay distance

When a source is gravitationally lensed by a foreground mass, the arrival time of photons traveling from the source to the observer depends on both the path length and the gravitational potential traversed by the light rays. For a single lens plane, the excess time delay of an image at a position $\bm{\theta}=(\theta_{1},\theta_{2})$ on the sky with a corresponding source position $\bm{\beta}=(\beta_{1},\beta_{2})$ relative to the case of no lensing is

[TABLE]

where $D_{\Delta t}$ is the time-delay distance and $\psi(\bm{\theta})$ is the lens potential. The time-delay distance $D_{\Delta t}$ (Refsdal, 1964; Schneider et al., 1992; Suyu et al., 2010) is defined111For historical reasons, the time-delay distance is written in terms of angular diameter distances. A more natural definition is $D_{\Delta t}\equiv{\hat{D}_{\mathrm{d}}\hat{D}_{\mathrm{s}}}/{\hat{D}_{\mathrm{ds}}}$ where $\hat{D}$ are the proper distances that the photons have travelled. as

[TABLE]

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

The time delay between two images, $i$ and $j$, of a lensed source is the difference of their excess time delays,

[TABLE]

where $\bm{\theta}_{i}$ and $\bm{\theta}_{j}$ are the positions of images $i$ and $j$, respectively, in the image plane. If the source is variable on short timescales (on the order of weeks to months), it is possible to monitor the lensed image fluxes at positions $\bm{\theta}_{i}$ and $\bm{\theta}_{j}$ and measure the time delay, $\Delta t_{ij}$, between them (e.g., Vanderriest et al., 1989; Schechter et al., 1997; Fassnacht et al., 1999; Fassnacht et al., 2002; Kochanek et al., 2006; Courbin et al., 2011). The lens potentials at the image positions, $\psi(\bm{\theta}_{i})$ and $\psi(\bm{\theta}_{j})$, as well as at the source position, $\bm{\beta}$, can be determined from by modeling the system. In this way, lenses with measured time delays and accurate lens models can constrain $D_{\Delta t}$, and in turn, $H_{0}$.

If there are multiple deflectors at different redshifts, the observed time delays depend on combinations of the angular diameter distances among the observer, the multiple deflectors, and the source. The observed image positions are determined by the multi-plane lens equation (e.g., Blandford & Narayan, 1986; Kovner, 1987; Schneider et al., 1992; Petters et al., 2001; Collett & Auger, 2014; McCully et al., 2014), but there is no longer a unique time-delay distance associated with the system. However, if the lensing is dominated by the mass in a single redshift plane, the observed time delays are mostly sensitive to the time-delay distance (Equation 2), with the deflector redshift set to the redshift of the main lens plane. This approximation is valid for WFI2033$-$4723 (see Appendix B), and thus our results can be interpreted in terms of the “effective" time-delay distance, $D_{\Delta t}(z_{\rm d},z_{\rm s})$. Hereafter, $D_{\Delta t}$ refers to the effective time-delay distance unless otherwise indicated.

A complicating factor in determining the time delay is the “microlensing time delay", an effect first described by Tie & Kochanek (2018). Stars and compact objects in the lens galaxy can act as microlenses, which causes a differential magnification of the accretion disk of the lensed quasar. Since the microlensing effect is different at the positions of the various lensed images and varies over time as the microlenses move, this may create an additional bias and scatter in the measurement of the time delay between different images. The microlensing time delay depends on a number of assumptions about the accretion disk size, its orientation and inclination, and the propagation of radiation through the disk. The effect tends to be small, of order $\sim$days or shorter, and can be modeled and accounted for under proper assumptions (Bonvin et al., 2018). This effect can also be mitigated by using the relative offsets between the measured time delays and those expected from lens modeling (Chen et al., 2018a).

Another difficulty is due to the fact that external perturbations from mass along the LOS can affect to the lens potential that light rays pass through. These perturbations not only can affect the lens model of the system, but also lead to additional focusing and defocusing of the light rays, which also affect the measured time delays (e.g., Seljak, 1994). If unaccounted for, these perturbers can lead to biased inferences of $D_{\Delta t}$. If the effects of LOS perturbers are small enough that higher-order terms are unimportant (Keeton, 2003; McCully et al., 2014), they can be approximated by an external convergence term in the lens plane. The true $D_{\Delta t}$ is related to the $D_{\Delta t}^{\rm model}$ inferred from a mass model by

[TABLE]

Here, $\kappa_{\rm ext}$ cannot, in general, be constrained from the lens model due to the mass-sheet degeneracy (e.g., Falco et al., 1985; Gorenstein et al., 1988; Saha, 2000), in which the addition of a uniform mass sheet and a rescaling of the source plane coordinates can affect the inferred $D_{\Delta t}$ but leaves other observables unchanged.

This degeneracy can be substantially mitigated by estimating the mass distribution along the LOS (e.g., Fassnacht et al., 2006; Momcheva et al., 2006, 2015; Williams et al., 2006; Wong et al., 2011) and assuming that the physical mass of the deflector profile goes to zero at large radius. However, perturbers that are very massive or projected very close to the lens may need to be included explicitly in the mass model since their higher-order effects need to be accounted for (McCully et al., 2017). In contrast, the lens profile is also degenerate with the time-delay distance in that the radial profile slope is tightly correlated with the time-delay distance (e.g., Kochanek, 2002; Wucknitz, 2002; Suyu, 2012). This degeneracy can affect models with the same form of mass density profile (e.g., a power-law density profile), as well as models with different forms of density profiles (described analytically or not). Furthermore, this degeneracy can mimic the effects of the mass-sheet degeneracy because different profiles can approximate or exactly match mass-sheet transformations of one form or another (e.g., Schneider & Sluse, 2013, 2014; Unruh et al., 2017). These degeneracies can be reduced by combining the lensing data with stellar kinematics information (e.g., Treu & Koopmans, 2002; Koopmans et al., 2003; Auger et al., 2010; Suyu et al., 2014; Yıldırım et al., 2019), and by making reasonable assumptions about the mass profile. Including a velocity dispersion measurement in the modeling helps constrain any internal uniform mass component from a local galaxy group that the dynamics is sensitive to (Koopmans, 2004).

2.2 Joint Inference

Our inference of $D_{\Delta t}$ generally follows that of previous H0LiCOW analyses (Suyu et al., 2013, H0LiCOW IV, IX). Our observational data are denoted by $\bm{d_{\mathrm{HST}}}$ for the HST imaging data, $\bm{\Delta t}$ for the time delays, $\bm{\sigma}$ for the velocity dispersion of the lens galaxy, and $\bm{d_{\mathrm{LOS}}}$ for the LOS mass distribution determined from our photometric and spectroscopic data. We want to determine the posterior probability distribution function (PDF) of the model parameters $\bm{\xi}$ given the data, $P(\bm{\xi}|\bm{d_{\rm HST},\Delta t,\sigma,d_{\mathrm{LOS}},A})$. The vector $\bm{\xi}$ includes the lens model parameters $\bm{\nu}$, the cosmological parameters $\bm{\pi}$, and nuisance parameters representing the external convergence ($\kappa_{\rm ext}$; Section 5) and anisotropy radius for the lens stellar velocity ellipsoid ($r_{\rm ani}$; Section 4.3). ${\bm{A}}$ denotes a discrete set of assumptions about the form of the model, which includes the data modeling region, the source reconstruction grid, the treatment of the various deflector mass distributions, etc. In general, ${\bm{A}}$ is not fully captured by continuous parameters. From Bayes’ theorem, we have

[TABLE]

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

[TABLE]

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

For each of our main lens models in Section 4.2.1 and Section 4.2.2, we have a range of systematics tests (Section 4.2.3) where we vary the content of $\bm{A}$ and repeat the inference of ${\bm{\xi}}$. These tests are important for checking the magnitude of various known but unmodeled systematic effects, but leave us with the question of how to combine the results. We follow H0LiCOW IX in using the Bayesian Information Criterion (BIC) to weight the various models in our final inference (Section 4.4). This effectively combines our various assumptions $\bm{A}$ using the BIC so that we obtain $P(\bm{\xi}|\bm{d_{\rm HST},\Delta t,\sigma,d_{\mathrm{LOS}}})$. We can further marginalise over the non-cosmological parameters ($\bm{\nu}$, $\kappa_{\rm ext}$, $r_{\rm ani}$) and obtain the posterior probability distribution of the cosmological parameters $\bm{\pi}$:

[TABLE]

In the lens model, we actually vary $H_{0}$, keeping other parameters fixed at $w=-1$, $\Omega_{\mathrm{m}}=0.3$, and $\Omega_{\Lambda}=0.7$. This assumes a fixed curvature of the expansion history of the Universe, but not the absolute scale (represented by $H_{0}$ or $D_{\Delta t}$). This is done because there is not a unique $D_{\Delta t}$ when accounting for multiple lens planes, but we convert this to an “effective" $D_{\Delta t}$ that is insensitive to assumptions of the cosmological model (see Appendix B). Specifically, given the lens/quasar redshifts and $\bm{\pi}$ (i.e., $H_{0}$ and the other fixed cosmological parameters), we can compute the effective time-delay distance $D_{\Delta t}(\bm{\pi},z_{\rm d},z_{\rm s})$ to obtain the posterior probability distribution of $D_{\Delta t}$, $P(D_{\Delta t}|\bm{d_{\rm HST},\Delta t,\sigma,d_{\mathrm{LOS}}})$.

3 Data

The data we use to infer $D_{\Delta t}$ consists of 1) the HST imaging used for lens modeling, which we present in Section 3.1; 2) the spectroscopy of the lensing galaxy, used to measure its stellar velocity dispersion, and 3) targeted spectroscopy of the LOS environment, both of which we present in Section 3.2; 4) wide-field multi-band imaging, which we present in Section 3.3 and we use to infer $\kappa_{\mathrm{ext}}$ in Section 5; and 5) the time delays measured by COSMOGRAIL, presented in Section 3.4.

3.1 HST Imaging

The HST images we use to model WFI2033$-$4723 consist of Wide Field Camera 3 (WFC3) F160W band observations (Program #12889; PI: Suyu), as well as archival Advanced Camera for Surveys (ACS) observations in the F814W filter (Program #9744; PI: Kochanek). The latter program also contains imaging in the F555W filter, which we do not use, because the signal-to-noise ratio from the lensed images is low and does not add much information to the lens model.

The details of the observations from Program #12889 are presented in H0LiCOW I. Using a combination of short (74 s) and long (599–699 s) exposures, we obtain the brightness distribution of the lens system covering a large dynamic range (of the bright lensed AGN, its much fainter host galaxy, and the foreground lens galaxy). The WFC3 images are drizzled using DrizzlePac222DrizzlePac is a product of the Space Telescope Science Institute, which is operated by AURA for NASA. to a final pixel scale of $0\aas@@fstack{\prime\prime}08$, whereas the ACS images are reduced using MultiDrizzle333MultiDrizzle is a product of the Space Telescope Science Institute, which is operated by AURA for NASA. onto a final pixel scale of 0$\aas@@fstack{\prime\prime}$05. More details of the reduction are presented in H0LiCOW IV.

We create cutouts of the reduced HST images and define an arcmask around the lens in each of the two filters, which encloses the region where we reconstruct the lensed arc from the extended quasar host galaxy. We expand the cutout to the west of the lens to include the nearby galaxy G2, which is a bright perturber at $z=0.7450$ whose light profile needs to be modeled, as it may contaminate the signal within the arcmask. The cutout region is $10\aas@@fstack{\prime\prime}4\times 6\aas@@fstack{\prime\prime}4$, which corresponds to a $208\times 128$ pixel cutout for the F814W image and a $130\times 80$ pixel cutout for the F160W image. These cutouts are shown in Figure 1.

The reconstruction of the point spread function (PSF) for each HST exposure, as well as of the weight images and bad pixel masking for each cutout, are analogous to the procedure described in H0LiCOW IV. As detailed in that paper, we note that in order to avoid biasing the modeling due to large residuals from PSF mismatch near the AGN image centers, we rescale the weights in those regions by a power law model such that a pixel originally given a noise value of $p_{i}$ is rescaled to a noise value of $A\times p_{i}^{b}$. The constants $A$ and $b$ are chosen for each band such that the normalized residuals in the AGN image regions are approximately consistent with the normalized residuals in the rest of the arc region. For completeness, we show the residuals for models using the weight images without this power-law weighting in Appendix A. The strong residuals in these images motivates our decision to adopt this rescaling.

We note that although the background noise for the WFC3 IR camera depends on the number of non-destructive reads, we check that the number of reads in the lensed arc region is the same as for the blank sky patch used for estimating the background noise, so this procedure is valid. Since most of the lens model constraints come from the parts of the lensing arcs away from the centers of the AGN images, we check that these arcs do not have pixels that were flagged as bad in too many exposures, which would otherwise affect our lens mass model.

3.2 Spectroscopic data

Our spectroscopic observations, presented in H0LiCOW X, reveal that the lens is part of a galaxy group at $z_{\mathrm{grp}}=0.6588$ with a velocity dispersion of $\sigma=500\pm 80~{}\mathrm{km~{}s}^{-1}$ measured from 22 member galaxies, which is independently confirmed by Wilson et al. (2016) based on a spectroscopic study by Momcheva et al. (2006, 2015).

We summarize hereafter the characteristics of the spectroscopic data used. A more exhaustive description of the data acquisition and analysis is provided in H0LiCOW X. WFI2033$-$4723 was observed with the ESO-MUSE integral field spectrograph (Bacon et al., 2010) during several observing runs between 2014-06-19 and 2016-07-20. The velocity dispersion measurement of the lensing galaxy was based on a total of 3$\times$2400 s exposures with the lensing galaxy located close to the center of the 1′$\times$1′ field of view (FOV). The data cubes are characterized by a $0\aas@@fstack{\prime\prime}2\times 0\aas@@fstack{\prime\prime}2$ spatial sampling, a wavelength coverage in the optical range from 4800 Å to 9350 Å, a spectral sampling of 1.25 Å per pixel, and a resolving power $R\sim 1800-3600$ (i.e. 2.5 Å spectral resolution; Richard et al., 2017). The analysis has been carried out on the combined datacube characterized by a median seeing of 1″. To deblend the lensing galaxy and the quasar images, we modeled each monochromatic slice with a model of the system composed of four Moffat (Moffat, 1969) components for the quasar lensed images, and one de Vaucouleurs (de Vaucouleurs, 1948) model for the lensing galaxy. After removing the quasar images from the datacube, we extracted the lensing galaxy spectrum within a square aperture of 9 pixels = 1$\aas@@fstack{\prime\prime}$8 side-length.

The velocity dispersion was obtained following the same procedure as Suyu et al. (2010, 2013), resulting in an inference of $\sigma_{\mathrm{LOS}}=250$ km s*-1* with a statistical uncertainty of $\approx 10~{}\mathrm{km~{}s^{-1}}$. The order of the polynomial continuum and spectral regions masked for the fit introduce additional systematic uncertainties. The various choices we made have been treated as nuisance parameters over which we have marginalised to derive our final velocity dispersion PDF (see H0LiCOW X). The overall uncertainty, accounting for the random and systematic errors, reaches $\sigma_{\sigma_{\mathrm{LOS}}}=19$ km s*-1*. We integrate this measurement in our cosmographic inference in Section 4.3.

In addition to ESO MUSE spectroscopy of the galaxies located in the vicinity of the lens, we have also obtained multi-object spectroscopy of the galaxies in the FOV with the ESO FORS (Appenzeller et al., 1998) and the Gemini GMOS (Hook et al., 2004) instruments. In total, we used 10 masks, with about 35 long-slits (6″ length) per mask positioned on targets located within 2′ from the lens. For each mask, we obtained 40 minutes long exposures, and used a setup allowing to cover most of the optical wavelength range (typically 4500-9000Å) with a resolving power of $\approx 440$ (FORS) / 1100 (GMOS).

3.3 Photometric data

Our photometric data consists of wide-field optical wavelength data from the Dark Energy Survey444https://www.darkenergysurvey.org (DES), ultraviolet data from the DES Camera (Flaugher et al., 2015) on the Blanco Telescope, VLT/HAWK-I (Pirard et al., 2004; Kissler-Patig et al., 2008) near-infrared data, and archival IRAC (Fazio et al., 2004) infrared data from the Spitzer Space Telescope. These data and their products, consisting of the galaxy-star classification, photometric redshifts and stellar masses of all galaxies with $i<23$ mag within a $120\arcsec$ radius around WFI2033$-$4723, are described in H0LiCOW X. In Section 5, where we measure the relative density of the environment of WFI2033$-$4723, we use a conservative cut of $i<22.5$ mag in order to ensure that the galaxy catalogue, with a $5\sigma$ limiting magnitude of $\sim 23.13$, is complete.555While shallow magnitude limits may bias the $\kappa_{\mathrm{ext}}$ distribution we determine in Section 5, Figure 6 in Collett et al. (2013) shows that the expected bias is at a level of $\sim 0.25\%$, which is acceptable given our goal of inferring $H_{0}$ with biases below the $1\%$ level. We show the $4\arcmin\times 4\arcmin$ FOV, with the galaxy catalogue overlapped, in Figure 2.

3.4 Time delays

3.4.1 Time-delay measurements

COSMOGRAIL XVIII presents the most comprehensive analysis of the time delays of WFI2033$-$4723 so far, with the analysis of four different data sets spanning across 14 years of monitoring, for a total of $\sim 447$ hours of observations. The data were acquired in the scope of the COSMOGRAIL collaboration, using three different telescopes in the Southern hemisphere; the C2 and ECAM instruments mounted on the 1.2m Leonhard Euler Swiss telescope and the WFI instrument mounted on the ESO/MPIA 2.2m telescope, both located at La Silla Observatory in Chile, and the 1.3m Small and Moderate Aperture Research System (SMARTS) at the Cerro Tololo Inter-American Observatory (CTIO) in Chile.

The data is split in four data sets, one per instrument (C2 and ECAM on the Euler telescope, WFI on the 2.2m telescope and SMARTS), each being reduced independently. The photometry of the four images of WFI2033$-$4723 is recovered using the MCS deconvolution scheme (Magain et al., 1998; Cantale et al., 2016). The light curves obtained are presented in Figure 2 of COSMOGRAIL Paper XVIII. For three of the four data sets (C2, ECAM and SMARTS), the deconvolution scheme is not able to properly resolve the flux coming from the A1 and A2 images. Thus, the A1 and A2 fluxes are summed into a virtual light curve A, under the assumption that the time delay between A1 and A2 is zero. The WFI data set being composed of exposures of better quality, the deconvolution scheme manages to properly resolve the A1 and A2 images. A virtual light curve A=A1+A2 is also constructed for WFI in order to compare it to the other data sets.

The time-delay measurements between each pair of light curves are made with the PyCS software (Tewes et al., 2013; Bonvin et al., 2016) and follow the formalism introduced by Bonvin et al. (2018). We use two different curve shifting techniques. Both techniques share a common framework to assess their own uncertainties, based on a statistical analysis of the residuals of the real data that prevents, by construction, involuntary fine-tuning of the curve-shifting technique parameters to recover a biased value of the time delays.

Each data set is analyzed independently. The time-delay estimates obtained are in good agreement with each other and a Bayes Factor analysis states that they can be combined without loss of consistency. In this work, we use the combined time-delay estimates with respect to image B. For our fiducial set of models, we use the B-A1 and B-A2 time delays estimated from the WFI data set (see Figure 4 of COSMOGRAIL XVIII), and the B-C time delay estimated by combining all the data sets together (labeled “PyCS-mult” on Figure 3 of COSMOGRAIL XVIII). They read $\Delta t_{\mathrm{B-A1}}=-36.2^{+1.6}_{-2.3}$, $\Delta t_{\mathrm{B-A2}}=-37.3^{+2.6}_{-3.0}$ and $\Delta t_{\mathrm{B-C}}=-59.4\pm 1.3$. Although using different time delays from different combinations of data sets might appear subjective, we recall that i) only the WFI data set is of good enough quality to resolve the A1 and A2 images, thus bringing an additional independent constraint to the modeling and solving the potential issue of where to anchor a time-delay estimate related to a virtual image A, and ii) all the time-delay estimates and combination of time-delay estimates are statistically consistent with each other.

3.4.2 Microlensing time-delay

Our time-delay measurements do not include the contribution from the microlensing time delay (Tie & Kochanek, 2018; Bonvin et al., 2019a), a time-dependent reweighting of the geometrical delay (originating from the extended spatial structure of the source) by the microlensing pattern affecting each image independently. As a result, an excess microlensing time delay adds to the excess cosmological time delay of each lensed image, and the measured time delays between pairs of images can deviate from the cosmological time delays by a noninegligible amount. The amplitude of the effect depends mainly on the mass of the central black hole of WFI2033$-$4723 (Sluse et al., 2012; Motta et al., 2017), and its estimation relies on the assumption that the accretion disk can be modeled as a thin-disk (Shakura & Sunyaev, 1973) - which, so far, is disfavored by the data (see e.g. Morgan et al., 2018) - and that the emission of the accretion disk follows an idealized lamp-post model (Cackett et al., 2007; Starkey et al., 2016).

In Figure 6 of COSMOGRAIL XVIII, we compute the amplitude of the microlensing time delay for various disk sizes. Although the measured time delays do not show any discrepancies that would be evidence for a microlensing time delay, it cannot be ruled out either. We thus chose to include it by default in our models, noting that the effect is much smaller than our other uncertainties. We follow the framework presented in Chen et al. (2018a) and assume the accretion disk size of Morgan et al. (2018) with $r=R_{0}$. We also test the effect of ignoring the microlensing time delay for one of our models, finding that it changes the $D_{\Delta t}$ accuracy by $<1\%$ (Section 4.2.3).

4 Lens Modeling

In this section, we describe our procedure to simultaneously model the images in the two HST bands, and the time delays, in order to infer the lens model parameters and $D_{\Delta t}$.

4.1 Overview

We perform our lens modeling using Glee, a software package developed by S. H. Suyu and A. Halkola (Suyu & Halkola, 2010; Suyu et al., 2012b). The lensing mass distribution is described by a parameterized profile. The extended host galaxy of the source is modeled separately on a $50\times 50$ pixel grid with curvature regularization (Suyu et al., 2006). The lensed quasar images are modeled as point sources on the image plane convolved with the PSF. The quasar image amplitudes are allowed to freely vary and are independent from the extended host galaxy light distribution to allow for variability due to microlensing, time delays, and substructure. The lens galaxy light distribution is modeled using either Sérsic profiles or Chameleon profiles. The Sérsic profile is defined as

[TABLE]

where $A$ is the amplitude, $k$ is a constant such that $r_{\rm eff}$ is the effective (half-light) radius, $q_{\mathrm{L}}$ is the axis ratio, and $n$ is the Sérsic index. The Chameleon profile (also known as the pseudo-Jaffe profile) is defined as the difference of two non-singular $r^{-2}$ elliptical profiles (Kassiola & Kovner, 1993; Dutton et al., 2011), which are a good approximation to Sérsic profiles.

We represent the galaxy light distribution as the sum of two Sérsic (or two Chameleon) profiles plus a point source (to account for possible AGN emission from the lens galaxy) with a common centroid. Since the light of G2 can also influence the model, we represent its light distribution as a single Sérsic profile plus a point source with a common centroid, although we mask its central regions since we only care about light from G2 that could affect the lens galaxy or arc light. There is a small nearby perturber (“X" in Figure 1), which we also represent as a single Sérsic profile plus a point source with a common centroid. Model parameters of the lens and source are constrained through Markov Chain Monte Carlo (MCMC) sampling.

In accounting for perturbers at different redshifts from the main lens galaxy, we use the full multi-plane lens equation (e.g., Blandford & Narayan, 1986; Kovner, 1987; Schneider et al., 1992; Petters et al., 2001; Collett & Auger, 2014; McCully et al., 2014) in our modeling. We vary $H_{0}$ directly in our models and use this distribution to calculate the effective model time-delay distance $D_{\Delta t}^{\rm model}$. In calculating $D_{\Delta t}^{\rm model}$, we assume $\Omega_{\rm m}=0.3$, $\Omega_{\Lambda}=0.7$, and $w=-1$. Relaxing these assumptions by allowing these cosmological parameters to vary freely shifts the resulting $D_{\Delta t}^{\rm model}$ distributions by $<1$% in previous analyses (see H0LiCOW IV), and we also verify that this is true for WFI2033$-$4723 (Appendix B). Thus, this approximation has no measurable effect on the inferred time delay distance, which can then be applied to constrain any arbitrary cosmology.

4.2 Mass Models

Our primary mass models for the lens galaxy are a singular power-law elliptical mass distribution (SPEMD; Barkana, 1998), and a model consisting of a baryonic component that traces the light distribution plus a separate dark matter component (hereafter the “composite" model). We also include an external shear in the strong lens plane. Non-linear couplings due to multi-plane effects are small and thus ignored.

We explicitly include the nearby perturber X in the lens model, linking its mass centroid to that of its light. Although we do not have a spectroscopic confirmation of the redshift of X, it is likely a satellite galaxy that is physically associated with the lens galaxy, given its small size and proximity. We also see evidence in the F160W image for possible tidal features emanating from X in the direction away from the lens galaxy, suggesting that it may be an infalling satellite. We therefore assume that X is at the same redshift as the lens and parameterize it as a singular isothermal sphere (SIS). In our models, X generally has a much smaller mass than the main lens, and therefore it has a minor influence on the potential, even if it is located at a different redshift.

We also explicitly include three nearby massive perturbing galaxies (denoted G2, G3, and G7, following the naming convention of Vuissoz et al., 2008) in Figure 3 that are projected close the lens. G2 is close enough that its influence may not be adequately described by external shear (H0LiCOW II; see also McCully et al., 2017), and H0LiCOW X showed that G3 ($z=0.6548$) and G7 ($z=0.6573$) may have a non-negligible higher-order influence on the model as well. Our updated estimation of the influence of these galaxies, computed in terms of the flexion shift considered in H0LiCOW X but taking into account the galaxy morphologies and velocity dispersions measured in that paper, shows that G2, with $\log M_{\star}\sim 11.15$, is in fact the only galaxy with significant impact on the modeling. Nonetheless, based on their proximity to the lensing galaxy, we choose to explicitly model G3 ($\log M_{\star}\sim 10.17$) and G7 ($\log M_{\star}\sim 11.16$) as well. G2 is modeled as a singular isothermal ellipsoid, which is a reasonable assumption since higher-order moments of the potential will have a small effect at the position of the main lens galaxy. G3 and G7 are modeled as SIS. The relative Einstein radii of G2, G3, and G7 are calculated from their measured velocity dispersions (H0LiCOW X), assuming isothermal profiles. The ratio of their Einstein radii is fixed, but with a global scaling allowed to vary freely, as in H0LiCOW IV, IX. This is done to prevent the model from optimizing the perturbers’ Einstein radii in a way that would be inconsistent with their measured redshifts and velocity dispersions. The centroid of G2 is linked to the centroid of its light distribution in the F160W band in the modeling, while the centroids of G3 and G7 are fixed to their measured positions in the F160W image. We set the redshifts of G3 and G7 equal to the lens redshift of $z=0.6575$ in our model, as their redshifts are consistent with this value within the range allowed by peculiar velocities. All masses are treated using the full multi-plane lens equation, as detailed by Suyu et al., in preparation.

Our constraints on the primary lens model include the positions of the lensed quasar images, the measured time delays, and the surface brightness of the pixels in the ACS/F814W and WFC3/F160W images that are fit simultaneously. The quasar positions are fixed to the positions of the point sources on the image plane (after they have stabilized) and are given a Gaussian uncertainty of width $0\aas@@fstack{\prime\prime}004$ to account for offsets due to substructure in the lens or LOS, which is small enough to satisfy astrometric requirements for cosmography (Birrer & Treu, 2019). The quasar flux ratios are not used as constraints, as they can be affected by microlensing. We first model the lens seperately in each band to iteratively update the respective PSFs using the lensed AGN images themselves, similar to Chen et al. (2016), but with the PSF corrections and source intensity reconstructed simultaneously in our case (H0LiCOW IV, IX) rather than separately. We keep these “corrected" PSFs fixed and use them in our final models that simultaneously use the surface brightness distribution in both bands as constraints. We then use the positions of the quasar images to align the images in the two HST bands. We do not enforce any similarity of pixel values at the same spatial position across different bands (i.e., the flux at any position in one band is independent of the other band). We also directly include the effect of microlensing time delays, as described in Section 3.4.2, although our tests show that this has a very small effect on our results (Section 4.2.3). In our MCMC sampling, we vary the light parameters of the lens galaxy, G2, X, and quasar images, the mass parameters of the lens galaxy, X, G2, G3, and G7, the external shear, and $H_{0}$. The quasar image positions are linked across both bands, but the other light parameters are allowed to vary independently.

4.2.1 Power-law Model

Our fiducial SPEMD model uses the double Sérsic parameterization for the lens galaxy light and has the additional free parameters:

1. (i)

position ($\theta_{1}$,$\theta_{2}$) of the centroid (allowed to vary independently from the centroid of the light distribution) 2. (ii)

Einstein radius $\theta_{\mathrm{E}}$ 3. (iii)

minor-to-major axis ratio $q$ and associated position angle $\theta_{q}$ 4. (iv)

3-dimensional slope of the power-law mass distribution $\gamma^{\prime}$ 5. (v)

position of X, linked to its light centroid 6. (vi)

Einstein radius of X 7. (vii)

position of G2, linked to its light centroid 8. (viii)

global scaling parameter that controls the Einstein radii of G2, G3, and G7 9. (ix)

minor-to-major axis ratio $q$ and associated position angle $\theta_{q}$ of G2 10. (x)

external shear $\gamma_{\rm ext}$ and associated position angle $\theta_{\gamma}$666$\theta_{\gamma}$ is defined to be the direction of the shear itself, i.e. orthogonal to the direction of the mass producing the shear. 11. (xi)

the Hubble constant, $H_{0}.$

We conservatively assume uniform priors on the model parameters over a wide physical range. Although the lens is not drawn from a random population, but rather with some selection function that could, in principle, bias the inferred time-delay distance, this selection function is not well known and these biases are negligible for this type of analysis (e.g., Collett & Cunnington, 2016). The parameters that are exceptions to our choice of uniform priors are that the global scaling parameter for the Einstein radii of the perturbers is given a Gaussian prior such that the expected mean and uncertainty of G2’s Einstein radius is constrained by its measured velocity dispersion, and that the position angle $\theta_{q}$ of G2’s is given a Gaussian prior based on the fit of its light profile. We anchor the scaling parameter to G2 as it is the perturber with the most precisely-measured velocity dispersion, and its proximity to the lens makes it the most significant of the three massive perturbing galaxies.

Figure 4 shows the data and the lens model results in both bands for our fiducial SPEMD model, as well as the source reconstructions. Our model reproduces the surface brightness structure of the lensed AGN and host galaxy in both bands simultaneously.

4.2.2 Composite Model

We follow Suyu et al. (2014) and H0LiCOW IV to construct the composite model, consisting of a baryonic component linked to the light profile of the lens galaxy, plus a dark matter component. The composite model assumes the double Chameleon light profile for the lens galaxy in the WFC3/F160W band scaled by an overall mass-to-light (M/L) ratio. We use the Chameleon light profiles for the composite model because it is straightforward to link the parameters describing the light distribution to those of the mass distribution, as they are fundamentally just a combination of isothermal profiles. We use the F160W band because it probes the rest-frame near-infrared and thus should be the best tracer of stellar mass. Although we include a point source in the light profile, we assume that this is due to low-level AGN emission from the lens galaxy, and do not associate it with a massive component in the model. This point source is roughly $\sim 2\%$ of the total light in the F160W band, so its inclusion would have a minor impact on our results. We keep the double Sérsic parameterization for the lens galaxy light in the F814W band to maintain consistency with the SPEMD models. The dark matter component is modeled as an elliptical NFW (Navarro et al., 1996) potential with the centroid linked to the light centroid in the F160W band, as non-contracted NFW profiles are a good representation of the dark matter halos of massive elliptical galaxies (Dutton & Treu, 2014).

Our fiducial composite model has the same free parameters (v) to (xi) as the SPEMD model in Section 4.2.1, as well as the additional parameters:

1. (a)

M/L ratio for the baryonic component 2. (b)

NFW halo normalization $\kappa_{\mathrm{0,h}}$ (defined as $\kappa_{\mathrm{0,h}}\equiv 4\kappa_{\mathrm{s}}$; Golse & Kneib, 2002) 3. (c)

NFW halo scale radius $r_{\mathrm{s}}$ 4. (d)

NFW halo minor-to-major axis ratio $q$ and associated position angle $\theta_{q}$

We set a Gaussian prior of $r_{\mathrm{s}}=11\aas@@fstack{\prime\prime}9\pm 1\aas@@fstack{\prime\prime}6$ based on the results of Gavazzi et al. (2007) for lenses in the Sloan Lens ACS Survey (SLACS; Bolton et al., 2006) sample, which encompasses the redshift and stellar mass of WFI2033$-$4723. All other parameters are given uniform priors, again with the exception of the Gaussian prior on global scaling parameter based on G2’s Einstein radius, as well as the Gaussian prior on G2’s position angle. The relative amplitudes of the two Chameleon profiles that represent the stellar light distribution of the lens galaxy can vary, but the relative amplitudes of these two components in the mass profiles are fixed. To account for this, we iteratively run a series of MCMC chains and update the relative amplitudes of the two mass components to match that of the light components after each chain. We iterate until the inferred $H_{0}$ stabilizes, then combine the chains after this point into a single distribution to represent the fiducial composite model. The other composite models use fixed relative amplitudes of the mass components based on the latest iteration of the fiducial composite model.

Figure 5 shows the data and the lens model results in both bands for the fiducial composite model described in this section, as well as the source reconstruction.

4.2.3 Systematics Tests

In this section we describe a range of tests of the effects of various systematics in our modeling, stemming from different assumptions in the way we constructed the model that might affect the posterior. In addition to the basic fiducial models described above, we perform inferences for both the SPEMD and composite models given the following sets of assumptions:

• •

A model with the arcmask region increased by one pixel on both the inner and outer edges. To compensate for the larger arcmask region, we increase the source plane resolution to $60\times 60$ pixels in all bands.

• •

A model where the region near the AGN images scaled by the power-law weighting is increased by one pixel around the outer edge. Increasing these regions by more pixels would start to greatly reduce the area of the arcmask where we fit the extended source.

• •

A model where the regions near the AGN images are given zero weight rather than being scaled by a power-law weighting.

• •

A model that includes the group at $z=0.6588$ (of which the lens galaxy is a member) as a spherical NFW halo. The halo centroid and mass are given Gaussian priors based on the calculations of H0LiCOW X. The scale radius is given a Gaussian prior of $r_{\mathrm{s,g}}=32\aas@@fstack{\prime\prime}0\pm 8\aas@@fstack{\prime\prime}0$ from a calculation of its virial mass and radius (H0LiCOW X) and a halo concentration based on the results of Diemer (2018). The redshift of the group is set to the lens redshift ($z=0.6575$), as the difference can be explained by peculiar velocity.

• •

A model that includes both the group at $z=0.6588$ (again set to the lens redshift) and a foreground group at $z=0.4956$ which may have a significant effect on the lens potential based on H0LiCOW X, who estimate its flexion shift (following the definition in McCully et al., 2014). The foreground group’s centroid, mass, and scale radius are given Gaussian priors in the same way as for the group at the lens redshift. The scale radius prior from H0LiCOW X and Diemer (2018) is $r_{\mathrm{s,gf}}=34\aas@@fstack{\prime\prime}8\pm 9\aas@@fstack{\prime\prime}3$.

In addition to the above models for both the SPEMD and composite models, we run one additional SPEMD model:

• •

A model where the light profile of the lens galaxy in both bands is represented by the sum of two Chameleon profiles rather than the sum of two Sérsic profiles.

As described in Section 4.4, we combine the MCMC chains from all of these tests, weighted by the BIC (e.g., H0LiCOW IX). We calculate the relative BIC for the SPEMD models and composite models separately, then give the combined distributions equal weight in the final inference so that we are not biased by the parameterization of the mass profile.

We also run a test to verify that the microlensing time delay does not significantly impact our results. We test our fiducial SPEMD model without including the microlensing time delay effect and compare the blinded effective time-delay distance to the model with this effect included, in Figure 6. We find that the microlensing time delay affects the inferred $D_{\Delta t}$ at $<1\%$, so its inclusion in our models, given our assumptions about the disk size, does not have an appreciable effect.

4.2.4 Comparison of Power Law and Composite Models

The marginalized parameter distributions of the SPEMD model are shown in Figure 7. We show the combined distributions of all SPEMD models where each model is given equal weight, as well as the BIC-weighted distribution. The parameter statistics for each model are given in Appendix C. There are some minor variations in the model parameters from model to model, but the $D_{\Delta t}$ distributions are generally consistent. We note that the model with Chameleon light profiles for the lens galaxy is somewhat offset toward a lower $D_{\Delta t}$ (see Section 6). This model is disfavored by our BIC weighting, so this has a minimal effect on our final results. This does not necessarily mean that the Chameleon profiles in general are a bad fit to the lens galaxy light, as the composite models (which use the Chameleon light profile by default in the F160W band) are not similarly offset.

The multi-modal distributions in some of the parameters arises primarily from differences in the posterior PDFs of different models corresponding to the various systematics tests, not from bimodality within individual lens models. We note that despite this multi-modal behavior, the effective $D_{\Delta t}$ distribution remains stable and unimodal, suggesting that the cosmological inference is robust to the various systematics tests.

The model that includes both group halos has the highest BIC weighting for both the SPEMD and composite models. To check that the addition of the $z=0.4956$ group contributes meaningful information to the modeling, we run a test where the centroid of this group is given a prior located at a similar distance but rotated by $90^{\circ}$ and $135^{\circ}$ on the sky relative to the lens. We compare the BIC weight values of these test models to that of the model with just the group at the lens redshift and the original model with both groups. These test cases show a lower BIC weight than the original model with both groups, suggesting that the addition of the foreground group with the actual centroid prior is contributing information, although the small BIC difference is within the typical BIC variance, so it is difficult to draw a firm conclusion. The $D_{\Delta t}$ distributions remain robust within the uncertainties for each of these test cases.

The offset between the mass centroid and the light centroid in the F160W band for the SPEMD model is typically $\sim 0\aas@@fstack{\prime\prime}02-0\aas@@fstack{\prime\prime}03$ (roughly $150-200$ pc for a flat $\Lambda$CDM cosmology with $h=0.7$ and $\Omega_{\rm m}=0.3$) such that the mass centroid is slightly southeast of the light centroid. This might be partially explained by the influence of object X, although we note that in our SPEMD models, the mass of X is consistent with zero. The centroids of the light profiles in F814W and F160W are consistent with each other at the $\sim 0\aas@@fstack{\prime\prime}002$ level for both models. The SPEMD models are able to fit the quasar positions to an rms of $\sim 0\aas@@fstack{\prime\prime}01$, while the composite models have a larger rms of $\sim 0\aas@@fstack{\prime\prime}025$. Despite these differences, the SPEMD and composite models’ $D_{\Delta t}$ distributions are not drastically different, and by weighting them equally in the final inference, we are accounting in part for the astrometric uncertainty.

We show the marginalized parameter distributions of the composite model in Figure 8. Again, we show the uniformly-combined distributions as well as the BIC-weighted composite model separately, and the parameter statistics for each model are given in Appendix C. As with the SPEMD model, there are small variations in the model parameters, but the $D_{\Delta t}$ inference is consistent.

We compare the physical parameters of our BIC-weighted SPEMD model to the composite model. The results are shown in Table 1, with the parameter statistics for all composite models given in Appendix C.

4.3 Kinematics

We compute the LOS stellar velocity dispersion of the strong lens galaxy through the spherical Jeans equation (see also Treu & Koopmans, 2002; Koopmans et al., 2003), similar to previous H0LiCOW analyses (e.g., Suyu et al., 2010, H0LiCOW IV). Yıldırım et al. (2019) recently showed that the assumption of spherical Jeans equation is applicable to time-delay cosmography with a single aperture-averaged lens velocity dispersion without significant bias, as in our case of WFI2033$-$4723. For a given lens model, we obtain the 3D mass profile of the lens galaxy by taking the spherical deprojection of the circularized surface mass density profile. The resulting 3D profile assumes analytical forms for both the SPEMD and the composite model. The 3D distribution of tracers is obtained by applying the same procedure to the surface brightness distribution of the lens galaxy, modeled as a Hernquist (1990) profile. We also tested a Jaffe (1983) profile, which has been shown to produce similar results (Suyu et al., 2010), and find that the results change by less than 1%. We parametrize the orbital anisotropy profile as an Osipkov-Merritt model (Osipkov, 1979; Merritt, 1985)

[TABLE]

where $\sigma_{\theta}$ and $\sigma_{r}$ are the tangential and radial velocity dispersions, respectively. Given values of the lens mass parameters in Section 4.2, the external convergence $\kappa_{\rm ext}$ in Section 5, and the anisotropy radius $r_{\rm ani}$, we then calculate the LOS velocity dispersion profile by numerically integrating the solutions of the spherical Jeans equation as given by Mamon & Łokas (2005). Finally, we calculate the integral over the spectroscopic slit of the seeing-convolved brightness-weighted LOS velocity dispersion $\sigma^{\rm P}$ (Equation (20) of Suyu et al., 2010) and compare to the measurements to calculate the likelihood of the kinematics data,

[TABLE]

where $\sigma_{\mathrm{LOS}}=250\,\rm km\,s^{-1}$ and $\sigma_{\sigma_{\mathrm{LOS}}}=19\,\rm km\,s^{-1}$ (H0LiCOW X). We adopt a uniform prior on $r_{\rm ani}$ in a range from 0.5 to 5 times the effective radius, $r_{\rm eff}$, which we calculate to be $r_{\rm eff}=1\aas@@fstack{\prime\prime}41$ from our lens light fitting in the F160W filter. We fit to the double Sérsic light profile, as the Chameleon profile does not provide an accurate representation of the galaxy light distribution at large radii (Dutton et al., 2011). The point source contributes a very small amount to the galaxy light, but not enough to impact this calculation. We note that the choice of filter affects $r_{\rm eff}$, but the impact is small and results in a negligible effect ($\lesssim 0.1\%$) on the final inference.

We use importance sampling (e.g., Lewis & Bridle, 2002) to simultaneously combine the velocity dispersion and external convergence distributions in Section 5 with the $D_{\Delta t}^{\rm model}$ inferred from our lens model. Specifically, for each set of lens and cosmological parameters $\{\bm{\nu},\bm{\pi}\}$ from our lens model MCMC chain, we draw a $\kappa_{\rm ext}$ sample from the distribution in Section 5 and a sample of $r_{\rm ani}$ from the uniform distribution [0.5,5]$r_{\rm eff}$. With these, we can then compute the kinematics likelihood in Equation (10) for the joint sample $\{\bm{\nu},\bm{\pi},\kappa_{\rm ext},r_{\rm ani}\}$ and use this to weight the joint sample. From the effective model time-delay distance computed from our multi-plane lensing ($D_{\Delta t}^{\rm model}$) and the external convergence ($\kappa_{\rm ext}$), we can then compute the effective time-delay distance ($D_{\Delta t}$) via Equation (4), keeping its absolute value blinded until we finalize our analysis. The resulting distribution of $D_{\Delta t}$ encapsulates the cosmological information from WFI2033$-$4723.

4.4 BIC Weighting

We weight our models using the BIC, defined as

[TABLE]

$n$ is the number of data points, which is the number of pixels in the image region across both bands that are outside the fiducial AGN mask (so that we are comparing equal areas), plus eight (for the four AGN image positions), plus three (for the time delays), plus one (for the velocity dispersion). $k$ is the number of free parameters, which is the number of parameters in the lens model that are given uniform priors, plus two (for the source position), plus one (for the anisotropy radius to predict the velocity dispersion). $\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 (which marginalizes over the source surface brightness pixel parameters and is thus the likelihood of the lens/cosmological parameters excluding the source pixel parameters; see Suyu & Halkola, 2010) times the likelihood of the lens model parameters within the image plane region that excludes the arcmask. We evaluate the BIC using the fiducial weight image and arcmask, as the majority of the models were optimized with these. This may penalize the model with a larger AGN mask and the 60x60 pixel source grid model with a larger arcmask, but choosing any other region would penalize the fiducial model and all the other models that used the same regions, so this choice is fair to the largest number of models. We note that our computation of the BIC described above uses all available data sets (lensing image, time delays and lens velocity dispersion) for model comparison.999Ranking the lensing mass models based on BIC values computed from only the lensing data would lead to insignificant changes to the final BIC-weighted $D_{\Delta t}$ distributions.. The kinematics have a relatively small impact in comparison to the other terms, and does not strongly favor either the SPEMD or composite model.

We estimate the variance in the BIC, $\sigma_{\mathrm{\mathrm{BIC}}}^{2}$, by running the fiducial model with source resolutions of [47, 48, 49, 50, 51, 52, 53, 54, 56, 58, 60] pixels on a side (the $50\times 50$ pixel case is just the original fiducial model), keeping the arcmask the same. Changing the source resolution in this way shifts the inferred $D_{\Delta t}$ values stochastically, but there is no overall trend with resolution, and the degree of the shifts are smaller than the scatter among the different models we run. We calculate the BIC for each of these models and take the variance of this set of models as the variance on the BIC, $\sigma_{\mathrm{BIC}}^{2}$. We find $\sigma_{\mathrm{BIC}}^{2}\sim 41$ for the SPEMD models and $\sigma_{\mathrm{BIC}}^{2}\sim 55$ for the composite models. In Appendix D, we show the BIC and BIC weight values for these source resolution tests.

To avoid biases due to our choice of lens model parameterization, we split the samples into the SPEMD and composite models and calculate the relative BIC and weighting for each set separately, similar to H0LiCOW IX. Specifically, we weight a model with a given BIC of $x$ by a function $f_{\mathrm{BIC}}(x)$, defined as the convolution

[TABLE]

where $\mathrm{BIC_{min}}$ is the smallest BIC value within a set of models (SPEMD or composite), and $h$ is a Gaussian centered on $x$ with a variance of $\sigma_{\mathrm{BIC}}^{2}$. We follow the analytic calculation of Yıldırım et al. (2019) in evaluating the convolution integral in Equation (12). Once we have a weighted $D_{\Delta t}$ distribution for the SPEMD models and another for the composite models, we combine these two with equal weight in the final inference.

4.5 Blind Analysis

We perform our analysis blindly using a similar procedure as for previous H0LiCOW analyses (Suyu et al., 2013, H0LiCOW IV, IX). In practice, this is done by subtracting the median of certain parameter PDFs from the distribution when displaying and analyzing results. In particular, we blind the values of $D_{\Delta t}$ or $H_{0}$. This blinding procedure still allows us to measure their precision and relative offsets, as well as their correlation with other lens model parameters, but without knowing their values. Blinding also eliminates confirmation bias and the tendency for experimenters to stop their analysis when they obtain a value consistent with an “expected" value, and forces us to be confident in our checks of systematic errors before finalizing our result. After completing our analysis, writing this paper draft with blinded $D_{\Delta t}$ distributions, and coming to an agreement among the coauthors to unblind the results on 7 May 2019, we unblinded and did not make any further changes to the models. There is also no iteration between the lens modeling and time-delay measurements. Throughout this paper, we show blinded $D_{\Delta t}$ distributions until Section 6, where we reveal the absolute $D_{\Delta t}$ values from our inference.

5 Estimating the external convergence

We estimate the external convergence $\kappa_{\mathrm{ext}}$ using the weighted number counts technique, introduced in Greene et al. (2013) and Rusu et al. (2017, hereafter H0LiCOW III) and reframed as approximate Bayesian computation in H0LiCOW IX. For details of the numerical implementation, and a justification of the applicability in cosmography, we refer the reader to H0LiCOW III. Briefly, we use the catalogue of galaxies and associated physical properties (redshifts and stellar masses) around WFI2033$-$4723 from Section 3.3 and H0LiCOW X, as well as a control catalogue of the same properties from non-contiguous regions of the sky, large enough to overcome sample variance, from the Canada-France-Hawaii Telescope Lensing Survey (CFHTLenS; Heymans et al., 2012). We compute relative galaxy number counts within the matching limiting magnitude and within the same apertures of $45\arcsec$- and $120\arcsec$-radii101010In order to ensure a fair comparison of the number counts between the lens fields and the control field, in H0LiCOW X we performed detections in the same $i$-filter, where the images have similar seeing. To account for the coarser pixel scale in the DES data compared to CFHTLenS, we used more aggressive deblending parameters. This has no noticeable effect, with the exceptions of regions around bright stars. We ignore the negligible differences between the $ugriz$ filter curves in DES and CFHTLenS., using physically motivated weights $\zeta_{q}$ introduced by Greene et al. (2013) and H0LiCOW III. Here $q$ stands for the redshift $z$, stellar mass $M_{\star}$, the inverse of the distance $r$ between each galaxy and the lens or the center of the aperture, etc. A full list of the weights and of the corresponding measured relative number counts is shown in Table 2.

We perform the calculation above three times. These correspond, first, to the fiducial mass model in Section 4.2, which incorporates the nearby galaxies G2, G3 and G7 from Figure 3. Since the effect of these galaxies has already been accounted for, we remove them from the input catalogue before running the computation. Second and third, we also remove one or both of the galaxy groups found in H0LiCOW X to impact the mass modeling beyond the tidal shear term, and therefore taken into account in the systematics tests presented in Section 4.2.3. However, since our spectroscopic completeness down to the limiting magnitude of $i<22.5$ is only $\sim 50\%$, and also non-uniform, decreasing with radius from the lens (see Figure 4 in H0LiCOW X), it is likely that there are other galaxies part of these groups, in addition to the ones spectroscopically confirmed. If we were to keep these galaxies in the number counts, our inferred $\kappa_{\mathrm{ext}}$ would be an overestimate, when coupled to the models from Section 4.2.3 which already include these galaxy groups. We use two different methods to account for these galaxies statistically. Briefly, in the first method we use the measured spectroscopic completeness and the total number of galaxies within the $120\arcsec$-radius aperture, as well as the number of confirmed group members inside the same aperture, and we apply Poisson statistics to infer the distribution of galaxy numbers we miss due to spectroscopic incompleteness. In the second method, we use the velocity dispersions of the two groups measured in H0LiCOW X, as well as the virial radii from Wilson et al. (2016), and we calculate the expected number of galaxies inside the virial radius, from the empirical relation in Andreon & Hurn (2010). Then, based on the projected distance between the group centroid and the lens measured in H0LiCOW X, we estimate the expected number of galaxies at the intersection of the sphere of virial radius and the $120\arcsec$-radius cylinder centered on the lens. Subtracting from this the number of galaxy members spectroscopically confirmed, we arrive at a distribution of the number of missing galaxies. We show the resulting distributions from both methods in Figure 9. For each group, the distributions from both methods overlap significantly, giving consistent results. The expected median number of galaxies missing from the group at $z=0.66$ is 6-8, and from the group at $z=0.49$ it is 3. Finally, we extract at random, 20 times, a number from these distributions, and remove these galaxies, picked at random from within our catalogue of galaxies around the lens, with photometric redshifts compatible with the group redshifts, before computing the weighted number counts. The resulting scatter is included in the values reported in Table 2.

Our measured relative weighted number counts, in particular for the fiducial model and inside the $120\arcsec$-radius aperture, show that the field of WFI2033$-$4723 is overdense. This was also remarked in the study by Fassnacht et al. (2011), where a non-weighted number count overdensity of 1.33 was obtained. While a direct comparison with this values would be biased because of the different limiting magnitude, detection filter, and the fact that Fassnacht et al. (2011) count all nearby galaxies without exception, our $\sim 1.44$ unweighted relative number count inside the same aperture of $45\arcsec$-radius is consistent, within $1.5\sigma$.

To convert the measured relative weighted number counts into a $\kappa_{\textrm{ext}}$ distribution based on these constraints, we follow Suyu et al. (2010, 2013); Greene et al. (2013); H0LiCOW III, IX; Chen et al. (2019), and use the results of ray-tracing by Hilbert et al. (2009) through the Millennium Simulation (MS; Springel et al., 2005), in the form of a convergence and shear map ($\kappa$,$\gamma$) covering each simulated sky location. Our technique is justified by the results of Suyu et al. (2010); Hilbert et al. (2009), which showed that the distribution of $\kappa$ for LOS containing strong lenses is very similar to that over all LOS. With a catalogue of galaxies painted on top of the dark matter halos from the MS, following the semi-analytical models in De Lucia & Blaizot (2007), and containing realistic simulated photometry, we follow a similar procedure and compute relative weighted number counts at each spatial location throughout the MS (see H0LiCOW III for details). Finally, we compute

[TABLE]

where we combine multiple weighted number count constraints $\zeta_{q}$, including from both $45\arcsec$- and $120\arcsec$-radius apertures. Following H0LiCOW III, we treat the external shear $\gamma_{\textrm{ext}}$ computed from the lens models in Section 4 at the location of the lens (in the case of the MS, at the center of each aperture) analogously to the weighted number count constraints.111111Here and in Chen et al. (2019) we modify the way we implement the $\gamma_{\textrm{ext}}$ constraint described in H0LiCOW III, in the sense that we no longer normalize by the number of LOS in each small division of the constraint range. This is because the shear values derived in Section 4 use a flat prior, and the distribution of $\gamma$ in the MS maps naturally introduces a cosmological prior. This effect is negligible, except for the case of very large $\gamma_{\textrm{ext}}$ uncertainties. In Section E we explore various combinations of constraints, show that our technique is free of biases, and describe several tests we ran. We settle on the combination of $P(\kappa_{\mathrm{ext}}|\zeta^{45\arcsec}_{\mathrm{1}},\zeta^{45\arcsec}_{\mathrm{1/r}},\zeta^{120\arcsec}_{\mathrm{1}},\zeta^{120\arcsec}_{\mathrm{1/r}},\gamma)$, which employs our most robust constraints. In Figure 10 we show the resulting distributions, corresponding to the various mass models explored in Section 4, and their associated shear values121212It is unexpected that the two distributions of the composite model, which correspond to the case where one or both galaxy groups are explicitly modeled, have larger medians than the standard composite model (by $\lesssim 0.2\sigma$ or at $\lesssim 1\%$ level), even though they are constrained by smaller values of shear and weighted counts. We attribute this to noise, and we have checked that the excess is consistent with variations between similar distributions for these models, employing a different choice of weighted count constraints..

6 Results

After conducting the analysis described in Section 4 and Section 5, and combining with the time delays from Section 3.4, we plot the final BIC-weighted $D_{\Delta t}$ distributions in Figure 11, with the blinded values shown on the bottom axis and the unblinded values shown on the top axis. We report the median and 68% quantiles of $D_{\Delta t}$ for each of the models Table 3, along with the $\Delta$BIC and associated weighting that each model receives. Our constraint on $D_{\Delta t}$ for WFI2033$-$4723 is $D_{\Delta t}=4784_{-248}^{+399}~{}\mathrm{Mpc}$, a $\sim 6.6\%$ precision measurement.

From this inferred $D_{\Delta t}$, we can calculate cosmological parameters for flat $\Lambda$CDM or other cosmologies. For flat $\Lambda$CDM with uniform priors on $H_{0}$ (within the range [0, 150] km s*-1* Mpc*-1*) and $\Omega_{\rm m}$ (within the range [0.05, 0.5]), this translates into a constraint on the Hubble constant of $H_{0}=71.6_{-4.9}^{+3.8}~{}\mathrm{km~{}s^{-1}~{}Mpc^{-1}}$. Within the uncertainties, our result is consistent with the previous measurements of $H_{0}$ from H0LiCOW. After B1608$+$656 (Suyu et al., 2010), this is the H0LiCOW lens producing the second smallest uncertainty on $H_{0}$, comparable to what is expected for the lensed supernovae “Refsdal” (Grillo et al., 2018). We note that weighting all of our models equally and ignoring the BIC weighting (as was done with some previous H0LiCOW lenses) gives $H_{0}=72.2_{-4.8}^{+4.3}~{}\mathrm{km~{}s^{-1}~{}Mpc^{-1}}$, which is a shift of $<1\%$ in the median value. Our result can be combined with the other five lenses in the H0LiCOW sample to give stronger constraints for a variety of cosmological models. The full cosmological analysis is presented in Wong et al. (2019). Our measurement is also consistent with recent ones from other techniques, not only from Planck and SH0ES, but also from the Carnegie Supernova Project (Burns et al., 2018), the Megamaser Cosmology project (Braatz et al., 2018), DES clustering and weak lensing + baryon acoustic oscillations + Big Bang nucleosynthesis experiments (Abbott et al., 2018), the inverse distance ladder (e.g., Aubourg et al., 2015; Macaulay et al., 2019), extragalactic background light attenuation (e.g., Domínguez et al., 2019), etc.

We note that WFI2033$-$4723 has been used in the past to measure $H_{0}$. Based on three years of monitoring and constrained by the relative quasar image positions measured from archival HST imaging, Vuissoz et al. (2008) infer $H_{0}=67_{-10}^{+13}~{}\mathrm{km~{}s^{-1}~{}Mpc^{-1}}$, using non-parametric modeling, in good agreement with our result, but with significantly larger uncertainties.

7 Conclusions

We have analyzed the gravitational lens WFI2033$-$4723, performing a blind cosmographic analysis to determine the time-delay distance of this system. We use deep HST imaging, precise time-delay measurements, a measurement of the lens galaxy’s velocity dispersion, and deep wide-area spectroscopic and photometric data to constrain the mass distribution along the LOS. By accurately modeling the lens and accounting for systematic uncertainties, we constrain the effective time-delay distance to be $D_{\Delta t}=4784_{-248}^{+399}~{}\mathrm{Mpc}$, a precision of $6.6\%$. This translates to a Hubble constant of $H_{0}=71.6_{-4.9}^{+3.8}~{}\mathrm{km~{}s^{-1}~{}Mpc^{-1}}$ in a flat $\Lambda$CDM cosmology with uniform priors on $H_{0}$ and $\Omega_{\rm m}$, a value consistent with measurements from other techniques, as well as previous H0LiCOW lenses. A joint analysis of all six H0LiCOW lenses and our constraints on different cosmologies is presented in Wong et al. (2019).

As with all galaxy-scale lenses where time-delay cosmography is applied, we expect that our results can be improved with future, higher resolution adaptive optics imaging (Chen et al., 2016), spatially resolved kinematics (Shajib et al., 2018), and a more tailored, non-statistical treatment of the external convergence (McCully et al., 2017).

Acknowledgements

We thank Adriano Agnello, Roger Blandford, Xuheng Ding, Georges Meylan, Danka Paraficz, Chiara Spiniello, Malte Tewes and Olga Tihhonova for their contributions to the H0LiCOW project. H0LiCOW and COSMOGRAIL are made possible thanks to the continuous work of all observers and technical staff obtaining the monitoring observations, in particular at the Swiss Euler telescope at La Silla Observatory. Euler is supported by the Swiss National Science Foundation. C.E.R. and C.D.F. are funded through the NSF grant AST-1312329, “Collaborative Research: Accurate cosmology with strong gravitational lens time delays,” and the HST grant GO-12889. 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. V.B. and F.C. acknowledge support from the Swiss National Science Foundation (SNSF) and through European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (COSMICLENS: grant agreement No 787866) D.S. acknowledges funding support from a Back to Belgium grant from the Belgian Federal Science Policy (BELSPO). S.H.S. thanks the Max Planck Society for support through the Max Planck Research Group. This research was supported in part by Perimeter Institute for Theoretical Physics. Research at Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development and by the Province of Ontario through the Ministry of Research, Innovation and Science. S.H. acknowledges support by the DFG cluster of excellence ‘Origin and Structure of the Universe’ (www.universe-cluster.de). C.D.F. and G.C.-F.C. acknowledge support for this work from the National Science Foundation under Grant No. AST-1715611. T.T. thanks the Packard Foundation for generous support through a Packard Research Fellowship, the NSF for funding through NSF grant AST-1450141, “Collaborative Research: Accurate cosmology with strong gravitational lens time delays". G. C.-F. Chen acknowledges support from the Ministry of Education in Taiwan via Government Scholarship to Study Abroad (GSSA) L.V.E.K. is supported in part through an NWO-VICI career grant (project number 639.043.308). P.J.M. acknowledges support from the U.S. Department of Energy under contract number DE-AC02-76SF00515.

This paper is based on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with Programs #12889 and #9744. Support for program #12889 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555.

This work is partly based on the analysis in (Sluse et al., 2019), which used public archival data from the Dark Energy Survey (DES). Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, the Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Científico e Tecnológico and the Ministério da Ciência, Tecnologia e Inovação, the Deutsche Forschungsgemeinschaft, and the Collaborating Institutions in the Dark Energy Survey. The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenössische Technische Hochschule (ETH) Zürich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciències de l’Espai (IEEC/CSIC), the Institut de Física d’Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig-Maximilians Universität München and the associated Excellence Cluster Universe, the University of Michigan, the National Optical Astronomy Observatory, the University of Nottingham, The Ohio State University, the OzDES Membership Consortium, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, and Texas A&M University. Based in part on observations at Cerro Tololo Inter-American Observatory, National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.

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.

Data analysis was in part carried out on common use data analysis computer system at the Astronomy Data Center, ADC, of the National Astronomical Observatory of Japan. This work made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration et al., 2013, 2018). Plots were produced with Matplotlib (Hunter, 2007).

Appendix A Residuals without power-law weighting

For completeness, we show the normalized residuals for the fiducial SPEMD model (Figure 12) and fiducial composite model (Figure 13) using the weight images without the power-law weighting in the region near the AGN images. We see that there are strong residuals due to the AGN images, which motivates our downweighting of these regions.

Appendix B Impact of Different Cosmologies

In multi-lens-plane modeling, we need to sample the cosmological parameters in order to carry out the ray tracing. For computational reasons, we directly vary $H_{0}$ but keep other cosmological parameters fixed ($\Omega_{\rm m}=0.3$, $\Omega_{\Lambda}=0.7$, $w=-1$). $D_{\Delta t}$ has a weak dependence on these other parameters in principle, but we found that varying these parameters affected the posterior $D_{\Delta t}$ distribution by $<1\%$ for HE 0435$-$1223 (H0LiCOW IV). We perform a similar check for WFI2033$-$4723 in which we run the fiducial SPEMD model while allowing either $\Omega_{\rm m}$ to vary, or allowing both $\Omega_{\rm m}$ and $w$ to vary. The resulting effective $D_{\Delta t}$ distributions are shown in Figure 14. The peaks of the distribution are consistent to within $1\%$ of the absolute value, which shows that the results are insensitive to these extra cosmological parameters at the level of accuracy that we are currently working at, similar to HE 0435$-$1223.

In the case of multiple lens planes (e.g., Blandford & Narayan, 1986; Kovner, 1987; Kochanek & Apostolakis, 1988; Schneider et al., 1992; Petters et al., 2001; Collett & Auger, 2014; McCully et al., 2014; Schneider, 2014), there is not a unique time-delay distance for the system, but instead multiple time-delay distances between planes $i$ and $j$,

[TABLE]

with $z_{i}$ being the redshift of plane $i$. The multi-plane time delay is given by

[TABLE]

where $s$ is the index of the source plane (starting with $i=1$ as the lowest-redshift lens plane and counting up towards the source plane). $\psi_{i}$ is the lens potential related to the scaled deflection angle via $\nabla\psi_{i}=\bm{\alpha}_{i}$, and $\beta$ is a ratio of angular diameter distances among deflector planes and the source,

[TABLE]

From Equation (15), we see that the time delay depends on the multiple time-delay distances and $\beta$ terms. In general, it is difficult to constrain all of these quantities independently, so we adopt specific cosmological models to compute the distances for the ray tracing, then compare the time-delay distance measurements from these different background cosmologies. For the case of WFI2033$-$4723, where G2 is at a different redshift from the main lens plane and is not strongly lensing the background source but merely perturbing it, the effect on the time delays is weak. The lack of sensitivity to $\Omega_{\rm m}$ and $w$ seen in Figure 14 suggests that WFI2033$-$4723 is not sensitive to the extra $\beta$ terms at an interesting level to probe it directly in the same way as a double source plane lens (e.g., Gavazzi et al., 2008; Collett & Auger, 2014). Since the time delays are mostly set by the main lens plane, we can measure the “effective" $D_{\Delta t}$ (which is $D_{\Delta t}^{is}$ with $i$ as the main lens plane and $s$ as the source plane) that is independent of assumptions on the background cosmology. This robust distance determination then permits us to constrain any reasonable cosmological model via the distance-redshift relation.

Appendix C Model Parameters

We show the marginalized parameter constraints for each of the SPEMD models in Table 4 and for each of the composite models in Table 5.

Appendix D Source resolution changes

In Table 6, we show the BIC and BIC weight values for the source resolution tests described in Section 4.4. All models are the fiducial models run with different source resolutions. The BIC weights have been renormalized.

Appendix E Further tests on the selection and implementation of the weighted number count constraints

E.1 The use of $\gamma_{\mathrm{ext}}$ constraints

In Section 4 we have shown that a large value of $\gamma_{\mathrm{ext}}\sim 0.10-0.14$ is required to model the lens system, even after the nearby perturbers are being accounted for. In addition, the different inner mass profiles of the lens, explored in Section 4, all require a large external shear. This means that the large-scale environment and/or LOS structures must be responsible for this effect.

The specific large-scale structures responsible for this effect have not been identified, since the mass models incorporating the two large galaxy groups identified in H0LiCOW X still require a large $\gamma_{\mathrm{ext}}$. This was also pointed out in previous studies, based on more limited ancillary data and less detailed modeling, both superseeded in our work. For example, Wong et al. (2011) compute an expected $\gamma_{\mathrm{ext}}=0.08\pm 0.03$ based on a spectroscopic galaxy catalogue, where galaxies are treated as SIS, and galaxy groups as NFW halos. They note that the orientation of this shear caused by the environment does not match the one obtained from the mass models of the lens. Vuissoz et al. (2008), based on mass models constrained only by the positions of the quasar images (and in some models by the measured time delays as well), find $\gamma_{\mathrm{ext}}$ as large as 0.3, but as small as 0.06 if they incorporate the galaxy group which includes the lens. However, our more complete spectroscopic catalogue shows that the group centroid is more distant from the lens ($\sim 30\arcsec$ compared to $\sim 10\arcsec$), leading to a larger amount of external shear necessary to model the system, after accounting for the group.

Based on the above, while we concede that we cannot find a mass model which explains most of the shear, we are confident that the external shear values measured for the various models in Section 4 are robust, and that we are justified to attribute the shear to the lens environment and/or LOS structures which are not captured in our mass models, but are captured in our statistical approach to computing $\kappa_{\mathrm{ext}}$. This justifies our use of the $\gamma_{\mathrm{ext}}$ constraints to infer $\kappa_{\mathrm{ext}}$ in Figure 10. As we will show in Appendix E.2, this constraint has a dominant effect on our inference, which is expected from results of ray-tracing through the MS (e.g., see Figure 8 in Collett & Cunnington, 2016), where $|\kappa|$ and $|\gamma|$ are found to correlate. We note that the lens RXJ1131$-$1231, which has also been modeled as part of H0LiCOW (Suyu et al., 2013; Chen et al., 2019), also has a fairly large measured shear of $\sim 0.08$, which was used to constrain its $\kappa_{\mathrm{ext}}$. Independently, its $\kappa_{\mathrm{ext}}$ was measured using a different methodology by McCully et al. (2017), which has found $\kappa_{\mathrm{ext}}$ to be offset to smaller, but nonetheless consistent values.

After unblinding our analysis, we checked what the impact of using $\kappa_{\mathrm{ext}}$ inferred without the shear constraint would have been on our analysis. We obtain $H_{0}=76_{-3.0}^{+2.9}~{}\mathrm{km~{}s^{-1}~{}Mpc^{-1}}$, a value 6% larger than the one in Section 6, and with significantly increased statistical precision due to the tighter $\kappa_{\mathrm{ext}}$ distributions (see Figure 15).

In order to avoid biases, we must ensure that we construct $P(\kappa_{\mathrm{ext}}|\gamma,...)$ such that it is consistent with the discussion above. In particular, we must ensure that when we select LOS with large shear from the MS, the shear is not due to galaxies very close to the LOS, or to galaxy groups/clusters so massive that we would incorporate them in our mass models, in the real data. We address each of these in the following.

E.1.1 Galaxies close to the LOS

In order to ensure that the weighted number counts are not dominated by the galaxies very close to the LOS, in Greene et al. (2013); H0LiCOW III, IX we have used a $5\arcsec$-radius inner mask, and set an upper limit to the weights incorporating $1/r$ of $1/10\arcsec$ for each galaxy. While the usage of a mask when computing weighted number counts is agnostic to the actual existence of galaxies inside of it, and therefore their contribution to $\kappa_{\mathrm{ext}}$ at that particular spatial location, here and in Chen et al. (2019) we only select LOS from the MS which have no galaxies inside the $5\arcsec$-radius. Such galaxies, if above the magnitude threshold, would be modeled explicitly in the real data. This radius corresponds to the inner $\gtrsim 4$ pixels of the $(\kappa,\gamma)$ map from Hilbert et al. (2009). The addition of this constraint has the effect of lowering $\kappa_{\mathrm{ext}}$ by approximately $0.1\kappa_{\mathrm{ext}}$.

E.1.2 Massive large-scale structures

For the FOV around WFI2033$-$4723, we used our spectroscopic catalogue to identify large galaxy groups, which we incorporated in our mass models, in cases where they exceeded the flexion shift threshold. We must take this fact into account when we infer $\kappa_{\mathrm{ext}}$ by selecting LOS from the MS. The easiest way to account for this in the MS is to use a complete catalogue of galaxy groups, and to test at each location of the $\kappa$ map whether their flexion shifts exceed the threshold. If so, those LOS are removed. In this way, we only use the $\kappa_{\mathrm{ext}}$ distribution free of the contribution of massive structures, which is what is needed.

In practice, we use a catalogue of galaxies from the MS, which identifies their parent halos and the masses of those halos. We impose the observational constraints from WFI2033$-$4723, that the groups (in this case the parent halos) we model have at least 5 galaxies within our magnitude threshold. For these haloes, we convert their masses into velocity dispersions (assuming the SIS mass profile, for simplicity), and compute the radius of the circle around the halo centroid inside which the flexion shift for each of these structures exceeds the threshold. We then remove the LOS inside those circles from the $\kappa_{\mathrm{ext}}$ computation. We have compared the $\kappa_{\mathrm{ext}}$ distributions with and without removing these groups, and have found them to be indistinguishable, given the rarity of these groups. We therefore safely ignore the effect of massive but rare large scale structures on our analysis.

E.2 The choice of conjoined constraints

Greene et al. (2013) and H0LiCOW III have explored different combinations of weighted number count ratio constraints, most of which incorporate unweighted number counts $\zeta_{1}$, as well as $\zeta_{1/r}$. These correspond to the most robust constraints which can be determined from imaging data, usually with the tightest uncertainties. In Figure 15 we explore the median and standard deviation of the $\kappa_{\mathrm{ext}}$ distributions for various combinations of constraints, typically incorporating the two constraints above. We explore combinations of constraints measured inside the same aperture, as well as measured inside both the $45\arcsec$- and $120\arcsec$-radii apertures. The latter is because the weighted number count ratios we measure in Table 2 appear larger inside the $120\arcsec$-radius aperture, suggesting that WFI2033$-$4723 is positioned on the outskirts of a galaxy overdensity, and this observation might include statistical information useful for tightening the $\kappa_{\mathrm{ext}}$ distribution.

Comparing the distributions constrained by $(\zeta^{45\arcsec}_{1},\zeta^{45\arcsec}_{1/r},\zeta^{45\arcsec}_{q})$ and $(\zeta^{120\arcsec}_{1},\zeta^{120\arcsec}_{1/r},\zeta^{120\arcsec}_{q})$ we see that the former produce lower $\kappa^{\mathrm{med}}_{\mathrm{ext}}$. This is partly because of the somewhat lower weighted number counts inside the smaller aperture, but mostly due to the fact that the same overdensity over a larger aperture implies a larger structure, and therefore larger convergence. This is demonstrated in Figure 16. In the case of constraints from both apertures, such as $(\zeta^{45\arcsec}_{1},\zeta^{45\arcsec}_{1/r},\zeta^{120\arcsec}_{1},\zeta^{120\arcsec}_{1/r},\zeta^{120\arcsec}_{q})$, $\kappa^{\mathrm{med}}_{\mathrm{ext}}$ is brought closer to the average value obtained from the two individual apertures. Amongst the various $\zeta_{q}$, $\kappa^{\mathrm{med}}_{\mathrm{ext}}$ varies by $\sim 0.01$, or at the $1\%$ level.

Once the large $\gamma_{\mathrm{ext}}$ is used as a constraint, $\kappa^{\mathrm{med}}_{\mathrm{ext}}$ reaches much larger values, as expected. All distributions including the $\zeta_{1/r}$ constraint are in agreement to within the $1\%$ level, even if we compare constraints from different aperture radii. Interestingly, once the constraints from different apertures are combined, the distributions for different $\zeta_{q}$ remain in agreement, but $\kappa^{\mathrm{med}}_{\mathrm{ext}}$ decreases by almost $\sim 2\%$, or $\sim 0.3\sigma$. We cannot fully explain this result, but we note that we have explored the $\kappa_{\mathrm{ext}}-\gamma_{\mathrm{ext}}$ plane constrained by either, as well as both apertures. It appears that this behavior is a result of both the very large $\gamma_{\mathrm{ext}}$ value and of the large weighted count constraints, and disappears for smaller values of either.

We do not use more than 5 conjoined constraints, because the number of MS LOS surviving the cut becomes too small, and the distributions are noisy. Our conclusions regarding the convergence distributions from different constraints are similar to those from H0LiCOW III, wherever a direct comparison is possible.

E.3 Tests for bias and scatter

In H0LiCOW III, we ran two types of simulations using the MS data, to check whether or not our combination of constraints biases $P(\kappa_{\mathrm{ext}})$. We ran these for the case of WFI2033$-$4723 as well, and found that $P(\kappa^{\mathrm{med}}_{\mathrm{ext}}-\kappa^{\mathrm{true}}_{\mathrm{ext}}|\zeta_{q},...)$ is centered on zero for all $\zeta_{q},...$ combinations in Figures 15 and 17. All our $P(\kappa_{\mathrm{ext}})$ distributions are, therefore, unbiased.

We also wish to determine which combinations of constraints produces the tightest $P(\kappa_{\mathrm{ext}})$. However, from Greene et al. (2013) as well as Figure 15 we see that the medians and standard deviations of $P(\kappa_{\mathrm{ext}})$ are always correlated, such that smaller $\kappa^{\mathrm{med}}_{\mathrm{ext}}$ implies smaller $\sigma_{\kappa}$. We therefore follow a different approach, where we use the scatter in $P(\kappa^{\mathrm{med}}_{\mathrm{ext}}-\kappa^{\mathrm{true}}_{\mathrm{ext}}|\zeta_{q},...)$. We normalize the scatter in $P(\kappa^{\mathrm{med}}_{\mathrm{ext}}-\kappa^{\mathrm{true}}_{\mathrm{ext}}|\zeta^{45\arcsec}_{1})$151515Using instead $P(\kappa^{\mathrm{med}}_{\mathrm{ext}}-\kappa^{\mathrm{true}}_{\mathrm{ext}}|\zeta^{120\arcsec}_{1})$ produces consistent results. to the unit value, and we fit the standard deviation of $P(\kappa_{\mathrm{ext}}|\zeta^{45\arcsec}_{1})$ as a function of its median with a linear function, while we vary $\zeta^{45\arcsec}_{1}$. Finally, we divide the standard deviation of $P(\kappa^{\mathrm{med}}_{\mathrm{ext}}-\kappa^{\mathrm{true}}_{\mathrm{ext}}|\zeta_{q},...)$ by the value of this linear function at the corresponding $\kappa^{\mathrm{med}}_{\mathrm{ext}}$, and we show the results, for selected combinations of constraints, in Figure 17. We can see that combining the results from the two different aperture radii typically results in scaled standard deviations reduced by up to $\sim 20\%$, whether the $\gamma_{\mathrm{ext}}$ constraint is used or not. While the figure shows a spike in value for the scaled standard deviation corresponding to the combination of weights we choose to consider as fiducial in this work, $(\zeta^{45\arcsec}_{1},\zeta^{45\arcsec}_{1/r},\zeta^{120\arcsec}_{1},\zeta^{120\arcsec}_{1/r},\gamma)$, most of the surrounding similar distributions (with only one of the $\zeta_{q}$ constraints replaced) show small values, and we know from Figure 15 that the resulting $P(\kappa_{\mathrm{ext}})$ for these constraints are consistent with each other.