footnote

# Black-hole ringdown as a probe of higher-curvature gravity theories

Hector O. Silva  Max Planck Institute for Gravitational Physics (Albert
Einstein Institute), Am Mühlenberg 1, Potsdam D-14476, Germany  Abhirup Ghosh
Max Planck Institute for Gravitational Physics (Albert Einstein Institute), Am
Mühlenberg 1, Potsdam D-14476, Germany  Alessandra Buonanno  Max Planck
Institute for Gravitational Physics (Albert Einstein Institute), Am Mühlenberg
1, Potsdam D-14476, Germany  Department of Physics, University of Maryland,
College Park, Maryland 20742, USA

###### Abstract

Detecting gravitational waves from coalescing compact binaries allows us to
explore the dynamical, nonlinear regime of general relativity and constrain
modifications to it. Some of the gravitational-wave events observed by the
LIGO-Virgo Collaboration have sufficiently high signal-to-noise ratio in the
merger, allowing us to probe the relaxation of the remnant black hole to its
final, stationary state — the so-called black-hole ringdown, which is
characterized by a set of quasinormal modes. Can we use the ringdown to
constrain deviations from general relativity, as predicted by several of its
contenders? Here, we address this question by using an inspiral-merger-
ringdown waveform model in the effective-one-body formalism, augmented with a
parametrization of the ringdown based on an expansion in the final black
hole’s spin. We give a prescription on how to include in this waveform model,
the quasinormal mode frequencies calculated on a theory-by-theory basis. In
particular, we focus on theories that modify general relativity by higher-
order curvature corrections, namely, Einstein-dilaton-Gauss-Bonnet, dynamical
Chern-Simons theories, and cubic- and quartic-order effective-field-theories
of general relativity. We use this parametrized waveform model to measure the
ringdown properties of the two loudest ringdown signals observed so far,
GW150914 and GW200129. We find that while the Einstein-dilaton-Gauss-Bonnet
theory cannot be constrained with these events, we can place upper bounds on
the fundamental lengthscale of cubic- ($\ell_{\rm cEFT}\leqslant 38.2$ km) and
quartic-order ($\ell_{\rm qEFT}\leqslant 51.3$ km) effective-field-theories of
general relativity, and of dynamical Chern-Simons gravity ($\ell_{\rm
dCS}\leqslant 38.7$ km). The latter result is a concrete example of a theory
presently unconstrained by inspiral-only analyses which, however, can be
constrained by merger-ringdown studies with current gravitational-wave data.

##  I Introduction

Since the first detection of gravitational waves (GWs) from a binary black-
hole (BBH) merger in 2015 Abbott  _et al._ (2016a), the LIGO Aasi  _et al._
(2015) and Virgo Acernese  _et al._ (2015) detectors have observed about 90 GW
events Abbott  _et al._ (2021a) from mergers of BHs, neutron stars (NSs)
Abbott  _et al._ (2017a, 2018, 2020) and their mixture Abbott  _et al._
(2021b). These results have been confirmed by independent analyses, which have
also identified a few new GW signals Nitz  _et al._ (2019a, b); Venumadhav
_et al._ (2020); Zackay  _et al._ (2021); Nitz  _et al._ (2021); Olsen  _et
al._ (2022).

The large number of GW observations has allowed us to infer relevant
astrophysical Abbott  _et al._ (2021c) and cosmological Abbott  _et al._
(2021d) information on the compact-object population in our local Universe,
and also to probe general relativity (GR) in the high-velocity, dynamical and
strong-field regime of gravity Yunes  _et al._ (2016); Abbott  _et al._
(2021e). The latter complement tests of GR in the low-velocity, quasistatic or
linear regimes available with Solar-System experiments Will (2014), binary-
pulsar Wex (2014); Kramer  _et al._ (2021) and galactic-center Abuter  _et
al._ (2018); Do  _et al._ (2019); Akiyama  _et al._ (2019) observations, and
cosmological measurements Clifton  _et al._ (2012).

The coalescence of two BHs in GR is characterized by a long inspiral stage,
during which the holes adiabatically and steadily come closer and closer to
each other, losing energy because of the emission of GWs. Then, they merge,
forming a common apparent horizon. Subsequently, during the ringdown stage,
the newly formed remnant object settles down to a Kerr BH emitting quasinormal
modes (QNMs) Vishveshwara (1970); Press (1971); Chandrasekhar and Detweiler
(1975). Because of the no-hair conjecture in GR Carter (1971); Israel (1967);
Hawking (1972); Robinson (1975), the QNM (complex) frequencies of
(electrically neutral) astrophysical BHs are only described by the BH’s mass
and spin. In GR, the QNM frequencies are labeled by the harmonic indices
$(\ell,m)$ and the overtone number $n$.

Several null tests have been proposed to probe the nature of gravity with GW
signals Abbott  _et al._ (2016b); Yunes  _et al._ (2016); Ghosh  _et al._
(2016, 2018); Abbott  _et al._ (2019a, b, 2021f, 2021e). They include tests of
GW generation Arun  _et al._ (2006); Yunes and Pretorius (2009a); Li  _et al._
(2012); Agathos  _et al._ (2014); Mehta  _et al._ (2022), where deviations in
the post-Newtonian (PN) coefficients in the inspiral, and phenomenological
coefficients in the plunge and merger stages can be bounded; tests of GW
propagation Will (1998); Mirshekari  _et al._ (2012), which allow us to set
upper limits on coefficients entering generalized dispersion relations,
including the Compton wavelength associated to the mass of the graviton; tests
of the polarization of gravitational radiation Eardley  _et al._ (1973); Will
(2014); Isi and Weinstein (2017); Abbott  _et al._ (2017b, 2019c); Wong  _et
al._ (2021), for which more than two GW detectors are needed to set
statistically significant bounds, and tests of the remnant properties Meidam
_et al._ (2014); Carullo  _et al._ (2018); Brito  _et al._ (2018); Carullo
_et al._ (2019); Isi  _et al._ (2019); Ghosh  _et al._ (2021); Carullo (2021)
in the postmerger stage. So far, none of these null tests have reported any
deviation from GR.

Probing the gravitational properties of the remnant object during the
ringdown, has attracted a lot of attention in the last twenty years. Reference
Dreyer  _et al._ (2004) proposed the idea of employing BH spectroscopy
Detweiler (1980) of the ringdown stage to rule out (or constrain) either
modified theories of gravity or exotic compact objects (in GR) rather than
BHs, thus testing the no-hair conjecture. Since then, many studies have
quantified the accuracy with which the QNM frequencies can be measured for GW
sources detectable with ground- and space-based detectors (see, e.g., Refs.
Berti  _et al._ (2006); Baibhav  _et al._ (2020); Bhagwat  _et al._ (2022);
Ota and Chirenti (2022)). Several analyses Carullo  _et al._ (2018); Brito
_et al._ (2018); Carullo  _et al._ (2019); Isi  _et al._ (2019); Ghosh  _et
al._ (2021); Carullo (2021) have used the GW observations of the LIGO-Virgo-
KAGRA (LVK) collaboration to set upper limits on deviations in the QNM
frequencies of BHs in GR. Others have claimed the measurements of QNMs beyond
the dominant $(2,2,0)$ mode Capano  _et al._ (2021), or overtones — for
example the $(2,2,1)$ mode Isi  _et al._ (2019); Isi and Farr (2022); Finch
and Moore (2022) , although the evidence for overtones can also be due to
noise Cotesta  _et al._ (2022). These studies have been pursued either using a
superposition of damped sinusoids Carullo  _et al._ (2019); Isi and Farr
(2021), in some cases augmented with QNM amplitudes calibrated to numerical-
relativity (NR) simulations, or with parametrized inspiral-merger-ringdown
(IMR) waveform models, where the QNM frequencies are not necessarily fixed to
the GR values for BHs, but kept free Brito  _et al._ (2018); Ghosh  _et al._
(2021).

Here, we will employ the parametrized IMR model of Ref. Ghosh  _et al._
(2021), constructed from a nonprecessing-spin effective-one-body (EOB)
waveform model Bohé  _et al._ (2017); Cotesta  _et al._ (2018); Mihaylov  _et
al._ (2021), to carry out theory-specific tests of the ringdown using four
high-curvature gravity theories. Previously, such parametrized waveform model
was employed in Refs. Abbott  _et al._ (2021f, e) for theory-independent tests
of the ringdown. More specifically, here we will focus on four modified
gravity theories, Einstein-dilaton-Gauss-Bonnet gravity, dynamical Chern-Simon
gravity, cubic and quartic effective-field theories (EFTs) of GR, and express
the QNM frequencies using the parametrized ringdown spin-expansion
coefficients (ParSpec) of Ref. Maselli  _et al._ (2020). In this framework,
the non-GR QNM frequencies are recast as deviations from the GR QNM values,
and are expressed in terms of a single free parameter, the fundamental
lengthscale $\ell_{\rm th}$ of the gravity theory under consideration, the GR
limit corresponding to $\ell_{\rm th}\rightarrow 0$.

With this formulation of the ringdown, we use the two loudest merger-ringdown
GW events, so far observed by the LVK collaboration, notably GW150914 and
GW200129, and use Bayesian-inference techniques to perform null tests. We find
no indication that GR is violated and, when possible, we place upper limits,
at $90\%$ credible level, on the lengthscale $\ell_{\rm th}$ of each theory.
In Table 1, we summarize our results, and compare them with existing
constraints.

Theory | Constraint | This work
---|---|---
EdGB |  $\ell_{\rm EdGB}\leqslant 1.18$ km (GW) Lyu  _et al._ (2022) | –
dCS |  $\ell_{\rm dCS}\leqslant 8.5$ km (EM+GW) Silva  _et al._ (2021a) |  $\ell_{\rm dCS}\leqslant 38.7$ km
cubic EFT | – |  $\ell_{\rm cEFT}\leqslant 38.2$ km
quartic EFT |  $\ell_{\rm qEFT}\leqslant 150$ km Sennett  _et al._ (2020) |  $\ell_{\rm qEFT}\leqslant 51.3$ km a

Table 1: Summary of the upper limits, at $90\%$ credible level, on the
lengthscale of the modified gravity theories under investigation. They were
obtained combining Bayesian-inference results from GW150914 and GW200129.
Current constraints on $\ell_{\rm th}$ are also listed.

  * a

We notice that our result for the quartic EFT of GR is only in marginal
tension with our hard cutoff scale for the validity of the theory, and hence
we do still quote it here (see Secs. IV.4 and V.4 for details).

The paper is organized as follows. In Sec. II, we briefly describe the four
higher-curvature modified gravity theories for which we perform the ringdown
test. In Sec. III we build our parametrized IMR model for nonprecessing-spin
compact-object binaries making use of the ParSpec framework. After reviewing
the Bayesian inference method, in Sec. IV, we motivate our selection of GW
events from the LVK catalog, and also discuss the range of validity of our
analyses. In particular, we discuss the impact on our results of the
assumptions underlying the ParSpec framework, and the fact that our modified
gravity theories have to be interpreted as EFTs. In Sec. V, we present our
results obtained by applying Bayesian analysis on the LVK data of GW150914 and
GW1200129, and discuss how we set the upper limits on the fundamental
lengthscales $\ell_{\rm th}$. Finally, in Sec. VI we summarize our findings,
and discuss how to make our framework more robust, in view of stronger GW
events in upcoming GW observational runs, by including physical effects
currently absent in our study (e.g., precessing-spins and eccentricity). In
the Appendix A we provide details in calculating the non-GR QNM frequencies,
when using ParSpec, for the modified gravity theories under consideration.
Henceforth, unless otherwise specified, we work in geometric units $G=1=c$.

##  II Overview of modified gravity theories

We will treat the modified theories of gravity as EFTs, and focus on finite-
size effects (see, e.g., Ref. Sennett  _et al._ (2020)). Thus, for each
gravity theory we impose that the fundamental lengthscale $\ell_{\rm
th}\lesssim GM/c^{2}$, where $M$ is the mass of the BH. This implies that
observable deviations from GR present in those theories arise from
modifications to the Kerr geometry of each individual BH. Those finite-size
effects can manifest themselves in the QNMs of the remnant produced by the
merger, but also in the GW phasing of the inspiral through corrections to the
GR spin-induced quadrupole and Love numbers. However, here we will not
consider the latter, instead, we will only study the impact of finite-size
effects on the QNMs of the remnant.

We start by briefly reviewing the modified gravity theories we consider in
this paper, what the current observational constraints are and what we know
about BH QNMs in each of these theories.

###  II.1 Einstein-dilaton-Gauss-Bonnet gravity

This theory belongs to the class of scalar-Gauss-Bonnet theories, which are
described by the action

$\displaystyle S_{\rm EdGB}=\frac{1}{16\pi}\int{\rm
d}^{4}x\,\sqrt{-g}\,\left[R-\tfrac{1}{2}(\partial\varphi)^{2}+\tfrac{1}{4}\ell^{2}_{\rm
EdGB}f(\varphi)\,\mathcal{G}\right],$

where $g\equiv\textrm{det}(g_{\mu\nu})$ is the metric determinant, $R$ is the
Ricci scalar, $\varphi$ is a dynamical scalar field, with kinetic term
$(\partial\varphi)^{2}=g^{\mu\nu}\partial_{\mu}\varphi\partial_{\nu}\varphi$,
which couples to the Gauss-Bonnet invariant
$\mathcal{G}=R^{\mu\nu\rho\sigma}R_{\mu\nu\rho\sigma}-4R^{\mu\nu}R_{\mu\nu}+R^{2}$,
and $R_{\mu\nu\rho\sigma}$ and $R_{\mu\nu}$ are the Riemann and Ricci tensors
respectively. By itself, $\int{\rm d}^{4}x\,\sqrt{-g}\,\mathcal{G}$ is a
boundary term in four dimensions and hence does not contribute to the field
equations Myers (1987). However, when coupled to $\varphi$, it can contribute
to the field equation through the coupling function $f(\varphi)$. The strength
of the coupling is set by $\ell_{\rm GB}$, with dimensions of length.

Different subclasses of this theory are determined by the function
$f(\varphi)$ and can be divided into two classes based on the properties of
their BH solutions. In the first class, the first derivative of the coupling
function $f^{\prime}(\varphi)\equiv{\rm d}f/{\rm d}\varphi$ is always nonzero
and BHs are known to always support scalar hair. This is the case of Einstein-
dilaton-Gauss-Bonnet (EdGB) gravity, for which $f(\varphi)=\exp(\varphi)$
Kanti  _et al._ (1996). In the second class, $f^{\prime}(\varphi)=0$ can
vanish for some constant $\varphi_{0}$. In this case, the theory admits the
same stationary, asymptotically flat BH solutions as GR Silva  _et al._ (2018)
and those of scalarized BHs Doneva and Yazadjiev (2018); Silva  _et al._
(2018); Macedo  _et al._ (2019); Dima  _et al._ (2020); Herdeiro  _et al._
(2021); Berti  _et al._ (2021). Examples include Gaussian
$f(\varphi)\propto\exp(-\varphi^{2})$ Doneva and Yazadjiev (2018) and the
quadratic $f(\varphi)\propto\varphi^{2}$ Silva  _et al._ (2018) coupling
functions.

BHs in EdGB gravity have scalar hair, to which we can associate a monopole
scalar charge, related to the asymptotic $r^{-1}$ fall-off of the scalar
field, where $r$ is the distance from the BH. This charge is not an
independent parameter, and depends on the BH’s mass and spin, thus being a
“secondary hair” Coleman  _et al._ (1992); Kanti  _et al._ (1996); Herdeiro
and Radu (2015). Since the scalar field is sourced by a curvature scalar, the
scalar charge is larger (smaller), the smaller (larger) the BH mass is.

These properties are not mere theoretical curiosities; they have important
observational consequences. First, the presence of the scalar charge implies
that when in binaries, BHs can source scalar-dipole radiation (see, e.g.,
Refs. Yagi  _et al._ (2012); Julié and Berti (2019); Shiralilou  _et al._
(2021, 2022); Julié  _et al._ (2022)) which affects the GW phase at $-1$PN
order (relative to the dominant quadrupolar GR contribution), with magnitude
proportional to the difference between the charges of binary components. This
makes EdGB gravity testable with GW observations of compact binaries where at
least one component is a BH. Indeed, Ref. Lyu  _et al._ (2022) placed the
bound $\ell_{\rm EdGB}\leqslant 7.1$ km using the NSBH binaries GW200105 and
GW200115 Abbott  _et al._ (2021b) while Refs. Nair  _et al._ (2019); Perkins
_et al._ (2021), obtained $\ell_{\rm EdGB}\leqslant 9.1$ km, by stacking the
posteriors of $\ell_{\rm EdGB}$ from a selection of BBHs from the GWTC-1 and
GWTC-2 catalogs Abbott  _et al._ (2019d, 2021g)111These bounds are, strictly
speaking, valid only when the scalar field $\varphi$ is small (i.e.,
$\varphi\ll 1$), and they take into account only the leading-order scalar
field interaction arising from the original dilatonic coupling (i.e.,
$f(\varphi)\approx\varphi$ in Eq. (II.1)). This results in what is often
referred to as shift-symmetric scalar-Gauss-Bonnet theory. In this theory, NSs
do not have scalar monopole charge Yagi  _et al._ (2016), while BHs do Yunes
and Stein (2011); Sotiriou and Zhou (2014a, b). Finally, we note that the
constant $\alpha_{\rm EdGB}$ used in Refs. Yagi  _et al._ (2012); Nair  _et
al._ (2019); Perkins  _et al._ (2021); Lyu  _et al._ (2022) is related to
$\ell_{\rm EdGB}$ as $\ell_{\rm EdGB}=4\pi^{1/4}|\alpha_{\rm EdGB}|^{1/2}$..
Secondly, the scalar field influences the response of BHs to linearized
perturbations and thus affects the BH’s QNM spectra. The coupling between
scalar field and the Gauss-Bonnet invariant, results in a coupling between
scalar perturbations and gravitational perturbations of polar parity which,
for instance, breaks the equivalence between the QNM spectra of polar and
axial gravitational perturbations Glampedakis  _et al._ (2017a); Lenzi and
Sopuerta (2021) (i.e., the isospectrality Chandrasekhar (1985)) of
Schwarzschild BHs in GR.

Using BH perturbation theory, the QNMs for nonrotating BHs in EdGB gravity
were first computed in Refs. Pani and Cardoso (2009); Blázquez-Salcedo  _et
al._ (2016) and were extended (in the polar sector) to leading-order in BH
spin in Ref. Pierini and Gualtieri (2021). See Ref. Bryant  _et al._ (2021)
for a study in the geometrical optics limit ($\ell\gg 1$) and Refs. Witek  _et
al._ (2019); Okounkova (2019) for NR studies.

###  II.2 Dynamical Chern-Simons gravity

This theory is described by the action Jackiw and Pi (2003); Alexander and
Yunes (2009),

$\displaystyle S_{\rm dCS}=\int{\rm
d}^{4}x\,\sqrt{-g}\,\left[\tfrac{1}{16\pi}\,R-\tfrac{1}{2}(\partial\vartheta)^{2}+\tfrac{1}{4}\ell^{2}_{\rm
dCS}\,\vartheta\,{}^{*}RR\right],$ (2)

where $\vartheta$ is a pseudoscalar field which couples to the Pontryagin
density ${}^{*}RR={}^{*}R_{\mu\nu}{}^{\rho\sigma}R^{\nu\mu}{}_{\rho\sigma}$,
where ${}^{*}R^{\mu\nu\rho\sigma}$ is the dual of the Riemann tensor defined
as
${}^{*}R^{\mu\nu\rho\sigma}=\epsilon^{\mu\nu}{}_{\gamma\delta}R^{\gamma\delta\rho\sigma}/2$,
and $\epsilon_{\mu\nu\gamma\delta}$ is the Levi-Civita tensor. The variation
of the Pontryagin density is a boundary term that does not contribute to the
field equations in four-dimensions Jackiw and Pi (2003). However, the
Pontryagin density can modify the field equations when coupled to $\vartheta$;
the strength of the coupling set by $\ell_{\rm dCS}$, with dimensions of
length.

The theory admits as a solution the garden-variety Schwarzschild BH of GR.
This is not the case when rotation is included and the Kerr metric is not a
solution of the theory Jackiw and Pi (2003). These rotating BHs support a
scalar field which falls off as $r^{-2}$ asymptotically, to which we can
associate a scalar dipole charge Yunes and Pretorius (2009b); Konno  _et al._
(2009) and the leading-order modification to the GW phase enters at 2PN Yagi
_et al._ (2012). Deviations from GR at this PN order are constrained with
present GW observations Abbott  _et al._ (2021e) only at the level of $\sim
50\%$ (see the constraint on the 2PN parameter $\varphi_{4}$ in Fig. 6 of Ref.
Abbott  _et al._ (2021e)). So far, analyses that used only the inspiral
portion of the BBH GW signals were not able to set meaningful bounds on such
deviation at 2PN order Nair  _et al._ (2019); Perkins  _et al._ (2021). These
works, as well as the analysis we do here, probe the effects of dCS in the
generation of GWs. Weak constraints of order $\sim 10^{3}$ km were obtained on
this theory by considering parity-violation propagation effects in GWs Yamada
and Tanaka (2020); Wang  _et al._ (2021); Okounkova  _et al._ (2022), which
show up as an amplitude birefringence between different GW polarizations
Jackiw and Pi (2003); Alexander and Martin (2005). Nonetheless, the theory has
been constrained in Ref. Silva  _et al._ (2021a), which found $\ell_{\rm
dCS}\leqslant 8.5$ km, by folding data from the X-ray observations of the
pulsar PSR J0030+0451 Lommen  _et al._ (2000); Arzoumanian  _et al._ (2018) by
NICER Riley  _et al._ (2019); Miller  _et al._ (2019) and from the GW
observation of the binary NS GW170817 Abbott  _et al._ (2017a, 2018), using
equation-of-state independent relations between NS moment of inertia and tidal
deformability Yagi and Yunes (2013a, b); Gupta  _et al._ (2018).

The QNMs of the Schwarzschild BH in dCS gravity were studied in Refs. Yunes
and Sopuerta (2008); Cardoso and Gualtieri (2009); Molina  _et al._ (2010),
which found that scalar perturbations couple to gravitational perturbations of
axial parity, in contrast with EdGB gravity, resulting in a breakdown of
isospectrality. The QNM spectra of slowly-rotating BHs in dCS gravity was
studied in Refs. Srivastava  _et al._ (2021); Wagle  _et al._ (2022). They
were also extracted from NR simulations of BH head-on collisions in Ref.
Okounkova  _et al._ (2019).

###  II.3 Effective-field-theory of General Relativity

Our last example of modified gravity theories are the so-called EFTs of GR
Endlich  _et al._ (2017); Sennett  _et al._ (2020); Cano and Ruipérez (2019);
Brandhuber and Travaglini (2020); Accettulli Huber  _et al._ (2021); de Rham
_et al._ (2020); Cano  _et al._ (2020, 2022). They are described by the action

$S_{\rm EFT}=\frac{1}{16\pi}\int{\rm d}^{4}x\sqrt{-g}\left[R+\sum_{n\geqslant
2}\ell_{\rm EFT}^{2n-2}\,L^{(2n)}\right]\,,$ (3)

where $\ell_{\rm EFT}$ is a lengthscale assumed to be small compared to the
lengthscale $M$ associated with a BH (i.e., $\ell_{\rm EFT}/M\ll 1$), and
$L^{(2n)}$ are corrections that introduce higher-order curvature tensors (with
$2n$ metric derivatives).

More specifically, we follow the notation of Refs. Cano  _et al._ (2020, 2022)
and consider up to dimension-eight operators ($n=4$),

$\displaystyle L^{(6)}$ $\displaystyle=\lambda_{\rm
e}R_{\mu\nu}{}^{\rho\sigma}R_{\rho\sigma}{}^{\gamma\delta}R_{\gamma\delta}{}^{\mu\nu}+\lambda_{\rm
o}R_{\mu\nu}{}^{\rho\sigma}R_{\rho\sigma}{}^{\gamma\delta}\tilde{R}_{\gamma\delta}{}^{\mu\nu},$
(4a) $\displaystyle L^{(8)}$
$\displaystyle=\varepsilon_{1}\mathcal{C}^{2}+\varepsilon_{2}\tilde{\mathcal{C}}^{2}+\varepsilon_{3}\mathcal{C}\tilde{\mathcal{C}},$
(4b)

where ${\cal C}=R_{\mu\nu\rho\sigma}R^{\mu\nu\rho\sigma}$, $\tilde{\cal
C}=R_{\mu\nu\gamma\delta}\epsilon^{\mu\nu}{}_{\rho\sigma}R^{\rho\sigma\gamma\delta}$,
and both $\lambda_{\rm o,e}$ and $\varepsilon_{i}$ (with $i=1,2,3$) are
dimensionless parameters. Due to the large number of free parameters in this
theory, we focus on a subset of the parameter space. In particular, we
consider dimension-six and dimension-eight operators separately. In addition,
in the dimension-six case we further assume that $\lambda_{\rm e}=\lambda_{\rm
o}=1$, leaving us with $\ell_{\rm cEFT}$ as our single free parameter.
Similarly, in the dimension-eight case, we set $\varepsilon_{1}=1$ and
$\varepsilon_{2}=\varepsilon_{3}=0$, as done in Ref. Sennett  _et al._ (2020).
This leaves us with $\ell_{\rm qEFT}$ as our single free parameter.

For the EFT of GR with dimension-eight operators, Ref. Sennett  _et al._
(2020) focused on the orbital effects (i.e., instead of finite-size effects)
and performed Bayesian model selection using the two lowest-mass BBHs events
of the second-observig run of the LIGO-Virgo Collaboration, notably GW151226
and GW170608. They found that the data disfavor the appearance of new physics
on distance scales around $\ell_{\rm qEFT}\sim 150$ km.

The QNMs of nonrotating BHs where calculated in Refs. de Rham  _et al._
(2020); Cano  _et al._ (2020) in the dimension-six EFT and in Refs. Cardoso
_et al._ (2018); Cano  _et al._ (2020) in the dimension-eight EFT. Reference
Cano  _et al._ (2020) calculated the leading-order BH spin corrections to the
QNM spectra.

##  III Methods

Having reviewed the modified gravity theories that we will consider, we now
present the sequential building blocks for the waveform model we use to test
these theories against GW observations. We start by reviewing the parametrized
ringdown spin-expansion coefficients (ParSpec) framework Maselli  _et al._
(2020) in Sec. III.1. Next, in Sec. III.2, we review our baseline parametrized
IMR waveform model Brito  _et al._ (2018); Ghosh  _et al._ (2021), and explain
how we extend it to include the ParSpec. Finally, in Sec. III.3, we show how
we can map theory-specific QNM calculations in modified gravity theories onto
the free coefficients in the ParSpec framework. Ultimately, this provides us
with an IMR waveform model, with the ringdown portion of the model informed by
QNM calculations in specific beyond-GR theories.

###  III.1 The parametrized ringdown spin expansion coefficients framework

A general procedure to describe deviations to the QNM frequencies $\omega^{\rm
GR}_{\ell mn}$ and damping times $\tau^{\rm GR}_{\ell mn}$ of BHs of GR is to
write,

$\displaystyle\omega_{\ell mn}$ $\displaystyle=\omega_{\ell mn}^{\rm
GR}\,(1+\delta\omega_{\ell mn}),$ (5a) $\displaystyle\tau_{\ell mn}$
$\displaystyle=\tau_{\ell mn}^{\rm GR}\,(1+\delta\tau_{\ell mn}),$ (5b)

where $\delta\omega_{\ell mn}$ and $\delta\tau_{\ell mn}$ are the fractional
deformation parameters, $\ell$ and $m$ are the multipole indices, and $n$ the
overtone number. This type of parametrization222See also Refs. Glampedakis
_et al._ (2017b); Glampedakis and Silva (2019); Silva and Glampedakis (2020)
and Refs. Cardoso  _et al._ (2019); McManus  _et al._ (2019); Kimura (2020)
for alternative parametrizations. was adopted, for instance, in Refs. Gossan
_et al._ (2012); Meidam  _et al._ (2014); Brito  _et al._ (2018); Carullo  _et
al._ (2018); Isi and Farr (2021); Ghosh  _et al._ (2021).

The current LVK tests of BH ringdown (see Sec.VII.A of Ref. Abbott  _et al._
(2021h) or Sec.VIII.A of Ref. Abbott  _et al._ (2021e)) take a flexible
theory-independent approach towards the inference of $\delta\omega_{\ell mn}$
and $\delta\tau_{\ell mn}$. These deviations are either assumed to occur
identically across all observed sources or belong to a generic underlying
Gaussian population. However, in reality, these parameters depend on the
source BH’s mass and spin. Ideally, one would like to explicitly reinstate
this dependence, by (i) introducing deformation parameters which can be
determined, once and for all, from a specific gravity theory (GR included),
and (ii) making it simpler to combine constraints coming from multiple
(independently observed) sources.

The ParSpec framework was introduced in Ref. Maselli  _et al._ (2020) and can
be used to our purpose. It is an observable-based bivariate expansion of Eq.
(5), given by

$\displaystyle\omega_{\ell mn}$ $\displaystyle=\frac{1}{M_{\rm
f}}\sum_{j=0}^{N_{\rm max}}\,\chi_{\rm f}^{j}\omega^{(j)}_{\ell
mn}\,\left(1+\gamma\,\delta\omega^{(j)}_{\ell mn}\right),$ (6a)
$\displaystyle\tau_{\ell mn}$ $\displaystyle=M_{\rm f}\sum_{j=0}^{N_{\rm
max}}\,\chi_{\rm f}^{j}\tau^{(j)}_{\ell
mn}\,\left(1+\gamma\,\delta\tau^{(j)}_{\ell mn}\right),$ (6b)

where $M_{\rm f}$ and $\chi_{\rm f}$ are the detector-frame final mass and
spin, respectively; the quantities $\omega_{\ell mn}^{(j)}$ and $\tau_{\ell
mn}^{(j)}$ are dimensionless coefficients of the expansion in spin for the
QNMs of BHs in GR, while $\delta\omega_{\ell mn}^{(j)}$ and $\delta\tau_{\ell
mn}^{(j)}$ are source-independent dimensionless coefficients that characterize
the corrections to the GR QNM at each spin-order, and $N_{\rm max}$ is the
order of the spin expansion. All source dependence due to a given modified
gravity theory is contained in the dimensionless parameter $\gamma$, which
reads

$\gamma=\left(\frac{\ell_{\rm th}}{M^{\rm s}_{\rm
f}}\right)^{p}=\left[\frac{\ell_{\rm th}c^{2}(1+z)}{GM_{\rm f}}\right]^{p},$
(7)

which depends on the lengthscale parameter $\ell_{\rm th}$ of the specific
gravity theory (non-GR modifications become important at distances
$\lesssim\ell_{\rm th}$), and the exponent $p$ is related to how the non-GR
modifications are added to the Einstein-Hilbert action. In Eq. (7) we made
$\gamma$ dimensionless by the lengthscale associated with the remnant BH
(i.e., its source-frame mass $M^{\rm s}_{\rm f}$333Throughout this paper, we
maintain the convention of using the upperscript “s” to denote source-frame
quantities and plain baseline symbols for detector-frame measurements.), which
we can also write in terms of the detector-frame mass $M_{\rm f}$ through the
redshift $z$ Krolak and Schutz (1987). Also, dividing by the factor $G/c^{2}$
allows us to express $\ell_{\rm th}$ in physically intuitive metric units.

In principle, a modification to GR would also affect $M_{\rm f}$ and
$\chi_{\rm f}$ and the expansion should be written in terms of the non-GR mass
and spin, say $\bar{M}_{\rm f}$ and $\bar{\chi}_{\rm f}$. If we assume that
the non-GR corrections are included perturbatively (as it is the case with all
the theories described in Sec. II), the modifications to the BH mass and spin
can be absorbed into the deviations parameters $\delta\omega^{(j)}_{\ell mn}$
and $\delta\tau^{(j)}_{\ell mn}$. This means we can identify the $M_{\rm f}$
and $\chi_{\rm f}$ with their _corresponding GR values_444For a more detailed
discussion, see Appendix A in Ref. Maselli  _et al._ (2020).. We will see in
Sec. IV that this assumption is indeed satisfied in our parameter estimation
studies (see, for instance, Fig. 7).

Finally, we remark that in the GR limit ($\gamma\to 0$) the series (6)
truncated at $N_{\rm max}=4$, reproduces with 1% accuracy the GR QNMs for BH’s
spins $\chi_{\rm f}\lesssim 0.7$. The values of the fitting coefficients
$\omega_{\ell mn}^{(j)}$ and $\tau_{\ell mn}^{(j)}$ can be found in Ref.
Maselli  _et al._ (2020). In Ref. Carullo (2021), the fitting coefficients
were calculated up to $N_{\rm max}=9$, which extend the validity of the spin-
expansion up to $\chi_{\rm f}\lesssim 0.99$. As we will discuss in Sec. IV,
the expansion to $N_{\rm max}=4$ is sufficient for our purposes. For
convenience we list the coefficients in the case of GR in Table 2.

###  III.2 The parametrized waveform model

We now describe the waveform model used in our paper to infer properties of a
BBH ringdown. As in Refs. Brito  _et al._ (2018); Ghosh  _et al._ (2021), we
use an IMR BBH waveform model where the complex-valued frequencies describing
the remnant object are left additionally free and estimated directly from the
data.

The GW signal from a quasicircular BBH can be described in GR by a unique set
of parameters $\bm{\theta}$, that includes the masses and spins of the two
BHs, $(m_{1},\,m_{2},\,\mathbf{S}_{1},\,\mathbf{S}_{2})$, the sky location
determined by the luminosity distance $D_{L}$, right ascension $\alpha$ and
declination $\delta$, and the orientation of the binary given by the
inclination $\iota$ and polarization $\psi$ angles. The set is completed by
the choice of a reference time $t_{0}$ and phase $\phi_{0}$. If we further
assume that the spins of the individual BHs are restricted to be aligned or
anti-aligned (for short, aligned) to the unit vector perpendicular to the
orbital plane ($\mathbf{\hat{L}}$), we reduce the six components of the spins
to just two, $\chi_{i}\equiv\mathbf{S}_{i}\cdot\mathbf{\hat{L}}/m_{i}^{2}$
with $i=1,2$, and our entire parameter set from 15 to 11. Let us also define
some additional parameters and set some conventions that will be useful in our
analysis later, namely, the total mass $M=m_{1}+m_{2}$, the chirp mass
$\mathcal{M}=(m_{1}m_{2})^{3/5}/(m_{1}+m_{2})^{1/5}$, the asymmetric mass
ratio $q=m_{1}/m_{2}$, with the convention $m_{1}\geqslant m_{2}$ (and thus
$q\geqslant 1$), and the symmetric mass ratio of the binary,
$\nu=m_{1}m_{2}/(m_{1}+m_{2})^{2}$. Note that for BHs
$-1\leqslant\chi_{i}\leqslant 1$.

For the polarizations of the GW signal (in the observer’s frame) we have

$h_{+}(\iota,\varphi_{0};t)-ih_{\times}(\iota,\varphi_{0};t)=\sum_{l,m}{}_{-\\!2}Y_{lm}(\iota,\varphi_{0})\,h_{lm}(t),$
(8)

where ${}_{-\\!2}Y_{lm}(\iota,\varphi_{0})$ are the $-2$ spin-weighted
spherical harmonics. As our baseline model, that is, the GR model upon which
non-GR modifications are added, we use the computationally efficient (time-
domain) multipolar waveform model for quasicircular spin-aligned BBHs
described in Ref. Mihaylov  _et al._ (2021), which contains the modes,
$(l,|m|)=(2,2),(2,1)$, $(3,3)$, $(4,4)$, and $(5,5)$. Such a model was built
by applying the post-adiabatic approximation Nagar and Rettegno (2019) to the
multipolar spin-aligned EOB waveform model of Refs. Bohé  _et al._ (2017);
Cotesta  _et al._ (2018) (henceforth we refer to our baseline model as
SEOBNR555This waveform model is available in LALSuite LIGO Scientific
Collaboration (2018) as the SEOBNRv4HM_PA waveform approximant.).

An accurate description of the merger is incorporated through calibration with
NR simulations, as described in Refs. Bohé  _et al._ (2017); Cotesta  _et al._
(2018), along with information for the merger and ringdown phases, from BH
perturbation theory. The merger-ringdown waveform, $h_{lm}^{\textrm{merger-
RD}}$, is then stitched to inspiral-plunge waveform, $h_{lm}^{\textrm{insp-
plunge}}$ at a certain time $t=t^{\textrm{match}}_{lm}$, as

$h_{lm}(t)=h_{lm}^{\textrm{insp-
plunge}}\,\Theta(t_{\textrm{match}}^{lm}-t)+h_{lm}^{\textrm{merger-
RD}}\,\Theta(t-t_{\textrm{match}}^{lm})\,,$ (9)

where $\Theta(t)$ is the Heaviside step function. The merger-ringdown waveform
is expressed as an exponentially damped sinusoid Bohé  _et al._ (2017);
Cotesta  _et al._ (2018)

$h_{lm}^{\textrm{merger-RD}}(t)=\nu\ \tilde{A}_{lm}(t)\
e^{i\tilde{\phi}_{lm}(t)}\ e^{-i\sigma_{lm0}(t-t^{\textrm{match}}_{lm})}\,,$
(10)

where

$\displaystyle\sigma_{lm0}\equiv{\rm Re}(\sigma_{lm0})+i\,{\rm
Im}(\sigma_{lm0})=\omega_{lm0}-\frac{i}{\tau_{lm0}}\,,$ (11)

are the complex frequencies of the fundamental ($0$-th overtone) QNMs of the
remnant BH. The functions $\tilde{A}_{lm}(t)$ and $\tilde{\phi}_{lm}(t)$ are
defined in Ref. Bohé  _et al._ (2017); Cotesta  _et al._ (2018).

In the SEOBNR model Bohé  _et al._ (2017); Cotesta  _et al._ (2018), the
complex frequencies $\sigma_{\ell m0}$ are computed by first determining the
final mass and spin from estimates of the initial masses and spins through NR-
fitting-formulas Taracchini  _et al._ (2014); Hofmann  _et al._ (2016), and
then converting them to the complex frequencies using BH perturbation theory-
inspired analytical fits Berti  _et al._ (2006, 2009). Hence,

$\displaystyle\omega_{\ell m0}^{\rm GR}$ $\displaystyle\equiv\omega^{\rm
GR}_{\ell m0}(m_{1},m_{2},\chi_{1},\chi_{2})\,,$ (12a)
$\displaystyle\tau_{\ell m0}^{\rm GR}$ $\displaystyle\equiv\tau^{\rm GR}_{\ell
m0}(m_{1},m_{2},\chi_{1},\chi_{2})\,,$ (12b)

where $(\omega^{\rm GR}_{\ell m0},\tau^{\rm GR}_{\ell m0})$ refer to the GR
QNM predictions in the baseline SEOBNR model. In this paper, we replace these
GR predictions with QNM frequencies defined through the ParSpec framework
introduced in Sec. III.1 (see Eqs. (6)). Hence,

$\displaystyle\omega_{\ell m0}$ $\displaystyle\equiv\omega_{\ell
m0}(m_{1},m_{2},\chi_{1},\chi_{2},\ell_{\rm th},\\{\delta\omega_{\ell
m0}^{(j)}\\}),$ (13a) $\displaystyle\tau_{\ell m0}$
$\displaystyle\equiv\tau_{\ell m0}(m_{1},m_{2},\chi_{1},\chi_{2},\ell_{\rm
th},\\{\delta\tau_{\ell m0}^{(j)}\\}),$ (13b)

where the mass and spin of the remnant object $(M_{\rm f},\chi_{\rm f})$ are
themselves functions of $(m_{1},\,m_{2},\,\chi_{1},\,\chi_{2})$ Taracchini
_et al._ (2014); Hofmann  _et al._ (2016), and we fix $p$ to a certain theory-
specific value. Additionally, the frequencies depend on the ParSpec
coefficients $\\{\delta\omega_{\ell m0}^{(j)},\delta\tau_{\ell m0}^{(j)}\\}$,
and the lengthscale $\ell_{\rm th}$.

Using this parametrized waveform model, which we call pSEOBNR, we infer bounds
on our non-GR parameter $\ell_{\rm th}$ for the specific cases of modified
gravity theories presented Sec. II. We detail our results in Sec. V.

###  III.3 From theory-independent to theory-specific QNM results

Let us now establish the connection between the theory-independent framework
of the pSEOBNR waveform model and the theory-specific QNM calculations. In
this paper we restrict ourselves to the leading and next-to-leading order
terms in the ParSpec expansion, as well as to the fundamental QNM
$(\ell,m,n)=(2,\,2,\,0)$. For this reason, for simplicity, we omit the
subscripts hereafter and rewrite $\omega_{\ell mn}$ and $\tau_{\ell mn}$,
given by Eqs. (6) as,

$\displaystyle M_{\rm f}\,\omega$
$\displaystyle=\gamma\left[\delta\omega^{(0)}\omega^{(0)}+\chi_{\rm
f}\,\delta\omega^{(1)}\omega^{(1)}\right]+\sum_{j=0}^{N_{\rm
max}}\chi^{j}_{\rm f}\,\omega^{(j)}\,,$ (14a) $\displaystyle\frac{\tau}{M_{\rm
f}}$ $\displaystyle=\gamma\left[\delta\tau^{(0)}\tau^{(0)}+\chi_{\rm
f}\,\delta\tau^{(1)}\tau^{(1)}\right]+\sum_{j=0}^{N_{\rm max}}\chi^{n}_{\rm
f}\,\tau^{(j)}\,,$ (14b)

where we pull out from the sum all non-GR corrections, restricting ourselves
to the nonspinning ($j=0$) and linear-order in spin ($j=1$) corrections to the
GR QNMs. The QNMs associated to the higher ($\ell,|m|$)-modes listed in Sec.
III.2 are kept with their GR values.

| GR Maselli  _et al._ (2020)
---|---
$\omega^{(0)}$ | 0.3737
$\tau^{(0)}$ | 11.2407
$\omega^{(1)}$ | 0.1258
$\tau^{(1)}$ | 0.2522

| EdGB ($p=4$) Blázquez-Salcedo  _et al._ (2016); Pierini and Gualtieri (2021) | dCS ($p=4$) Wagle  _et al._ (2022) | cubic EFT ($p=4$) Cano  _et al._ (2022) | quartic EFT ($p=6$) Cano  _et al._ (2022)
---|---|---|---|---
$\delta\omega^{(0)}$ | 0.0107 | 3.1964 |  $-$0.5813 |  $-$0.2114
$\delta\tau^{(0)}$ | 0.0044 | 6.3619 |  $-$0.2114 |  $-$0.6070
$\delta\omega^{(1)}$ |  $-$0.2480 | 41.199 | 6.4439 |  $-$1.5263
$\delta\tau^{(1)}$ |  $-$1.1014 | 794.66 | 265.12 | 171.35

Table 2: Summary of theory-specific QNM calculations. We summarize each theory
we consider together with: the exponent $p$ at which their QNM-modification
enters, the corresponding modifications to the oscillation frequency
$\delta\omega^{(n)}_{220}$ and decay time $\delta\tau^{(n)}_{220}$, and the
references from which we used the results from. We also include for comparison
the GR coefficients, up to $j=1$, obtained in Ref. Maselli  _et al._ (2020).
The remaining GR coefficients for $1<j\leqslant 4$, for which their non-GR
counterpart cannot be determined as of yet for the theories under
consideration, can be found in Table I of Ref. Maselli  _et al._ (2020).

How can we determine the beyond-GR corrections? In GR, comparison between the
numerically determined Kerr QNMs against the fitting formula (14) fixes the GR
expansion coefficients $\omega^{(j)}$ and $\tau^{(j)}$. We can proceed in a
similar way with QNMs calculated in the context of a non-GR theory. In
particular, in the literature, we can already find fitting formulas relating
the QNMs to the BHs mass, spin and lengthscale $\ell_{\rm th}$, the latter
being specific to each theory, up to $j=1$ in the spin expansion (see Table
2). The idea is then to compare these formulas against Eq. (14) to fix $p$,
$\delta\omega^{(j)}$, and $\delta\tau^{(j)}$. Because QNMs of rotating BHs in
modified gravity theories are not known to all spin values, we can expect that
the $j=1$ coefficients to change as calculations beyond-leading order in spin
are accomplished in the future. That is not the case for the $j=0$
coefficients and the situation is the same as in GR, in which the $j=0$
coefficients are simply the QNMs of the Schwarzschild BH.

In the end, the pSEOBNR waveform model with theory-specific QNMs has only
$\ell_{\rm th}$ as a free beyond-GR parameter. We emphasize that our procedure
is different from that of Ref. Carullo (2021) which, for a given value of $p$,
varied all $\ell_{\rm th}$, $\delta\omega^{(j)}$, and $\delta\tau^{(j)}$
parameters, and then proceeded to use the posteriors on $\ell_{\rm th}$,
considering up to $j=2$ in the GR deformation coefficients, and remaining
agnostic about the underlying theory which would predict the modifications to
the QNMs. We will see in Sec. V that adding theory-specific information to the
ParSpec coefficients _can lead to different interpretations of the bounds on
$\ell_{\rm th}$, even for different theories that predict the same value of
the exponent $p$._

As we have seen in Sec. II, QNMs of slowly rotating BHs in modified gravity
theories can belong to two families depending on how they behave under a
parity transformation: axial and polar. Which one do we use to match with Eqs.
(14)? To answer this question one has to work with a chosen theory and perform
a translation between the metric perturbations $h_{\mu\nu}$ in the Regge-
Wheeler-Zerilli gauge Regge and Wheeler (1957); Zerilli (1970) and connect it
with the transverse-traceless gauge used to described GWs (see, for instance,
Ref. Maggiore (2018), Chapter 12), In GR, both axial and polar QNMs are
isospectral and hence which QNM we use to model the ringdown makes no
difference. In beyond-GR theories, isospectrality is in general broken (see in
Ref. Hui  _et al._ (2021) for a counterexample). Thus, how axial and polar
gravitational QNMs appear in the GW signal has to be answered on a theory-by-
theory basis. This is outside the scope of this paper and here we take the
more pragmatic approach of simply choosing the _least damped gravitational
mode between the two parities._ Underlying this choice, are the assumptions
that either (i) the least-damped QNM is also the one excited with largest
amplitude or (ii) that QNMs of both parities are excited with comparable
amplitudes, and one of the modes decays sufficiently fast to not appear in the
ringdown. We performed the mapping between theory-specific QNM calculations
and the ParSpec framework under the hypothesis above, for the theories listed
in Sec. II. We summarize our results in Table 2 and leave the details of our
calculations to Appendix A.

Figure 1: Gravitational-wave signal with GW150914-like parameters for both GR
(solid line) and cubic EFT of GR (dashed line) with the leading-order $n=0$
modifications to the fundamental QNM, for $\ell_{\rm cEFT}=65$ km. The former
is computed with the SEOBNR model, while the latter with the pSEOBNR model,
with ringdown modifications according to the results in Table 2. Top panel:
the $+$ polarization $h_{+}(t)$. Bottom panel: the GW amplitude
$|\tilde{h}(t)|$ (left axes) and the instantaneous frequency $f(t)$ (right
axes).

In Fig. 1 we show an illustrative waveform for GR (solid line; using the
SEOBNR model) and in the cubic EFT of gravity (dashed line; using the pSEOBNR
model), including the leading-order $j=0$ deformations to the fundamental QNM
(see Table 2). We choose binary parameters similar to GW150914, but
nonspinning, with detector-frame masses $m_{1}=39~{}{\rm M}_{\odot}$ and
$m_{2}=31~{}{\rm M}_{\odot}$. The top panel shows the $+$ GW polarization
$h_{+}$ in both theories, while the bottom panel shows the amplitude
$|h|=(h_{+}^{2}+h_{\times}^{2})^{1/2}$ and instantaneous frequency $f$, all as
functions of time $t$. The signals are identical up to the merger,
specifically, the time $t_{\textrm{match}}$ defined in Sec. III.2, after which
they differ during the ringdown. By construction, the ringdown lasts longer
for the cubic EFT of GR waveform and with a smaller instantaneous frequency
(see the bottom panel) due to the negative value of the $\delta\omega^{(0)}$
coefficient in this theory.

##  IV Parameter inference and validity of our bounds

In this section, we first provide a basic outline of the Bayesian formalism
that we use to infer the properties of the underlying GW signal; then, we
identify the most promising events from the catalog of LVK GW observations to
base our analyses on. Finally, we discuss how we can interpret our results
after taking into account the region of validity of the non-GR theories that
we are considering, which are EFTs.

###  IV.1 Bayesian formalism

If we assume that a GW signal observed in detector data $d$ is accurately
described by our waveform model pSEOBNR, we can infer the parameters of the
model, $\bm{\lambda}$, given the hypothesis $\mathcal{H}$, using Bayes’
theorem,

$P(\bm{\lambda}|d,\mathcal{H})=\frac{p(\bm{\lambda}|\mathcal{H})\,\mathcal{L}(d|\bm{\lambda},\mathcal{H})}{E(d|\mathcal{H})}\,,$
(15)

where $P(\bm{\lambda}|d,\mathcal{H})$ is the posterior probability
distribution, $p(\bm{\lambda}|\mathcal{H})$ the prior,
$\mathcal{L}(d|\bm{\lambda},\mathcal{H})$ the likelihood, and
$E(d|\mathcal{H})$ the evidence. The set of parameters, $\bm{\lambda}$ is a
union of the GR waveform model parameters $\bm{\theta}$ (see Sec. III.2) and
$\ell_{\rm th}$, the only non-GR parameter in this problem which, we recall,
sets the characteristic lengthscale in which deviations from GR become
relevant in each of the theories described in Sec. II.

Assuming stationary Gaussian noise, we can write the (log) likelihood function
as,

$\ln\mathcal{L}(d|\bm{\lambda},\mathcal{H})\propto-\tfrac{1}{2}\langle
d-h(\bm{\lambda})|d-h(\bm{\lambda})\rangle\,,$ (16)

with the noise-weighted inner product $\langle\cdot|\cdot\rangle$ defined as,

$\langle A|B\rangle=2\int_{f_{\rm low}}^{f_{\rm high}}\,{\rm
d}f\,\frac{\tilde{A}^{\ast}(f)\tilde{B}(f)+\tilde{A}(f)\tilde{B}^{\ast}(f)}{S_{n}(f)}\,,$
(17)

where $\tilde{A}(f)$ is the Fourier transform of $A(t)$, the asterisk denotes
complex conjugation, $S_{n}(f)$ is the power spectrum density of the detector,
and $[f_{\rm low},f_{\rm high}]$ span the detector sensitivity frequency band.
Assuming a specific prior distribution for our parameters (discussed further
in the next section), we stochastically sample over the parameter space using
a Markov-Chain Monte Carlo algorithm as implemented in LALInferenceMCMC Rover
_et al._ (2006); van der Sluys  _et al._ (2008), as part of the LALInference
software suite Veitch  _et al._ (2015); LIGO Scientific Collaboration (2018).
We subsequently marginalize over the remaining parameters to obtain the
posterior probability distribution function (PDF) on $\ell_{\rm th}$ i.e.,
$P(\ell_{\rm th}|d,\mathcal{H})$, our main parameter of interest.

For $N$ independent GW observations $\\{d_{j}\\}$, $j=1,...,N$, each
characterized by a PDF $P_{j}(\ell_{\rm th}|d_{j},\mathcal{H})$, the joint
posterior can be written as:

$P(\ell_{\rm th}|\\{d_{j}\\},\mathcal{H})=p(\ell_{\rm
th})\prod_{j=1}^{N}\frac{P_{j}(\ell_{\rm
th}|d_{j},\mathcal{H})}{p_{j}(\ell_{\rm th}|\mathcal{H})}\,.$ (18)

where $p_{j}(\ell_{\rm th}|\mathcal{H})$ are the priors used for each
observation, $p(\ell_{\rm th})$ is an overall prior, and we assume that the
value of $\ell_{\rm th}$ is shared among all events. Since we assume a uniform
prior on $\ell_{\rm th}$, the joint posterior is equal to the joint
likelihood. Hereafter, we will drop the explicit usage of $\mathcal{H}$.

###  IV.2 Priors

The prior distribution functions on the GR parameters are assumed to be
uniform over the component masses, $(m_{1},m_{2})$, isotropically distributed
on a sphere in the sky for the source location with $p(D_{L})\propto
D_{L}^{2}$, and isotropic on the binary orientation,
$p(\iota,\psi,\phi_{0})\propto\sin\iota$. For the spins $(\chi_{1},\chi_{2})$,
we assume a prior uniform and isotropic in the spin magnitudes666This spin-
prior choice can be specified in LALInference using the option alignedspin-
zprior..

Among our non-GR parameters $\\{\ell_{\rm
th},\delta\omega^{(j)}\delta\tau^{(j)}\\}$, as already mentioned in the
previous section, we hold $\\{\delta\omega^{(j)},\delta\tau^{(j)}\\}$ fixed to
theory-specific predictions, and only allow $\ell_{\rm th}$ to vary freely. We
assume a uniform prior on $\ell_{\rm th}$, which ranges between $\ell_{\rm
th}=0$ km and $\ell_{\rm th}\sim 100-300$ km, the specific value chosen to
ensure that the marginalized posterior distributions on this parameter do not
rail against the prior’s maximum value. The lower limit is set by the fact
that the modified gravity theories we consider all have $p$ even and hence we
can assume $\ell_{\rm th}>0$ without loss of generality.

###  IV.3 Events selection

The pSEOBNR model, as described in Sec. III.2, is an IMR model that infers the
properties of the underlying GW signal, including (independently) its ringdown
properties, using the Bayesian formalism above. Naturally, the most promising
candidates for our analyses are high-mass _and_ loud GW observations with a
significant signal-to-noise ratio (SNR) in the inspiral and post-merger stages
to break the degeneracy between the total mass and the QNM frequencies. The
latest LVK GW catalog Abbott  _et al._ (2021a) reported 90 observed signals
not all of which are relevant for our BH ringdown analysis. In fact, in the
accompanying paper Abbott  _et al._ (2021e) on tests of GR, the pSEOBNRv4HM
Brito  _et al._ (2018); Ghosh  _et al._ (2021) analysis 777See, in particular,
Sec.VIII A.2 in Ref. Abbott  _et al._ (2021e), which is most similar to the
pSEOBNR model presented in this paper, identified two events which provided
the strongest bounds on the measurements of the dominant $(220)$ QNM: GW150914
Abbott  _et al._ (2016a) and GW200129 Abbott  _et al._ (2021a). These two
events, with a total (source-frame) mass of $65M_{\odot}$ and $63.4M_{\odot}$
respectively, are extremely similar in their source properties. These are also
two of the loudest BBH signals observed to date with a total network SNR of 24
and 26.8, respectively. Moreover, and what is more relevant for our analysis,
are their post-inspiral (merger-ringdown) SNRs which are both $\approx 16$
(see the columns for $\rho_{\text{post-insp}}$ in Table III of Ref. Abbott
_et al._ (2019b) and Table IV of Ref. Abbott  _et al._ (2021e)). In this
paper, we are going to focus on these two GW events as our probes of the BH
ringdown in modified theories of gravity.

The parameter inference in this paper follows configurations identical to the
ones used on these events for the pSEOBNRv4HM analysis in Ref. Abbott  _et
al._ (2021e). GW150914 was a 2-detector (Hanford-Livingston) event while
GW200129 was 3-detector (Hanford-Livingston-Virgo). We consequently use the
same strain data $h(t)$, detector power-spectral-densities $S_{n}(f)$ and
calibration envelopes as were used for the analyses in Ref. Abbott  _et al._
(2021e).

In Sec. V, we enumerate through the different theories and outline the main
results. Whenever possible, we also combine results from both events to obtain
the strongest possible bound on $\ell_{\rm th}$.

###  IV.4 EFT interpretation of our results

There are two conditions that we must verify before we can confidently claim
to have placed a constraint on $\ell_{\rm th}$. First, as we have explained in
Sec. II, all theories that we consider must be interpreted as an EFT, meaning
that they should be considered valid only below an energy scale, or
equivalently, above a lengthscale. As a cutoff lengthscale for the validity of
the EFT we use,

$\Lambda_{\rm
EFT}(\varepsilon,\mathfrak{m})=\varepsilon\,\frac{G\mathfrak{m}}{c^{2}}\,,$
(19)

where $\varepsilon$ is a dimensionless number and $\mathfrak{m}$ is the median
value of one of the mass scales involved in the problem. We note that
$\Lambda_{\rm EFT}$ has dimensions of length and hence can be compared to each
theory’s fundamental lengthscale $\ell_{\rm th}$. Here we explore the range
$\varepsilon\in[0,1]$, but following Refs. Nair  _et al._ (2019); Perkins  _et
al._ (2021); Lyu  _et al._ (2022) we quote our final results using
$\varepsilon=1/2$, but we stress that there is no fundamental justification
for this choice.

Under these assumptions, we will say that a bound has been placed on
$\ell_{\rm th}$, if most of the PDF $P(\ell_{\rm th}|d)$ support is in the
interval $[0,\,\Lambda_{\rm EFT}(1/2,\mathfrak{m})]$. In practice, this can be
quantified through the cumulative distribution function (CDF) associated with
the marginalized posterior distribution $P(\ell_{\rm th}|d)$, namely

$P(\ell_{\rm th}\leqslant\ell^{\rm max}_{\rm th}|d)=\int_{0}^{\,\ell^{\rm
max}_{\rm th}}{\rm d}\ell^{\prime}\,P(\ell^{\prime}|d).$ (20)

For instance, we require that for a bound at 90% credible level to be placed
on $\ell_{\rm th}$ that

$P(\ell_{\rm th}\leqslant\Lambda_{\rm EFT}|d)\geqslant 0.9,\quad\textrm{(EFT
bound)},$ (21)

where we let $\ell^{\rm max}_{\rm th}=\Lambda_{\rm EFT}$ in Eq. (20), and
likewise for other credibility percentiles.

Second, as already emphasized in Ref. Maselli  _et al._ (2020), the ParSpec
formalism is by construction perturbative. This means that the non-GR
deformation parameters are small, that is,

$\gamma\,\delta\omega^{(j)}\ll
1,\,\,\textrm{and}\,\,\gamma\,\delta\tau^{(j)}\ll 1,\quad\textrm{(ParSpec
bound)},$ (22)

for all orders $j$ in the expansion in dimensionless spin $\chi_{\rm f}$ and
where $\gamma$ was defined in Eq. (7). We also construct posterior
distributions for these parameters and check if most of their support is
concentrated to a domain with values much smaller than unity.

Another question we must consider is the following: what is the mass
$\mathfrak{m}$ that we should use in Eq. (19)? In Refs. Nair  _et al._ (2019);
Perkins  _et al._ (2021); Lyu  _et al._ (2022), which attempted to constrain
dCS and EdGB theories with the _inspiral_ part of the GW signal alone, it was
natural to choose the secondary’s mass $m_{2}$ as the most conservative
choice, since it is by definition the smaller component mass and hence places
the lowest cutoff scale $\Lambda_{\rm EFT}$ for the validity of either of
these theories as an EFT.

In our problem, the answer is not as clear. On the one hand, since we are
interested in the ringdown part of the signal, it is natural to use the final
mass $M_{\rm f}$ to compute $\Lambda_{\rm EFT}$. On the other hand, one may
argue that the modified gravity theory under consideration should be able to
predict a full inspiral, merger, and ringdown of the BBH before we can even
make such a test, and thus the same, more conservative choice
$\mathfrak{m}=m_{2}$ should be used. Here we adopt a pragmatic approach to
this issue and consider _both_ masses, $m_{2}$ and $M_{\rm f}$, to determine
$\Lambda_{\rm EFT}$. Specifically, we will use the median value of the
marginalized PDF of these masses. We then compare how different assumptions
yield to different interpretations of the results of our parameter estimation.

##  V Results using LIGO-Virgo events

###  V.1 Einstein-dilaton-Gauss-Bonnet gravity

We start with EdGB gravity. In Fig. 2 we show the marginalized PDFs of the
coupling constant $\ell_{\rm GB}$, for GW150914 (top panel) and GW200129
(middle panel), with and without the spin corrections to the $(2,2,0)$ QNM.
The bottom panel shows the joint posterior obtained by combining both events.
We see, for both events, the $N_{\rm max}=0$ posteriors are characterized by a
peak away from zero. This does not mean that we are inferring a deviation from
GR. We recall that the deviations from GR in the ParSpec framework are
controlled by the dimensionless parameter $\gamma$, which here reads,

$\gamma_{\rm EdGB}=\left(\frac{c^{2}\ell_{\rm EdGB}}{GM^{\rm s}_{\rm
f}}\right)^{4}\,.$ (23)

As shown in Fig. 3, $\gamma_{\rm EdGB}$ does indeed have a posterior
distribution with largest support at zero, indicating consistency between the
underlying signal and GR. We also observe that the inclusion of the spin
corrections (i.e., the curves with $N_{\rm max}=1$) displaces the posteriors
distributions towards smaller values of $\ell_{\rm EdGB}$ (see Fig. 2), and
larger values of $\gamma_{\rm EdGB}$ (see Fig. 3).

Figure 2: Posterior distribution function on coupling constant $\ell_{\rm
EdGB}$ in EdGB gravity for GW150914 (top panel), GW200129 (middle panel) and
from combining results (bottom panel). In all panels, different line colors
correspond to the inclusion ($N_{\rm max}=1$) or not ($N_{\rm max}=0$) of the
linear-in-spin QNM correction. The joint posteriors are shown for illustrative
purposes only. As we explain in Fig. 4 and in the main text, our analysis of
these events fails to satisfy the EFT bound (21). We mark with solid vertical
lines the $90\%$ upper credible intervals, while the dashed and dot-dashed
lines correspond to the EFT bound, $\Lambda_{\rm EFT}(0.5,m_{2})$ and
$\Lambda_{\rm EFT}(0.5,M_{\rm f})$, respectively. We see that the $90\%$ upper
credible level lines are located to right of the EFT cutoffs. This means that
the EFT bound is not satisfied.  Figure 3: Posterior distribution function on
the dimensionless parameter $\gamma_{\rm EdGB}$, defined in Eq. (23).
Different line colors distinguish between events, while different line styles
distinguish between different $N_{\rm max}$. We see that this parameter which
controls the ParSpec expansion in EdGB gravity does have maximum support at
$\gamma_{\rm EdGB}=0$. This shows that our model is consistent with GR. The
right panel shows the PDFs in the range $0\leqslant\gamma_{\rm EdGB}\leqslant
1$. Note how the curves are flat in this range for both events and when
$N_{\rm max}=1$. Figure 4: The CDF evaluated at the cutoff $\Lambda_{\rm
EFT}(\varepsilon,\mathfrak{m})$ for EdGB gravity as a function of the
parameter $\varepsilon$ for both (dashed curves) $\mathfrak{m}=m_{2}$, the
secondary’s source mass, and (solid curves) $\mathfrak{m}=M_{\rm f}$, the
remnant’s source mass, without (black curves) and with (purple curves) linear
in spin QNM corrections. The horizontal lines mark the $90\%$ credible levels.
Having set the maximum value of $\varepsilon$ to be 1/2, we see that in no
situation the curves pass through the 0.90 lines. This means that no bound on
$\ell_{\rm EdGB}$ can be placed with the events we analyzed.

As we have emphasized in Sec. IV.4, we must first check whether the “EFT” (21)
and “ParSpec” (22) bounds are satisfied, before drawing any conclusions on the
allowed values for $\ell_{\rm EdGB}$ from our parameter estimations. We check
the validity of the EFT bound in Fig. 4. In the top (bottom) panel we show the
CDF of the $\ell_{\rm EdGB}$ posteriors for GW150914 (GW200129), obtained by
evaluating the integral (20) with $\ell^{\rm max}_{\rm th}=\Lambda_{\rm
EFT}(\varepsilon,\mathfrak{m})$, with the mass scale set by the secondary’s
mass (i.e., $\mathfrak{m}=m_{2}$, dashed lines) or the remnant’s mass
($\mathfrak{m}=M_{\rm f}$, solid lines), while varying $\varepsilon$ between 0
and 1. For GW150914, we see that for the $N_{\rm max}=0$ curves, that the CDF
never goes past $0.2$, regardless of the mass scale $\mathfrak{m}$ used and
even at $\varepsilon=1$, at which the EFT description of the theory would not
be valid anyway. This shows that the “EFT bound” given by Eq. (21) is never
met to a significant credible level and thus that we cannot place a bound on
$\ell_{\rm EdGB}$. The situation is similar for GW200129 with $N_{\rm max}=0$
and does not change for either event when we add spin corrections to the EdGB
QNM. For the case with $N_{\rm max}=1$, we find that the “EFT bound” is
satisfied only for $\varepsilon\approx 0.8$ and $\approx 1$ for GW150914 and
GW200129, respectively. However, we set the maximum value of $\varepsilon$ to
be 1/2, thus, taken together we are led to _conclude that we cannot constrain
EdGB gravity with our present model_. We summarize our findings in Table 3.

$N_{\rm max}$ | Event | EFT | ParSpec | Constraint
---|---|---|---|---
|  | bound? | bound? | ($\mathfrak{m}=M_{f}$)
| GW150914 | No | Yes | –
0 | GW200129 | No | Yes | –
| Combined | – | Yes | –
| GW150914 | No | Yes | –
1 | GW200129 | No | Yes | –
| Combined | – | Yes | –
Table 3: Detailed summary of our results for EdGB gravity for GW150914,
GW200129, and combined events using $\mathfrak{m}=M_{\rm f}$,
$\varepsilon=1/2$ and quoting only 90% credible results. We find that we
cannot place any bound on $\ell_{\rm EdGB}$ with our waveform model from
either GW event.

We can compare this conclusion with that of Ref. Carullo (2021), which found
that $p=4$ modifications (such as the case of EdGB gravity) are constrained to
$\ell\lesssim 35$ km, _but not_ including theory-specific QNM information on
$\delta\omega^{(j)}$ and $\delta\tau^{(j)}$. Furthermore, Ref. Carullo (2021)
did not impose the EFT bound that we imposed. Our results provide a concrete
example of the importance of including theory-specific QNM calculations
information into the parameter estimation and how this can dramatically change
the outcome of the results.

Let us also contrast our results with those of Refs. Nair  _et al._ (2019);
Perkins  _et al._ (2021); Lyu  _et al._ (2022) which relied on the BBH
inspiral to constrain $\ell_{\rm EdGB}$, as discussed in Sec. II.1. We see
that EdGB gravity provides an example of a theory in which, with current GW
events, the inspiral portion of the signal can be more constraining than the
ringdown portion of the signal. Two reasons together can explain our negative
results. First, as observed by Ref. Blázquez-Salcedo  _et al._ (2016), the
QNMs of EdGB BHs only differ slightly from their Schwarzschild counterparts.
Second, the larger mass $M_{\rm f}$ of the remnant BH, suppresses scalar
field’s charge relatively to the initial binary components.

###  V.2 Dynamical Chern-Simons gravity

Figure 5: Similar to Fig. 2, but for dCS gravity. We stress that posteriors
obtained including $n=1$ corrections, violate conditions (22) and therefore
should not be used to draw meaningful conclusions. We show them for
illustrative purposes and also to emphasize the importance of taking
conditions (21) and (22) simultaneously into consideration when analyzing the
results of the parameter estimation.  Figure 6: Similar to Fig. 4, but for dCS
gravity. We see that the CDF curves for $\mathfrak{m}=M_{\rm f}$ are above 90%
for $\varepsilon=1/2$ for both events, with and without including the $n=1$
dCS corrections to the dominant QNM.

We now consider dCS gravity where the main results are summarized in Fig. 5.
As with EdGB gravity (see Sec. V.1) although the PDF of $\ell_{\rm dCS}$ is
peaked away from 0, this does not signify a deviation from GR, as we have
verified that,

$\gamma_{\rm CS}=\left(\frac{c^{2}\ell_{\rm dCS}}{GM^{\rm s}_{\rm
f}}\right)^{4}\,,$ (24)

does indeed peak at zero indicating consistency with GR, similarly to what is
shown in Fig. 3 for $\gamma_{\rm EdGB}$. We also see that in both cases the
inclusion of leading-order–in-spin correction to the QNM displaces the
posteriors toward smaller values of $\ell_{\rm dCS}$. This can be seen more
evidently by looking at the location of posterior peaks. Finally, in the
bottom panel, we show the combined result for both events.

In Fig. 6 we show the CDF for GW150914, we see that with $\mathfrak{m}=m_{2}$,
Eq. (21) is not satisfied unless $\varepsilon\approx 0.9$ (with only $j=0$
corrections) and $\varepsilon\approx 0.7$ (with both $j=0$ and $1$
corrections). The situation is different if we use $\mathfrak{m}=M_{\rm f}$.
In this case, we find that with or without spin corrections Eq. (21) can be
satisfied with $\varepsilon\leqslant 1/2$ (i.e., below the criteria used Refs.
Nair  _et al._ (2019); Perkins  _et al._ (2021); Lyu  _et al._ (2022)). This
means that with our model’s assumptions and using the remnant’s source mass
$M_{\rm f}$ to set the cutoff scale that we can claim an upper bound

$\ell_{\rm dCS}\leqslant 41.9~{}\textrm{km}\quad\textrm{at 90\% credible
level},$ (25)

on dCS gravity. This result would constitute _the strongest bound to date on
this theory with GW observations alone, and also the first bound using GW
generation effects._

We can draw qualitatively similar conclusions from the GW200129 event. In
particular, we find,

$\ell_{\rm dCS}\leqslant 35.8~{}\textrm{km}\quad\textrm{at 90\% credible
level}.$ (26)

These stronger bounds are a consequence of the larger support for $\ell_{\rm
dCS}\lessapprox 15$ km for GW200129 (compare the top and middle panels in Fig.
5), and in part due to the smaller median remnant ($M_{\rm f}\approx
59.5~{}{\rm M}_{\odot}$ versus $M_{\rm f}\approx 61.8~{}{\rm M}_{\odot}$ for
GW150914). We also found for both GW events, that the perturbative-conditions
(22) required by the ParSpec is violated for the $\gamma_{\rm
dCS}\,\delta\tau^{(1)}$ coefficient. This means that we cannot use this
posterior to infer any meaningful bound on dCS gravity and that is why we
quoted only the $N_{\rm max}=0$ bound above.

Finally, since both events individually lead to a bound on $\ell_{\rm dCS}$
(assuming a cutoff scale for $\mathfrak{m}=M_{\rm f}$ and $N_{\rm max}=0$), we
can combine the posteriors to obtain the cumulative bound,

$\ell_{\rm dCS}\leqslant 38.7~{}\textrm{km}\quad\textrm{ at 90\% credible
level},$ (27)

which is the main result of this section. This bound is approximately a factor
of four weaker than that placed by Ref. Silva  _et al._ (2021a), but it (i)
relies only on GW observations, and (ii) suggests that a ringdown analysis can
potentially place constraints on theories that, with current GW events, can
evade GR tests using inspiral information alone, such as the case of dCS
gravity Nair  _et al._ (2019); Perkins  _et al._ (2021); Lyu  _et al._ (2022).
In Table 4 we summarize our findings of this section.

Figure 7: Corner plot showing that the inferred final spin $\chi_{\rm f}$, remnant mass $M_{\rm f}$ for GW150914, using the same waveform model, but without (purple contours) and with the non-GR parameters different from zero (blue contours), for dCS gravity and $N_{\rm max}=0$. The contours represent 90% credible levels. We see that the introduction of the non-GR parameters does not bias the inference on the source parameters as required by the ParSpec.  $N_{\rm max}$ | Event | EFT | ParSpec | Constraint ($\mathfrak{m}=M_{f}$)
---|---|---|---|---
|  | bound? | bound? |
| GW150914 | Yes | Yes |  $\ell_{\rm dCS}\leqslant 41.9$ km
0 | GW200129 | Yes | Yes |  $\ell_{\rm dCS}\leqslant 35.8$ km
| Combined |  |  |  $\ell_{\rm dCS}\leqslant 38.7$ km
| GW150914 | Yes | No | –
1 | GW200129 | Yes | No | –
| Combined |  |  | –
Table 4: Detailed summary of our results for dCS gravity for GW150914,
GW200129, and combined events using $\mathfrak{m}=M_{\rm f}$,
$\varepsilon=1/2$ and quoting only 90% credible bounds. We found that while
our posteriors satisfy the condition (21) (with $\varepsilon=1/2$), they do
not obey the condition (22) for $N_{\rm max}=1$. This means that our results
for $N_{\rm max}=0$ are the only ones we can confidently quote. The combined
bound, which is also quoted in Table 2, is $\ell_{\rm dCS}\leqslant 38.7$ km
at 90% credible level.

As an additional check, to verify the robustness of our constraint, we show in
Fig. 7, the final spin $\chi_{\rm f}$ and remnant mass $M_{\rm f}$ for
GW150914 for GR and dCS gravity. We see that our pSEOBNR waveform model does
not introduce substantial changes to the GR estimates on these parameters, as
required by the ParSpec expansion (see discussion in Sec. III.1). In fact, we
observed no bias on the estimation of $M_{\rm f}$ and $\chi_{\rm f}$ for all
theories considered here. For completeness, in Appendix B we also show how all
other intrinsic parameters remain unbiased.

###  V.3 Cubic effective-field-theory of general relativity

We now consider the cubic EFT of GR. In Fig. 8 we show the marginalized
posterior distributions functions of $\ell_{\rm cEFT}$ for GW150914 (top
panel) and GW200129 (middle panel), with different curve colors corresponding
to different $N_{\rm max}$ in the spin expansion. We find that in this theory,
the posterior distributions are mostly uniform for $\ell_{\rm cEFT}\lesssim
40$ km (contrast this with the EdGB and dCS gravity cases in Figs. 2 and 5).
For values $\ell_{\rm cEFT}\gtrsim 40$ km, the posteriors smoothly approach
zero.

In Fig. 9 we show the CDF for both events, calculated in the same way as
already described for the EdGB and dCS theories. We see that curves are very
similar to those of dCS gravity for GW200129 (see bottom panel in Fig. 9).
Moreover, we find that the EFT (21) and ParSpec (22) bounds are satisfied for
both events both when $\mathfrak{m}=M_{\rm f}$, $\varepsilon=1/2$, and $N_{\rm
max}=0$. This allows us to place the combined bound of

$\ell_{\rm cEFT}\leqslant 38.2~{}\textrm{km},\quad\textrm{at 90\% credible
level}\,.$ (28)

As also happened for our study for dCS, the find that, for the cubic EFT, the
ParSpec bound is violated by the $N_{\rm max}=1$ corrections to the QNMs,
meaning that we cannot use this case to draw any meaningful constraint on this
parameter. We summarize our results in Table 5.

$N_{\rm max}$ | Event | EFT | ParSpec | Constraint ($\mathfrak{m}=M_{f}$)
---|---|---|---|---
|  | bound? | bound? |
| GW150914 | Yes | Yes |  $\ell_{\rm cEFT}\leqslant 38.2$ km
0 | GW200129 | Yes | Yes |  $\ell_{\rm cEFT}\leqslant 42.5$ km
| Combined |  |  |  $\ell_{\rm cEFT}\leqslant 38.2$ km
| GW150914 | Yes | No | –
1 | GW200129 | Yes | No | –
| Combined |  |  | –
Table 5: Detailed summary of our results the cubic EFT of GR for GW150914,
GW200129, and combined events using $\mathfrak{m}=M_{\rm f}$,
$\varepsilon=1/2$ and quoting only 90% credible results.  Figure 8: Similar to
Fig. 2, but for the cubic EFT of GR. We show our results for GW150914 (top
panel), GW200129 (middle panel) and combined events (bottom panel). The colors
distinguish different $N_{\rm max}$ in the spin expansion. Once again, the
solid vertical lines mark the 90% upper credible intervals, while the dashed
and dot-dashed lines correspond to the EFT bound, $\Lambda_{\rm
EFT}(0.5,m_{2})$ and $\Lambda_{\rm EFT}(0.5,M_{\rm f})$, respectively.  Figure
9: Similar to Fig. 4, but for the cubic EFT of GR. We see that the CDF curves
for $\mathfrak{m}=M_{\rm f}$ are above 90% for $\varepsilon=1/2$ for both
events, with and without including the $j=1$ corrections to the dominant QNM.

###  V.4 Quartic effective-field-theory of general relativity

Let us now consider the quartic EFT of GR, as our final example. In Fig. 10 we
show the posteriors on $\ell_{\rm qEFT}$ for GW150914 (top panel), GW200129
(middle panel) for $N_{\rm max}=0$, which are qualitatively similar to the
cubic EFT of GR. We find that while the Parspec bound is satisfied, the EFT
bound is only marginally so, As shown in Fig. 11, the 90% credible level is
reached for $\varepsilon\approx 0.58$ (in the case of GW15094) and for
$\varepsilon\approx 0.64$ (in the case of GW200129). Having in mind that the
cut off $\varepsilon=1/2$ is not fundamental, but to keep consistency across
our analysis, our final result

$\ell_{\rm qEFT}\leqslant 51.3~{}\textrm{km},$ (29)

at 90% credible level should be taken lightly. However, we _can_ claim the
validity of the bound above, but at a lower, 68% credible level.

$N_{\rm max}$ | Event | EFT | ParSpec | Constraint ($\mathfrak{m}=M_{\rm f}$)
---|---|---|---|---
|  | bound? | bound? |
| GW150914 | Yes | Yes |  $\ell_{\rm qEFT}\leqslant 51.7$ km
0 | GW200129 | Yes | Yes |  $\ell_{\rm qEFT}\leqslant 54.8$ km
| Combined |  |  |  $\ell_{\rm qEFT}\leqslant 51.3$ km
Table 6: Detailed summary of our results the quartic EFT of GR for GW150914,
GW200129, and combined events using $\mathfrak{m}=M_{\rm f}$,
$\varepsilon=1/2$, and $N_{\rm max}=0$. The quoted result correspond to 90%
credible values _if_ we allow for a more flexible cutoff $\varepsilon\lesssim
0.65$. However, the result is robust for the cut off $\varepsilon=1/2$, at 65%
credible level.

In this theory, we have considered only $N_{\rm max}=0$. We find that the
addition of spin corrections (while maintaining the same prior ranges on
$\ell_{\rm qEFT}$ as used in the $N_{\rm max}=0$ study) can result in
waveforms that can have a ringdown segment larger (sometimes seconds long)
than the inspiral-plunge segment in the detectors’ frequency band, making the
parameter estimation challenging. To overcome this issue we have lowered the
value of $\ell^{\rm max}_{\rm qEFT}$, but by doing so we have obtained
posteriors which were flat, just as our prior, and were thus uninformative.
Hence, we do not quote any results for $N_{\rm max}=1$. Table 6 summarizes our
findings for the quartic EFT of GR.

Figure 10: Similar to Fig. 2, but for the quartic EFT of GR. We show our
results for GW150914 (top panel), GW200129 (middle panel) and combined events
(bottom panel). As in Fig. 2, we mark with the solid vertical lines the 90%
upper credible intervals. Here we calculated the values of $\Lambda_{\rm
EFT}(\varepsilon,\mathfrak{m})$ using $\varepsilon=0.58$ (for GW150914) and
$\varepsilon=0.64$ (for GW200129) as discussed in main text.  Figure 11:
Similar to Fig. 4, but for the quartic EFT of GR. We find that similarly to
what happens in dCS gravity and the cubic EFT of GR, we can place a constraint
on $\ell_{\rm qEFT}$ when $\mathfrak{m}=M_{\rm f}$ for both events at 68%
credibility.

##  VI Conclusions

We presented an unified framework that combines the ParSpec framework to model
deviations to the GR QNMs Maselli  _et al._ (2020) with the pSEOBNR waveform
model Brito  _et al._ (2018); Ghosh  _et al._ (2021). We showed with concrete
examples, how theory-specific QNM calculations of slowly rotating BHs in
modified gravity theories can be mapped onto the non-GR parameters of the
ParSpec formalism. The resulting pSEOBNR waveform model does not bias
(relative to GR) the inference of the intrinsic binary parameters, as required
by ParSpec (see, in particular, Fig. 7 and Fig. 12 in Appendix B), Put
together this allowed us to test four modified gravity theories (EdGB, dCS,
cubic, and quartic EFTs of GR) using observational data from the LVK events
GW150914 and GW200129. Our results are summarized in Table 1.

In particular, we found, that within the interpretation of these theories as
EFTs and the region of validity of the ParSpec framework, the fundamental
lengthscale of dCS gravity is bound as $\ell_{\rm dCS}\leqslant 34.5$ km, at
90% credible level, when stacking the posteriors of GW150914 and GW200129.
This is the strongest constrain to date on this theory with GW observations
alone. In contrast, we could not place any bounds on the fundamental
lengthscale of EdGB gravity $\ell_{\rm EdGB}$. This dichotomy between the two
theories has a counterpart with works that considered the inspiral part of the
GW signal alone Nair  _et al._ (2019); Perkins  _et al._ (2021); Lyu  _et al._
(2022). Using data of the LVK BBHs, it was found that the posterior
distributions for deviations from GR were uninformative in dCS gravity, but
not in EdGB gravity. We emphasize that both those theories (and the cubic EFT
of GR also studied here) all predict the same exponent $p$ in ParSpec. Hence,
our results show how the inclusion of theory-specific information into the
ParSpec framework can result in different outcomes for different theories,
even if they predict the same value of $p$.

Let us discuss some avenues for future work. First, we could implement a high-
spin version of the GR fitting coefficients to the ParSpec formulas. This has
already been done in Ref. Carullo (2021) extending the validity of the ParSpec
formulas up to spins of $\chi_{\rm f}\approx 0.99$. For the events analyzed
here, the original fit by Ref. Maselli  _et al._ (2020) was sufficient, but it
might not be the case with upcoming GW observation campaigns. Second, it would
be important to incorporate additional effects, such as spin-precession and
eccentricity to pSEOBNR (see, for instance, Refs. Ossokine  _et al._ (2020);
Ramos-Buades  _et al._ (2022)). Third, it will be interesting, to perform
tests of modified theories of gravity using IMR waveform models that include,
during the inspiral stage, finite-size effects induced by the non-GR geometry
around the BHs — for example the ones due to spin-induced quadrupole, tidal
deformability and absorption, and also orbital effects due to non-GR
gravitational interactions between the BHs. Those effects could be included
using the flexible theory-independent method Mehta  _et al._ (2022), as done
in Ref. Sennett  _et al._ (2020) or the TIGER code Li  _et al._ (2012);
Agathos  _et al._ (2014). However, setting bounds on deviations from GR caused
by orbital effects requires a different EFT interpretation than what we
adopted in Sec. IV.4 (see also Sec. IIC in Ref. Sennett  _et al._ (2020)).
Indeed, in this case one would need to analyze the data considering that the
modified theory of gravity is valid for $\ell_{\rm th}\gtrsim M$, but
$\ell_{\rm th}\lesssim\mathcal{R}$, being $\mathcal{R}$ the binary’s
separation. Fourth, to test the robustness of the results obtained in this
paper, it will be very useful to employ NR waveforms produced in some of the
non-GR theories under consideration, as synthetic signals, and carry out
Bayesian analysis to recover the binary’s parameters, including the non-GR
ones during the ringdown. As today, there are only a small number of such BBH
NR simulations, for a given theory Healy  _et al._ (2012); Witek  _et al._
(2019); Okounkova  _et al._ (2017); Hirschmann  _et al._ (2018); Okounkova
(2020, 2020); Silva  _et al._ (2021b); East and Ripley (2021a, b); Figueras
and França (2022). Last, but not least, it will be very beneficial to
calculate the QNMs (complex) frequencies of rapidly rotating BHs in modified
gravity theories (with BH perturbation theory). This is a challenging problem,
but certainly necessary, also in the context of ParSpec, if higher spin
corrections would need to be included to make the framework robust.

## Acknowledgments

We thank Emanuele Berti, Andrea Maselli, Caio F. B. Macedo, Deyan Mihaylov,
Serguei Ossokine, Scott E. Perkins, and Helvi Witek for discussions. We also
thank Anuradha Gupta and Gregorio Carullo for comments on this manuscript. We
also thank the referees for the detailed review of this work and for the
suggestions that improved its presentation. We are grateful for the
computational resources provided by the Max Planck Institute for Gravitational
Physics in Potsdam, specifically the high-performace computing cluster Hypatia
and to Steffen Grunewald for assistance. We acknowledge funding from the
Deutsche Forschungsgemeinschaft (DFG) - project number: 386119226. The
material presented in this manuscript is based upon work supported by NSF’s
LIGO Laboratory, which is a major facility fully funded by the NSF. This
research has made use of data, software and/or web tools obtained from the
Gravitational Wave Open Science Center (https://www.gw-openscience.org), a
service of LIGO Laboratory, the LIGO Scientific Collaboration and the Virgo
Collaboration, and Zenodo (https://zenodo.org/record/5172704). The authors
would like to thank everyone at the frontline of the Covid-19 pandemic.

##  Appendix A Details of the determination of the theory-specific ParSpec
coefficients

Here, for the theories described in Sec. II, we use QNM calculations from the
literature and determine the coefficients in the ParSpec, which we have
summarized in Table 2. We consider only the fundamental QNM
$(\ell,m,n)=(2,2,0)$, hence we omit the QNM subscript “$(2,2,0)$” for brevity
and, likewise, the subscript “f” for final BH’s spin and mass.

###  A.1 Einstein-dilaton-Gauss-Bonnet gravity

We start by considering EdGB gravity and focus on Refs. Blázquez-Salcedo  _et
al._ (2016); Pierini and Gualtieri (2021) to determine the ParSpec
coefficients for this theory. In particular, Ref. Blázquez-Salcedo  _et al._
(2016) found that the damping time of the dominant _axial_ gravitational-led
mode increases as the lengthscale $\ell_{\rm EdGB}$ is increased. The leading-
order spin corrections to the _polar_ -parity QNMs was studied in Ref. Pierini
and Gualtieri (2021). Hence, according to the prescription of Sec. III.3, we
select the axial-parity branch of QNMs. For the nonrotating QNMs we use the
numerical data of Ref. Blázquez-Salcedo  _et al._ (2016) and generate a new
linear fit in $\gamma_{\rm EdGB}$ using numerical QNM data valid for small
values of the $\gamma_{\rm EdGB}$ (see, in particular, Eq. (27) and Fig. 1 of
Ref. Blázquez-Salcedo  _et al._ (2016)). We find,

$\displaystyle M{\rm Re}(\sigma)_{\rm EdGB}$ $\displaystyle=M{\rm
Re}(\sigma)_{\rm GR}\,\left(1+0.0107\,\gamma_{\rm EdGB}\right),$ (30a)
$\displaystyle M{\rm Im}(\sigma)_{\rm EdGB}$ $\displaystyle=M{\rm
Im}(\sigma)_{\rm GR}\,\left(1-0.0044\,\gamma_{\rm EdGB}\right).$ (30b)

The small values of the numerical prefactors of $\gamma_{\rm EdGB}$ are a
consequence of the how weakly the QNMs of BHs in EdGB gravity deviate from
their GR counterparts, even at moderately large values of $\gamma_{\rm
EdGB}\approx 0.3$.

The spin-corrections to the _polar_ gravitational-led modes were calculated in
Ref. Pierini and Gualtieri (2021) (see, in particular, their Eqs. (51) and
(52)). For consistency with our previous discussion, we truncate these
equations at leading-order in $\gamma_{\rm EdGB}$, but we emphasize that we
are being inconsistent in mixing results valid for modes of different
parities. We still do so, simply to explore what the rotational corrections to
EdGB gravity QNMs _might_ tells us in our ringdown analysis and the results of
Ref. Pierini and Gualtieri (2021) are our best presently available guide.

We can expand the resulting formula in $\gamma_{\rm EdGB}$ and the
coefficients $\delta\omega^{(i)}$ and $\delta\tau^{(i)}$, $i=1,2$ can be read-
off by comparison against Eqs. (14), where for the damping time we use the
relation ${\rm Im}(\sigma)_{\rm EdGB}=-1/\tau_{\rm EdGB}$ and reexpand in
$\ell_{\rm EdGB}$ and $\chi$. These steps yield for $p_{\rm EdGB}=4$,

$\displaystyle\delta\omega_{\rm EdGB}^{(0)}=0.0107,\quad\delta\tau_{\rm
EdGB}^{(0)}=-0.2480,$ (31)

for the $j=0$ coefficients and

$\delta\omega_{\rm EdGB}^{(1)}=-0.2480,\quad\delta\tau_{\rm
EdGB}^{(1)}=-1.1014.$ (32)

for the $j=1$ coefficients.

###  A.2 Dynamical Chern-Simons gravity

For dCS gravity, we follow Ref. Wagle  _et al._ (2022), which numerically
calculated the QNMs of slowly rotating BHs, and found that for the axial
gravitational-led modes the damping time increases, as we increase $\ell_{\rm
dCS}$, at constant, small BH spin. Hence, according to the prescription of
Sec. III.3, this is the branch of QNMs we choose to work with.

We then proceed to determine $\delta\omega^{(j)}$ and $\delta\tau^{(j)}$ as
follows. Using the fitting formula Eq. (54a) of Ref. Wagle  _et al._ (2022),
namely,

$M{\rm Re}(\sigma)_{\rm
dCS}=c_{1}+c_{2}\kappa\zeta+(c_{3}+c_{4}\kappa\zeta)\,(1-\chi_{\rm
f})^{c_{5}+c_{6}\kappa\zeta},$ (33)

and similarly for the imaginary part, ${\rm Im}(\sigma)_{\rm dCS}=-1/\tau_{\rm
dCS}$. Here $\kappa=1/(16\pi)$, $\zeta=\ell_{\rm dCS}^{4}/(M_{\rm
s}^{4}\kappa)$, thus $\kappa\zeta=\gamma_{\rm dCS}$ and where $c_{i}$ (with
$i=1,\dots,6$) are fitting coefficients which can be found in Table II of Ref.
Wagle  _et al._ (2022),

We now expand Eq. (33) to leading orders in $\chi$ and $\gamma_{\rm dCS}$, and
gather the terms proportional to $\gamma_{\rm dCS}$. We obtain

$\displaystyle M\omega_{\rm dCS}=\left(0.3722+1.1945\gamma_{\rm
dCS}\right)+\left(0.1861+5.1828\gamma_{\rm dCS}\right)\,\chi,$

where we make use of the numerical values of the coefficients $c_{i}$. We find
(reassuringly) that the nonrotating GR part of the expression above agrees
with $\omega^{(0)}$ of Ref. Maselli  _et al._ (2020) to 0.5% relative error.
The same estimate leads to a larger relative error ($\approx 20$%) for the
linear-in-spin coefficient (i.e., 0.1861 in comparison to 0.1258 of Ref.
Maselli  _et al._ (2020)). We attribute this difference to Ref. Wagle  _et
al._ (2022) having fitted Eq. (33) to QNM data computed to linear-order in
spin, whereas Maselli  _et al._ (2020) fitted Eq. (6) to Kerr QNM valid to all
orders in spin.

We can now isolate the dCS corrections from Eq. (LABEL:eq:omega_dcs_interm)
and compare against Eq. (14), to find $p_{\rm dCS}=4$,

$\delta\omega^{(0)}_{\rm dCS}=3.1964\,,\quad\delta\omega^{(1)}_{\rm
dCS}=41.199.$ (35)

We can carry the same steps for $\tau_{\rm dCS}=-1/{\rm Im}(\sigma)_{\rm dCS}$
and find

$\delta\tau^{(0)}_{\rm dCS}=6.3619\,,\quad\delta\tau^{(1)}_{\rm dCS}=794.66,$
(36)

which completes the set of fixed non-GR parameters in the ringdown of the
pSEOBNR waveform model for this theory. We remark that the alarmingly large
values of $\delta\omega_{\rm dCS}^{(1)}$ and $\delta\tau_{\rm dCS}^{(1)}$ are
compensated by the assumptions that $\gamma_{\rm dCS}$ and $\chi$ are much
less than unity, which are indeed the assumptions used in Ref. Wagle  _et al._
(2022) to compute the QNMs.

###  A.3 Effective-field-theory of general relativity

The QNMs of slowly rotating BHs in both cubic and quartic EFT of GR where
calculated in Ref. Cano  _et al._ (2022). For the cubic EFT, we use their Eq.
(67), in the particular case of $\lambda_{\rm e}=\lambda_{\rm o}=1$. We then
linearize the resulting expression in $\chi$ and consider $m=2$ the harmonic.
As an outcome, we find that the fundamental axial-parity QNM is the least
damped one, and it is the one we use. Direct comparison with Eqs. (14) results
in $p_{\rm cEFT}=4$,

$\displaystyle\delta\omega_{\rm cEFT}^{(0)}=-0.5813,\quad\delta\tau_{\rm
cEFT}^{(0)}=2.6469,$ $\displaystyle\delta\omega_{\rm
cEFT}^{(1)}=-3.8620,\quad\delta\tau_{\rm cEFT}^{(1)}=265.12,$ (37)

for this theory.

We proceed in the same away for the quartic EFT. Here we use Eq. (68) (with
$\epsilon_{1}=1$ and $\epsilon_{2}=0$) and Eq. (70) of Ref. Cano  _et al._
(2022), In this case we find that both axial and polar modes reduce the
damping time of the fundamental QNM mode relative to GR. Hence, we choose the
axial-parity mode for this which reduction is the smallest. This time we then
find that $p_{\rm qEFT}=6$,

$\displaystyle\delta\omega_{\rm qEFT}^{(0)}=-0.2114,\quad\delta\tau_{\rm
qEFT}^{(0)}=-0.6070,$ $\displaystyle\delta\omega_{\rm
qEFT}^{(1)}=-1.5263,\quad\delta\tau_{\rm qEFT}^{(1)}=171.35,$ (38)

for this theory. As in the case of the previous theories, the large values of
some of these coefficients are compensated by the assumptions of weak coupling
and small spin used to calculate the QNM frequencies.

##  Appendix B The estimation of intrinsic binary parameters in General
Relativity and modified theories of gravity

We show in Fig. 12 a corner plot for all the intrinsic binary parameters from
our parameter-estimation study of GW150914, using the pSEOBNR waveform models
for GR and dCS. We also included $\ell_{\rm dCS}$ in the latter case. We see
that the medians of the posterior distributions are not affected substantially
by the inclusion of the non-GR parameters. This is particularly important for
the $M_{\rm f}$ and $\chi_{\rm f}$ parameters, since the ParSpec framework
assumes that the non-GR theory induces small deviations from GR. Indeed, since
pSEOBNR introduces only minimal modifications to the plunge-merger and because
the remnant BH parameters are estimated according to GR predictions using the
binary’s component masses and spins, the fact that only small biases are
introduced on $M_{\rm f}$ and $\chi_{\rm f}$ is, to some extend, expected. We
obtain qualitative similar results for GW200129 and the other modified gravity
theories considered in our work. We also remark that the fact that the
posterior on $\chi_{\rm f}$ has most support around $\approx 0.7$ justifies
our use of the fitting coefficients in the ParSpec formulas of Ref. Maselli
_et al._ (2020).

Figure 12:  Corner plot showing that the inferred source binary parameters and
$\ell_{\rm dCS}$ for GW150914. We used the same waveform model pSEOBNR setting
(purple contours) or not (blue contours) the non-GR parameters different from
zero. In the latter case, we considered dCS gravity and $N_{\rm max}=0$ as an
example. Here, $\chi_{\rm eff}$ is the dimensionless effective-spin parameter,
related to the individual spins $\chi_{i}$ and masses $m_{i}$ of each binary
component as $\chi_{\rm
eff}\equiv(m_{1}\chi_{1}+m_{2}\chi_{2})/(m_{1}+m_{2})$. All contours
correspond to 90% credible levels. We see that the addition of the non-GR
parameters does not introduce biases in the inference of the source
parameters. We found the same qualitative behavior in the posteriors
distributions of the source binary parameters for the other modified gravity
theories studied in the main text. The same conclusions apply for the other GW
event studied here, GW200129.

## References

  * Abbott  _et al._ (2016a) B. P. Abbott _et al._ (LIGO Scientific, Virgo), “Observation of Gravitational Waves from a Binary Black Hole Merger,” Phys. Rev. Lett. 116, 061102 (2016a), arXiv:1602.03837 [gr-qc] .
  * 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] .
  * Abbott  _et al._ (2021a) R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), “GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run,” (2021a), arXiv:2111.03606 [gr-qc] .
  * Abbott  _et al._ (2017a) B. P. Abbott _et al._ (LIGO Scientific, Virgo), “GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral,” Phys. Rev. Lett. 119, 161101 (2017a), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ (2018) B. P. Abbott _et al._ (LIGO Scientific, Virgo), “GW170817: Measurements of neutron star radii and equation of state,” Phys. Rev. Lett. 121, 161101 (2018), arXiv:1805.11581 [gr-qc] .
  * Abbott  _et al._ (2020) B. P. Abbott _et al._ (LIGO Scientific, Virgo), “GW190425: Observation of a Compact Binary Coalescence with Total Mass $\sim 3.4M_{\odot}$,” Astrophys. J. Lett. 892, L3 (2020), arXiv:2001.01761 [astro-ph.HE] .
  * Abbott  _et al._ (2021b) R. Abbott _et al._ (LIGO Scientific, KAGRA, VIRGO), “Observation of Gravitational Waves from Two Neutron Star–Black Hole Coalescences,” Astrophys. J. Lett. 915, L5 (2021b), arXiv:2106.15163 [astro-ph.HE] .
  * Nitz  _et al._ (2019a) A. H. Nitz, C. Capano, A. B. Nielsen, S. Reyes, R. White, D. A. Brown, and B. Krishnan, “1-OGC: The first open gravitational-wave catalog of binary mergers from analysis of public Advanced LIGO data,” Astrophys. J. 872, 195 (2019a), arXiv:1811.01921 [gr-qc] .
  * Nitz  _et al._ (2019b) A. H. Nitz, T. Dent, G. S. Davies, S. Kumar, C. D. Capano, I. Harry, S. Mozzon, L. Nuttall, A. Lundgren, and M. Tápai, “2-OGC: Open Gravitational-wave Catalog of binary mergers from analysis of public Advanced LIGO and Virgo data,” Astrophys. J. 891, 123 (2019b), arXiv:1910.05331 [astro-ph.HE] .
  * Venumadhav  _et al._ (2020) T. Venumadhav, B. Zackay, J. Roulet, L. Dai, and M. Zaldarriaga, “New binary black hole mergers in the second observing run of Advanced LIGO and Advanced Virgo,” Phys. Rev. D 101, 083030 (2020), arXiv:1904.07214 [astro-ph.HE] .
  * Zackay  _et al._ (2021) B. Zackay, L. Dai, T. Venumadhav, J. Roulet, and M. Zaldarriaga, “Detecting gravitational waves with disparate detector responses: Two new binary black hole mergers,” Phys. Rev. D 104, 063030 (2021), arXiv:1910.09528 [astro-ph.HE] .
  * Nitz  _et al._ (2021) A. H. Nitz, S. Kumar, Y.-F. Wang, S. Kastha, S. Wu, M. Schäfer, R. Dhurkunde, and C. D. Capano, “4-OGC: Catalog of gravitational waves from compact-binary mergers,” (2021), arXiv:2112.06878 [astro-ph.HE] .
  * Olsen  _et al._ (2022) S. Olsen, T. Venumadhav, J. Mushkin, J. Roulet, B. Zackay, and M. Zaldarriaga, “New binary black hole mergers in the LIGO-Virgo O3a data,” Phys. Rev. D 106, 043009 (2022), arXiv:2201.02252 [astro-ph.HE] .
  * Abbott  _et al._ (2021c) R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), “Constraints on the cosmic expansion history from GWTC-3,” (2021c), arXiv:2111.03604 [astro-ph.CO] .
  * Abbott  _et al._ (2021d) R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), “The population of merging compact binaries inferred using gravitational waves through GWTC-3,” (2021d), arXiv:2111.03634 [astro-ph.HE] .
  * 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. D 94, 084002 (2016), arXiv:1603.08955 [gr-qc] .
  * Abbott  _et al._ (2021e) R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), “Tests of General Relativity with GWTC-3,” (2021e), arXiv:2112.06861 [gr-qc] .
  * Will (2014) C. M. Will, “The Confrontation between General Relativity and Experiment,” Living Rev. Rel. 17, 4 (2014), arXiv:1403.7377 [gr-qc] .
  * Wex (2014) N. Wex, “Testing Relativistic Gravity with Radio Pulsars,” (2014), arXiv:1402.5594 [gr-qc] .
  * Kramer  _et al._ (2021) M. Kramer _et al._ , “Strong-Field Gravity Tests with the Double Pulsar,” Phys. Rev. X 11, 041050 (2021), arXiv:2112.06795 [astro-ph.HE] .
  * Abuter  _et al._ (2018) R. Abuter _et al._ (GRAVITY), “Detection of the gravitational redshift in the orbit of the star S2 near the Galactic centre massive black hole,” Astron. Astrophys. 615, L15 (2018), arXiv:1807.09409 [astro-ph.GA] .
  * Do  _et al._ (2019) T. Do _et al._ , “Relativistic redshift of the star S0-2 orbiting the Galactic center supermassive black hole,” Science 365, 664–668 (2019), arXiv:1907.10731 [astro-ph.GA] .
  * Akiyama  _et al._ (2019) K. Akiyama _et al._ (Event Horizon Telescope), “First M87 Event Horizon Telescope Results. IV. Imaging the Central Supermassive Black Hole,” Astrophys. J. Lett. 875, L4 (2019), arXiv:1906.11241 [astro-ph.GA] .
  * Clifton  _et al._ (2012) T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis, “Modified Gravity and Cosmology,” Phys. Rept. 513, 1–189 (2012), arXiv:1106.2476 [astro-ph.CO] .
  * Vishveshwara (1970) C. V. Vishveshwara, “Stability of the Schwarzschild metric,” Phys. Rev. D 1, 2870–2879 (1970).
  * Press (1971) W. H. Press, “Long Wave Trains of Gravitational Waves from a Vibrating Black Hole,” Astrophys. J. Lett. 170, L105–L108 (1971).
  * Chandrasekhar and Detweiler (1975) S. Chandrasekhar and S. L. Detweiler, “The quasi-normal modes of the Schwarzschild black hole,” Proc. Roy. Soc. Lond. A 344, 441–452 (1975).
  * Carter (1971) B. Carter, “Axisymmetric Black Hole Has Only Two Degrees of Freedom,” Phys. Rev. Lett. 26, 331–333 (1971).
  * Israel (1967) W. Israel, “Event horizons in static vacuum space-times,” Phys. Rev. 164, 1776–1779 (1967).
  * Hawking (1972) S. W. Hawking, “Black holes in general relativity,” Commun. Math. Phys. 25, 152–166 (1972).
  * Robinson (1975) D. C. Robinson, “Uniqueness of the Kerr black hole,” Phys. Rev. Lett. 34, 905–906 (1975).
  * Abbott  _et al._ (2016b) B. P. Abbott _et al._ (LIGO Scientific, Virgo), “Tests of general relativity with GW150914,” Phys. Rev. Lett. 116, 221101 (2016b), [Erratum: Phys.Rev.Lett. 121, 129902 (2018)], arXiv:1602.03841 [gr-qc] .
  * Ghosh  _et al._ (2016) A. Ghosh _et al._ , “Testing general relativity using golden black-hole binaries,” Phys. Rev. D 94, 021101(R) (2016), arXiv:1602.02453 [gr-qc] .
  * Ghosh  _et al._ (2018) A. Ghosh, N. K. Johnson-Mcdaniel, A. Ghosh, C. K. Mishra, P. Ajith, W. Del Pozzo, C. P. Berry, A. B. Nielsen, and L. London, “Testing general relativity using gravitational wave signals from the inspiral, merger and ringdown of binary black holes,” Class. Quant. Grav. 35, 014002 (2018), arXiv:1704.06784 [gr-qc] .
  * Abbott  _et al._ (2019a) B. P. Abbott _et al._ (LIGO Scientific, Virgo), “Tests of General Relativity with GW170817,” Phys. Rev. Lett. 123, 011102 (2019a), arXiv:1811.00364 [gr-qc] .
  * Abbott  _et al._ (2019b) 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 (2019b), arXiv:1903.04467 [gr-qc] .
  * Abbott  _et al._ (2021f) 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 (2021f), arXiv:2010.14529 [gr-qc] .
  * Arun  _et al._ (2006) 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–L43 (2006), arXiv:gr-qc/0604018 .
  * Yunes and Pretorius (2009a) N. Yunes and F. Pretorius, “Fundamental Theoretical Bias in Gravitational Wave Astrophysics and the Parameterized Post-Einsteinian Framework,” Phys. Rev. D 80, 122003 (2009a), arXiv:0909.3328 [gr-qc] .
  * 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] .
  * Mehta  _et al._ (2022) A. K. Mehta, A. Buonanno, R. Cotesta, A. Ghosh, N. Sennett, and J. Steinhoff, “Tests of General Relativity with Gravitational-Wave Observations using a Flexible–Theory-Independent Method,” (2022), arXiv:2203.13937 [gr-qc] .
  * Will (1998) C. M. Will, “Bounding the mass of the graviton using gravitational wave observations of inspiralling compact binaries,” Phys. Rev. D 57, 2061–2068 (1998), arXiv:gr-qc/9709011 .
  * Mirshekari  _et al._ (2012) S. Mirshekari, N. Yunes, and C. M. Will, “Constraining Generic Lorentz Violation and the Speed of the Graviton with Gravitational Waves,” Phys. Rev. D 85, 024041 (2012), arXiv:1110.2720 [gr-qc] .
  * Eardley  _et al._ (1973) D. M. Eardley, D. L. Lee, and A. P. Lightman, “Gravitational-wave observations as a tool for testing relativistic gravity,” Phys. Rev. D 8, 3308–3321 (1973).
  * Isi and Weinstein (2017) M. Isi and A. J. Weinstein, “Probing gravitational wave polarizations with signals from compact binary coalescences,” (2017), arXiv:1710.03794 [gr-qc] .
  * Abbott  _et al._ (2017b) B. P. Abbott _et al._ (LIGO Scientific, Virgo), “GW170814: A Three-Detector Observation of Gravitational Waves from a Binary Black Hole Coalescence,” Phys. Rev. Lett. 119, 141101 (2017b), arXiv:1709.09660 [gr-qc] .
  * Abbott  _et al._ (2019c) B. P. Abbott _et al._ (LIGO Scientific, Virgo), “Tests of General Relativity with GW170817,” Phys. Rev. Lett. 123, 011102 (2019c), arXiv:1811.00364 [gr-qc] .
  * Wong  _et al._ (2021) I. C. F. Wong, P. T. H. Pang, R. K. L. Lo, T. G. F. Li, and C. Van Den Broeck, “Null-stream-based Bayesian Unmodeled Framework to Probe Generic Gravitational-wave Polarizations,” (2021), arXiv:2105.09485 [gr-qc] .
  * Meidam  _et al._ (2014) J. Meidam, M. Agathos, C. Van Den Broeck, J. Veitch, and B. S. Sathyaprakash, “Testing the no-hair theorem with black hole ringdowns using TIGER,” Phys. Rev. D 90, 064009 (2014), arXiv:1406.3201 [gr-qc] .
  * Carullo  _et al._ (2018) G. Carullo _et al._ , “Empirical tests of the black hole no-hair conjecture using gravitational-wave observations,” Phys. Rev. D 98, 104020 (2018), arXiv:1805.04760 [gr-qc] .
  * Brito  _et al._ (2018) R. Brito, A. Buonanno, and V. Raymond, “Black-hole Spectroscopy by Making Full Use of Gravitational-Wave Modeling,” Phys. Rev. D 98, 084038 (2018), arXiv:1805.00293 [gr-qc] .
  * Carullo  _et al._ (2019) G. Carullo, W. Del Pozzo, and J. Veitch, “Observational Black Hole Spectroscopy: A time-domain multimode analysis of GW150914,” Phys. Rev. D 99, 123029 (2019), [Erratum: Phys.Rev.D 100, 089903 (2019)], arXiv:1902.07527 [gr-qc] .
  * Isi  _et al._ (2019) M. Isi, M. Giesler, W. M. Farr, M. A. Scheel, and S. A. Teukolsky, “Testing the no-hair theorem with GW150914,” Phys. Rev. Lett. 123, 111102 (2019), arXiv:1905.00869 [gr-qc] .
  * Ghosh  _et al._ (2021) A. Ghosh, R. Brito, and A. Buonanno, “Constraints on quasinormal-mode frequencies with LIGO-Virgo binary–black-hole observations,” Phys. Rev. D 103, 124041 (2021), arXiv:2104.01906 [gr-qc] .
  * Carullo (2021) G. Carullo, “Enhancing modified gravity detection from gravitational-wave observations using the parametrized ringdown spin expansion coeffcients formalism,” Phys. Rev. D 103, 124043 (2021), arXiv:2102.05939 [gr-qc] .
  * Dreyer  _et al._ (2004) O. Dreyer, B. J. Kelly, B. Krishnan, L. S. Finn, D. Garrison, and R. Lopez-Aleman, “Black hole spectroscopy: Testing general relativity through gravitational wave observations,” Class. Quant. Grav. 21, 787–804 (2004), arXiv:gr-qc/0309007 .
  * Detweiler (1980) S. L. Detweiler, “Black holes and gravitational waves. III. The resonant frequencies of rotating holes,” Astrophys. J. 239, 292–295 (1980).
  * Berti  _et al._ (2006) E. Berti, V. Cardoso, and C. M. Will, “On gravitational-wave spectroscopy of massive black holes with the space interferometer LISA,” Phys. Rev. D 73, 064030 (2006), arXiv:gr-qc/0512160 .
  * Baibhav  _et al._ (2020) V. Baibhav, E. Berti, and V. Cardoso, “LISA parameter estimation and source localization with higher harmonics of the ringdown,” Phys. Rev. D 101, 084053 (2020), arXiv:2001.10011 [gr-qc] .
  * Bhagwat  _et al._ (2022) S. Bhagwat, C. Pacilio, E. Barausse, and P. Pani, “Landscape of massive black-hole spectroscopy with LISA and the Einstein Telescope,” Phys. Rev. D 105, 124063 (2022), arXiv:2201.00023 [gr-qc] .
  * Ota and Chirenti (2022) I. Ota and C. Chirenti, “Black hole spectroscopy horizons for current and future gravitational wave detectors,” Phys. Rev. D 105, 044015 (2022), arXiv:2108.01774 [gr-qc] .
  * Capano  _et al._ (2021) C. D. Capano, M. Cabero, J. Westerweck, J. Abedi, S. Kastha, A. H. Nitz, A. B. Nielsen, and B. Krishnan, “Observation of a multimode quasi-normal spectrum from a perturbed black hole,” (2021), arXiv:2105.05238 [gr-qc] .
  * Isi and Farr (2022) M. Isi and W. M. Farr, “Revisiting the ringdown of GW150914,” (2022), arXiv:2202.02941 [gr-qc] .
  * Finch and Moore (2022) E. Finch and C. J. Moore, “Searching for a ringdown overtone in GW150914,” Phys. Rev. D 106, 043005 (2022), arXiv:2205.07809 [gr-qc] .
  * Cotesta  _et al._ (2022) R. Cotesta, G. Carullo, E. Berti, and V. Cardoso, “Analysis of Ringdown Overtones in GW150914,” Phys. Rev. Lett. 129, 111102 (2022), arXiv:2201.00822 [gr-qc] .
  * Isi and Farr (2021) M. Isi and W. M. Farr, “Analyzing black-hole ringdowns,” (2021), arXiv:2107.05609 [gr-qc] .
  * Bohé  _et al._ (2017) A. Bohé _et al._ , “Improved effective-one-body model of spinning, nonprecessing binary black holes for the era of gravitational-wave astrophysics with advanced detectors,” Phys. Rev. D 95, 044028 (2017), arXiv:1611.03703 [gr-qc] .
  * Cotesta  _et al._ (2018) R. Cotesta, A. Buonanno, A. Bohé, A. Taracchini, I. Hinder, and S. Ossokine, “Enriching the Symphony of Gravitational Waves from Binary Black Holes by Tuning Higher Harmonics,” Phys. Rev. D 98, 084028 (2018), arXiv:1803.10701 [gr-qc] .
  * Mihaylov  _et al._ (2021) D. P. Mihaylov, S. Ossokine, A. Buonanno, and A. Ghosh, “Fast post-adiabatic waveforms in the time domain: Applications to compact binary coalescences in LIGO and Virgo,” Phys. Rev. D 104, 124087 (2021), arXiv:2105.06983 [gr-qc] .
  * Maselli  _et al._ (2020) A. Maselli, P. Pani, L. Gualtieri, and E. Berti, “Parametrized ringdown spin expansion coefficients: a data-analysis framework for black-hole spectroscopy with multiple events,” Phys. Rev. D 101, 024043 (2020), arXiv:1910.12893 [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] .
  * Silva  _et al._ (2021a) H. O. Silva, A. M. Holgado, A. Cárdenas-Avendaño, and N. Yunes, “Astrophysical and theoretical physics implications from multimessenger neutron star observations,” Phys. Rev. Lett. 126, 181101 (2021a), arXiv:2004.01253 [gr-qc] .
  * Sennett  _et al._ (2020) N. Sennett, R. Brito, A. Buonanno, V. Gorbenko, and L. Senatore, “Gravitational-Wave Constraints on an Effective Field-Theory Extension of General Relativity,” Phys. Rev. D 102, 044056 (2020), arXiv:1912.09917 [gr-qc] .
  * Myers (1987) R. C. Myers, “Higher Derivative Gravity, Surface Terms and String Theory,” Phys. Rev. D 36, 392 (1987).
  * 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–5058 (1996), arXiv:hep-th/9511071 .
  * 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] .
  * 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] .
  * Macedo  _et al._ (2019) C. F. B. Macedo, J. Sakstein, E. Berti, L. Gualtieri, H. O. Silva, and T. P. Sotiriou, “Self-interactions and Spontaneous Black Hole Scalarization,” Phys. Rev. D 99, 104041 (2019), arXiv:1903.06784 [gr-qc] .
  * Dima  _et al._ (2020) A. Dima, E. Barausse, N. Franchini, and T. P. Sotiriou, “Spin-induced black hole spontaneous scalarization,” Phys. Rev. Lett. 125, 231101 (2020), arXiv:2006.03095 [gr-qc] .
  * Herdeiro  _et al._ (2021) C. A. R. Herdeiro, E. Radu, H. O. Silva, T. P. Sotiriou, and N. Yunes, “Spin-induced scalarized black holes,” Phys. Rev. Lett. 126, 011103 (2021), arXiv:2009.03904 [gr-qc] .
  * Berti  _et al._ (2021) E. Berti, L. G. Collodel, B. Kleihaus, and J. Kunz, “Spin-induced black-hole scalarization in Einstein-scalar-Gauss-Bonnet theory,” Phys. Rev. Lett. 126, 011104 (2021), arXiv:2009.03905 [gr-qc] .
  * Coleman  _et al._ (1992) S. R. Coleman, J. Preskill, and F. Wilczek, “Quantum hair on black holes,” Nucl. Phys. B 378, 175–246 (1992), arXiv:hep-th/9201059 .
  * Herdeiro and Radu (2015) C. A. R. Herdeiro and E. Radu, “Asymptotically flat black holes with scalar hair: a review,” Int. J. Mod. Phys. D 24, 1542014 (2015), arXiv:1504.08209 [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] .
  * Julié and Berti (2019) F.-L. Julié and E. Berti, “Post-Newtonian dynamics and black hole thermodynamics in Einstein-scalar-Gauss-Bonnet gravity,” Phys. Rev. D 100, 104061 (2019), arXiv:1909.05258 [gr-qc] .
  * Shiralilou  _et al._ (2021) B. Shiralilou, T. Hinderer, S. Nissanke, N. Ortiz, and H. Witek, “Nonlinear curvature effects in gravitational waves from inspiralling black hole binaries,” Phys. Rev. D 103, L121503 (2021), arXiv:2012.09162 [gr-qc] .
  * Shiralilou  _et al._ (2022) B. Shiralilou, T. Hinderer, S. M. Nissanke, N. Ortiz, and H. Witek, “Post-Newtonian gravitational and scalar waves in scalar-Gauss–Bonnet gravity,” Class. Quant. Grav. 39, 035002 (2022), arXiv:2105.13972 [gr-qc] .
  * Julié  _et al._ (2022) F.-L. Julié, H. O. Silva, E. Berti, and N. Yunes, “Black hole sensitivities in Einstein-scalar-Gauss-Bonnet gravity,” Phys. Rev. D 105, 124031 (2022), arXiv:2202.01329 [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] .
  * Perkins  _et al._ (2021) S. E. Perkins, R. Nair, H. O. Silva, and N. Yunes, “Improved gravitational-wave constraints on higher-order curvature theories of gravity,” Phys. Rev. D 104, 024060 (2021), arXiv:2104.11189 [gr-qc] .
  * Abbott  _et al._ (2019d) B. P. Abbott _et al._ (LIGO Scientific, Virgo), “GWTC-1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs,” Phys. Rev. X 9, 031040 (2019d), arXiv:1811.12907 [astro-ph.HE] .
  * Abbott  _et al._ (2021g) R. Abbott _et al._ (LIGO Scientific, Virgo), “GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run,” Phys. Rev. X 11, 021053 (2021g), arXiv:2010.14527 [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] .
  * 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] .
  * Sotiriou and Zhou (2014a) T. P. Sotiriou and S.-Y. Zhou, “Black hole hair in generalized scalar-tensor gravity,” Phys. Rev. Lett. 112, 251102 (2014a), arXiv:1312.3622 [gr-qc] .
  * Sotiriou and Zhou (2014b) T. P. Sotiriou and S.-Y. Zhou, “Black hole hair in generalized scalar-tensor gravity: An explicit example,” Phys. Rev. D 90, 124063 (2014b), arXiv:1408.1698 [gr-qc] .
  * Glampedakis  _et al._ (2017a) K. Glampedakis, A. D. Johnson, and D. Kennefick, “Darboux transformation in black hole perturbation theory,” Phys. Rev. D 96, 024036 (2017a), arXiv:1702.06459 [gr-qc] .
  * Lenzi and Sopuerta (2021) M. Lenzi and C. F. Sopuerta, “Darboux covariance: A hidden symmetry of perturbed Schwarzschild black holes,” Phys. Rev. D 104, 124068 (2021), arXiv:2109.00503 [gr-qc] .
  * Chandrasekhar (1985) S. Chandrasekhar,  _The mathematical theory of black holes_ (1985).
  * Pani and Cardoso (2009) P. Pani and V. Cardoso, “Are black holes in alternative theories serious astrophysical candidates? The Case for Einstein-Dilaton-Gauss-Bonnet black holes,” Phys. Rev. D 79, 084031 (2009), arXiv:0902.1569 [gr-qc] .
  * Blázquez-Salcedo  _et al._ (2016) J. L. Blázquez-Salcedo, C. F. B. Macedo, V. Cardoso, V. Ferrari, L. Gualtieri, F. S. Khoo, J. Kunz, and P. Pani, “Perturbed black holes in Einstein-dilaton-Gauss-Bonnet gravity: Stability, ringdown, and gravitational-wave emission,” Phys. Rev. D 94, 104024 (2016), arXiv:1609.01286 [gr-qc] .
  * Pierini and Gualtieri (2021) L. Pierini and L. Gualtieri, “Quasi-normal modes of rotating black holes in Einstein-dilaton Gauss-Bonnet gravity: the first order in rotation,” Phys. Rev. D 103, 124017 (2021), arXiv:2103.09870 [gr-qc] .
  * Bryant  _et al._ (2021) A. Bryant, H. O. Silva, K. Yagi, and K. Glampedakis, “Eikonal quasinormal modes of black holes beyond general relativity. III. Scalar Gauss-Bonnet gravity,” Phys. Rev. D 104, 044051 (2021), arXiv:2106.09657 [gr-qc] .
  * Witek  _et al._ (2019) H. Witek, L. Gualtieri, P. Pani, and T. P. Sotiriou, “Black holes and binary mergers in scalar Gauss-Bonnet gravity: scalar field dynamics,” Phys. Rev. D 99, 064035 (2019), arXiv:1810.05177 [gr-qc] .
  * Okounkova (2019) M. Okounkova, “Stability of Rotating Black Holes in Einstein Dilaton Gauss-Bonnet Gravity,” Phys. Rev. D 100, 124054 (2019), arXiv:1909.12251 [gr-qc] .
  * Jackiw and Pi (2003) R. Jackiw and S. Y. Pi, “Chern-Simons modification of general relativity,” Phys. Rev. D 68, 104012 (2003), arXiv:gr-qc/0308071 .
  * Alexander and Yunes (2009) S. Alexander and N. Yunes, “Chern-Simons Modified General Relativity,” Phys. Rept. 480, 1–55 (2009), arXiv:0907.2562 [hep-th] .
  * Yunes and Pretorius (2009b) N. Yunes and F. Pretorius, “Dynamical Chern-Simons Modified Gravity. I. Spinning Black Holes in the Slow-Rotation Approximation,” Phys. Rev. D 79, 084043 (2009b), arXiv:0902.4669 [gr-qc] .
  * Konno  _et al._ (2009) K. Konno, T. Matsuyama, and S. Tanda, “Rotating black hole in extended Chern-Simons modified gravity,” Prog. Theor. Phys. 122, 561–568 (2009), arXiv:0902.4767 [gr-qc] .
  * Yamada and Tanaka (2020) K. Yamada and T. Tanaka, “Parametrized test of parity-violating gravity using GWTC-1 events,” PTEP 2020, 093E01 (2020), arXiv:2006.11086 [gr-qc] .
  * Wang  _et al._ (2021) Y.-F. Wang, R. Niu, T. Zhu, and W. Zhao, “Gravitational Wave Implications for the Parity Symmetry of Gravity in the High Energy Region,” Astrophys. J. 908, 58 (2021), arXiv:2002.05668 [gr-qc] .
  * Okounkova  _et al._ (2022) M. Okounkova, W. M. Farr, M. Isi, and L. C. Stein, “Constraining gravitational wave amplitude birefringence and Chern-Simons gravity with GWTC-2,” Phys. Rev. D 106, 044067 (2022), arXiv:2101.11153 [gr-qc] .
  * Alexander and Martin (2005) S. Alexander and J. Martin, “Birefringent gravitational waves and the consistency check of inflation,” Phys. Rev. D 71, 063526 (2005), arXiv:hep-th/0410230 .
  * Lommen  _et al._ (2000) A. N. Lommen, A. Zepka, D. C. Backer, M. McLaughlin, J. C. Cordes, Z. Arzoumanian, and K. Xilouris, “New pulsars from an Arecibo drift scan search,” Astrophys. J. 545, 1007 (2000), arXiv:astro-ph/0008054 .
  * Arzoumanian  _et al._ (2018) Z. Arzoumanian _et al._ (NANOGrav), “The NANOGrav 11-year Data Set: High-precision timing of 45 Millisecond Pulsars,” Astrophys. J. Suppl. 235, 37 (2018), arXiv:1801.01837 [astro-ph.HE] .
  * Riley  _et al._ (2019) T. E. Riley _et al._ , “A NICER View of PSR J0030+0451: Millisecond Pulsar Parameter Estimation,” Astrophys. J. Lett. 887, L21 (2019), arXiv:1912.05702 [astro-ph.HE] .
  * Miller  _et al._ (2019) M. C. Miller _et al._ , “PSR J0030+0451 Mass and Radius from NICER Data and Implications for the Properties of Neutron Star Matter,” Astrophys. J. Lett. 887, L24 (2019), arXiv:1912.05705 [astro-ph.HE] .
  * Yagi and Yunes (2013a) K. Yagi and N. Yunes, “I-Love-Q,” Science 341, 365–368 (2013a), arXiv:1302.4499 [gr-qc] .
  * Yagi and Yunes (2013b) K. Yagi and N. Yunes, “I-Love-Q Relations in Neutron Stars and their Applications to Astrophysics, Gravitational Waves and Fundamental Physics,” Phys. Rev. D 88, 023009 (2013b), arXiv:1303.1528 [gr-qc] .
  * Gupta  _et al._ (2018) T. Gupta, B. Majumder, K. Yagi, and N. Yunes, “I-Love-Q Relations for Neutron Stars in dynamical Chern Simons Gravity,” Class. Quant. Grav. 35, 025009 (2018), arXiv:1710.07862 [gr-qc] .
  * Yunes and Sopuerta (2008) N. Yunes and C. F. Sopuerta, “Perturbations of Schwarzschild Black Holes in Chern-Simons Modified Gravity,” Phys. Rev. D 77, 064007 (2008), arXiv:0712.1028 [gr-qc] .
  * Cardoso and Gualtieri (2009) V. Cardoso and L. Gualtieri, “Perturbations of Schwarzschild black holes in Dynamical Chern-Simons modified gravity,” Phys. Rev. D 80, 064008 (2009), [Erratum: Phys.Rev.D 81, 089903 (2010)], arXiv:0907.5008 [gr-qc] .
  * Molina  _et al._ (2010) C. Molina, P. Pani, V. Cardoso, and L. Gualtieri, “Gravitational signature of Schwarzschild black holes in dynamical Chern-Simons gravity,” Phys. Rev. D 81, 124021 (2010), arXiv:1004.4007 [gr-qc] .
  * Srivastava  _et al._ (2021) M. Srivastava, Y. Chen, and S. Shankaranarayanan, “Analytical computation of quasinormal modes of slowly rotating black holes in dynamical Chern-Simons gravity,” Phys. Rev. D 104, 064034 (2021), arXiv:2106.06209 [gr-qc] .
  * Wagle  _et al._ (2022) P. Wagle, N. Yunes, and H. O. Silva, “Quasinormal modes of slowly-rotating black holes in dynamical Chern-Simons gravity,” Phys. Rev. D 105, 124003 (2022), arXiv:2103.09913 [gr-qc] .
  * Okounkova  _et al._ (2019) M. Okounkova, L. C. Stein, M. A. Scheel, and S. A. Teukolsky, “Numerical binary black hole collisions in dynamical Chern-Simons gravity,” Phys. Rev. D 100, 104026 (2019), arXiv:1906.08789 [gr-qc] .
  * Endlich  _et al._ (2017) S. Endlich, V. Gorbenko, J. Huang, and L. Senatore, “An effective formalism for testing extensions to General Relativity with gravitational waves,” JHEP 09, 122 (2017), arXiv:1704.01590 [gr-qc] .
  * Cano and Ruipérez (2019) P. A. Cano and A. Ruipérez, “Leading higher-derivative corrections to Kerr geometry,” JHEP 05, 189 (2019), [Erratum: JHEP 03, 187 (2020)], arXiv:1901.01315 [gr-qc] .
  * Brandhuber and Travaglini (2020) A. Brandhuber and G. Travaglini, “On higher-derivative effects on the gravitational potential and particle bending,” JHEP 01, 010 (2020), arXiv:1905.05657 [hep-th] .
  * Accettulli Huber  _et al._ (2021) M. Accettulli Huber, A. Brandhuber, S. De Angelis, and G. Travaglini, “From amplitudes to gravitational radiation with cubic interactions and tidal effects,” Phys. Rev. D 103, 045015 (2021), arXiv:2012.06548 [hep-th] .
  * de Rham  _et al._ (2020) C. de Rham, J. Francfort, and J. Zhang, “Black Hole Gravitational Waves in the Effective Field Theory of Gravity,” Phys. Rev. D 102, 024079 (2020), arXiv:2005.13923 [hep-th] .
  * Cano  _et al._ (2020) P. A. Cano, K. Fransen, and T. Hertog, “Ringing of rotating black holes in higher-derivative gravity,” Phys. Rev. D 102, 044047 (2020), arXiv:2005.03671 [gr-qc] .
  * Cano  _et al._ (2022) P. A. Cano, K. Fransen, T. Hertog, and S. Maenaut, “Gravitational ringing of rotating black holes in higher-derivative gravity,” Phys. Rev. D 105, 024064 (2022), arXiv:2110.11378 [gr-qc] .
  * Cardoso  _et al._ (2018) V. Cardoso, M. Kimura, A. Maselli, and L. Senatore, “Black Holes in an Effective Field Theory Extension of General Relativity,” Phys. Rev. Lett. 121, 251105 (2018), arXiv:1808.08962 [gr-qc] .
  * Glampedakis  _et al._ (2017b) K. Glampedakis, G. Pappas, H. O. Silva, and E. Berti, “Post-Kerr black hole spectroscopy,” Phys. Rev. D 96, 064054 (2017b), arXiv:1706.07658 [gr-qc] .
  * Glampedakis and Silva (2019) K. Glampedakis and H. O. Silva, “Eikonal quasinormal modes of black holes beyond General Relativity,” Phys. Rev. D 100, 044040 (2019), arXiv:1906.05455 [gr-qc] .
  * Silva and Glampedakis (2020) H. O. Silva and K. Glampedakis, “Eikonal quasinormal modes of black holes beyond general relativity. II. Generalized scalar-tensor perturbations,” Phys. Rev. D 101, 044051 (2020), arXiv:1912.09286 [gr-qc] .
  * Cardoso  _et al._ (2019) V. Cardoso, M. Kimura, A. Maselli, E. Berti, C. F. B. Macedo, and R. McManus, “Parametrized black hole quasinormal ringdown: Decoupled equations for nonrotating black holes,” Phys. Rev. D 99, 104077 (2019), arXiv:1901.01265 [gr-qc] .
  * McManus  _et al._ (2019) R. McManus, E. Berti, C. F. B. Macedo, M. Kimura, A. Maselli, and V. Cardoso, “Parametrized black hole quasinormal ringdown. II. Coupled equations and quadratic corrections for nonrotating black holes,” Phys. Rev. D 100, 044061 (2019), arXiv:1906.05155 [gr-qc] .
  * Kimura (2020) M. Kimura, “Note on the parametrized black hole quasinormal ringdown formalism,” Phys. Rev. D 101, 064031 (2020), arXiv:2001.09613 [gr-qc] .
  * Gossan  _et al._ (2012) S. Gossan, J. Veitch, and B. S. Sathyaprakash, “Bayesian model selection for testing the no-hair theorem with black hole ringdowns,” Phys. Rev. D 85, 124056 (2012), arXiv:1111.5819 [gr-qc] .
  * Abbott  _et al._ (2021h) 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 (2021h), arXiv:2010.14529 [gr-qc] .
  * Krolak and Schutz (1987) A. Krolak and B. F. Schutz, “Coalescing binaries – Probe of the universe,” Gen. Rel. Grav. 19, 1163–1171 (1987).
  * Nagar and Rettegno (2019) A. Nagar and P. Rettegno, “Efficient effective one body time-domain gravitational waveforms,” Phys. Rev. D 99, 021501 (2019), arXiv:1805.03891 [gr-qc] .
  * LIGO Scientific Collaboration (2018) LIGO Scientific Collaboration, “LIGO Algorithm Library - LALSuite,” free software (GPL) (2018).
  * Taracchini  _et al._ (2014) A. Taracchini _et al._ , “Effective-one-body model for black-hole binaries with generic mass ratios and spins,” Phys. Rev. D 89, 061502 (2014), arXiv:1311.2544 [gr-qc] .
  * Hofmann  _et al._ (2016) F. Hofmann, E. Barausse, and L. Rezzolla, “The final spin from binary black holes in quasi-circular orbits,” Astrophys. J. Lett. 825, L19 (2016), arXiv:1605.01938 [gr-qc] .
  * Berti  _et al._ (2009) E. Berti, V. Cardoso, and A. O. Starinets, “Quasinormal modes of black holes and black branes,” Class. Quant. Grav. 26, 163001 (2009), arXiv:0905.2975 [gr-qc] .
  * Regge and Wheeler (1957) T. Regge and J. A. Wheeler, “Stability of a Schwarzschild singularity,” Phys. Rev. 108, 1063–1069 (1957).
  * Zerilli (1970) F. J. Zerilli, “Effective potential for even parity Regge-Wheeler gravitational perturbation equations,” Phys. Rev. Lett. 24, 737–738 (1970).
  * Maggiore (2018) M. Maggiore,  _Gravitational Waves. Vol. 2: Astrophysics and Cosmology_ (Oxford University Press, 2018).
  * Hui  _et al._ (2021) L. Hui, A. Podo, L. Santoni, and E. Trincherini, “Effective Field Theory for the perturbations of a slowly rotating black hole,” JHEP 12, 183 (2021), arXiv:2111.02072 [hep-th] .
  * Rover  _et al._ (2006) C. Rover, R. Meyer, and N. Christensen, “Bayesian inference on compact binary inspiral gravitational radiation signals in interferometric data,” Class. Quant. Grav. 23, 4895–4906 (2006), arXiv:gr-qc/0602067 .
  * van der Sluys  _et al._ (2008) M. van der Sluys, V. Raymond, I. Mandel, C. Rover, N. Christensen, V. Kalogera, R. Meyer, and A. Vecchio, “Parameter estimation of spinning binary inspirals using Markov-chain Monte Carlo,” Class. Quant. Grav. 25, 184011 (2008), arXiv:0805.1689 [gr-qc] .
  * Veitch  _et al._ (2015) J. Veitch _et al._ , “Parameter estimation for compact binaries with ground-based gravitational-wave observations using the LALInference software library,” Phys. Rev. D 91, 042003 (2015), arXiv:1409.7215 [gr-qc] .
  * Ossokine  _et al._ (2020) S. Ossokine _et al._ , “Multipolar Effective-One-Body Waveforms for Precessing Binary Black Holes: Construction and Validation,” Phys. Rev. D 102, 044055 (2020), arXiv:2004.09442 [gr-qc] .
  * Ramos-Buades  _et al._ (2022) A. Ramos-Buades, A. Buonanno, M. Khalil, and S. Ossokine, “Effective-one-body multipolar waveforms for eccentric binary black holes with nonprecessing spins,” Phys. Rev. D 105, 044035 (2022), arXiv:2112.06952 [gr-qc] .
  * Healy  _et al._ (2012) J. Healy, T. Bode, R. Haas, E. Pazos, P. Laguna, D. Shoemaker, and N. Yunes, “Late Inspiral and Merger of Binary Black Holes in Scalar-Tensor Theories of Gravity,” Class. Quant. Grav. 29, 232002 (2012), arXiv:1112.3928 [gr-qc] .
  * Okounkova  _et al._ (2017) M. Okounkova, L. C. Stein, M. A. Scheel, and D. A. Hemberger, “Numerical binary black hole mergers in dynamical Chern-Simons gravity: Scalar field,” Phys. Rev. D 96, 044020 (2017), arXiv:1705.07924 [gr-qc] .
  * Hirschmann  _et al._ (2018) E. W. Hirschmann, L. Lehner, S. L. Liebling, and C. Palenzuela, “Black Hole Dynamics in Einstein-Maxwell-Dilaton Theory,” Phys. Rev. D 97, 064032 (2018), arXiv:1706.09875 [gr-qc] .
  * Okounkova (2020) M. Okounkova, “Numerical relativity simulation of GW150914 in Einstein dilaton Gauss-Bonnet gravity,” Phys. Rev. D 102, 084046 (2020), arXiv:2001.03571 [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] .
  * East and Ripley (2021a) W. E. East and J. L. Ripley, “Evolution of Einstein-scalar-Gauss-Bonnet gravity using a modified harmonic formulation,” Phys. Rev. D 103, 044040 (2021a), arXiv:2011.03547 [gr-qc] .
  * East and Ripley (2021b) W. E. East and J. L. Ripley, “Dynamics of Spontaneous Black Hole Scalarization and Mergers in Einstein-Scalar-Gauss-Bonnet Gravity,” Phys. Rev. Lett. 127, 101102 (2021b), arXiv:2105.08571 [gr-qc] .
  * Figueras and França (2022) P. Figueras and T. França, “Black hole binaries in cubic Horndeski theories,” Phys. Rev. D 105, 124004 (2022), arXiv:2112.15529 [gr-qc] .