# Galaxy-Scale Test of General Relativity with Strong Gravitational Lensing

Xiao-Hui Liu  Department of Physics, College of Sciences, Northeastern
University, Shenyang 110819, China  Zhen-Hua Li  Department of Physics,
College of Sciences, Northeastern University, Shenyang 110819, China  Jing-
Zhao Qi  Department of Physics, College of Sciences, Northeastern University,
Shenyang 110819, China  Xin Zhang  Department of Physics, College of Sciences,
Northeastern University, Shenyang 110819, China  Frontiers Science Center for
Industrial Intelligence and Systems Optimization, Northeastern University,
Shenyang 110819, China  Key Laboratory of Data Analytics and Optimization for
Smart Industry (Northeastern University), Ministry of Education, Shenyang
110819, China

###### Abstract

Although general relativity (GR) has been precisely tested at the solar system
scale, precise tests at a galactic or cosmological scale are still relatively
insufficient. Here, in order to test GR at the galactic scale, we use the
newly compiled galaxy-scale strong gravitational lensing (SGL) sample to
constrain the parameter $\gamma_{\rm{PPN}}$ in the parametrized post-Newtonian
(PPN) formalism. We employ the Pantheon sample of type Ia supernovae
observation to calibrate the distances in the SGL systems using the Gaussian
Process method, which avoids the logical problem caused by assuming a
cosmological model within GR to determine the distances in the SGL sample.
Furthermore, we consider three typical lens models in this work to investigate
the influences of the lens mass distributions on the fitting results. We find
that the choice of the lens models has a significant impact on the constraints
on the PPN parameter $\gamma_{\rm{PPN}}$. By using the minimum $\chi^{2}$
comparison and the Bayesian information criterion as evaluation tools to make
comparisons for the fitting results of the three lens models, we find that the
most reliable lens model gives the result of
$\gamma_{\rm{PPN}}=1.065^{+0.064}_{-0.074}$, which is in good agreement with
the prediction of $\gamma_{\rm{PPN}}=1$ by GR. As far as we know, our 6.4%
constraint result is the best result so far among the recent works using the
SGL method.

Strong gravitational lensing (1643), General relativity (641), Cosmology (343)

##  1 Introduction

It is well known that Einstein’s theory of general relativity (GR), one of the
cornerstones of modern physics, is an extremely successful description of
gravity, which has been confirmed in all current experiments, especially in
some famous experiments, e.g. Dyson et al. (1920); Pound & Rebka (1960);
Shapiro (1964); Taylor et al. (1979). In recent years, the continued pursuit
for tests of GR at a much higher precision has yielded many accomplishments,
especially at the solar system scale. However, the long-range nature of
gravity on the extragalactic scale is still poorly understood. Any possible
violation of GR would have the potential to advance our understanding of many
important problems in fundamental physics, for which a good example is the
accelerating expansion of the Universe.

In the framework of GR, an exotic component with negative pressure, called
dark energy, was introduced to explain the cosmic accelerated expansion.
However, up to now, not only the nature of dark energy is still poorly
understood, but also any current phenomenological models of it are more or
less burdened with some theoretical problems (Weinberg, 1989; Carroll, 2001;
Copeland et al., 2006). Another perspective for understanding the cosmic
acceleration is modifying Einstein’s theory of gravity. Indeed, the true
gravity theory may deviate GR at cosmological scale (Sotiriou & Faraoni,
2010). Some geometric effects from modifying gravity could also successfully
explain the cosmic accelerated expansion without introducing exotic dark
energy (Sotiriou & Faraoni, 2010; Qi et al., 2017, 2016; Pan et al., 2021;
Bengochea & Ferraro, 2009; Zhang et al., 2014; Li et al., 2015; Geng et al.,
2015; Zhao et al., 2018). Therefore, clarifying whether it is dark energy or
modified gravity that is responsible for the accelerating expansion of the
Universe requires an accurate test of GR on cosmological scales. At
cosmological scales, GR with dark energy could be distinguished from modified
gravity by accurately measuring the growth of structure (Zhang et al., 2007;
Reyes et al., 2010; Song et al., 2007). Moreover, the test of GR at the
galactic scale could also be well established, which is the focus of this
paper, and its details will be described in the following.

The validity of GR can be verified by constraining the parametrized post-
Newtonian (PPN) parameter $\gamma_{\rm{PPN}}$, which describes how much space
curvature is provided by the unit rest mass of the objects along or near the
path of particles (Thorne & Will, 1971; Bertotti et al., 2003). The PPN
parameter $\gamma_{\rm{PPN}}=1$ for GR. Thus, any deviation from
$\gamma_{\rm{PPN}}=1$ on large scales indicates a possible violation of GR.

Strong gravitational lensing (SGL) could provide an effective way to constrain
the PPN parameter on the galactic scale. This method was first adopted in
(Bolton et al., 2006b). They used 15 lenses sample from Sloan Lens ACS (SLACS)
Survey (Bolton et al., 2006a), giving a result of $\gamma_{\rm{PPN}}=0.98\pm
0.07$ (Bolton et al., 2006b), based on priors on galaxy structure taken from
local observations. Subsequently, many works have provided effective tests for
GR with different and larger SGL samples (Smith, 2009; Schwab et al., 2010;
Cao et al., 2017; Collett et al., 2018; Yang et al., 2019). More recently, a
larger SGL sample including 161 galaxy-scale strong lensing systems has been
released (Chen et al., 2019), which makes it possible to further explore the
validity of GR on galactic scale at higher precision.

In this paper, we will not only improve the constraints on the PPN parameter
with the larger SGL sample in exploring the validity of GR, but also refine
the test of GR from two aspects compared to the previous works motioned above
(Bolton et al., 2006b; Smith, 2009; Schwab et al., 2010; Cao et al., 2017).
Firstly, it should be emphasized that the influence of lens model selection
cannot be ignored in this method for testing GR, since the theoretical value
of light deflection is determined by the mass distribution model of lens
galaxies and gravitational theory together. In the previous works, most of
them only considered a single lens mass model. We will fully consider the
different lens mass models and investigate their impacts on the results.
Secondly, in previous related works, the distance information required to
constrain the PPN parameter is conventionally provided by the prediction of
the standard $\Lambda$ cold dark matter ($\Lambda$CDM) model that is derived
within the framework of GR, which more or less leads to a logic problem in
testing gravity. To avoid this problem, in this work, by using the method of
Gaussian Process, we will construct a smooth distance-redshift curve directly
from the type Ia supernovae (SN Ia) observation to calibrate the distances in
the SGL sample.

This paper is organized as follows. In Section 2, we describe the methodology
of testing GR by SGL and Gaussian Process. In addition, the SGL sample and the
latest Pantheon SN Ia sample are also introduced. In Section 3, the constraint
results in different lens mass models, and the model comparison based on the
information criterion are presented. Finally, the conclusions and discussions
are given in Section 4.

##  2 Methodology and Observational Data

In the limit of a weak gravitational field, the PPN form of Schwarzschild
metric can be expressed as

$\displaystyle ds^{2}=$ $\displaystyle
c^{2}dt^{2}\left(1-2GM/c^{2}r\right)-dr^{2}\left(1+2\gamma_{\rm{PPN}}GM/c^{2}r\right)$
(1) $\displaystyle-r^{2}d\Omega^{2},$

where $\gamma_{\rm{PPN}}$ is the PPN parameter, and $\gamma_{\rm{PPN}}=1$
stands for GR. $M$ is the mass of the central object, in this case the mass of
the lens galaxy.

###  2.1 Strong Gravitational Lensing Systems

For a SGL system, the Einstein radius $\theta_{\mathrm{E}}$ reflecting the
angular separations between multiple images is corresponding to the mass
contained within a cylinder of Einstein radius $M_{\mathrm{E}}$,

$\theta_{\mathrm{E}}=\sqrt{\frac{1+\gamma_{\rm{PPN}}}{2}}\left(\frac{4GM_{\mathrm{E}}}{c^{2}}\frac{D_{\rm{ds}}}{D_{\rm{s}}D_{\rm{d}}}\right)^{1/2},$
(2)

where $D_{\rm{s}}$ is the angular diameter distance between observer and
source, $D_{\rm{d}}$ is the angular diameter distance between observer and
lens, and $D_{\rm{ds}}$ is the angular diameter distance between lens and
source. By introducing the radius perpendicular to the line of sight
$R_{\rm{E}}=D_{\rm{d}}\theta_{\rm{E}}$ and rearranging Eq. (2), we can obtain

$\frac{GM_{\mathrm{E}}}{R_{\mathrm{E}}}=\frac{2}{(1+\gamma_{\rm{PPN}})}\frac{c^{2}}{4}\frac{D_{\rm{s}}}{D_{\rm{ds}}}\theta_{\mathrm{E}}.$
(3)

It should be noted that the mass $M_{\mathrm{E}}$ is derived from the
gravitational theory and the observations of angular separations between
multiple images, and thus it could also be called the gravitational mass,
denoted as $M_{\mathrm{grl}}^{\mathrm{E}}$. On the other hand, by observing
the velocity dispersion of the lens galaxy and assuming a lens mass model, the
dynamical mass enclosed within the Einstein radius,
$M_{\mathrm{dyn}}^{\mathrm{E}}$, can also be inferred. In theory, these two
masses should be equivalent, i.e.,
$M_{\mathrm{grl}}^{\mathrm{E}}=M_{\mathrm{dyn}}^{\mathrm{E}}$ (Chen et al.,
2019).

To estimate the dynamical mass from the velocity dispersion $\sigma$ of lens
galaxy, here we choose a general mass model based on power-law density
profiles to describe the mass distribution of lens galaxy (Koopmans, 2006):

$\left\\{\begin{array}[]{l}\rho(r)=\rho_{0}\left(r/r_{0}\right)^{-\gamma}\\\
v(r)=v_{0}\left(r/r_{0}\right)^{-\delta}\\\
\beta(r)=1-\sigma_{\theta}^{2}/\sigma_{r}^{2},\end{array}\right.$ (4)

where $\rho(r)$ is the total mass density distribution including luminous
matter and dark matter, and $v(r)$ represents the density distribution of
luminous mass. The parameter $\beta(r)$ characterizes the anisotropic
distribution of velocity dispersion pattern, where $\sigma_{\theta}$ and
$\sigma_{r}$ are the tangential and radial components of the velocity
dispersion, respectively. Here $\gamma$ and $\delta$ are power law indices;
note that $\gamma$ should not be confused with the PPN parameter
$\gamma_{\rm{PPN}}$. If $\gamma=\delta=2$ and $\beta=0$, this mass model will
reduce to the singular isothermal sphere (SIS) model commonly used in studies
of SGL. It is important to note that the total mass density slope $\gamma$ is
probably dependent on redshift (Ruff et al., 2011; Bolton et al., 2012; Cao et
al., 2016; Cui et al., 2017). Therefore, three parametrized models of $\gamma$
commonly used will be considered in this work, namely

$\displaystyle(\rm{i})$ $\displaystyle P_{1}:\gamma=\gamma_{0},$
$\displaystyle(\rm{ii})$ $\displaystyle
P_{2}:\gamma=\gamma_{0}+\gamma_{z}z_{\rm{d}},$ $\displaystyle(\rm{iii})$
$\displaystyle
P_{3}:\gamma=\gamma_{0}+\gamma_{z}z_{\rm{d}}+\gamma_{s}\log\tilde{\Sigma},$

where $\gamma_{0}$, $\gamma_{z}$ and $\gamma_{s}$ are arbitrary constants. For
$P_{1}$ model, the dependences on redshift and mass distribution are not
considered, and this model is also frequently used in previous related works
(Smith, 2009; Schwab et al., 2010; Cao et al., 2017). $P_{2}$ model assumes
that the total mass density distribution of galaxies evolves with redshift.
$z_{d}$ is the redshift of lens galaxy. In $P_{3}$ model, in addition to the
dependence of redshift, the dependence on the surface mass density in the mass
distribution is also considered. Here, $\tilde{\Sigma}$ represents the
normalized surface mass density, expressed as (Chen et al., 2019)

$\tilde{\Sigma}=\frac{\left(\sigma_{0}/100\mathrm{~{}km}\mathrm{~{}s}^{-1}\right)^{2}}{R_{\mathrm{eff}}/10h^{-1}\mathrm{kpc}},$
(6)

where $h=H_{0}/(100~{}\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}})$ and $H_{0}$ is the
Hubble constant. $\sigma_{0}$ is the observed velocity dispersion, and
$R_{\mathrm{eff}}$ is the half-light radius of the lens galaxy.

Based on the well-known spherical Jeans equation and the mass density profiles
in Eq. (4), the mass enclosed within a sphere with radius $r$ is (Schwab et
al., 2010)

$M(r)=\frac{2}{\sqrt{\pi}}\frac{\Gamma(\gamma/2)}{\Gamma\left(\frac{\gamma-1}{2}\right)}\left(\frac{r}{R_{\mathrm{E}}}\right)^{3-\gamma}M_{\mathrm{dyn}}^{\mathrm{E}},$
(7)

where $\Gamma(x)$ is Euler’s Gamma function. The radial velocity distance
dispersion of the luminous mass $\sigma_{r}$ could be expressed as

$\sigma_{r}^{2}(r)=\frac{G\int_{r}^{\infty}\mathrm{d}r^{\prime}r^{\prime
2\beta-2}v\left(r^{\prime}\right)M\left(r^{\prime}\right)}{r^{2\beta}v(r)}.$
(8)

By taking $\xi=\gamma+\delta-2$ and assuming $\beta$ is independent of the
radius $r$, the combination of Eqs. (7) and (8) could infer the expression of
$\sigma_{r}^{2}$,

$\sigma_{r}^{2}(r)=\frac{2}{\sqrt{\pi}}\frac{GM_{\mathrm{dyn}}^{\mathrm{E}}}{R_{\mathrm{E}}}\frac{1}{\xi-2\beta}\frac{\Gamma(\gamma/2)}{\Gamma\left(\frac{\gamma-1}{2}\right)}\left(\frac{r}{R_{\mathrm{E}}}\right)^{2-\gamma}.$
(9)

In the observation of SGL, the actual measurement of velocity dispersion is
the component of luminosity-weighted along the line of sight (LOS) and over
the effective spectrometer aperture. Moreover, for the spectroscopic data, the
velocity dispersion $\sigma_{\rm{ap}}$ measured within certain apertures
$\theta_{\rm{ap}}$ should be normalized to the one within typical physical
aperture via
$\sigma_{0}=\sigma_{\mathrm{ap}}\left[\theta_{\mathrm{eff}}/\left(2\theta_{\mathrm{ap}}\right)\right]^{\eta}$,
where $\theta_{\mathrm{eff}}$ is the effective angular radius of the lens
galaxy, and the value of the correction factor $\eta$ we adopt is
$\eta=-0.066\pm 0.033$ (Jorgensen et al., 1995). The detailed derivation and
description could be found in Koopmans (2006); Chen et al. (2019). Taking all
of these into account, the theoretical value of the velocity dispersion could
be expressed as (Chen et al., 2019)

$\sigma_{0}=\sqrt{\frac{c^{2}}{2\sqrt{\pi}}\frac{2}{(1+\gamma_{\rm{PPN}})}\frac{D_{\mathrm{s}}}{D_{\mathrm{ds}}}\theta_{\mathrm{E}}f(\gamma,\delta,\beta)\left(\frac{\theta_{\mathrm{eff}}}{2\theta_{\mathrm{E}}}\right)^{(2-\gamma)}},$
(10)

where

$\displaystyle f(\gamma,\delta,\beta)=$
$\displaystyle\frac{(2\sqrt{\pi})(3-\delta)}{(\xi-2\beta)(3-\xi)}$ (11)
$\displaystyle\times\left[\frac{\Gamma[(\xi-1)/2]}{\Gamma(\xi/2)}-\beta\frac{\Gamma[(\xi+1)/2]}{\Gamma[(\xi+2)/2]}\right]$
$\displaystyle\times\frac{\Gamma(\gamma/2)\Gamma(\delta/2)}{\Gamma[(\gamma-1)/2]\Gamma[(\delta-1)/2]}.$

Considering the extra mass contribution from matters along the LOS, the
fractional uncertainty of velocity dispersion at the level of 3% will be taken
as the systematic error (Jiang & Kochanek, 2007), and the total uncertainty of
$\sigma_{0}$ from the observation is

$\displaystyle\left(\Delta\sigma_{0}\right)^{2}$
$\displaystyle=\left(\Delta\sigma_{0}^{\mathrm{stat}}\right)^{2}+\left(\Delta\sigma_{0}^{\mathrm{AC}}\right)^{2}+\left(\Delta\sigma_{0}^{\mathrm{sys}}\right)^{2}$
(12) $\displaystyle=$
$\displaystyle\left[\frac{\Delta\sigma_{\rm{ap}}^{2}}{\sigma_{\rm{ap}}^{2}}+0.03^{2}+\left(\ln\left(\frac{\theta_{\rm{eff}}}{2\theta_{\rm{ap}}}\right)\Delta\eta\right)^{2}\right]\sigma_{0}^{2},$

where $\Delta\sigma_{0}^{\mathrm{stat}}$ is the statistical uncertainty
propagated from the measurement error of $\sigma_{\rm{ap}}$, and
$\Delta\sigma_{0}^{\mathrm{AC}}$ is caused by the uncertainty of $\eta$.

###  2.2 SGL Data Sample

The SGL sample we use in this paper includes 161 galaxy-scale SGL systems
compiled by Chen et al. (2019) from the following surveys: the Lens Structure
and Dynamics (LSD) survey (Koopmans & Treu, 2002, 2003; Treu & Koopmans, 2002,
2004), the Strong Lensing Legacy Survey (Ruff et al., 2011; Sonnenfeld et al.,
2013, 2015), the Sloan Lens ACS survey (Bolton et al., 2008; Auger et al.,
2009, 2010; Shu et al., 2015, 2017), and the Baryon Oscillation Spectroscopic
Survey (BOSS) emission-line lens survey (Brownstein et al., 2012; Shu et al.,
2016a, b), which are selected with strict criteria to ensure the validity of
the assumption of spherical symmetry on the lens galaxy. Specifically, the
lens galaxies should be early-type galaxies with E/S0 morphologies and should
not have significant substructures or close massive companions. The
information of these selected 161 SGL systems could be found in Appendix of
Chen et al. (2019).

As mentioned above, the power-law density profiles expressed by Eq. (4) could
describe the mass distribution of lens galaxy well, in which the slope of the
luminosity density profile $\delta$ for a specific lens galaxy can be obtained
by observation actually, although it is commonly treated as a universal
parameter for all lens galaxies. By fitting the two-dimensional power-law
luminosity profile over a circle of radius $\theta_{\rm{eff}}/2$ for lens
galaxies with the high-resolution Hubble Space Telescope imaging data, the
measured values of $\delta$ for the 130 lens galaxies in the full sample were
obtained (Chen et al., 2019). By constraining cosmological parameters with the
full sample ($\delta$ as a universal parameter) and the truncated sample
($\delta$ as an observable for each SGL system), respectively, Chen et al.
(2019) concluded that the intrinsic scatter of $\delta$ among the lenses
should be taken into account to get an unbiased cosmological estimate. For the
detailed analysis and descriptions, one can refer to Chen et al. (2019) and
references therein.

In this paper, therefore, the SGL sample we adopt is the truncated sample
including 130 SGL systems with the intrinsic scatter of $\delta$, for which
the redshift range of lenses is $0.0624\leq z_{\rm{d}}\leq 0.7224$ and the
redshift range of sources is $0.1970\leq z_{\rm{s}}\leq 2.8324$. However, in
the following, we will use the Gaussian process to reconstruct a distance-
redshift curve directly from the SN Ia data and then to calibrate the
distances of the SGL systems. Therefore, the maximum redshift of the SN Ia
sample determines the number of SGL systems that can be calibrated. Actually,
only 120 SGL data are available in this paper.

###  2.3 Gaussian Process

To constrain on the PPN parameter $\gamma_{\rm{PPN}}$, it is necessary to
determine the distances $D_{\rm{d}}$ and $D_{\rm{ds}}$, which are
conventionally derived from the $\Lambda$CDM model (Cao et al., 2017; Schwab
et al., 2010). However, theoretically, the logic of this approach is
problematic because the $\Lambda$CDM is derived within the framework of GR.
Therefore, in this paper, we will calibrate the distances of SGL via a
cosmological model-independent method, Gaussian Process (GP), which could make
us reconstruct a smooth distance-redshift curve from SN Ia observation without
assuming any parametrization forms. This reconstruction method has been widely
used in cosmology (Wang et al., 2020; Qi et al., 2019a; Seikel et al., 2012a,
b; Wang et al., 2021a; Liao et al., 2019).

In GP, the values of the reconstructed function at different points $z$ and
$\tilde{z}$ are connected by a covariance function $k(z,\tilde{z})$. According
to the analysis in Seikel & Clarkson (2013), the squared exponential form with
the Matérn $(\nu=9/2)$ covariance function can lead to more reliable results
than all others. So we take it here, and its expression is

$\displaystyle k(z,\tilde{z})$ $\displaystyle=$
$\displaystyle\sigma_{f}^{2}\exp(-\frac{3|z-\tilde{z}|}{\ell})$ (13)
$\displaystyle\times$
$\displaystyle(1+\frac{3|z-\tilde{z}|}{\ell}+\frac{27(z-\tilde{z})^{2}}{7\ell^{2}}$
$\displaystyle+$
$\displaystyle\frac{18|z-\tilde{z}|^{3}}{7\ell^{3}}+\frac{27(z-\tilde{z})^{4}}{35\ell^{4}}),$

where $\sigma_{f}$ and $\ell$ are hyperparameters which can be optimized by
the GP itself via the observational data.

Figure 1: Reconstruction of the luminosity distance $D_{L}(z)$ from the
Pantheon sample of SN Ia observation using the Gaussian Process method. Here,
the purple dashed line represents the best fit and the shaded area represents
the 68% confidence level.

###  2.4 Type Ia Supernovae

The SN Ia sample we choose is the Pantheon sample from the Pan-STARRS1 Medium
Deep Survey (Scolnic et al., 2018), which consists of 1048 SN Ia data covering
redshift range $0.001<z<2.3$. This means that only the systems with
$z_{\rm{s}}<2.3$ in the SGL sample can be calibrated by SN Ia. Thus, only 120
SGL systems are actually available.

For SN Ia as standard candle, the observational data are connected with the
distance modulus $\mu_{\rm{obs}}$ by the SALT2 light curve fitter (Guy et al.,
2010),

$\mu_{\mathrm{obs}}(z)=m_{B}(z)+\alpha_{\rm{SN}}\cdot
X_{1}-\beta_{\rm{SN}}\cdot\mathcal{C}-M_{B},$ (14)

where $m_{\rm{B}}$ represents the rest frame B-band peak magnitude, $X_{1}$
and $\mathcal{C}$ denote the time stretch of light curve and the supernova
color at maximum brightness, respectively. $\alpha_{\rm{SN}}$ and
$\beta_{\rm{SN}}$ are two nuisance parameters. The parameter $M_{B}$ is the
absolute B-band magnitude, which is exactly degenerate with $H_{0}$.
Theoretically, if the value of $M_{B}$ or $H_{0}$ is determined, absolute
luminosity distance can be obtained from SN Ia. In the context of Hubble
tension (Riess et al., 2019; Di Valentino et al., 2021; Zhang, 2019; Guo et
al., 2019a, b, 2020; Feng et al., 2020; Gao et al., 2021; Wang et al., 2021b;
Dainotti et al., 2021), the selection of $H_{0}$ often has an effect on the
results (Qi et al., 2021). However, in this paper, the velocity dispersion we
focus is just sensitive to the distance ratio as seen from Eq. (10), which
eliminates the influence of the choice for $H_{0}$.

The stretch-luminosity parameter $\alpha_{\rm{SN}}$ and the color-luminosity
parameter $\beta_{\rm{SN}}$ could be calibrated to zero through a method
called BEAMS with bias corrections (Kessler & Scolnic, 2017). The expression
of the observed distance modulus can be simplified as

$\mu_{\mathrm{obs}}(z)=m_{B}(z)-M_{B}.$ (15)

The theoretical distance modulus is defined as

$\mu_{\mathrm{th}}=5\log_{10}\left(D_{L}/\mathrm{Mpc}\right)+25.$ (16)

The reconstruction of the luminosity distance from the Pantheon sample by
using GP is shown in Figure 1.

###  2.5 The Likelihood Function

By using the emcee Python module (Foreman-Mackey et al., 2013) based on the
Markov Chain Monte Carlo (MCMC) method to maximize the likelihood
$\mathcal{L}\propto e^{-\chi^{2}/2}$, the PPN parameter ($\gamma_{\rm{PPN}}$)
and the parameters of lens models ($\gamma_{0}$, $\gamma_{z}$ and
$\gamma_{s}$) are fitted simultaneously. Specifically, the $\chi^{2}$ function
is defined as

$\chi^{2}\left(\boldsymbol{p},\gamma_{\rm{PPN}}\right)=\sum_{i=1}^{120}\frac{\left(\sigma_{0}^{\mathrm{th}}\left(z_{i};\gamma_{\rm{PPN}}\right)-\sigma_{0}^{\mathrm{obs}}\left(z_{i};\boldsymbol{p}\right)\right)^{2}}{\Delta\sigma_{0}^{\rm{tot}}\left(z_{i}\right)^{2}},$
(17)

where $\boldsymbol{p}$ represents the parameters of lens models. It should be
noted that the uncertainty $\Delta\sigma^{\rm{tot}}_{0}$ is contributed by
both uncertainties from the observation, $\Delta\sigma^{\rm{SGL}}_{0}$ given
by Eq. (12), and from the distance calibration by SN Ia
($\Delta\sigma^{\rm{SN}}_{0}$), and thus we have

$\left(\Delta\sigma^{\rm{tot}}_{0}\right)^{2}=\left(\Delta\sigma^{\rm{SGL}}_{0}\right)^{2}+\left(\Delta\sigma^{\rm{SN}}_{0}\right)^{2}.$
(18)

By using Eq. (10), the fractional uncertainty of $\sigma_{0}$ from the
distance calibration by SN Ia could be estimated as

$\Delta\sigma^{\rm{SN}}_{0}=\sigma_{0}^{\rm{th}}\frac{\Delta
D_{\rm{ratio}}}{2D_{\rm{ratio}}},$ (19)

where $D_{\rm{ratio}}$ is a convenient notation for the ratio of distances in
Eq. (10), i.e., $D_{\rm{ratio}}\equiv D_{\rm{s}}/D_{\rm{ds}}$, and its
uncertainty is $\Delta D_{\rm{ratio}}$. Through the smooth luminosity
distance-redshift curve reconstructed from the SN Ia data, the angular
diameter distance values of $D_{\rm{d}}$, $D_{\rm{s}}$ and their uncertainties
($\Delta D_{\rm{d}}$, $\Delta D_{\rm{s}}$) of each SGL system can be obtained
by using the distance duality relation (Etherington, 1933; Qi et al., 2019b;
Liao et al., 2016). The ratio of distances $\Delta D_{\rm{ratio}}$ expressed
by luminosity distances with superscript $L$ is

$D_{\rm{ratio}}=\frac{1}{1-\frac{D_{\rm{d}}^{L}\left(1+z_{\rm{d}}\right)}{D_{\rm{s}}^{L}\left(1+z_{\rm{s}}\right)}},$
(20)

and its uncertainty is given by

$\displaystyle\Delta D_{\rm{ratio}}^{2}=$ $\displaystyle
D_{\rm{ratio}}^{4}\left(\frac{(1+z_{\rm{d}})}{D_{\rm{s}}^{L}\left(1+z_{\rm{s}}\right)}\Delta
D_{\rm{d}}^{L}\right)^{2}$ (21)
$\displaystyle+D_{\rm{ratio}}^{4}\left(\frac{D_{\rm{d}}^{L}(1+z_{\rm{d}})}{(D_{\rm{s}}^{L})^{2}(1+z_{\rm{s}})}\Delta
D_{\rm{s}}^{L}\right)^{2}.$

In addition, based on the previous studies (Gerhard et al., 2001; Bolton et
al., 2006b; Schwab et al., 2010; Cao et al., 2017; Chen et al., 2019), we will
marginalize the anisotropy parameter $\beta$ using a Gaussian prior with
$\beta=0.18\pm 0.13$ over the range of
$\left[\bar{\beta}-2\sigma_{\beta},\bar{\beta}+2\sigma_{\beta}\right]$, where
$\bar{\beta}=0.18$ and $\sigma_{\beta}=0.13$.

##  3 Results and Discussion

Table 1: Constraint results (68.3% confidence level) for $\gamma_{\rm{PPN}}$ and the SGL parameters ($\gamma_{0}$, $\gamma_{z}$, $\gamma_{\mathrm{s}}$) in different lens models. Lens model | $\gamma_{\rm{PPN}}$ | $\gamma_{0}$ | $\gamma_{z}$ | $\gamma_{\mathrm{s}}$ | $\chi^{2}_{\rm{min}}$ | BIC
---|---|---|---|---|---|---
$P_{1}:\gamma=\gamma_{0}$ | $1.179_{-0.096}^{+0.081}$ | $2.031\pm{0.027}$ | - | - | 334.896 | 349.259
$P_{2}:\gamma=\gamma_{0}+\gamma_{z}z_{d}$ | $1.455_{-0.127}^{+0.154}$ | $2.241\pm{0.039}$ | $-0.629\pm{0.101}$ | - | 292.534 | 311.684
$P_{3}:\gamma=\gamma_{0}+\gamma_{z}z_{d}+\gamma_{\mathrm{s}}\ln\tilde{\Sigma}$ | $1.065_{-0.074}^{+0.064}$ | $1.265\pm{0.103}$ | $-0.246\pm{0.091}$ | $0.281\pm{0.028}$ | 186.196 | 210.133
$P_{3}$ in $\Lambda$CDM | $1.077_{-0.080}^{+0.061}$ | $1.293\pm 0.098$ | $-0.247\pm 0.077$ | $0.272\pm 0.026$ | 188.742 | 213.080
Figure 2: 1D and 2D marginalized probability distributions at 68.3% and 95.4%
confidence level for the PPN parameter $\gamma_{\rm{PPN}}$ and the lens model
parameter $\gamma_{0}$ in the $P_{1}$ model. The dashed line represents
$\gamma_{\rm{PPN}}=1$ predicted by GR. Figure 3: 1D and 2D marginalized
probability distributions at 68.3% and 95.4% confidence level for the PPN
parameter $\gamma_{\rm{PPN}}$ and the lens model parameters ($\gamma_{0}$,
$\gamma_{z}$) in the $P_{2}$ model. Figure 4: 1D and 2D marginalized
probability distributions at 68.3% and 95.4% confidence level for the PPN
parameter $\gamma_{\rm{PPN}}$ and the lens model parameters ($\gamma_{0}$,
$\gamma_{z}$, $\gamma_{\mathrm{s}}$) in the $P_{3}$ model based on two methods
of calibrating distances using the $\Lambda$CDM model and the Gaussian
process, respectively. The dashed line represents $\gamma_{\rm{PPN}}=1$
predicted by GR.

The one-dimensional and two-dimensional marginalized probability distributions
of the PPN parameter and SGL parameters in three lens models are shown in
Figures 2-4, and the fitting values of parameters are presented in Table 1.

For $P_{1}$ model, according to the constraint result,
$\gamma_{\rm{PPN}}=1.179_{-0.096}^{+0.081}$, we find that the inferred value
of the PPN parameter deviates from $\gamma_{\rm{PPN}}=1$ at 1$\sigma$
confidence level, but is still consistent with $\gamma_{\rm{PPN}}=1$ at
2$\sigma$ confidence level, as shown in Figure 2, which means that there is no
significant signal indicating the violation of GR. The constraint accuracy of
$\gamma_{\rm{PPN}}$ is about 7.5%. For the SGL parameter $\gamma_{0}$
reflecting the total mass density distribution of lens galaxy, the constraint
result of $\gamma_{0}=2.031\pm{0.027}$ is consistent with a spherically
symmetric total mass distribution of lens galaxy ($\gamma=2$) which is
predicted by the SIS model.

In Figure 3, we show the constraint results for $P_{2}$ model. It can be seen
clearly that the change of the lens model has a significant impact on the
results. The constraint result of the PPN parameter is
$\gamma_{\rm{PPN}}=1.455^{+0.154}_{-0.127}$ indicating that the GR is excluded
above 2$\sigma$ confidence level. Moreover, the constraint accuracy of
$\gamma_{\rm{PPN}}$ is also reduced to 9.6%. This reduction should not be
caused alone by the addition of one parameter compared to $P_{1}$ model.
Because from the results of $P_{3}$ model with additional three parameters in
Table 1, the constraint accuracy of $\gamma_{\rm{PPN}}$ becomes improved
instead. Therefore, the reduction in constraint accuracy of
$\gamma_{\rm{PPN}}$ here is likely to be caused by the lens model, which will
be discussed in detail in the later analysis of the Bayesian information
criterion. Concerning the parameters of the lens model, the constraint result
of $\gamma_{0}=2.241\pm{0.039}$ indicates that the total mass density
distribution of lens galaxy is not spherically symmetric. In addition, the
parameter $\gamma_{z}$ represents the dependence of the total density profile
of lens galaxy on the redshift. From the result $\gamma_{z}=-0.629\pm{0.101}$,
we find that the mass distribution of lens galaxies, usually early-type
galaxies, is likely evolving with cosmic time.

For $P_{3}$ model, as shown in Figure 4 and Table 1, the inferred value of the
PPN parameter with $\gamma_{\rm{PPN}}=1.065_{-0.074}^{+0.064}$ is in good
agreement with the prediction of $\gamma_{\rm{PPN}}=1$ by GR, and its
constraint accuracy is improved to about 6.4%, although more parameters are
considered in the fitting. For the parameters of lens model, the fit value of
$\gamma_{0}=1.265\pm{0.103}$ is much smaller than those of the above two lens
models, consistent with the previous results of Chen et al. (2019) obtained by
constraining $\Lambda$CDM model with the same SGL sample used here. The other
two parameters of lens model are constrained to be
$\gamma_{z}=-0.246\pm{0.091}$ and $\gamma_{\mathrm{s}}=0.281\pm{0.028}$, both
excluding zero at above the 2$\sigma$ confidence level, indicating that the
dependence of the total density profile on redshift and surface mass density
cannot be ignored.

As we can see, the choice of lens model leads to rather different results for
constraints on the PPN parameter $\gamma_{\rm{PPN}}$, which urges us to
explore which models are better matched to the observational data. In general,
the minimum $\chi^{2}$ comparison is a sufficient tool for comparing different
models with the same number of parameters. A smaller value of
$\chi^{2}_{\rm{min}}$ means a favored model by observations. In Table 1, the
corresponding $\chi^{2}_{\rm{min}}$ for each lens model is listed. However, in
the case of models with different number of parameters, more parameters
usually result in a smaller $\chi^{2}_{\rm{min}}$ for the model. To be fair,
in this work, we also consider using the Bayesian information criterion (BIC)
(Schwarz, 1978) as an evaluation tool to compare different models. The BIC is
expressed as

$\mathrm{BIC}=-2\ln\mathcal{L}_{\max}+k\ln N,$ (22)

where $\mathcal{L}_{\max}$ is the maximum likelihood, $k$ is the number of
model parameters, and $N$ is the number of data points. A model with a smaller
value of BIC is believed to be more favored by observations. In the model
comparison, the absolute value of BIC is not important, and what we are truly
interested in are the relative values of BIC between different models.
Usually, $\Delta\rm{BIC}=2$ is regarded as positive evidence, and
$\Delta\rm{BIC}\geq 6$ is considered to be strong evidence against the model
with higher value of BIC (Mukherjee et al., 1998; Liddle, 2004). The
corresponding BIC value for each lens model is listed in Table 1. It can be
seen that the BIC and $\chi^{2}_{\rm{min}}$ values of $P_{3}$ model are much
smaller than those of the other two models, which strongly indicates that
$P_{3}$ model is the most compatible lens model with the observation.

On the other hand, as mentioned above, to avoid the logical problem caused by
assuming a cosmological model within GR adopted in the previous works
commonly, we use the GP method with the latest Pantheon SN Ia sample to
calibrate the distances of the SGL systems. Here, we wish to investigate
whether using a cosmological model to calibrate distances influences the
inference of $\gamma_{\rm{PPN}}$. For the comparison, in Figure 4 and Table 1,
we present the constraint results of $P_{3}$ model assuming a flat
$\Lambda$CDM model with parameters taken from the Planck 2018 results (Aghanim
et al., 2018). It could be seen clearly that the inferred value of the PPN
parameter with $\gamma_{\rm{PPN}}=1.077_{-0.080}^{+0.061}$ is also consistent
with $\gamma_{\rm{PPN}}=1$ at $1\sigma$ confidence level. Compared with the
results obtained by the GP method, the best-fit value of $\gamma_{\rm{PPN}}$
has a very slight change, and their constraint precisions are almost the same
order, which means that the method of calibrating distances by the
cosmological model does not have a substantial impact on the inference of
$\gamma_{\rm{PPN}}$. For the constraint results of lens model parameters in
the framework of $\Lambda$CDM model, their best-fit values have few changes
compared to the results obtained by the GP method, and their constraint
precision is slightly improved because the distances provided by the
cosmological model have no error, whereas the distances provided by the GP
method contain. Specifically, in Eq. (18), the term of
$\left(\Delta\sigma^{\rm{SN}}_{0}\right)$ from the distance calibration by SN
Ia is absent in the framework of $\Lambda$CDM model. On the other hand, this
also leads to a larger denominator of Eq. (17) in the fitting process for this
method, resulting in larger $\chi^{2}_{\rm{min}}$ and BIC values compared to
that obtained by GP method.

##  4 Conclusion

Strong gravitational lensing provides us with an effective way to test GR by
constraining the PPN parameter at the galactic scale. In this paper, we
consider to improve the constraints on the PPN parameter in the following
three aspects. (i) A larger SGL sample obtained using high-resolution imaging
and stellar dynamical data is considered, in which the luminosity density
slope $\delta$ for each lens is treated as an observational quantity rather
than a universal parameter for all lenses. (ii) To determine the distances in
the SGL sample, we use the GP method to construct a smooth distance-redshift
curve directly from the latest Pantheon SN Ia sample. The logical problem
caused by the introduction of cosmological model in GR can thus be avoided in
such a way of from observation to observation. (iii) The construction of lens
mass model has a great influence on the constraints on the PPN parameter. Here
we consider three commonly used lens models and investigate their impacts on
the fit results. The most favored lens model can be selected by a model
comparison using BIC and $\chi^{2}_{\rm{min}}$.

For the three lens models, we obtain rather different constraint results for
the PPN parameter. First, for $P_{1}$ model, the inferred value of
$\gamma_{\rm{PPN}}$ deviates from $\gamma_{\rm{PPN}}=1$ at 1$\sigma$
confidence, but the result is still consistent with GR at 2$\sigma$ confidence
level. Next, for $P_{2}$ model, GR is excluded at above 2$\sigma$ confidence
level. Finally, for $P_{3}$ model, the constraint result of
$\gamma_{\rm{PPN}}=1.065_{-0.074}^{+0.064}$ is in good agreement with the
prediction of GR. Therefore, we conclude that the choice of lens models has a
great influence on the constraints on the PPN parameter as expected.

Better characterization and modeling for the mass-dynamical structure of
lensing galaxies are necessary for strong gravitational lensing to be used as
an important probe to test the validity of GR. The ongoing and future massive
surveys with high-quality imaging and spectroscopy such as the Dark Energy
Survey (Frieman et al., 2013), the Euclid Survey (Pocino et al., 2021), the
Wide Field Infrared Survey Telescope (Spergel et al., 2015), and the Large
Synoptic Survey Telescope (Abell et al., 2009) will be very helpful in this
study.

To make a model selection by comparing the three lens models, we use the BIC
to assess the fitting results. We find that the BIC value of the $P_{3}$ model
is much smaller than those of the other two models, and thus we conclude that
the $P_{3}$ model is the most compatible model with the observation among the
three lens models, which is also confirmed by another indicator of comparing
the models, $\chi^{2}_{\rm{min}}$.

In previous related works, the distances of SGL systems are generally
calibrated assuming a specific cosmological model within GR, which more or
less leads to a logic problem in testing GR. To investigate the impact of this
problem on the inference of $\gamma_{\rm{PPN}}$, we compare the results in the
$P_{3}$ model based on two methods of calibrating distance using the
$\Lambda$CDM model and the GP method, and we find that using a cosmological
model within GR to calibrate the distances does not have a significant impact
on inferring $\gamma_{\rm{PPN}}$.

The result of $\gamma_{\rm{PPN}}=1.065_{-0.074}^{+0.064}$ in the $P_{3}$ model
is at a 6.4% accuracy, which is the most precise constraint result obtained so
far among the recent tests on GR using the similar method. Here we make a
comparison with some recent representative results. For instance, considering
uncertainties of total mass profile, velocity anisotropy, and luminous matter
profile as the same in this work, but using a different SGL sample, Cao et al.
(2017) obtained the constraint accuracy of the PPN parameter to be 25%. The
nearby SGL system, ESO 325-G004, provided a constraint on the PPN parameter,
$\gamma_{\rm{PPN}}=0.97\pm 0.09$ at 68% confidence level, with a 9.3% accuracy
(Collett et al., 2018). Using four time-delay measurements of SGL from $H_{0}$
Lenses in COSMOGRAIL’s Wellspring collaboration, Yang et al. (2020) obtained a
constraint with an accuracy of 20%, and they forecasted that the joint
constraint from 40 time-delay strong lensing could reach the accuracy of 7.7%
for the PPN parameter in the next few years. Our results show that SGL could
be used as an effective tool to test GR at the galactic scale, and more
promising studies can be expected with the help of more large survey
telescopes in the future.

## Acknowledgments

We thank Yun Chen for very helpful discussions and for sharing the SGL data
with us. We are grateful to the referee for useful comments. This work was
supported by the National Natural Science Foundation of China (Grants Nos.
11975072, 11835009, and 11875102), the Liaoning Revitalization Talents Program
(Grant No. XLYC1905011), the Fundamental Research Funds for the Central
Universities (Grant Nos. N2005030 and N2105014), the National 111 Project of
China (Grant No. B16009), and the China Manned Space Project (Grant No. CMS-
CSST-2021-B01).

## References

  * Abell et al. (2009) Abell, P. A., et al. 2009.  https://arxiv.org/abs/0912.0201
  * Aghanim et al. (2018) Aghanim, N., et al. 2018.  https://arxiv.org/abs/1807.06209
  * Auger et al. (2009) Auger, M. W., Treu, T., Bolton, A. S., et al. 2009, Astrophys. J., 705, 1099, doi: 10.1088/0004-637X/705/2/1099
  * Auger et al. (2010) —. 2010, Astrophys. J., 724, 511, doi: 10.1088/0004-637X/724/1/511
  * Bengochea & Ferraro (2009) Bengochea, G. R., & Ferraro, R. 2009, Phys. Rev. D, 79, 124019, doi: 10.1103/PhysRevD.79.124019
  * Bertotti et al. (2003) Bertotti, B., Iess, L., & Tortora, P. 2003, Nature, 425, 374, doi: 10.1038/nature01997
  * Bolton et al. (2008) Bolton, A. S., Burles, S., Koopmans, L. V. E., et al. 2008, Astrophys. J., 682, 964, doi: 10.1086/589327
  * Bolton et al. (2006a) Bolton, A. S., Burles, S., Koopmans, L. V. E., Treu, T., & Moustakas, L. A. 2006a, Astrophys. J., 638, 703, doi: 10.1086/498884
  * Bolton et al. (2006b) Bolton, A. S., Rappaport, S., & Burles, S. 2006b, Phys. Rev. D, 74, 061501, doi: 10.1103/PhysRevD.74.061501
  * Bolton et al. (2012) Bolton, A. S., Brownstein, J. R., Kochanek, C. S., et al. 2012, Astrophys. J., 757, 82, doi: 10.1088/0004-637X/757/1/82
  * Brownstein et al. (2012) Brownstein, J. R., et al. 2012, Astrophys. J., 744, 41, doi: 10.1088/0004-637X/744/1/41
  * Cao et al. (2016) Cao, S., Biesiada, M., Yao, M., & Zhu, Z.-H. 2016, Mon. Not. Roy. Astron. Soc., 461, 2192, doi: 10.1093/mnras/stw932
  * Cao et al. (2017) Cao, S., Li, X., Biesiada, M., et al. 2017, Astrophys. J., 835, 92, doi: 10.3847/1538-4357/835/1/92
  * Carroll (2001) Carroll, S. M. 2001, Living Rev. Rel., 4, 1, doi: 10.12942/lrr-2001-1
  * Chen et al. (2019) Chen, Y., Li, R., Shu, Y., & Cao, X. 2019, Mon. Not. Roy. Astron. Soc., 488, 3745, doi: 10.1093/mnras/stz1902
  * Collett et al. (2018) Collett, T. E., Oldham, L. J., Smith, R. J., et al. 2018, Science, 360, 1342
  * Copeland et al. (2006) Copeland, E. J., Sami, M., & Tsujikawa, S. 2006, Int. J. Mod. Phys. D, 15, 1753, doi: 10.1142/S021827180600942X
  * Cui et al. (2017) Cui, J.-L., Li, H.-L., & Zhang, X. 2017, Sci. China Phys. Mech. Astron., 60, 080411, doi: 10.1007/s11433-017-9054-5
  * Dainotti et al. (2021) Dainotti, M. G., De Simone, B., Schiavone, T., et al. 2021, Astrophys. J., 912, 150, doi: 10.3847/1538-4357/abeb73
  * Di Valentino et al. (2021) Di Valentino, E., Mena, O., Pan, S., et al. 2021, Class. Quant. Grav., 38, 153001, doi: 10.1088/1361-6382/ac086d
  * Dyson et al. (1920) Dyson, F. W., Eddington, A. S., & Davidson, C. 1920, Phil. Trans. Roy. Soc. Lond. A, 220, 291, doi: 10.1098/rsta.1920.0009
  * Etherington (1933) Etherington, I. 1933, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 15, 761
  * Feng et al. (2020) Feng, L., He, D.-Z., Li, H.-L., Zhang, J.-F., & Zhang, X. 2020, Sci. China Phys. Mech. Astron., 63, 290404, doi: 10.1007/s11433-019-1511-8
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, Publications of the Astronomical Society of the Pacific, 125, 306
  * Frieman et al. (2013) Frieman, J., Collaboration, D. E. S., et al. 2013, in American Astronomical Society Meeting Abstracts# 221, Vol. 221, 335–01
  * Gao et al. (2021) Gao, L.-Y., Zhao, Z.-W., Xue, S.-S., & Zhang, X. 2021, JCAP, 07, 005, doi: 10.1088/1475-7516/2021/07/005
  * Geng et al. (2015) Geng, J.-J., Guo, R.-Y., He, D.-Z., Zhang, J.-F., & Zhang, X. 2015, Front. Phys. (Beijing), 10, 109501, doi: 10.1007/s11467-015-0507-8
  * Gerhard et al. (2001) Gerhard, O., Kronawitter, A., Saglia, R. P., & Bender, R. 2001, Astron. J., 121, 1936, doi: 10.1086/319940
  * Guo et al. (2019a) Guo, R.-Y., Zhang, J.-F., & Zhang, X. 2019a, JCAP, 02, 054, doi: 10.1088/1475-7516/2019/02/054
  * Guo et al. (2020) —. 2020, Sci. China Phys. Mech. Astron., 63, 290406, doi: 10.1007/s11433-019-1514-0
  * Guo et al. (2019b) Guo, R.-Y., Zhang, L., Zhang, J.-F., & Zhang, X. 2019b, Sci. China Phys. Mech. Astron., 62, 30411, doi: 10.1007/s11433-018-9278-1
  * Guy et al. (2010) Guy, J., et al. 2010, Astron. Astrophys., 523, A7, doi: 10.1051/0004-6361/201014468
  * Jiang & Kochanek (2007) Jiang, G., & Kochanek, C. S. 2007, Astrophys. J., 671, 1568, doi: 10.1086/522580
  * Jorgensen et al. (1995) Jorgensen, I., Franx, M., & Kjaergaard, P. 1995, Mon. Not. Roy. Astron. Soc., 276, 1341
  * Kessler & Scolnic (2017) Kessler, R., & Scolnic, D. 2017, Astrophys. J., 836, 56, doi: 10.3847/1538-4357/836/1/56
  * Koopmans & Treu (2002) Koopmans, L. V., & Treu, T. 2002, The Astrophysical Journal Letters, 568, L5
  * Koopmans (2006) Koopmans, L. V. E. 2006, EAS Publ. Ser., 20, 161, doi: 10.1051/eas:2006064
  * Koopmans & Treu (2003) Koopmans, L. V. E., & Treu, T. 2003, Astrophys. J., 583, 606, doi: 10.1086/345423
  * Li et al. (2015) Li, Y.-H., Zhang, J.-F., & Zhang, X. 2015, Phys. Lett. B, 744, 213, doi: 10.1016/j.physletb.2015.03.063
  * Liao et al. (2016) Liao, K., Li, Z., Cao, S., et al. 2016, Astrophys. J., 822, 74, doi: 10.3847/0004-637X/822/2/74
  * Liao et al. (2019) Liao, K., Shafieloo, A., Keeley, R. E., & Linder, E. V. 2019, The Astrophysical Journal Letters, 886, L23
  * Liddle (2004) Liddle, A. R. 2004, Mon. Not. Roy. Astron. Soc., 351, L49, doi: 10.1111/j.1365-2966.2004.08033.x
  * Mukherjee et al. (1998) Mukherjee, S., Feigelson, E. D., Babu, G. J., et al. 1998, Astrophys. J., 508, 314, doi: 10.1086/306386
  * Pan et al. (2021) Pan, Y., He, Y., Qi, J., et al. 2021, Astrophys. J., 911, 135, doi: 10.3847/1538-4357/abebe0
  * Pocino et al. (2021) Pocino, A., et al. 2021.  https://arxiv.org/abs/2104.05698
  * Pound & Rebka (1960) Pound, R. V., & Rebka, Jr., G. A. 1960, Phys. Rev. Lett., 4, 337, doi: 10.1103/PhysRevLett.4.337
  * Qi et al. (2017) Qi, J.-Z., Cao, S., Biesiada, M., Zheng, X., & Zhu, H. 2017, Eur. Phys. J. C, 77, 502, doi: 10.1140/epjc/s10052-017-5069-1
  * Qi et al. (2019a) Qi, J.-Z., Cao, S., Zhang, S., et al. 2019a, Mon. Not. Roy. Astron. Soc., 483, 1104, doi: 10.1093/mnras/sty3175
  * Qi et al. (2019b) Qi, J.-Z., Cao, S., Zheng, C., et al. 2019b, Phys. Rev. D, 99, 063507, doi: 10.1103/PhysRevD.99.063507
  * Qi et al. (2016) Qi, J.-Z., Yang, R.-J., Zhang, M.-J., & Liu, W.-B. 2016, Res. Astron. Astrophys., 16, 022, doi: 10.1088/1674-4527/16/2/022
  * Qi et al. (2021) Qi, J.-Z., Zhao, J.-W., Cao, S., Biesiada, M., & Liu, Y. 2021, Mon. Not. Roy. Astron. Soc., 503, 2179, doi: 10.1093/mnras/stab638
  * Reyes et al. (2010) Reyes, R., Mandelbaum, R., Seljak, U., et al. 2010, Nature, 464, 256, doi: 10.1038/nature08857
  * Riess et al. (2019) Riess, A. G., Casertano, S., Yuan, W., Macri, L. M., & Scolnic, D. 2019, Astrophys. J., 876, 85, doi: 10.3847/1538-4357/ab1422
  * Ruff et al. (2011) Ruff, A. J., Gavazzi, R., Marshall, P. J., et al. 2011, Astrophys. J., 727, 96, doi: 10.1088/0004-637X/727/2/96
  * Schwab et al. (2010) Schwab, J., Bolton, A. S., & Rappaport, S. A. 2010, Astrophys. J., 708, 750, doi: 10.1088/0004-637X/708/1/750
  * Schwarz (1978) Schwarz, G. 1978, Annals Statist., 6, 461
  * Scolnic et al. (2018) Scolnic, D. M., et al. 2018, Astrophys. J., 859, 101, doi: 10.3847/1538-4357/aab9bb
  * Seikel & Clarkson (2013) Seikel, M., & Clarkson, C. 2013.  https://arxiv.org/abs/1311.6678
  * Seikel et al. (2012a) Seikel, M., Clarkson, C., & Smith, M. 2012a, JCAP, 06, 036, doi: 10.1088/1475-7516/2012/06/036
  * Seikel et al. (2012b) Seikel, M., Yahya, S., Maartens, R., & Clarkson, C. 2012b, Phys. Rev. D, 86, 083001, doi: 10.1103/PhysRevD.86.083001
  * Shapiro (1964) Shapiro, I. I. 1964, Phys. Rev. Lett., 13, 789, doi: 10.1103/PhysRevLett.13.789
  * Shu et al. (2015) Shu, Y., et al. 2015, Astrophys. J., 803, 71, doi: 10.1088/0004-637X/803/2/71
  * Shu et al. (2016a) Shu, Y., Bolton, A. S., Kochanek, C. S., et al. 2016a, The Astrophysical Journal, 824, 86, doi: 10.3847/0004-637x/824/2/86
  * Shu et al. (2016b) Shu, Y., Bolton, A. S., Mao, S., et al. 2016b, The Astrophysical Journal, 833, 264
  * Shu et al. (2017) Shu, Y., et al. 2017, Astrophys. J., 851, 48, doi: 10.3847/1538-4357/aa9794
  * Smith (2009) Smith, T. L. 2009.  https://arxiv.org/abs/0907.4829
  * Song et al. (2007) Song, Y.-S., Hu, W., & Sawicki, I. 2007, Phys. Rev. D, 75, 044004, doi: 10.1103/PhysRevD.75.044004
  * Sonnenfeld et al. (2013) Sonnenfeld, A., Treu, T., Gavazzi, R., et al. 2013, Astrophys. J., 777, 98, doi: 10.1088/0004-637X/777/2/98
  * Sonnenfeld et al. (2015) Sonnenfeld, A., Treu, T., Marshall, P. J., et al. 2015, Astrophys. J., 800, 94, doi: 10.1088/0004-637X/800/2/94
  * Sotiriou & Faraoni (2010) Sotiriou, T. P., & Faraoni, V. 2010, Rev. Mod. Phys., 82, 451, doi: 10.1103/RevModPhys.82.451
  * Spergel et al. (2015) Spergel, D., et al. 2015.  https://arxiv.org/abs/1503.03757
  * Taylor et al. (1979) Taylor, J. H., Fowler, L. A., & McCulloch, P. M. 1979, Nature, 277, 437, doi: 10.1038/277437a0
  * Thorne & Will (1971) Thorne, K. S., & Will, C. M. 1971, Astrophys. J., 163, 595, doi: 10.1086/150803
  * Treu & Koopmans (2002) Treu, T., & Koopmans, L. 2002, Astrophys. J., 575, 87, doi: 10.1086/341216
  * Treu & Koopmans (2004) Treu, T., & Koopmans, L. V. E. 2004, Astrophys. J., 611, 739, doi: 10.1086/422245
  * Wang et al. (2020) Wang, B., Qi, J.-Z., Zhang, J.-F., & Zhang, X. 2020, Astrophys. J., 898, 100, doi: 10.3847/1538-4357/ab9b22
  * Wang et al. (2021a) Wang, G.-J., Ma, X.-J., & Xia, J.-Q. 2021a, Mon. Not. Roy. Astron. Soc., 501, 5714, doi: 10.1093/mnras/staa4044
  * Wang et al. (2021b) Wang, L.-F., He, D.-Z., Zhang, J.-F., & Zhang, X. 2021b.  https://arxiv.org/abs/2102.09331
  * Weinberg (1989) Weinberg, S. 1989, Rev. Mod. Phys., 61, 1, doi: 10.1103/RevModPhys.61.1
  * Yang et al. (2020) Yang, T., Birrer, S., & Hu, B. 2020, Mon. Not. Roy. Astron. Soc., 497, L56, doi: 10.1093/mnrasl/slaa107
  * Yang et al. (2019) Yang, T., Hu, B., Cai, R.-G., & Wang, B. 2019, Astrophys. J., 880, 50, doi: 10.3847/1538-4357/ab271e
  * Zhang et al. (2014) Zhang, J.-F., Li, Y.-H., & Zhang, X. 2014, Phys. Lett. B, 739, 102, doi: 10.1016/j.physletb.2014.10.044
  * Zhang et al. (2007) Zhang, P., Liguori, M., Bean, R., & Dodelson, S. 2007, Phys. Rev. Lett., 99, 141302, doi: 10.1103/PhysRevLett.99.141302
  * Zhang (2019) Zhang, X. 2019, Sci. China Phys. Mech. Astron., 62, 110431, doi: 10.1007/s11433-019-9445-7
  * Zhao et al. (2018) Zhao, M.-M., Zhang, J.-F., & Zhang, X. 2018, Phys. Lett. B, 779, 473, doi: 10.1016/j.physletb.2018.02.042