key: cord-0668522-zu369ww0
authors: Silva, Hector O.; Ghosh, Abhirup; Buonanno, Alessandra
title: Black-hole ringdown as a probe of higher-curvature gravity theories
date: 2022-05-10
journal: nan
DOI: nan
sha: 2cfb39d6d0fd4fac4dcdc1396b19160fca83fa42
doc_id: 668522
cord_uid: zu369ww0

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 (EdGB), dynamical Chern-Simons (dCS) theories, and cubic- and quartic-order effective-field-theories (EFT) of general relativity. We use this parametrized waveform model to measure the ringdown properties of the two loudests ringdown signals observed so far, GW150914 and GW200129. We find that while EdGB theory cannot be constrained with these events, we can place upper bounds on the fundamental length-scale of cubic- ($ell_{rm cEFT} leqslant 38.2$ km) and quartic-order ($ell_{rm qEFT} leqslant 51.3$ km) EFTs of general relativity, and of dCS 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.

Since the first detection of gravitational waves (GWs) from a binary black-hole (BBH) merger in 2015 [1] , the LIGO [2] and Virgo [3] detectors have observed about 90 GW events [4] from mergers of BHs, neutron stars (NSs) [5] [6] [7] and their mixture [8] . These results have been confirmed by independent analyses, which have also identified a few new GW signals [9] [10] [11] [12] [13] [14] .

The large number of GW observations has allowed us to infer relevant astrophysical [15] and cosmological [16] 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 [17, 18] . The latter complement tests of GR in the low-velocity, quasi-static or linear regimes available with Solar-System experiments [19] , binary-pulsar [20, 21] and galactic-center [22] [23] [24] observations, and cosmological measurements [25] .

The coalescence of two BHs in GR is characterized by a long inspiral stage, during which the holes adiabatically and steadly come closer and closer to each other, loosing 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) [26] [27] [28] . Because of the no-hair conjecture in GR [29] [30] [31] [32] , 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 ( , m) and the overtone number n.

Several null tests have been proposed to probe the nature of gravity with GW signals [17, 18, [33] [34] [35] [36] . They include tests of GW generation [37] [38] [39] [40] [41] , where deviations in the post-Newtonian (PN) coefficients in the inspiral, and phenomeno-logical coefficients in the plunge and merger stages can be bounded; tests of GW propagation [42] , 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 [19] , for which more than two GW detectors are needed to set statistically significant bounds, and tests of the remnant properties [43] [44] [45] [46] [47] [48] [49] in the post-merger 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 [50] proposed the idea of employing BH spectroscopy [51] 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 nohair 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. [52] [53] [54] [55] ). Several analyses [44] [45] [46] [47] [48] [49] 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 [56] , or overtones -for example the (2, 2, 1) mode [47, 57, 58] . These studies have been pursued either using a superposition of damped sinusoids [46, 59] , in some cases augmemted with QNM amplitudes calibrated to numerical-relativity (NR) simulations, or with parameterized inspiral-merger-ringdown (IMR) waveform models, where the QNM frequencies are not necessarily fixed to the GR values for BHs, but kept free [45, 48] .

Here, we will employ the parameterized IMR model of Ref. [48] , constructed from a nonprecessing-spin effective-onebody (EOB) waveform model [60] [61] [62] , to carry out theory- I. Summary of the upper limits, at 90% credible level, on the length-scale of the modified gravity theories under investigation. They were obtained combining Bayesian-inference results from GW150914 and GW200129. Current constraints on th are also listed. a We notice that our result for the quartic EFT of GR is only in marginal tension with our hard cut-off scale for the validity of the theory, and hence we do still quote it here (see Secs. IV D and V D for details). specific tests of the ringdown using four high-curvature gravity theories. Previously, such parameterized waveform model was employed in Refs. [18, 36] 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 effectivefield theories (EFTs) of GR, and express the QNM frequencies using the parametrized ringdown spin-expansion coefficients (ParSpec) of Ref. [63] . 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 length-scale th of the gravity theory under consideration, the GR limit corresponding to th → 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 length-scale th of each theory. In Table I , we summarize our results, and compare them with existing constraints.

The paper is organized as follows. In Sec. II, we briefly describe the four higher-curvature modified gravity theories for which we peform the ringdown test. In Sec. III we build our parameterized IMR model for nonprecessing-spin compactobject 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 length-scales 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., precessingspins and eccentricity). In the Appendix A we provide details in calculating the non-GR QNM frequencies, when using Par-Spec, for the modified gravity theories under consideration.

Henceforth, unless otherwise specified, we work in geometric units G = 1 = c.

We will treat the modified theories of gravity as EFTs, and focus on finite-size effects (see, e.g., Ref. [66] ). Thus, for each gravity theory we impose that the fundamental length-scale th 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.

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

where g ≡ det(g µν ) is the metric determinant, R is the Ricci scalar, ϕ is a dynamical scalar field, with kinetic term (∂ϕ) 2 = g µν ∂ µ ϕ∂ ν ϕ, which couples to the Gauss-Bonnet invariant G = R µνρσ R µνρσ − 4R µν R µν + R 2 , and R µνρσ and R µν are the Riemann and Ricci tensors respectively. By itself, d 4 x √ −g G is a boundary term in four dimensions and hence does not contribute to the field equations [67] . However, when coupled to ϕ, it can contribute to the field equation through the coupling function f (ϕ). The strength of the coupling is set by GB , with dimensions of length. Different subclasses of this theory are determined by the function f (ϕ) 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 (ϕ) ≡ d f /dϕ 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 (ϕ) = exp(ϕ) [68] . In the second class, f (ϕ) = 0 can vanish for some constant ϕ 0 . In this case, the theory admits the same stationary, asymptotically flat BH solutions as GR [69] and those of scalarized BHs [69] [70] [71] [72] [73] [74] . Examples include Gaussian f (ϕ) ∝ exp(−ϕ 2 ) [70] and the quadratic f (ϕ) ∝ ϕ 2 [69] 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" [68, 75, 76] . 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. [77] [78] [79] [80] [81] ) which affects the GW phase at −1PN 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. [64] placed the bound EdGB 7.1 km using the NSBH binaries GW200105 and GW200115 [8] while Refs. [82, 83] , obtained EdGB 9.1 km, by stacking the posteriors of EdGB from a selection of BBHs from the GWTC-1 and GWTC-2 catalogs [84, 85] 1 . 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 [90, 91] (i.e., the isospectrality [92] ) of Schwarzschild BHs in GR.

Using BH perturbation theory, the QNMs for nonrotating BHs in EdGB gravity were first computed in Refs. [93, 94] and were extended (in the polar sector) to leading-order in BH spin in Ref. [95] . See Ref. [96] for a study in the geometrical optics limit ( 1) and Refs. [97, 98] for NR studies.

This theory is described by the action [99, 100] ,

where ϑ is a pseudoscalar field which couples to the Pontryagin density * RR = * R µν ρσ R νµ ρσ , where * R µνρσ is the dual of the Riemann tensor defined as * R µνρσ = µν γδ R γδρσ /2, and µνγδ 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 [99] . However, the Pontryagin density can modify the field equations when coupled to ϑ; the strength of the coupling set by 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 1 These bounds are, strictly speaking, valid only when the scalar field ϕ is small (i.e., ϕ 1), and they take into account only the leading-order scalar field interaction arising from the original dilatonic coupling (i.e., f (ϕ) ≈ ϕ in Eq. (2.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 [86] , while BHs do [87] [88] [89] . Finally, we note that the constant α EdGB used in Refs. [64, 77, 82, 83] is related to EdGB as EdGB = 4π 1/4 |α EdGB | 1/2 . is included and the Kerr metric is not a solution of the theory [99] . These rotating BHs support a scalar field which fall off as r −2 asymptotically, to which we can associate a scalar dipole charge [101, 102] and the leading-order modification to the GW phase enters at 2PN [77] . Deviations from GR at this PN order are constrained with present GW observations [18] only at the level of ∼ 50% (see the constraint on the 2PN parameter ϕ 4 in Fig. 6 of Ref. [18] ). 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 [82, 83] . Nonetheless, the theory has been constrained in Ref. [65] , which found dCS 8.5 km, by folding data from the X-ray observations of the pulsar PSR J0030+0451 [103, 104] by NICER [105, 106] and from the GW observation of the binary NS GW170817 [5, 6] , using equation-of-state independent relations between NS moment of inertia and tidal deformability [107] [108] [109] .

The QNMs of the Schwarzschild BH in dCS gravity were studied in Refs. [110] [111] [112] , 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. [113, 114] . They were also extracted from NR simulations of BH head-on collisions in Ref. [115] .

Our last example of modified gravity theories are the socalled EFTs of GR [66, [116] [117] [118] [119] [120] [121] [122] . They are described by the action

where EFT is a length-scale assumed to be small compared to the length-scale M associated with a BH (i.e., EFT /M 1), and L (2n) are corrections that introduce higher-order curvature tensors (with 2n metric derivatives).

More specifically, we follow the notation of Refs. [121, 122] and consider up to dimension-eight operators (n = 4),

where C = R µνρσ R µνρσ ,C = R µνγδ µν ρσ R ρσγδ , and both λ o,e and ε 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 λ e = λ o = 1, leaving us with cEFT as our single free parameter. Similarly, in the dimension-eight case, we set ε 1 = 1 and ε 2 = ε 3 = 0, as done in Ref. [66] . This leaves us with qEFT as our single free parameter.

For the EFT of GR with dimension-eight operators, Ref. [66] focused on the orbital effects (i.e., instead of finite-size effects) and performed Bayesian model selection using the two lowestmass 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 qEFT ∼ 150 km.

The QNMs of nonrotating BHs where calculated in Refs. [120, 121] in the dimension-six EFT and in Refs. [121, 123] in the dimension-eight EFT. Reference [121] calculated the leading-order BH spin corrections to the QNM spectra.

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 spinexpansion coefficients (ParSpec) framework [63] in Sec. III A. Next, in Sec. III B, we review our baseline parametrized IMR waveform model [45, 48] , and explain how we extend it to include the ParSpec. Finally, in Sec. III C, 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.

A general procedure to describe deviations to the QNM frequencies ω GR mn and damping times τ GR mn of BHs of GR is to write,

where δω mn and δτ mn are the fractional deformation parameters, and m are the multipole indices, and n the overtone number. This type of parametrization 1 was adopted, for instance, in Refs. [43-45, 48, 130] . The current LVK tests of BH ringdown (see Sec.VII.A of Ref. [131] or Sec.VIII.A of Ref. [18] ) take a restrictive theoryindependent approach towards the inference of δω mn and δτ mn . These deviations are either assumed to occur identically across all observed sources or belong to a generic underlying Gaussian population. A more realistic scenario, however, is that 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. 1 See also Refs. [124] [125] [126] and Refs. [127] [128] [129] for alternative parametrizations.

The ParSpec framework was introduced in Ref. [63] and can be used to our purpose. It is an observable-based bivariate expansion of Eq. (3.1), given by

where M f and χ f are the detector-frame final mass and spin, respectively; the quantities ω ( j) mn and τ ( j) mn are dimensionless coefficients of the expansion in spin for the QNMs of BHs in GR, while δω ( j) mn and δτ ( j) mn are source-independent dimensionless coefficients that characterize the corrections to the GR QNM at each spin-order, and N max is the order of the spin expansion. All source dependence due to a given modified gravity theory is contained in the dimensionless parameter γ, which reads

which depends on the length-scale parameter th of the specific gravity theory (non-GR modifications become important at distances th ), and the exponent p is related to how the non-GR modifications are added to the Einstein-Hilbert action. In Eq. (3.3) we made γ dimensionless by the length-scale associated with the remnant BH (i.e., its source-frame mass M s f 2 ), which we can also write in terms of the detector-frame mass M f through the redshift z [132] . Also, dividing by the factor G/c 2 allows us to express th in physically intuitive metric units.

In principle, a modification to GR would also affect M f and χ f and the expansion should be written in terms of the non-GR mass and spin, sayM f andχ 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 δω ( j) mn and δτ ( j) mn . This means we can identify the M f and χ f with their corresponding GR values 3 . 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 (γ → 0) the series (3.2) truncated at N max = 4, reproduces with 1% accuracy the GR QNMs for BH's spins χ f 0.7. The values of the fitting coefficients ω ( j) mn and τ ( j) mn can be found in Ref. [63] . In Ref. [49] , the fitting coefficients were calculated up to N max = 9, which extend the validity of the spin-expansion up to χ f 0.99. As we will discuss in Sec. IV, the expansion to N max = 4 is sufficient for our purposes. For convenience we list the coefficients in the case of GR in Table II .

We now describe the waveform model used in our paper to infer properties of a BBH ringdown. As in Refs. [45, 48] , 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 θ, that includes the masses and spins of the two BHs, (m 1 , m 2 , S 1 , S 2 ), the sky location determined by the luminosity distance D L , right ascension α and declination δ, and the orientation of the binary given by the inclination ι and polarization ψ angles. The set is completed by the choice of a reference time t 0 and phase φ 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 (L), we reduce the six components of the spins to just two, χ i ≡ S i ·L/m 2 i 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 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 m 2 (and thus q 1), and the symmetric mass ratio of the binary, ν = m 1 m 2 /(m 1 + m 2 ) 2 . Note that for BHs −1 χ i 1.

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

where −2 Y lm (ι, ϕ 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. [62] , 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 [133] to the multipolar spin-aligned EOB waveform model of Refs. [60, 61] (henceforth we refer to our baseline model as SEOBNR 1 ). An accurate description of the merger is incorporated through calibration with NR simulations, as described in Refs. [60, 61] , along with information for the merger and ringdown phases, from BH perturbation theory. The mergerringdown waveform, h

is the Heaviside step function. The mergerringdown waveform is expressed as an exponentially damped sinusoid [60, 61] h merger-RD lm where

are the complex frequencies of the fundamental (0-th overtone) QNMs of the remnant BH. The functionsà lm (t) andφ lm (t) are defined in Ref. [60, 61] . In the SEOBNR model [60, 61] , the complex frequencies σ m0 are computed by first determining the final mass and spin from estimates of the initial masses and spins through NR-fittingformulas [135, 136] , and then converting them to the complex frequencies using BH perturbation theory-inspired analytical fits [52, 137] . Hence,

where (ω GR m0 , τ GR 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 A (see Eqs. (3.2)). Hence,

where the mass and spin of the remnant object (M f , χ f ) are themselves functions of (m 1 , m 2 , χ 1 , χ 2 ) [135, 136] , and we fix p to a certain theory-specific value. Additionally, the frequencies depend on the ParSpec coefficients {δω ( j) m0 , δτ ( j) m0 }, and the length-scale th .

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

Let us now establish the connection between the theoryindependent 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 ( , m, n) = (2, 2, 0). For this reason, for simplicity, we omit the subscripts hereafter and rewrite ω mn and τ mn , given by Eqs. (3.2) as,

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. We summarize each theory we consider together with: the exponent p at which their QNM-modification enters, the corresponding modifications to the oscillation frequency δω (n) 220 and decay time δτ (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. [63] . The remaining GR coefficients for 1 < j 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. [63] .

How can we determine the beyond-GR corrections? In GR, comparison between the numerically determined Kerr QNMs against the fitting formula (3.10) fixes the GR expansion coefficients ω ( j) and τ ( 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 length-scale th , the latter being specific to each theory, up to j = 1 in the spin expansion (see Table II ). The idea is then to compare these formulas against Eq. (3.10) to fix p, δω ( j) , and δτ ( 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 theoryspecific QNMs has only th as a free beyond-GR parameter. We emphasize that our procedure is different from that of Ref. [49] which, for a given value of p, varied all th , δω ( j) , and δτ ( j) parameters, and then proceeded to use the posteriors on 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 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. (3.10)? To answer this question one has to work with a chosen theory and perform a translation between the metric perturbations h µν in the Regge-Wheeler-Zerilli gauge [138, 139] and connect it with the transverse-traceless gauge used to described GWs (see, for instance, Ref. [140] , 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. [141] for a counterexample). Thus, how axial and polar gravitational QNMs appear in the GW signal has to be answered on a theoryby-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.

The justification is that this is the mode (if excited during the merger) which is the most likely to appear in the signal. 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 II and leave the details of our calculations to Appendix A. 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 leadingorder j = 0 deformations to the fundamental QNM (see Table II ). We choose binary parameters similar to GW150914, but non-spinning, with detector-frame masses m 1 = 39 M and m 2 = 31 M . The top panel shows the + GW polarization h + in both theories, while the bottom panel shows the amplitude |h| = (h 2 + + h 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 match defined in Sec. III B, 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 δω (0) coefficient in this theory.

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.

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, λ, given the hypothesis H, using Bayes' theorem,

where P(λ|d, H) is the posterior probability distribution, p(λ|H) the prior, L(d|λ, H) the likelihood, and E(d|H) the evidence. The set of parameters, λ is a union of the GR waveform model parameters θ (see Sec. III B) and th , the only non-GR parameter in this problem which, we recall, sets the characteristic length-scale 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,

with the noise-weighted inner product ·|· defined as,

whereÃ( 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 low , f 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 [142, 143] , as part of the LALInference software suite [134, 144] . We subsequently marginalize over the remaining parameters to obtain the posterior probability distribution function (PDF) on th (i.e., P( th |d, H)), our main parameter of interest. For N independent GW observations {d j }, j = 1, ..., N, each characterized by a PDF P j ( th |d j , H), the joint posterior can be written as:

where p j ( th |H) are the priors used for each observation, p( th |H) is an overall prior, and we assume that the value of th is shared among all events. Since we assume a uniform prior on th , the joint posterior is equal to the joint likelihood.

Hereafter, we will drop the explicit usage of H.

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 ) ∝ D 2 L , and isotropic on the binary orientation, p(ι, ψ, φ 0 ) ∝ sin ι. For the spins (χ 1 , χ 2 ), we assume a prior uniform and isotropic in the spin magnitudes 1 .

Among our non-GR parameters { th , δω ( j) δτ ( j) }, as already mentioned in the previous section, we hold {δω ( j) , δτ ( j) } fixed to theory-specific predictions, and only allow th to vary freely. We assume an uniform prior on th between appropriate ranges. The lower limit is set by the fact that the modified gravity theories we consider all have p even and hence we can assume th > 0 without loss of generality.

The pSEOBNR model, as described in Sec. III B, 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 post-merger part of the signal. The latest LVK GW catalog [4] reported 90 observed signals not all of which are relevant for our BH ringdown analysis. In fact, in the accompanying paper [18] on tests of GR, the pSEOBNRv4HM [45, 48] analysis 2 , 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 [1] and GW200129 [4] . These two events, with a total (source-frame) mass of 65M and 63.4M 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 (mergerringdown) SNRs which are both ≈ 16 (see the columns for ρ post-insp in Table III of Ref. [35] and Table IV of Ref. [18] ). 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. [18] . 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. [18] .

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 th .

There are two conditions that we must verify before we can confidently claim to have placed a constraint on 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, a length-scale. As a cut-off length-scale for the validity of the EFT we use,

where ε is a dimensionless number and m is the median value of one of the mass scales involved in the problem. We note that Λ EFT has dimensions of length and hence can be compared to each theory's fundamental length-scale th . Here we explore the range ε ∈ [0, 1], but following Refs. [64, 82, 83] we quote our final results using ε = 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 th , if most of the PDF P( th |d) support is in the interval [0, Λ EFT (1/2, m)]. In practice, this can be quantified through the cumulative distribution function (CDF) associated with the marginalized posterior distribution P( th |d), namely For instance, we require that for a bound at 90% credible level to be placed on th that

where we let max th = Λ EFT in Eq. (4.6), and likewise for other credibility percentiles.

Second, as already emphasized in Ref. [63] , the ParSpec formalism is by construction perturbative. This means that the non-GR deformation parameters are small, that is, γ δω ( j) 1, and γ δτ ( j) 1, (ParSpec bound), (4.8) for all orders j in the expansion in dimensionless spin χ f and where γ was defined in Eq. (3.3). 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 m that we should use in Eq. (4.5)? In Refs. [64, 82, 83] , 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 cut-off scale Λ 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 f to compute Λ 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 m = m 2 should be used. Here we adopt a pragmatic approach to this issue and consider both masses, m 2 and M f , to determine Λ EFT . We then compare how different assumptions yield to different interpretations of the results of our parameter estimation.

We start with EdGB gravity. In Fig. 2 we show the marginalized PDFs of the coupling constant 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 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 γ, which here reads,

As shown in Fig. 3 , γ 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 max = 1) displaces the posteriors distributions towards smaller values of EdGB (see Fig. 2 ), and larger values of γ EdGB (see Fig. 3 ).

As we have emphasized in Sec. IV D, we must first check whether the "EFT" (4.7) and "ParSpec" Fig. 4 and in the main text, our analysis of these events fails to satisfy the EFT bound (4.7). We mark with vertical lines the 90% upper credible intervals.

regardless of the mass scale m used and even at ε = 1, at which the EFT description of the theory would not be valid anyway. This shows that that the "EFT bound" given by Eq. (4.7) is never met to a significant credible level and thus that we cannot place a bound on EdGB . The situation is similar for GW200129 with N max = 0 and does not change for either event when we add spin corrections to the EdGB QNM. For the case with N max = 1, we find that the "EFT bound" is satisfied only for ε ≈ 0.8 and ≈ 1 for GW150914 and GW200129, respectively. However, we set the maximum value of ε 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 III . We can compare this conclusion with that of Ref. [49] , which found that p = 4 modifications (such as the case of EdGB gravity) are constrained to 35 km, but not including theoryspecific QNM information on δω ( j) and δτ ( j) . Furthermore, Ref. [49] 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 contrasts our results with those of Refs. [64, 82, 83] which relied on the BBH inspiral to constrain EdGB , as discussed in Sec. II A. 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. [94] , the QNMs of EdGB BHs only differ slightly from their Schwarzschild counterparts. Second, the larger mass M f of the remnant BH, suppresses scalar field's charge relatively to the initial binary components.

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

does indeed peak at zero indicating consistency with GR, similarly to what is shown in Fig. 3 for γ EdGB . We also see that in both cases the inclusion of leading-order-in-spin correction In Fig. 6 we show the CDF for GW150914, we see that with m = m 2 , Eq. (4.7) is not satisfied unless ε ≈ 0.9 (with only j = 0 corrections) and ε ≈ 0.7 (with both j = 0 and 1 corrections). The situation is different if we use m = M f . In this case, we find that with or without spin corrections Eq. (4.7) can be satisfied with ε 1/2 (i.e., below the criteria used Refs. [64, 82, 83] ). This means that with our model's assumptions and using the remnant's source mass M f to set the cut-off scale that we can claim an upper bound dCS 41.9 km at 90% credible level,

on dCS gravity and would constitute the first bound on this theory with GW observations alone. We can draw qualitatively similar conclusions from the GW200129 event. In particular, we find, Fig. 2 , but for dCS gravity. We stress that posteriors obtained including n = 1 corrections, violate conditions (4.8) 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 (4.7) and (4.8) simultaneously into consideration when analyzing the results of the parameter estimation.

GW150914). We also found for both GW events, that the perturbative-conditions (4.8) required by the ParSpec is violated for the γ dCS δτ (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 max = 0 bound above. Finally, since both events individually lead to a bound on dCS (assuming a cut-off scale for m = M f and N max = 0), we can combine the posteriors to obtain the cumulative bound, dCS 38.7 km at 90% credible level, (5.5) which is the main result of this section. This bound is approximately a factor of four weaker than that placed by Ref. [65] , 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 [64, 82, 83] .

In Table IV we summarize our findings of this section.

As an additional check, to verify the robustness of our constraint, we show in Fig. 7 , the final spin χ f and remnant mass M 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 A). In fact, we observed no bias on the estimation of M f and χ f for all theories considered here. For completness, in Appendix B we also show how all other intrinsic parameters remain unbiased. TABLE IV . Detailed summary of our results for dCS gravity for GW150914, GW200129, and combined events using m = M f , ε = 1/2 and quoting only 90% credible bounds. We found that while our posteriors satisfy the condition (4.7) (with ε = 1/2), they do not obey the condition (4.8) for N max = 1. This means that our results for N max = 0 are the only ones we can confidently quote. The combined bound, which is also quoted in Table II , is dCS 38.7 km at 90% credible level.

We now consider the cubic EFT of GR. In Fig. 8 we show the marginalized posterior distributions functions of cEFT for GW150914 (top panel) and GW200129 (middle panel), with different curve colors corresponding to different N max in the spin expansion. We find that in this theory, the posterior distributions are mostly uniform for cEFT 40 km (contrast this with the EdGB and dCS gravity cases in Figs. 2 and 5) . For values cEFT 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 65 70 GW200129 (see bottom panel in Fig. 9 ). Moreover, we find that the EFT (4.7) and ParSpec (4.8) bounds are satisfied for both events both when m = M f , ε = 1/2, and N max = 0. This allows us to place the combined bound of cEFT 38.2 km, at 90% credible level .

(5.6)

As also happened for our study for dCS, the find that, for the cubic EFT, the ParSpec bound is violated by the N 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 V 

Let us now consider the quartic EFT of GR, as our final example. In Fig. 10 for GW150914 (top panel), GW200129 (middle panel) for N 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 ε ≈ 0.58 (in the case of GW15094) and for ε ≈ 0.64 (in the case of GW200129). Having in mind that the cut off ε = 1/2 is not fundamental, but to keep consistency across our analysis, our final result qEFT 51.3 km, (5.7) 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.

In this theory, we have considered only N max = 0. We find that the addition of spin corrections (while maintaining the same prior ranges on qEFT as used in the N 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 max 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 max = 1. Table VI summarizes our findings for the quartic EFT of GR.

We presented an unified framework that combines the Par-Spec framework to model deviations to the GR QNMs [63] with the pSEOBNR waveform model [45, 48] . 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 I. In particular, we found, that within the interpretation of these theories as EFTs and the region of validity of the ParSpec framework, the fundamental length-scale of dCS gravity is bound as dCS 34.5 km, at 90% credible level, when stacking the posteriors of GW150914 and GW200129. This is the first GW-alone constrain on this theory. In contrast, we could not place any bounds on the fundamental length-scale of EdGB gravity EdGB . This dichotomy between the two theories has a counterpart with works that considered the inspiral part of the GW signal alone [64, 82, 83] . 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 theoryspecific 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. [49] extending the validity of the ParSpec formulas up to spins of χ f ≈ 0.99. For the events analyzed here, the original fit by Ref. [63] 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. [145, 146] ). 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 [41] , as done in Ref. [66] or the TIGER code [39, 40] . However, setting bounds on deviations from GR caused by orbital effects requires a different EFT interpretation than what we adopted in Sec. IV D (see also Sec. IIC in Ref. [66] ). Indeed, in this case one would need to analyze the data considering that the modified theory of gravity is valid for th M, but th R, being 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 [97, [147] [148] [149] [150] [151] [152] [153] [154] . 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.

and to Steffen Grunewald for assistance. 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 II . We consider only the fundamental QNM ( , 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.

We start by considering EdGB gravity and focus on Refs. [94, 95] to determine the ParSpec coefficients for this theory. In particular, Ref. [94] found that the damping time of the dominant axial gravitational-led mode increases as the length-scale EdGB is increased. The leading-order spin corrections to the polar-parity QNMs was studied in Ref. [95] . Hence, according to the prescription of Sec. III C, we select the axial-parity branch of QNMs. For the nonrotating QNMs we use the numerical data of Ref. [94] and generate a new linear fit in γ EdGB using numerical QNM data valid for small values of the γ EdGB (see, in particular, Eq. (27) and Fig. 1 of Ref. [94] ). We find,

The small values of the numerical prefactors of γ 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 γ EdGB ≈ 0.3. The spin-corrections to the polar gravitational-led modes were calculated in Ref. [95] (see, in particular, their Eqs. (51) and (52)). For consistency with our previous discussion, we truncate these equations at leading-order in γ 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. [95] are our best presently available guide.

We can expand the resulting formula in γ EdGB and the coefficients δω (i) and δτ (i) , i = 1, 2 can be read-off by comparison against Eqs. (3.10) , where for the damping time we use the relation Im(σ) EdGB = −1/τ EdGB and re-expand in EdGB and χ. These steps yield for p EdGB = 4, for the j = 1 coefficients.

For dCS gravity, we follow Ref. [114] , 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 dCS , at constant, small BH spin. Hence, according to the prescription of Sec. III C, this is the branch of QNMs we choose to work with.

We then proceed to determine δω ( j) and δτ ( j) as follows. Using the fitting formula Eq. (54a) of Ref. [114] , namely, MRe(σ) dCS = c 1 + c 2 κζ + (c 3 + c 4 κζ) (1 − χ f ) c 5 +c 6 κζ , (A4) and similarly for the imaginary part, Im(σ) dCS = −1/τ dCS . Here κ = 1/(16π), ζ = 4 dCS /(M 4 s κ), thus κζ = γ dCS and where c i (with i = 1, . . . , 6) are fitting coefficients which can be found in Table II of Ref. [114] ,

We now expand Eq. (A4) to leading orders in χ and γ dCS , and gather the terms proportional to γ dCS . We obtain Mω dCS = (0.3722 + 1.1945γ dCS ) + (0.1861 + 5.1828γ dCS ) χ,

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 ω (0) of Ref. [63] to 0.5% relative error. The same estimate leads to a larger relative error (≈ 20%) for the linear-in-spin coefficient (i.e., 0.1861 in comparison to 0.1258 of Ref. [63] ). We attribute this difference to Ref. [114] having fitted Eq. (A4) to QNM data computed to linear-order in spin, whereas [63] fitted Eq. (3.2) to Kerr QNM valid to all orders in spin.

We can now isolate the dCS corrections from Eq. (A5) and compare against Eq. (3.10), to find p dCS = 4,

We can carry the same steps for τ dCS = −1/Im(σ) dCS and find δτ (0) dCS = 6.3619 , δτ (1) dCS = 794.66,

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 δω (1) dCS and δτ (1) dCS are compensated by the assumptions that γ dCS and χ are much less than unity, which are indeed the assumptions used in Ref. [114] to compute the QNMs.

The QNMs of slowly-rotating BHs in both cubic and quartic EFT of GR where calculated in Ref. [122] . For the cubic EFT, we use their Eq. (67) , in the particular case of λ e = λ o = 1. We then linearize the resulting expression in χ 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. (3.10) results in p cEFT = 4, δω (0) cEFT = −0.5813, δτ (0) cEFT = 2.6469, δω (1) cEFT = −3.8620, δτ (1) cEFT = 265.12,

for this theory. We proceed in the same away for the quartic EFT. Here we use Eq. (68) (with 1 = 1 and 2 = 0) and Eq. (70) of Ref. [122] , 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 qEFT = 6, 

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 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 f and χ 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 f and χ 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 χ f has most support around ≈ 0.7 justifies our use of the fitting coefficients in the ParSpec formulas of Ref. [63] . . Corner plot showing that the inferred source binary parameters for GW150914, using the same waveform model pSEOBNR, but with (purple contours) and without setting the non-GR parameters different from zero (blue contours), having as an illustrative example dCS gravity and N max = 0. Here, χ eff is the dimensionless effective-spin parameter, related to the individual spins χ i and masses m i of each binary component as χ eff ≡ (m 1 χ 1 + m 2 χ 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.

Observation of Gravitational Waves from a Binary Black Hole Merger

Advanced LIGO

Advanced Virgo: a secondgeneration interferometric gravitational wave detector

GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run

GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral

GW170817: Measurements of neutron star radii and equation of state

GW190425: Observation of a Compact Binary Coalescence with Total Mass ∼ 3.4M

Observation of Gravitational Waves from Two Neutron Star-Black Hole Coalescences

1-OGC: The first open gravitational-wave catalog of binary mergers from analysis of public Advanced LIGO data

2-OGC: Open Gravitational-wave Catalog of binary mergers from analysis of public Advanced LIGO and Virgo data

New binary black hole mergers in the second observing run of Advanced LIGO and Advanced Virgo

Detecting Gravitational Waves With Disparate Detector Responses: Two New Binary Black Hole Mergers

4-OGC: Catalog of gravitational waves from compact-binary mergers

New binary black hole mergers in the LIGO-Virgo O3a data

Constraints on the cosmic expansion history from GWTC-3

The population of merging compact binaries inferred using gravitational waves through GWTC-3

Theoretical Physics Implications of the Binary Black-Hole Mergers GW150914 and GW151226

The Confrontation between General Relativity and Experiment

Testing Relativistic Gravity with Radio Pulsars

Strong-Field Gravity Tests with the Double Pulsar

Detection of the gravitational redshift in the orbit of the star S2 near the Galactic centre massive black hole

Relativistic redshift of the star S0-2 orbiting the Galactic center supermassive black hole

First M87 Event Horizon Telescope Results. IV. Imaging the Central Supermassive Black Hole

Modified Gravity and Cosmology

Stability of the Schwarzschild metric

Long Wave Trains of Gravitational Waves from a Vibrating Black Hole

The quasi-normal modes of the Schwarzschild black hole

Axisymmetric Black Hole Has Only Two Degrees of Freedom

Event horizons in static vacuum space-times

Black holes in general relativity

Uniqueness of the Kerr black hole

Tests of general relativity with GW150914

Tests of General Relativity with GW170817

Tests of General Relativity with the Binary Black Hole Signals from the LIGO-Virgo Catalog GWTC-1

Tests of General Relativity with Binary Black Holes from the second LIGO-Virgo Gravitational-Wave Transient Catalog

Testing post-Newtonian theory with gravitational wave observations

Fundamental Theoretical Bias in Gravitational Wave Astrophysics and the Parameterized Post-Einsteinian Framework

Towards a generic test of the strong field dynamics of general relativity using compact binary coalescence

TIGER: A data analysis pipeline for testing the strong-field dynamics of general relativity with gravitational wave signals from coalescing compact binaries

Tests of General Relativity with Gravitational-Wave Observations using a Flexible-Theory-Independent Method

Bounding the mass of the graviton using gravitational wave observations of inspiralling compact binaries

Testing the no-hair theorem with black hole ringdowns using TIGER

Empirical tests of the black hole no-hair conjecture using gravitational-wave observations

Black-hole Spectroscopy by Making Full Use of Gravitational-Wave Modeling

Observational Black Hole Spectroscopy: A time-domain multimode analysis of GW150914

Testing the no-hair theorem with GW150914

Constraints on quasinormal-mode frequencies with LIGO-Virgo binary-blackhole observations

Enhancing modified gravity detection from gravitational-wave observations using the parametrized ringdown spin expansion coeffcients formalism

Black hole spectroscopy: Testing general relativity through gravitational wave observations

Black holes and gravitational waves. III. The resonant frequencies of rotating holes

On gravitational-wave spectroscopy of massive black holes with the space interferometer LISA

LISA parameter estimation and source localization with higher harmonics of the ringdown

The landscape of massive black-hole spectroscopy with LISA and Einstein Telescope

Black hole spectroscopy horizons for current and future gravitational wave detectors

Observation of a multimode quasi-normal spectrum from a perturbed black hole

On the detection of ringdown overtones in GW150914

Revisiting the ringdown of GW150914

Black Hole Ringdown: The Importance of Overtones

Improved effective-one-body model of spinning, nonprecessing binary black holes for the era of gravitational-wave astrophysics with advanced detectors

Enriching the Symphony of Gravitational Waves from Binary Black Holes by Tuning Higher Harmonics

Fast post-adiabatic waveforms in the time domain: Applications to compact binary coalescences in LIGO and Virgo

Parametrized ringdown spin expansion coefficients: a data-analysis framework for black-hole spectroscopy with multiple events

Constraints on Einsteindilation-Gauss-Bonnet gravity from black hole-neutron star gravitational wave events

Astrophysical and theoretical physics implications from multimessenger neutron star observations

Gravitational-Wave Constraints on an Effective Field-Theory Extension of General Relativity

Higher Derivative Gravity, Surface Terms and String Theory

Dilatonic black holes in higher curvature string gravity

Spontaneous scalarization of black holes and compact stars from a Gauss-Bonnet coupling

New Gauss-Bonnet Black Holes with Curvature-Induced Scalarization in Extended Scalar-Tensor Theories

Self-interactions and Spontaneous Black Hole Scalarization

Spininduced black hole spontaneous scalarization

Spin-induced scalarized black holes

Spin-induced black-hole scalarization in Einstein-scalar-Gauss-Bonnet theory

Quantum hair on black holes

Asymptotically flat black holes with scalar hair: a review

Post-Newtonian, Quasi-Circular Binary Inspirals in Quadratic Modified Gravity

Post-Newtonian dynamics and black hole thermodynamics in Einstein-scalar-Gauss-Bonnet gravity

Nonlinear curvature effects in gravitational waves from inspiralling black hole binaries

Post-Newtonian gravitational and scalar waves in scalar-Gauss-Bonnet gravity

Black hole sensitivities in Einstein-scalar-Gauss-Bonnet gravity

Fundamental Physics Implications for Higher-Curvature Theories from Binary Black Hole Signals in the LIGO-Virgo Catalog GWTC-1

Improved gravitational-wave constraints on higher-order curvature theories of gravity

GWTC-1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs

GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run

Challenging the Presence of Scalar Charge and Dipolar Radiation in Binary Pulsars

Non-Spinning Black Holes in Alternative Theories of Gravity

Black hole hair in generalized scalar-tensor gravity

Black hole hair in generalized scalar-tensor gravity: An explicit example

Darboux transformation in black hole perturbation theory

Darboux covariance: A hidden symmetry of perturbed Schwarzschild black holes

The mathematical theory of black holes

Are black holes in alternative theories serious astrophysical candidates? The Case for Einstein-Dilaton-Gauss-Bonnet black holes

Perturbed black holes in Einstein-dilaton-Gauss-Bonnet gravity: Stability, ringdown, and gravitational-wave emission

Quasi-normal modes of rotating black holes in Einstein-dilaton Gauss-Bonnet gravity: the first order in rotation

Eikonal quasinormal modes of black holes beyond general relativity

Scalar Gauss-Bonnet gravity

Black holes and binary mergers in scalar Gauss-Bonnet gravity: scalar field dynamics

Stability of Rotating Black Holes in Einstein Dilaton Gauss-Bonnet Gravity

Chern-Simons modification of general relativity

Chern-Simons Modified General Relativity

Dynamical Chern-Simons Modified Gravity. I. Spinning Black Holes in the Slow-Rotation Approximation

Rotating black hole in extended Chern-Simons modified gravity

New pulsars from an Arecibo drift scan search

The NANOGrav 11-year Data Set: High-precision timing of 45 Millisecond Pulsars

A NICER View of PSR J0030+0451: Millisecond Pulsar Parameter Estimation

PSR J0030+0451 Mass and Radius from NICER Data and Implications for the Properties of Neutron Star Matter

I-Love-Q

I-Love-Q Relations in Neutron Stars and their Applications to Astrophysics, Gravitational Waves and Fundamental Physics

I-Love-Q Relations for Neutron Stars in dynamical Chern Simons Gravity

Perturbations of Schwarzschild Black Holes in Chern-Simons Modified Gravity

Perturbations of Schwarzschild black holes in Dynamical Chern-Simons modified gravity

Gravitational signature of Schwarzschild black holes in dynamical Chern-Simons gravity

Analytical computation of quasinormal modes of slowly rotating black holes in dynamical Chern-Simons gravity

Quasinormal modes of slowly-rotating black holes in dynamical Chern-Simons gravity

Numerical binary black hole collisions in dynamical Chern-Simons gravity

An effective formalism for testing extensions to General Relativity with gravitational waves

Leading higher-derivative corrections to Kerr geometry

On higher-derivative effects on the gravitational potential and particle bending

From amplitudes to gravitational radiation with cubic interactions and tidal effects

Black Hole Gravitational Waves in the Effective Field Theory of Gravity

Ringing of rotating black holes in higher-derivative gravity

Gravitational ringing of rotating black holes in higher-derivative gravity

Black Holes in an Effective Field Theory Extension of General Relativity

Post-Kerr black hole spectroscopy

Eikonal quasinormal modes of black holes beyond General Relativity

Eikonal quasinormal modes of black holes beyond general relativity. II. Generalized scalartensor perturbations

Parametrized black hole quasinormal ringdown: Decoupled equations for nonrotating black holes

Parametrized black hole quasinormal ringdown. II. Coupled equations and quadratic corrections for nonrotating black holes

Note on the parametrized black hole quasinormal ringdown formalism

Bayesian model selection for testing the no-hair theorem with black hole ringdowns

Tests of general relativity with binary black holes from the second LIGO-Virgo gravitational-wave transient catalog

Coalescing binaries -Probe of the universe

Efficient effective one body timedomain gravitational waveforms

Effective-one-body model for black-hole binaries with generic mass ratios and spins

The final spin from binary black holes in quasi-circular orbits

Quasinormal modes of black holes and black branes

Stability of a Schwarzschild singularity

Effective potential for even parity Regge-Wheeler gravitational perturbation equations

Astrophysics and Cosmology

Effective Field Theory for the perturbations of a slowly rotating black hole

Bayesian inference on compact binary inspiral gravitational radiation signals in interferometric data

Parameter estimation of spinning binary inspirals using Markovchain Monte Carlo

Parameter estimation for compact binaries with ground-based gravitational-wave observations using the LAL-Inference software library

Multipolar Effective-One-Body Waveforms for Precessing Binary Black Holes: Construction and Validation

Effective-one-body multipolar waveforms for eccentric binary black holes with nonprecessing spins

Late Inspiral and Merger of Binary Black Holes in Scalar-Tensor Theories of Gravity

Numerical binary black hole mergers in dynamical Chern-Simons gravity: Scalar field

Black Hole Dynamics in Einstein-Maxwell-Dilaton Theory

Numerical relativity simulation of GW150914 in Einstein dilaton Gauss-Bonnet gravity

Dynamical Descalarization in Binary Black Hole Mergers

Evolution of Einstein-scalar-Gauss-Bonnet gravity using a modified harmonic formulation

Dynamics of Spontaneous Black Hole Scalarization and Mergers in Einstein-Scalar-Gauss-Bonnet Gravity

Black Hole Binaries in Cubic Horndeski Theories

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 are grateful for the computational resources provided by the Max Planck Institute for Gravitational Physics in Potsdam, specifically the high-performace computing cluster Hypatia