# A Tale of Two Molecules: The Underprediction of CO2 and Overprediction of
PH3 in Late T and Y Dwarf Atmospheric Models

Samuel A. Beiler Ritter Astrophysical Research Center, Department of Physics &
Astronomy, University of Toledo, 2801 W. Bancroft St., Toledo, OH 43606, USA
Sagnick Mukherjee Department of Astronomy and Astrophysics, University of
California, Santa Cruz, 1156 High Street, Santa Cruz, CA 95064, USA
Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD,
USA
Michael C. Cushing Ritter Astrophysical Research Center, Department of Physics
& Astronomy, University of Toledo, 2801 W. Bancroft St., Toledo, OH 43606, USA
J. Davy Kirkpatrick IPAC, Mail Code 100-22, Caltech, 1200 E. California
Boulevard, Pasadena, CA 91125, USA  Adam C. Schneider United States Naval
Observatory, Flagstaff Station, 10391 West Naval Observatory Road, Flagstaff,
AZ 86005, USA  Harshil Kothari Ritter Astrophysical Research Center,
Department of Physics & Astronomy, University of Toledo, 2801 W. Bancroft St.,
Toledo, OH 43606, USA  Mark S. Marley Lunar and Planetary Laboratory,
University of Arizona, 1629 E. University Boulevard, Tucson, AZ 85721, USA
Channon Visscher Department of Chemistry and Planetary Sciences, Dordt
University, Sioux Center, IA, USA  Center for Exoplanetary Systems, Space
Science Institute, Boulder, CO, USA

###### Abstract

The sensitivity and spectral coverage of JWST is enabling us to test our
assumptions of ultracool dwarf atmospheric chemistry, especially with regards
to the abundances of phosphine (PH3) and carbon dioxide (CO2). In this paper,
we use NIRSpec PRISM spectra ($\sim$0.8$-$5.5 $\mu$m, $R\sim$100) of four late
T and Y dwarfs to show that standard substellar atmosphere models have
difficulty replicating the 4.1$-$4.4 $\mu$m wavelength range as they predict
an overabundance of phosphine and an underabundance of carbon dioxide. To help
quantify this discrepancy, we generate a grid of models using PICASO based on
the Elf Owl chemical and temperature profiles where we include the abundances
of these two molecules as parameters. The fits to these PICASO models show a
consistent preference for orders of magnitude higher CO2 abundances and a
reduction in PH3 abundance as compared to the nominal models. This tendency
means that the claimed phosphine detection in UNCOVER$-$BD$-$3 could instead
be explained by a CO2 abundance in excess of standard atmospheric model
predictions; however the signal-to-noise of the spectrum is not high enough to
discriminate between these cases. We discuss atmospheric mechanisms that could
explain the observed underabundance of PH3 and overabundance of CO2, including
a vertical eddy diffusion coefficient ($K_{\mathrm{zz}}$) that varies with
altitude, incorrect chemical pathways, or elements condensing out in forms
such as NH4H2PO4. However, our favored explanation for the required CO2
enhancement is that the quench approximation does not accurately predict the
CO2 abundance, as CO2 remains in chemical equilibrium with CO after CO
quenches.

Brown dwarfs(185), Carbon dioxide(196), Chemical abundances(224), Near
infrared astronomy(1093), Spectroscopy(1558), Y dwarfs(1827), James Webb Space
Telescope(2291)

##  1 Introduction

As we enter the 3rd year since the James Webb Space Telescope’s (JWST) launch
(Rigby et al., 2023), we are starting to gain a new perspective on the
atmospheres of ultracool objects. In particular, we have seen some of the
first JWST observations of the 3.5$-$6 $\mu$m region, where the spectral
energy distributions of late T and Y dwarfs peak (e.g. Burgasser et al., 2023;
Luhman et al., 2023; Lew et al., 2024). This region, which we will refer to as
the 5 $\mu$m peak, is a general feature of cool substellar objects (Marley et
al., 1996). This region contains strong absorption bands of several molecules,
including $\mathrm{H}_{2}\mathrm{O}$, $\mathrm{CH}_{4}$, $\mathrm{CO}$, and
$\mathrm{CO}_{2}$ and has a long pathlength through the atmosphere owing to
low background opacity. This long pathlength allows lower abundance molecules
to be observed in this portion of the spectrum, for example, the detection of
CO, PH3, and GeH4 in Jupiter (Beer, 1975; Ridgway et al., 1976; Fink et al.,
1978). As a result, the spectral morphology of this region is especially
sensitive to variations in atmospheric abundances and chemistry.

Atmospheric mixing is a key factor in the chemistry of these cold atmospheres
(Barshay & Lewis, 1978; Saumon et al., 2000). Mixing can keep the atmosphere
from achieving chemical equilibrium if the mixing timescale is shorter than
the chemical timescale. In late T and Y dwarfs, for example, CH4 is chemically
favored in the upper, cooler layers of atmospheres over CO, but CO is still
observed because it is being dredged up from deeper, hotter layers at a faster
rate than it can be converted to CH4 (Lodders & Fegley, 2002; Visscher &
Moses, 2011). Disequilibrium chemistry has often been observed in brown dwarfs
across several molecules (e.g. Noll et al., 1997; Sorahana & Yamamura, 2012;
Miles et al., 2020), and is an important consideration in the modeling of
their atmospheres. As such, most atmospheric models include disequilibrium
chemistry via vertical mixing, parameterized by a vertical eddy diffusion
coefficient $K_{\mathrm{zz}}$ (e.g. Hubeny & Burrows, 2007; Saumon & Marley,
2008; Phillips et al., 2020; Mukherjee et al., 2022; Lacy & Burrows, 2023),
with a higher value corresponding to stronger mixing.

Models that include disequilibrium chemistry generally fit observations well,
but they still have difficulty fitting the 5 $\mu$m peak (Beiler et al., 2023;
Luhman et al., 2023; Leggett & Tremblin, 2023). This is in part due to the
sheer number of molecular bands present whose abundances impact this portion
of the spectrum. In this paper, we show that the two molecules that are
causing the most problems for the forward model fits are CO2 and PH3. The
models generally predict significant phosphine absorption, and minimal carbon
dioxide absorption, which is the opposite of what has been seen in
observations (Beiler et al., 2023; Luhman et al., 2023). The absence of
observed PH3 absorption in late-type brown dwarfs (with the exception of
Burgasser et al. (2023) claiming a phosphine detection in the galactic disk
object UNCOVER$-$BD$-$3) is especially strange given its strong features in
the spectra Jupiter and Saturn (Ridgway et al., 1976; Prinn & Lewis, 1975;
Gillett & Forrest, 1974).

Previous work has noted that the CO2 absorption band seen in AKARI spectra of
late-type T dwarfs is poorly fit by atmospheric models (Yamamura et al.,
2010). The CO2 absorption band at 4.2 $\mu$m implies an abundance that
disequilibrium models could not replicate while simultaneously fitting the CO
and CH4 features. Tsuji et al. (2011) and Sorahana & Yamamura (2014) both
proposed that enhanced C and O abundances could account for this mismatch.
However, after the completion of the AKARI mission in 2011, the CO2 features
were unobservable until JWST began observation in 2022. As we enter the era of
JWST, understanding the chemistry of CO2 and PH3 is of increased importance.

In this paper, we demonstrate how a variety of standard model grids over-
predict and under-predict PH3 and CO2 abundances, respectively, and explore
how modifying the abundances of PH3 and CO2 affect the model fits. In $\S$2 we
discuss the sample selection and observations, whose reduction we discuss in
$\S$3. In $\S$4 we discuss fitting three commonly used model grids to the
sample and the difficulties in fitting the 5 $\mu$m peak. In $\S$5 we detail
the generation of a grid of models which include parameters to modify the
abundances of phosphine and carbon dioxide, and compare how this grid of
models fit the observations. In $\S$6 we discuss what atmospheric processes
could result in these modified CO2 and PH3 abundances.

##  2 Sample and Observations

Our sample includes two previously published JWST Near Infrared Spectrograph
(NIRSpec, Böker et al., 2023) PRISM spectra of WISE J035934.06$-$540154.6 (Y0,
Beiler et al., 2023) and UNCOVER$-$BD$-$3 (T8 or T9, Burgasser et al., 2023),
as well as two other ultracool dwarfs to more thoroughly probe the spectral
sequence: WISE J210200.15$-$442919.5 (T9) and WISEPC J205628.90+145953.3 (Y0),
hereafter WISE 2102$-$44 and WISE 2056+14, respectively, both observed in GO
Program #2302. The properties for these two objects are in Table 3. Even
though the WISE objects also have MIRI LRS data available, we include only the
NIRSpec PRISM spectra so that a clearer comparison can be made between the
WISE objects and UNCOVER$-$BD$-$3, which only has NIRSpec PRISM data.

We observed WISE 2056+14 and WISE 2102$-$44 using JWST’s NIRSpec, collecting a
low-resolution spectrum for both objects. The instrument was used in fixed-
slit mode with the CLEAR/PRISM filter and the S200A1 slit ($0\farcs 2\times
3\farcs 2$). This filter allowed us to obtain a spectrum from 0.6–5.3 $\mu$m
with a spectral resolving power ($R={\lambda}/{\Delta\lambda}$) of $\sim$100\.
The observations were completed with the NRSRAPID readout pattern and a 5
point dither, with both object having a total exposure time of 38.95 s (5
integrations of 7.79 s).

##  3 Data Reduction

We use the JWST pipeline (Version 1.12.1) for the data reduction, using the
11.17.0 CRDS (Calibration Reference Data System) version and 1132.pmap CRDS
context to assign the reference files. We modified the “extract_1D” step of
the Stage 3 pipeline to extract from the center of the point spread function
as we did in Beiler et al. (2023). We also ignore data shortward of the $Y$
band after the first pixel with negative flux. This occurs at approximately
0.7 $\mu$m for WISE 2102$-$44 and 0.8 $\mu$m for WISE 2056+14.

The pre-flight goal for the precision of the absolute flux calibration of JWST
spectra was $\sim 10$% (Gordon et al., 2022). We improved the absolute
calibration by using Spitzer/IRAC Channel 2 ([4.5], 4.5 $\mu$m) photometry
(Kirkpatrick et al., 2012), which can be found in Table 3. From these flux
densities we calculate the scaling factor needed to convert our spectra and
their errors to absolute units of Janskys (Reach et al., 2005; Cushing et al.,
2006).

Table 1: Astrometric and Photometric Properties

Property | WISE 2102$-$44 | WISE 2056+14 | Ref.
---|---|---|---
Spectral Type | T9 | Y0 | 1,2
Parallax $\varpi_{\mathrm{abs}}$ (mas) | $92.9\pm 1.9$ | $140.8\pm 2.0$ | 3,4
Distance (pc) | $10.76\pm 0.22$ | $7.1\pm 0.1$ | 3,4
IRAC [3.6] (mag) | $16.325\pm 0.036$ | $16.031\pm 0.030$ | 3,3
IRAC [4.5] (mag) | $14.223\pm 0.019$ | $13.923\pm 0.018$ | 3,3

References. — (1) Kirkpatrick et al. (2012),(2) Cushing et al. (2011), (3)
Kirkpatrick et al. (2019), (4) Tinney et al. (2014).

##  4 Comparing the Spectra to Standard Atmospheric Models

We fit this sample of NIRSpec spectra with a variety of standard atmospheric
models to demonstrate the difficulty in replicating the spectra at the 5
$\mu$m peak. We used the Elf Owl (Mukherjee et al., 2024), ATMO2020++
(Phillips et al., 2020; Meisner et al., 2023), and LOW-Z models (Meisner et
al., 2021), all of which include disequilibrium chemistry. The Elf Owl and
LOW-Z models parameterize this with a constant vertical eddy diffusion
coefficient ($K_{\mathrm{zz}}$). The ATMO2020++ models do not have an
independent $K_{\mathrm{zz}}$ parameter, but instead defines a vertically
uniform $K_{\mathrm{zz}}$ as a function of log $g$ [$\mathrm{cm~{}s}^{-2}$]
such that log $K_{\mathrm{zz}}$$=5+2\times(5-$log $g$) [$\mathrm{cm}^{2}\
\mathrm{s}^{-1}$]. The parameter ranges we use for each model grid can be
found in Table 2. The increments on the log $g$ and $T_{\mathrm{eff}}$
parameters change both between and within models. The increments for the Elf
Owl models can be found in Mukherjee et al. (2024) (Table 1). The ATMO2020++
model grid increments $T_{\mathrm{eff}}$ at 25 K from 250-300 K, 50 K from
300-500 K, and 100 K elsewhere, while the LOWZ model grid increments at 50 K
below 1000 K, and 100 K above 1000 K. For log $g$, the ATMO2020++ and LOWZ
model grids both increment at 0.5 dex, with the exception of LOWZ including
models where log $g$=5.25 [cm s-2].

Table 2: Atmospheric Model Parameter Ranges Model Name |  $T_{\mathrm{eff}}$ (K) |  log $g$ [cm s-2] | log $K_{\mathrm{zz}}$ [cm2 s-1] | [M/H] | C/O
---|---|---|---|---|---
Sonora Elf Owl | 2400$-$275 | 3.25$-$5.5 | 2, 4, 7, 8, 9 |  $-$1.0, $-$0.5, 0.0, 0.5, 0.7, 1.0 | 0.23, 0.46, 0.92, 1.15 aaThese values are often labeled as (C/O)/(C/$\mathrm{O}_{\odot}$), where C/O⊙ = 0.458 (i.e. 0.5, 1, 2, and 2.5), but we have converted them for ease of comparison.
ATMO2020++ | 1200$-$250 | 2.5$-$5.5 |  $5+2\times(5-$log $g$) |  $-$1.0, $-$0.5, 0.0 | 1
LOWZ | 1600$-$500 | 3.5$-$5.25 | 0, 2, 10 |  $-$2.5 to 1.0 | 0.1, 0.55, 0.85

The published version of the Elf Owl model grid handles the phosphine
abundance differently than the model’s general disequilibrium chemistry
scheme. Instead PH3 is assumed to be in chemical equilibrium, due to the non-
detection of phosphine in most brown dwarf atmospheres. Since we are
interested in the current disequilibrium chemistry paradigm and its prediction
of phosphine abundances, we use an earlier version of the Elf Owl models which
uses a phosphine abundance consistent with the general disequilibrium
chemistry scheme. This model gird is publicly available at
https://doi.org/10.5281/zenodo.11370830.

To compare the models to the observed spectra, we first need to convolve the
model spectra to the resolving power of the NIRSpec PRISM spectra at each
wavelength due to the non-uniform resolving power of the spectrum, which we
calculated as detailed in Beiler et al. (2023). The best-fit model for each
grid is the one that minimizes $\chi^{2}$, defined as:

$\chi^{2}=\sum_{i}\left(\frac{f_{\lambda,i}-CM_{\lambda,i}}{\sigma_{\lambda,i}^{2}}\right)^{2},$
(1)

where $M_{\lambda,i}$, $f_{\lambda,i}$, and $\sigma_{\lambda,i}$ are
respectively the model flux density, observed flux density, and the observed
flux density uncertainty in Janskys at each data point $i$. $C$ is the scaling
factor for the model spectrum that minimizes $\chi^{2}$, defined by:

$C=\frac{\sum\limits_{i}M_{\lambda,i}f_{\lambda,i}/\sigma_{\lambda,i}^{2}}{\sum\limits_{i}{M_{\lambda,i}^{2}/\sigma_{\lambda,i}^{2}}}.$
(2)

The best-fit model parameters for these three nominal model grids are reported
in Table 3, and Figure 1 shows the best-fit model spectra plotted against all
four objects. For UNCOVER–BD–3, we find a different best-fit LOWZ model than
Burgasser et al. (2023) (which is $T_{\mathrm{eff}}$=550 K, log $g$=5.25,
[M/H] = $-$0.5), but the $\Delta\chi^{2}$ between the two best-fit models is
only 17 and is likely caused by our choice to convolve the model to NIRSpec’s
non-uniform resolving power instead of assuming a uniform resolving power. The
best-fit models of the three model grids provide reasonable fits to these four
objects, with similar temperatures and surface gravities to the other models.
The range where the models do have difficulty replicating these observations
is from 4–4.5 $\mu$m due to predicting strong absorption (see Figure 2, a
zoomed in version of Figure 1).

Table 3: Best-Fit Nominal Atmospheric Model Parameters  Object Name | Model Name |  $T_{\mathrm{eff}}$ (K) |  log $g$ [cm s-2] | log $K_{\mathrm{zz}}$ [cm2 s-1] | [M/H] | C/O | $\chi^{2}$ | $\chi^{2}_{r}$
---|---|---|---|---|---|---|---|---
WISE 2102$-$44 | Sonora Elf Owl | 600 | 3.75 | 2 | 0.0 | 0.23 aaThese models report C/O as (C/O)/(C/$\mathrm{O}_{\odot}$), where C/O⊙ = 0.458, but we have converted them for ease of comparison. | 5206.31 | 13.38 |  |
| ATMO2020++ | 600 | 4.5 | 6 | 0.0 | 0.5 | 10197.29 | 26.21 |  |
| LOWZ | 500 | 4.5 | 10 |  $-$0.25 | 0.85 | 19695.41 | 50.63 |  |
WISE 2056+14 | Sonora Elf Owl | 475 | 3.25 | 4 | 0.0 | 0.23 aaThese models report C/O as (C/O)/(C/$\mathrm{O}_{\odot}$), where C/O⊙ = 0.458, but we have converted them for ease of comparison. | 11246.58 | 30.07 |  |
| ATMO2020++ | 450 | 4.0 | 7 | 0.0 | 0.5 | 3881.94 | 10.38 |  |
| LOWZ | 500 | 3.5 | 10 |  $-$1.0 | 0.85 | 16547.60 | 44.24 |  |
WISE 0359$-$54 | Sonora Elf Owl | 450 | 4.5 | 9 |  $-$1.0 | 1.15 aaThese models report C/O as (C/O)/(C/$\mathrm{O}_{\odot}$), where C/O⊙ = 0.458, but we have converted them for ease of comparison. | 4220.04 | 12.34 |  |
| ATMO2020++ | 450 | 3.0 | 9 |  $-$1.0 | 0.5 | 7661.16 | 22.04 |  |
| LOWZ | 500 | 4.5 | 10 |  $-$1.0 | 0.85 | 10455.49 | 30.57 |  |
UNCOVER$-$BD$-$3 | Sonora Elf Owl | 550 | 4.25 | 7 |  $-$1.0 | 1.15 aaThese models report C/O as (C/O)/(C/$\mathrm{O}_{\odot}$), where C/O⊙ = 0.458, but we have converted them for ease of comparison. | 749.66 | 1.86 |  |
| ATMO2020++ | 500 | 4.5 | 6 | 0.0 | 0.5 | 764.05 | 1.89 |  |
| LOWZ | 500 | 4.5 | 10 |  $-$0.5 | 0.85 | 814.40 | 2.02 |  |
Figure 1: The four objects in this sample (black) plotted alongside the best-
fit model for each of the three nominal model grids: ATMO2020++ (blue), Sonora
Elf Owl (purple), LOWZ (green). The models generally fit well, but fit the
4–4.5 $\mu$m region poorly where PH3 and CO2 affect the spectrum.
Uncertainties on the data are shown as grey bars on the data. The two newly
published spectra shown in this figure are available as the Data behind the
Figure. Figure 2: The same plot as Figure 1, but zoomed in to clearly show the
4–4.5 $\mu$m region where PH3 and CO2 affect the spectrum. For all of the
models, the PH3 absorption is too strong in at least one of the $Q$ and $R$
branches, and the CO2 is generally too weak.

To get a sense are what is causing the poor fits to this region, in Figure 3
we plot WISE 2056+14 alongside the absorption coefficients
($\alpha$=$n\sigma$, where $n$ is the number density and $\sigma_{i}$ is the
cross section) for important molecules (H2O, CH4, CO, NH3, PH3, H2S, and CO2).
The number density are from the best-fit Sonora Elf Owl model at $P$ = 1 bar
and $T$ = 745 K, and the cross sections are from PICASO at $P$ = 1 bar and the
closest temperature available, $T$ = 725 K.

Figure 3: Top: The absorption coefficents ($\alpha_{i}$=$n_{i}\sigma_{i}$) for
important molecules in the 3.5–5.3 $\mu$m range for important molecules (H2O,
CH4, CO, NH3, PH3, H2S, and CO2). The abundances ($n_{i}$) are from the best-
fit Sonora Elf Owl model for WISE 2056+14 at P = 1 bar. The cross-sections
($\sigma_{i}$) are from PICASO at P = 1 bar and T = 725 K, the closest
temperature to the best-fit Elf Owl model at 1 bar (T = 745 K). Note the
y-axis is inverted for ease of comparison to the spectrum. Bottom: A zoom-in
on the 3.4–5.3 $\mu$m region of the WISE 2056+14 portion of Figure 1, showing
clearly how the best-fit Elf Owl model does not fit the CO2 and PH3 features
of the observed spectrum.

The 5 $\mu$m peak is predominately shaped by the $\nu_{2}$ water band, the
$\nu_{3}$ methane band, and the 1-0 X${}^{1}\Sigma^{+}$-X${}^{1}\Sigma^{+}$
band of carbon monoxide. These molecular features are as a whole well fit with
some discrepancy in the height of the water features past 4.8 $\mu$m. However,
the worst fitting region is 4–4.5 $\mu$m, where $\nu_{3}$ CO2 band and the
overlapping $\nu_{1}$ and $\nu_{3}$ PH3 bands play a strong role in shaping
the spectra. It is clear from these fits that the models have an overabundance
of phosphine, seen most clearly in the deep $Q$-branch of phosphine at 4.3
$\mu$m, but also in the $P$\- and $R$-branches that bracket the $Q$-branch.

This overabundance of phosphine arises as an effect of the non-zero
$K_{\mathrm{zz}}$ values that are required to fit the prominent CO and
$\nu_{3}$ CH4 band at effective temperatures below $\sim$600 K. This
enhancement of CO has been seen many times in the spectra of L, T and Y dwarfs
and even Jupiter (e.g. Beer, 1975; Noll et al., 1997; Oppenheimer et al.,
1998; Sorahana & Yamamura, 2012; Miles et al., 2020). For phosphorus, vertical
mixing enhances PH3 over oxidized P-bearing species expected to dominate in
chemical equilibrium at these temperatures (Wang et al., 2016; Visscher,
2020). PH3 absorption has been seen in the spectra of Jupiter and Saturn
(Gillett & Forrest, 1974; Prinn & Lewis, 1975; Ridgway et al., 1976), but has
not been found outside of our solar system. For example, in the best-fit
parameters for WISE 2056+14 the phosphine abundance at 1 bar increases by
$\approx$100 when going from log $K_{\mathrm{zz}}$ = 0 to log
$K_{\mathrm{zz}}$ = 4 [cm2 s-1] (Mukherjee et al., 2023b). Unlike CO, PH3 is
not expected to appear in the spectra of hotter L and T dwarfs due to the
preference for other species such as P2 in the low pressure regions of these
objects (Visscher et al., 2006).

Even though the large $K_{\mathrm{zz}}$ values that are selected to fit the CO
spectral feature result in an overabundance of PH3, these actually improve the
$\chi^{2}$ compared to low $K_{\mathrm{zz}}$ values by compensating for the
underabundance of CO2 in the models. Like PH3 and CO, the CO2 abundance is set
by the $K_{\mathrm{zz}}$ parameter, which dredges up CO2 from deeper
atmospheric layers. In the case of the Elf Owl best-fit model for WISE
2056+14, the log $K_{\mathrm{zz}}$ = 4 [cm2 s-1] increases the CO and CO2
abundance over chemical equilibrium by a factor of $\approx 10^{8}$ and
$10^{6}$, respectively. However, this increased CO2 abundance is still too low
in models to match our sample.

The lone exception to the standard models over-abundance of phosphine in this
sample is the best-fitting Elf Owl model of WISE 2102$-$44\. A combination of
a low $K_{\mathrm{zz}}$, higher $T_{\mathrm{eff}}$, and low C/O ratio are
required to fit the carbon monoxide and methane features, and this represses
the phosphine in the model. However, even in this case we see the
underabundance of carbon dioxide in the model.

##  5 Scaling PH3 and CO2 Abundances with PICASO

To test the extent to which the phosphine and carbon dioxide abundances are
creating the 4–4.5 $\mu$m discrepancy, we generate a grid of models where we
include the abundances of these molecules as parameters. PH3 and CO2 are both
trace species, and so changing their abundances will not have a significant
effect on the pressure/temperature $(P/T)$ profile. As such, we generate our
new grid by taking $P/T$ and chemical profiles from the Elf Owl models,
multiplying the PH3 and CO2 abundance profiles by various factors, and using
the PICASO framework (v3.0, Batalha et al., 2019; Marley et al., 2021;
Mukherjee et al., 2023a) to generate emission spectra. We use the same
opacities as Mukherjee et al. (2023b) (Table 2) at a native spectral
resolution of $R$=60,000 (Batalha et al., 2020). We will refer to this grid as
the PICASO grid hereafter. The parameter range of the grid can be seen in
Table 4. Adding two parameters greatly increases the number of models, so we
reduce the run of parameters compared to the Elf Owl range. In total, we
generated 86,400 models. These models, convolved to the resolution of our
objects, can be found at https://doi.org/10.5281/zenodo.11370830.

Table 4: PICASO Atmospheric Model Parameter Range  Parameter | Ranges
---|---
$T_{\mathrm{eff}}$ (K) | 400–600 (25 K steps)
log $g$ [cm s-2] | 3.25–5.5 (0.25 steps)
log $K_{\mathrm{zz}}$ [cm2 s-1] | 2, 4, 7, 9
[M/H] |  $-$1.0, 0.0
C/O | 0.23, 0.46, 1.15 aaThese values are often labeled as (C/O)/(C/$\mathrm{O}_{\odot}$), where C/$\mathrm{O}_{\odot}$ = 0.458, but we have converted them for ease of comparison.
$\times$ CO2 | 1, 50, and 10n and 5$\times$10n for n = 2–4
$\times$ PH3 | 1, 0.5, 0.2, 0.1, 0.01

We follow the same fitting procedure as Section 4, resulting in the best-fit
parameters listed in Table 5 and the best-fit PICASO model spectrum plotted in
Figure 4 alongside the Elf Owl best-fitting model spectrum. The PICASO models
are a significant improvement over the Elf Owl models in every case. The
PICASO grid is not the preferred model for every object, as WISE 2056+14 is
better fit by a ATMO2020++ model, but the PICASO grid is a better fit than the
other standard models for the rest of the sample.

Table 5: Best-Fit PICASO Atmospheric Model Parameters  Object Name |  $T_{\mathrm{eff}}$ (K) |  log $g$ [cm s-2] | log $K_{\mathrm{zz}}$ [cm2 s-1] | [M/H] | C/O |  $\times$ CO2 |  $\times$ PH3 | $\chi^{2}$ | $\chi^{2}_{r}$
---|---|---|---|---|---|---|---|---|---
WISE 2102$-$44 | 600 | 3.5 | 2 | 0.0 | 0.23aaThese values are often labeled as (C/O)/(C/O⊙) where C/O⊙ = 0.458, but we have converted them for ease of comparison. | 100 | 1 | 4504.33 | 11.58 |
WISE 2056+14 | 475 | 4  | 9 |  $-$1.0 | 1.15aaThese values are often labeled as (C/O)/(C/O⊙) where C/O⊙ = 0.458, but we have converted them for ease of comparison. | 5000 | 0.5 | 6457.95 | 17.27 |
WISE 0359$-$54 | 425 | 3.25 | 9 |  $-$1.0 | 0.46aaThese values are often labeled as (C/O)/(C/O⊙) where C/O⊙ = 0.458, but we have converted them for ease of comparison. | 1000 | 0.01 | 1726.82 | 5.05 |
UNCOVER$-$BD$-$3 | 525 | 3.75 | 7 |  $-$1.0 | 1.15aaThese values are often labeled as (C/O)/(C/O⊙) where C/O⊙ = 0.458, but we have converted them for ease of comparison. | 5000 | 0.5 | 712.59 | 1.76 |
Figure 4: The best-fit Elf Owl (purple) and PICASO (blue) models for the
objects in our sample. The PICASO grid, which parameterizes the CO2 and PH3
abundances, grid provides the most consistently good fit across all four
objects with a consistent preference for an elevated CO2 abundance, and a
decreased PH3 abundance.

For all of the objects, there is a strong preference for elevated CO2
abundances, up to a factor of 5000 above the base Elf Owl model with the same
parameters. This abundance is high enough that the $\nu_{3}$ CO2 band is now
one of the dominant absorption features at 4.25 $\mu$m matching observations,
which can be seen in Figure 5 and 6. Similar figures for the other two objects
can be found in the Appendix. With such a large increase in CO2 abundance, it
might be expected that the C/O ratio would be impacted, however for the best-
fit models it changes by less than 0.1% due to CO2 being at such low overall
abundances.

Figure 5: Top: The same plot as Figure 3, but this time using the abundances
from the PICASO best-fit model at P = 1 bar, where T = 469 K, and the cross-
sections at T = 475 K. Note the dramatic increase in the CO2 absorption
coefficent. Bottom: The best-fit Elf Owl (purple) and PICASO(blue) models from
3.4–5.3 $\mu$m. Note the particular improvement of the fit in the range that
CO2 and PH3 dominate (4.1–4.4 $\mu$m). Figure 6: The same plot as Figure 5 but
with the fits for UNCOVER–BD–3. The abundances are from the PICASO best-fit
model at P = 1 bar, where T = 568 K, and the cross sections at T = 575 K.

We do note that all of the objects which require such high CO2 increases are
also best fit by low metallicity PICASO models ([M/H] = $-$1.0). Of these
models, only WISE 2056+14 has solar metallicity models in their top 20 best-
fitting models, almost half of which are solar metallicity. To achieve a good
fit, these PICASO solar metallicity models all require lower
$T_{\mathrm{eff}}$ (475 or 450 K) and log $K_{\mathrm{zz}}$ (4 [cm2 s-1])
values to try and balance the CO and CH4 abundances, however these models are
unfavored as the CH4 absorption is still too strong. All of these factors lead
to solar abundance models only requiring a 100-500$\times$ increase of CO2,
compared to the 5000$\times$ increase of the sub-solar models.

All of these objects are best fit with PICASO models whose log $g$$\leq 4$,
implying the somewhat unlikely case that they all have low masses (M $<$ 5
$M_{\mathrm{Jup}}$) and young ages ($\tau<$ 0.3 Gyr) Marley et al. (2021).
None of the objects have a model with log $g$$\geq 5$ in their 250 best-fit
models. This tendency is inherited from the Elf Owl models, whose preference
for low gravity fits is seen in this paper, and was also seen in fitting the
full spectrum of WISE 0359$-$54 to Elf Owl models (Beiler et al., 2023).
Kothari et al. (2024) showed that the retrieved gravity and temperature for
WISE 0359$-$54 is more in line with the 1-10 Gyr age estimate for field
objects, and as such we do not believe the log $g$ values presented here to be
accurate. Fortunately, the best-fitting models with log $g$$>4$ have similar
parameters to the overall best-fit models, so while the preference for low
gravities is interesting, it does not impact our conclusions.

The PICASO model fits also show a preference for decreasing the abundance of
phosphine by a factor of 2-100$\times$ in three of our four objects. The
exception is WISE 2102$-$44, whose phosphine abundance was already too low to
be seen in the Elf Owl model for the reasons explained above (high
$T_{\mathrm{eff}}$, low $K_{\mathrm{zz}}$). The change in abundance is not as
large for phosphine as it is for carbon dioxide, but it still has a
significant impact on the best-fit model spectra. At these lower abundances,
the $P$\- and $R$-branches of phosphine are no longer discernible in the model
spectra. However, while the 4.25-4.4 $\mu$m region where the PH3 $Q$-branch
sits is now well fit in WISE 0359$-$54 and UNCOVER–BD–3, the best-fit mode1
for WISE 2056+14 still shows too much PH3 absorption. That said, this region
is better reproduced than in the Elf Owl model even with the large increase in
$K_{\mathrm{zz}}$ that improves the fit redward of 4.6 $\mu$m. It might be
expected that a further reduction in PH3 would result in a better fit, but
without the $P$\- and $R$-branches acting to suppress the 4 and 4.4 $\mu$m
peaks there is too much flux at these peaks to make it the best-fitting model.
The 2nd best fit model is identical to the best-fit model, except with
$\times$0.2 PH3 instead of $\times$0.5, so the true best fit is likely
somewhere between these two values.

There is a further complication with identifying a PH3 absorption feature at
4.3 $\mu$m. Methane and water both have local peaks that are aligned at this
wavelength, which when combined can sometimes be equivalent to or dominate
over phosphine in the standard models. This overlap increases the difficulty
in identifying phosphine in these objects. With this in mind, it is clear
there is no detection of phosphine in the three WISE objects.

It is more difficult to rule out a phosphine detection for UNCOVER–BD–3.
Burgasser et al. (2023) claims the absorption at 4.2 $\mu$m is due to the
$R$-branch of phosphine. This claim was made assuming the abundances of their
LOWZ best-fit model was correct, except for CO2 where they adopted
$N_{\mathrm{CO_{2}}}/N_{\mathrm{CO}}=10^{-3}$ from Yamamura et al. (2010),
which is based off model fits to AKARI spectra (2.6-5 $\mu$m) of mid- and
late-T dwarfs. However, we have shown that standard atmospheric models over-
predict the phosphine abundance, and the best-fit PICASO model gives
$N_{\mathrm{CO_{2}}}/N_{\mathrm{CO}}=3.8\times 10^{-3}$ at 1 bar (T = 567 K).
These factors combined to cause CO2 to dominate over PH3 for the best-fit
PICASO model. Yamamura et al. (2010) actually estimated
$N_{CO_{2}}/N_{CO}\approx 10^{-3}$ as a lower bound, which would be consistent
with our findings and points to enhanced CO2 as another possible explanation
for the 4.2 $\mu$m absorption feature.

To further demonstrate the lack of clarity in the phosphine abundance in
UNCOVER–BD–3, we can look at the 15 best-fitting PICASO models. Of these
models there is a roughly even distribution across the phosphine parameter:
three models have no change to the phosphine, six have phosphine halved, four
have it reduced by a factor of 5, and two have it reduced by a factor of 10.
All of these models have roughly the same temperature (525 or 550 K), gravity
(3.25$-$4.25 [cm s-2]), and CO2 factor (5000$\times$ or 10000$\times$), and
all the other parameters are constant. At the signal-to-noise of this spectrum
($\sim 8$) it is not possible to definitively determine the abundance of
phosphine from models.

With a higher signal-to-noise spectrum of UNCOVER–BD–3 it might be possible to
confirm a detection of phosphine; it is the only ultracool object with a
published JWST NIRSpec spectrum that potentially has observable phosphine.
However, UNCOVER–BD–3’s great distance from Earth makes the exposure time for
high signal-to-noise spectrum prohibitive. Burgasser et al. (2023) predicts
that phosphine might be a hallmark of sub-solar metallicity, given that
UNCOVER–BD–3 is a disk object, making nearby low-metallicity objects a
possibility for confirmation. However, WISE 2056+14 and WISE 0359$-$54 are
best-fit by models with similarly low metallicities, and both lack phosphine
absorption features. A JWST Cycle 3 proposal (GO #4668, PI A. Burgasser) has
been accepted to observe 32 L and T subdwarfs, which will help determine if
phosphine absorption is a hallmark of sub-solar metallicity.

##  6 Possible Explanations for These Differences in Abundances

There are several possible explanations for the overabundance of phosphine and
underabundance of carbon dioxide in standard atmospheric models. For CO2, one
possibility is that vertical mixing is function of altitude, as all the models
in this paper assume $K_{\mathrm{zz}}$ is constant with altitude. The deep
atmosphere $K_{\mathrm{zz}}$ is expected to set the quench pressures for CH4,
CO, and NH3, whereas the upper atmosphere $K_{\mathrm{zz}}$ should set the
quenched abundance of CO2 (Zahnle & Marley, 2014; Mukherjee et al., 2022,
2023b; Phillips et al., 2020). We perform a test of the effect of
$K_{\mathrm{zz}}$ varying with altitude where we take a self-consistent Elf
Owl model for the best-fit set of parameters and see if a higher or lower
$K_{\mathrm{zz}}$ in the upper atmosphere can explain the difference in CO2
abundances. We do this test in an ad-hoc manner, where we use the $P/T$
profile of the Elf Owl model, instead of a fully self-consistent approach. We
tried to enhance and decrease $K_{\mathrm{zz}}$ in the upper atmosphere by a
factor relative to the deeper atmosphere $K_{\mathrm{zz}}$ and find that
changes in the upper atmosphere $K_{\mathrm{zz}}$ cannot by themselves explain
the much larger enhancements in CO2 abundance required to explain the data.

Another possibility is that our disequilibrium approximations are insufficient
for CO2. In the Elf Owl models, the pressure at which a molecule quenches is
set by its abundance at the pressure where its chemical timescale becomes
slower than the atmospheric mixing timescale. This approximation works well
for most molecules, but does not account for the fact that CO quenches at a
higher pressure than CO2. After CO quenches with respect to CH4 (by quenching
of CO$\rightleftharpoons$CH4 reaction pathways), CO2 remains in equilibrium
with the quenched CO abundance due to relatively faster
CO$\rightleftharpoons$CO2 reaction pathways. This results in the CO2 abundance
increasing with altitude above the CO quench level, until CO2 itself quenches
to a fixed abundance (Visscher et al., 2010).

To demonstrate this behavior, we model the CO, CO2, and CH4 abundances with
the Elf Owl framework assuming chemical equilibrium, a quench approximation,
and a full 1D chemical kinetics model (Photochem, Wogan et al., 2024) with
parameters of 500 K, log $g$=3.5 [cm s-1], [M/H]=$-1.0$, and solar C/O, which
roughly represents the best-fitting PICASO models for this sample. These
abundance profiles can be seen in Figure 7. The kinetics and quench
approximation models have log $K_{\mathrm{zz}}$= 8 [cm2 s-1]. With these
parameters, we see a $\sim$200 fold increase in CO2 in the kinetics model
(solid line) compared to the quench approximation model (dotted line). While
this does not represent the full enhancement ($\sim$1000$\times$) we see in
the PICASO best-fit models, it is meaningful increase and is therefore a
strong candidate for the mechanism responsible for the CO2 enhancement needed
to fit our spectra.

Figure 7: The chemical profiles of CO2 (red), CO (pink), and CH4 (green)
assuming various chemical schemes within the Elf Owl framework at 500 K, log
$g$ = 4.5 [cm s-1], [M/H]=$-1.0$, and solar C/O ratio. The dashed lines are
from a model assuming equilibrium chemistry (meaning $K_{\mathrm{zz}}$=0), the
dotted lines are from a model with disequilibrium chemistry implemented via a
quenching timescale approximation, and the solid lines are from a kinetics
model. Both nonequilibrium models have a log $K_{\mathrm{zz}}$ of 8 [cm2 s-1].
The kinetic model enhances CO2 by a factor of $\sim$200 over the base Elf Owl
disequilibrium models, due to CO$\rightleftharpoons$CO2 reactions keeping CO2
in equilibrium with the quenched CO. This enhancement factor is not as large
as the best-fit PICASO enhancement of CO2 ($\sim$1000$\times$), but does move
the CO2 meaningfully closer to the preferred value. Due to CO2 being a trace
species, CO and CH4 are not strongly impacted by this change.

For phosphorus, it could be that our understanding of the chemical pathways
are incomplete. If the chemical timescale for PH3 to be removed by oxidation
or condensation (e.g. P4O6, P4O10, H3PO4, NH4H2PO4) is shorter than is
currently believed (Visscher et al., 2006; Wang et al., 2016; Morley et al.,
2018; Visscher, 2020), the appearance of PH3 absorption in model spectra would
begin at higher $K_{\mathrm{zz}}$’s and/or lower $T_{\mathrm{eff}}$’s. Even in
WISE 2102$-$44, the atmospheric model is sufficiently hot that with a low
$K_{\mathrm{zz}}$, PH3 absorption is weak and does not need to be reduced, so
a small change in the phosphorus chemistry could cause the PH3 features to
disappear.

At a more basic level, it is possible there is no longer enough phosphorus in
gaseous form to create the phosphine absorption feature. The condensate
NH4H2PO4 could act as sink for phosphorus, as discussed in Visscher et al.
(2006) and Morley et al. (2018). The condensation curve of NH4H2PO4, assuming
it condenses directly from PH3, crosses the PICASO $P/T$ profiles between 1
and 0.5 bar for all our objects except for WISE 2102$-$44, where it crosses at
0.1 bar. This would be consistent with out observation that the PH3 mismatch
gets worse for colder objects where NH4H2PO4 is more likely to condense out.
As an aside, the clear detection of PH3 on Jupiter and Saturn makes the P
chemistry on warmer objects that much more enigmatic. It is possible that the
super-solar metallicities of these gas giants (Mahaffy et al., 2000; Orton et
al., 2000) are responsible for these PH3 detections. Observing cold metal-rich
objects may prove fruitful in our search for a clear phosphine detection, and
would be an good test our understand of phosphorus chemistry. Clearly the
chemistry of carbon dioxide and phosphine in these ultracool atmospheres are
exciting areas in need of further modeling efforts.

##  7 Conclusion

We show comparisons between three forward model grids and four NIRSpec PRISM
spectra of late T and Y dwarfs, two of which are being presented for the first
time in this paper. These comparisons reveal that the current generation of
forward models generally overestimate the PH3 absorption from 4$-$4.5$\mu$m,
and underestimate the CO2 absorption in this same range. As a further test, we
generate a grid of models using PICASO by taking the Elf Owl $P/T$ and
chemical profiles and modify the abundances of CO2 and PH3, and then fit our
four objects with this grid. These PICASO models show a preference for
enhanced CO2 abundance by at least a factor of a 100 and reduced the PH3
abundance for three the coldest objects compared to the Elf Owl models. The
warmest object shows minimal PH3 absorption in the base best-fit Elf Owl
model, and so does not need PH3 reduction. We present several ways one or both
of these chemical changes can be explained. These include a $K_{\mathrm{zz}}$
parameter that is variable with altitude, condensation of phosphorus in the
form of NH4H2PO4, and possibly an incomplete understanding of phosphorus and
carbon chemical pathways. However, our favored explanation for the required
CO2 enhancement is that the quench approximation does not accurately predict
the CO2 abundance, as CO2 remains in chemical equilibrium with CO after CO
quenches. This results in a $\sim$200 enhancement of CO2 for parameters
similar to the PICASO best fits. More modeling work is needed to understand
the phosphine and carbon dioxide chemistry in ultracool atmospheres.

##  8 acknowledgments

This work was inspired by conversations between some of the authors at The
First Year of JWST Science Conference, and we thank the organizers for the
this collaboration space. This work is based [in part] on observations made
with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from
the Mikulski Archive for Space Telescopes at the Space Telescope Science
Institute, which is operated by the Association of Universities for Research
in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These
observations are associated with program #2302. The specific observations
analyzed can be accessed via https://doi.org/10.17909/ntwg-k441 (catalog
10.17909/ntwg-k441). Support for program #2302 was provided by NASA through a
grant from the Space Telescope Science Institute, which is operated by the
Association of Universities for Research in Astronomy, Inc., under NASA
contract NAS 5-03127. This research has benefited from the Y Dwarf Compendium
maintained by Michael Cushing at
https://sites.google.com/view/ydwarfcompendium, and also from conversations
with the Brewster retrieval community. This research has made use of the
Spanish Virtual Observatory (https://svo.cab.inta-csic.es) project funded by
MCIN/AEI/10.13039/501100011033/ through grant PID2020-112949GB-I00. SM and CV
acknowledge support from JWST cycle 1 GO AR theory program PID-2232.

##  Appendix A Absorption Coefficient Plots

For completeness, we present the plots of absorption coefficients for
important molecules (H2O, CH4, CO, NH3, PH3, H2S, and CO2) from the best-fit
PICASO models for WISE 2102$-$44 and WISE 0359$-$54.

Figure 8: The same plot as Figure 5 but with the fits for WISE 2102$-$44\. The
abundances are from the PICASO best-fit model at P = 1 bar, where T = 861 K,
and the cross sections at T = 850 K. Figure 9: The same plot as Figure 5 but
with the fits for WISE 0359$-$54\. The abundances are from the PICASO best-fit
model at P = 1 bar, where T = 548 K, and the cross sections at T = 550 K.

## References

  * Barshay & Lewis (1978) Barshay, S. S., & Lewis, J. S. 1978, Icarus, 33, 593, doi: 10.1016/0019-1035(78)90192-6
  * Batalha et al. (2020) Batalha, N., Freedman, R., Gharib-Nezhad, E., & Lupu, R. 2020, Resampled Opacity Database for PICASO, Zenodo, doi: 10.5281/zenodo.6928501
  * Batalha et al. (2019) Batalha, N. E., Marley, M. S., Lewis, N. K., & Fortney, J. J. 2019, The Astrophysical Journal, 878, 70, doi: 10.3847/1538-4357/ab1b51
  * Beer (1975) Beer, R. 1975, The Astrophysical Journal, 200, L167, doi: 10.1086/181923
  * Beiler et al. (2023) Beiler, S. A., Cushing, M. C., Kirkpatrick, J. D., et al. 2023, The Astrophysical Journal Letters, 951, L48, doi: 10.3847/2041-8213/ace32c
  * Burgasser et al. (2023) Burgasser, A. J., Gerasimov, R., Bezanson, R., et al. 2023, UNCOVER: JWST Spectroscopy of Three Cold Brown Dwarfs at Kiloparsec-scale Distances, arXiv.  http://arxiv.org/abs/2308.12107
  * Böker et al. (2023) Böker, T., Beck, T. L., Birkmann, S. M., et al. 2023, Publications of the Astronomical Society of the Pacific, 135, 038001, doi: 10.1088/1538-3873/acb846
  * Cushing et al. (2006) Cushing, M. C., Roellig, T. L., Marley, M. S., et al. 2006, The Astrophysical Journal, 648, 614, doi: 10.1086/505637
  * Cushing et al. (2011) Cushing, M. C., Kirkpatrick, J. D., Gelino, C. R., et al. 2011, The Astrophysical Journal, 743, 50, doi: 10.1088/0004-637X/743/1/50
  * Fink et al. (1978) Fink, U., Larson, H. P., & Treffers, R. R. 1978, Icarus, 34, 344, doi: 10.1016/0019-1035(78)90172-0
  * Gillett & Forrest (1974) Gillett, F. C., & Forrest, W. J. 1974, The Astrophysical Journal, 187, L37, doi: 10.1086/181387
  * Gordon et al. (2022) Gordon, K. D., Bohlin, R., Sloan, G. C., et al. 2022, The Astronomical Journal, 163, 267, doi: 10.3847/1538-3881/ac66dc
  * Hubeny & Burrows (2007) Hubeny, I., & Burrows, A. 2007, The Astrophysical Journal, 669, 1248, doi: 10.1086/522107
  * Kirkpatrick et al. (2012) Kirkpatrick, J. D., Gelino, C. R., Cushing, M. C., et al. 2012, The Astrophysical Journal, 753, 156, doi: 10.1088/0004-637X/753/2/156
  * Kirkpatrick et al. (2019) Kirkpatrick, J. D., Martin, E. C., Smart, R. L., et al. 2019, The Astrophysical Journal Supplement Series, 240, 19, doi: 10.3847/1538-4365/aaf6af
  * Kothari et al. (2024) Kothari, H., Cushing, M. C., Burningham, B., et al. 2024, Probing the Heights and Depths of Y Dwarf Atmospheres: A Retrieval Analysis of the JWST Spectral Energy Distribution of WISE J035934.06$-$540154.6, doi: 10.48550/arXiv.2406.06493
  * Lacy & Burrows (2023) Lacy, B., & Burrows, A. 2023, The Astrophysical Journal, 950, 8, doi: 10.3847/1538-4357/acc8cb
  * Leggett & Tremblin (2023) Leggett, S. K., & Tremblin, P. 2023, The First Y Dwarf Data From JWST Show That Dynamic and Diabatic Processes Regulate Cold Brown Dwarf Atmospheres, arXiv, doi: 10.48550/arXiv.2309.14567
  * Lew et al. (2024) Lew, B. W. P., Roellig, T., Batalha, N. E., et al. 2024, The Astronomical Journal, 167, 237, doi: 10.3847/1538-3881/ad3425
  * Lodders & Fegley (2002) Lodders, K., & Fegley, B. 2002, Icarus, 155, 393, doi: 10.1006/icar.2001.6740
  * Luhman et al. (2023) Luhman, K. L., Tremblin, P., de Oliveira, C. A., et al. 2023, JWST/NIRSpec Observations of the Coldest Known Brown Dwarf, arXiv.  http://arxiv.org/abs/2311.17316
  * Mahaffy et al. (2000) Mahaffy, P. R., Niemann, H. B., Alpert, A., et al. 2000, Journal of Geophysical Research: Planets, 105, 15061, doi: 10.1029/1999JE001224
  * Marley et al. (1996) Marley, M. S., Saumon, D., Guillot, T., et al. 1996, Science, 272, 1919, doi: 10.1126/science.272.5270.1919
  * Marley et al. (2021) Marley, M. S., Saumon, D., Visscher, C., et al. 2021, The Astrophysical Journal, 920, 85, doi: 10.3847/1538-4357/ac141d
  * Meisner et al. (2023) Meisner, A. M., Leggett, S. K., Logsdon, S. E., et al. 2023, The Astronomical Journal, 166, 57, doi: 10.3847/1538-3881/acdb68
  * Meisner et al. (2021) Meisner, A. M., Schneider, A. C., Burgasser, A. J., et al. 2021, arXiv:2106.01387 [astro-ph].  http://arxiv.org/abs/2106.01387
  * Miles et al. (2020) Miles, B. E., Skemer, A. J. I., Morley, C. V., et al. 2020, The Astronomical Journal, 160, 63, doi: 10.3847/1538-3881/ab9114
  * Morley et al. (2018) Morley, C. V., Skemer, A. J., Allers, K. N., et al. 2018, The Astrophysical Journal, 858, 97, doi: 10.3847/1538-4357/aabe8b
  * Mukherjee et al. (2023a) Mukherjee, S., Batalha, N. E., Fortney, J. J., & Marley, M. S. 2023a, The Astrophysical Journal, 942, 71, doi: 10.3847/1538-4357/ac9f48
  * Mukherjee et al. (2022) Mukherjee, S., Fortney, J. J., Batalha, N. E., et al. 2022, The Astrophysical Journal, 938, 107, doi: 10.3847/1538-4357/ac8dfb
  * Mukherjee et al. (2023b) Mukherjee, S., Fortney, J., Morley, C., et al. 2023b, The Sonora Substellar Atmosphere Models. IV. Elf Owl: Atmospheric Mixing and Chemical Disequilibrium with Varying Metallicity and C/O Ratios (Y- type Models), Zenodo, doi: 10.5281/zenodo.10381250
  * Mukherjee et al. (2024) Mukherjee, S., Fortney, J. J., Morley, C. V., et al. 2024, The Sonora Substellar Atmosphere Models. IV. Elf Owl: Atmospheric Mixing and Chemical Disequilibrium with Varying Metallicity and C/O Ratios, doi: 10.48550/arXiv.2402.00756
  * Noll et al. (1997) Noll, K. S., Geballe, T. R., & Marley, M. S. 1997, The Astrophysical Journal, 489, L87, doi: 10.1086/310954
  * Oppenheimer et al. (1998) Oppenheimer, B. R., Kulkarni, S. R., Matthews, K., & Kerkwijk, M. H. v. 1998, The Astrophysical Journal, 502, 932, doi: 10.1086/305928
  * Orton et al. (2000) Orton, G. S., Serabyn, E., & Lee, Y. T. 2000, Icarus, 146, 48.  https://www.sciencedirect.com/science/article/pii/S0019103500963905
  * Phillips et al. (2020) Phillips, M. W., Tremblin, P., Baraffe, I., et al. 2020, doi: 10.1051/0004-6361/201937381
  * Prinn & Lewis (1975) Prinn, R. G., & Lewis, J. S. 1975, Science, 190, 274, doi: 10.1126/science.190.4211.274
  * Reach et al. (2005) Reach, W., Megeath, S., Cohen, M., et al. 2005, Publications of the Astronomical Society of the Pacific, 117, 978, doi: 10.1086/432670
  * Ridgway et al. (1976) Ridgway, S. T., Wallace, L., & Smith, G. R. 1976, The Astrophysical Journal, 207, 1002, doi: 10.1086/154570
  * Rigby et al. (2023) Rigby, J., Perrin, M., McElwain, M., et al. 2023, Publications of the Astronomical Society of the Pacific, 135, 048001, doi: 10.1088/1538-3873/acb293
  * Saumon et al. (2000) Saumon, D., Geballe, T. R., Leggett, S. K., et al. 2000, The Astrophysical Journal, 541, 374, doi: 10.1086/309410
  * Saumon & Marley (2008) Saumon, D., & Marley, M. S. 2008, The Astrophysical Journal, 689, 1327, doi: 10.1086/592734
  * Sorahana & Yamamura (2012) Sorahana, S., & Yamamura, I. 2012, The Astrophysical Journal, 760, 151, doi: 10.1088/0004-637X/760/2/151
  * Sorahana & Yamamura (2014) —. 2014, The Astrophysical Journal, 793, 47, doi: 10.1088/0004-637X/793/1/47
  * Tinney et al. (2014) Tinney, C. G., Faherty, J. K., Kirkpatrick, J. D., et al. 2014, The Astrophysical Journal, 796, 39, doi: 10.1088/0004-637X/796/1/39
  * Tsuji et al. (2011) Tsuji, T., Yamamura, I., & Sorahana, S. 2011, The Astrophysical Journal, 734, 73, doi: 10.1088/0004-637X/734/2/73
  * Visscher (2020) Visscher, C. 2020, Journal of Geophysical Research: Planets, 125, e2020JE006526, doi: 10.1029/2020JE006526
  * Visscher et al. (2006) Visscher, C., Lodders, K., & Fegley, Jr., B. 2006, The Astrophysical Journal, 648, 1181, doi: 10.1086/506245
  * Visscher & Moses (2011) Visscher, C., & Moses, J. I. 2011, The Astrophysical Journal, 738, 72, doi: 10.1088/0004-637X/738/1/72
  * Visscher et al. (2010) Visscher, C., Moses, J. I., & Saslow, S. A. 2010, Icarus, 209, 602, doi: 10.1016/j.icarus.2010.03.029
  * Wang et al. (2016) Wang, D., Lunine, J. I., & Mousis, O. 2016, Icarus, 276, 21, doi: 10.1016/j.icarus.2016.04.027
  * Wogan et al. (2024) Wogan, N. F., Batalha, N. E., Zahnle, K. J., et al. 2024, The Astrophysical Journal, 963, L7, doi: 10.3847/2041-8213/ad2616
  * Yamamura et al. (2010) Yamamura, I., Tsuji, T., & Tanabe, T. 2010, The Astrophysical Journal, 722, 682, doi: 10.1088/0004-637X/722/1/682
  * Zahnle & Marley (2014) Zahnle, K. J., & Marley, M. S. 2014, The Astrophysical Journal, 797, 41, doi: 10.1088/0004-637X/797/1/41


# A robust test of general relativity at the galactic scales by combining
strong lensing systems and gravitational wave standard sirens

Tonghua Liu  School of Physics and Optoelectronic Engineering, Yangtze
University, Jingzhou, 434023, China;  Marek Biesiada  National Centre for
Nuclear Research, Pasteura 7, PL-02-093 Warsaw, Poland;  Shuxun Tian
Department of Astronomy, Beijing Normal University, 100875, Beijing, China;
Kai Liao<EMAIL_ADDRESS>School of Physics and Technology, Wuhan
University, Wuhan 430072, China;

###### Abstract

The measurement of the parametrized post-Newtonian parameter
$\gamma_{\rm{PPN}}$ is a robust test of general relativity (GR). In some
modified theories of gravity, $\gamma_{\rm{PPN}}$ may evolve with the redshift
and deviate from one at high redshifts. This means that precise constraints on
$\gamma_{\rm{PPN}}$ acquired in the solar system experiments could not be
sufficient to test such theories and it is necessary to constrain
$\gamma_{\rm{PPN}}$ with high precision at high redshifts. However, in many
approaches aimed at extragalactic tests of GR, the results might be biased due
to entanglement of various factors, such as cosmic curvature, cosmic opacity,
and the Hubble constant. Strong lensing systems naturally provide a laboratory
to test $\gamma_{\rm{PPN}}$ at galactic scales and high redshifts, but there
is degeneracy between measured strength of gravity and cosmic distances in the
lensing system. Gravitational waves (GWs) from binary neutron star mergers
(standard sirens) provide a direct way to break this degeneracy by providing
self-calibrated measurements of the luminosity distance. We investigate the
possibility of estimating $\gamma_{\rm{PPN}}$ by combining well measured
strongly lensed systems with GW signals from coalescing neutron stars. Such
combination provides a cosmological-model independent, relatively pure and
unbiased method for the inference of $\gamma_{\rm{PPN}}$ parameter, avoiding
the influence of the above factors and the mass-sheet degeneracy in the lens.
Based on the simulated future 55 lensed quasar systems we demonstrated that
the precision of $\gamma_{\rm{PPN}}$ parameter obtained by our method could be
of order of $\sim 10^{-2}$. One may reasonably expect that our approach will
play an increasingly important role in precise testing the validity of general
relativity at galactic scales and high redshifts.

##  I Introduction

Over the past few decades, observations of Type Ia supernovae (SNe Ia) have
revealed that the expansion of the Universe is accelerating Riess et al.
(1998); Perlmutter et al. (1999). Based on Einstein’s theory of General
Relativity (GR) Wald (1984), and assuming that the Universe is homogeneous and
isotropic on large scales (Weinberg, 1972), it is generally accepted that the
so called $\Lambda$CDM model correctly describes our Universe. Although this
model is currently taken as the concordance cosmological scenario and
supported by the vast majority of astronomical observations (Planck
Collaboration et al., 2016), it still faces some problems. First, there is a
significant tension in the Hubble constant $H_{0}$ between the cosmic
microwave background (CMB) measurement obtained within GR$+\Lambda$CDM using
the _Planck_ data, which yielded $H_{0}=67.4\pm
0.5\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$ at the $68\%$ confidence level (CL)
(Planck Collaboration et al., 2020) and the value of $H_{0}=73.2\pm
1.3\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$ at the $68\%$ CL reported by SH0ES
(Supernova $H_{0}$ for the Equation of State) collaboration using SNe Ia
calibrated by local Cepheid variable stars (Riess et al., 2021). The $H_{0}$
Lenses in COSMOGRAIL’s Wellspring (H0LiCOW) collaboration recently reported
the measurement $H_{0}=73.3^{+1.7}_{-1.8}$
$\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$ with a joint analysis of six lensed
quasars for a flat $\Lambda$CDM model (Wong et al., 2020) using a technique
independent of above mentioned methods. Further, the same happens with
measurements of cosmic curvature (Di Valentino et al., 2020; Handley, 2021)
and the growth of structure $S_{8}$ parameter, which also reveal tension (see
a recent review (Abdalla et al., 2022) for more details) between alternative
approaches and analyses. These inconsistencies are the fairly challenging
problem in astrophysics and cosmology. From the theoretical point of view,
cosmological constant $\Lambda$ responsible for the accelerating expansion of
Universe should be invoked in the framework of GR, but inconsistency between
its observed value and theoretical predictions from the quantum field theory
is considerable (Weinberg, 1989).

Considering the fact that, the values of $H_{0}$ inferred from _Planck_ data
and H0LiCOW collaboration are based on GR plus $\Lambda$CDM model, it opens a
discussion of whether the GR could fail at larger, cosmological scales.
Although GR has passed with a very high precision all tests at the millimeter
scale in the laboratory, up to to the solar system scales (Ashby, 2002;
Bertotti et al., 2003), and the detection of gravitational waves also provided
the possibility for testing the validity of GR on very extreme scales
(Weisberg & Huang, 2016; Isi et al., 2019; Shao, 2021; Ding et al., 2021; Gair
et al., 2013), the long-range nature of gravity on the extragalactic or
cosmological scale is still relatively insufficiently tested. At present,
precise constraints regarding the parametrized post-Newtonian parameter
$\gamma_{\rm{PPN}}$ can be obtained at the Solar System or stellar scales in
the local Universe, e.g., $\Delta\gamma_{\rm{PPN}}\sim 10^{-5}$ given by
Cassini mission (Bertotti et al., 2003). However, the extragalactic
constraints on this parameter are much weaker, e.g.,
$\gamma_{\rm{PPN}}=0.97\pm 0.09$ on kiloparsec scale by using a nearby lens,
ESO 325-G004 given by (Collett et al., 2018). Recent reviews of the progress
in experimental testing of GR can be found in (Will, 2014; Ishak, 2019; Uzan,
2010). The issue of whether general relativity breaks down on larger
cosmological scales should be further validated (Koyama, 2016).

Independent alternative techniques could provide new perspectives. As one of
the most important predictions of GR, strong gravitational lensing, especially
strong lensing time-delay measurements, have become a powerful tool for
studying cosmology and gravity. In addition, lenses with measured time-delay
have recently been found to be more powerful cosmological probes with
capability to measure the angular diameter distances to the lenses $D^{A}_{d}$
by combining the measurements of time-delays, image configuration and stellar
kinematics (spectroscopy of lens galaxy). The possible deviations from GR
enter into the lens formula, which further affects the measured $D^{A}_{d}$.
However, the possible effects of GR deviation and measured distances are
highly degenerate. In other words, the gravitational mass is inferred from the
Einstein radius of the lens, but the observed Einstein radius also depends on
the distances in the optical system. Fortunately, gravitational wave (GW)
signals from inspiraling and merging compact binaries provide absolute
distance (luminosity distances $D^{L}$) (Schutz, 1986) and do not suffer from
interstellar extinction effects, which can help to break aforementioned
degeneracy. Such combination of time-delay lenses and GWs opens a new
possibility to test deviations from GR at galactic scales and high redshifts
in the limit of weak gravitational field. This letter is organized as follows:
In Sec. 2, we present the methodology and observational data. The results and
discussion are given in Sec. 3. Finally, we summarize our findings in Sec. 4.
The natural units of $c=G=1$ are adopted throughout this letter. The fiducial
cosmological model, a flat $\Lambda$CDM with $\Omega_{m}=0.30$ and $H_{0}=70$
$\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$ is assumed for simulating the lensing
data and GW data.

##  II Methodology

###  II.1 Background

In the weak field limit, background FLRW metric perturbed by the presence of
the lens can be written as

$ds^{2}=-\big{(}1+2\Phi\big{)}dt^{2}+a(t)^{2}\big{(}1-2\Psi\big{)}d\mathbf{r}^{2},$
(1)

where $a(t)$ is the scale factor, $\Phi$ is the Newtonian gravitational
potential, and $\Psi$ is the spatial curvature potential. We remark here that
$\Phi$ and $\Psi$ may be different in some modified gravity theories, i.e.,
$\Phi$ is gravitational potential that responds to the motion of non-
relativistic matter such as baryons and dark matter. The relationship of
$\Psi$ to the Newtonian potential is determined by the matter content (e.g.,
usual matter species in the concordance model plus possible extra scalar
fields) and gravity itself. The deviation from GR is quantified by the ratio
$\gamma_{\rm{PPN}}=\Psi/\Phi$ (denoted as the PPN parameter), with
$\gamma_{\rm{PPN}}=1$ in the framework of GR. In modified gravity theories,
$\Phi$ and $\Psi$ could be unequal with each other (Will, 2014), and even
their ratio could be a function of space (Hu & Sawicki, 2007; Saaidi et al.,
2011; Thomas et al., 2023) and time (Tian & Zhu, 2019). In this letter, we aim
to construct an unbiased method to constrain $\gamma_{\rm{PPN}}$ at galactic
scales and high redshifts. For simplicity, we assume that $\gamma_{\rm{PPN}}$
is independent of the space coordinates, and only depends on time. Such case
can be realized in many modified gravity theories, e.g., the nonlocal Gauss-
Bonnet gravity.

Strong lensing systems enables to perform two types of mass measurements:
gravitational mass of the lens inferred from the lensed images of the source,
and dynamical mass obtained from spectroscopic measurements of stellar
kinematics of the lensing galaxy. The dynamical mass is sensitive to the
Newtonian gravitational potential $\Psi$ only. On the other hand,
gravitational mass is sensitive to both potentials, i.e. to the Weyl potential
(defined as $\Phi_{+}=\frac{\Phi+\Psi}{2}=\frac{1+\gamma_{\rm{PPN}}}{2}\Phi$).
Thus, strong lensing provides a natural laboratory to test gravity and further
measure the PPN parameter $\gamma_{\rm{PPN}}$ by directly comparing the
difference between the dynamical mass and gravitational mass or comparing
$\Phi$ and $\Phi_{+}$.

###  II.2 Angular diameter distance from strong lensing

Combination of time delay and stellar kinematics measurements, can provide the
measurement of the angular diameter distance $D_{d}^{A}$ to the lens
(deflector). Let us briefly outline the standard procedure for determination
of $D_{d}^{A}$ pursued in the H0LiCOW program. For a typical strong lensing
system, with quasar at the redshift $z_{s}$ acting as a background source,
lensed by a foreground elliptical galaxy (at the redshift $z_{d}$), multiple
bright images of the active galactic nucleus (AGN) are formed together with
the arcs of its host galaxy. Time delays between multiple images can be
measured from monitoring of variability of the AGN light curves. From
theoretical point of view, lensing time delay is determined by both the
geometry of the Universe (different paths of rays forming different images) as
well as the Shapiro effects through (Shapiro, 1964)

$\Delta t_{\rm AB}=D_{\Delta t}\left[\phi(\theta_{\rm
A},\beta)-\phi(\theta_{\rm B},\beta)\right]=D_{\Delta t}\Delta\phi_{\rm
AB}(\xi_{\rm lens}),$ (2)

where $\phi(\theta,\beta)=\left[{(\theta-\beta)^{2}}/{2}-\psi(\theta)\right]$
is the Fermat potential at images, $\beta$ is the source position, $\psi$ is
effective lensing potential (the integral of the Weyl potential along the
line-of-sight) obeying the Poisson equation $\nabla^{2}\psi=2\kappa$, where
$\kappa$ (so called convergence) is the surface mass density of the lens in
units of critical density $\Sigma_{\rm crit}=D^{A}_{\rm s}/(4\pi D^{A}_{\rm
d}D^{A}_{\rm ds})$, and $\xi_{\rm lens}$ denotes the lens model parameters.
The cosmological background is reflected in the so-called ”time delay
distance” $D_{\Delta t}=(1+z_{\rm d})\frac{D^{A}_{\rm d}D^{A}_{\rm
s}}{D^{A}_{\rm ds}}$, where subscripts $d$ and $s$ stand for lens (deflector)
and source, and superscript $A$ denotes the angular diameter distance. The key
point here is that the Fermat potential difference $\Delta\phi_{\rm
AB}(\xi_{\rm lens})$ can be reconstructed by high-resolution lensing imaging
from space telescopes.

In the framework of parametrized post-Newtonian (PPN) formalism, the inferred
lensing mass parameters are rescaled by a factor of $(1+\gamma_{\rm{PPN}})/2$.
Hence, we denote the actually inferred lens model parameters in the Fermat
potential as $\xi^{{}^{\prime}}_{\rm lens}$. Adopting the notations in work
Yang et al. (2020), one can write the time-delay distance as

$D_{\Delta t}=(1+z_{\rm d})\frac{D^{A}_{\rm d}D^{A}_{\rm s}}{D^{A}_{\rm
ds}}=\frac{\Delta t_{\rm AB}}{\Delta\phi_{\rm AB}(\xi^{{}^{\prime}}_{\rm
lens})}\,.$ (3)

Thus, the time delay distance can be obtained from both the measurements of
time delay and the actual reconstructed Fermat potential with parameter
$\xi^{{}^{\prime}}_{\rm lens}$. The only difference in this equation is that
in the case of $\gamma_{\rm{PPN}}\neq 1$ under the PPN framework the inferred
lens model parameters are $\xi^{{}^{\prime}}_{\rm lens}$ but not the original
$\xi_{\rm lens}$ under GR.

On the other hand, assuming some explicit model of the lens such as the
simplest Singular Isothermal Sphere (SIS) model (or its extensions like
singular ellipsoid SIE or power-law model), observations regarding stellar
kinematics (dynamical mass determination) allow to obtain the following
distance ratio

$\frac{D^{\rm A}_{\rm s}}{D^{\rm A}_{\rm
ds}}=\frac{\sigma_{v}^{2}}{c^{2}J(\xi_{\rm{lens}},\xi_{\rm{light}},\beta_{\rm{ani}})}\;,$
(4)

where $\sigma_{v}$ is the line of sight (LOS) projected stellar velocity
dispersion of the lens galaxy. This distance ratio provides a valuable extra
constraint. The function $J$ captures all of the model components calculated
from the lensed images and luminosity-weighted projected velocity dispersion
(from the spectroscopy). Since the radial velocity dispersion $\sigma_{v}$ can
be modelled via the anisotropic Jeans equation, arguments of $J$ function
comprise lens model parameters $\xi_{\rm{lens}}$, parameters related to
luminosity distribution of the lens $\xi_{\rm{light}}$ and anisotropy
parameter $\beta_{\rm{ani}}$. More details concerning modelling issues related
to the function $J$ can be found in Section 4.6 of (Birrer et al., 2019).
Stellar kinematics of lensing galaxies is sensitive to the Newtonian potential
only, which is independent of PPN parameter. Thus, the inferred distance ratio
$D^{A}_{s}/D^{A}_{ds}$ combined with the well-measured velocity dispersion is
independent of the cosmological model and time delays, but still relies on the
lens model $\xi_{\rm{lens}}$ (Birrer et al., 2016, 2019). The lens model
parameter in $J$ is the ”unrescaled” $\xi_{\rm{lens}}$. If we replace
$\xi_{\rm{lens}}$ with $\xi^{{}^{\prime}}_{\rm{lens}}$, the resulted distance
ratio shall also be rescaled, correspondingly

$\frac{2}{1+\gamma_{\rm{PPN}}}\frac{D^{\rm A}_{\rm s}}{D^{\rm A}_{\rm
ds}}=\frac{\sigma_{v}^{2}}{c^{2}J(\xi^{\prime}_{\rm{lens}},\xi_{\rm{light}},\beta_{\rm{ani}})}\;.$
(5)

The PPN parameter is introduced here explicitly, and this formula was widely
used in the works on constraining the GR parameter for strong gravitational
lensing systems Cao et al. (2017); Liu et al. (2022); Wei et al. (2022); Liu &
Liao (2024). Furthermore, we can define $D^{\prime
A}_{d}=\frac{1+\gamma_{\rm{PPN}}}{2}D_{d}^{A}$, and by combining Eqs. (3) and
(5), the angular diameter distance to lens can be expressed as

$\displaystyle D^{\prime A}_{\rm d}=\frac{1}{1+z_{\rm d}}\frac{c\Delta t_{\rm
AB}}{\Delta\phi_{\rm
AB}(\xi^{\prime}_{\rm{lens}})}\frac{c^{2}J(\xi^{\prime}_{\rm{lens}},\xi_{\rm{light}},\beta_{\rm{ani}})}{\sigma_{v}^{2}}\,.$
(6)

Let us stress that this distance is unaffected by cosmic opacity111In the case
of gravitational lensing observables, cosmic opacity can change the absolute
intensity (magnitude) of images but not the relative intensity, thus not
biasing the distance determination in strong lensing system. Besides, the
velocity dispersion based on spectroscopic measurements are also unaffetced by
the opacity.. In fact, one could use the time delay distance (Eq. 3) to test
gravity theory, but this will inevitably introduce cosmic curvature parameter,
because $D^{A}_{\rm ds}$ is not directly observable. More importantly, one of
the main of obstacle for lensing mass modelling is the mass-sheet degeneracy,
which is completely circumvented by testing gravity using angular diameter
distances. Many factors such as the Hubble constant, dark energy model and
cosmic opacity (e.g. caused by light extinction by intergalactic dust) affect
determination of cosmological distances. This would bias most straightforward
tests of GR and reliable test of GR should be pure enough to ensure its
validity and being unbiased. Therefore, we point out in this work that GR
testing by gravitational lensing could be independent of these factors.

The current and future programs for time-delay with lensed quasars have great
progress, such as H0LiCOW Wong et al. (2020), COSMOGRAIL Eigenbrod et al.
(2005), STRIDES Treu et al. (2018), and SHARP, which have now combined in the
TDCOSMO collaboration Birrer & Treu (2021); Shajib et al. (2023). Although our
proposed method has many advantages, yet the existing sample of time-delay
lenses is still insufficient, so we turn to a new generation of wide and deep
sky observations and make appropriate simulations. Regarding the assessment of
the uncertainties of $D^{A}_{d}$, for each lens one should followed the
process presented in (Yıldırım et al., 2020), which resulted in 1.8% precision
of $D_{d}^{A}$ for RXJ1131-like system. In a realistic situation, the analysis
of each lens is very complicated, and such simulations are beyond this work.
Instead, we follow the work of (Jee et al., 2016) and consider that the main
sources of uncertainties give the expected uncertainty level for $D^{A}_{d}$
determination. These uncertainties, including the measurement of time delays,
the mass fluctuation along the LOS, recovering stellar kinematics from
measured velocity dispersions, the lens mass parameterisation from highly
resolved imaging will result in a few percent level of $D_{d}^{A}$
determination accuracy. We take the 5% uncertainty on the measurement of
$D_{d}^{A}$ as the best case scenario, which was also also used in (Jee et
al., 2016; Linder, 2011) and as the aim of the H0LiCOW program (Suyu et al.,
2017). In the conservative scenario, we also adopt 10% uncertainty on
$D_{d}^{A}$ for comparison. This conservative uncertainty is estimated from
the (Shajib et al., 2020), which used simulations to obtain realistic error
estimates with current/upcoming instruments, e.g. OH Suppressing InfraRed
Imaging Spectrograph on Keck (Larkin et al., 2006), and Near-Inrared
Spectrograph on James Webb Space Telescope (Birkmann et al., 2022), and
InfraRed Imaging Spectrograph on the Thirty Metre Telescope (Wright et al.,
2016).

Figure 1:  Redshift distribution of lenses and sources in a sample of 55
strong lensing systems expected to be observed in current and upcoming
projects.

In the near future, the next generation of wide and deep sky surveys, with
improved depth, area and resolution, may increase the current galactic-scale
lens sample sizes by orders of magnitude. For example, the upcoming Legacy
Survey of Space and Time (LSST) conducted in Vera Rubin Observatory, is
expected to find ten thousand lensed quasars (Ivezić et al., 2019). However,
considering that the metric efficiency is about 20%, the Time Delay Challenge
(TDC) program showed that only 400 well-measured time delays are available
(Liao et al., 2015). More precisely, LSST will discover 10000 lensed quasars,
but it has been estimated that only 2000 lensed quasar systems would have a
long enough (10$<\Delta t<$ 120 days) time delays measurements. The TDC
expects to be able to make time delay measurements with the high accuracy in
at least 20% efficiency of an LSST sample of 2000 lenses. This would
correspond to a well-measured sample of around 400 lensed quasars. In this
simulation, we use the OM10 catalog of mock
lenses222https://github.com/drphilmarshall/OM10 (Oguri & Marshall, 2010),
which has the distribution of time-delay lenses expected from the LSST.
Besides, one would need auxiliary data comprising high resolution imaging and
spectroscopy from instruments such as the Hubble Space Telescope (HST) and
ground-based observatories, in order to have all ingredients needed for
accurate lens mass modeling and determination of $D_{d}^{A}$. It should be
pointed out that LSST will discover a number of fainter, smaller-separation
lenses where it is not clear that the same level of precision of lens model
reconstruction can be reached. Therefore, we follow the work of (Jee et al.,
2016) and set the following criteria: 1) the lensed quasar image separation
should be greater than $1$ arcsec; 2) the third brightest quasar image has
i-band magnitude $m_{i}<21$ mag; 3) the lens galaxy has $m_{i}<22$ mag; (4)
the lensed quasar image is quadruple imaging lens systems (this type of system
provide more information to break the Source-Position Transformation
(Schneider & Sluse, 2013, 2014)). After applying these conditions to the OM10
catalog, 55 high-quality quadruple lens systems are kept (Jee et al., 2016).
Fig. 1 shows the source and the lens redshift distributions of 55 quadruple
lenses that match the selection criteria.

###  II.3 Luminosity distance from gravitational waves

From Eq. (6), it is clear that the PPN parameter and the angular diameter
distance are entangled. In order to break this degeneracy and complete the
test of the gravity theory, one needs independent, complementary distance
measurements. For this purpose, we turn to standard sirens.

Simultaneous detections of GW signal (Abbott et al., 2017a) from the binary
neutron star (NS-NS) merger and the electromagnetic (EM) counterparts (Abbott
et al., 2017b) from the same transient source opened a new era of multi-
messenger astronomy. GW signals from coalescing compact binary systems, i.e.
binary black holes (BH-BH) (Abbott et al., 2016), neutron stars (NS-NS) or
possible mixed neutron star - black hole (NS-BHs) systems provide us direct
measurements of luminosity distances $D^{L}$. Hence, they are called standard
sirens (Schutz, 1986). Standard sirens are self-calibrating, meaning that the
luminosity distances can be directly inferred from the detected waveforms
using matched-filter method. Redshift of the source $z$ is unfortunately non-
measurable in GW domain alone, and one has to identify the host galaxy in
order to obtain $z$. There are two advantages of standard sirens: firstly, GW
signals are unaffected by cosmic opacity, and propagate through the Universe
without any absorption and dissipation. Secondly, GW signals from standard
sirens provide direct luminosity distance measurements instead of the relative
distances as in case of SNe Ia. Thus $D^{L}$ from standard sirens do not need
to be anchored to the cosmic distance ladder.

Einstein Telescope (ET)333 The Einstein Telescope Project, https://www.et-
gw.eu/et/. (Dalal et al., 2006; Taylor, et al., 2012) is the third-generation
gravitational wave detector, which is designed to have a fantastic sensitivity
in the frequency range of $1-10^{4}$ Hz. See (Abernathy et al., 2011) for the
details of the ET conceptual design study. We simulate GW signals based on the
foreseen performance of the ET. Compared with the advanced detectors such as
AdLIGO and AdVirgo, such proposed third-generation detector aims for a
broadband factor of 10 sensitivity improvement, especially for the
characteristic distance sensitivity (Taylor & Gair, 2012). In this simulation,
we follow the work by (Califano et al., 2023a, b), and assume that the GW
sources are caused by binary merger of the NS with the NS, which can generate
the intense short gamma ray bursts (SGRBs) with measurable source redshifts,
or binary merger of the BH with the BH, which the redshift information comes
from a statistical analysis over a catalogue of potential host galaxies.

Following the strategy proposed in (Califano et al., 2023a, b), the
uncertainty for luminosity distance obtained from the GW signal from merger
event is

$\sigma^{2}_{\rm tot}=\sigma^{2}_{\rm inst}+\sigma^{2}_{\rm
len}+\sigma_{pec}^{2},$ (7)

where $\sigma_{\rm inst}$ is the instrumental measurement uncertainty and
calculated by

$\sigma_{\rm
inst}\simeq\sqrt{\bigg{\langle}\frac{\partial\mathcal{H}(f)}{\partial
D^{L,GW}},\frac{\partial\mathcal{H}(f)}{\partial
D^{L,GW}}\bigg{\rangle}^{-1}},$ (8)

where $\mathcal{H}(f)$ is the Fourier transform of strain $h(t)$ of a chirp
waveform in frequency domain. Since the strain $h$ and $\mathcal{H}$ is
inversely proportional to $D^{L,GW}$, we obtain $\sigma_{\rm inst}\simeq
D^{L,GW}/\rho$, where $\rho$ denotes the combined signal-noise-ratio (SNR) for
the network of independent interferometers, determined by the square root of
the inner product of $\mathcal{H}(f)$. Considering the uncertainty from the
inclination angle $\iota$ would also affect the SNR, the maximal inclination
angle effect from $\iota=0^{\circ}$ to $\iota=90^{\circ}$ on the SNR, we
choose to double the instrumental uncertainty of $D^{L}$ as the upper limit,
as was proposed in the recent study of the third-generation GW detector (ET)
based on Fisher information matrix (Cai & Yang, 2017), i.e.,

$\sigma_{\rm inst}\simeq\frac{2D^{L,GW}}{\rho}.$ (9)

The second uncertainty term $\sigma_{\rm len}$ is the weak gravitational
lensing effect caused by the large scale structure, and ignoring it would
result in the biased distance estimation. Following the work (Califano et al.,
2023a, b), this uncertainty is modeled as Hirata et al. (2010); Tamanini et
al. (2016)

$\sigma_{lens}=0.066\left(\frac{1-(1+z)^{-0.25}}{0.25}\right)^{1.8}D^{L}(z)F_{delens}(z),$
(10)

where $F_{delens}(z)=1-\frac{0.3}{\pi/2}\arctan{{z}/{z_{*}}}$, with
$z_{*}=0.073$ (Tamanini et al., 2016). The latter factor takes into account
the possibility to reduce the uncertainty due to weak lensing with the future
detectors such as the Extremely Large Telescope (Speri et al., 2021). The
final uncertainty term $\sigma_{pec}$ is caused by the peculiar velocity of
the host galaxy (Hjorth et al., 2017; Howlett & Davis, 2020; Mukherjee et al.,
2021), and can be approximated by a fitting formula (Kocsis et al., 2006)

$\sigma_{pec}=\left[1+\frac{c(1+z)^{2}}{H(z)D^{L}(z)}\right]\frac{\sqrt{\langle
v^{2}\rangle}}{c}D^{L}(z)\,,$ (11)

where the averaged peculiar velocity $\langle v^{2}\rangle$ is set as $500$
km/s, in agreement with the observed values in galaxy catalogs (Cen &
Ostriker, 2000).

Although the next generation ground-based detector ET is expected to detect
tens or hundreds of thousand NS-NS inspiral events per year up to the redshift
$z\sim 2$ and NS-BH mergers up to $z\sim 5$ with $\rho>8$, yet SGRBs necessary
for $z$ measurements are strongly beamed, and only the nearly face-on
configurations of NS-NS or NS-BH mergers are useful. Probability of their
occurence is $\sim 10^{-3}$. Assuming that the redshift distribution of
simulated GW sources follows the cosmic star formation history (Schneider et
al., 2001), we sample the the mass of neutron star and black hole within [1,
2] $M_{\odot}$ and [3, 10] $M_{\odot}$, respectively. Thus, we simulate 1000
GW events observable in the ET and their accurate redshifts up to $z_{GW}=5$.
Matching these events with the redshifts of the lens $z_{d}$ needed for
testing the PPN parameter, reduces the sample to 300 sources within redshift
$z_{GW}<1.25$.

##  III Results and discussion

Table 1:  Summary of the constraints on the $\gamma_{\rm{PPN}}$ parameter using the combination of strong lensing systems plus GW events. $\sigma_{D^{A}}/{D^{A}}=5\%$ | $\rm\gamma_{0}$ | $\rm\gamma_{1}$ | $\rm\gamma_{2}$
---|---|---|---
$\gamma_{\rm{PPN}}=\gamma_{0}$ | $1.003^{+0.018}_{-0.018}$ | $-$ | $-$
$\gamma_{\rm{PPN}}=1+\gamma_{1}*z$ | $-$ | $0.006^{+0.030}_{-0.031}$ | $-$
$\gamma_{\rm{PPN}}=1+\gamma_{2}*a$ | $-$ | $-$ | $0.008^{+0.028}_{-0.028}$
$\sigma_{D^{A}}/{D^{A}}=10\%$ | $\rm\gamma_{0}$ | $\rm\gamma_{1}$ | $\rm\gamma_{2}$
$\gamma_{\rm{PPN}}=\gamma_{0}$ | $1.003^{+0.032}_{-0.033}$ | $-$ | $-$
$\gamma_{\rm{PPN}}=1+\gamma_{1}*z$ | $-$ | $0.007^{+0.050}_{-0.050}$ | $-$
$\gamma_{\rm{PPN}}=1+\gamma_{2}*a$ | $-$ | $-$ | $0.004^{+0.050}_{-0.049}$

Figure 2:  Posterior distributions for PPN parameters obtained by the
combination of the future LSST survey and ET detector yields. Three panels
illustrate the following parametrizations: $\gamma_{\rm{PPN}}=\gamma_{0}$
constant case (left), $\gamma_{\rm{PPN}}=1+\gamma_{1}*z$ varying with redshift
(middle), and $\gamma_{\rm{PPN}}=1+\gamma_{1}*a$ varying with the scale factor
(right) .

To demonstrate the performance of our method to detect the possible deviation
of GR, we use three parameterized forms of $\gamma_{\rm{PPN}}$: 1)
$\gamma_{\rm{PPN}}=\gamma_{0}$ being constant; 2)
$\gamma_{\rm{PPN}}=1+\gamma_{1}*z$ changing with redshift; 3)
$\gamma_{\rm{PPN}}=1+\gamma_{2}*a$ changing with the scale factor $a=1/(1+z)$,
which unlike the redshift is the gravitational degree of freedom. The above
mentioned PPN parameters and their uncertainties are fitted by minimizing the
$\chi^{2}$ objective function defined in the following way:

$\chi^{2}=\sum^{55}_{i=1}\bigg{[}\frac{(D^{\prime
A}_{d,i}-D^{L}_{GW,i}(1+z)^{-2})^{2}}{\sigma_{D^{L}_{GW,i}}^{2}+\sigma_{D^{\prime
A}_{d,i}}^{2}}\bigg{]},$ (12)

where $\gamma_{\rm{PPN}}$ parameter is implicitly present in $D^{\prime
A}_{d}$. The total uncertainty consists of observational uncertainties from GW
and strong lensing contributions (see Section 2. B and C). The uncertainty
$\sigma_{D^{L}_{GW}}=\sigma_{\rm tot}/(1+z)^{2}$ is calculated using
$\sigma_{\rm tot}$ given by Eq. (7), strong lensing related uncertainties are
given by $\sigma_{D^{\prime A}_{d}}=5\%D^{\prime A}_{d}$ and
$\sigma_{D^{\prime A}_{d}}=10\%D^{\prime A}_{d}$ for the best case scenario
and the conservative scenario, respectively. We simulated ten thousand
realizations of the data with different random seeds and repeated the
minimization process to get an unbiased estimation for the PPN parameters. We
emphasize that in this work we do not attempt to constrain PPN parameter from
real data, but propose an unbiased approach and discuss its ability to explore
the presence of possible evolution or deviations from GR on simulated data
based on future observational forecasts.

The final Probability Distribution Functions (PDFs) for $\gamma_{\rm{PPN}}$
are reported in Table 1 and displayed in Fig. 2. Since the PDFs are
approximately Gaussian-like, we calculate the standard deviations as the
$1\sigma$ uncertainty levels. For the 5% uncertainty (optimistic case) on the
$D_{d}^{A}$ measurements from strong lensing systems, our results regarding
$\gamma_{\rm{PPN}}=\gamma_{0}$ constant case suggest that the uncertainty on
$\gamma_{0}$ is 0.018, corresponding to 1.8% precision. The precision obtained
by our method is $>10$ times better than the one obtained from current four
well-measured H0LiCOW quasars lenses, and almost 3 times better that the one
from simulated future 40 lensed quasar systems (Yang et al., 2020). Meanwhile,
the precision obtained by our method is $>4$ times better than expected from
the simulated future 10 strongly lensed fast radio bursts (FRBs) systems (Gao
et al., 2022). This result demonstrates the superiority of our method. There
are two reasons for this. Firstly, only one parameter $\gamma_{0}$ is
estimated in this work, whereas the previous work was simultaneously
estimating $H_{0}$ and PPN parameters, and the distance information was
derived there by assuming a concrete cosmological model. This affects the
ability to constrain PPN parameter alone. Secondly, one of the main of
obstacle for lensing mass modelling is the mass-sheet degeneracy444It is a
mathematical degeneracy that leaves the lensing observables unchanged, while
rescaling the absolute time delay, and thus the inferred $H_{0}$., which is
considered the dominant source of residual modeling error in time-delay
cosmography. Our method circumvents the mass-sheet degeneracy by using angular
diameter distances to perform the test of gravity, which greatly reduces
systemic uncertainty and further improves the precision on PPN parameter in
our work.

Some modified gravitational theories, such as the nonlocal Gauss-Bonnet
gravity (Tian & Zhu, 2019) require that PPN parameter may evolve with redshift
or the scale factor. With this motivation, we further assume a parameterized
form of $\gamma_{\rm{PPN}}$ changing with redshift or the scale factor. The
results are also shown in Fig. 2 and reported in Table 1. We obtain
$\Delta\gamma_{1}=0.030$ and $\Delta\gamma_{2}=0.028$ in the optimistic
scenario. Compared with PPN parameter treated as a constant, the precision of
PPN parameters evolving with redshift or the scale factor is noticeably
reduced. In Fig. 3, we plot the reconstruction of $\gamma_{\rm{PPN}}$ for the
different parameterizations with their uncertainties, and a prediction from a
nonlocal Gauss-Bonnet theory with specific parameters (Capozziello et al.,
2009). It is clearly seen that if GR does not hold and is replaced by this
nonlocal gravity, our method will distinguish them with very high precision.
Here, for the purple dashed line in Fig. 3, we adopt the conventions in (Tian
& Zhu, 2019) and set the model parameter $\alpha=5.6\times
10^{-6}c^{4}/H_{0}^{4}$, the integral constant $C_{1}=0$, and the matter
density $\rho\propto a^{-3}$. Note that, generally nonlocal Gauss-Bonnet
gravity predicts $|1-\gamma_{\rm{PPN}}|=\mathcal{O}(1)$ (Tian & Zhu, 2019).
The $\alpha$ we use here is too small to explain the cosmic late-time
acceleration, which is the initial motivation to propose this theory
(Capozziello et al., 2009). Therefore, this purple line corresponds to a toy
model. This illustration is just to show that $\gamma_{\rm{PPN}}$ in modified
gravity can be very close to 1 at low redshifts and deviate from 1 at high
redshifts. Such possible behavior of $\gamma_{\rm{PPN}}$ highlights the
importance of constraining it at high redshifts. In addition, the method
proposed in this paper is independent of any specific cosmological models and
gravitational theories.

The main purpose of this work is to quantify the ability of upcoming TDSL plus
GW to test GR against modified gravity theories. The key question now arises:
Is this combination of observations sufficient to detect possible effects of
deviation from GR? They can hardly become competitive to high-precision
measurements of $\gamma_{\rm{PPN}}=1+(2.1\pm 2.3)\times 10^{-5}$ in the solar
system made by Cassini mission. Considering the 10% uncertainty (conservative
case) on the $D_{d}^{A}$ measurements from strong lensing systems, the
precision on PPN parameter is almost half of that for the optimistic case.
This suggests that the main source of uncertainty in our method is the
observational error from the angular diameter distance measured by the
gravitational lensing system. This is reasonable, considering that only 55
gravitational lensing systems have been used in this work, and the uncertainty
of the PPN parameters caused by systematic errors is much greater than that
caused by statistical errors. If we take a more optimistic estimation
presented by (Yıldırım et al., 2020), the future uncertainty on $D_{d}^{A}$
can even decrease dramatically to the order of 1.8%, which indicates a more
promising future of TDSL plus GW technique for testing general relativity.
With the upcoming sky surveys such as LSST, the number of lensing systems will
increase dramatically. Hence, one can expect that our approach will greatly
improve the measurement precision on $\gamma_{\rm{PPN}}$ at the galactic
scale.

Figure 3: The reconstructed $\gamma_{\rm{PPN}}$ for the three
parameterizations with their $1\sigma$ uncertainties. The purple dashed line
represents the nonlocal Gauss-Bonnet theory. The green, pink, and blue bands
correspond to the $1\sigma$ regions of the best-fitting in the constant,
evolution with the redshift, and evolution with the scale factor cases,
respectively.

##  IV Conclusion

In this work, we investigated, for the first time, the possibility of
estimating post-Newtonian parameter $\gamma_{\rm{PPN}}$ by combining well
measured strongly lensed systems with gravitational wave signals. The
combination of strong lensing and self calibrating standard sirens observed in
GWs enables to avoid the possible bias that comes with assuming a particular
cosmological model. Meanwhile, such combination provides a relatively pure and
unbiased method for GR testing at the galactic scales and high redshifts, and
brings various benefits. It can reduce possible bias on GR testing induced by
cosmic curvature, cosmic opacity, the dark energy, the Hubble constant (the
$D^{A}_{d}$ provided by strong lensing and the $D^{L}$ provided by
gravitational waves are both absolute distances and are therefore unaffected
by the Hubble constant), and the mass-sheet degeneracy inside the lens.
Combining measurements of time delay, stellar velocity dispersion, high
resolution images, LOS environment modeling from upcoming LSST survey and GWs
acting standard siren from future ET detector, we obtained the precision of
$\sim$ 1.8% (optimistic case) and 3.2% (conservative case) in the case of
$\gamma_{\rm{PPN}}$ being a constant. Secondly, we considered the case of
$\gamma_{\rm{PPN}}$ parameter displaying evolution with redshift or the scale
factor and studied the precision with which this effect can be revealed. This
combination of data also allows us to infer post-Newtonian parameter and
detect possible deviation of GR at different redshifts. Although the precision
of testing GR at the galactic scale can hardly achieve the results competitive
to the measurements within the solar system, yet our technique would be able
to distinguish the departures from the GR quite precisely. We also indicated a
more promising future of TDSLs plus GWs for testing general relativity.

As a final remark, there are many potential ways to improve our method. For
instance, current and future surveys like the Dark Energy Survey (DES) (Treu
et al., 2018), the Hyper SuprimeCam Survey (More et al., 2017), and the Legacy
Survey of Space and Time (LSST) (Oguri & Marshall, 2010; Collett, 2015) will
bring us hundreds of thousands of lensed quasars in the most optimistic
discovery scenario. Even if only some small fraction of them will have precise
measurements of time delays between multiple images, the resulting statistics
will outshine current catalogs. With high-quality auxiliary observations, one
can use high-cadence, high-resolution and multi-filter imaging of the resolved
lensed images, to derive an accurate determination of the Fermat potential,
which will increase the precision of time delay distance by an order of
magnitude. On the other hand, the first strongly lensed supernova “SN Refsdal”
Kelly et al. (2015) with multiple images opened a new window for astrophysics
and cosmology with the concept of “lensed transients”. We also expect that the
discovery of various lensed transients explosive sources like gamma ray bursts
(GRBs), fast radio bursts (FRBs) and even GWs will give us new advantages over
the traditional targets in studying the Universe Liao et al. (2022, 2017). At
the same time, we also expect in the future, the synergies between GW and EMW
observations in various bands will yield reliable cosmological probe. It is
reasonable to expect that our approach will play an increasingly important
role in precise testing the validity of general relativity.

###### Acknowledgements.

This work was supported by National Natural Science Foundation of China under
Grant Nos. 12222302, 12203009, 11973034. Tonghua Liu was supported by Chutian
Scholars Program in Hubei Province (X2023007); Kai Liao was supported by Funds
for the Central Universities (Wuhan University 1302/600460081). Shuxun Tian
was supported by the Fundamental Research Funds for the Central Universities
(Beijing Normal University 10800-310400209522). Marek Biesiada was supported
by Hubei Province Foreign Expert Project (2023DJC040).

## References

  * Riess et al. (1998) Riess, A. G., Filippenko, A. V., Challis, P., et al. 1998, Astron. J., 116, 1009.
  * Perlmutter et al. (1999) Perlmutter, S., Aldering, G., Goldhaber, G., et al. 1999, Astrophys. J. , 517, 565.
  * Wald (1984) Wald R. 1984 ”General Relativity.”
  * Weinberg (1972) Weinberg, S. 1972, Gravitation and Cosmology: Principles and Applications of the General Theory of Relativity, by Steven Weinberg, pp. 688. ISBN 0-471-92567-5. Wiley-VCH , July 1972., 688
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, Astron. Astrophys., 594, A13.
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, Astron. Astrophys., 641, A6.
  * Riess et al. (2021) Riess, A.G., Casertano, S., Yuan, W., Bowers, J.B., Macri, L., Zinn, J.C. and Scolnic, D. 2021, Astrophys. J. 908, L6.
  * Wong et al. (2020) Wong, K. C., Suyu, S. H., Chen, G. C.-F., et al. 2020, Mon. Not. Roy. Astron. Soc., 498, 1420.
  * Di Valentino et al. (2020) Di Valentino, E., Melchiorri, A., & Silk, J. 2020, Nature Astronomy, 4, 196.
  * Handley (2021) Handley, W. 2021, Phys. Rev. D, 103, L041301.
  * Abdalla et al. (2022) Abdalla, E., Abellán, G. F., Aboubrahim, A., et al. 2022, Journal of High Energy Astrophysics, 34, 49.
  * Weinberg (1989) Weinberg, S. 1989, Rev. Modern Phys., 61, 1
  * Ashby (2002) Ashby, N. 2002, Phys. Today, 55, 41
  * Bertotti et al. (2003) Bertotti, B., Iess, L., & Tortora, P. 2003, Nature, 425, 374
  * Weisberg & Huang (2016) Weisberg, J. M. & Huang, Y. 2016, Astrophys. J. , 829, 55.
  * Shao (2021) Shao, L. 2021, Physics Online Journal, 14, 173.
  * Isi et al. (2019) Isi, M., Chatziioannou, K., & Farr, W. M. 2019, Phys. Rev. Lett. , 123, 121101.
  * Ding et al. (2021) Ding, H., Deller, A. T., Fonseca, E., et al. 2021, Astrophys. J. Lett., 921, L19.
  * Gair et al. (2013) Gair, J. R., Vallisneri, M., Larson, S. L., et al. 2013, Living Reviews in Relativity, 16, 7.
  * Collett et al. (2018) Collett, T. E., Oldham, L. J., Smith, R. J., et al. 2018, Science, 360, 1342.
  * Will (2014) Will, C. M. 2014, Living Reviews in Relativity, 17, 4
  * Ishak (2019) Ishak, M. 2019, Living Reviews in Relativity, 22, 1.
  * Uzan (2010) Uzan, J.-P. 2010, General Relativity and Gravitation, 42, 2219.
  * Koyama (2016) Koyama, K. 2016, Rept. Prog. Phys., 79, 046902
  * Schutz (1986) Schutz, B. F., 1986, Nature, 323, 310
  * Hu & Sawicki (2007) Hu, W. & Sawicki, I. 2007, Phys. Rev. D, 76, 064004.
  * Saaidi et al. (2011) Saaidi, K., Mohammadi, A., & Sheikhahmadi, H. 2011, Phys. Rev. D, 83, 104019.
  * Thomas et al. (2023) Thomas, D. B., Clifton, T., & Anton, T. 2023, JCAP, 2023, 016.
  * Tian & Zhu (2019) Tian, S. X. & Zhu, Z.-H. 2019, Phys. Rev. D, 99, 064044.
  * Shapiro (1964) Shapiro, I. I. 1964, Phys. Rev. Lett. , 13, 789.
  * Yang et al. (2020) Yang, T., Birrer, S., & Hu, B. 2020, Mon. Not. Roy. Astron. Soc., 497, L56.
  * Birrer et al. (2019) Birrer, S., Treu, T., Rusu, C. E., et al. 2019, Mon. Not. Roy. Astron. Soc., 484, 4726.
  * Birrer et al. (2016) Birrer, S., Amara, A., & Refregier, A. 2016, JCAP, 2016, 020.
  * Cao et al. (2017) Cao, S., Li, X., Biesiada, M., et al. 2017, Astrophys. J. , 835, 92.
  * Liu et al. (2022) Liu, X.-H., Li, Z.-H., Qi, J.-Z., et al. 2022, Astrophys. J. , 927, 28.
  * Wei et al. (2022) Wei, J.-J., Chen, Y., Cao, S., et al. 2022, Astrophys. J. Lett., 927, L1.
  * Liu & Liao (2024) Liu, T. & Liao, K. 2024, Mon. Not. Roy. Astron. Soc., 528, 1354. doi:10.1093/mnras/stae119
  * Eigenbrod et al. (2005) Eigenbrod, A., Courbin, F., Vuissoz, C., et al. 2005, Astron. Astrophys., 436, 25.
  * Treu et al. (2018) Treu, T., Agnello, A., Baumer, M. A., et al. 2018, Mon. Not. Roy. Astron. Soc., 481, 1041.
  * Birrer & Treu (2021) Birrer, S. & Treu, T. 2021, Astron. Astrophys., 649, A61.
  * Shajib et al. (2023) Shajib, A. J., Mozumdar, P., Chen, G. C.-F., et al. 2023, Astron. Astrophys., 673, A9.
  * Yıldırım et al. (2020) Yıldırım, A., Suyu, S. H., & Halkola, A. 2020, Mon. Not. Roy. Astron. Soc., 493, 4783.
  * Jee et al. (2016) Jee, I., Komatsu, E., Suyu, S. H., et al. 2016, JCAP, 04, 031.
  * Linder (2011) Linder, E. V. 2011, Phys. Rev. D, 84, 123529.
  * Suyu et al. (2017) Suyu, S. H., Bonvin, V., Courbin, F., et al. 2017, Mon. Not. Roy. Astron. Soc., 468, 2590.
  * Shajib et al. (2020) Shajib, A. J., Birrer, S., Treu, T., et al. 2020, Mon. Not. Roy. Astron. Soc., 494, 6072.
  * Larkin et al. (2006) Larkin, J., Barczys, M., Krabbe, A., et al. 2006, New Astronomy Reviews, 50, 362.
  * Birkmann et al. (2022) Birkmann, S. M., Ferruit, P., Giardino, G., et al. 2022, Astron. Astrophys., 661, A83.
  * Wright et al. (2016) Wright, S. A., Walth, G., Do, T., et al. 2016, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 9909, 990905.
  * Ivezić et al. (2019) Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, Astrophys. J. , 873, 111.
  * Liao et al. (2015) Liao, K., Treu, T., Marshall, P., et al. 2015, Astrophys. J. , 800, 11.
  * Oguri & Marshall (2010) Oguri, M. & Marshall, P. J. 2010, Mon. Not. Roy. Astron. Soc., 405, 2579.
  * Schneider & Sluse (2013) Schneider, P. & Sluse, D. 2013, Astron. Astrophys., 559, A37.
  * Schneider & Sluse (2014) Schneider, P. & Sluse, D. 2014, Astron. Astrophys., 564, A103.
  * Abbott et al. (2017a) Abbott, B. P., et al. [LIGO Scientific Collaboration, the Virgo Collaboration], 2017a, Phys. Rev. Lett. , 119, 161101
  * Abbott et al. (2017b) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017b, Astrophys. J. Lett., 848, L12.
  * Abbott et al. (2016) Abbott, B. P., et al. [LIGO Scientific and Virgo Collaborations], 2016, Phys. Rev. Lett. , 116, 061102
  * Dalal et al. (2006) Dalal, N., Holz, D. E., Hughes, S.A., Jain, B. 2006, PRD, 74, 063006
  * Taylor, et al. (2012) Taylor, S. R., et al. 2012, Phys. Rev. D, 85, 023535
  * Abernathy et al. (2011) Abernathy, M., et al., “Einstein gravitational wave Telescope: Conceptual Design Study”, 2011 [available from European Gravitational Observatory, document number ET-0106A-10].
  * Taylor & Gair (2012) Taylor, S. R. & Gair, J. R. 2012, Phys. Rev. D, 86, 023502
  * Califano et al. (2023a) Califano, M., de Martino, I., Vernieri, D., et al. 2023a, Mon. Not. Roy. Astron. Soc., 518, 3372.
  * Califano et al. (2023b) Califano, M., de Martino, I., Vernieri, D., et al. 2023b, Phys. Rev. D, 107, 123519.
  * Cai & Yang (2017) Cai, R.-G. & Yang, T. 2017, Phys. Rev. D, 95, 044024
  * Hirata et al. (2010) Hirata, C. M., Holz, D. E., & Cutler, C. 2010, Phys. Rev. D, 81, 124046.
  * Tamanini et al. (2016) Tamanini, N., Caprini, C., Barausse, E., et al. 2016, JCAP, 2016, 002.
  * Speri et al. (2021) Speri, L., Tamanini, N., Caldwell, R. R., et al. 2021, Phys. Rev. D, 103, 083526.
  * Hjorth et al. (2017) Hjorth, J., Levan, A. J., Tanvir, N. R., et al. 2017, Astrophys. J. Lett., 848, L31.
  * Howlett & Davis (2020) Howlett, C. & Davis, T. M. 2020, Mon. Not. Roy. Astron. Soc., 492, 3803.
  * Mukherjee et al. (2021) Mukherjee, S., Lavaux, G., Bouchet, F. R., et al. 2021, Astron. Astrophys., 646, A65. 2001, Mon. Not. Roy. Astron. Soc., 324, 797.
  * Kocsis et al. (2006) Kocsis, B., Frei, Z., Haiman, Z., et al. 2006, Astrophys. J. , 637, 27.
  * Cen & Ostriker (2000) Cen, R. & Ostriker, J. P. 2000, Astrophys. J. , 538, 83.
  * Schneider et al. (2001) Schneider, R., Ferrari, V., Matarrese, S., et al. 2001, Mon. Not. Roy. Astron. Soc., 324, 797.
  * Gao et al. (2022) Gao, R., Li, Z., & Gao, H. 2022, Mon. Not. Roy. Astron. Soc., 516, 1977.
  * Capozziello et al. (2009) Capozziello, S., Elizalde, E., Nojiri, S., et al. 2009, Physics Letters B, 671, 193.
  * More et al. (2017) More, A., Lee, C.-H., Oguri, M., et al. 2017, Mon. Not. Roy. Astron. Soc., 465, 2411.
  * Collett (2015) Collett, T. E. 2015, Astrophys. J. , 811, 20.
  * Liao et al. (2022) Liao, K., Biesiada, M., & Zhu, Z.-H. 2022, Chinese Physics Letters, 39, 119801.
  * Liao et al. (2017) Liao, K., Fan, X.-L., Ding, X., et al. 2017, Nature Communications, 8, 1148.
  * Kelly et al. (2015) Kelly, P. L., Rodney, S. A., Treu, T., et al. 2015, Science, 347, 1123.