Hubble Constant Measurement from Three Large-separation Quasars Strongly Lensed by Galaxy Clusters

Tension between cosmic microwave background–based and distance ladder–based determinations of the Hubble constant H 0 motivates the pursuit of independent methods that are not subject to the same systematic effects. A promising alternative, proposed by Refsdal in 1964, relies on the inverse scaling of H 0 with the delay between the arrival times of at least two images of a strongly lensed variable source such as a quasar. To date, Refsdal’s method has mostly been applied to quasars lensed by individual galaxies rather than by galaxy clusters. Using the three quasars strongly lensed by galaxy clusters (SDSS J1004+4112, SDSS J1029+2623, and SDSS J2222+2745) that have both multiband Hubble Space Telescope data and published time delay measurements, we derive H 0, accounting for the systematic and statistical sources of uncertainty. While a single time delay measurement does not yield a well-constrained H 0 value, analyzing the systems together tightens the constraint. Combining the six time delays measured in the three cluster-lensed quasars gives H 0 = 74.1 ± 8.0 km s−1 Mpc−1. To reach 1% uncertainty in H 0, we estimate that a sample size of order of 620 time delay measurements of similar quality as those from SDSS J1004+4112, SDSS J1029+2623, and SDSS J2222+2745 would be needed. Improving the lens modeling uncertainties by a factor of two and a half may reduce the needed sample size to 100 time delays, potentially reachable in the next decade.


INTRODUCTION
The Hubble parameter H 0 , which describes the current expansion rate of the Universe, has been sought since the discovery in the 1920s that the Universe is expanding (Lemaître 1927;Hubble 1929).At the turn of the last century, measurements of H 0 started converging around H 0 = 70 km s −1 Mpc −1 .However, as H 0 measurements have become increasingly precise, the so-called 'Hubble Tension' has arisen between the estimates from early-and late-Universe probes.The Planck Collaboration reported H 0 = 67.4± 0.5 km s −1 Mpc −1 (Planck Collaboration et al. 2020).They used density fluctuations encoded in the Cosmic Microwave Background (CMB) at the surface of last scattering to determine H at that epoch and then used a spatially flat cosmological model to extrapolate to H 0 .By contrast, the "Supernovae, H 0 , for the Equation of State of Dark Energy" (SH0ES) collaboration combined Gaia parallaxes and multiband HST photometry of Milky Way Cepheids to calibrate the extragalactic distance scale and derive H 0 = 73.04 ± 1.04 km s −1 Mpc −1 (Riess et al. 2022).The Planck and SH0ES values, which respectively capture the early and late-time physics of the Universe, differ by 5σ.Freedman (2021) applied an updated Tip of the Red Giant Branch (TRGB) calibration to a distant sample of Type Ia supernovae from the Carnegie Supernova Project and obtained H 0 = 69.8± 0.6 (stat) ± 1.6 (sys) km s −1 Mpc −1 , consistent with the CMB value, and within 2σ of the SH0ES value, owing to the larger uncertainties.The discrepancy between different H 0 methods may indicate a deviation from the standard Λ Cold Dark Matter (ΛCDM) model, and therefore new physics, or the presence of unknown or underestimated systematics.Either way, this tension remotivates the pursuit of other H 0 determination methods that are not prone to the same systematics.
An alternative H 0 determination method, proposed by Refsdal (1964), uses the fact that H 0 scales inversely with the delay between the arrival times of at least two images of a strongly-lensed variable source, such as a quasar or a supernova.Due to the rarity of galaxy clusters lensing quasars or supernovae, the Refsdal H 0 technique has primarily been sought with galaxy-scale lenses (see e.g., the recent reviews by Moresco et al. 2022;Birrer et al. 2022;Treu et al. 2022).
Of the >300 known lensed quasars, the vast majority are lensed by individual galaxies (Lemon et al. 2019(Lemon et al. , 2022)).Quasars lensed by individual galaxies have been used to obtain H 0 .For example, the H 0 Lenses in COS-MOGRAIL's Wellspring (H0LiCOW) collaboration performed a joint analysis of six galaxy-lensed quasars to obtain H 0 = 73.3+1.7  −1.8 km s −1 Mpc −1 (Wong et al. 2020).This value is consistent with the Cepheid-calibrated measurement from the SH0ES collaboration.
The main uncertainty with galaxy-scale lenses is the mass-sheet degeneracy which does not allow the total density profile to be precisely constrained from the lensing information alone.Kochanek (2020) argued that the H0LiCOW H 0 measurement uncertainty was underestimated due to artificially restricting the mass-sheet degeneracy through specific mass model assumptions.To completely account for the mass-sheet degeneracy in the uncertainty, Birrer et al. (2020) permitted additional freedom in the lens model and constrained the mass profile with stellar kinematics.Unsurprisingly, this led to an increase in the uncertainty level to ∼8%, but the mean H 0 value remained at ∼73 km s −1 Mpc −1 .To reduce the uncertainty, Birrer et al. (2020) included further information on the galaxy mass profile from the Sloan Lens ACS sample and obtained H 0 = 67.4± 3.7 km s −1 Mpc −1 .Resolved JWST and Extremely Large Telescope observations of the stellar kinematics in the lens galaxies may significantly reduce these sources of systematic errors (Birrer & Treu 2021).
What has remained largely unexplored until now is determining H 0 by using quasars that are strongly lensed by galaxy clusters.For several reasons, cluster-lensed quasars can potentially overcome some of the difficulties faced by individual galaxy lenses.First, since galaxy clusters have deeper potential wells than galaxies, cluster lenses produce longer time delays of order months to years compared to typically a month in galaxy lenses.Consequently, the observationally measured time delay values will have smaller fractional uncertainty, which then will propagate to reduced uncertainty in H 0 due to the inverse scaling of H 0 with time delays.Second, the light curves of cluster-lensed quasars are less likely to be affected by microlensing from stars in the lens plane, because the mass distribution is dominated by dark matter at the projected radius at which the images appear.Third, galaxy cluster mass distributions are less affected by complex baryonic physics than those of galaxy lenses; the complex baryonic surface density of galaxy-scale lenses may be a significant source of systematic uncertainty.Cluster-scale lenses are less affected by the mass-sheet degeneracy due to having multiple images of lensed sources from various source plane redshifts behind the lens (Grillo et al. 2020).A challenge that must be contended with, however, is the complexity of cluster lenses.
Two inputs are necessary to use cluster-lensed quasars to determine H 0 .The first is an observational measurement of the time delay between the multiple quasar images, and the second is an accurate mapping of the projected density of the dark and luminous mass at the cluster core.High accuracy lens models require space-based resolution and spectroscopic follow-up.Of the seven published cluster-lensed quasars to date (Inada et al. 2003(Inada et al. , 2006;;Dahle et al. 2013;Shu et al. 2018Shu et al. , 2019;;Martinez et al. 2022;Napier et al. 2023), only three have the necessary data available to determine H 0 : SDSS J1004+4112, SDSS J1029+2623, and SDSS J2222+2745.In this paper, we use the available archival HST data and the latest measurements of time delay and spectroscopic redshifts of background sources from the literature to obtain an independent measurement of H 0 from these three systems.
This paper is organized as follows: In Section 2, we outline the theory of observational gravitational lensing time delay and its dependence on H 0 .In Section 3 we detail the lens modeling procedure.In Sections 4, 5, and 6, we give an overview of the three cluster-lensed quasar systems used in this H 0 analysis and provide details about their HST and spectroscopic data, time delays, and lens models.In Section 7, we present our constraints on H 0 .We conclude in Section 8 with a discussion of our H 0 result and the future prospects of the time delay H 0 method.
Throughout the paper, we adopt the standard ΛCDM flat cosmological model with Ω m = 0.3 and Ω Λ = 0.7.

TIME DELAY ANALYSIS
The Refsdal H 0 method is possible due to the measurable delay between the arrival time of two or more images of a variable source such as a quasar.Under the thin lens approximation, a packet of light that travels from the source to the observer will be delayed by time t given by the arrival time surface (Schneider 1985): (1) where z l is the redshift of the lens, d l , d s , and d ls are the angular diameter distances from the observer to the lens, to the source, and between the lens and the source, respectively; ⃗ θ is the image position in the image plane; ⃗ β is the unobserved source position; and ψ( ⃗ θ) is the gravitational lensing potential.The arrival time t is a combination of the path length and the gravitational time delay (t = t geometric + t grav ).The last term, τ (θ; , is also known as the Fermat potential.The multiple images of a stronglylensed source appear in the stationary points of the arrival time surface, that is, in the minima, maxima, and saddle points.H 0 is incorporated in Eq. 1 because of its inverse scaling with the angular diameter distances: for flat cosmology, where The matter density and vacuum energy density parameters are Ω m and Ω Λ , respectively.Conveniently, H 0 is disentangled from the other cosmological parameters in the angular diameter distance equation (Eq.2).After substituting Eq. 2 into d l d s /d ls in Eq. 1, the time delay is determined by applying Eq. 1 for two image positions corresponding to the same source position and taking the difference.The time delay between the images thus becomes: The first term on the right-hand side of the time delay equation gives the Hubble parameter; the second term contains the dependence on cosmological parameters other than H 0 ; and the third term is solved by the strong gravitational lens model.We neglect the higher order effects of the cosmological parameters and fix the second term on the right-hand side in Eq. 3 using the fiducial cosmology.The left-hand side of the equation is the measurement of the time delay, e.g., from monitoring and comparing the observed light curves of two images of the variable source.
Once we compute a model of the lensing mass distribution (see Section 3), we determine the model-predicted excess arrival time surface (Eq. 3) with respect to one of the quasar images.Assuming our lens model is a correct description of the matter distribution, we then leverage the fact that the time delay scales inversely with H 0 .We compare the model-predicted time delays between images to the observational measurements of the time delays to obtain H 0 via: where H 0,model is the H 0 value used to generate the Fermat potential from the lensing analysis (we assumed H 0,model = 70 km s −1 Mpc −1 ), ∆t model is the modelpredicted time delay between the quasar images, and ∆t measured is the observational measurement of the time delay between the pair of quasar images.

LENS MODELING
We computed the lens models with the publicly available software Lenstool (Jullo et al. 2007).Lenstool is a 'parametric' modeling algorithm which describes the lensing mass distribution as a linear combination of galaxy-scale, group-scale, and cluster-scale halos, each of which is parameterized as a pseudo-isothermal ellipsoidal mass distribution (PIEMD, also called dPIE; Elíasdóttir et al. 2007).The three-dimensional density distribution of the PIEMD, has seven parameters whose values can either be fixed or varied: position (x, y); ellipticity e = (a 2 -b 2 )/(a 2 +b 2 ), where a and b are the semi-major and semi-minor axes, respectively; position angle θ; core radius r c ; truncation radius r cut ; and effective velocity dispersion σ 0 .The parameters of the group-scale and cluster-scale halos are typically allowed to vary.The exception is r cut for the cluster-scale halos as this radius usually occurs outside the region where strong lensing evidence is found, and thus, cannot be constrained.
Lenstool uses a Markov Chain Monte Carlo (MCMC) sampling of parameter space.The best-fit model is identified as the one that minimizes the scatter between the model-predicted and observed image locations in the image plane ("image plane minimization") or minimizes the scatter between the predicted source locations of multiple images in the source plane ("source plane minimization").The lens models employ the strong lensing information of multiply-imaged galaxies (arcs), whose positions and redshifts are used as model constraints.The availability of lensing constraints strongly affects the accuracy of lens models, as they are used as local solutions of the lensing equations and constrain the projected mass density distribution at the cluster's core.The mass distribution and magnification are sensitive to the accurate identifications and positions of multiple images and to the redshifts of the lensed galaxies.It is necessary to include a few spectroscopic redshifts in the lens model to avoid incorrect results (Johnson & Sharon 2016).As described in more detail in the following sections, all of the multiply-imaged galaxies used as constraints in the lens models for SDSS J1004+4112, SDSS J1029+2623, and SDSS J2222+2745 were robustly determined; hence the image identifications are not a significant source of systematic uncertainty.
To select cluster-member galaxies, we followed the procedure of Gladders & Yee (2000), by selecting galaxies that fall on the cluster red sequence in a colormagnitude diagram.For SDSS J1029+2623 we also incorporated spectroscopic redshift information (see Section 5).The galaxy-scale halos' positional parameters (x, y, e, θ) are measured with Source Extractor (Bertin & Arnouts 1996) and fixed.Ther c , r cut , and σ 0 of the galaxy-scale halos are scaled to their observed luminosity following the scaling relations in Limousin et al. (2005).A potential misidentification of cluster member galaxies is unlikely to be a significant source of systematic uncertainty.For example, Napier et al. (2023) found that excluding a subsample of cluster members from the lens model did not significantly change the predicted time delays.
We modeled SDSS J1004+4112 using one cluster-scale halo, one brightest cluster galaxy (BCG)-scale halo, and a galaxy-scale halo for each of the cluster member galaxies, 168 in total, five (including the BCG) of which have their slope parameters optimized instead of adopting the scaling relations from Limousin et al. (2005).We optimized the parameters for these five galaxies because of either their proximity to the quasar images or their necessity for reproducing the lensed image configuration.
We modeled the cluster using both source-plane minimization and image-plane minimization, and evaluated the quality of the models obtained by each approach.While formally the image-plane minimization resulted in a better image-plane scatter, these models produced additional quasar images that are not observed.Therefore, we proceeded with the source-plane minimization for SDSS J1004+4112 for the remainder of the analysis.We note that the best-fit lens model produced large scatter between the observed and model-predicted positions in the image plane for quasar image C. In our results, we checked what happens when image C is removed from the H 0 measurement.
The model consists of 27 free parameters and 78 constraints.The HST data and the lens model for SDSS J1004+4112 are shown in Figure 1.The redshifts of the arcs in our lens model are the same as those used by Forés-Toribio et al. (2022).The strong lensing mass model parameters are reported in Table 3.
The measured time delay between images A and B (∆t AB = -38.4± 2.0 days) was first published in Fohlmeister et al. (2007).In this notation, a positive value of the time delay means image A leads the other image.In addition to reporting a refined value of ∆t AB = -40.6 ± 1.8 days, Fohlmeister et al. (2008) 4).

SDSS J1029+2623
SDSS J1029+2623 is a cluster at z = 0.588 that is strongly lensing a quasar at z = 2.1992 into three images (Inada et al. 2006;Oguri et al. 2008).The quasar images are in a naked cusp configuration with a maximum image separation of 22. ′′ 5 (Table 1).Acebron et al. (2022) reported spectroscopic redshifts of several galaxies in the field, based on Multi Unit Spectroscopic Explorer (MUSE) spectroscopy from the Very Large Telescope.They refined the redshift measurement of the quasar to z = 2.1992 (formerly reported as z = 2.197, Inada et al. (2006)).The other spectroscopically confirmed objects from MUSE include a doublyimaged galaxy at z =2.1812, a septuply-imaged galaxy at z =3.0275, a quadruply-imaged galaxy at z =3.0278, a doubly-imaged galaxy at z =1.0232, and a quadruplyimaged galaxy at z =5.0622 (Acebron et al. 2022) (Table 2).
We used archival HST multi-color imaging from GO-12195 (PI: Oguri): WFC3/F160W (2 orbits), ACS/F814W (3 orbits), and ACS/F475W (2 orbits).These data were originally proposed to identify multiply-imaged galaxies to construct a mass model that could be used to better understand the anomalous flux ratios between two of the quasar images and the dynam-ical state of the cluster (Oguri et al. 2013).These HST data also enabled a weak lensing analysis and a morphology study of the quasar host galaxy (Oguri et al. 2013).
Our lens model, which builds on the results from Acebron et al. ( 2022) and Oguri et al. (2013), contains 48 constraints and 33 free parameters.All of the model constraints are taken from Acebron et al. (2022).The model includes two cluster-scale dark matter halos that were allowed to vary in position around the two BCGs as well as two galaxy-scale halos that were fixed to the BCGs' positions.Additionally, a foreground galaxy (z =0.5111 from MUSE) and a background galaxy (z =0.6735 from MUSE) along the line of sight are both modeled at the cluster redshift since Lenstool does not yet implement a multi-plane lensing framework.This approach improves the accuracy of the lensing analysis outputs compared to omitting these interlopers from the model (Raney et al. 2020).The lens model includes 204 galaxy-scale halos.
Our lens model differs from Acebron et al. (2022) in the following ways.Whereas Acebron et al. (2022) include a model (Model 1) with an external shear component, we opted to not include this component as its inclusion does not significantly improve the rms scatter in the source plane.Additionally, for consistency with the other clusters modeled in this paper, our galaxy-scale halos have ellipticities, whereas Acebron et al. (2022) use spherical halos.We constructed our galaxy catalog as described in Section 3, taking into account the MUSE spectroscopy to determine the red sequence (see Sharon et al. 2022).We used the ACS F814W vs. F475W for selection.We identified the red sequence by fitting a line to the spectroscopic members in this phase space, with four iterations of sigma clipping.
We found that the source-plane minimization did a better job at predicting the quasar image positions in this cluster than the image-plane minimization, possibly due to the close proximity of quasar images B and C. Once a best-fit model was obtained, we examined the posterior distribution of image predictions and rejected from the MCMC sampling steps that did not produce this lensing configuration, i.e., not producing two separate images for A and B on either side of the critical curve.Since these two images lie very close to the critical curve, some parameter combinations produce solutions in which these two images merge and only image A of the quasar remains, in contrast to the observed lensing evidence.
The HST data and the lens model for SDSS J1029+2623 are shown in Figure 1.The strong lensing mass model parameters are reported in Table 5. Fohlmeister et al. (2013) published the time delay measurement between images A and B (∆t AB = 744 ± 10 days) based on photometric monitoring campaign at the FLWO 1.2m.
Spectroscopy of other lensed galaxies was obtained by the Gemini North Telescope.These data include triplyimaged and doubly-imaged knots from a galaxy at z = 4.562 and a doubly-imaged galaxy at z = 2.2963 (Sharon et al. 2017).
We used archival HST multi-color imaging from GO-13337 (PI: Sharon): WFC3/F160W, F110W (1 orbit) and ACS/F814W, F606W, F435W (6 orbits).These data were originally proposed to detect any additional quasar images and to compute a mass model (Sharon et al. 2017).Additionally, these HST data have enabled a spatially resolved study of the Lyman-alpha emission in the quasar host galaxy (Bayliss et al. 2017).
We adopted the lens model from Sharon et al. (2017) with 32 constraints and 31 free parameters.SDSS J2222+2745 is modeled using image plane minimization with one cluster-scale halo and 172 galaxyscale halos.Sharon et al. (2017) included as constraints triply-imaged and doubly-imaged knots at the quasar's redshift of z = 2.805, and triply-imaged and doublyimaged knots from a galaxy at z = 4.562.Two separate triply-imaged galaxies had their redshifts left as free pa-rameters, with priors of 2.0 ≤ z ≤ 4.0 and 3.8 ≤ z ≤ 5.0, respectively, based on photometric redshift analysis.The HST data and the lens model for SDSS J2222+2745 are shown in Figure 1.Table 5 of Sharon et al. (2017) lists the strong lensing mass model parameters.Dahle et al. (2015) first published the time delay measurements between images A and B (∆t AB = 47.7 ± 6.0 days) and A and C (∆t AC = -722 ± 24 days).Then Dyrland ( 2019) reported updated values for the time delays between images A and B (∆t AB = 42.44 +1.36  −1.44 days) and images A and C (∆t AC = -696.65+2.00 −2.10 days).These measurements were based on data from a monitoring campaign at the 2.5m Nordic Optical Telescope.
In the analysis that follows, we used the most up-to-date time delay values for SDSS J1004+4112, SDSS J1029+2623, and SDSS J2222+2745 which are listed in Table 4. Constraints on H0 from three clusterlensed quasars, SDSS J1004+4112, SDSS J1029+2623, and SDSS J2222+2745.The histograms are created from 500 random models sampled from the MCMC.Overplotted are Gaussian fits to the distributions.Whereas individual time delay measurements produce H0 values with an average of 41% error, the error is decreased to 11% when the systems are analyzed together.The inverse-variance weighted mean of H0 is 74.1 km s −1 Mpc −1 (solid gray line) and the standard error of the weighted mean is 8.0 km s −1 Mpc −1 .

RESULTS
Using the outputs of the lens models described in the previous sections, we computed the model-predicted time delay values for each of the quasar images in each cluster field with respect to image A of the quasar (Equation 3 and Table 6).
The time delay is a sensitive function of the positions of the source ( ⃗ β) and its multiple images ( ⃗ θ 1 , ⃗ θ 2 ).The unobservable source position and the locations of its multiple images are strongly coupled to the time delay, since stationary points in the arrival time surface determine the image-plane positions of multiple images of any given source-plane position (see Section 2).It is therefore important to measure time delays selfconsistently, by obtaining the time delay at the image positions predicted by the same lensing potential.Lens models are never perfect, and small scatter between the observed and predicted positions is expected.To maintain this self-consistency, we calculated the source position ⃗ β by ray-tracing the observed position of image A ( ⃗ θ A ) through the lens equation, and used the same lens model to predict the image-plane positions of its counter images ( ⃗ θ 2 , ⃗ θ 3 ,...).The time delay was then calculated from these predicted positions, rather than the observed positions, which may be slightly shifted from the stationary points in the Fermat potential.The scatter in the image or source plane contributes to the error budget through the MCMC exploration of the parameter space.An alternative approach to determining the source position is averaging the predicted source locations from all the quasar images, and calculating the predicted image locations of the average source.We found that this alternative approach of determining the source position produced an H 0 value consistent with the method described above, and hence, the choice of one approach vs. the other is not a significant source of systematic uncertainty.In what follows, we describe the results using the former method.
Using Equation 4, we computed the H 0 value corresponding to each independent published time delay value and corresponding predicted time delays.To generate the 1σ uncertainties in H 0 , we used 500 random models from the MCMC sampling of the parameter space for each cluster to ensure smooth sampling of the H 0 posteriors.
The number of measured time delays in each field determines the number of H 0 measurements derived from each cluster: three from SDSS J1004+4112, one from SDSS J1029+2623, and two from SDSS J2222+2745, for a total of six H 0 measurements.Table 7 lists the derived H 0 values and uncertainties, obtained for the "best" lens model, i.e., the one producing the smallest scatter, and for the full posterior distribution.We find that the fractional uncertainty on H 0 calculated using the lensing analysis is consistent with the analytical propagation of errors due to the source plane uncertainty, following Birrer & Treu (2019) Eq. 16 and the examples therein, giving δH 0 (δθ)/H 0 = 15 − 50%.
The resulting H 0 measurement from each quasar pair has large uncertainties due to the complexity of the lens and systematic uncertainties in the lens modeling process.Grillo et al. (2020) found that systematic effects such as adding a constant mass sheet at the cluster redshift and using a power-law profile for the mass density of the cluster's primary halo did not introduce a significant bias on the inferred H 0 value and that the statistical uncertainties dominate the total error budget.Given that SDSS J1004+4112, SDSS J1029+2623, and SDSS J2222+2745reside in the same universe, they all must have the same H 0 ; we can leverage these three independent lines of sight, with six time delays, to obtain a tighter constraint than what is possible from a single time delay, and marginalize over systematics related to line of sight effects such as intervening structure.We combine the results from the six time delays by taking the inverse-variance weighted mean of the six H 0 measurements, sampled from their posterior distributions.We accounted for the correlations between measurements made in the same cluster lens by using the same model realization from the MCMC to generate a linked set of H 0 values.For example, since SDSS J1004+4112 has three measured time delays, the same model realization was used to calculate the H 0 values from ∆t AB , ∆t AC , and ∆t AD .
We note that the observational time delay measurement uncertainties are negligible compared to the lens modeling uncertainties.The inverse-variance weighted mean and the standard error of the weighted mean of H 0 is 74.1 ± 8.0 km s −1 Mpc −1 (Fig. 2).Combining the H 0 values derived from multiple time delay values improves the constraints on H 0 , decreasing the uncertainty from ∼41% for an individual H 0 measurement to 11% for the sample.If SDSS J1004+4112's quasar image C is excluded from the analysis (see Section 4), we obtain H 0 = 77.4± 9.9 km s −1 Mpc −1 .This value is consistent with the H 0 measurement calculated with all six time delays.
Increasing the number of systems used for a combined time-delay measurement of H 0 will improve this method's competitiveness with CMB-based and distance ladder-based methods.Although four other clusterlensed quasars are published in the literature, none has all the necessary time delays measurements, spaceresolution imaging, and spectroscopic redshifts of secondary arcs for a measurement of H 0 .All four of the other published cluster-lensed quasars have ongoing photometric monitoring campaigns to measure their time delays.Additionally, one of the other four systems, COOL J0542-2125 (Martinez et al. 2022) will be observed by HST in Cycle 30 (GO-17243; PI: Napier).
To estimate the improvement in the H 0 constraint from a sample of twice as many time delay measurements, we simulated H 0 distributions guided by the existing sample, as follows.We sampled integer H 0 values from a Gaussian distribution centered at our combined H 0 value of 74.1 km s −1 Mpc −1 .We then randomly assigned to these six H 0 values the standard deviation of one of the six H 0 distributions (Table 7), and produced the corresponding Gaussian distributions.We repeated this simulation process 100,000 times.Incorporating these new six H 0 distributions for a total of 12 constraints, and averaging the 100,000 iterations, gave a standard error of the weighted mean of 5.5 km s −1 Mpc −1 .Therefore, doubling the number of systems results in a ∼35-40% improvement in the constraint on H 0 , reducing the uncertainty on H 0 from 11% to ∼7.4%.
A 1% uncertainty measurement of H 0 from clusterlensed quasars would be competitive with the current precision level of CMB and distance ladder methods.Extending the simulation described above to a larger number of systems, we estimated that 620 delay measurements from cluster-lensed quasars would achieve a 1% uncertainty level on H 0 from cluster lensed-quasars.Based on SDSS J1004+4112, SDSS J1029+2623, and SDSS J2222+2745 each having an average of two time delay measurements, a sample size of 310 cluster-lensed quasars would be needed to produce 620 time delay measurements.Future surveys are expected to detect of order ∼50 such systems in the next decade (Robertson et al. 2020).Therefore, this increase in sample size alone will not achieve 1% uncertainty in H 0 ; to reach 1% with of order of 50 systems (100 time delays) will require a decrease in the lens modeling uncertainties by about a factor of two and a half, on average.Future work will explore whether this decrease in the uncertainties is feasible.

ACKNOWLEDGMENTS
Based on observations made with the NASA/ESA Hubble Space Telescope, obtained from the Multimission Archive at the Space Telescope Science Institute (MAST) at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS 5-26555.These archival observations are associated with programs GO-10509, GO-9744, GO-10793, GO-12195, and GO-13337.Support for HST program AR-16150, which enabled this work, was provided through grants from the STScI under NASA contract NAS5-26555.Co-author GM received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 896778.We thank Ana Acebron for her useful discussions about SDSS J1029+2623.We thank the anonymous referee for their constructive feedback and comments that improved the manuscript.

Figure 1 .
Figure1.Hubble Space Telescope imaging of the three cluster-lensed quasars used to derive H0.We computed the lens models for SDSS J1004+4112 and SDSS J1029+2623.SDSS J2222+2745 is reproduced fromSharon et al. (2017).The positions of the quasar images are denoted with the cyan letters.The critical curves, the loci of maximum magnification at a specified source redshift, are generated at the quasar redshifts -z = 1.734, z = 2.1992, and z = 2.805, for SDSS J1004+4112, SDSS J1029+2623, and SDSS J2222+2745, respectively, and are plotted in red.±2.1 days) and set a lower limit of ∆t AD > 1250 days.After the completion of a 14.5 year monitoring campaign at the 1.2m Fred Lawrence Whipple Observatory (FLWO),Muñoz et al. (2022) presented new light curves for the four brightest images in SDSS J1004+4112, resulting in updated time delay values of ∆t AB = -43.01± 0.27, ∆t AC = -825.23± 0.46 days, and ∆t AD = 1633.23± 0.97 days (Table4).
Figure 2.Constraints on H0 from three clusterlensed quasars, SDSS J1004+4112, SDSS J1029+2623, and SDSS J2222+2745.The histograms are created from 500 random models sampled from the MCMC.Overplotted are Gaussian fits to the distributions.Whereas individual time delay measurements produce H0 values with an average of 41% error, the error is decreased to 11% when the systems are analyzed together.The inverse-variance weighted mean of H0 is 74.1 km s −1 Mpc −1 (solid gray line) and the standard error of the weighted mean is 8.0 km s −1 Mpc −1 .

Table 1 .
The quasar image positions and redshifts.Also included are the model-predicted median magnifications at the observed positions of the quasar images.