# Are Parametrized Tests of General Relativity with Gravitational Waves Robust
to Unknown Higher Post-Newtonian Order Effects?

Scott Perkins<EMAIL_ADDRESS>Illinois Center for Advanced Study of the
Universe & Department of Physics, University of Illinois at Urbana-Champaign,
Champaign, Illinois 61820 USA  Center for AstroPhysical Surveys, National
Center for Supercomputing Applications, Urbana, Illinois 61801 USA  Nicolás
Yunes<EMAIL_ADDRESS>Illinois Center for Advanced Study of the Universe
& Department of Physics, University of Illinois at Urbana-Champaign,
Champaign, Illinois 61820 USA

###### Abstract

Gravitational wave observations have great potential to reveal new information
about the fundamental nature of gravity, but extracting that information can
be difficult. One popular technique is the parametrized inspiral test of
general relativity (a realization of the parametrized post-Einsteinian
framework), where the gravitational waveform, as calculated in Einstein’s
theory as a series expansion in the orbital velocity, is parametrically
deformed at a given set of orders in velocity. However, most current
approaches usually only analyze the data while considering a single, specific
modification at a time. Are then constraints placed with a single modification
robust to our ignorance of higher post-Newtonian order corrections? We show
here that for a wide class of theories, specifically those that admit a post-
Newtonian expansion, single-parameter tests are indeed robust. In particular,
through a series of full Bayesian parameter estimation studies on several
different sets of synthetic data, we show that single-parameter constraints
are not degraded but rather are improved by the inclusion of multiple
parameters, provided one includes information about the mathematical structure
of the series. We then exemplify this with a specific theory of gravity,
shift-symmetric scalar Gauss-Bonnet theory, where the waveform has been
calculated to higher post-Newtonian orders than leading. We show that the
inclusion of these higher order terms strengthens single-parameter
constraints, instead of weakening them, and that the strengthening is very
mild. This analysis therefore provides strong evidence that single-parameter
post-Einsteinian tests of general relativity are robust to ignorance of high
post-Newtonian order terms in the general relativistic deformations.

##  I Introduction

As the field of gravitational wave (GW) astrophysics matures, the growing
amount of data is being used in a variety of ways and in different contexts to
solve problems in topics like fundamental physics, astrophysics, and
cosmology. One critical subject that has drastically benefited by the
observation of GWs is the study of the gravitational interaction. General
relativity (GR), our best theory of gravity to date, has some shortcomings, if
we may call them that [1, 2]. These include the incompatibility of GR with
quantum mechanics (c.f. Ref. [3]) and the singularity problem [4, 5], as
examples. Aside from these shortcomings, the theory also struggles to
adequately explain certain phenomena observed in our universe without the
inclusion of additional structure, like the late time acceleration of the
universe (unless one includes an “unnaturally” small cosmological constant)
[6, 7], the rotation curves of galaxies (unless one includes dark matter
particles that have so far been undetected via direct experiment) [8, 9], and
the matter-antimatter asymmetry of the universe (unless one includes new
parity-violating interactions that satisfy the Sakharov conditions) [10].
Modifications to GR in the ultraviolet or in the infrared may potentially
resolve some of these theoretical and observational issues, which is why they
have attracted attention lately in light of new GW data.

If we are searching for or attempting to constrain “alternatives” to GR, GW
observations are an excellent place to look [11, 12]. The dynamics of the
system that generates them are highly relativistic (reaching characteristic
speeds up to a substantial fraction of the speed of light), and involve large
gravitational potentials and fields. However, deducing information about
fundamental physics from GW data is a nuanced business, where the information
encoded inside the waveform can, in some cases, be very faint and buried in
noise. One possible philosophy is to take a maximally “ignorant” approach,
where one seeks to involve as little outside information as possible. This
approach has its appeal in terms of its perceived robustness (uninfluenced by
any external, _a priori_ information), but the exact method with which one
encodes this “ignorance” is not a trivial or well defined matter. Furthermore,
a maximally uninformed approach can minimize and sometimes even erase the
strength of any inferences about fundamental physics one might hope to extract
[13]. While it would be appealing to conduct these analyses independently of
our priors, this seems to be a quixotic approach.

For example, one maximally uninformed approach would be to reconstruct the
signal with an orthogonal basis of wavelets, as is normally done with
BayesWave [14, 15]. The fact that this method is totally agnostic to what
physics produced the signal allows one to reconstruct the signal to almost
arbitrary precision, which is useful in many contexts. However, such an
approach usually provides little information about the theory of gravity that
led to the signal that is being reconstructed, or more specifically about the
fundamental physics involved in the generation and propagation of GW in our
universe (although there have been recent attempts to tie in to certain topics
in fundamental physics [16]). Given this, it’s useful to remember that general
tests like “residual tests,” in which one studies the signal-to-noise ratio
contained in the difference between the data and the best re-constructed
signal [17, 18], are inherently a consistency check of the reconstruction
procedure, and should not be interpreted as direct tests of the underlying
physics. The value of such tests is only as great as the accuracy of the
underlying reconstruction models, regardless of what physics those models are
trying to represent.

Meanwhile, a less uninformative approach, though still agnostic, might attempt
to model the signal as a GW produced in GR, but with some small deformation
that is represented parametrically. This type of approach defines the
parametrized post-Einsteinian (ppE) formalism [19, 17, 20, 21], which is the
general framework that the “parametrized inspiral tests of GR,” used by the
LIGO/Virgo collaboration (LVC), is built from [22, 23, 24, 25, 26, 27, 28, 29,
30, 31, 13]. In this approach, deformations from GR are incorporated directly
at the level of the model for the signal (i.e. the GW template or GW model)
that one will use to compare against the data. Moreover, these deformations
are represented parametrically through a particular basis (i.e. in the
inspiral, a polynomial in orbital velocity), with the polynomial exponent of
the leading post-Newtonian (PN) order term encoding the type of modification
one is considering, and the polynomial coefficient the strength of the
deformation.

The resulting functional form of the deformation of the signal is not
arbitrary. In fact, one can easily show that such a representation derives
from parametrically deforming the binding energy of the binary, the rate of
change of this binding energy, or the dispersion relation of the propagating
GW through a polynomial in velocity [19]. Such an analysis also reveals that
there is a mapping between deformations in these physical quantities and those
that are introduced in the waveform, thus allowing for constraints on the
waveform deformations to then be mapped back to constraints on theoretical
physics [19, 12, 32, 33]. Moreover, a polynomial in velocity is a natural
basis for expansions in the inspiral phase, where the PN expansion is expected
to hold to approximate the solution to the field equations. Obviously, such a
representation need not apply in the merger, and would have to be changed in
the ringdown, as discussed extensively in [19], but we are here focusing on
inspiral tests only.

The initial and simplest ppE proposal was to include a single deformation at a
time [19]. As explained above, the deformation is represented by adding a term
of the form $\beta\,v^{b}$ in the GW Fourier phase (and a similar term in the
Fourier amplitude), where $\beta$ is a real parameter one is attempting to
estimate, $v$ is the orbital velocity of the binary (which is connected to the
GW frequency via Kepler’s third law), and $b$ is a fixed constant. One
attempts to measure or constrain $\beta$ (sometimes called the ppE amplitude
coefficient) because it determines the strength of the GR deformation, and it
is connected to the coupling constant of modified GR theories. The quantity
$b$ (sometimes called the ppE exponent coefficient) is not a parameter, but
rather it is a real number that determines the type of modification one is
considering; therefore, one chooses a value of $b$ a priori, before carrying
out parameter estimation or model selection. By picking a set of such $b$
constants, one can then derive posterior probability distributions for each
$\beta$ (corresponding to each value of $b$ that one chose), to then construct
the so-called “violin plots” the LVC generates in their testing GR papers [13,
31, 24, 25].

However, as explained when the ppE framework was first developed [19],
deformations to GR in the inspiral should not consist of a single PN
modification. Rather, one expects modifications to enter at some leading PN
order (a so-called “controlling factor” of an asymptotic expansion),
multiplied by an entire PN series in velocity, i.e.
$\beta_{0}\,v^{b}(1+\sum_{i=n}\beta_{n}v^{n})$. One way to incorporate this
idea into an inspiral test in GR would be to truncate the series at some PN
order, allow every $\beta_{n}$ coefficient to be independent, and attempt to
estimate all of them simultaneously [23, 34]. As argued in [20], and then
verified in [13], this approach is doomed to fail because of large covariances
between the $\beta_{n}$ parameters that prevents one from estimating any one
of them with any precision (i.e. the marginalized posteriors of $\beta_{n}$
that one recovers are very similar to the priors one chooses for them). If
allowing for multiple, independent deformation terms at once erases all useful
information one may glean from the GW signal, are constraints one obtains with
a single deformation at a time robust?

At first sight, one may think the answer to this question is no. After all, we
know that a modified gravity signal will not contain just a single
modification in the phase. Therefore, forcing a template to be of this form
ought to introduce uncontrolled systematic error. Moreover, the velocities at
merger are close to $0.3$–$0.5$ the speed of light, so PN terms at higher-
than-leading order need not be negligible. Therefore, by neglecting these
higher order terms one could be ignoring “degeneracies” that could deteriorate
constraints if they were properly modeled. If this were the case, single-
parameter ppE tests or single-parameter parametrized inspiral tests could
over-confidently predict and place constraints that are more stringent than
what we have any right in claiming.

We are of course not the first to study the concerns mentioned above. To our
knowledge, one of the first studies to do this for aLIGO observations appeared
in an Appendix [12], where the authors investigated how constraints on Brans-
Dicke theory are improved as higher PN order terms are added to the inspiral
phase in the extreme mass-ratio limit. The authors found that constraints on
the Brans-Dicke coupling parameter are very stable to the inclusion of these
higher PN order terms, with relative fractional changes of at most a few tens
of percent. A more recent study of this same topic was by [35], where
constraints on a specific theory of gravity, shift-symmetric scalar Gauss-
Bonnet (ssGB), were studied with the O2 and O3 GW catalog. In this work, the
authors modeled the higher PN order terms as unknown, and then marginalized
over them; they found that the inclusion of these higher order terms only
changed the constraints again by only tens of percent.

In spite of these studies, there has not yet been a dedicated study of the
concerns mentioned above, which is what motivated this paper. As we will show
in detail here, and in agreement with the previous work described above, the
pessimistic viewpoint expressed above regarding single- versus multiple-
parameter tests of GR is completely unwarranted. For a multitude of reasons,
which we will elaborate on in this work, constraints obtained by using single-
parameter models are both robust and reliable. In particular, we will show
that although systematic error is incurred when one uses a single-parameter
ppE model, the error is small and tends to predict a less stringent constraint
than what one would obtain if one included higher PN order corrections subject
to a reasonable PN prior; in this sense, single-parameter ppE and parametrized
inspiral tests are conservative. Moreover, we will show that the improvement
one gains by including higher PN order terms in the GR deformation makes
constraints only slightly more stringent (of at most one order of magnitude in
the most extreme cases).

We arrive at these conclusions through a detailed and full Bayesian analysis
of a multitude of synthetic signals. We first carry out a fully Bayesian
Markov-chain Monte-Carlo (MCMC) parameter estimation study, in which we inject
a synthetic GR signal in stationary and Gaussian noise, and then we extract it
with single-parameter ppE models, as well as multi-parameter models. For the
latter, we develop a new prior based on the theoretical foundations of PN
theory [36, 37], which ensures that the modifications one introduces lead to a
non-GR PN series that has similar or better convergence properties than that
of GR. If the modified theory of gravity one is considering accepts a PN
expansion, then the prior we develop is guaranteed to be valid. It is the
development and use of this prior that restricts the impact of the covariances
found in previous studies [13, 23, 34].

We then examine our conclusions under the light of a specific example:
constraints on ssGB theory. Higher-than-leading PN order waveforms have been
recently calculated in this theory [38], thus allowing us to verify our prior.
With this at hand, we carry out again a fully Bayesian test of GR, injecting
synthetic GR signals in stationary and Gaussian noise, and extracting them
with the new, high-PN-order ssGB model, the leading-PN-order ssGB model, a
single-parameter ppE model, and a multiple-parameter ppE model with the PN-
inspired prior developed above. In all cases, we map constraints on the
parameters we search over to constraints on the coupling constant of ssGB
gravity. We find that the the leading PN order ssGB model, the single-
parameter ppE model, and the multiple-parameter ppE model with the PN-inspired
prior all lead to roughly the same constraints, in agreement with previous
work [35]. The new, high-PN-order ssGB model leads to constraints that are
stronger (not weaker) than the other constraint, but only by a factor of
roughly 3 (at 90% confidence), because this model adds new physical
information (at higher PN orders). Therefore, our analysis proves that, at
least in this theory (and very likely in all theories of this type), higher-PN
order corrections to the ssGB waveform improve constraints instead of
deteriorating them, and the improvement is only very mild.

The conclusions we arrive at in this paper are admittedly strong, so they come
with a couple of caveats. One of them is that the signal one analyzes is
dominated by the inspiral (and not the merger) phase of the coalescence. If
one can only hear the last few cycles of coalescence and the ringdown, then
one is outside the regime of validity of the PN approximation, rendering some
of the arguments presented above invalid. In practice, this implies that our
conclusions apply only to binaries of sufficiently small total mass, with the
maximum mass allowed dependent on the characteristics of the instrument (and
in particular, the low frequency seismic wall of the noise spectral density of
the detector); for the Advanced LIGO detectors at design sensitivity, the
validity of our conclusions require that at least the total mass of the binary
be less than roughly $40M_{\odot}$, so that at least the frequency of the
innermost stable circular orbit (of the effective one-body problem) be above
in the sensitivity bucket of aLIGO at 100 Hz. Another caveat is that we
consider modifications to GR that lead to “persistent” effects in the
inspiral, as opposed to modifications that turn on suddenly during the
inspiral. The latter can occur in theories with additional length scales that
lead to screening (like massive gravity theories [39]), or theories with
additional fields that activate during the inspiral due to non-linear effects
(like dynamical scalarization [40, 41, 42, 43] or vectorization [44, 45]). As
shown in [46], however, even for such theories a single-parameter ppE test is
sufficient to detect such sudden deviations (at the cost of deteriorating the
effectiveness of the test).

The rest of this paper will present the details of the calculations that have
led us to the conclusions we described above, and it is structured as follows.
We begin by discussing some of basics of Bayesian inference with GWs in Sec.
II.1. This is followed by a discussion of the current methodology for testing
GR in the inspiral phase of GW binaries in Sec. II.2. We continue with a
discussion of a proposed improvement on those methods in Sec. III, in which we
outline a framework informed by the PN formalism. After this, we expand on the
experimental design we use to test this new framework in Sec. IV. With the
experiment summarized, we move on to discuss the results and implications of
those experiments in Sec. V. In Sec. VI, we investigate alternative
parametrizations to determine the robustness of our conclusions. From here, we
focus on a specific theory of gravity, ssGB, to determine how realistic our
conclusions are in the context of actual modifications to GR in Sec. VII.
Finally, we summarize our findings and propose future research in Sec. VIII.
Throughout this work, we will work in geometric units, where $G=1=c$.

##  II Testing GR with GWs

Inferences about fundamental physics from GW data begins with matched
filtering and Bayesian signal analysis, which we review in Sec. II.1. We then
move on to one of the currently implemented methodologies used to test GR with
GW data (the LVC implementation of the ppE framework, which they dub a
parametrized inspiral test) in Sec. II.2.

###  II.1 Bayesian Analysis of GW Data

The foundation of any Bayesian analysis is, of course, Bayes’ theorem. In the
context of GW analysis, this can be written as

$p(\bm{\theta}|D)=\frac{p(\bm{\theta})p(D|\bm{\theta})}{p(D)}\,,$ (1)

where $p(\bm{\theta}|D)$ is the posterior probability of the source parameters
$\bm{\theta}=\\{\theta_{1},\theta_{2},\ldots,\theta_{N}\\}$ given some dataset
$D$, $p(\bm{\theta})$ is the prior probability of $\bm{\theta}$,
$p(D|\bm{\theta})$ is the likelihood that the data $D$ is described by the
model and parameter vector $\bm{\theta}$, and $p(D)$ is the evidence for the
model. The posterior distribution is the quantity of interest in this study,
as we are primarily interested in deducing the values additional, beyond-GR
parameters can take while still being consistent with data observed in our
Universe. The evidence can also be of interest when calculating Bayes’ factors
between models, but we will not be focused on model selection in this study.
The information gained from the current experiment at hand is encoded in the
likelihood distribution, while the beliefs about the different parameters held
before the experiment began is embodied by the prior distribution.

In the case of inference with GW data, the likelihood function is calculated
by matched-filtering. With this method, a template response function defined
by a model and a set of parameters is compared against the data. This function
can be succinctly written in the Fourier domain (assuming Gaussian, stationary
noise) as

$\log
p(D|\bm{\theta})=-\frac{1}{2}\sum_{i}^{N}\left(D_{i}-h_{i}|D_{i}-h_{i}\right)\,,$
(2)

for $N$ detectors, where $h_{i}$ is the response template of the $i$-detector
given a GW metric perturbation (the perturbation template contracted on to the
$i$-th detector $h_{i}=h_{+}F_{+,i}+h_{\times}F_{\times,i}$) and the inner
product is defined as

$\left(A|B\right)=4\operatorname{Re}\int_{0}^{\infty}\frac{AB^{\ast}}{S(f)}df\,,$
(3)

where $({}^{\ast})$ denotes complex conjugation and $S(f)$ is the one-sided
power spectral density of the noise in the detector.

As a general rule in Bayesian inference, one should be careful in how the
prior information is specified. Either in the case of encoding information
from past experiments or through mathematical intuition, or in how one
attempts to encode one’s “ignorance” about the parameter in question, an
overly restrictive prior can lead to biased posterior distributions. Standard
choices used in GW science include a uniform prior in volume, a uniform prior
on the component masses and the spin magnitudes, with a prior for the spin
directions, sky location, and binary orientation uniform on the unit two-
sphere. These choices of prior encode a certain amount of basic information
through certain assumptions we have made about these parameters (e.g. that the
universe is uniform and isotropic on large scales, or that the masses of
compact objects cannot be negative), but are considered to be “uninformative.”

Ultimately, we seek to determine the form of the posterior defined in Eq. (1),
because the degree to which the final marginalized posterior differs from the
prior determines how much new information we have gained by performing the
experiment. Once the posterior has been calculated, we can deduce many useful
quantities, such as the confidence regions and maximum likelihood values for
various parameters. However, without making any simplifying assumptions, this
function is not analytic. Therefore, we are typically forced to explore these
distributions through sampling techniques like Markov-Chain Monte-Carlo (MCMC)
algorithms. Analyses such as these produce samples from the distribution of
interest, which can then be binned into a histogram representation of the
distribution.

Throughout this work, we will sample the likelihood, and therefore construct
the posterior, through an MCMC exploration of the parameter space using
software written for and implemented in previous works [35]. The Bayesian
exploration of parameter space is based on the Metropolis-Hasting algorithm,
and implements a variety of techniques to optimize the sampling efficiency,
which we will briefly outline here. A critical component of the sampler is the
use of parallel-tempering [47, 48], in which multiple MCMC chains are run in
parallel. The first chain in the series samples from an unmodified posterior
distribution, but all of the other chains in the series sample from a
“tempered” distribution, in which the likelihood is raised to some power,
$1/T_{i}$, where $T_{i}\geq 1$. The modified posterior is given by

$p_{i}(\bm{\theta}|D)\propto p(\bm{\theta})p(D|\bm{\theta})^{1/T_{i}}\,,$ (4)

for some “temperature” $T_{i}$. Periodically, the chains are allowed to
exchange information by proposing a swap of the parameter coordinates of two
of the chains. Each temperature ladder has a certain set of $\\{T_{i}\\}$
comprised of 15 (30) separate chains for the 2g (3g) injections (where the 2g
and 3g injection parameters are discussed in Sec. IV). The temperatures run
from $T_{0}=1$ (untempered) to $T_{14~{}(29)}\approx\infty$ (effectively
sampling the prior). To increase efficiency, we run 14 (5) copies of the above
set of chains in parallel, so that we have multiple chains sampling at each
temperature, $T_{i}$. All chains in the ensemble are allowed to swap with all
other chains to improve convergence. Samples are then harvested from the $T=1$
chains, which are sampling from the unmodified, target posterior. While this
does introduce some correlations between the samples from the cold chains, we
harvest more than enough samples to mitigate the impact of this, as discussed
below. Furthermore, correlations of this type are analogous to those that
would be present with other ensemble samplers, which base proposal
distributions on the position of other chains in the ensemble.

In terms of the proposals we have implemented, we used two types: steps drawn
from a random normal distribution centered at the current location and steps
along the eigenvectors of the Fisher information matrix. We “burn-in”
approximately 80,000 (28,000) samples per chain, during which we allow the
sampler to optimize itself to the target distribution by tuning the step
widths, continuously updating the Fisher matrix, and optimizing the
temperature spacing of the various chains using the algorithm by Vousden
et.al. [49]. After this initial burn-in phase, we then fix all the parameters
of the sampler and continue to run for an additional 100 independent samples
per untempered chain, as determined by the auto-correlation length, to make
sure that indeed the chains have exited the burn-in phase. Then, we proceed to
harvest samples that we will use to build the distribution, while thinning out
by the autocorrelation length. We sample until we have harvested approximately
4,000 independent samples per untempered chain. This produces on the order of
56,000 (20,000) total, independent samples per 2g (3g) injections. With the
samples drawn, we calculate the 1-$\sigma$ confidence intervals by calculating
the $16\%$ and $84\%$ quantiles from the data using Numpy routines [50], and
taking the difference then dividing by 2.

###  II.2 Current Parametric Tests of GR

With the general framework outlined above, we can now discuss some of the
current methodologies used in searches for physics beyond GR in GW data. In
the course of this paper, we will always work with small deformations away
from GR, consistent with our experience that GR describes gravity extremely
well up to (and possibly beyond) the most extreme energy scales we have
observed so far. In the spirit of this philosophy, we use a base GR waveform
model (IMRPhenomD [51, 52]) and append modifications on top through a
phenomenological framework. We could have used a more advanced base GR model,
like IMRPhenomPv3 [53], but doing so will probably not impact our conclusions
qualitatively.

The LVC parametrized inspiral tests are a realization of the ppE framework
[12], restricted to single-parameter deformations and deformations that enter
at half-integer PN orders, and thus, there are two “standard” ways to
implement them. In the first implementation, the modifications are added as
fractional deviations at each PN order directly to the GW phase. Explicitly,
let us write the GW Fourier phase of the dominant mode as as [51, 52, 36]

$\displaystyle\Psi_{\rm GW}(f)$
$\displaystyle=\frac{3}{128\eta}v^{-5}\left(1+\sum_{i=2}^{7}p_{i}(\bm{\Xi})v^{i}\right)$
$\displaystyle+\frac{3}{128\eta}v^{-5}\left(\sum_{i=5}^{6}p_{i,l}(\bm{\Xi})\log\\!\left(v^{3}\right)\,v^{i}\right)\,,$
(5)

where $v\equiv(\pi mf)^{1/3}$, $m\equiv m_{1}+m_{2}$ is the total mass for a
binary with component masses $m_{1}$ and $m_{2}$, and $f$ is the GW frequency
of the dominant mode, while $p_{i}\equiv p_{i,\mbox{\tiny
GR}}(1+\delta_{i}^{\bar{i}}\,\delta p_{i})$ and $p_{i,l}\equiv
p_{i,l,\mbox{\tiny GR}}(1+\delta_{i}^{\bar{j}}\,\delta p_{i,l})$. The
parameters $p_{i,\mbox{\tiny GR}}$ are the coefficients of the PN expansion in
GR at the $(i/2)$-th PN order, and $p_{i,l,\mbox{\tiny GR}}$ are the
coefficients proportional to logarithmic terms; these coefficients are purely
functions of the source parameters, $\bm{\Xi}$, and not of frequency. The
parameters $\delta p_{i}$ and $\delta p_{i,l}$ are the non-GR deformation
parameters that one attempts to estimate or constrain [28, 23, 22, 30, 29].
The constants $\bar{i}$ and $\bar{j}$ determine the PN order of the single
deformation that is turned on, with the Kronecker delta guaranteeing that all
other deformations are turned off. This is the parametrization that has been
largely adopted by the LVC in their parametrized inspiral tests of GR (c.f.
[24, 25, 26, 27, 31, 13]).

A second implementation is to write the GW Fourier phase of the dominant mode
as

$\displaystyle\Psi_{\rm GW}(f)$
$\displaystyle=\frac{3}{128\eta}v^{-5}\left(1+\sum_{i=2}^{7}p_{i,\mbox{\tiny
GR}}(\bm{\Xi})\,v^{i}\right)$
$\displaystyle+\frac{3}{128\eta}v^{-5}\left(\sum_{i=5}^{6}p_{i,l,\mbox{\tiny
GR}}(\bm{\Xi})\log\\!\left(v^{3}\right)\,v^{i}\right)$
$\displaystyle+\beta\,\mathcal{U}^{(b-5)}\,,$ (6)

where $\mathcal{U}=(\pi\mathcal{M}f)^{1/3}$, $\mathcal{M}=\eta^{3/5}m$ is the
chirp mass of the binary, with $\eta=m_{1}m_{2}/m^{2}$ the symmetric mass
ratio, $\beta$ is a deformation parameter and $b$ is a real number one fixes
before carrying out parameter estimation. This is closer to the original ppE
framework [19] and it is employed by researchers because it allows for a
slightly more direct mapping to constraints on specific modified theories [17,
20, 21, 32]. As it should be painstakingly clear by comparing Eqs. (5) and
(6), these two implementations are equivalent to each other with
$\bar{i}=b-5$, $\bar{j}=0$, a reparametrization between $\delta p_{i}$ and
$\beta$, and the transformation of their priors, as shown in detail in [12].

Regardless of the preferred implementation, both of the ones presented above
focus on a single modification at a time (e.g. a single $\beta$ or $\delta
p_{i}$ or $\delta p_{i,l}$ at a time). This is more a choice made out of
practicality than one made because of a physically-motivated reason. In
reality, solving for the phase of the gravitational waveform, even in GR,
leads to a series with an infinite number of terms when using the PN
formalism, as shown in Eqs. (5) and (6). As these generic parametrizations
represent unknown corrections to these terms in the PN series (whether in the
GW phase directly or in the energy flux, binding energy or dispersion
relations), there should also be an infinite number of deformations when the
calculation is done for any specific modified theory of gravity. One would
therefore expect a Fourier waveform phase of the form of Eq. (5) but without
the Kronecker deltas in the definitions of $p_{i}$ and $\delta p_{i,l}$ for
the first implementation. For the second implementation shown in Eq. (6), the
natural generalization would be

$\displaystyle\Psi_{\rm GW}(f)$
$\displaystyle=\frac{3}{128\eta}v^{-5}\left(1+\sum_{i=2}^{7}p_{i}(\bm{\Xi})v^{i}\right)$
$\displaystyle+\frac{3}{128\eta}v^{-5}\left(\sum_{i=5}^{6}p_{i,l}(\bm{\Xi})\log\left(\pi
mf\right)v^{i}\right)$
$\displaystyle+\beta\,\mathcal{U}^{(b-5)}\left(1+\sum_{i=1}^{7}\beta_{i}\mathcal{U}^{i}\right)\,,$
(7)

where we now include $\beta$ along with higher terms in the PN expansion,
$\beta_{i}$

To be maximally robust, one would then think that one should simultaneously
fit for all the $\delta p_{i}$ and $\delta p_{i,l}$ or all the $\beta_{i}$ and
$\beta$ parameters at all known PN orders (for waveforms used in current
parameter estimation analysis, this would go up to 3.5PN order). However, when
this was attempted in the past, the constraints on any one $\delta p_{i}$ or
$\beta_{i}$ parameter quickly degraded to the point that the posterior
devolved into the prior [13], and no new information could be gleaned about
the underlying physics. This was predicted back in [20] and then verified
recently in [13]. As is commonly understood, this is due to overfitting [23,
34, 54]. By increasing the dimensions of the parameter space and allowing a
large enough prior range that significant degeneracies surface, any insight
into the individual modifications is washed out by covariances. But if
information is lost when including multiple deformations, are the constraints
found while using a single deformation robust, as these covariances are
completely neglected in this restricted case? If all one wishes to establish
is whether a deviation exists, then previous work has shown that indeed a
single-term deformation is sufficient [20]. Here, however, we are not
concerned with detecting a deviation, but rather determining the accuracy to
which deformation coefficients can be estimated, and thus, constraints can be
placed.

##  III Improved Parametric Tests of GR

From the standpoint of testing for new physics, is the situation truly so
bleak? There are some critical details that are completely ignored by the
overly agnostic approach of varying over all deformation parameters
simultaneously, and which one could may be use to perform robust parametric
tests. We discuss the critical, additional information we will incorporate
into our analysis in Sec. III.1 and the details of the exact implementation in
Sec. III.2.

###  III.1 Restrictions Based on Mathematical Structure of Modified Theories

There exists at least two pieces of basic information that one can infer and
use about the exact form of the deformations introduced into the GW Fourier
phase by real modified theories of gravity.

The first piece of information concerns the PN structure of the deformations.
Let us consider theories in which the early inspiral waveform can be
calculated with the PN approximation, as we do in GR, where the phase is
written as a series in the orbital velocity. If this is the case, there is a
hierarchy in the magnitude of the modifications at each PN order, so that the
PN series is convergent in certain regimes of parameter space. For example,
using the parametrization of Eq. (5)

$\delta p_{i}p_{i,{\rm GR}}v^{i}\gg\delta p_{i+1}p_{i+1,{\rm GR}}v^{i+1}\,,$
(8)

while using the parametrization of Eq. (7)

$\beta_{i}v^{i}\gg\beta_{i+1}v^{i+1}\,.$ (9)

These conditions then lead to the asymptotic inequalities $\delta
p_{i+1}\,v\ll\delta p$ and $\beta_{i+1}\,v\ll\beta_{i}$. The implication here
is that allowing the $\delta p_{i}$ or the $\beta_{i}$ parameters at each PN
order to vary over completely independent ranges is not representative of any
PN-compatible theory of gravity beyond GR. For example, if the 2 PN
deformation coefficient $\delta p_{4}$ or $\beta_{4}$ had a prior range
$(-10^{7},+10^{7})$, then there would be many choices for which the 2PN terms
would completely dominate over the 1PN and the Newtonian term, rendering the
PN approximation inaccurate in the regime of velocities of interest to ground-
based detectors.

Whether the PN series is convergent, asymptotic, or divergent has not been
formally established due to the lack of an exact solution one can compare
against. For comparable masses, Ref. [55, 36] computed the radius of
convergence of the PN series using the PN binding energy, and up to 2PN order
the author found it to be $v/c\approx 1/2$. In the extreme mass-ratio limit,
Ref. [56] have shown the PN series is asymptotic, and determined the optimal
radius of convergence as a function of PN order to be $v/c\approx(1/5,1/3)$
depending on the PN order studied. Either way, the growth of the coefficients
of the PN approximation cannot grow with PN order too rapidly if the solution
is to remain accurate in the late inspiral (at velocities of roughly $1/3c$).
This is true not just in GR, but also in modified theories of gravity, since
it is a mathematical statement about the validity of a series solution.

The second piece of information one can use relates to the fact that the
mathematical structure of all physical theories guarantee that each modified
GR term will be a function of the coupling constants of the theory. Let us
restrict attention to modified theories that have a continuous GR limit, i.e.
they reduce to GR continuously as their coupling constants tend to a specific
value. Then, for small deformation from GR, each modified GR term can be
Taylor expanded in the coupling constants of the theory about their GR values.
To make this concrete, let us focus on theories with a single coupling
constant, and let this constant be called $\ell$. Then, Taylor expanding the
deformation coefficients at any given PN order, we have

$\delta p_{i}(\ell,\bm{\Xi})=\delta\bar{p}_{i}(\bm{\Xi})\,\ell^{p}\,,$ (10)

or

$\beta_{i}(\ell,\bm{\Xi})=\bar{\beta}_{i}(\bm{\Xi})\,\ell^{p}\,,$ (11)

where $\delta\bar{p}_{i}(\bm{\Xi})$ and $\bar{\beta}_{i}(\bm{\Xi})$ are
functions of the source parameter vector $\bm{\Xi}$ only (like the masses,
spins, etc), and $p$ is a positive real number. Note that although the
functional dependence on the source parameters may change with PN order, the
functional dependence on the coupling parameter will not in the small
deformation limit (as otherwise the Taylor expansion of the full series would
not be well-defined). One can now introduce a new parameter $\gamma=\ell^{p}$
to write all deformation coefficients as linear functions in $\gamma$ to
leading-order in small deformations. This conclusion is not only valid for
theories with a single coupling constant, but rather, it can be generalized to
theories with multiple couplings if there is one that dominates due to others
being already constrained by other experiments or observations.

In summary, the PN expansion of most theories of gravity should give a GW
Fourier phase $\Psi_{\rm GW}$ that takes the form

$\displaystyle\Psi_{\rm GW}(f,\gamma,\bm{\Xi})$ $\displaystyle=\Psi_{\rm
GR}(f,\bm{\Xi})+\gamma\,\delta\psi_{0}(\bm{\Xi})\,v^{b}$
$\displaystyle\times\left[1+\delta\psi_{2}(\bm{\Xi})\,v^{2}+\delta\psi_{3}(\bm{\Xi})\,v^{3}+\ldots\right]\,,$
(12)

where we use the new parameter $\delta\psi_{i}$ to denote the deformation
amplitudes because we now choose to use the PN expansion parameter $v=(\pi
mf)^{1/3}$; this is a slight deviation from the notation of the standard ppE
PN expansion, which uses the parameter $\mathcal{U}=(\pi\mathcal{M}f)^{1/3}$.
Here, $\Psi_{\rm GR}$ is the GW phase in GR, and $\gamma$ generically
represents the coupling constant in a given theory of gravity. The subscripts
$i$ in $\delta\psi_{i}$ represent the order in $v$ or the $i/2$ PN order at
which the deformation enters the GW Fourier phase relative to the leading PN
order term in the non-GR sector. The 0PN modification could enter the phase at
any PN order relative to GR, as controlled by the $b$ parameter. As required
by the first piece of information presented above, the range of the
$\delta\psi_{i}$ functions above must satisfy
$\delta\psi_{i+1}\,v\ll\delta\psi_{i}$. As required by the second piece of
information, the coupling constant of the modified theory shows up as a pre-
factor that multiplies all deformations, so that in the GR limit, the non-GR
terms vanish. When explicitly attempting to search for physics beyond GR, the
real information comes from the first term in the series, where the coupling
constant $\gamma$ appears. Therefore, estimates on our ability to discern new
information about physics should be quantified purely in terms of constraints
on this constant, regardless of the number of higher order terms we include.

One of the first things one notices from Eq. (III.1) is that the
$\delta\psi_{i}$ are not constants but rather functions of the source
parameters $\bm{\Xi}$. Given a specific modified theory, one knows what these
functions are, and therefore, the $\delta\psi_{i}$ are not new parameters. In
fact, in this case, the only new parameter to the model is simply the coupling
constant $\gamma$. If one is attempting to carry out an agnostic test of GR,
however, one is not privy to the functional form of $\delta\psi_{i}$. In this
case, one can replace Eq. (III.1) via

$\displaystyle\Psi_{\rm GW}(f,\gamma,\bm{\Xi})$
$\displaystyle=\Psi_{\mbox{\tiny
GR}}(f,\bm{\Xi})+\gamma\,\delta\bar{\psi}_{0}\,v^{b}$
$\displaystyle\times\left[1+\delta\bar{\psi}_{2}\,v^{2}+\delta\bar{\psi}_{3}\,v^{3}+\ldots\right]\,,$
(13) $\delta\bar{\psi}_{i}=\frac{\int\delta{\psi}_{i}\;d^{N}{{\Xi}}}{\int
d^{N}{{\Xi}}}\,,$ (14)

where $N$ is the dimensionality of the source parameter space and $d^{N}{\Xi}$
is the parameter space volume factor. Clearly then, the coefficients
$\delta\bar{\psi}_{i}$ are now constants that become new independent
parameters to search over. Moreover, since $\gamma$ enters multiplied by
$\delta\bar{\psi}_{0}$, these two parameters are 100% degenerate, so one must
employ the reparametrization $\bar{\gamma}=\gamma\,\delta\bar{\psi}_{0}$. The
waveform parameter space is then enlarged to $\bm{\Xi}\cup\bm{\Theta}$, where
$\bm{\Theta}=\\{\bar{\gamma},\delta\bar{\psi}_{2},\delta\bar{\psi}_{3},\ldots\\}$.
As before, the GR deviations are all still controlled by a single parameter
$\bar{\gamma}$, while the range of the deformation coefficients is restricted
by $\delta\bar{\psi}_{i+1}\,v\ll\delta\bar{\psi}_{i}$.

###  III.2 Implementation of Improvements
in Parametric Inspiral Tests

Using the representation of Eq. (III.1) and enforcing convergence of the
series, we can now attempt to include multiple deformations at once. However,
care must be taken in the way one enforces convergence, as the criteria
$\delta\bar{\psi}_{i+1}\,v\ll\delta\bar{\psi}_{i}$ depends on the orbital
velocity $v$.

To do this, we will implement a specialized prior that ensures the series of
$\delta\bar{\psi}_{i}$ are convergent at a “reasonable” value of $v$. That is,
we choose the prior to be

$|\delta\bar{\psi}_{i}|>|\delta\bar{\psi}_{(i+1)}\left(v_{\text{eval}}\right)|\,,$
(15)

for $i\in[2,3,4,6,7]$, where $v_{\text{eval}}=(\pi mf_{\text{eval}})^{1/3}$;
the exact value of $f_{\text{eval}}$ and thus $v_{\text{eval}}$ that we choose
are discussed below. For the 1PN term, the prior must be handled separately,
and we choose it to be

$|\delta\bar{\psi}_{2}\left(v_{\text{eval}}\right)^{2}|<1\,,$ (16)

as this is the next-to leading order (NLO) term due to our omission of the
$0.5$PN term.

We must now choose the value of the velocity, or equivalently the frequency,
at which to evaluate this prior. We do so by choosing the GW frequency and
orbital velocity that would correspond to the waves emitted by a binary at an
orbital separation of $r_{12,\rm eval}=100m$ at leading-order in PN theory.
Note that $r_{12,\rm eval}$ is _not_ the formal radius of convergence of the
series, assuming the series were convergent. Rather, this quantity represents
the orbital separation at which we are confident that adding higher order PN
terms improves the PN approximation to the (unknown) exact answer. More
aggressive (smaller) choices of the orbital separation could be viable options
as well, but we will remain conservative in this work, as we will explain
below. Using the Newtonian version of Kepler’s third law for quasi-circular
orbits, we then have

$\displaystyle v_{\rm eval}$ $\displaystyle=\left(\frac{m}{r_{12,\rm
eval}}\right)^{1/2}=10^{-1}\,,$ (17) $\displaystyle f_{\text{e}val}$
$\displaystyle=\frac{(r_{12,\rm eval}/m)^{-3/2}}{(\pi m)}\approx\frac{3\times
10^{-4}}{m}\,.$ (18)

This choice in frequency is inspired by the convergence properties of the GW
phase in GR, as studied in [56] and revisited below. GR is the obvious model
for us to base our new prior on, as it is much better understood than any
modified theory and it is our null hypothesis.

As we now show, with this choice of separation (or equivalently frequency or
velocity) the coefficients of the PN expansion in GR present a convergent
structure for most mass ratios. The phase in GR (via the Taylor F2 waveform
[57, 58, 59], accurate to 3.5PN) was already presented in Eq. (5), and the
exact forms of the $p_{i}$ and $p_{i,l}$ parameters can be found in various
places in the literature (c.f. [36, 51, 52]). For a non-spinning binary in GR,
the ratios $p_{i}/(p_{i+1}v_{\rm eval})$ only depend on the symmetric mass
ratio $\eta$ alone. If these ratios are larger than unity, then the PN
coefficient exhibit a convergent structure111This is in the sense that $v_{\rm
eval}$ would then be within the radius of convergence of the series if the
series were convergent. If the series is asymptotic, then this would be an
indication that one is evaluating the series in the regime where the
asymptotic series is a good approximation to the exact answer.. Figure 1
presents this ratio as a function of $\eta$, where we observe that our choice
of $v_{\rm eval}$ leads to a convergent structure of the series for all
$\eta>0.09$. The spikes at around $\eta\sim 0.053$ are because $p_{6}$
vanishes, but this occurs outside the comparable-mass regime, which is what we
focus on here.

Figure 1:  Ratio of the PN coefficients of the Fourier GW phase (Eq. (5)) for
a non-spinning binary as functions of the symmetric mass ratio $\eta$. The
ratios are evaluated at an orbital separation of $r_{\rm eval}=100m$,
corresponding to $v_{\rm eval}=0.1$. Each line corresponds to a different pair
of coefficients, and the dotted black line simply identifies the threshold we
use for convergence ($|p_{i}/(p_{i+1}u_{\rm eval}|>1$). The spikes shown here
correspond to when the coefficient $p_{6}=0$ at $\eta\sim 0.053$. Observe that
all ratios are above unity for symmetric mass ratios $\eta>0.09$, indicating
that our choice of $v_{\rm eval}$ has led to a convergent structure of PN
coefficients.

As the parametrization in Eq. (III.1) has not been used extensively in
previous literature, we will also study a related parametrization, closely
associated with the standard ppE formalism. Let us then denote absolute
deviations as $\bar{\Delta}_{i}$, and relative deviations as
$\delta\bar{\psi}_{i}$ (as introduced in Eq. (III.1)) for the $i/2$-th PN
order. These new $\bar{\Delta}_{i}$ parameters then enter the GW phase at
various PN orders via

$\displaystyle\Psi_{\rm GW}(f,\gamma,\bm{Xi})$ $\displaystyle=\Psi_{{\rm
GW},GR}(f,\bm{\Xi})+\bar{\gamma}\,v^{b}+$
$\displaystyle\bar{\Delta}_{2}\,v^{(2+b)}+\bar{\Delta}_{3}\,v^{(3+b)}+\ldots\,,$
(19)

where each modification ($\bar{\gamma},\bar{\Delta}_{i}$) now explicitly
depends on the coupling constant of the theory. The mapping between Eq.
(III.2) and Eq. (III.1) is simply given by

$\bar{\Delta}_{i}=\delta\bar{\psi}_{i}\times\bar{\gamma}\,,$ (20)

and $\bar{\gamma}$ is unchanged.

The priors we have defined above, although simple to state mathematically, are
non-trivial. To illustrate this, we sampled from the prior directly using a
variation of the sampling methods outlined in Sec. II.1. In effect, we drew
each of the deformation parameters from a uniform distribution $U(-B,+B)$ with
some boundaries $-B$ and $+B$. For the leading order and NLO deformation,
$|B|=100$, while for the other, higher order parameters, $|B|$ is set to a
value large enough to ensure the final distributions are not affected by
$|B|$. With these random draws, we excluded any sample that violated Eq. (15).
Although each parameter began as a uniform distribution, our prior in Eq. (15)
imposed very non-trivial structure to the histogramed distributions, as the
boundary of the prior for each deformation depends on the value of other
deformations in the series.

The priors on $\delta\bar{\psi}_{i}$ with 6 deformations and beginning at 0PN
order relative to GR (proportional to $v^{-5}$) are shown in the left panels
of Fig. 2, while those transformed to the absolute deformations
$\bar{\Delta}_{i}$ are in the right panel. Observe that the prior on the
leading PN order term in the series (which encodes most of the information we
are interested in) is indeed flat, because Eq. (15) does not affect it.
Observe also that the larger PN order the relative deviation, the larger the
range of the prior. This range, however, is indeed finite, since if it were
any larger, then the relative deviations would begin to affect the convergent
structure of the PN series. Finally, observe that although the priors on
$\delta\bar{\psi}_{i}$ seem to be pushing modifications away from GR
($\delta\bar{\psi}_{i}=0$), this is just an artifact of looking at relative
deviations; the priors on the absolute modifications $\bar{\Delta}_{i}$ are
fully consistent with GR ($\bar{\Delta}_{i}=0$).

Figure 2:  The one dimensional, marginalized prior on a set of deformations.
The leading order deformation enters at Newtonian order relative to GR, and
each of the subsequent deformations are labeled as the $i/2$PN order relative
to this leading order deformation. The priors for this parametrization are
non-trivial, and can be misleading. They appear to disallow the GR limit
$\delta\bar{\psi}_{i}=0$ for several of the different PN orders. This is
resolved when looking at the transformed priors of the right panel, rewritten
in terms of the absolute modifications, $\bar{\Delta}_{i}$.

##  IV Experimental Design

To determine the impact of these higher order deformations on our inferences
about fundamental physics, we now conduct a series of parameter estimation
analyses on synthetic signals (injections) with a variety of different
recovery models. Each of the recovery models has the same base GR waveform
(IMRPhenomD), but they differ in the different number of phase deformations
overlaid on top of that base model. For each model, the modifications are
added in ascending order beyond leading, beginning with the 1PN term (relative
to the leading order deformation), and increasing by 0.5PN order up to 3.5PN
order relative to GR, or relative to the leading order term, which ever
criteria is met first. We skip terms of 2.5PN order relative to GR in all
models, as they are 100% degenerate with the coalescence phase (which is an
arbitrary constant). When skipping terms in the series due to these total-
degeneracies, we also skip them in the enforcement of the prior. For example,
if we were to skip the $\delta\bar{\psi}_{i}$ term in the series because it
was degenerate with the coalescence phase, we would update our prior for
$\delta\bar{\psi}_{i-1}$ and $\delta\bar{\psi}_{i+1}$ to satisfy
$\delta\bar{\psi}_{i-1}\geq\delta\bar{\psi}_{i+1}v^{2}$. Therefore, every
model is uniquely defined by the PN order of the leading deformation (relative
to GR) and the number of subsequent terms.

We will specifically target leading order GR deformations that enter at
$[-1,0,1,1.5,2,3,3.5]$PN orders, relative to Newtonian order in GR, and we
include up to 5 additional deformations (again, never exceeding 3.5PN order
relative to Newtonian order in GR). For example, one model might have a
leading order deformation at 1PN order relative to Newtonian order in GR, with
3 higher PN order deformations at $[2,3,3.5]$PN orders (again relative to
Newtonian order in GR). As another example, another model might have a leading
order deformation that enters at 1.5PN order relative to the Newtonian term in
GR, with only a single higher PN order deformation. This would leave us with
deformations at $[1.5,3]$PN order relative to Newtonian order in GR.

As expected, all inferences will depend on the astrophysical source that
produced the GW we are assuming has been detected, as well as the detector
network that observes the signals. To attempt to explore these aspects as much
as possible (while still keeping the scope of this work computational
tractable), we focus on two sources: a “heavy” source and a “light” source,
whose properties are listed in Table 1. The spin configurations also differ
slightly, with one binary having both spins aligned with the orbital angular
momentum, while the other binary has one spin anti-aligned. In both cases, we
focus on binary black hole inspirals, and do not consider neutron star
inspirals or mixed binaries. We expect the qualitative conclusions we will
find to also hold for those systems.

For the detector networks, we also focus on two configurations. For the first
configuration, and as a proxy for a 2g network, we use a network comprised of
LIGO Hanford [60], LIGO Livingston [60], and Virgo [61]. For their
sensitivities, we use analytic approximations to the aLIGO design sensitivity
[62] and the first phase of Advanced Virgo’s sensitivity estimate [62]. For
the second configuration, and as a proxy for a 3g network, we use a network
comprised of Cosmic Explorer (CE) [63] and the Einstein Telescope (ET) [64].
For their sensitivities, we use the first phase of the CE noise curve [65] and
the ET-D configuration of the ET noise curve [66]. To reduce the impact of the
uncertainties concerning the noise curves we use, the luminosity distances of
the injected sources are all scaled such that the SNR is exactly 20 (when
observed by the entire 2g detector network). The distances are then left fixed
at these values when we transitioned to the 3g network, so as to isolate the
impact of a pure boost in SNR.

With all these considerations in mind, we then have 21 separate analyses (one
for each combination of multi-deformation-parameter ppE model) for each
detector network and source combination , for which there are four
combinations. This comes out to a total of 84 separate Bayesian parameter
estimation studies, whose results we summarize next.

Source Identifier |  $m_{1}$ ($M_{\odot}$) |  $m_{2}$ ($M_{\odot}$) | $\chi_{1}$ | $\chi_{2}$ |  $D_{L}$ (Mpc) | SNR2g | SNR3g
---|---|---|---|---|---|---|---
Heavy | 25 | 25 | 0.3 | -0.1 | 1613 | 20 | 517.798
Light | 5 | 5 | 0.3 | 0.1 | 466 | 20 | 496.355
Table 1:  Choices of GW sources for injection campaigns. The source parameters
are the following: $m_{1}$ and $m_{2}$ are the masses of the larger and
smaller black holes, respectively, $D_{L}$ is the luminosity distance from
Earth to the source, $\chi_{i}$ is the aligned, dimensionless spin of the
$i$-th black hole, SNR2g and SNR3g are the SNRs of the source as measured by a
2g and 3g detector network respectively. The “heavy” source and the light
“source” both have an SNR of 20 (for the 2g detector network), but have
different total masses and different spin configurations. All source
parameters were kept the same between the analyses involving the 2g and the 3g
detector networks.

##  V Bayesian Results

Let us now present the results of our experiments. Because the leading PN
order deformation contains all the relevant constants controlling the
magnitude of GR deviations, we focus on constraints on this leading-order
term, presenting results both visually and in tabular form.

Figure 3: Marginalized $1\sigma$ constraints on the leading PN order
deformations for the “heavy” source (left) and the “light” source (right)
using the 2g network (see Sec. IV for details) as a function of the number of
sub-leading PN order terms added in the modified sector (in ascending order).
Each line corresponds to modifications that start at different leading PN
order. The functional form of GW phase is given in Eq. (III.1), and the priors
used are those presented in Sec. III. The top panels show the $1\sigma$
constraint on the leading PN order deformation, while the lower panel shows
the strengthening factor, as defined in Eq. (21). Observe how the lines in the
upper panels trend downward, which means that constraints with just a single
parameter deviation are conservative (i.e. adding higher PN order corrections
strengthens the constraint obtained with a single parametric deviation).
Moreover, observe how the slope of the lines are small, which means that the
strengthening of the constraint is mild, with improvements of at most roughly
one order of magnitude.

Terms/LO | -1 | 0 | 1 | 1.5 | 2
---|---|---|---|---|---
1 | $4.1\times 10^{-4}$ | $6.6\times 10^{-3}$ | $9.7\times 10^{-2}$ | $4.8\times 10^{-1}$ | $3.1$
2 | $2.4\times 10^{-4}$ | $5.3\times 10^{-3}$ | $1.6\times 10^{-1}$ | $2.6\times 10^{-1}$ | $2.5$
3 | $2.2\times 10^{-4}$ | $4.2\times 10^{-3}$ | $4.8\times 10^{-2}$ | $7.1\times 10^{-2}$ |
4 | $1.4\times 10^{-4}$ | $3.0\times 10^{-3}$ | $1.1\times 10^{-2}$ |  |
5 | $6.1\times 10^{-5}$ | $8.5\times 10^{-4}$ |  |  |
6 | $3.6\times 10^{-5}$ | $2.5\times 10^{-4}$ |  |  |

Terms/LO | -1 | 0 | 1 | 1.5 | 2
---|---|---|---|---|---
1 | $8.3\times 10^{-6}$ | $6.8\times 10^{-3}$ | $3.3\times 10^{-2}$ | $2.5\times 10^{-1}$ | $2.7$
2 | $8.9\times 10^{-6}$ | $2.0\times 10^{-3}$ | $3.4\times 10^{-2}$ | $1.1\times 10^{-1}$ | $1.3$
3 | $9.4\times 10^{-6}$ | $1.2\times 10^{-3}$ | $2.5\times 10^{-2}$ | $8.0\times 10^{-2}$ |
4 | $7.8\times 10^{-6}$ | $8.2\times 10^{-4}$ | $9.9\times 10^{-3}$ |  |
5 | $6.4\times 10^{-6}$ | $3.0\times 10^{-4}$ |  |  |
6 | $5.1\times 10^{-6}$ | $1.6\times 10^{-4}$ |  |  |

Table 2:  Marginalized $1\sigma$ constraints on the leading PN order
deformations for the “heavy” source (left) and the “light” source (right)
using the 2g network (see Sec. IV for details). The columns represent the PN
order (relative to the Newtonian term in GR) at which the GR deformation is
first introduced. The rows corresponds to the number of PN corrections that
are added on top of the leading PN order one in ascending PN order. The
functional form of the phase is given in Eq. (III.1), where the priors are
presented in Sec. III. Observe the (order of magnitude) consistency of the
constraints as one includes more and more sub-leading PN order deviations
(i.e. as one moves down the column for a fixed leading-order term). This
indicates that higher PN order corrections to the current modified gravity ppE
waveforms will not invalidate current bounds placed on modified theories with
only leading PN order deformations.

The results for the 2g network injections are shown in Fig. 3 and Table 2 for
the “heavy” and “light” sources. The top panels show the marginalized
$1\sigma$ constraint on the leading PN order deformation as a function of the
number of terms included in the series, beginning with 1 (only the leading PN
order deformation) and up to 6 total terms. The lower panel shows the
strengthening factor, which we define here as

$(\text{strengthening factor})=\frac{\sigma_{\rm Single}}{\sigma_{\rm
Multiple}}\,.$ (21)

If this number is larger (smaller) than unity, then adding higher PN order
deformations strengthens (weakens) the constraint one obtains by using a
single parameter deformation. The bottom panels of Fig. 3 show the
strengthening factor as a function of PN order.

Two main conclusions can be drawn from these figures and tables. First,
observe how all curves in the top panels trend downwards (have a negative
slope). This implies that adding higher PN order corrections to the modified
sector strengthens the constraints one would have gotten if one included only
a single parametric deformation. This is confirmed in the bottom panels, which
shows the strengthening factor is always greater than one. In this sense,
single-parameter constraints are therefore conservative. Second, observe how
the slopes of the curves are small. This implies that the strengthening one
obtains is somewhat mild, with improvements of only up to a factor of 40 in
the most extreme case, as shown in the bottom panels of the figure.

These results and conclusions may seem counter-intuitive and, in fact,
opposite to what one would expect. Adding additional, independent parameters
to a model has been shown to increase degeneracies, and therefore deteriorate
our ability to estimate any given parameter [23, 34, 20, 13]. The reason this
does not happen here is because of the series structure of the deformations in
the model. Equation (III.1) is different from how other work has modeled
deviations, because the strength of all deviations is here controlled by the
leading PN order term (as expected from most modified theories of gravity). If
an MCMC chain keeps the leading PN order deformation small, then the overall
modification at each PN order
($\bar{\Delta}=\delta\bar{\psi}_{i}\times\bar{\gamma}$) can remain small even
if the chain visits large values of $\delta\bar{\psi}_{i}$. This opens up more
of the prior volume for chain exploration, which can thus be preferred more by
the posterior distribution. In other words, there are more “states” for the
model to take with small leading PN order deformation than with large leading
PN order deformations, putting “pressure” on the leading coefficient to remain
small.

But if the higher PN order terms were allowed to explore infinitely large
values, then this would put infinite “pressure” on the leading PN order
coefficient, forcing its posterior to be artificially tight around zero.
Whether this happens or not depends on the prior one chooses on the higher PN
order parameters. Naively, one may think that the most conservative prior is
one that is infinitely wide and flat (this would correspond to our PN-based
priors but evaluated at an infinite orbital separation.) Choosing such a flat,
uniform prior, however, would allow the magnitude of the higher PN order
deformation parameters to increase indefinitely. In turn, this would put
infinite “pressure” on the leading PN order coefficient and push its posterior
to zero. The result would be an overly confident or overly aggressive
constraint that is generated not by the information contained in the data, but
by the prior choice. Our PN-based prior prevents this from happening. By
ensuring the higher order deformations are not infinitely large (as otherwise
the PN series would break down), the prior volume remains a reasonable size,
and the “pressure” on the leading PN order deformation, sourced by this prior
volume, is therefore kept to a reasonable amount. This then means that the
constraint on the leading PN order parameter is governed by the information
contained in the data (and the mathematical requirement that the PN
approximation be valid in the inspiral for the modified theory), and not by
the imposition of an overly restrictive prior.

Given how important the PN prior is to prevent overly aggressively constraints
on the leading PN order term, one may wonder whether these are robust to
variations of the details associated with this prior. The main quantity we can
vary to change the PN prior is the choice of velocity $v_{\rm eval}$ (or
orbital separation $r_{12,eval}$) at which the prior is evaluated. Figure 4
shows the constraints on $\bar{\gamma}$ and the strengthening factor as a
function of the number of higher PN order terms kept in the modified sector,
but this time obtained with PN-based priors evaluated at 3 different orbital
radii. Observe that in all cases the constraints improve as you add higher PN
order terms in the modified sector. Therefore, the first main conclusion of
our paper (i.e. that single-parameter constraints are conservative) is robust
to modifications in the PN-based prior.

Figure 4 also shows that the degree of importance of the higher PN order terms
does depend on the PN-based prior. Indeed, the strengthening factor is larger,
the larger the value of $r_{12,\rm eval}$ that one chooses for the PN-based
prior. This is consistent with our explanation above that as $r_{12,\rm
eval}\to\infty$, then the prior on the higher PN order terms becomes flat and
infinite, and, thus, $\Delta\bar{\gamma}\to 0$ because of the infinite
“pressure” created by the higher PN order terms. The choice of $r_{12,\rm
eval}$ should then be the smallest value of the orbital separation (or the
orbital velocity) at which one expects the PN approximation to still be valid
in the inspiral of compact binaries in the modified theory. Without specifying
a particular modified theory, one can therefore not choose $r_{12,\rm eval}$
precisely. This is why we used the properties of the PN series in GR to set
$r_{12,\rm eval}=100m$ in this paper, since we are sure that for such large
values of $r_{12,\rm eval}$, the PN series in GR is still a good
approximation. Had we chosen a smaller value, our conclusions about the
importance of the higher PN order terms would have been even stronger (i.e. we
would have concluded the higher PN order terms are even less important than
stated so far).

Figure 4:  Marginalized $1\sigma$ constraint on $\bar{\gamma}$ for a
modification that first enters at 0PN order as a function of number of PN
terms kept in the modified sector, using 3 different choices of $r_{12,\rm
eval}$ for the PN prior. In all cases, we here focus on the light source and
the 2g detector network. Observe that in all cases the constraint on
$\bar{\gamma}$ becomes stronger the more PN order terms are added. Observe
also that the larger we choose $r_{12,\rm eval}$ to be, the stronger the
“pressure” on $\bar{\gamma}$ and thus the stronger the constraint.

Our conclusions thus far are the following. First, single-parameter tests of
GR are conservative, and would become stronger if higher PN order terms are
included in the modified sector. Second, the improvement of these higher PN
order terms is relatively mild, with enhancements of at most 1–2 orders of
magnitude, depending on the inspiral signal observed. Given these conclusions,
we then infer that current constraints on modified theories of gravity are
robust to uncertainties in waveform modeling related to unknown, higher PN
order corrections. But are these conclusions robust also to an improvement in
the detector’s sensitivity, or alternatively, to an increase in the signal-to-
noise ratio of the signals detected? Motivated by this, we repeated the
analysis described above for the proxy of a 3g network described in Sec. IV,
and we arrived at very similar conclusions, as shown in Fig. 5 and Table 3.
Observe that the increased SNR of the injections does not change the trends we
have described above for the less sensitive 2g networks. While the bounds
become stronger because of the more sensitive detectors, the inferences made
with single-parameter models are still robust to uncertainties in the higher
PN order terms.

Figure 5:  Same as Fig. 3, but for a 3g network, using a “heavy” source (left)
and a “light” source (right). Observe that the trends found with the 2g
detector networks continue when considering 3g detector networks.

Terms/LO | -1 | 0 | 1 | 1.5 | 2
---|---|---|---|---|---
1 | $1.4\times 10^{-7}$ | $1.9\times 10^{-4}$ | $7.8\times 10^{-4}$ | $6.0\times 10^{-3}$ | $7.8\times 10^{-2}$
2 | $1.3\times 10^{-7}$ | $7.1\times 10^{-5}$ | $7.9\times 10^{-4}$ | $3.3\times 10^{-3}$ | $3.7\times 10^{-2}$
3 | $1.6\times 10^{-7}$ | $3.6\times 10^{-5}$ | $7.1\times 10^{-4}$ | $2.7\times 10^{-3}$ |
4 | $1.5\times 10^{-7}$ | $3.0\times 10^{-5}$ | $2.5\times 10^{-4}$ |  |
5 | $1.6\times 10^{-7}$ | $1.3\times 10^{-5}$ |  |  |
6 | $1.6\times 10^{-7}$ | $3.2\times 10^{-6}$ |  |  |

Terms/LO | -1 | 0 | 1 | 1.5 | 2
---|---|---|---|---|---
1 | $5.5\times 10^{-9}$ | $5.8\times 10^{-5}$ | $9.0\times 10^{-4}$ | $7.7\times 10^{-3}$ | $8.2\times 10^{-2}$
2 | $5.3\times 10^{-9}$ | $2.0\times 10^{-5}$ | $7.4\times 10^{-4}$ | $7.0\times 10^{-3}$ | $5.2\times 10^{-2}$
3 | $5.0\times 10^{-9}$ | $1.5\times 10^{-5}$ | $5.2\times 10^{-4}$ | $3.6\times 10^{-3}$ |
4 | $5.4\times 10^{-9}$ | $1.1\times 10^{-5}$ | $3.6\times 10^{-4}$ |  |
5 | $5.7\times 10^{-9}$ | $7.6\times 10^{-6}$ |  |  |
6 | $5.5\times 10^{-9}$ | $3.6\times 10^{-6}$ |  |  |

Table 3:  Same as Table 2, but for a 3g network, using the “heavy” source
(top) and the “light” source (bottom). Once more, the trends found with the 2g
detector networks continue when considering a 3g detector network.

##  VI Alternative Parametrizations

One may wonder how robust the conclusions presented in Sec. III are to the
exact form of the parametrization we implemented in that section. Therefore,
before continuing on to a specific theory of gravity, let us briefly examine
two other reasonable parameterizations to determine how our conclusions are
affected. In what follows, we will re-analyze the light system with the 2g
network, exactly as defined earlier, using 6 deformations and starting with
the Newtonian-order term relative to GR, i.e. terms at [0,1,1.5,2,3,3.5]PN
orders relative to GR. The difference, however, will be in the exact form of
the parametrization of the deformations at each of those orders and how the
prior for those deformations is imposed.

Figure 6:  The one dimensional, marginalized prior on a set of deformations in
the three models considered in Sec. VI. The leading-order deformation enters
at Newtonian order relative to GR, and each of the subsequent deformations are
labeled as the $i/2$PN order relative to this leading order deformation.
Included in this figure is the original model described in Sec. III, shown as
a solid blue line. The first alternative model uses the $\delta\bar{\psi}_{i}$
parametrization described in Eq. (III.1) but with simple, uniform ranges with
fixed (increasingly larger with higher PN order) boundaries, shown as a dotted
orange line. The second alternative model is described by a series of
$\bar{\Delta}_{i}$ parameters with uniform distributions with fixed
(increasingly larger with higher PN order) boundaries, defined in Eq. (III.2),
shown as a dashed green line. Note that the range of the prior on the
$\delta\bar{\psi}_{i}$ parameters in the second alternative parametrization
extends beyond the frame of the figure, but the range was restricted for
visual purposes.

The first model uses the same functional form as the majority of our work
here, defining the deformations in the phase as a series of coefficients
exactly as shown in Eq. (III.1), but with a different prior. We modify the
prior to have simple, fixed boundaries, as opposed to the more complicated
prior used in the main body of this work (previously, we enforced the
convergence criteria at every point in parameter space through Eq. (15) and
Eq. (16)). In particular, the prior for the new model will have fixed
boundaries on $\bar{\gamma}$ and the $\delta\bar{\psi}_{i}$’s with
progressively larger ranges, namely

$\displaystyle|\bar{\gamma}|<10^{-1}\,,$
$\displaystyle|\delta\bar{\psi}_{2}|<10^{2}\,,$
$\displaystyle|\delta\bar{\psi}_{3}|<10^{3}\,,$
$\displaystyle|\delta\bar{\psi}_{4}|<10^{4}\,,$
$\displaystyle|\delta\bar{\psi}_{6}|<10^{6}\,,$
$\displaystyle|\delta\bar{\psi}_{7}|<10^{7}\,.$ (22)

This prior ensures that the convergence criteria is _generally_ satisfied but
not guaranteed, while removing some of the complicated structure from our
original prior (shown in Fig 2). By reanalyzing our synthetic data with this
new model, we will quantify how strong of an impact our enforcement of a
strict convergence criteria had on our results above.

We note in passing that the prior range on $\bar{\gamma}$ never played much of
a role in our analysis, as the prior bounds on the higher-order parameters,
$\delta\bar{\psi}_{i}$’s, were independent of $\bar{\gamma}$. As we will see
below, that will no longer be the case, and we therefore updated our prior
bound on $\bar{\gamma}$ to a much more reasonable range of
$|\bar{\gamma}|<10^{-1}$ (informed by our original experiments). We have
verified that we obtain the same results with this new prior range and with
the original parametrization.

The second parametrization we investigate is written in the form of Eq.
(III.2), where our model is parametrized by $\bar{\gamma}$ and the series of
$\bar{\Delta}_{i}$’s. In other words, we are moving from parametrizing our
deformations as a relative series of terms to working with a model described
by the _absolute_ deformations. With this parametrization, we enforce similar
priors as our first alternative parametrizations, with fixed boundaries of
successively larger sizes, namely

$\displaystyle|\bar{\gamma}|<10^{-1}\,,$
$\displaystyle|\bar{\Delta}_{2}|<10^{1}\,,$
$\displaystyle|\bar{\Delta}_{3}|<10^{2}\,,$
$\displaystyle|\bar{\Delta}_{4}|<10^{3}\,,$
$\displaystyle|\bar{\Delta}_{6}|<10^{5}\,,$
$\displaystyle|\bar{\Delta}_{7}|<10^{6}\,.$ (23)

With this model, we will investigate the impact that our choice of
parameterization had on our conclusions of the previous section.

As the prior distribution is the major reason for testing these alternative
parametrizations, we show all three prior distributions on $\bar{\gamma}$,
$\delta\bar{\psi}_{i}$, and $\bar{\Delta}_{i}$ in Fig. 6. From this figure, we
can see that the first alternative parametrization produces the strongest
prior on the absolute deformations, $\bar{\Delta}_{i}$, for low orders, but
the original prior still places the most prior weight at small values of
$\bar{\Delta}_{i}$ for large PN orders. Our first alternative parametrization
does alleviate the “pressure” on the $\delta\bar{\psi}_{i}$ parameters seen in
the original parametrization, which was pushing low PN deformations away from
zero (and away from GR) and pushing high PN deformations towards zero. The
second alternative parametrization has a uniform distribution in the absolute
deformations, or the $\bar{\Delta}_{i}$ parameters, giving a variation of an
uninformative prior. From this figure, we would expect comparable results from
the original parametrization and the first alternative model, as the prior
distributions are relatively similar. The second alternative parametrization,
uniform in $\bar{\Delta}_{i}$, will most likely result in the most degradation
between the six-deformation and single-deformation model, as the prior is the
least informative and most similar to the current, state-of-the art
methodology.

We now use these two additional parametrizations to analyze the “light”
source, as seen by the 2g detector network, and compare the constraints on
$\bar{\gamma}$ to the posterior coming from a model that has a single
deformation at Newtonian order and to the constraint coming from the original
parametrization with six deformations. As all three models (the original model
and the two alternative parametrizations) reduce to the same form as the
number of deformations is taken to one, this will provide a consistent metric
for comparison. The results for those analyses are shown in Fig. 7, and as you
can observe, they are consistent with the expectations presented above. The
two models described by a relative series of deformations at higher orders
(the original model and the first alternative model) produce comparable
results. The second alternative model (uniform in $\bar{\Delta}_{i}$) produces
a constraint that is mildly weaker than the original, single deformation
model. We therefore conclude that the conclusions presented in Sec. III are
robust to our parameterization and choice of priors.

Figure 7:  The final, marginalized posterior distributions on $\bar{\gamma}$
coming from the three models discussed in Sec. VI using the synthetic data
from the “light” source as observed by the 2g detector network. The top panel
shows a larger range of $\bar{\gamma}$ and the lower panel zooms-in the range
of $\bar{\gamma}$ to present details in the two narrower distributions. The
probabilities on the y-axis are normalized to the shown range of
$\bar{\gamma}$ in both panels. The distribution coming from the model using a
single deformation parameter at Newtonian order is shown as a solid blue line.
The distribution coming from the model described by the first alternative
parametrization with six deformations (with a prior uniform in
$\delta\bar{\psi}_{i}$) is shown as the dotted orange line. The distribution
coming from the model described by the second alternative parametrization with
six deformations(with a prior uniform in $\bar{\Delta}_{i}$) is shown as the
dashed green line. The distribution coming from the model described by the
original parametrization with six deformations (with a prior that strictly
enforces our notion of convergence) is shown as the dotted-dashed red line.
Note that the first alternative parametrization and the original
parametrization are almost totally overlapping in the top panel.

##  VII A Concrete Example: ssGB theory

The methodology we have proposed here seems to accomplish the purpose it was
designed for: allowing for realistic uncertainty in our limited modeling while
illustrating the current robustness of the bounds placed on modified theories
of gravity. A critical question, however, is the accuracy to which it actually
relates to known, interesting theories of gravity actively being studied. One
particular theory currently of interest is scalar Gauss Bonnet, or sGB, which
is inspired by low-energy limits of higher energy theories [67, 68, 69]. sGB
contains a scalar field, $\phi$, that couples to a curvature invariant called
the Gauss Bonnet invariant, $\mathcal{G}$. The Gauss Bonnet invariant is
defined by
$\mathcal{G}=R^{2}-4R_{\mu\nu}R^{\mu\nu}+R_{\mu\nu\rho\sigma}R^{\mu\nu\rho\sigma}$,
where $R_{\mu\nu\rho\sigma}$, $R_{\mu\nu}$, and $R$ are the Riemann tensor,
Ricci tensor, and Ricci scalar, respectively. The full action can be written
as [38]

$S=\frac{1}{16\pi}\int
d^{4}x\sqrt{-g}\left[R-2\partial^{\mu}\phi\partial_{\mu}\phi+\alpha
f(\phi)\mathcal{G}\right]\,,$ (24)

where $g$ is the determinant of the metric, $\alpha$ is the (dimensionful)
coupling constant of the theory, and $f(\phi)$ is a coupling function. The
theory reduces to GR minimally coupled to the scalar field in the $\alpha\to
0$. Therefore, constraints on this theory typically relate to the magnitude of
$\alpha$. Certain choices of the coupling function $f(\phi)$ result in
different flavors of sGB. In particular, an interesting choice motivated by
string theory is Einstein-dilaton-Gauss-Bonnet, where $f(\phi)=e^{2\phi}/4$
[67]. Another option, which will be the focus of this paper, is shift-
symmetric sGB or ssGB, obtained by making the choice of $f(\phi)=2\phi$ [70,
71, 72], which results from a small-field expansion of Einstein-dilaton-Gauss-
Bonnet gravity [71].

The GW phase in ssGB was recently calculated beyond the leading PN order
modification (at -1PN order, relative to Newtonian order in GR) in [38],
giving us an example theory to explore. In their work, they classified the
inspiral into two distinct regimes: that which is dominated by the emission of
dipolar radiation by the scalar field (which is absent in GR), and that which
is dominated by quadrupolar radiation. The boundary between dipolar and
quadrupolar radiation is typically well before the binary reaches a separation
of $100m$, so the regime of interest to ground-based GW detectors is the
quadrupolar driven one. To evaluate the convergence properties of the series,
we take the expression for the GW phase derived in Eq. 93–96 of [38] and
isolate terms involving deviations from GR, i.e. terms that involve $\alpha$,
the coupling constant in ssGB. This gives the following form for the phase
modification

$\Psi_{\rm GW}-\Psi_{\rm{GW},\rm{GR}}=\delta\psi_{\rm LO,\rm
ssGB}v^{-7}+\delta\psi_{\rm NLO,\rm ssGB}v^{-5}+\delta\psi_{\rm NNLO,\rm
ssGB}v^{-3}\,,$ (25)

where the coefficients take the form

$\displaystyle\delta\psi_{\rm LO,\rm ssGB}$
$\displaystyle=-\frac{5}{7168}\zeta\frac{\left(4\eta-1\right)}{\eta^{5}}\,,$
(26) $\displaystyle\delta\psi_{\rm NLO,\rm ssGB}$
$\displaystyle=-\frac{5}{688128}\zeta\frac{\left(685-3916\eta+2016\eta^{2}\right)}{\eta^{5}}\,,$
(27) $\displaystyle\delta\psi_{\rm NNLO,\rm ssGB}$
$\displaystyle=\frac{5}{387072}\zeta\frac{\left(1-2\eta\right)^{2}\left(995+952\eta\right)}{\eta^{5}}\,,$
(28)

where $\zeta\equiv{\alpha^{2}}/{m^{4}}$ is the dimensionless coupling constant
of the theory222The expressions presented here are not valid in the $\eta\ll
1$ regime, because then the curvature of the small black hole becomes very
large, and the effective field theory treatment used to derive these
expressions breaks down.. We have verified that this expression reduce exactly
to the leading PN order results first obtained in [73, 12] in the non-spinning
limit333Note that the results of [38] are formally complete only to 1PN order,
implying that the $v^{-3}$ coefficient has a 2PN correction coming from the
dipolar term that has not yet been computed.. Observe that the phase
deformation has the exact structure that we anticipated in Eq. (III.1).
Namely, the above phase deformation can be rewritten as

$\displaystyle\Psi_{\rm GW}-\Psi_{\rm GW,\rm GR}$
$\displaystyle=\zeta\,\delta\psi_{\rm LO,ssGB}(\eta)\,v^{-7}$
$\displaystyle\times\left[1+\delta\psi^{\prime}_{\rm NLO,\rm
ssGB}(\eta)\,v^{2}+\delta\psi^{\prime}_{\rm NNLO,\rm
ssGB}(\eta)\,v^{4}\right]\,,$ (29)

with the following re-definitions of the coefficients

$\displaystyle\delta\psi^{\prime}_{\rm NLO,\rm ssGB}$
$\displaystyle=\frac{\delta\psi_{\rm NLO,\rm ssGB}}{\delta\psi_{\rm
LO,ssGB}}=\frac{685-3916\eta+2016\eta^{2}}{96\left(4\eta-1\right)}\,$ (30)
$\displaystyle\delta\psi^{\prime}_{\rm NNLO,\rm ssGB}$
$\displaystyle=\frac{\delta\psi_{\rm NNLO,\rm ssGB}}{\delta\psi_{\rm
LO,ssGB}}=-\frac{\left(1-2\eta\right)^{2}(995+952\eta)}{54\left(4\eta-1\right)}\,.$
(31)

The phase modification has therefore been written as a series in $v$, all
proportional to the coupling constant $\alpha$ appearing in the coefficient of
the overall controlling factor of the series $\delta\psi_{\rm LO,ssGB}$.

With this in hand, we can now compare the above results to the general
framework we developed in the previous section when studying the convergence
properties of the PN series in modified theories. To do so, we need to look at
the ratios of these coefficients, which we show below. For the LO and NLO
terms, we arrive at the following expression for their ratio

$\frac{\delta\psi_{\rm LO,\rm ssGB}}{\delta\psi_{\rm NLO,\rm ssGB}v_{\rm
eval}^{2}}=\frac{96(4\eta-1)}{\left[685+4\eta\left(-979+504\eta\right)\right]v_{\rm
eval}^{2}}\,.$ (32)

For the NLO and NNLO terms, we evaluate their ratio to

$\frac{\delta\psi_{\rm NLO,\rm ssGB}}{\delta\psi_{\rm NNLO,\rm ssGB}v_{\rm
eval}^{2}}=\frac{9\left[685+4\eta\left(-979+504\eta\right)\right]}{16\left(1-2\eta\right)^{2}\left(995+952\eta\right)v_{\rm
eval}^{2}}\,.$ (33)

As we are only considering non-spinning binaries, these ratios are only a
function $\eta$, so we can easily plot them to determine their convergence
properties. The absolute deviations (divided by $\zeta$), as shown in Eqs.
(26), (27), and (28) and including appropriate values of $v_{\rm eval}$, are
shown in the left panel of Fig. 8, while the ratios of these coefficients, as
defined in Eqs. (32) and (33), are shown in the right panel.

Figure 8: Absolute value of the terms in the phase deformation at different PN
order (left) and ratio of terms (right) in ssGB as a function of symmetric
mass ratio, evaluated at $v_{\rm eval}=0.1$. Observe that, for this choice of
$v_{\rm eval}$, the leading PN order term is larger than the next-to-leading
order one, which in turn is larger than the next-to-next-to-leading order,
except at $\eta\sim 0.19$, where the 1PN term vanishes identically. Observe
also that the ratios are all larger than unity, except again for a specific
value of $\eta$ at which the 1PN term vanishes.

Figure 8 illustrates that the PN expansion of ssGB theory conforms, on
average, to the criteria we have outlined in Sec. III. This is most easily
evidenced by the right panel of this figure, which shows that the ratio of
coefficients is always above unity, except for a specific value of $\eta$ at
which the 1PN term $\delta\psi_{\rm NLO,\rm ssGB}$ vanishes. Since this
happens at a single point, averaging over all values of $\eta$, it is clear
that the convergence criteria is satisfied.

But are our conclusions from Sec. V still valid in this specific theory? To
answer this question, we completed another series of Bayesian studies, but now
within the context of ssGB theory. We used the “light” injection from Table 1,
and recovered with a IMRPhenomD model with one of the following phase
deformations appended to its inspiral phase:

  * •

Model 1: only the leading PN order deformation in ssGB theory (defined in Eq.
(26)),

  * •

Model 2: all the phase deformation terms in ssGB theory (defined in Eq. (27)
and Eq. (28))

  * •

Model 3: a single ppE deformation term (as in Eq. (III.1) but with
$\delta\bar{\psi}_{i>0}=0$),

  * •

Model 4: 6-parameter ppE deformations (as in Eq. (III.1) but with
$\delta\bar{\psi}_{i>6}=0$).

For models 3 and 4, after carrying out a Bayesian parameter estimation study,
we mapped the constraints on $\bar{\gamma}$ to constraints on $\alpha$ through

$\bar{\gamma}=\delta\psi_{\rm LO,ssGB},\qquad b=-7\,,$ (34)

to enable comparisons with the results using models 1 and 2. For reference,
the posteriors for $\bar{\gamma}$ for models 3 and 4 are shown in Fig. 9. This
figure simply reiterates in more detail the conclusions from Sec. V (Fig. 3
and Table 2). Namely, the difference between using a single deformation and
six deformations within the model established in Sec. III is a mild
enhancement of constraints.

Figure 9:  For reference, we show the two distributions on $\bar{\gamma}$
coming from the two generic deformation models (models 3 and 4). This figure
shows that the distributions on the generic parameter $\bar{\gamma}$ for these
two models are comparable, with a slight enhancement when using the 6
deformation variation. This same information is shown graphically in Fig. 3
and tabulated in Table 2. When compared with Fig. 10, we see that this
relation persists. Even once the posteriors are transformed to theory-specific
constants, constraints are not meaningfully changed when using these two
models. Furthermore, note that the posterior on $\bar{\gamma}$ is fully
consistent with GR (the model that truly describes the injected data), which
is further evidence that the deviation away from GR in Fig 10 for the generic
models is a artifact of the transformation.

Marginalized posterior distributions of $\sqrt{\alpha}$ obtained using these
four models are shown in Fig. 10. While the transformation leaves certain
artifacts in the posterior distributions because of singularities in the
transformation (as discussed in great detail in [35] and [74]), the upper
limits on $\sqrt{\alpha}$ obtained with the four models are consistent with
each other. The $90\%$ confidence upper limit on $\sqrt{\alpha}$ with a single
deformation to the phase (the red curve in Fig. 10) is $5.0$ km, while that
obtained using the full, three term modification to the GW phase (the green
curve in Fig. 10) is $1.5$ km, leading to a strengthening factor of about 3.
The additional information incorporated into the waveform through the higher
order deformation only serve to improve the constraint on $\sqrt{\alpha}$.
This comparison provides further evidence that constraints derived from
leading PN order deformations are robust to future work on deriving higher PN
order corrections.

Figure 10:  Marginalized posterior distributions on $\sqrt{\alpha}$ for a GR
injection extracted with the four models described in the text. We here
considered the “light” source, defined by Table 1, and a 2g detector network.
When carrying out the Bayesian studies with the multi-parameter ppE model, we
employed the PN prior with $r_{12,\rm eval}=100$, as done in the rest of this
paper. The constraints on $\bar{\gamma}$ obtained with models 3 and 4 were
mapped to constraints on $\sqrt{\alpha}$ to enable comparisons with the
results obtained with models 1 and 2. Observe that clearly the posterior
distributions are all consistent with each other (modulo singularities in the
transformation at $\alpha=0$, discussed in detail in [74, 35]). This shows
clearly that in ssGB theory, leading PN order constraints are conservative and
sufficient to place bounds on the theory.

The bias in the distribution on $\sqrt{\alpha}$ coming from the generic
parametrizations is a known issue [74, 35], and is not of much concern in the
present context. The issue fundamentally lies with the Jacobian of the
transformation between the two parametrizations, and causes the derived prior
distribution on $\sqrt{\alpha}$ to go to zero in the GR limit when mapping the
generic parametrization to ssGB. When transforming distributions between two
different basis, one must account for the Jacobian of the transformation as

$p(\bm{\theta}_{1})=p(\bm{\theta}_{2})\frac{d\bm{\theta}_{2}}{d\bm{\theta}_{1}}\,,$
(35)

In the present case, we are using one parametrization
$\bm{\theta}_{1}\equiv[\Xi\cup\sqrt{\alpha}]$ in ssGB and one parametrization
$\bm{\theta}_{2}\equiv[\Xi\cup\bar{\gamma}\cup\sum_{i}\delta\bar{\psi}_{i}]$
in the generic framework. From Eq. (35), we know we need the Jacobian,
$d\bm{\theta}_{1}/d\bm{\theta}_{2}$, to transform from parametrization 2 to
parametrization 1. Here lies the issue, as the component of the Jacobian
related to $\sqrt{\alpha}$ and $\bar{\gamma}$ is given as

$\frac{d\bar{\gamma}}{d\sqrt{\alpha}}\propto(\alpha)^{3/2}\,,$ (36)

using the relation in Eq. (26). As this expression goes to 0 in the
$\sqrt{\alpha}\rightarrow 0$ limit, the prior will have zero weight for the GR
limit. Now the discrepancy is clear: the lack of agreement with the posterior
on $\sqrt{\alpha}$ coming from the generic parametrization is not an
indication that GR is lacking, but instead an inherent flaw of the
parametrization. While this is an indication of a failing of the
parametrization, the fact that the prior seems to be widening the constraint
on $\sqrt{\alpha}$ instead of artificially shrinking it indicates that any
constraints derived from this method related to _upper limits on
$\sqrt{\alpha}$_ are actually _conservative_.

Graphically, this is illustrated in Fig. 11 where the posteriors on
$\sqrt{\alpha}$ coming from mapping constraints from the generic models are
plotted along side samples from the full, derived prior on $\sqrt{\alpha}$.
The samples from the prior were drawn from a uniform distribution for
$\bar{\gamma}$ and the usual priors for the source parameters, then mapped to
$\sqrt{\alpha}$ in the same way as the analysis using models 3 and 4. The
original posterior and prior are completely consistent with GR, as shown in
this figure, so the lack of support in the GR limit for $\sqrt{\alpha}$ must
come from the transformation itself.

Figure 11:  Marginalized posterior distributions on $\sqrt{\alpha}$ for a GR
injection extracted with the two generic models described in the text. One of
the models uses a single deformation while the second model uses a series of
six deformations. We here considered the “light” source, defined by Table 1,
and a 2g detector network. These two data sets are models 3 and 4, described
in the text and shown in Fig. 10. For comparison, we also show the derived
prior on $\sqrt{\alpha}$, calculated by taking a uniform prior on
$\bar{\gamma}$ and the usual priors on source parameters, then mapping it to
$\sqrt{\alpha}$ with the same prescription as was used for transforming the
two generic model constraints. This figure illustrates the interesting
behavior around $\sqrt{\alpha}\rightarrow 0$ (the GR limit). The fact that the
derived prior on $\sqrt{\alpha}$, as transformed from a uniform prior on
$\bar{\gamma}$, disallows $\sqrt{\alpha}=0$ explains the bias in the
posteriors for models 3 and 4. The issue is related to the Jacobian of the
transformation, discussed in the text, and is not of serious concern in
analysis such as these.

##  VIII Conclusions

We have here studied whether the inclusion of higher PN order terms in the
modified gravity deformations to the GW phase of inspiraling binaries affects
the constraints one can place on these theories. First, we focus on a
particular implementation of multi-parameter phase deformation, in which an
overall controlling factor contains all of the coupling-constant information
of the modified theory. Then, we develop a novel, PN-based prior to impose on
parametric deviations, which ensures that the terms added obey certain
convergence criteria, as they must if they derive from a PN expansion in the
inspiral phase (even in a modified theory).

Our analysis shows that the constraints placed on modified theories with
single-parameter ppE waveforms are robust and reliable. More specifically, the
inclusion of higher PN order terms in the inspiral phase do not weaken the
constraints we can place with single-parameter ppE models. In fact, the
inclusion of these terms actually improves the constraints on coupling
constants of modified theories, and typically the strengthening of the bound
is mild. We verified that these conclusions with an array of Bayesian
parameter estimation studies, in which we injected synthetic GR signals and
extracted with a variety of single- and multi-parameter ppE models. We further
considered a specific theory, ssGB gravity, to exemplify our findings with a
concrete set of deformations to the phase. As expected from our generic
analysis, the inclusion of higher PN order terms in the ssGB inspiral phase
does not weaken bounds obtained with leading order waveforms. Instead, the
higher PN order terms improve the bounds on the ssGB coupling parameter, but
only by about a factor of 3. These results are consistent with a very recent
analysis of ssGB theory with the same higher PN order model we use here on
real aLIGO/Virgo data [75].

One can compare our methodology to other techniques to measure multiple phase
deformations simultaneously. One such methodology is through the use of a
principal component analysis (PCA) on tests of GR [76, 77, 78]. In this
method, the phase deformations are combined through certain linear
transformations, so that covariances are minimized. To do this, samples are
first drawn from the posterior distribution with a waveform model that
includes multiple ppE phase deformations, and then this posterior distribution
is decomposed into its eigenvectors. Constraints are placed on linear
combinations of deformations that produce the tightest constraints (the
eigenvectors with the largest eigenvalues). Effectively, this just equates to
finding a parametrization that optimizes the constraint you can place on these
deformations, but it comes with serious drawbacks related to tests of
fundamental physics. All information is lost about the physical meaning of
these deformations, as this basis has no connection to calculations performed
in modified gravity. This can still be an effective consistency test of GR (in
the same category of “residual tests” discussed in Sec. I), but it does not
provide information about fundamental physics beyond this specific type of
test. Our approach, instead, attempts to stay close to theoretical physics,
using physical insight to concentrate our attention on the “relevant”
deformation, the controlling factor $\bar{\gamma}$ that encodes the coupling
constants of the theory, at the expense of losing less relevant information
about the higher order corrections (which are just functions of the system
parameters). With this, we still retain the mapping between theories of
modified gravity and $\bar{\gamma}$, allowing constraints in the latter to
lead to constraints on the coupling constants of modified theories.

These results encourage the use of leading PN order deformations in the
inspiral phase to constrain theories of gravity beyond GR. The results will
certainly improve with the inclusion of higher PN order deformations, but the
bounds will not degrade through their inclusion. This does, however, come with
the caveat that we have only considered theories of gravity that allow for
well-behaved, series solutions to the field equations, both in the PN
expansion and the expansion in the coupling. There are theories that do not
conform to these criteria, such as theories of gravity exhibiting spontaneous
scalarization [46, 79, 80, 42, 43] or other abrupt, discontinuous
transformations. Future work could focus on those theories to attempt to
develop a generic framework that can also be applied to them.

###### Acknowledgements.

The authors would like to thank Walter del Pozzo and Emanuele Berti for useful
comments on the first draft of this manuscript. We would also like to
especially thank Jonathan Gair for detailed comments on the first draft, which
led us to add Sec. VI to the final version. This work was partially supported
by the Center for AstroPhysical Surveys (CAPS) at the National Center for
Supercomputing Applications (NCSA), University of Illinois at Urbana-
Champaign. S. P. acknowledges support from the Illinois Center for Advanced
Studies of the Universe (ICASU)/CAPS Graduate Fellowship. N. Y. acknowledges
support from the Simons Foundation through Grant No. 896696. This work made
use of the Illinois Campus Cluster, a computing resource that is operated by
the Illinois Campus Cluster Program (ICCP) in conjunction with the National
Center for Supercomputing Applications (NCSA) and which is supported by funds
from the University of Illinois at Urbana- Champaign.

## References

  * Will [2014] C. M. Will, The Confrontation between General Relativity and Experiment, Living Rev. Rel. 17, 4 (2014), arXiv:1403.7377 [gr-qc] .
  * Berti  _et al._ [2015] E. Berti _et al._ , Testing General Relativity with Present and Future Astrophysical Observations, Class. Quant. Grav. 32, 243001 (2015), arXiv:1501.07274 [gr-qc] .
  * Shomer [2007] A. Shomer, A Pedagogical explanation for the non-renormalizability of gravity, arXiv:0709.3555 [hep-th] (2007).
  * Senovilla and Garfinkle [2015] J. M. M. Senovilla and D. Garfinkle, The 1965 Penrose singularity theorem, Class. Quant. Grav. 32, 124008 (2015), arXiv:1410.5226 [gr-qc] .
  * Penrose [1965] R. Penrose, Gravitational collapse and space-time singularities, Phys. Rev. Lett. 14, 57 (1965).
  * Riess  _et al._ [1998] A. G. Riess _et al._ (Supernova Search Team), Observational evidence from supernovae for an accelerating universe and a cosmological constant, Astron. J. 116, 1009 (1998), arXiv:astro-ph/9805201 .
  * Perlmutter  _et al._ [1999] S. Perlmutter _et al._ (Supernova Cosmology Project), Measurements of $\Omega$ and $\Lambda$ from 42 high redshift supernovae, Astrophys. J. 517, 565 (1999), arXiv:astro-ph/9812133 .
  * Sofue and Rubin [2001] Y. Sofue and V. Rubin, Rotation curves of spiral galaxies, Ann. Rev. Astron. Astrophys. 39, 137 (2001), arXiv:astro-ph/0010594 .
  * Bertone and Hooper [2018] G. Bertone and D. Hooper, History of dark matter, Rev. Mod. Phys. 90, 045002 (2018), arXiv:1605.04909 [astro-ph.CO] .
  * Canetti  _et al._ [2012] L. Canetti, M. Drewes, and M. Shaposhnikov, Matter and Antimatter in the Universe, New J. Phys. 14, 095012 (2012), arXiv:1204.4186 [hep-ph] .
  * Carson and Yagi [2020] Z. Carson and K. Yagi,  _Testing General Relativity with Gravitational Waves_ (2020) arXiv:2011.02938 [gr-qc] .
  * Yunes  _et al._ [2016] N. Yunes, K. Yagi, and F. Pretorius, Theoretical Physics Implications of the Binary Black-Hole Mergers GW150914 and GW151226, Phys. Rev. D94, 084002 (2016), arXiv:1603.08955 [gr-qc] .
  * Abbott  _et al._ [2016] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Tests of general relativity with GW150914, Phys. Rev. Lett. 116, 221101 (2016), [Erratum: Phys.Rev.Lett. 121, 129902 (2018)], arXiv:1602.03841 [gr-qc] .
  * Cornish and Littenberg [2015] N. J. Cornish and T. B. Littenberg, BayesWave: Bayesian Inference for Gravitational Wave Bursts and Instrument Glitches, Class. Quant. Grav. 32, 135012 (2015), arXiv:1410.3835 [gr-qc] .
  * Cornish  _et al._ [2021] N. J. Cornish, T. B. Littenberg, B. Bécsy, K. Chatziioannou, J. A. Clark, S. Ghonge, and M. Millhouse, BayesWave analysis pipeline in the era of gravitational wave observations, Phys. Rev. D 103, 044006 (2021), arXiv:2011.09494 [gr-qc] .
  * Chatziioannou  _et al._ [2021] K. Chatziioannou, M. Isi, C.-J. Haster, and T. B. Littenberg, Morphology-independent test of the mixed polarization content of transient gravitational wave signals, Phys. Rev. D 104, 044005 (2021), arXiv:2105.01521 [gr-qc] .
  * Cornish  _et al._ [2011] N. Cornish, L. Sampson, N. Yunes, and F. Pretorius, Gravitational Wave Tests of General Relativity with the Parameterized Post-Einsteinian Framework, Phys. Rev. D 84, 062003 (2011), arXiv:1105.2088 [gr-qc] .
  * Vallisneri [2012] M. Vallisneri, Testing general relativity with gravitational waves: a reality check, Phys. Rev. D 86, 082001 (2012), arXiv:1207.4759 [gr-qc] .
  * Yunes and Pretorius [2009] N. Yunes and F. Pretorius, Fundamental Theoretical Bias in Gravitational Wave Astrophysics and the Parameterized Post-Einsteinian Framework, Phys. Rev. D 80, 122003 (2009), arXiv:0909.3328 [gr-qc] .
  * Sampson  _et al._ [2013] L. Sampson, N. Cornish, and N. Yunes, Gravitational Wave Tests of Strong Field General Relativity with Binary Inspirals: Realistic Injections and Optimal Model Selection, Phys. Rev. D 87, 102001 (2013), arXiv:1303.1185 [gr-qc] .
  * Chatziioannou  _et al._ [2012] K. Chatziioannou, N. Yunes, and N. Cornish, Model-Independent Test of General Relativity: An Extended post-Einsteinian Framework with Complete Polarization Content, Phys. Rev. D 86, 022004 (2012), [Erratum: Phys.Rev.D 95, 129901 (2017)], arXiv:1204.2585 [gr-qc] .
  * Mishra  _et al._ [2010] C. K. Mishra, K. G. Arun, B. R. Iyer, and B. S. Sathyaprakash, Parametrized tests of post-Newtonian theory using Advanced LIGO and Einstein Telescope, Phys. Rev. D 82, 064010 (2010), arXiv:1005.0304 [gr-qc] .
  * Arun  _et al._ [2006a] K. G. Arun, B. R. Iyer, M. S. S. Qusailah, and B. S. Sathyaprakash, Testing post-Newtonian theory with gravitational wave observations, Class. Quant. Grav. 23, L37 (2006a), arXiv:gr-qc/0604018 .
  * Abbott  _et al._ [2019a] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Tests of General Relativity with the Binary Black Hole Signals from the LIGO-Virgo Catalog GWTC-1, Phys. Rev. D 100, 104036 (2019a), arXiv:1903.04467 [gr-qc] .
  * Abbott  _et al._ [2021a] R. Abbott _et al._ (LIGO Scientific, Virgo), Tests of general relativity with binary black holes from the second LIGO-Virgo gravitational-wave transient catalog, Phys. Rev. D 103, 122002 (2021a), arXiv:2010.14529 [gr-qc] .
  * Abbott  _et al._ [2020] R. Abbott _et al._ (LIGO Scientific, Virgo), Tests of General Relativity with Binary Black Holes from the second LIGO-Virgo Gravitational-Wave Transient Catalog, arXiv:2010.14529 [gr-qc] (2020).
  * Abbott  _et al._ [2021b] R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), Tests of General Relativity with GWTC-3, arXiv:2112.06861 [gr-qc] (2021b).
  * Li  _et al._ [2012] T. G. F. Li, W. Del Pozzo, S. Vitale, C. Van Den Broeck, M. Agathos, J. Veitch, K. Grover, T. Sidery, R. Sturani, and A. Vecchio, Towards a generic test of the strong field dynamics of general relativity using compact binary coalescence, Phys. Rev. D 85, 082003 (2012), arXiv:1110.0530 [gr-qc] .
  * Agathos  _et al._ [2014] M. Agathos, W. Del Pozzo, T. G. F. Li, C. Van Den Broeck, J. Veitch, and S. Vitale, TIGER: A data analysis pipeline for testing the strong-field dynamics of general relativity with gravitational wave signals from coalescing compact binaries, Phys. Rev. D 89, 082001 (2014), arXiv:1311.0420 [gr-qc] .
  * Meidam  _et al._ [2018] J. Meidam _et al._ , Parametrized tests of the strong-field dynamics of general relativity using gravitational wave signals from coalescing binary black holes: Fast likelihood calculations and sensitivity of the method, Phys. Rev. D 97, 044033 (2018), arXiv:1712.08772 [gr-qc] .
  * Abbott  _et al._ [2019b] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Tests of General Relativity with GW170817, Phys. Rev. Lett. 123, 011102 (2019b), arXiv:1811.00364 [gr-qc] .
  * Tahura and Yagi [2018] S. Tahura and K. Yagi, Parameterized Post-Einsteinian Gravitational Waveforms in Various Modified Theories of Gravity, Phys. Rev. D 98, 084042 (2018), [Erratum: Phys.Rev.D 101, 109902 (2020)], arXiv:1809.00259 [gr-qc] .
  * Perkins  _et al._ [2021a] S. E. Perkins, N. Yunes, and E. Berti, Probing Fundamental Physics with Gravitational Waves: The Next Generation, Phys. Rev. D 103, 044024 (2021a), arXiv:2010.09010 [gr-qc] .
  * Arun  _et al._ [2006b] K. G. Arun, B. R. Iyer, M. S. S. Qusailah, and B. S. Sathyaprakash, Probing the non-linear structure of general relativity with black hole binaries, Phys. Rev. D 74, 024006 (2006b), arXiv:gr-qc/0604067 .
  * Perkins  _et al._ [2021b] S. E. Perkins, R. Nair, H. O. Silva, and N. Yunes, Improved gravitational-wave constraints on higher-order curvature theories of gravity, arXiv:2104.11189 [gr-qc] (2021b).
  * Blanchet [2014] L. Blanchet, Gravitational Radiation from Post-Newtonian Sources and Inspiralling Compact Binaries, Living Rev. Rel. 17, 2 (2014), arXiv:1310.1528 [gr-qc] .
  * Damour and Nagar [2016] T. Damour and A. Nagar, The effective-one-body approach to the general relativistic two body problem, in  _Astrophysical Black Holes_, edited by F. Haardt, V. Gorini, U. Moschella, A. Treves, and M. Colpi (Springer International Publishing, Cham, 2016) pp. 273–312.
  * Shiralilou  _et al._ [2021] B. Shiralilou, T. Hinderer, S. Nissanke, N. Ortiz, and H. Witek, Post-Newtonian Gravitational and Scalar Waves in Scalar-Gauss-Bonnet Gravity, arXiv:2105.13972 [gr-qc] (2021).
  * de Rham [2014] C. de Rham, Massive Gravity, Living Rev. Rel. 17, 7 (2014), arXiv:1401.4173 [hep-th] .
  * Damour and Esposito-Farese [1993] T. Damour and G. Esposito-Farese, Nonperturbative strong field effects in tensor - scalar theories of gravitation, Phys. Rev. Lett. 70, 2220 (1993).
  * Damour and Esposito-Farese [1996] T. Damour and G. Esposito-Farese, Tensor - scalar gravity and binary pulsar experiments, Phys. Rev. D 54, 1474 (1996), arXiv:gr-qc/9602056 .
  * Doneva and Yazadjiev [2018] D. D. Doneva and S. S. Yazadjiev, New Gauss-Bonnet Black Holes with Curvature-Induced Scalarization in Extended Scalar-Tensor Theories, Phys. Rev. Lett. 120, 131103 (2018), arXiv:1711.01187 [gr-qc] .
  * Silva  _et al._ [2018] H. O. Silva, J. Sakstein, L. Gualtieri, T. P. Sotiriou, and E. Berti, Spontaneous scalarization of black holes and compact stars from a Gauss-Bonnet coupling, Phys. Rev. Lett. 120, 131104 (2018), arXiv:1711.02080 [gr-qc] .
  * Ramazanoğlu [2017] F. M. Ramazanoğlu, Spontaneous growth of vector fields in gravity, Phys. Rev. D 96, 064009 (2017), arXiv:1706.01056 [gr-qc] .
  * Silva  _et al._ [2021a] H. O. Silva, A. Coates, F. M. Ramazanoğlu, and T. P. Sotiriou, The ghost of vector fields in compact stars, arXiv:2110.04594 [gr-qc] (2021a).
  * Sampson  _et al._ [2014] L. Sampson, N. Cornish, and N. Yunes, Mismodeling in gravitational-wave astronomy: The trouble with templates, Phys. Rev. D 89, 064037 (2014), arXiv:1311.4898 [gr-qc] .
  * Swendsen and Wang [1986] R. H. Swendsen and J.-S. Wang, Replica monte carlo simulation of spin-glasses, Phys. Rev. Lett. 57, 2607 (1986).
  * Earl and Deem [2005] D. J. Earl and M. W. Deem, Parallel tempering: theory, applications, and new perspectives., Physical chemistry chemical physics : PCCP 7 23, 3910 (2005).
  * Vousden  _et al._ [2016] W. D. Vousden, W. M. Farr, and I. Mandel, Dynamic temperature selection for parallel tempering in Markov chain Monte Carlo simulations, MNRAS 455, 1919 (2016), arXiv:1501.05823 [astro-ph.IM] .
  * Harris  _et al._ [2020] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, Array programming with NumPy, Nature 585, 357 (2020).
  * Husa  _et al._ [2016] S. Husa, S. Khan, M. Hannam, M. Pürrer, F. Ohme, X. Jiménez Forteza, and A. Bohé, Frequency-domain gravitational waves from nonprecessing black-hole binaries. I. New numerical waveforms and anatomy of the signal, Phys. Rev. D 93, 044006 (2016), arXiv:1508.07250 [gr-qc] .
  * Khan  _et al._ [2016] S. Khan, S. Husa, M. Hannam, F. Ohme, M. Pürrer, X. Jiménez Forteza, and A. Bohé, Frequency-domain gravitational waves from nonprecessing black-hole binaries. II. A phenomenological model for the advanced detector era, Phys. Rev. D 93, 044007 (2016), arXiv:1508.07253 [gr-qc] .
  * Chatziioannou  _et al._ [2017] K. Chatziioannou, A. Klein, N. Yunes, and N. Cornish, Constructing Gravitational Waves from Generic Spin-Precessing Compact Binary Inspirals, Phys. Rev. D 95, 104004 (2017), arXiv:1703.03967 [gr-qc] .
  * Gupta  _et al._ [2020] A. Gupta, S. Datta, S. Kastha, S. Borhanian, K. G. Arun, and B. S. Sathyaprakash, Multiparameter tests of general relativity using multiband gravitational-wave observations, Phys. Rev. Lett. 125, 201101 (2020), arXiv:2005.09607 [gr-qc] .
  * Blanchet [2002] L. Blanchet, On the accuracy of the postNewtonian approximation, in  _25th Johns Hopkins Workshop on Current Problems in Particle Theory: 2001: A Relativistic Spacetime Odyssey. Experiments and Theoretical Viewpoints on General Relativity and Quantum Gravity_ (2002) pp. 411–430, arXiv:gr-qc/0207037 .
  * Yunes and Berti [2008] N. Yunes and E. Berti, Accuracy of the post-Newtonian approximation: Optimal asymptotic expansion for quasicircular, extreme-mass ratio inspirals, Phys. Rev. D 77, 124006 (2008), [Erratum: Phys.Rev.D 83, 109901 (2011)], arXiv:0803.1853 [gr-qc] .
  * Damour  _et al._ [2001] T. Damour, B. R. Iyer, and B. S. Sathyaprakash, A Comparison of search templates for gravitational waves from binary inspiral, Phys. Rev. D 63, 044023 (2001), [Erratum: Phys.Rev.D 72, 029902 (2005)], arXiv:gr-qc/0010009 .
  * Damour  _et al._ [2002] T. Damour, B. R. Iyer, and B. S. Sathyaprakash, A Comparison of search templates for gravitational waves from binary inspiral - 3.5PN update, Phys. Rev. D 66, 027502 (2002), arXiv:gr-qc/0207021 .
  * Arun  _et al._ [2005] K. G. Arun, B. R. Iyer, B. S. Sathyaprakash, and P. A. Sundararajan, Parameter estimation of inspiralling compact binaries using 3.5 post-Newtonian gravitational wave phasing: The Non-spinning case, Phys. Rev. D 71, 084008 (2005), [Erratum: Phys.Rev.D 72, 069903 (2005)], arXiv:gr-qc/0411146 .
  * Aasi  _et al._ [2015] J. Aasi _et al._ (LIGO Scientific), Advanced LIGO, Class. Quant. Grav. 32, 074001 (2015), arXiv:1411.4547 [gr-qc] .
  * Acernese  _et al._ [2015] F. Acernese _et al._ (VIRGO), Advanced Virgo: a second-generation interferometric gravitational wave detector, Class. Quant. Grav. 32, 024001 (2015), arXiv:1408.3978 [gr-qc] .
  * O’Reilly  _et al._ [2020] B. O’Reilly, M. Branchesi, S. Haino, and G. Gemme,  _LIGO Document T2000012-v1_, Tech. Rep. (2020).
  * Dwyer  _et al._ [2015] S. Dwyer, D. Sigg, S. W. Ballmer, L. Barsotti, N. Mavalvala, and M. Evans, Gravitational wave detector with cosmological reach, Phys. Rev. D 91, 082001 (2015), arXiv:1410.0612 [astro-ph.IM] .
  * Punturo  _et al._ [2010] M. Punturo _et al._ , The third generation of gravitational wave observatories and their science reach, Class. Quant. Grav. 27, 084007 (2010).
  * Cosmic Explorer [2020] Cosmic Explorer, https://cosmicexplorer.org/researchers.html (2020).
  * Hild  _et al._ [2011] S. Hild _et al._ , Sensitivity Studies for Third-Generation Gravitational Wave Observatories, Class. Quant. Grav. 28, 094013 (2011), arXiv:1012.0908 [gr-qc] .
  * Kanti  _et al._ [1996] P. Kanti, N. E. Mavromatos, J. Rizos, K. Tamvakis, and E. Winstanley, Dilatonic black holes in higher curvature string gravity, Phys. Rev. D 54, 5049 (1996), arXiv:hep-th/9511071 .
  * Boulware and Deser [1985] D. G. Boulware and S. Deser, String Generated Gravity Models, Phys. Rev. Lett. 55, 2656 (1985).
  * Gross and Sloan [1987] D. J. Gross and J. H. Sloan, The quartic effective action for the heterotic string, Nuclear Physics B 291, 41 (1987).
  * Yunes and Stein [2011] N. Yunes and L. C. Stein, Non-Spinning Black Holes in Alternative Theories of Gravity, Phys. Rev. D 83, 104002 (2011), arXiv:1101.2921 [gr-qc] .
  * Yagi  _et al._ [2016] K. Yagi, L. C. Stein, and N. Yunes, Challenging the Presence of Scalar Charge and Dipolar Radiation in Binary Pulsars, Phys. Rev. D 93, 024010 (2016), arXiv:1510.02152 [gr-qc] .
  * Sotiriou and Zhou [2014] T. P. Sotiriou and S.-Y. Zhou, Black hole hair in generalized scalar-tensor gravity: An explicit example, Phys. Rev. D 90, 124063 (2014), arXiv:1408.1698 [gr-qc] .
  * Yagi  _et al._ [2012] K. Yagi, L. C. Stein, N. Yunes, and T. Tanaka, Post-Newtonian, Quasi-Circular Binary Inspirals in Quadratic Modified Gravity, Phys. Rev. D 85, 064022 (2012), [Erratum: Phys.Rev.D 93, 029902 (2016)], arXiv:1110.5950 [gr-qc] .
  * Nair  _et al._ [2019] R. Nair, S. Perkins, H. O. Silva, and N. Yunes, Fundamental Physics Implications for Higher-Curvature Theories from Binary Black Hole Signals in the LIGO-Virgo Catalog GWTC-1, Phys. Rev. Lett. 123, 191101 (2019), arXiv:1905.00870 [gr-qc] .
  * Lyu  _et al._ [2022] Z. Lyu, N. Jiang, and K. Yagi, Constraints on Einstein-dilation-Gauss-Bonnet gravity from black hole-neutron star gravitational wave events, Phys. Rev. D 105, 064001 (2022), arXiv:2201.02543 [gr-qc] .
  * Saleem  _et al._ [2021] M. Saleem, S. Datta, K. G. Arun, and B. S. Sathyaprakash, Parametrized tests of post-Newtonian theory using principal component analysis, arXiv:2110.10147 [gr-qc] (2021).
  * Shoom  _et al._ [2021] A. A. Shoom, P. K. Gupta, B. Krishnan, A. B. Nielsen, and C. D. Capano, Testing GR with the Gravitational Wave Inspiral Signal GW170817, arXiv:2105.02191 [gr-qc] (2021).
  * Pai and Arun [2013] A. Pai and K. G. Arun, Singular value decomposition in parametrised tests of post-Newtonian theory, Class. Quant. Grav. 30, 025011 (2013), arXiv:1207.1943 [gr-qc] .
  * Khalil  _et al._ [2019] M. Khalil, N. Sennett, J. Steinhoff, and A. Buonanno, Theory-agnostic framework for dynamical scalarization of compact binaries, Phys. Rev. D 100, 124013 (2019), arXiv:1906.08161 [gr-qc] .
  * Silva  _et al._ [2021b] H. O. Silva, H. Witek, M. Elley, and N. Yunes, Dynamical Descalarization in Binary Black Hole Mergers, Phys. Rev. Lett. 127, 031101 (2021b), arXiv:2012.10436 [gr-qc] .