Parameter free Hubble constant from the quadruply lensed quasar SDSS J 1004 + 4112

,


Introduction
The SDSS J1004+4112 galaxy cluster is a gravitational lens at a redshift of z = 0.68, which produces five images of a background quasar at redshift z = 1.734.This gravitationally lensed quasar was found in the Sloan Digital Sky Survey by Inada et al. (2003), when the large separation between its multiple images was noticed.Following the discovery of SDSS J1029+2623 (Inada et al. 2006), currently SDSS J1004+4112 is the second lensed quasar with the greatest reported separation.
There is a broad range of previously published works based on this lens system and many of them have been focused on observational constraints, for instance, cluster substructure (Oguri et al. 2004), central quasar image (faintest image E; Naohisa et al. 2005), background lensed sources (Sharon et al. 2005;Liesenborgs et al. 2009;Oguri 2010), and time delays between the four brightest images (ABCD) of the quasar (Fohlmeister et al. 2007(Fohlmeister et al. , 2008;;Muñoz et al. 2022).
Most recently, following the method proposed by Refsdal (1964), Napier et al. (2023) used the Lenstool software to measure the Hubble constant from six measured time delays (including three independent delays of SDSS J1004+4112) and parametric mass models.They obtained H 0 = (71.5 ± 61) cm s −1 Mpc −1 .In addition, Liu et al. (2023) focused on SDSS J1004+4112, using the GLAFIC software to estimate H 0 from its measured delays and 16 parametric models with equal weight.The resulting Hubble constant is H 0 = 67.5 +14.5  −8.9 km s −1 Mpc −1 .
In this work, we explore Refsdal's method in the lensed quasar with the most precise time delays among those with a large number of observation constraints for the lensing mass and image separation exceeding 10/arcsec.A large number of observation constraints are required to successfully reconstruct the lensing mass from a non-parametric method, so galaxyscale lensed quasars are usually studied from parametric models.At present, there are no non-parametric lens mass reconstructions that use the three well-measured time delays of SDSS J1004+4112 to estimate the Hubble constant.Thus, the current work complements the two recent studies of the system (see: Liu et al. 2023;Napier et al. 2023) to achieve a more comprehensive perspective.
For the lens mass reconstructions, the redshifts and image positions of seven lensed galaxies were considered, taken from Oguri (2010).The lensed sources are found at four different redshifts, which means that we do not need to be concerned about the mass-sheet degeneracy (Diego et al. 2007).In the case of sources at different redshifts, different models tend to find similar mass profiles that are not biased, demonstrating that the degeneracy does not exist for multiple lensed sources (Meneghetti et al. 2017).With regard to the measured time delays, as in the work of Forés-Toribio et al. (2022), their formal errors (see Muñoz et al. 2022) are enlarged by a factor of 5. We also remark that quasar image flux ratios at near-infrared and radio wavelengths (which are less sensitive to microlensing effects) are used to weight mass models.The method used for this cluster can be applied to future data with similar measurements.
The paper is organized as follows: Sect. 2 gives a brief description of the gravitational lensing theory, proposing the mathematical formulation of the problem and briefly explains the WSLAP+ strong lensing inversion code (Diego et al. 2005(Diego et al. , 2007) ) used to derive the free-form lens model.Section 3 presents 21 different reconstructed lens models, used to account for the uncertainty in the predicted time delays from the lens model.In Sect.4, the likelihood functions of H 0 and the weight of each model are discussed.Finally, in Sect.5, the obtained estimation of H 0 is shown, followed by a discussion of the results in Sect.6.
In this work, we assume a cosmology with parameters Ω M = 0.3 and Ω Λ = 0.7.The mass models derived in this work scale as 1/H 0 , which depends very weakly on the values of Ω M and Ω Λ within the currently accepted ranges of these parameters.The time delay also scales as 1/H 0 which is by far the most relevant cosmological dependence.

Lens modeling
The reconstructed lens models of SDSS J1004+4112 described in this article are obtained with the free-form code WSLAP+ (Diego et al. 2005(Diego et al. , 2007)).Details are given in the references.Here, we simply give a brief description of the algorithm.The lensing problem in its most basic form is formulated linearly as with θ corresponding to the observed positions of the lensed images, β the unknown positions of the background galaxies, M(θ) the unknown mass distribution of the lens and α(θ, M(θ) the deflection angle that corresponds to the observed position θ.
The deflection angle is the contribution to the deflection of all the mass distribution, that is obtained by integrating where D ls , D l , and D s are the angular distances from the lens to the source galaxy, from the observer to the lens, and from the observer to the source correspondingly.
The problem is discretized assuming the lens is divided into N c cells, in each of which the mass is nearly constant, assigning a mass m i to each cell.The deflection angle Eq. ( 2) is then approximated by a sum this is, the deflection is approximated by the contribution of N c discrete masses m i (contained in the array of masses M), in positions θ i in the direction of θ − θ i and magnitude m i /|θ − θ i |.The matrix γ hence contains the deflection field at position θ j from a grid point at position θ i with mass m i = 1.The galaxy cluster SDSS J1004+4112 shows a distribution of tangential strong lensing arcs over the image.Discretizing these lensed arc positions into N θ points, the lens equation can be formulated as a linear problem between unknown positions β and the masses m i where β and θ are 2N θ vectors containing the x and y unknown positions in the source plane and the image positions in the image plane, while M is the mass vector containing the N c mass cells.γ is a 2N θ × N c whose i-th row γ i gives the linear relation between the mass distribution contained in the vector M and the deflection angle in the corresponding θ i position contained in the i-th row of the θ vector.The mass distribution can also contain a contribution from a compact component associated with member galaxies.This component can be divided into different layers that from the point of view of Eq. ( 4) behave as additional grid points.The number of layers must be at least N l = 1 if all galaxies are placed in the same layer and up to N l = N gal if we assign one layer per galaxy.In this work, we adopt one single layer for all galaxies, that is, N l = 1.
The system has 2N s + N c + N l unknowns (lens masses and source positions), that can be reduced assuming the N s sources can be approximated by point sources, this is, forcing the θ positions of the arcs correspond to the N s source positions.This can be formalized as where Γ is a 2N θ ×(N c +N l +2N s ) known matrix given by the grid configuaration and positions of lensed galaxies and M contains all the unknowns, namely, the mass elements and positions of the lensed galaxies in the source plane.
As the problem may be ill-conditioned and too big to use matrix inversions, approximate numerical methods are used, for instance, the bi-conjugate gradient method (Press et al. 1997) can produce a fast solution to the system of equations given by Eq. ( 5).However, models obtained with this approach often return unphysical solutions where elements in the mass vector, M, can adopt negative values.This is solved by imposing the constraint that the mass must be always positive.This type of problem can be solved by somewhat slower, but still very fast quadratic programming algorithms (see Diego et al. 2005).
A first solution is obtained where the lens plane is initially discretized with a regular 16 × 16 grid.This is similar to assuming no prior knowledge about the distribution of the mass since all grid points play the same role in the system of equations given by Eq. ( 5).After this first solution, a multi-resolution grid is performed, which is based on the previously obtained lens model.This new grid increases the number of grid points in the areas where the matter concentration of the original solution is greater by making the cells' size inversely proportional to the matter density.Figure 1 shows a multi-resolution grid with 199 points, which concentrates more grid points in the center as a higher mass density is expected.

Lens model results
The uncertainty in the derived Hubble constant is dominated by the uncertainty in the lens model.In order to capture this uncertainty, we derived a range of models, all of them consistent with the observed set of linear constraints (positions of lensed galaxies and QSO).The observed time delays or relative magnifications, are not included in our set of constraints.In the case of WSLAP+, a major source of uncertainty in the lens model originates in the particular choice of the grid used to describe the smooth component of the mass.The grid configuration is a free variable in WSLAP+.For the number of constraints used in this work, grid configurations with 150-300 points can produce satisfactory results in terms of reproducing the position of lensed galaxies with relatively smooth critical curves.A larger number of grid points results in unstable and unreliable solutions and smaller numbers result in a lack of resolution.The cases presented should capture the expected range of uncertainty.
Another source of variability comes from the starting point of the minimization.It is possible to choose among an infinite number of initial conditions.After each minimization, each solution differs slightly from other solutions obtained with a different seed.We varied both the grid and initial seed to explore the range of valid models.

Range of lens models
We derived a total of 21 lens models: one with a uniform grid of 16×16 points, nine a dynamical grid of 199 points, six with a different dynamical grid of 248 points, and an additional five models with a different dynamical grid of 299 points.However, only the 20 multiresolution grids were used to compute the estimation of the Hubble constant since the solution obtained with the regular 16 × 16 grid lacks resolution in the central region.
Figure 2 shows the convergence of 10 of the 21 models used, including the initial model obtained with the uniform grid 16 × 16 = 256 grid points.The first number in the legend indicates the number of grid points, while Model1, Model2, and Model3 are different realizations (different initial conditions) for the same grid configuration.Models that are not shown have similar profiles.These profiles are consistent with the ones derived by Liu et al. (2023) using a parametric model.Despite the relative similarity between profiles, the small changes between models are the main source of uncertainty in the predicted time delay, and hence in the derived value for the Hubble constant.Figure 3 shows the magnification map of the lens model, obtained with the dynamic grid of 199 points in Fig. 1 plotted in red on top of the galaxy cluster shown in blue.

Predicted time-delay for image E
Each of the lens models described in the previous subsection gives a prediction for the image E time-delay for a canonical Hubble constant H 0 = 70 km s −1 Mpc −1 .Thus, by re-scaling the predicted time delays to the value of the Hubble constant measured, H 0 = 74 km s −1 Mpc −1 (see Sect. 5), we can obtain a predicted time-delay with respect to image C of 3200 ± 200 days.This is approximately 1 yr longer than the expected value reported in Forés-Toribio et al. ( 2022), who obtained a predicted time delay of about 8 yr.
A robust determination of the time delay of image E (with respect to image C) would play an important role since this measurement would significantly improve the precision in the estimation of the Hubble constant.However, image E is very  faint.This means that the typical magnification ratio between images D and E is µ D /µ E ∼ 30 (with a great variation among different models) and, consequently, it is reasonable to assume that image E is about 30 times fainter than image D at optical wavelengths (neglecting possible extinction and microlensing effects).Hence, if the r-band magnitude of D is ∼21 (see Table 1 of Muñoz et al. 2022), the expected r-band magnitude of D would be ∼24.7 mag.Additionally, E is located in a bright and crowded sky region, so it seems extremely difficult to build its optical light curve from current facilities.

Constraints on H 0
It is common to express the Hubble constant H 0 as: where h = H 0 /100 is referred to as the reduced Hubble constant and it is a dimensionless physical constant.
A187, page 3 of 7 In order to measure H 0 , Refsdal (1964) proposed that the gravitational lensing effect could be used.In fact, when a light ray is deflected by a gravitational lens, its trajectory is modified and therefore its traveling time, due to the light's finite propagation speed.The expression that describes the time delay of the perturbed ray with respect to the one that follows a straight path is displayed in Eq. ( 6): with z l indicating the redshift of the lens.
The three angular distances involved in Eq. ( 6) are inversely proportional to the Hubble constant, H 0 , therefore, we say that the time delay scales as 1/h.In addition, the magnitude expressed in Eq. ( 6) is not measurable; instead, we chose to measure the time delay relative to one of the images.
Using the lens models derived in the previous section, we are able to confront the predictions of the time delay of the multiple images of the QSO with the observations, both relative to image C. Our lens models predict the time delay for a canonical value of H 0 = 70 km s −1 Mpc −1 , which is a reduced Hubble constant of h = 0.7.Hence, these time delays are re-scaled to a new reduced Hubble constant (h ) by multiplying the predicted time delays, x = h /0.7

Likelihood
For each lens model, we can derive a probability for the presence of h by minimizing the weighted difference between the predicted and observed time delays (∆T M i and ∆T obs , respectively): where the index j runs over the four positions, A, B, C, and D. In addition, the proportionality sign just expresses the lack of a normalization constant.The values of σ j are computed at each of these four positions as the dispersions of the time delays predicted by the lens models.Image E has no observable time delay and hence is not included in the definition of the likelihood.
Since we have different lens models derived under different assumptions, we need to compute the joint probability from all models.We adopted the conservative approach that the combined probability from all lens models is the sum of probabilities of the models (see, e.g., Vega-Ferrero et al. 2018).A more aggressive approach would be to consider that the joint probability is the product of the individual probabilities, but this implies all models are independent from each other.This is not true since all lens models are derived using the same algorithm, with some assumptions common to several models (e.g., the grid configuration, definition of the compact component, or lensing constraints).
The joint probability is then defined as: where the sum runs over the N lens models, p(M i , h) is the probability given in Eq. ( 7) and the term w i is the weight assigned to each model.Once again, the proportionality sign just indicates the lack of a normalization constant.This weight is determined by how well a particular model reproduces the observed flux ratios between the different images of the QSO.We define this weight as (dropping the subscript i); In the above definition, we have ignored image E since this is very close to the center of the BCG galaxy and extremely faint.In addition, initially, we did not include the flux ratio, F A /F C , in w Eq. ( 9) because the A image is very close to a critical line and, thus, the magnification ratio, µ A /µ C , in some models might be a strongly biased estimator of F A /F C .Hence, in an abundance of caution, image A was ignored for respect to the definition of the weight.We computed the flux ratios as the average of three flux ratios from our own measurements in the IR band F814W, measurements published in Oguri (2010) and radio flux in Hartley et al. (2021).The values of σ CB , and σ CD are the combination of the observed dispersion of flux ratios between different measurements and the expected variation due to the intrinsic variability: For the intrinsic variability in the F814W HST filter, we adopted a value of σ int = 0.2 for image B and σ int = 0.16 for image D from Li et al. (2018) and Morganson et al. (2014).The measured time delays, flux ratios measured by the indicated authors, and estimated intrinsic variability are summarized in Table 1.The flux ratios from Oguri and this work correspond to the F814W HST filter and Hartley's in the 5 GHz radio band.

H 0 results
Each of the models obtained in this work has a different contribution to the prediction of the value of the Hubble constant.
Figure 4 shows the contribution of each of the models to the final likelihood, this is, it shows the function defined by each weighted probability in the likelihood function in Eq. (8).
Our main results are incorporated in Fig. 5.In this figure, we show the probability for the reduced Hubble constant using three different schemes.The standard scheme relies on Eqs. ( 7)-( 9) and the 20 lens models, leading to the probability distribution shown in blue.Selecting the 11 lens models with the highest resolutions (248 and 299 grid points), we were also able to obtain the distribution in red.Additionally, we considered the 20 lens models and a modified version of Eq. ( 9), including flux A187, page 4 of 7

Conclusions
These results are in agreement with recent estimates reported in Liu et al. (2023) and Napier et al. (2023).The present paper, along with these two works, show that the Hubble constant must lie in the 60−80 km s −1 Mpc −1 range.The uncertainty in H 0 can be reduced in the future with deeper observations of this cluster, which would reveal more lensed galaxies and their corresponding redshifts.This would take advantage of lensing by galaxy clusters, which typically can have over an order of magnitude more lensing constraints than classic galaxy-QSO lensing.Also, future observations of lensed SNe in galaxy clusters will provide competitive constraints on H 0 , similar to those obtained with SN Refsdal (Vega-Ferrero et al. 2018;Kelly et al. 2023), or SN H0pe (Frye, in prep.).
To conclude, we analyze the contribution of this work to the Hubble tension problem.The Hubble tension refers to the 4σ to 6σ discrepancy between estimations of H 0 from early universe measurements and late universe results.The most precise predictions of the Hubble constant using early universe results come from the CMB Plank Collaboration (Planck Collaboration XII 2020) that obtained H 0 = 67.49± 0.53 km s −1 Mpc −1 with 68% confidence level (CL).In the late universe, however, the most precise results come from different techniques: using distance ladder Cepheid calibration in the Large Magellanic, Riess et al. ( 2019) estimated H 0 = 74.03± 1.42 km s −1 Mpc −1 ; the SNe Ia distance calibration by the SH0ES collaboration (Riess et al. 2022) that obtained H 0 = 73.30± 1.04 km s −1 Mpc −1 (68% CL); and gravitational lensing on quasars by H0LiCOW collaboration (Wong et al. 2019) that obtained H 0 = 73.3+1.7 −1.8 km s −1 Mpc −1 (68% CL) and TDCOSMO with H 0 = 77.1 +7.3 −7.1 km s −1 Mpc −1 (Shajib et al. 2023) in their latest result.There is a chance that the tension might be fictitious (e.g., the resulting effects of the selection of the SNe Ia or use of inadequate lens mass models) or real, in which case a modified cosmological model would be required (e.g., Dainotti et al. 2021;Goicoechea & Shalyapin 2023;Shalyapin et al. 2023).Nonetheless, our results do not show any tension with any of the results in either the early universe or late universe, although this is largely due to the lens model uncertainties.

Appendix A: Additional estimated parameters
Table A.1 includes the convergence, κ; shear, γ; and shear direction, φ, for the quasar's five images for the 20 lens models used in this work for the estimation of the Hubble constant, together with the mean value for these in the last row of the Table.

Fig. 1 .
Fig. 1.Example of a dynamic grid with 199 points used for the smooth dark matter distribution.

Fig. 3 .
Fig. 3. Magnification map of the reconstructed lens of SDSS J1004+4112 using the dynamic grid with 199 points.

Fig. 4 .
Fig. 4. Contribution of the seven models with the highest weights to the final likelihood.

Fig. 5 .
Fig. 5. Probability distribution of h obtained from time-delay predicted by the lens models using different analyses.

Table 1 .
Observed time delay, flux ratios, and adopted intrinsic flux ratio variation for the QSO images of SDSS J1004.

Table A .
1. Convergence κ, shear γ and shear direction φ of the multi-resolution lens models for the different QSO images.