Direct Estimate of the Post-Newtonian Parameter and Cosmic Curvature from Galaxy-scale Strong Gravitational Lensing

Einstein's theory of general relativity (GR) has been precisely tested on solar system scales, but extragalactic tests are still poorly performed. In this work, we use a newly compiled sample of galaxy-scale strong gravitational lenses to test the validity of GR on kiloparsec scales. In order to solve the circularity problem caused by the preassumption of a specific cosmological model based on GR, we employ the distance sum rule in the Friedmann-Lema\^{\i}tre-Robertson-Walker metric to directly estimate the parameterized post-Newtonian (PPN) parameter $\gamma_{\rm PPN}$ and the cosmic curvature $\Omega_k$ by combining observations of strong lensing and Type Ia supernovae. This is the first simultaneous measurement of $\gamma_{\rm PPN}$ and $\Omega_k$ without any assumptions about the contents of the universe or the theory of gravity. Our results show that $\gamma_{\rm PPN}=1.11^{+0.11}_{-0.09}$ and $\Omega_{k}=0.48^{+1.09}_{-0.71}$, indicating a strong degeneracy between the two quantities. The measured $\gamma_{\rm PPN}$, which is consistent with the prediction of 1 from GR, provides a precise extragalactic test of GR with a fractional accuracy better than 9.0\%. If a prior of the spatial flatness (i.e., $\Omega_{k}=0$) is adopted, the PPN parameter constraint can be further improved to $\gamma_{\rm PPN}=1.07^{+0.07}_{-0.07}$, representing a precision of 6.5\%. On the other hand, in the framework of GR (i.e., $\gamma_{\rm PPN}=1$), our results are still marginally compatible with zero curvature ($\Omega_k=-0.12^{+0.48}_{-0.36}$), supporting no significant deviation from a flat universe.


INTRODUCTION
Einstein's theory of general relativity (GR) is one of the major pillars of modern physics. Any possible violation of GR would have far-reaching consequences for our understanding of fundamental physics; testing GR at a much higher precision has therefore been one of the most enduring pursuits of scientists. At the post-Newtonian level, the validity of GR can be tested by constraining the parameterized post-Newtonian (PPN) parameter γ PPN , since GR predicts exactly γ PPN ≡ 1 (Thorne & Will 1971;Will , 2014. Here, γ PPN stands for the amount of space-curvature generated by a unit rest mass. On solar system scales, tests of GR through numerical values of γ PPN have reached high precision. By measuring the arrival-time delay of radar signals passing close to the Sun, the Cassini spacecraft yielded an agreement with GR to 10 −3 %, i.e., γ PPN = 1 + (2.1 ± 2.3) × 10 −5 (Bertotti et al. 2003). However, current extragalactic tests of GR are much less precise. On scales of 10-100 Mpc, only ∼ 20% precision on the constraints of γ PPN has been obtained using the joint measurements of weak gravitational lensing and redshift-space distortions (Song et al. 2011;Simpson et al. 2013;Blake et al. 2016). On megaparsec scales, γ PPN has been limited to just 30% precision by analyzing the mass profiles of galaxy clusters (Wilcox et al. 2015;Pizzuti et al. 2016).
On kiloparsec scales, strong gravitational lensing (SGL) systems, combined with stellar dynamical data of lensing galaxies, provide an effective tool to verify the weak-field metric of gravity. For a specific SGL system with the foreground galaxy acting as a lens, multiple images, arcs, or even an Einstein ring can form with angular separations close to the so-called Einstein radius (Chakraborty & SenGupta 2017). In theory, the Einstein radius is related to the mass of the lens, the PPN parameter γ PPN , and a ratio of three angular diameter distances (i.e., the distances from the observer to the lens and the source, D l and D s , and the distance between the lens and the source D ls ) (Cao et al. 2015). With the required angular diameter distances and measurements of the lens mass and the Einstein radius, one can therefore constrain γ PPN and test whether GR is a suitable theory of gravity on the corresponding scales. This method was first performed on 15 lensing galaxies from the Sloan Lens ACS Survey by Bolton et al. (2006), which yielded γ PPN = 0.98 ± 0.07 based on prior assumptions on galaxy structure from local observations. Subsequently, different SGL samples have been used to test the accuracy of GR (Smith 2009;Schwab et al. 2010;Cao et al. 2017;Collett et al. 2018;Yang et al. 2020;Liu et al. 2021). In most previous studies, the distance information required to constrain the PPN parameter γ PPN is provided by the prediction of the standard ΛCDM cosmological model. It should, however, be emphasized that ΛCDM is established based on the framework of GR. Thus, there is a circularity problem in testing GR . To overcome this problem, one has to determine the lensing distance ratio in a cosmology-independent way.
The circularity problem can be alleviated by determining the two distances D l and D s through observations of Type Ia supernovae (SNe Ia). But, the distance D ls cannot be determined directly from the observations. In the Friedmann-Lemaître-Robertson-Walker (FLRW) metric, these three distances are related via the distance sum rule (DSR), which depends on the curvature parameter of the universe Ω k . Turning this around, supposing that the universe is described by the FLRW metric, we can use combined observations of strong lensing and SNe Ia to estimate not only γ PPN but also Ω k independently of the cosmological model (Cao et al. 2017). Based on the DSR in the FLRW metric, and assuming that GR is valid (i.e., γ PPN = 1), model-independent constraints on the cosmic curvature Ω k have been implemented by combining SGL systems with other distance indicators (Räsänen et al. 2015;Liao et al. 2017;Xia et al. 2017;Denissenya et al. 2018;Li et al. 2018aLi et al. ,b, 2019Cao et al. 2019Cao et al. , 2021Collett et al. 2019;Liao 2019;Qi et al. 2019a;Liu et al. 2020;Qi et al. 2019bQi et al. , 2021Wang et al. 2020;Wei & Melia 2020;Zhou & Li 2020;Dhawan et al. 2021). Without the prior assumption on GR, Cao et al. (2017) proposed that this cosmology-independent method could be extended to study the degeneracy between the PPN parameter γ PPN and the curvature parameter Ω k . They used the simulated strong-lensing data to estimate both γ PPN and Ω k . We will now for the first time apply such a method to real data.
We should note that a recent work by Liu et al. (2021) used strong lensing and SNe Ia to obtain model-independent constraints on γ PPN within the framework of the flat FLRW metric (i.e., Ω k = 0). However, Cao et al. (2017) proved that there exists a significant degeneracy between γ PPN and Ω k by simulation. Obviously, a simple flatness assumption may lead to a biased estimate of γ PPN , even if the real curvature is tiny. Therefore, it would be better to simultaneously optimize γ PPN and Ω k , as we do in this work.
The outline of this work is as follows. In Section 2, we introduce the gravitational lensing theory and the DSR method. In Section 3, we describe the observational data used for our analysis. Model-independent constraints on γ PPN and Ω k are presented in Section 4. Finally, a brief summary and discussions are given in Section 5.

METHODOLOGY
In the limit of a weak gravitational field, the general form of the Schwarzschild metric for a point mass M can be written as (1) where γ PPN is the PPN parameter and Ω is the angle in the invariant orbital plane. In GR, γ PPN is predicted to be 1.
can be inferred from the spectroscopic measurement of the lens velocity dispersion. Here we adopt a general mass model with power-law density profiles for the lensing galaxy (Koopmans 2006;Cao et al. 2016): where r is the spherical radial coordinate from the lens center, ρ(r) is the total (i.e., luminous plus dark matter) mass density, and ν(r) denotes the luminosity density of stars. The parameter β(r) represents the anisotropy of the stellar velocity dispersion, which relates to the velocity dispersions, σ 2 t and σ 2 r , in the tangential and radial directions. Also, α and δ are the slopes of the power-law density profiles. It is worth noting that the total mass density slope α is significantly dependent on both the lens redshift z l and the surface mass density (e.g., Sonnenfeld et al. 2013;Chen et al. 2019). Chen et al. (2019) proved that the most compatible lens mass model is where α 0 , α z , and α s are free parameters. HereΣ denotes the normalized surface mass density of the lensing galaxy, which is given byΣ = (σ0/100 km s −1 ) 2 Reff/10 h −1 kpc , where σ 0 is the observed velocity dispersion, h = H 0 /(100 km s −1 Mpc −1 ) is the reduced Hubble constant, and R eff is the half-light radius of the lensing galaxy. In the literature, the velocity anisotropy parameter β is usually assumed to be independent of r (e.g., Koopmans et al. 2006;Treu et al. 2010). From a wellstudied sample of nearby elliptical galaxies (Gerhard et al. 2001), the posterior probability of β is found to be characterized by a Gaussian distribution, β = 0.18 ± 0.13, that is extensively adopted in previous works (e.g., Bolton et al. 2006;Schwab et al. 2010;Cao et al. 2017;Chen et al. 2019;Liu et al. 2021). Following these previous works, we will marginalize the anisotropy parameter β using a Gaussian prior of β = 0.18 ± 0.13 over the range of β − 2σ β ,β + 2σ β , whereβ = 0.18 and σ β = 0.13.
Based on the radial Jeans equation in spherical coordinate, the radial velocity dispersion of luminous matter in earlytype lens galaxies can be expressed as where M(r) is the total mass contained within a spherical radius r. With the mass density profiles in Equation (5), we can derive the relation between the dynamical mass M dyn E enclosed within the Einstein ring radius R E and M(r) as (see Koopmans 2006;Chen et al. 2019 for the detailed derivation) where λ(x) = Γ x−1 2 /Γ x 2 stands for the ratio of two respective Gamma functions. By substituting Equations (8) and (5) into Equation (7), one can have where ξ = α + δ − 2.
The actual velocity dispersion of the lensing galaxy is effectively averaged by line-of-sight luminosity and measured over the effective spectroscopic aperture R A , which can be expressed as (see Chen et al. 2019 for the detailed derivation) where Lastly, with the relations expressed in Equations (2) and (4), Equation (10) can be rewritten as From the spectroscopic data, one can measure the lens velocity dispersion σ ap inside the circular aperture with the angular radius θ ap . In practice, the luminosity-weighted average of the line-of-sight velocity dispersion σ ap measured within a certain aperture should be normalized to a typical physical aperture with the radius θ eff /2, where θ eff = R eff /D l is the effective angular radius of the lensing galaxy. Following Chen et al. (2019), we adopt the value of the correction factor η = −0.066 ± 0.035 from Cappellari et al. (2006). Then, we can calculate the total uncertainty of σ obs 0 using the expression where ∆σ stat 0 is the statistical uncertainty propagated from the measurement error of σ ap . The uncertainty caused by the aperture correction, ∆σ AC 0 , is propagated from the error of η. The extra mass contribution from other matters (outside of the lensing galaxy) along the line of sight in the estimation of M grl E can be treated as a systematic uncertainty ∆σ sys 0 , which contributes an uncertainty of ∼ 3% to the velocity dispersion (Jiang & Kochanek 2007).
With Equation (12), the theoretical value of the velocity dispersion within the radius θ eff /2 takes the form (Koopmans 2006) For the case of α = δ = 2 and β = 0, the mass model is reduced to the singular isothermal sphere (SIS) model, and the theoretical value of the velocity dispersion is simplified as σ SIS = c 2 4π 2 (1+γPPN) Ds D ls θ E . By comparing the observational values of the velocity dispersions (Equation (13)) with the corresponding theoretical ones (Equation (15)), one can place constraints on the PPN parameter γ PPN . For this purpose, it is also necessary to know the distance ratio D s /D ls , which is conventionally calculated in the context of flat ΛCDM (Schwab et al. 2010;Cao et al. 2017). However, a circularity problem exists in this approach because the standard ΛCDM cosmological model is built on the framework of GR . In order to avoid the circularity problem, we will apply a cosmology-independent method to constrain γ PPN . This method is based on the sum rule of distances along null geodesics of the FLRW metric.

Distance Sum Rule
If space is exactly homogeneous and isotropic, the FLRW metric can be used to describe the spacetime geometry of the universe. In the FLRW metric, the dimensionless comoving distance d(z l , z s ) ≡ (H 0 /c)(1 + z s )D A (z l , z s ) is given by where Ω k is the curvature parameter and E(z) = H(z)/H 0 is the dimensionless Hubble parameter. Also, sinn(x) = sinh (x) for Ω k > 0 and sinn(x) = sin(x) for Ω k < 0. For a flat universe with Ω k = 0, Equation (16) reduces to a linear function of the integral. For an SGL system with the notations d(z) ≡ d(0, z), d l ≡ d(0, z l ), d s ≡ d(0, z s ), and d ls ≡ d(z l , z s ), a simple sum rule of distances in the FLRW framework can be easily derived as (Peebles 1993;Bernstein 2006;Räsänen et al. 2015): This relation is very general because it only assumes that geometrical optics holds and that light propagation is described with the FLRW metric. Once the derived Ω k from the three distances (d l , d s , and d ls ) is observationally found to be different for any two pairs of (z l , z s ), we can rule out the FLRW metric. Given independent measurements of d l and d s on the right side of Equation (17), we are able to access the dimensionless distance ratio d ls /d s , 1 depending only on the curvature parameter Ω k (Geng et al. 2020;Liu et al. 2020;Zheng et al. 2021). Therefore, we can directly determine γ PPN and Ω k from Equations (15) and (17)  In order to obtain model-independent estimate of γ PPN and Ω k via Equations (15) and (17), we need to know the distances d l and d s on the right-hand-side terms of Equation (17). In principle, we can use different kinds of distance indicators such as standard candles, sirens, and rulers for providing these two distances. Here, we use SN Ia observations to obtain d l and d s . Scolnic et al. (2018) released the largest combined sample of SNe Ia called Pantheon, which contains 1,048 SNe in the redshift range 0.01 < z < 2.3. Generally, the observed distance modulus of each SN is given by µ SN = m B + κ · X 1 − ω · C − M B , where m B is the observed peak magnitude in the rest-frame B band, X 1 and C are the light-curve stretch factor and the SN color at maximum brightness, respectively, and M B is a nuisance parameter that represents the absolute B-band magnitude of a fiducial SN. Here, κ and ω are two light-curve parameters, which could be calibrated to zero through a method called BEAMS with Bias Corrections (BBC; Kessler & Scolnic 2017). With the BBC method, Scolnic et al. (2018) reported the corrected apparent magnitudes m corr = µ SN + M B for all SNe. Therefore, the observed distance moduli µ SN can be directly obtained by subtracting M B from m corr .
As proposed in Räsänen et al. (2015), we determine the dimensionless distances d l and d s by fitting a polynomial to the Pantheon SN Ia data. Here, we parameterize the dimensionless distance function as a third-order polynomial with initial conditions d(0) = 0 and d ′ (0) = 1, i.e., d(z) = z + a 1 z 2 + a 2 z 3 , where a 1 and a 2 are two free parameters that need to be optimized along with the absolute magnitude M B . We find that higher-order polynomials do not improve the fitting performance, taking into account the larger number of free parameters. That is, a simple third-order polynomial is flexible enough to fit the SN Ia data. Given a vector of distance residuals of the Pantheon SN sample that may be expressed as ∆μ =μ SN −μ model , wherê µ SN (μ model ) is the observed (model) vector of distance moduli, the likelihood for the model fit is defined by where Cov is a covariance matrix that includes both statistical and systematic uncertainties of SNe. Here the observed vectorμ SN is given by µ SN,i = m corr,i − M B , and the model vec- According to the analysis in Section 2.1, one can learn that the underlying method requires the following observational information of each SGL system, including the source redshift z s , the lens redshift z l , the Einstein angle θ E , the halflight angular radius of the lensing galaxy θ eff , the spectroscopic aperture angular radius θ ap , and the lens velocity dispersion σ ap measured within θ ap .
Recently, Chen et al. (2019) compiled a sample of 161 galaxy-scale SGL systems with gravitational lensing and stellar velocity dispersion measurements. In this sample, the slopes of the luminosity density profile δ of 130 SGL systems were measured by fitting the two-dimensional power-law luminosity profile convolved with the instrumental point spread function to imaging data over a circle of radius θ eff /2 centered on the lens galaxies. By constraining the cosmological parameter Ω m separately with the entire sample of 161 SGL systems (treating δ as a universal parameter for all lenses) and the truncated sample of 130 systems (treating δ as an observable for each lens), Chen et al. (2019) suggested that the intrinsic scatter δ among the lenses should be considered in order to get an unbiased estimate of Ω m . Therefore, we adopt this truncated sample of 130 SGL systems with δ measurements for the analysis demonstrated in this paper. The redshift ranges of lens and source galaxies of these 130 SGL systems are 0.0624 ≤ z l ≤ 0.7224 and 0.1970 ≤ z s ≤ 2.8324, respectively.
One of the limitations we must deal with in using the SGL data, however, is that the SN Ia measurements extend only to z = 2.3. As such, only a subset of the SGL sample that overlaps with the SN Ia catalog is actually available. Our analysis will therefore be based only on the 120 SGL systems with z s < 2.3. The likelihood function for strong-lensing data is then constructed as The third-order polynomial modeling the distance function d(z) has two free parameters (a 1 and a 2 ). The absolute magnitude M B enters into the SN likelihood as a nuisance parameter. The PPN parameter γ PPN and the lens model parameters (α 0 , α z , and α s ) enter into the SGL likelihood as four free parameters. In addition, the d ls /d s given by Equation (17) Figure 1. 1D and 2D marginalized probability distributions with 1σ and 2σ confidence contours for the PPN parameter γPPN and cosmic curvature Ω k . The dashed lines correspond to a flat universe with the validity of GR (Ω k = 0, γPPN = 1).
involves the curvature parameter Ω k , making it eight free parameters in total.
By marginalizing the lens model parameters (α 0 , α z , and α s ), the polynomial coefficients (a 1 and a 2 ), and the SN absolute magnitude M B , we obtain the 1D and 2D marginalized probability distributions with 1σ − 2σ confidence regions for γ PPN and Ω k , which are presented in Figure 1. These contours show that, whereas Ω k = 0.48 +1.09 −0.71 is weakly constrained, we can set a good limit of γ PPN = 1.11 +0.11 −0.09 at the 68% confidence level. The inferred value of the PPN parameter is compatible with the prediction of γ PPN = 1 from GR. The constraint accuracy of γ PPN is about 9.0%. As shown in Table 1, the lens model parameters are constrained to be α 0 = 1.266 +0.105 −0.105 , α z = −0.332 +0.169 −0.188 , and α s = 0.656 +0.065 −0.065 at the 68% confidence level, which are consistent with the results of Chen et al. (2019). We find that α z = 0 is ruled out at ∼ 2σ level and α s = 0 is ruled out at ∼ 10σ level, confirming the significant dependencies of the total mass density slope α on both the lens redshift and the surface mass density.
If a prior of flatness (i.e., Ω k = 0) is adopted, the resulting posterior probability distribution for γ PPN is shown in Figure 2. The result γ PPN = 1.07 +0.07 −0.07 (1σ confidence level) is in good agreement with γ PPN = 1 predicted by GR, and its constraint accuracy is improved to about 6.5%. If we instead assume GR holds (i.e., γ PPN = 1), and allow Ω k to be a free parameter, we obtain the marginalized probability distribution for Ω k , as illustrated in Figure 3. The curvature parameter is constrained to be Ω k = −0.12 +0.48 −0.36 , consistent with a flat uni-  verse. The corresponding results for all parameters are summarized in lines 1-3 of Table 1 for the cases with no priors, the prior of Ω k = 0, and the prior of γ PPN = 1, respectively. The comparison among these three cases indicates that the nuisance parameters (α 0 , α z , α s , a 1 , a 2 , and M B ) have little effect on the PPN parameter γ PPN and cosmic curvature Ω k .

SUMMARY AND DISCUSSIONS
Galaxy-scale strong-lensing systems with measured stellar velocity dispersions provide an excellent extragalactic test of GR by constraining the PPN parameter (γ PPN ). Measuring γ PPN in this manner, however, one has to know the lensing distances (the angular diameter distances between the source, lens, and observer), which are conventionally calculated within the standard ΛCDM cosmological model. Because ΛCDM itself is built on the theoretical framework of GR, these distance calculations would involve a circularity problem. In this work, aiming to overcome the circularity problem, we have applied the DSR in the FLRW metric to obtain cosmology-independent constraints on both γ PPN and the cosmic curvature parameter Ω k . Though the DSR method has been used to directly infer the value of Ω k by confronting observations of SGL systems with SN Ia luminosity distances, the simultaneous measurement of Ω k and γ PPN has not yet been achieved by the community in the literature.
Combining 120 well-measured SGL systems at z s < 2.3 with the latest Pantheon SN Ia observations, we have simultaneously placed limits on γ PPN and Ω k without any assumptions about the contents of the universe or the theory of gravity. This analysis suggests that the PPN parameter is constrained to be γ PPN = 1.11 +0.11 −0.09 , representing a precision of 9.0%, consistent with the prediction of 1 from GR at a 68% confidence level. Meanwhile, the optimized curvature parameter is Ω k = 0.48 +1.09 −0.71 . If using the spatial flatness as a prior, we find γ PPN = 1.07 +0.07 −0.07 , representing an agreement with GR to 6.5%. Assuming GR is valid and allowing Ω k to be a free parameter, we infer that Ω k = −0.12 +0.48 −0.36 . This cosmic curvature value does not significantly deviate from a flat universe.
Previously, Cao et al. (2017) obtained a 25% precision on the determination of γ PPN by analyzing a sample of 80 lenses in the flat ΛCDM model. Under the assumption of the fiducial ΛCDM cosmology with parameters taken from Planck observations, Collett et al. (2018) estimated γ PPN on scales around 2 kpc to be 0.97 ± 0.09 (representing a 9.3% precision measurement) by using a nearby SGL system, ESO 325-G004. Yang et al. (2020) derived γ PPN = 0.87 +0.19 −0.17 (representing a precision of 21%) for flat ΛCDM using a sample of four time-delay lenses. Within the framework of the flat FLRW metric, Liu et al. (2021) used 120 strong-lensing data to obtain a model-independent constraint of γ PPN = 1.065 +0.064 −0.074 (representing a precision of 6.5%) by implementing Gaussian processes to extract the SN distances. Despite not assuming a specific cosmological model, the uncertainties in our constraints are comparable to these previous results. Most importantly, our method offers a new cosmology-independent way of simultaneously constraining both γ PPN and Ω k .
Forthcoming lens surveys such as the Large Synoptic Survey Telescope, with improved depth, area, and resolution, will be able to increase the current galactic-scale lens sample sizes by orders of magnitude (Collett 2015). With such abundant observational information in the future, the mass-dynamical structure of the lensing galaxies will be better characterized, and model-independent constraints on the PPN parameter γ PPN and cosmic curvature Ω k , as discussed in this work, will be considerably improved.
Finally, we investigated whether the approximation of the dimensionless distance function d(z) (as a linear polynomial; see Equation 18) affects the inference of γ PPN . To probe the  Table 1), it is clear that the linear polynomial function provides a good approximation of d(z) and the adoption of the exact expression for d(z) in the flat ΛCDM model only has a minimal influence on these results.