Determining cosmological-model-independent $H_0$ and post-Newtonian parameter with time-delay lenses and supernovae

Strong gravitational lensing provides a natural opportunity to test General Relativity (GR). We propose a model-independent method for simultaneous constraining on Hubble constant ($H_0$) and post-Newtonian parameter (${\gamma_{\rm{PPN}}}$) using strong lensing systems and observational SNe Ia. The time-delay measurements from strong lesning can directly determine the Hubble constant, and the lens distance inferred from the spectroscopic measurement of the stellar kinematics of the deflector galaxy can help us to constrain the post-Newtonian parameter. We seek the Pantheon dataset and reconstruct unanchored distances using Gaussian process regression to achieve the cosmological model-independent GR testing instead of assuming a specific model, which can reduce possible bias on GR testing and measurement of Hubble constant. Combining the reconstructed unanchored distances and the four H0LiCOW lenses datasets, our results are $H_0=72.9^{+2.0}_{-2.3} {\mathrm{~km~s^{-1}~Mpc^{-1}}}$ and ${\gamma_{\rm{PPN}}}=0.89^{+0.17}_{-0.15}$. All the lenses show that there is no obvious evidence to support GR deviation within observational uncertainties. In the subsequent analysis, we consider a ratio of distance ${D_{\Delta t}}/{D^{'}_{d}}$ method to further avoid the influence of $H_0$ on GR testing. The results show that, except J1206 within the $\sim1.2\sigma$ observational uncertainty, the remaining 3 lenses support GR holds within the $1\sigma$ observational uncertainties.


INTRODUCTION
The modern theory of cosmology is based on two pillars, Einstein's theory of General Relativity (GR) and the cosmological principle.The former was the first to equate the gravitational field with the curvature of space-time and was extremely successful in describing the gravitational interaction between matter.The latter describes that our Universe is homogeneous and isotropic on large scales.The Lambda cold dark matter (ΛCDM) model is consistent with most popular observational evidences, such as the observations of type Ia supernovae (SNe Ia) (Riess et al. 2007), Cosmic Microwave Background Radiation (CMB) (Planck Collaboration et al. 2016), and regarded as the standard cosmological model.However, there is a recent tension of 4.4σ or more between the Hubble constant inferred by CMB observations under the assumption of a flat ΛCDM model (Planck Collaboration et al. 2020) and its value measured through Cepheid-calibrated distance ladder by the Su-⋆ E-mail:liaokai@whu.edu.cnpernova H0 for the Equation of State collaboration (SH0ES) (Riess et al. 2019).Such inconsistency could be caused by unknown systematic errors in astrophysical observations or reveal new physics significantly different from ΛCDM model.In a recent work, Brout et al. (2022) reanalyzed and recalibrated the photometric systems in the Pantheon+sample of SN Ia, including those in the SH0ES distance-ladder measurement of H0, and suggested that supernova calibration is not currently capable of resolving Hubble tension.
As an important prediction of GR, gravitational lensing is a powerful tool to study the velocity dispersion function of early-type galaxies (Matsumoto & Futamase 2008;Geng et al. 2021;Chae 2007), the distribution of dark matter (Cao et al. 2022;Ruff et al. 2011;Mellier et al. 1993;Newman et al. 2009), and cosmological parameters (Suyu et al. 2014;Bonvin et al. 2017;Liu et al. 2019;Chen et al. 2019b).More importantly, the time-delay measurements between multiple images from strong gravitational lensin provide a valuable opportunity for determination of H0 (Refsdal 1964).In the milestone work of Wong et al. (2020), the six gravitationally lensed quasars with well-measured time delays were jointly used to constrain the H0 by the H0 Lenses in COSMOGRAIL's Wellspring (H0LiCOW) collaboration.Assuming a flat ΛCDM model, the H0LiCOW collaboration reported the results on H0 with these six lensed quasars, H0 = 73.3+1.7 −1.8 km s −1 Mpc −1 , which is consistent with local measurements from SN Ia, but in 3.1σ tension with CMB observations.However, it needs to be stressed that both the Planck and H0LiCOW inferred H0 values are based on GR plus ΛCDM model, that inspires us to investigate the validity of GR with cosmological model-independent way.The validity of GR can be verified by constraining the parametrized post-Newtonian (PPN) parameter γPPN (since GR predicts exactly γPPN ≡ 1), which describes the spatial curvature generated by an unit rest mass.In recent years, especially on the solar system scale, many tests on GR have made great achievements in extremely high precision (see review (Will 2014) for more works about testing GR).However, testing GR at extra-galactic scale is still not very precise.For instance, there is only ∼ 20% precision on γPPN at 10 − 100 Mpc scales by using the joint measurements of weak gravitational lensing and redshift-space distortions (Simpson et al. 2013;Blake et al. 2016).At cosmological scales, strong gravitational lensing systems provide an effective way to probe deviation of GR (Bolton et al. 2006;Cao et al. 2017;Liu et al. 2022;Wei et al. 2022).Recently, the work by Collett et al. (2018) used a nearby strong gravitational lens ESO 325-G004 to test GR and reported the constrained results on γPPN = 0.97 ± 0.09 at 1σ confidence level (C.L).In further research, Yang et al. (2020) proposed a new methodology through time-delay measurements combined with the stellar kinematics of the deflector lens from strong lensing to simultaneously constrain H0 and γPPN, and showed that there is no obvious deviation from the GR with the result γPPN = 0.87 +0.19  −0.17 .However, it should be emphasized here that these work are cosmological model-dependent (in ΛCDM model).The testing GR should be done without invoking any particular background cosmology model in order to reduce potential bias from the forms of the parametric or model assumptions.
Inspired by above, we propose a cosmological modelindependent method to constrain the Hubble constant and PPN parameter simultaneously at cosmological scales using four strongly lensed quasars published by H0LiCOW with both time delay distance and lens distance.For the distance information required to constrain the PPN parameter, we seek for unanchored (or relative) distances from SN Ia observations using a Gaussian process (GP) regression.Intuitively, the Time-Delay Strong Lensing (TDSL) measurements can directly determine the Hubble constant, and the lens distance inferred from the spectroscopic measurement of the stellar kinematics of the deflector galaxy can help us to constrain the post-Newtonian parameter.This paper is organized as follows: In Section 2 we introduce the methodology and H0LiCOW lensing data including reconstructed H0D L using GP regression.The constrained results on H0 and γPPN and discussions are given in Section 3. We conclusion our result in Section 4. The natural units of c = G = 1 are adopted throughout this work.

Distances inferred from H0LiCOW program
In the limit of a weak gravitational field, the Schwarzschild line element of space-time can be expressed as where Ψ is the Newtonian potential and Φ represents the spatial curvature potential, and Ω is the angle in the invariant orbital plane.The ratio γPPN = Φ/Ψ denotes as the PPN parameter, which describes the spatial curvature generated per unit mass.It should be emphasized that the PPN parameter γPPN is predicted to be one or Ψ = Φ.In this work, we assume that γPPN is a constant at the lens galaxy scales.
The crucial idea of using strong lensing systems to test gravity is through two different mass measurements, i.e., gravitational mass inferred from the lensing image, and dynamical mass obtained from the spectroscopic measurement of stellar kinematics of the deflector galaxy.The motion of non-relativistic matters (usually made up of baryonic matter and dark matter) is governed by the Newtonian gravitational potential Ψ, which obeys Poisson equation.However, the motion of relativistic particles is sensitive to both potentials.Testing γPPN requires observing the motion of relativistic and non-relativistic particles around the same massive object.Thus, strong lensing systems provide a natural laboratory to test gravity and further measure the PPN parameter γPPN.The difference between dynamical mass and lensing mass can be directly used for difference between Ψ and Ψ+ = Ψ+Φ 2 = ( 1+γ PPN

2
)αGR(θ), and effective lensing potential (the integral of the Weyl potential along the line-of-sight) rescales as ψ+ = ( 1+γ PPN 2 )ψ, as well as convergence field κ ′ = ( 1+γ PPN 2 )κ.It is worth noting that the deflector angle is directly related to the cosmological distance, so the PPN parameter and the cosmological distance are highly degenerated.This is one of the main limitations of using strong lensing systems to test gravity.To break this degeneracy, additional data needs to be taken into account, either the cosmological or the gravitational one.The time delay measurements are able to break the degeneracy alone.
For a given strong lensing system, quasar acting as background source, time delays between multiple images can be measured from variable AGN light curves, and determined by both the geometry of the Universe as well as the gravitational potential of the lensing galaxy through (Shapiro 1964) where φ(θ, β) = (θ − β)2 /2 − ψ(θ) is the Fermat potential at images, β is the source position, ξ lens is the lens model parameter.The cosmological background is reflected in socalled "time delay distance" D∆t = (1+z d )D d Ds/D ds , which is inverse proportional to H0.The key here is the Fermat potential difference ∆φAB(ξ lens ), which can be reconstructed by high-resolution lensing imaging from space telescopes.As we mentioned above, in the PPN framework, the inferred mass parameters are rescaled by a factor of (1 + γPPN)/2.Therefore, we denote the actually inferred lens model parameters in the Fermat potential as ξ ′ lens .We rewrite the time-delay distance as .
(3) This is the first distance we need.It can be obtained from both the measurements of time delay and the Fermat potential reconstructed with parameter ξ ′ lens .On the other hand, the stellar kinematics of lensing galaxies are only sensitive to the Newtonian potential Ψ, which are independent of PPN parameters.It can also be obtained from the modeling of the kinematic observable in lensing galaxies.The modeling of the stellar kinematic in lensing galaxies by the H0LiCOW collaboration is quite mature, and here we give a brief introduction, thus provide the reader with a detailed background.The dynamics of stars with the luminosity density distribution of lenses ρ * (r) in the gravitational potential Ψ follows the Jeans equation.In the limit of a relaxed (vanishing time derivatives) and spherically symmetric system, with the only distinction between radial (σr) and tangential (σt) dispersions, the anisotropic Jeans equation is where βani(r) ≡ 1 − is the stellar distribution anisotropy.The luminosity-weighted projected velocity dispersion σs is Suyu et al. 2010), where I(R) is the projected light distribution and R is the projected radius.Considering observational conditions, the luminosity-weighted line-of-sight velocity dispersion σv within an aperture A is real observations, which is given by σ A [I(R) * P]dA , where P is point spread function (PSF) convolution of the seeing.The prediction of the stellar kinematics requires a three-dimensional stellar density ρ * and mass M (r) profile.In terms of imaging data, the information about the parameters of the lens mass surface density with parameters mass ξ lens and the surface brightness of the lens with parameters ξ light can be extract.Finally, the prediction of any σv from any model can be decomposed into cosmological part Ds/D ds and stellar kinematics part J(ξ lens , ξ light , βani) (Birrer et al. 2019).The function J captures all of the model components calculated from the sky angle (from the imaging data) and the anisotropy distribution of the stellar orbit (from the spectroscopy).This allows us to obtain the distance ratio Ds/D ds from the wellmeasured velocity dispersion, independent of the cosmological model and time delays, but still relies on the lens model ξ lens (not the ξ ′ lens under PPN) (Birrer et al. 2016(Birrer et al. , 2019) The lens model parameter in J is the "unrescaled" ξ lens .
Here, we use ξ ′ lens to replace ξ lens , the corresponding distance ratio shall also be rescaled, as  3) and ( 6), we obtain (Birrer et al. 2016(Birrer et al. , 2019) This is the second distance we need.More details for these two distance obtained from strong lensing systems refer to work (Wong et al. 2020) and references therein.Thanks to the H0LiCOW collaboration, four lenses (namely RXJ1131-1231 (Suyu et al. 2013(Suyu et al. , 2014)), PG1115+080 (Chen et al. 2019a), B1608+6561 (Suyu et al. 2010;Jee et al. 2019), J1206+4332 (Birrer et al. 2019)) posterior distributions with both angular diameter distances of lens D d and time delay distances D∆t have been published, making it easier to use them.There are available on the website of H0LiCOW 2 .The redshifts of both lens and source, the time delay distances, and the angular diameter distance to the lenses for these lensed quasar systems are summarized in Table 2 of Wong et al. (2020).For more relevant work by using these lensing systems, we refer the reader to see the literature (Ding et al. 2021;Liu et al. 2022;Bag et al. 2022;Sonnenfeld 2021;Liao et al. 2015;Rathna Kumar et al. 2015;Liao et al. 2020;Liu et al. 2022).
The next problem is the acquisition of these two distance information.However, since D d depend on the cosmological model, D∆t also change in different cosmological models.The traditional method is to assume a particular cosmological model, such as a flat ΛCDM model, but it is cosmological model-dependent.In this work, we follow the reconstructed method used in work (Liao et al. 2019), and seek for current SN Ia observations to determine distances, even though the observations of SN Ia are anchored relative distances.

Unanchored distance from observations of SNe Ia using GP regression
As the most explosive variable source in Universe, due to the nature of SNe Ia (as standard candles), they were regarded as powerful cosmological probes.It was the observations of SNe Ia that led to the discovery of the accelerating expansion of Universe.We use recent Pantheon dataset, consisting of 1048 SNe Ia spanning the redshift range 0.01 < z < 2.3 (Scolnic et al. 2018).To combine the Pantheon SNe and the H0LiCOW strong lenses datasets, we generate samples of the unanchored luminosity distance H0D L from the posterior of the Pantheon dataset.In order to perform the posterior sampling in a cosmological model-independent way, the GP regression (Holsclaw et al. 2010,b;Joudaki et al. 2018;Shafieloo et al. 2012Shafieloo et al. , 2013;;Keeley et al. 2019) is considered here by using the code GPHist3 (Kirkby & Keeley 2017).GP regression is powerful tool for reconstruction of function since the regression occurs in an infinite-dimensional function space without overfitting problem (Joudaki et al. 2018).GP regression works by generating a large sample of functions γ(z) determined by the covariance function.The covariance between these functions can be described by a kernel function.We adopt a squared-exponential kernel to parameterize the covariance with hyperparameters σ f and ℓ that are marginalized over.The γ(z) is a random function inferred from the distribution defined by the covariance, and we adopt γ(z) = ln([H fid (z)/H0]/[H(z))/H0]) to generate expansion histories H(z)/H0 by using Pantheon dataset.Here H fid (z)/H0 is chosen to be the best fit ΛCDM model for the Pantheon data and serves the role of the mean function for GP regression.Since the final reconstruction result is not completely independent of the mean function, it has some influence on the final reconstruction result because the value of the hyperparameter helps to track the deviation from the mean function.And the true model should be very close to flat ΛCDM, hence, our choice for the mean function is very reasonable (Shafieloo et al. 2012(Shafieloo et al. , 2013;;Aghamousa et al. 2017).
To highlight the purpose of our work, i.e., the Hubble constant and γPPN, we assume a flat universe here.With the reconstructed expansion history H(z)/H0, the unanchored SN luminosity distances can be calculated The 1000 unanchored luminosity distance curves H0D L (z) reconstructed from the SN data are shown in Fig. 1.It shows the shape of the cosmic distance-redshift relation of Pantheon data very well.It should be noted that the redshift of SN dataset well covers range of four strong lensing system redshifts, so the redshift range needn't extrapolation.

Simultaneous measurements on Hubble constant and PPN parameters
In summary, combining the three measurements, namely optical lensing image, deflector spectroscopies as well as time delays, two distances (D∆t, D ′ d ) can be inferred, simultaneously, based on Eqs. ( 3) and ( 7).It should be stressed that there is only deflector galaxy distance (D ′ d ) carried the information of the PPN parameter, while time delay distance (D∆t) is sensitive to H0.Thus, the combination of D∆t and d directly provides a new way for simultaneous measurement of H0 and γPPN.
The steps for simultaneous constraining H0 and γPPN are summarized as follows: (i) Draw 1000 unanchored luminosity distances curves H0D L from SN data, and convert to unanchored angular diameter distances H0D A for serving the strong lensing systems.
(ii) Calculate 1000 values of H0D d at the lens from the 1000 H0D A curves for serving four posterior distributions of lens from H0LiCOW program; Adopt same procedure for source redshifts of the four strong lens systems to calculate 1000 values of H0Ds; Combine these H0D d and H0Ds to calculate H0D∆t for each system using H0D∆t = (1 + z d )(H0D d )(H0Ds)/(H0D ds )4 ; (iii) Compute the likelihood, for each of the 1000 realizations, from the H0LiCOW's D∆t and D d data for each lens system for many values of H0 and γPPN; (iv) Multiply the four likelihoods to form the full likelihood for each realization, for each values of H0 and γPPN; (v) Marginalize over the realizations to form the posterior distributions of H0 and γPPN.

RESULTS AND DISCUSSIONS
Working on four posteriors of D∆t and D d published by H0LiCOW program and combining reconstructed unanchored SN dataset, we obtain the final distributions for Hubble constant H0 and PPN parameter γPPN.Our modelindependent constraints are H0 = 72.9+2.0 −2.3 km s −1 Mpc −1 and γPPN = 0.89 +0.17 −0.15 (median value plus the 16 th and 84 th percentiles around this) for combination of all lenses.The one-dimensional posterior distributions are shown in Fig. 2. The numerical constraint results for four individual lenses can be found in Table 1.One can see that constraint results of all the lenses show that GR is supported within observational uncertainties, and there is no obvious evidence of GR deviation.This can be compared to the results of Yang et al. (2020) of H0 = 73.65 1.95  −2.26 km s −1 Mpc −1 and γPPN = 0.87 +0.19  −0.17 within a flat ΛCDM for combining all lenses, our results are agreement with their results, which support the robustness of our method.Our results on H0 are consistent with the time delay and time delay plus SN results (Collett et al. 2019;Wong et al. 2020;Taubenberger et al. 2019) made within specific cosmological models or with polynomial fitting of the distance relation.However, it should be stressed that, comparing with assuming a specific model, combinations of strong lensing and current astronomical probe such as SN Ia can reduce possible bias in our work.More importantly, there is no significant increase in uncertainty.
On the other hand, Hubble constant and PPN parameter are in fact completely degenerate together, this can also be seen in work of Yang et al. (2020).From a theoretical point of view, PPN parameter is encoded in the inferred angular diameter distance to the lens,      2. For combination of all lenses, the γPPN parameter is even closer to one.However, for the lens J1206, the corresponding constraint result is γPPN = 1.78 +0.87  −0.67 .Although GR is somewhat deviated within the 1σ confidence level, GR is still supported within the ∼ 1.2σ confidence level.As pointed out in work of Millon et al. (2020), the dispersion measurements do not play a significant role in the H0 estimation of the H0LiCOW analysis, except for the lens J1206.In other word, in the analysis of the distance ratios, since we have not given any value for H0, the degeneracy between H0 and γPPN causes a shift in γPPN.In addition, for the lens B1608, the constraint of γPPN using the distance ratio does not change at all, the D∆t and D d of this lens are completely independent, because of the lack of blind analysis of this lens.Since no new data has been added, we do not expect any improvement in precision, and our results show this point.

CONCLUSION
In this work, we propose a model-independent method for simultaneous constraining on Hubble constant and post-Newtonian parameter using time-delay strong lensing sys-tems and observational SNe Ia.To match the lensing and source redshifts of the four strong lensing systems analyzed by H0LiCOW program, we use GP regression to reconstruct distance from observations of SNe Ia instead of assuming a specific model.Although the observations of SNe Ia provide the unanchored or relative distance, strong lensing systems encode absolute distance.Thus, such dataset combinations can anchor cosmological distances.
Firstly, we directly use four posteriors of D∆t (inferred from lensing mass) and D d (inferred from dynamic mass) published by H0LiCOW program, combining with reconstructed unanchored SN dataset.For combination of all lenses, we find that the constraint result on PPN parameter is γPPN = 0.89 +0.17 −0.15 which demonstrates that GR is supported within observational uncertainties.The result on Hubble constant is H0 = 72.9+2.0 −2.3 km s −1 Mpc −1 , which is consistent with the time delay and time delay plus SN results (Collett et al. 2019;Wong et al. 2020;Taubenberger et al. 2019) made within specific cosmological models or with polynomial fitting of the distance relation.
Secondly, in order to avoid the influence of H0 on GR testing, we do not assume any value for H0, and consider a ratio of distance D∆t/D ′ d method to test PPN parameter.This method can independently give a constraint on γPPN, because the distance ratio contains only velocity dispersion measurements and the dispersion measurements do not play a significant role in the H0.The mean values of γPPN for four individual lenses have a little changes (though not significant) using ratio of distance method.However, the constraint result using the lens J1206 is γPPN = 1.78 +0.87 −0.67 , which no longer supports that GR holds within observational uncertainty.
As a final remark, we point out that time-delay lenses plus observations of SN provide a quite promising and model independent method to test General Relativity.The uncertainty of model-independent analyses with such dataset combination can be comparable to the uncertainty of assuming specific models, while reducing possible biases.We also look forward to a large amount of future data, not only from strong lensing system, but also from the SN Ia, allowing us to further improve the measurements of H0 and γPPN.In the future, current surveys such as Dark Energy Survey (DES) (Treu et al. 2018) and Hyper SuprimeCam Survey (HSC) (More et al. 2017), and the upcoming wide-area and deep surveys like the Large Synoptic Survey Telescope (LSST) (Oguri & Marshall 2010) and Euclid and WFIRST satellites (Barnacka 2018;Petrushevska et al. 2018) with much wider field of view and higher sensitivity will be able to discover and precisely localize a large number of lensed quasars and even other lensed sources, which will have well-measured time delays.The High-resolution imaging from space telescopes such as HST or ground-based adaptive optics will help to better model the stellar kinematic in lensing galaxies.In addition, SN Ia data will continue to improve as well, playing an important role as a dense sampler of cosmic expansion history over a wide range of redshifts.All these help us to carry out more accurate analysis of General Relativity, Hubble constant, and cosmology in our subsequent work.

ACKNOWLEDGMENTS
This work was supported by National Natural Science Foundation of China under Grant Nos.12203009, 12222302, 11973034.Liu.T.-H was supported by Chutian Scholars Program in Hubei Province (X2023007).Liao.K was supported by Funds for the Central Universities (Wuhan University 1302/600460081).

DATA AVAILABILITY STATEMENTS
The data underlying this article will be shared on reasonable request to the corresponding author.

Figure 1 .
Figure 1.The unanchored luminosity distance curves H 0 D L (z) reconstructed from the Pantheon SN Ia data for a representative sample of the 1000 GP realizations.

Figure 2 .
Figure 2. The simultaneous constraints of H 0 (left panel) and γ PPN (right panel) from four of the H0LiCOW lenses.The dashed line is γ PPN = 1 predicted by GR.
PPN parameter and cosmological distance (related to Hubble constant) are degenerated.For instance, if γPPN increase, it means the enhancement of gravitational force.One also can keep gravity unmodified, but change the corresponding distances.To try to break this degeneracy and do not assume any value for H0, we can also consider ratios of distances D∆t/D ′ d = 2 1+γ PPN Ds/D ds , which are independent of H0.The final results are displayed in Fig 3. We see that the mean values of γPPN for four individual lenses have a little changes (though not significant), the numerical results are shown in Table

Figure 3 .
Figure 3.The one-dimensional posterior distributions of γ PPN using distance ratio method with four H0LiCOW lenses.The dashed line is γ PPN = 1 predicted by GR.

Table 1 .
Summary of the constraints on the Hubble constant H 0 and γ PPN from four of the H0LiCOW lenses.

Table 2 .
Summary of the constraints on γ PPN using the distance ratio method from four H0LiCOW lenses.