The gravitational waves emitted in the ringdown phase of binary black hole coalescence are a unique probe of strong gravity. At late times in the ringdown these waves can be described by quasinormal modes whose frequencies encode the mass and spin of the remnant as well as the theory of gravity in play. Understanding precisely how deviations from general relativity affect the quasinormal mode frequencies of ringing black holes however is extremely challenging as it requires solving highly coupled and sometimes higherorder partial differential equations. We here extend a novel approach metric perturbations with spectral methods METRICS to study the gravitational metric perturbations and the quasinormal mode frequencies of ringing black holes in modified gravity. We first derive the asymptotic behavior of gravitational perturbations at the event horizon and spatial infinity for rotating black holes beyond general relativity. We then extend the eigenvalueperturbation theory approach of METRICS to allow us to compute the leadingorder beyond general relativity corrections to the quasinormalmode frequencies and metric perturbations. As an example we apply METRICS to black holes with moderate spins in scalarGaussBonnet gravity. Without decoupling or simplifying the linearized field equations in this theory we compute the leadingorder corrections to the quasinormal frequencies of the axial and polar perturbations of the inlineformulammlmath displayinlinemmlminmmlmimmlmilmmlmimmlmimmmlmimmlmommlmommlmn022mmlmnmmlmathinlineformula 021 and 033 modes of black holes with dimensionless spin inlineformulammlmath displayinlinemmlmiammlmimmlmommlmommlmn0.85mmlmnmmlmathinlineformula. The numerical accuracy of the METRICS frequencies is inlineformulammlmath displayinlinemmlmo formprefixmmlmommlmsupmmlmn10mmlmnmmlmrowmmlmommlmommlmn5mmlmnmmlmrowmmlmsupmmlmathinlineformula when inlineformulammlmath displayinlinemmlmiammlmimmlmommlmommlmn0.6mmlmnmmlmathinlineformula inlineformulammlmath displayinlinemmlmommlmommlmsupmmlmn10mmlmnmmlmrowmmlmommlmommlmn4mmlmnmmlmrowmmlmsupmmlmathinlineformula when inlineformulammlmath displayinlinemmlmn0.6mmlmnmmlmoltmmlmommlmiammlmimmlmommlmommlmn0.7mmlmnmmlmathinlineformula and inlineformulammlmath displayinlinemmlmommlmommlmsupmmlmn10mmlmnmmlmrowmmlmommlmommlmn2mmlmnmmlmrowmmlmsupmmlmathinlineformula when inlineformulammlmath displayinlinemmlmn0.7mmlmnmmlmoltmmlmommlmiammlmimmlmommlmommlmn0.85mmlmnmmlmathinlineformula for all modes studied. We fit the frequencies with a polynomial in spin whose coefficients up to second order in spin are consistent with those obtained in previous slowrotating approximations. These results are the first accurate computations of the gravitational quasinormalmode frequencies of rapidly rotating black holes of inlineformulammlmath displayinlinemmlmiammlmimmlmommlmommlmn0.85mmlmnmmlmathinlineformula in scalarGaussBonnet gravity.
['General relativity', 'alternative theories of gravity', 'General Relativity and Quantum Cosmology', 'Astrophysics - Cosmology and Nongalactic Astrophysics', 'Astrophysics - High Energy Astrophysical Phenomena', 'Astrophysics - Instrumentation and Methods for Astrophysics']
Quasinormal mode frequencies and gravitational perturbations of black holes with any subextremal spin in modified gravity through METRICS The scalarGaussBonnet gravity case
{'Quasinormal mode frequencies and gravitational perturbations of black holes with any subextremal spin in modified gravity through METRICS: The scalar-Gauss-Bonnet gravity case': 'Adrian Ka-Wai Chung 1, ∗ and Nicol´as Yunes 1  \n1 Illinois Center for Advanced Studies of the Universe & Department of Physics, University of Illinois Urbana-Champaign, Urbana, Illinois 61801, USA (Dated: September 10, 2024)  \nThe gravitational waves emitted in the ringdown phase of binary black hole coalescence are a unique probe of strong gravity. At late times in the ringdown, these waves can be described by quasinormal modes, whose frequencies encode the mass and spin of the remnant, as well as the theory of gravity in play. Understanding precisely how deviations from general relativity affect the quasinormal mode frequencies of ringing black holes, however, is extremely challenging, as it requires solving highly coupled and sometimes higher-order partial differential equations. We here extend a novel approach, Metric pErTuRbations wIth speCtral methodS (METRICS), to study the gravitational metric perturbations and the quasinormal mode frequencies of ringing black holes in modified gravity. We first derive the asymptotic behavior of gravitational perturbations at the event horizon and spatial infinity for rotating black holes beyond general relativity. We then extend the eigenvalue perturbation theory approach of METRICS to allow us to compute the leading-order beyond general relativity corrections to the quasinormal-mode frequencies and metric perturbations. As an example, we apply METRICS to black holes with moderate spins in scalar-Gauss-Bonnet gravity. Without decoupling or simplifying the linearized field equations in this theory, we compute the leading-order corrections to the quasinormal frequencies of the axial and polar perturbations of the nlm = 022, 021 and 033 modes of black holes with dimensionless spin a ≤ 0 . 85. The numerical accuracy of the METRICS frequencies is ≤ 10 -5 when a ≤ 0 . 6, ≲ 10 -4 when 0 . 6 < a ≤ 0 . 7, and ≲ 10 -2 when 0 . 7 < a ≤ 0 . 85 for all modes studied. We fit the frequencies with a polynomial in spin, whose coefficients (up to second order in spin) are consistent with those obtained in previous slow-rotating approximations. These results are the first accurate computations of the gravitational quasinormal-mode frequencies of rapidly rotating black holes (of a ∼ 0 . 85) in scalar-Gauss-Bonnet gravity.', 'I. INTRODUCTION': "The detection of gravitational waves (GWs) allows us to probe the rich physics of the Universe, such as that in play during binary black hole (BBH) coalescence [113]. Broadly speaking, a BBH coalescence consists of three phases. The first phase is the inspiral, during which the BHs orbit around each other with decreasing orbital radius, as a result of GW emission. The second phase is the merger, during which the BHs merge to form a single remnant BH. The third phase is the ringdown, during which the newly formed, dynamical, and distorted remnant BH relaxes into a stationary configuration by emitting GWs that are characterized by a discrete set of complex quasinormal-mode (QNM) frequencies.  \nThe validity of general relativity (GR) in the strongfield regime is one aspect about our Universe that GWs can probe. GR has passed numerous experimental tests [7, 11, 14-23] and, thus far, represents our best understanding of spacetime and gravity. Nonetheless, GR exhibits theoretical and observational anomalies. Theoretically, GR predicts that gravitational collapse inevitably leads to the formation of a spacetime singularity, where GR cannot further describe nature. Observationally, GR  \nmay require additional parity-violating physics in the early Universe to explain the matter-antimatter asymmetry [24-27], a cosmological constant that is finely tuned [28, 29] to explain the acceleration of the Universe at late times [30, 31], and cold dark matter to describe the rotation curves of galaxies [32, 33]. These anomalies have led to the proposal and development of modified gravity theories, each of which amends different aspects of GR in an attempt to remedy these (and other) anomalies. To avoid repetition, we refer the reader to [34, 35] for a survey of the motivation of modified gravity theories and relevant GW tests.  \nThe coalescence of BBHs is a powerful laboratory in which to test these modified theories. In many such theories, BHs (or compact objects more generally) still exist, but their spacetime geometry may be different from that of their GR counterpart. Moreover, the field equations that describe the dynamics of these compact objects and the behavior of metric perturbations are generically also different from that which arise from the Einstein equations. These differences force the GWs emitted during BBH coalescence in these theories to also be different from those in GR. By comparing GW observations against GW predictions in GR, we can then test the validity of the latter and probe for GR deviations in the data directly.  \nThe ringdown phase of BBH coalescence has unique features that make it ideal for testing GR. First, the GWs  \nemitted in the ringdown phase usually have the highest frequencies among all three phases. Loosely speaking, the GW frequency gives a measure of the relative rate of change in perturbations of spacetime. By this standard, the spacetime changes most rapidly during the ringdown phase and might activate effects that could not otherwise be observed in other phases. Second, the QNM frequencies depend only on the properties of the remnant BH and are progenitor-independent, which makes them easy to characterize. Third, the QNM response of a BH involves perturbations of the BH geometry near the BH's horizon, the edge of the observable spacetime, where the field strength is the strongest, and where nature is perhaps most likely to be different from GR predictions. All these unique features make the ringdown phase a unique probe of the validity of GR in the strong-field regime.  \nHowever, many ringdown tests of gravity have mostly been model agnostic [22, 36-42]. This is understandable because computing the QNM frequencies of a BH in a particular gravity theory has been extremely challenging; this is because of the computational complexity required to solve the dynamical evolution equation of metric perturbations in these theories, which is typically a set of highly coupled partial differential equations. These computations are even more challenging when one takes into account BH backgrounds that do not spin slowly, as is the case for BH remnants of BBH coalescence.  \nTwo approaches have been recently developed to deal with these difficulties. The first approach works with curvature perturbations, expressing first the field equations and the Bianchi identities in terms of spinor coefficients, Weyl scalars, and differential operators. Then, the field equations in terms of these variables are linearized and, if the background spacetime is of Petrov-type D, then the linearized field equations can be simplified into a single master equation, known as the Teukolsky equation, for a certain Teukolsky master function [43-46]. Very recently, this formalism has been extended to modified gravity theories with leading-order deviations from GR [47-53], whose spacetime is not Petrov-type D. In this modified Teukolsky formalism, one then solves the modified Teukolsky equation, subject to appropriate boundary conditions, to compute the QNM frequencies of the modified Teukolsky function, which can be related to the frequencies of the emitted GWs at future null infinity. If one wants to reconstruct the metric perturbations from the Teukolsky master function everywhere outside the perturbed BH, one has to use a Hertz potential in the so-called ChrzanowskiKogen-Kegeles (CKK) approach [54-59], which amounts to solving ∼ 10 partial differential equations.  \nThe second approach is the Metric pErTuRbations wIth speCtral methodS (METRICS) that we developed in [34, 35]. This approach works directly with metric perturbations and perturbations of any other fields that may be present in the theory. This approach begins by calculating the asymptotic behavior of the field perturbations at the event horizon and at spatial infinity. Using this asymptotic behavior, one can then construct an asymp-  \ntotic factor that regulates the divergent behavior of the field perturbations at the event horizon and spatial infinity. The metric perturbations can be power decomposed through the product of the asymptotic factor and a finite but unknown correction function. Substituting this decomposition of the field perturbations into the field equations, one can then obtain a set of linearized field equations for the unknown functions, which are highly coupled and complicated partial differential equations. By performing a spectral expansion of the finite correction functions, and evoking the orthogonality of spectral functions, we transform the linearized field equations into a set of linear, homogeneous, algebraic equations of the spectral coefficients of the finite correction functions. The algebraic equations can then be solved as an eigenvalue problem, whose eigenvalues are the QNM frequencies. METRICS can accurately compute the QNM frequencies without decoupling and simplifying the linearized field equations into several master equations. Since METRICS works directly with metric perturbations, one can swiftly reconstruct the metric perturbations by simply reading the eigenvector corresponding to the QNM frequencies, without undergoing the CKK approach.  \nIn this paper, we extend METRICS to a wide class of modified gravity theories. We begin by reviewing the field equations of this class of theories (Sec. II A), background metrics that represent BHs in these theories, and the background scalar-field profile (Sec. II B). Then, we specify the perturbation ansatz of the metric tensor and the scalar field in the Regge-Wheeler gauge (Sec. II C), a gauge that is known to exist in many modified gravity theories. We then derive the asymptotic behavior of the field perturbations at the event horizon and spatial infinity for BHs in modified gravity (Sec. III). The asymptotic behavior at the horizon, derived using the properties of the Killing vector and the horizon, is the first important result of our paper, as it is valid for perturbations that fall through the horizon of a general BH following a null geodesic, regardless of the details of the gravity theory. The asymptotic behavior at the horizon and at spatial infinity allows us to construct the asymptotic factor needed to apply METRICS.  \nWe then continue to apply METRICS to this class of modified gravity theories by product decomposing the field perturbations as a product of the asymptotic factor and a finite correction function. We substitute this form of perturbations into the field equations to derive the linearized field equations for the finite correction functions (Sec. IV A). By expanding the finite correction functions into a product of Chebyshev polynomials along a compactified spatial coordinate and the associated Legendre polynomial along the azimuthal angle, we obtain a set of linear homogeneous algebraic equations of the spectral coefficients (Sec. IV B). Since existing tests of GR indicate that deviations from GR must be small, and since the background metric of BHs in this class of modified gravity theories is constructed only up to leading order in deformations from GR, we develop an eigenvalue per-  \n<!-- image -->  \n<!-- image -->  \nFIG. 1. Trajectory in the complex plane of the quasinormal frequencies of the nlm = 022, 033, and 021 polar (blue) and axial (red) modes in sGB gravity with a dimensionless coupling of ζ = 0 . 1 for a sequence of black holes with various spins. The symbols correspond to the frequencies computed with METRICS, while the lines correspond to polynomial fits to these frequencies. For comparison, the frequencies in GR are also presented (black symbols and gray lines). Observe how the sGB and GR frequency trajectories are close to each other at small spins, but they separate as the spin increases, specially for the axial modes, which couple more strongly with sGB corrections.  \n<!-- image -->  \nturbation theory for METRICS (Sec. V). This eigenvalue perturbation scheme is the second key result of this paper, as it allows us to accurately and consistently compute the leading-order modifications to the QNM frequencies and the metric and scalar-field perturbations.  \nTo exemplify the applications of METRICS to a modified gravity theory, we use METRICS to compute the QNM frequencies of rotating BHs in scalar-Gauss-Bonnet (sGB) gravity with zero potential and a shift-symmetric coupling (Sec. VI). We focus on sGB gravity for the following reasons. First, it is well-motivated as the low-energy limit of certain string theories [60, 61]. Moreover, it has recently been shown that sGB gravity is well-posed [62, 63] for small GR deformations, laying the foundation for the first numerical simulations of BBH coalescence in this theory [64-66]. Computing the QNM frequencies of BHs in sGB gravity is critical to establish a waveform model in this theory, which ultimately leads to an accurate modelspecific test with GW observations. Second, sGB gravity belongs to a wide class of modified gravity theories, known as quadratic gravity theories. The Lagrangian density of quadratic gravity theories consists of quadratic products of the curvature tensor, which can be viewed as the next leading-order correction to the Lagrangian density if one generalizes it into an infinite power series of the curvature tensor [67]. By focusing on sGB gravity, we can learn lessons that will be needed before we apply METRICS to more complicated theories.  \nSpecifically, we apply METRICS to compute the leading-order modifications to the complex frequency of the nlm = 022, 033, and 021 modes of rotating BHs of dimensionless spin a ≲ 0 . 85 in sGB gravity (Sec. VI C). Figure 1 presents a summary of our results through the trajectories of the frequencies in the complex plane for a sequence of black holes with different dimensionless spins. We immediately observe the breaking of isospectrality, i.e. the QNM frequencies of the axial and polar perturbations are different in this theory, as expected [68-76].  \nWe observe also that the sGB frequencies become more and more different from the GR frequencies as the spin increases, specially for the axial modes. For all the QNMs we study, the numerical uncertainty of the leading-order correction to the QNM frequencies is is smaller than 10 -5 for a ≤ 0 . 6, smaller than 10 -4 for 0 . 6 < a ≤ 0 . 7, and smaller than 10 -2 for 0 . 7 < a ≤ 0 . 85 (see also Fig. 4). We construct a polynomial fitting function to the QNM frequencies we calculate numerically (also shown in Fig. 1), which allows for rapid and accurate evaluation. The numerical values and the fitting functions are the third key results of this paper. These fitting functions can now be used to analyze the detected ringdown signals and test sGB gravity. Apart from the QNM frequencies, we also examine the properties of the leading-order modifications to the metric and scalar-field perturbations in sGB gravity, leading to deeper insight into metric perturbations in sGB gravity (Secs. VI E and VI F).  \nThis work is the first to compute the QNM frequencies of BHs that are not slowly rotating and possess scalar hair to the accuracy described above. We note that the QNM frequencies of a rotating BH of a = 0 . 7 in sGB gravity were estimated in [77] using a slow-spin expansion valid only up to second order in dimensionless spin. In what follows, we include sGB corrections to the metric up to 40 orders in dimensionless spin, depending on the dimensionless spin of the BH, to ensure that the effects of the higher spin-order terms in the metric corrections are accurately accounted for when computing the QNM frequencies. Also, in [77], Pad'e resummation is used to help improve the accuracy of the QNM frequencies for a ≥ 0 . 7, whereas here we compute the QNM frequencies at a = 0 . 7 without using any resummation or extrapolation technique. Moreover, QNM frequencies of the axial perturbations are not explicitly computed in [77], whereas, in this work, we also compute the axial-mode frequencies for dimensionless spin a ≤ 0 . 7. Apart from the semianalytical studies described above, the QNM response of a  \nremnant BH formed by binary BH coalescence in sGB gravity were also simulated in [78-80] and beyond an order-reduction scheme [62, 65]. In principle, one could also extract the QNM frequencies of BHs in sGB gravity from these simulations, but the accuracy of this extraction would be limited by numerical errors in the simulations. Hence, our work greatly extends the important work of [62, 65, 77-80], allowing us to obtain the first accurate computation of the QNM frequencies of rapidly spinning BHs in sGB gravity. Our work also extends other important work that focused on QNM frequencies of rapidly rotating BHs using the modified Teukolsky formalism [50-53] in modified gravity theories without BH scalar hair.  \nThe remainder of this paper describes the computational details of the work summarized above. In Sec. II, we review the field equations (Sec. II A) and the BH background metric (Sec. II B) of the wide class of modified gravity theories we consider. We also describe the ansatz of the gravitational and scalar perturbations to BHs in these gravity theories in Sec. II C. In Sec. III, we derive the asymptotic behavior of the gravitational and scalar perturbations of modified BHs at spatial infinity and the event horizon. The highlight of this section is the derivation of the asymptotic behavior of the perturbations at the event horizon through two separate arguments. The asymptotic behavior at the horizon is the first key result of our paper because it is valid for a general BH regardless of the nature of gravity, provided that the BH admits a stationary Killing horizon, and the perturbations follow null geodesics. In Sec. IV, we review the METRICS approach, which we developed in [34, 35], and apply it to transform the linearized field equations in sGB gravity into a set of linear homogeneous algebraic equations. In Sec. V, we develop an eigenvalue perturbation theory to solve the algebraic equations for the leading-order modification to the QNM frequencies due to the activation of modifications to gravity. The eigenvalue perturbation theory is the second important result of this paper in its own right. In Sec. VI, using METRICS and the eigenvalue perturbation theory, we compute the QNM frequencies of rotating BHs in sGB gravity for dimensionless spin up to 0.85. The QNM frequencies of these modes are presented in Table I, II and III of Appendix B. The coefficients (and their uncertainty) of the fitting polynomials of the frequencies are presented in Table V (and VI) of Appendix B.  \nHenceforth, following [34, 35], we assume the following conventions: x µ = ( x 0 , x 1 , x 2 , x 3 ) = ( t, r, χ, ϕ ), where χ = cos θ and θ is the polar angle; the signature of the metric tensor is ( -, + , + , +); gravitational QNMs are labeled by nlm or ( n, l, m ), where n is the principal-mode number, l is the azimuthal-mode number 1 and m is the magnetic-mode number of the QNM; Greek letters in in-  \nsts stand for spacetime coordinates; unless otherwise specified, we used geometric units in which c = G = 1; for the convenience of the reader, we have presented a list of all definitions and symbols in Appendix A.", 'II. MODIFIED GRAVITY THEORIES': 'In this section, we define the class of modified gravity theories we will consider. We begin by describing the field equations and their solution for stationary and axisymmetric black hole spacetimes. We then discuss how we perturb these backgrounds.', 'A. The field equations': "The Lagrangian density of a wide class of modified gravity theories can be written as [20, 81, 82]  \n16 π L = R -1 2 ∇ µ Φ ∇ µ Φ -V (Φ) + αf (Φ) Q , (1)  \nwhere Φ is a scalar field to which the BH in the modified gravity couples, V (Φ) is the potential of Φ, α is a coupling constant, which characterizes the strength of the modifications to gravity and has dimensions of length squared in geometric units, and f (Φ) is a function of Φ only. In the Lagrangian, Q is a scalar constructed from the curvature tensor. For example, in sGB gravity [67, 83], Q is given by  \nQ sGB = G = R 2 -4 R αβ R αβ + R αβγδ R αβγδ , (2)  \nand in dynamical Chern-Simons (dCS) gravity [84-86],  \nQ dCS = P = 1 4 R νµρσ ∗ R µνρσ , (3)  \nwhere ∗ R µνρσ is the dual Riemann tensor  \n∗ R µνρσ = 1 2 ϵ ρσαβ R µν αβ , (4)  \nand ϵ ρσαβ is the Levi-Civita tensor, defined as  \nϵ ρσαβ = 1 √ -g ˜ ϵ ρσαβ , (5)  \ng is the determinant of g µν , and ˜ ϵ ρσαβ is the totally asymmetric Levi-Civita symbol.  \nThrough this work, to exemplify the application of METRICS to study spinning BHs in modified gravity, we focus on the cases of zero potential and a shift-symmetric coupling function,  \nV (Φ) = 0 , f (Φ) = Φ . (6)  \nWe consider a shift-symmetric coupling function because this is the small-coupling approximation (or limit, i.e. when α ≪ 1 such that the background Φ is only slightly displaced from its ground-state value Φ 0 ) of a general coupling function when Q is a topological invariant, as in the sGB and dCS cases. To see the correspondence, we can Taylor expand a general f (Φ) around its ground-state potential Φ 0 ,  \nf (Φ) = f (Φ 0 ) + f ' (Φ 0 )(Φ -Φ 0 ) + ... . (7)  \nThe constant term f (Φ 0 ) -f ' (Φ 0 )Φ 0 has no physical effects when Q is a topological invariant; this is because f (Φ 0 ) -f ' (Φ 0 )Φ 0 pulls out of the integral and Q has zero variation, so the field equations are not modified. The coefficient f ' (Φ 0 ) can just be absorbed into the definition of the coupling constant, and one obtains identical field equations as if one had chosen a shift-symmetric coupling function.  \nUsing the Lagrangian and the least action principle, one can derive the field equations of the modified gravity theory in vacuum, which can be schematically expressed as  \nR µ ν + ζ ( A µ ν -T µ ν ) = 0 , (8)  \n□ ϑ + A ϑ = 0 , (9)  \nwhere ζ is a dimensionless coupling parameter,  \nζ = α 2 M 4 , (10)  \nwhere M is the mass of the background BH we consider in this paper. In the small-coupling approximation, ζ ≪ 1 because of existing tests of GR. The quantities A µ ν and A ϑ are a tensor and a scalar respectively, which represent modifications to GR, while ϑ is a rescaled scalar field, such that Φ = αϑ , and  \nT µ ν ≡ 1 2 ( ∇ µ ϑ ) ( ∇ ν ϑ ) + 1 2 ζ δ µ ν V (Φ) , (11)  \nis the trace-reversed energy-momentum tensor of the rescaled scalar field. In general, A µ ν satisfies the contracted Bianchi identity ∇ ν A µ ν = 0, and may involve higher-than-second-order derivatives of the metric tensor and derivatives of ϑ . For example, in sGB gravity [67, 83],  \nA µ ν ≡ δ νσαβ µλγδ R γδ αβ ∇ λ ∇ σ ϑ -1 2 δ µ ν δ ησαβ ηλγδ R γδ αβ ∇ λ ∇ σ ϑ, (12)  \nwhere δ νσαβ µλγδ is the generalized Kronecker δ , defined as  \nδ ν 1 ν 2 ν 3 ν 4 µ 1 µ 2 µ 3 µ 4 = det     δ ν 1 µ 1 δ ν 1 µ 2 δ ν 1 µ 3 δ ν 1 µ 4 δ ν 2 µ 1 δ ν 2 µ 2 δ ν 2 µ 3 δ ν 2 µ 4 δ ν 3 µ 1 δ ν 3 µ 2 δ ν 3 µ 3 δ ν 3 µ 4 δ ν 4 µ 1 δ ν 4 µ 2 δ ν 4 µ 3 δ ν 4 µ 4     , (13)  \nwhile for dCS gravity [70-73, 84-86],  \nA µν ≡ ( ∇ σ ϑ ) ϵ σδα ( µ | ∇ α R | ν ) δ +( ∇ σ ∇ δ ϑ ) ∗ R δ ( µν ) σ . (14)  \nA ϑ usually involves the product of the curvature tensor. For example, in sGB gravity, A ϑ = G , whereas in dCS gravity, A ϑ = P .", 'B. Background spacetime and scalar field': "Equation (8) allows us to construct a rotating BH spacetime (i.e. a stationary, axisymmetric and vacuum spacetime) as an expansion in powers of the dimensionless spin a in a given modified gravity theory. Since existing constraints indicate that ζ ≪ 1 [20, 82], we here focus on solutions within the small-coupling approximation. We observe that, whereas Eq. (8) explicitly depends on ζ , Eq. (9) does not. Thus, we can solve Eq. (8) in the following iterative manner [67, 87-92]. We first use the GR Kerr metric to compute A ϑ , with which we solve the second equation of Eq. (8) for ϑ . Since it is very difficult to solve the scalar-field equation exactly, one can instead solve for ϑ as a power series of a . Doing so, one finds that the solution takes the form [67, 87, 88]  \nϑ ( r, χ ) = ∑ k =0 N r ( K ) ∑ p =0 N χ ( K ) ∑ q =0 ϑ i,k,p,q a k χ q r p , (15)  \nwhere we recall that χ = cos θ and where ϑ i,k,p,q are constant, while N r ( K ) and N χ ( K ) are also constants that depend on the order in a to which one expands.  \nThen, we use ϑ in the form of Eq. (15) to solve the first equation of Eq. (8) for the background metric. Specifically, we solve the metric in the Boyer-Lindquist coordinates of the following form [50, 67],  \nds 2 = g (0) µν dx µ dx ν = -( 1 -2 Mr Σ -ζH 1 ( r, χ ) ) dt 2 -[1 + ζH 2 ( r, χ )] 4 M 2 ar Σ (1 -χ 2 ) dϕdt +[1 + ζH 3 ( r, χ )] ( Σ ∆ dr 2 + Σ 1 -χ 2 dχ 2 ) +[1 + ζH 4 ( r, χ )] (1 -χ 2 ) × [ r 2 + M 2 a 2 + 2 M 3 a 2 r Σ (1 -χ 2 ) ] dϕ 2 , (16)  \nwhere  \nΣ = r 2 + M 2 a 2 χ 2 , ∆ = ( r -r + )( r -r -) , r ± = M (1 ± b ) , b = √ 1 -a 2 , (17)  \nand where r ± are the GR radial location of the inner and outer event horizons respectively. The parameters  \nM and 0 ≤ a < 1 2 are the observable BH mass, which will be set to unity when we numerically compute the QNM frequencies (see Sec. VI), and the observable dimensionless spin parameter, related to the observable spin angular momentum through J = M 2 a . We denote the Kerr-BH metric in GR by g (GR) µν = g (0) µν ( ζ = 0), and thus, H i =1 ,..., 4 ( r, χ ) are the corrections to the Kerr background due to modified gravity, which can similarly be solved iteratively as series in a [50, 67], with the boundary condition that g (0) µν is asymptotically flat at spatial infinity (i.e. r → + ∞ ). Doing so, one finds that the solution takes the form  \nH i ( r, χ ) = K ∑ k =0 N ' r ( K ) ∑ p =0 N ' χ ( K ) ∑ q =0 h i,k,p,q a k χ q r p , (18)  \nwhere h i,k,p,q are again all numbers, while N ' r ( k ) and N ' χ ( k ) are also numbers that depend on the truncation order of the series in a (that is, on K ). Note that within this parametrization, H i =1 ,..., 4 ( r, χ ) do not explicitly depend on ζ or α . One crucial advantage of this form of metric is that the radial coordinate of the outer and inner event horizons, corresponding to the radial roots of  \ng rr = 1 1 + ζH 3 ( r, χ ) ∆ Σ = 0 , (19)  \nare not changed by the modifications to gravity,  \nr ± = M (1 ± b ) . (20)  \nTwo quantities related to the BH background that we will often use in our work are the horizon angular velocity, Ω H , which can be computed using the tϕ and ϕϕ components of the metric via  \nΩ H = -g tϕ g ϕϕ ∣ ∣ ∣ r = r + , (21)  \nand the surface gravity of the horizon, κ , which can be computed via  \nκ = lim r → r + ∂ r ξ 2 √ ξ 2 /g rr , (22)  \nwhere ξ 2 is the modulus square of the Killing vector of the BH horizon [67], which will be defined in more detail  \nin Sec. III B 1, but for now we define via  \nlim r → r + ξ 2 = [ g tt +2Ω H g tϕ +Ω 2 H g ϕϕ ] r = r + = ( g tt -g 2 tϕ g ϕϕ ) r = r + . (23)  \nBoth can be written as  \nΩ H = Ω (0) H + ζ Ω (1) H , κ = κ (0) + ζκ (1) , (24)  \nwhere Ω (0) H and κ (0) are, respectively, the GR horizon angular velocity and surface gravity  \nΩ (0) H = a 2 Mb , κ (0) = b 2 M (1 + b ) , (25)  \nand Ω (1) H and κ (1) are the leading order in ζ correction to the horizon angular velocity and surface gravity respectively, each of which is a series in a [67]. For the reference of the reader, the explicit power series of Ω (1) H , κ (1) , ϑ and H i =1 ,..., 4 as a power series up to the 40th order of a are given in a Mathematica notebook as Supplemental Material.", 'C. Perturbations of fields': "We now consider both metric and scalar perturbations of a BH in modified gravity theory. The linear metric perturbations can be written as  \ng µν = g (0) µν + ϵ h µν , (26)  \nwhere g (0) µν is the background metric of Eq. (16), h µν is the metric perturbations, and ϵ is a bookkeeping parameter for the perturbations. The metric perturbations h µν can be written as  \nh µν ( t, r, χ, ϕ ) = e imϕ -iωt     h 1 ( r, χ ) h 2 ( r, χ ) -im (1 -χ 2 ) -1 h 5 ( r, χ ) (1 -χ 2 ) ∂ χ h 5 ( r, χ ) ∗ h 3 ( r, χ ) -im (1 -χ 2 ) -1 h 6 ( r, χ ) (1 -χ 2 ) ∂ χ h 6 ( r, χ ) ∗ ∗ ( 1 -χ 2 ) -1 h 4 ( r, χ ) 0 ∗ ∗ ∗ ( 1 -χ 2 ) h 4 ( r, χ ) .     (27)  \nthe modified extremal limit is not strictly relevant for parameter estimation of astrophysical ringdown signals.  \nHere we have made use of the Regge-Wheeler gauge [93, 94], a gauge that many modified gravity theories have sufficient residual degrees of freedom to enforce [70-76]. Apart from the metric perturbations, we also need to consider perturbations of the scalar field to which the BH couples,  \nϑ ( r, χ ) = ϑ (0) ( r, χ ) + ϵe imϕ -iωt h 7 ( r, χ ) , (28)  \nwhere ϑ (0) ( r, χ ) is the unperturbed rescaled scalar field (which is just Eq. (15)).  \nLet us briefly discuss some of the variables introduced in this decomposition of the metric and scalar perturbations. The ω in e -iωt that multiplies h i =1 ,..., 7 in Eqs. (27) and (28) is the same quantity, but this does not mean that we are imposing isospectrality, i.e. that axial, polar, and scalar perturbations have the same frequency. As we shall see in the later sections, we will develop an algorithm to isolate axial modes, the polar modes and the scalar modes. That is, we will develop a method (essentially by picking the right 'parity-isolating' initial guesses) that will ensure that only the h i of a given parity or scalar/metric type are turned on, while the others are suppressed. Doing so ensures that ω we solve for corresponds to the parity or type we have intended to isolate. Thus, Eqs. (27) and (15) are still general enough for the computation of the QNM frequencies when isospectrality is not preserved, even though the frequency of all types of perturbations is labeled by the same ω .  \nBefore proceeding, let us discuss the parity content of Eq. (27). The axial (also known as 'odd' or 'magnetic') metric perturbations are defined to satisfy [95, 96]  \nˆ P [ h (A) µν ] = -( -1) l h (A) µν . (29)  \nHere ˆ P is the parity reversal operator [ ˆ Pf ( χ, ϕ ) = f ( -χ, π + ϕ )] and l is the azimuthal mode number of the QNM. The polar (also known as 'even' and 'electric') metric perturbations are defined to satisfy  \nˆ P [ h (P) µν ] = ( -1) l h (P) µν . (30)  \nThese two conditions do not, in general, imply that h i =1 ,..., 4 are polar and h i =5 , 6 are axial. This statement is only true when a = 0. To illustrate this important observation, let us use the associated Legendre polynomials as an example angular spectral basis for the angular representation of the metric perturbations. These polynomials obey  \nˆ P [ P m ℓ ( χ )] = ( -1) ℓ + m P m ℓ ( χ ) , (31)  \nwhere ℓ and m are the degree and order of the associated Legendre polynomial (and note that the degree of the Legendre polynomial ℓ is not to be confused with the QNM number l ). When a = 0, the perturbations of different ℓ decouple (so that ℓ and l are the same), but when a > 0, the perturbations of different ℓ are coupled.  \nFor spinning BHs then, the metric perturbations h i =1 ,..., 4 and h i =5 , 6 will 'take turns' being of polar and axial type. For example, for the l = 2 modes, the metric perturbations h i =1 ,..., 4 are polar and h i =5 , 6 are axial when ℓ = 2 , 4 , 6 , ... ; however, when ℓ = 3 , 5 , 7 , ... , because of Eq. (31), the metric perturbations h i =1 ,..., 4 are axial and h i =5 , 6 are polar. Therefore, one cannot simply set h i =5 , 6 = 0 to study polar metric perturbations or h i =1 ,..., 4 = 0 to study axial metric perturbations for BH backgrounds that are spinning.", 'III. ASYMPTOTIC BEHAVIOR OF THE PERTURBATIONS OF THE FIELDS': 'Having specified our ansatz of field perturbations around a BH in modified gravity, in this section, we determine the asymptotic behavior of the field perturbations near spatial infinity and the event horizon. In particular, using Boyer-Lindquist coordinates, we derive the asymptotic behavior of a massless field (i.e. the wavefront of the field travels along a null geodesic) near the Killing horizon of a stationary and axisymmetric BH (in or outside GR).', 'A. Behavior near spatial infinity': 'The asymptotic behavior of h i ( r, χ ) can be determined by studying the asymptotic form of the background spacetime near spatial infinity. Near spatial infinity, the background spacetime (Eq. (16)) is asymptotic to  \nds 2 ∼-[ 1 -2 M r ( 1 -1 2 ζH (0) 3 )] dt 2 -4 Ma r ( 1 -1 2 ζH (0) 3 ) (1 -χ 2 ) dtdϕ + [ 1 + ζH (0) 3 ( 1 -M r )] dr 2 1 -2 M r + ( 1 + ζH (0) 3 ) [ r 2 dχ 2 1 -χ 2 + r 2 (1 -χ 2 ) dϕ 2 ] , (32)  \nwhere  \nH (0) 3 = lim r → + ∞ H 3 ( r ) . (33)  \nThis asymptotic form of the background metric allows us to derive the asymptotic limit of the EddingtonFinkelstien coordinate r ∗ . Setting dχ = dϕ = 0, we can write the asymptotic form of the line element near spatial infinity as  \nds 2 ∼ [ 1 -2 M r ( 1 -1 2 ζH (0) 3 )] ( dt 2 -dr 2 ∗ ) , (34)  \nwhere r ∗ has been defined by  \ndr ∗ =   1 + ζH (0) 3 ( 1 -M r ) [ 1 -2 M r ( 1 -1 2 ζH (0) 3 )] ( 1 -2 M r )   1 2 dr. (35)  \nThis choice of r ∗ is such that the null cone has slope of unity near spatial infinity, i.e. dt/dr ∗ = 1 when ds = 0. Expanding the defining equation for r ∗ to first order in ζ and integrating it, we obtain the asymptotic form of r ∗ near spatial infinity as a function of r , namely  \nr ∗ ∼ ( 1 + 1 2 ζH (0) 3 ) r +2 M log ( r -2 M ) . (36)  \nWith this in hand, outgoing null waves near spatial infinity must travel along null cones in ( t, r ∗ ) coordinates, which means that null waves must be proportional to  \ne iωr ∗ ∼ e iM (1+ 1 2 ζH (0) 3 ) ωr r 2 iMω . (37)  \nObserve that metric perturbations do not need to propagate at the speed of light near spatial infinity. This is because the radial Boyer-Lindquist coordinate in modified gravity does not coincide with the usual radial coordinate, R , which measures the surface area of the two-sphere of a constant radius [67]. Instead, the modified gravity Boyer-Lindquist coordinate r and the GR Boyer-Lindquist coordinate R near spatial infinity are related by  \nr ∼ R ( 1 -1 2 ζH (0) 3 ) -1 2 ζH (1) 3 , (38)  \nwhere H (0) 3 and H (1) 3 are defined by the asymptotic behavior of H 3 near spatial infinity,  \nH 3 = H (0) 3 + H (1) 3 r + O ( 1 r 2 ) . (39)  \nThus, the asymptotic behavior of both the metric and scalar-field perturbations near spatial infinity can be written as  \nh k ( r, χ ) = e i ( 1+ 1 2 ζH (0) 3 ) ωr r 2 iMω + ρ ( k ) ∞ ∞ ∑ p =0 a p r p , (40)  \nwhere a p are constants, and ρ ( k ) ∞ is another k -dependent parameter controlling the divergent behavior of h k near spatial infinity (which will be specified later, in Eq. (87)).', 'B. Behavior near the event horizon': 'The asymptotic behavior of h k ( r, χ ) near the event horizon can be determined using the properties of the horizon. Since the derivation of the asymptotic behavior near the horizon is more involved than that near spatial infinity, we will derive it through two different sets of arguments and then cross-check the conclusions to make sure the final result is correct. Anticipating the final result, we will find that the asymptotic behavior of h k near the event horizon, for a fixed χ , must be  \nlim r → r + h k ( r, χ ) ∝ lim r → r + e -iωv + imφ  \n∼ ( r -r + ) -i ω -m Ω H 2 κ -ρ ( k ) H ∞ ∑ p =0 b p ( r -r + ) p , (41)  \nwhere b p are constants and ρ ( k ) H is a k -dependent parameter controlling the divergent behavior of h k at r = r + [35].  \nBefore presenting the details of the derivation of the above result, let us make two remarks. First, having the same asymptotic boundary condition for all metric perturbations seems to contradict the Petrov type of sGB BHs. An sGB BH is of Petrov-type I, meaning that there are four distinct principal null directions, two of which are ingoing and the remaining two are outgoing [88, 97]. However, we notice that the sGB corrections to the Kerr principal null directions found in [97] are finite, and thus, they are suppressed relative to the Kerr principal null directions, which diverge at the event horizon. In other words, the divergent part of the Kerr principal null directions still dominates the asymptotic behavior of massless field perturbations. Second, the asymptotic form of dv and dφ that we will derive below may be useful to perform field quantization near the event horizon, thereby facilitating the study of the Hawking radiation of rotating BHs in modified gravity in the future.', '1. Null geodesic approach': 'Near the event horizon, massless perturbations propagate along null geodesics that are ingoing at the event horizon 3 . These geodesics are more suitably described by ingoing Kerr null coordinates ( v, r, θ, φ ), where  \nv = t + r ∗ , φ = ϕ + ˜ r.  \n(42)  \nHere v is the advanced time, and r ∗ and ˜ r are two coordinates. At the event horizon, these coordinates are the coordinates traced by a congruence generated by a Killing vector at the event horizon (see below). For GR BHs, r ∗ and ˜ r can be determined as a function of r , valid for r ∈ [ r + , + ∞ ), by studying principal null geodesics of the spacetime. But for BHs in modified gravity, r ∗ and ˜ r cannot be determined in this way because some constants of the equations of motions, such as the Carter constant, may not exist, and the (natively second-order) geodesic equations cannot be decoupled into a set of first-order differential equations. Nonetheless, if the asymptotic dependence of v and φ on r as r → r + is known, we can still determine the asymptotic behavior of h k ( r, χ ) as r → r + .  \nOne way to determine the asymptotic form of v as a function of r is using an argument that relies on Killing vectors. First, we note that ˜ ξ µ = t µ + Ω H ϕ µ is a Killing vector in stationary and axisymmetric spacetimes, where t µ and ϕ µ are the Killing vectors associated with stationarity and axisymmetry. Naturally, then, ξ µ = -˜ ξ µ = -( t µ +Ω H ϕ µ ) is also a Killing vector, and geometrically, since ξ µ is null, tangent, and normal to the event horizon, so is -ξ µ . Analytically, the defining equations of a Killing vector, L ξ µ g αβ = 0 (where L ξ µ is the Lie derivative along ξ µ ) or ∇ ( α ξ β ) = 0, are invariant under sign reversal.  \nWith this at hand, we then consider the infinitesimal displacement of a congruence whose four-velocity (defined with respect to advance time v ) is the Killing vector ξ µ , (see, e.g. page 204 of [99]),  \ndx µ = ξ µ dv. (43)  \nTo obtain a differential equation that defines v , we perform the inner product of dx µ with ξ µ ,  \nξ µ dx µ = ξ µ ξ µ dv = ξ 2 dv, (44)  \nwhere ξ 2 = ξ µ ξ µ is the modulus square of the Killing vector at the event horizon, which is [see, e.g. Eq. (4.7) of [67]]  \nξ 2 = g tt g ϕϕ -g 2 tϕ g ϕϕ = O ( r -r + ) , (45)  \nwhere the second equality holds near the event horizon.  \nNext, let us consider the defining equation for the BH  \nsurface gravity,  \nξ µ ∇ µ ξ ν = -κξ ν , (46)  \nwhere κ > 0 is the surface gravity of the BH. Using this equation together with the Killing equation ∇ ( µ ξ ν ) = 0, we then have that [see, e.g., Eq. (5.78) of [99] and Eq. (4.13) of [67]]  \nξ µ = 1 2 κ ∂ µ ξ 2 . (47)  \nNote that this equation, strictly speaking, is valid only on the event horizon because Eq. (46) is only valid at the event horizon [see Sec 4 of [67]]. Note also that Eqs. (46) and (47) differ from those in the literature (e.g., Eq. (5.78) of [99] and Eq. (4.13) of [67]) only by a minus sign, because we have chosen to work with ξ µ instead of ˜ ξ µ and because κ > 0.  \nWith all of this background, we can now proceed to evaluate the asymptotic form of the Killing vector ξ µ . First, we note that since the background metric is a function of χ and r only, so is ξ 2 . Thus, at the horizon, ξ t = ξ ϕ = 0, and explicit calculations show that ξ χ = 0. Hence, near the event horizon, we must have that  \nξ t , ξ χ , ξ ϕ at most ∼ O ( r -r + ) . (48)  \nTechnically, the falloff away from the horizon need not be an integer, but the metric of the background spacetime contains only rational functions of r and χ , so noninteger values will not arise. As for ξ r , we use Eq. (47) to find that near the event horizon  \nξ r = 1 2 κ ∂ r ξ 2 + O ( r -r + ) . (49)  \nNote that the argument that ξ r = 0 implies ξ r = g rα ξ α = 0 does not apply at the event horizon because g rr is singular there.  \nWith all of this preliminary work behind us, we can now finally derive an expression for advanced time as a function of the modified Boyer-Lindquist coordinate r . We begin by substituting the above equation into Eq. (44) and using Eq. (43) to find  \nξ 2 dv = ξ r ( ξ r dv ) = ξ r dr = 1 2 κ ∂ r ξ 2 dr + O ( r -r + ) . (50)  \nRearranging, we have  \ndv = 1 2 κ 1 ξ 2 ∂ r ξ 2 dr + O ( ( r -r + ) 0 ) = 1 2 κ ∂ r log ξ 2 dr + O ( ( r -r + ) 0 ) . (51)  \n̸  \nFrom Eq. (45), we note that ξ 2 has a nonrepeated zero at r = r + , which means we can write ξ 2 = ( r -r + ) × δ , where δ = 0 at r = r + . Hence, asymptotically, at a constantt hypersurface (i.e. dt = 0) and as r → r + , we have  \nv -v 0 ∼ 1 2 κ log( r -r + ) + O ( r -r + ) , (52)  \nNote that the coordinate time t has thus far not appeared in our calculations because we are working on a constanttime hypersurface, so that dt = 0. But even if t were not a constant, its contribution to the asymptotic behavior will be suppressed relative to the logarithm divergence of log( r -r + ). where we have absorbed the δ -dependent term from expanding the logarithm into v 0 , which will then contribute a constant phase factor to h k . Note that while r + has not modified gravity correction (per the choice of coordinate system adopted to find the modified BH metric), the surface gravity κ does admit a ζ -dependent correction, which should be taken into account as we compute the QNM frequencies. Equation (52) indeed gives the correct asymptotic form of the advanced time of the GR Schwarzschild, Reissner-Nordstrom, Kerr and Kerr-Newman BHs [100] when ζ = 0, and it is consistent with the discussion in [101] (see page 28 in that reference).  \nTo find the relation between φ and the modified BoyerLindquist coordinate r , we first recall that, in the GR Kerr background, φ is defined as  \ndφ = dϕ + a ∆ dr. (53)  \nNear the event horizon, dφ is then asymptotic to  \ndφ ∼ dϕ + Ω H 2 κ dr r -r + , (54)  \nwhich is a differential equation relating ϕ and r . We seek a similar equation for a rotating BH in modified gravity, so we go back to Eq. (44), take µ = ϕ , and recall that we are choosing ξ µ = -( t µ +Ω H ϕ µ ) to obtain  \ndϕ = ξ ϕ dv = -Ω H dv ∼ -Ω H 2 κ dr r -r + , (55)  \nwhere we have made use of the asymptotic behavior of dv in Eq. (52). The asymptotic behavior of dϕ suggests that  \ndϕ + Ω H 2 κ dr r -r + ∼ O [ ( r -r + ) p ] , (56)  \nwhere p > -1. In principle, p can assume a noninteger value between -1 and 0. However, as the BH background metric (Eq. (16)) contains only a rational function of r , p can only be a non-negative integer. This implies that dϕ +[Ω H / (2 κ )][ dr/ ( r -r + )] is finite (and thus regular) as r → r + , and one may define a φ coordinate in exactly the same way as done in Eq. (54). When ζ = 0, Eq. (56) reduces to the dφ of the GR Kerr BH (Eq. (54)).  \nUsing the asymptotic form of v and φ , we can determine the asymptotic behavior of h k near the horizon to be  \nh k ( r → r + ) ∝ e -iωv + imφ ∝ ( r -r + ) -i ω -m Ω H 2 κ e imϕ -iωt , (57)  \nwhich is exactly what we anticipated in Eq. (41).', "2. d'Alembertian operator approach": "Another way to determine the asymptotic behavior of h k ( r, χ ) at the event horizon is to carry out an asymptotic analysis of the wave equation of a massless scalar field near the event horizon. Analytically, metric perturbations obey massless wave equations, which are similar to the wave equation governing a massless scalar field. Physically, both the metric and scalar perturbations of a BH follow null geodesics into the BH. Thus, deriving the asymptotic behavior of scalar-field perturbations near the horizon allows us to learn about the leadingorder divergent behavior of the perturbations along the radial coordinate, offering us insight to fully determine the asymptotic behavior of metric perturbations.  \nThe wave equation governing massless scalar-field perturbations is simply  \n✷ Φ = 0 , (58)  \nwhere ✷ is the d'Alembertian operator in curved spacetime,  \n✷ Φ = 1 √ -g ∂ µ ( √ -gg µν ∂ ν Φ ) = 0 , (59)  \nand Φ is a scalar field. Let us rewrite this equation as  \n✷ Φ = ∂ µ ( g µν ∂ ν Φ) + ( ∂ µ log √ -g ) g µν ∂ ν Φ , (60)  \ng  \nand consider the first and second contractions separately.  \nIn Boyer-Lindquist coordinates, the first contraction can be written as  \n∂ µ ( g µν ∂ ν Φ) = g tt ∂ 2 Φ ∂t 2 +2 g tϕ ∂ 2 Φ ∂t∂ϕ + g ϕϕ ∂ 2 Φ ∂ϕ 2 + ∂ ∂r ( g rr ∂ Φ ∂r ) + ∂ ∂χ ( g χχ ∂ Φ ∂χ ) , (61)  \nwhere we recall that the explicit expressions for the components of the inverse of the metric g µν are  \ng g  \ng tt = ϕϕ g tt g ϕϕ -g 2 tϕ , g tϕ = -tϕ g tt g ϕϕ -g 2 tϕ , (62) ϕϕ = g tt g tt g ϕϕ -g 2 tϕ , g rr = 1 g rr , g χχ = 1 g χχ . (63)  \nUsing Eq. (45) evaluated at the event horizon,  \n[ g tt g ϕϕ -g 2 tϕ ] r = r + = g ϕϕ ( r = r + ) ξ 2 , (64)  \nand since ξ 2 ∼ O ( r -r + ) near the event horizon, we then have that  \ng tt ( r = r + ) = [ g 2 tϕ g ϕϕ ] r = r + . (65)  \nWe can now use this expression on the components of the inverse metric to find that  \ng tt ( r = r + ) = 1 ξ 2 , g tϕ ( r = r + ) = -[ g tϕ g ϕϕ ] r = r + 1 ξ 2 = Ω H ξ 2 , g ϕϕ ( r = r + ) = [ g tt g ϕϕ ] r = r + 1 ξ 2 = [ g 2 tϕ g 2 ϕϕ ] r = r + 1 ξ 2 = Ω 2 H ξ 2 , (66)  \nwhere we have assumed that Ω H ≥ 0, since Ω H = -g tϕ ( r = r + ) /g ϕϕ ( r = r + ). Using what we have learned about the behavior of the components of the inverse metric near the event horizon, and multiplying both sides of ✷ Φ = 0 by ξ 2 , we can transform this equation into the following form:  \n∂ 2 Φ ∂t 2 +2Ω H ∂ 2 Φ ∂t∂ϕ +Ω 2 H ∂ 2 Φ ∂ϕ 2 + ξ 2 ∂ ∂r ( g rr ∂ Φ ∂r ) + ξ 2 ∂ ∂χ ( g χχ ∂ Φ ∂χ ) + ξ 2 ( ∂ µ log √ -g ) g µν ∂ ν Φ = 0 . (67)  \nLet us now focus on the second contraction of Eq. (60), namely ( ∂ µ log √ -g ) g µν ∂ ν Φ. Expanding this contraction, we have that  \n( ∂ µ log √ -g ) g µν ∂ ν Φ = ( ∂ r log √ -g ) g rr ∂ r Φ+ ( ∂ χ log √ -g ) g χχ ∂ χ Φ . (68)  \nNow, note that g χχ and ∂ µ log √ -g are at most ∼ O (1) near the event horizon. Thus, ( ∂ χ log √ -g ) g χχ ∂ χ Φ is also at most ∼ O (Φ). As for ( ∂ r log √ -g ) g rr ∂ r Φ, near the event horizon we have that  \ng rr ∼ O ( r -r + ) , ∂ Φ ∂r at most ∼ O [ ( r -r + ) -1 Φ ] . (69)  \nHence ( ∂ r log √ -g ) g rr ∂ r Φ is also at most ∼ O (Φ).  \nWith this in hand, we can now combine what we have learned about the two contractions that define Eq. (60). From the above arguments, we have that  \nξ 2 ∂ ∂χ ( g χχ ∂ Φ ∂χ ) + ξ 2 ( ∂ µ log √ -g ) g µν ∂ ν Φ ∼ O [( r -r + )Φ] . (70)  \nHence, as r → r + , the equation ✷ Φ = 0 is asymptotic to  \n∂ 2 Φ ∂t 2 +2Ω H ∂ 2 Φ ∂t∂ϕ +Ω 2 H ∂ 2 Φ ∂ϕ 2 + ξ 2 ∂ ∂r ( g rr ∂ Φ ∂r ) = 0 . (71)  \nThis equation can be further simplified using the asymptotic limit of g rr at the horizon, which can be computed using Eq. (47) and taking µ = r [67],  \nξ r = 1 2 κ ∂ r ξ 2 . (72)  \nSince ξ t = ξ χ = ξ ϕ = 0 at r = r + , ξ 2 can be computed using just ξ r and g rr at the event horizon,  \nξ 2 ∣ ∣ r + = -g rr ( ξ r ) 2 = -g rr ( ∂ r ξ 2 2 κ ) 2 . (73)  \nWe include the negative sign here because by Eq. (45),  \nlim r -r + → 0 + ξ 2 = 0 -, (74)  \nwhereas g rr > 0 for r > r + . Rearranging, we have the asymptotic value of g rr at the event horizon,  \ng rr | r + = -ξ 2 ( 2 κ ∂ r ξ 2 ) 2 . (75)  \nUsing this asymptotic limit of g rr at the event horizon, we then have that  \nξ 2 ∂ ∂r ( g rr ∂ Φ ∂r ) = -ξ 2 ∂ ∂r [ ξ 2 ( 2 κ ∂ r ξ 2 ) 2 ∂ Φ ∂r ] (76) = -2 κξ 2 ∂ r ξ 2 ∂ ∂r ( 2 κξ 2 ∂ r ξ 2 ∂ Φ ∂r ) -ξ 2 ( ∂ ∂r 2 κ ∂ r ξ 2 ) 2 κξ 2 ∂ r ξ 2 ∂ Φ ∂r , (77)  \nwhere in the second line we used the product rule. As r → r + , ξ 2 ∼ O ( r -r + ), ∂ r Φ is at most O [ ( r -r + ) -1 Φ ] ,  \nand  \nξ 2 ∂ r ξ 2 = 1 ∂ r log ξ 2 = 1 ∂ r log( r -r + ) ∼ O ( r -r + ) , ∂ ∂r 1 ∂ r ξ 2 at most ∼ O (1) , ⇒ ξ 2 ( ∂ ∂r 2 κ ∂ r ξ 2 ) 2 κξ 2 ∂ r ξ 2 ∂ Φ ∂r ∼ O [( r -r + )Φ] . (78)  \nHence, asymptotically, the last term can be discarded, so that ✷ Φ = 0 is asymptotic to,  \n∂ 2 Φ ∂t 2 +2Ω H ∂ 2 Φ ∂t∂ϕ +Ω 2 H ∂ 2 Φ ∂ϕ 2 -2 κξ 2 ∂ r ξ 2 ∂ ∂r ( 2 κξ 2 ∂ r ξ 2 ∂ Φ ∂r ) = 0 . (79)  \nEquation (79) motivates us to define a new radial coordinate, r ∗ , whose asymptotic form as r → r + is  \n∂ ∂r ∗ = 2 κξ 2 ∂ r ξ 2 ∂ ∂r = 2 κ ∂ r log ξ 2 ∂ ∂r . (80)  \nSince ξ 2 ∼ O ( r -r + ), integrating this differential equation yields  \nr ∗ ∼ 1 2 κ log( r -r + ) , (81)  \nas r → r + . Note that this is exactly the same expression we found in Eq. (52) when we considered the derivation in terms of coordinates adapted to ingoing null geodesics. In terms of r ∗ , Eq. (79) becomes  \n∂ 2 Φ ∂t 2 +2Ω H ∂ 2 Φ ∂t∂ϕ +Ω 2 H ∂ 2 Φ ∂ϕ 2 -∂ 2 Φ ∂r 2 ∗ = 0 . (82)  \nIf Φ ∝ e imϕ -iωt , Eq. (82) becomes  \n∂ 2 Φ ∂r 2 ∗ +( ω -m Ω H ) 2 Φ = 0 , (83)  \nwhich implies that the Φ that is ingoing at the horizon is asymptotic to  \nΦ ∝ e -i ( ω -m Ω H ) r ∗ + imϕ -iωt ∼ ( r -r + ) -i ω -m Ω H 2 κ e imϕ -iωt , (84)  \nwhich is the asymptotic behavior we anticipated in Eq. (41).", 'C. Ansatz for the metric and scalar perturbations': "The results above motivate us to 'resum and peel off' the asymptotic behaviors of all of the h k ( r, χ ) perturbations through the following product decomposition  \nh k ( r, χ ) = A k ( r ) u k ( r, χ ) , (85)  \nwhere u k ( r, χ ) is a finite and bounded 'correction factor,' while A k ( r ) is an 'asymptotic factor.' The latter is chosen to be  \nA k ( r ) = e i ( 1+ 1 2 ζH (0) 3 ) ωr r 2 iMω + ρ ( k ) ∞ ( r -r + r ) -i ω -m Ω H 2 κ -ρ ( k ) H (86)  \nso that the correction factor goes to unity both at the event horizon and at spatial infinity. The parameters ρ ( k ) H and ρ ( k ) ∞ in Eq. (86) control the divergence of the metric function at the event horizon and spatial infinity. In principle, these parameters could depend on the index k , the dimensionless spin a , and the coupling constant ζ . For the time being, we shall assume that ρ ( k ) H and ρ ( k ) ∞ are not modified by a and ζ , so that they assume the values for a GR Kerr BH,  \n̸  \nρ ( k ) H =      2 , for k = 3, 1 , for k = 2 or 6, 0 , otherwise, ρ ( k ) ∞ =      2 , for k = 4 nor 7, 1 , for k = 4, -1 , for k = 7. (87)  \nWe determined ρ (7) H and ρ (7) ∞ by inspecting the asymptotic behavior of the solutions to the scalar D'Alembertian equation in a Kerr background [43, 44]. In Sec. VI B, we shall justify the correctness of these ρ ( k ) H and ρ ( k ) ∞ choices through the exponential convergence of the QNM frequencies we calculate.  \nLet us now use that the correction function u k ( r, χ ) is bounded and finite to represent it as a spectral expansion. In order to achieve this, we must change coordinates so that they all have a finite domain. The radial coordinate, in particular, is semi-infinite, as its domain is r ∈ [ r + , ∞ ]. Following [35, 102, 103], let us then define the compactified spatial coordinate  \nz ( r ) = 2 r + r -1 , (88)  \nwhich maps the semi-infinite domain of r to the finite domain z ∈ [ -1 , 1]. As in [34, 35], let us then expand u k ( r, χ ) as a linear combination of the product of Chebyshev and associated Legendre polynomials,  \nu k ( z, χ ) = ∞ ∑ n =0 ∞ ∑ ℓ = | m | v nℓ k T n ( z ) P | m | ℓ ( χ ) , (89)  \nwhere v nℓ k are constant coefficients. As pointed out in [34, 35], we could have chosen a different basis for the spectral expansion, as long as the basis is complete and orthogonal. As we will show below, the product of Chebyshev and associated Legendre polynomials suffices for our purposes.  \n,  \nLet us then compose the ansatz for the metric and field perturbations to obtain  \nh k ( z, χ ) = A k ( z ) ∞ ∑ n =0 ∞ ∑ ℓ = | m | v nℓ k T n ( z ) P | m | ℓ ( χ ) . (90)  \nThis equation provides a formal expression for the full spectral expansion of all field perturbations in the angular coordinate χ and the compactified spatial coordinate z as a double infinite series over the chosen basis functions. In practice, infinite series cannot be used in numerical computations, so we truncate them at some N z for the Chebyshev sum and some N χ for the associated Legendre sum, leading to  \nh k ( r, χ ) = A k ( r ) N z ∑ n =0 N χ + | m | ∑ ℓ = | m | v nℓ k T n [ z ( r )] P | m | ℓ ( χ ) , (91)  \nMoreover, we can also, in principle, choose N z and N χ to be different from each other, but we shall further choose to set N z = N χ = N , so that  \nh k ( r, χ ) = A k ( r ) N ∑ n =0 N + | m | ∑ ℓ = | m | v nℓ k T n [ z ( r )] P | m | ℓ ( χ ) . (92)  \nWhen numerically calculating the QNM frequencies of BHs in modified gravity, we will use this last truncated spectral expansion.", 'IV. LINEARIZED FIELD EQUATIONS ABOUT THE BACKGROUND SPACETIME AS AN ALGEBRAIC PROBLEM IN THE METRICS APPROACH': 'In this section, we will derive the linearized field equations with the field perturbation ansatz described in the previous section. We will then convert these equations into an algebraic problem.', 'A. The linearized field equations': "We begin by substituting the perturbed metric and scalar field of Eq. (90) into the field equations [Eq. (8)]. Linearizing the latter, one finds a system of 11 coupled, partial differential equations for the seven unknown functions h k ( r, χ ). The process of linearization, however, can be simplified slightly as follows. When substituting the asymptotic factor into the field equations for the metric, i.e. R µ ν + ζ ( A µ ν -T µ ν ) = 0, we need to include the explicit ζ modifications to Ω H and κ in the asymptotic factor, because the equation is satisfied up to first order in ζ . However, when considering the scalar-field equation, i.e. □ ϑ + A ϑ = 0, we can set ζ to zero when computing Ω H and κ , because A ϑ is computed in the GR Kerr background. Nonetheless, this is not to say that the scalar-field  \nequation is not modified by ζ , because the d'Alambertian depends on h k ( r, χ ), and thus, the equation is coupled to the equations for the metric perturbations.  \nWe can now cast the linearized field equations into a form that can also be spectrally expanded. We note that, even with the corrections due to modified gravity, the components of the background metric tensor g (0) µν in Boyer-Lindquist coordinates are still rational functions of r and χ [c.f. Eq. (16)]. Therefore, the coefficient functions multiplying the field perturbations h k in the linearized field equations must also be rational functions of r and χ , since they can only depend on background quantities and their derivatives [34, 35]. Thus, we can always cast the k -th linearized field equation, after factorization and multiplication through common denominators, in the following form:  \n7 ∑ j =1 1 ∑ η =0 ∑ α,β =0 ∑ γ =0 d r ∑ δ =0 d χ ∑ σ =0 G k,η,γ,δ,σ,α,β,j ζ η ω γ r δ χ σ ∂ α r ∂ β χ h j = 0 ,  \n(93)  \nwhere G k,η,γ,δ,σ,α,β,j is a complex function of M , m , and a only. The upper limit of the sums over α + β and over γ depends on the nature of the modified gravity theory considered. For sGB gravity, α + β ≤ 3 and γ ≤ 2 because A µ ν involves at most second-order derivatives of the metric and the scalar field [c.f. Eq. (12)]. On the other hand, for dCS gravity, α + β ≤ 4 and γ ≤ 3 because derivatives of the Ricci tensor are involved [c.f. Eq. (14)]. Note that the upper limit of the sum over j has been increased to 7 to accommodate the scalarfield perturbation. The constants d r and d χ are the degree in r and χ of the coefficient of a given term in the linearized field equations respectively, which depend on the component of the linearized field equation we are focusing on, and, thus, on the summation indices ( α, β, k, j ). We sum only from η = 0 to 1, because the metric modifications and the field equations, and thus, all calculations that follow, are valid up to first order in ζ only. After factorizing each of the linearized field equations to obtain common denominators, we find that there are also prefactors, such as powers of 1 -χ 2 , ∆, and Σ, multiplying the equations. These prefactors contain no metric perturbation functions and are nonzero except at r = r + , and χ = ± 1. As these common factors are nonzero in the computational domain (except at the boundaries), we divide the equations by them to simplify them and improve their numerical stability. When ζ = 0, Eq. (93) reduces to the equations that correspond to the metric and scalar perturbations of a Kerr BH in GR.  \nEquation (93) represents a system of coupled and (twodimensional) partial differential equations. As observed in [34, 35], the modulus of G k,η,γ,δ,σ,α,β,j can change by ∼ 10 orders of magnitude across different γ, δ, σ, α, β , and j in one equation, and the largest modulus of G k,η,γ,δ,σ,α,β,j of different equations can also change by ∼ 20 orders of magnitude. As in [34, 35], we prevent overflow by normalizing every partial differential equation, such that  \nthe largest modulus of the coefficient of each equation is one. This is allowed by the homogeneity of the linearized field equations.  \nWe can now substitute the truncated spectral expansion of the field perturbation functions (Eq. (90)) into Eq. (93) to transform the latter into a system of linear algebraic equations. Since r is a rational function of z (see Eq. (88)), the coefficient functions of the linearized equations of u k can only be rational functions of z [34]. Therefore, when we substitute Eq. (85) into Eq. (93), we can factorize the k -th partial differential equation as  \n7 ∑ j =1 1 ∑ η =0 ∑ α,β =0 ∑ γ =0 d z ∑ δ =0 d χ ∑ σ =0 K k,η,γ,δ,σ,α,β,j ζ η ω γ z δ χ σ ∂ α z ∂ β χ u j = 0 ,  \n(94)  \nwhere d z and d χ are the degree of z and χ of the coefficient of the partial derivative ∂ α z ∂ β χ { ... } in the equations respectively, while K k,η,γ,δ,σ,α,β,j are complex functions of M,m , a , ζ , ρ ( k ) H and ρ ( k ) ∞ only (for every value of the summation indices α, β, γ, δ, σ , and j ).", 'B. Converting the linearized field equations into algebraic equations': "Let us now convert the linearized field equations into an algebraic system of equations using our spectral expansion. We first substitute the truncated spectral expansion of the u k functions into Eq. (94),  \n7 ∑ j =1 1 ∑ η =0 ∑ α,β =0 ∑ γ =0 d z ∑ δ =0 d χ ∑ σ =0 K k,η,γ,δ,σ,α,β,j ζ η ω γ z δ χ σ × ∂ α z ∂ β χ { N z ∑ n =0 N χ + | m | ∑ ℓ = | m | v nℓ j T n ( z ) P | m | ℓ ( χ ) } = 0 . (95)  \nThese equations can be further simplified by using the defining equations for the Chebyshev polynomials and associated Legendre polynomials, namely  \nd 2 T n dz 2 = 1 1 -z 2 ( z dT n dz -n 2 T n ) , d 2 P | m | ℓ dχ 2 = 1 1 -χ 2 ( 2 χ dP | m | ℓ dχ -ℓ ( ℓ +1) P | m | ℓ -m 2 1 -χ 2 P | m | ℓ ) . (96)  \nThe above allows us to obtain more factors of 1 -χ 2 , 1 -z , or 1 + z when factorizing the linearized field equations, further simplifying Eq. (94).  \nWe then express the left-hand side of Eq. (94) as a linear combination of the Chebyshev and associated Legendre  \npolynomials,  \nN z ∑ n =0 N χ + | m | ∑ ℓ = | m | w nℓ k T n ( z ) P | m | ℓ ( χ ) = 0 , (97)  \nwhere w nℓ k are independent of z and χ , but depend on M , a , n , ℓ , m , and 4 ω , and k ∈ [1 , 11]. By the orthogonality of T n ( z ) and of P | m | ℓ ( χ ), the linearized field equations are satisfied when w nℓ k = 0 for every k, n and ℓ because they are homogenous. Comparing Eq. (94) to Eq. (95), we notice that w nℓ k depends on v nℓ k linearly, because  \nw nℓ k = 7 ∑ j =1 N z ∑ n ' =0 N χ + | m | ∑ ℓ ' = | m | [ D nℓ,n ' ℓ ' ( ω )] kj v n ' ℓ ' j , (98)  \nwhere D nℓ,n ' ℓ ' ( ω ) are 11 × 7 matrices, whose elements are polynomials in ω and can be obtained by evaluating the inner product given by Eq. (41) of [35].  \nLet us now introduce some notation to simplify the resulting expressions. As in [35], we introduce the following (Euclidean) vectors [35],  \nv nℓ = ( v nℓ 1 , v nℓ 2 , v nℓ 3 , v nℓ 4 , v nℓ 5 , v nℓ 6 , v nℓ 7 ) T , w nℓ = ( w nℓ 1 , w nℓ 2 , ..., w nℓ 11 ) T , (99)  \nso that Eq. (98) can be written as  \nw nℓ = N z ∑ n ' =0 N χ + | m | ∑ ℓ ' = | m | D nℓ,n ' ℓ ' ( ω ) v n ' ℓ ' = 0 , (100)  \nwhere the D nℓ,n ' ℓ ' matrix is now 'dotted' into v n ' ℓ ' (with a flat, Euclidean metric). Let us further define the vectors v and w as  \nv = { v T 00 , v T 01 , ..., v T 0 N χ , ..., v T 1 N χ , ..., v T ( N z +1)( N χ +1) } T , w = { w T 00 , w T 01 , ..., w T 0 N χ , ..., w T 1 N χ , ..., w T ( N z +1)( N χ +1) } T , (101)  \nwhich store all v nℓ and w nℓ respectively. Note that v is a 7( N z +1)( N χ +1)-dimensional vector, whereas w is a 11( N z +1)( N χ +1)-dimensional vector. Then, Eq. (98) can be more compactly written as  \nw = ˜ D ( ω ) v = [ ∑ γ =0 ˜ D γ ω γ ] v = 0 , (102)  \nwhere the ˜ D γ =0 , 1 , 2 ,.... matrices are constant, 11( N z + 1)( N χ + 1) × 7( N z + 1)( N χ + 1) rectangular matrices, which are all linear in ζ . The QNM frequencies of the modified BH in modified gravity correspond to the ω such that Eq. (102) admits a nontrivial solution v .", 'V. FIRST ζ -ORDER EIGENVALUE PERTURBATION IN THE METRICS APPROACH': 'The modified gravity theories we are considering are effective, and thus, their defining actions ought to be understood as curvature (or derivative) expansions. The field equations that result from varying such an approximate action are therefore also approximate. Any solutions of these approximate field equations are, of course, also approximate by construction. More concretely for the problem at hand, the background metric (Eq. (16)) satisfies the field equations only up to the first order in ζ , because the field equations are only valid to this order. For consistency, then, we now compute the first-order-inζ modification to the QNM frequencies. To this end, we develop an eigenvalue-perturbation scheme, which allows us to estimate the first ζ -order modifications to the QNM frequencies.  \nWe begin by expanding D k , ω , and v as a power series in ζ ,  \nD k = ∑ j =0 ζ j D ( j ) k , ω = ∑ j =0 ζ j ω ( j ) , v = ∑ j =0 ζ j v ( j ) , (103)  \nand then truncate the series at j = 1,  \nD k = D (0) k + ζ D (1) k , ω = ω (0) + ζω (1) , v = v (0) + ζ v (1) . (104)  \nSubstituting these perturbed quantities into Eq. (102), we have the zeroth-order-inζ equation,  \n˜ D (0) ( ω (0) ) v (0) = 0 , (105)  \nwhich is the same as the METRICS equation for the QNM frequencies for Kerr black holes in GR [34]. The first-order parts in ζ of Eq. (102) is  \nω (1) ∂ ˜ D (0) ( ω ) ∂ω ∣ ∣ ∣ ω = ω (0) · v (0) + ˜ D (0) ( ω (0) ) · v (1) + ˜ D (1) ( ω (0) ) · v (0) = 0 , (106)  \nwhere we have defined  \n˜ D ( p ) ( ω ( q ) ) = ∑ k =0 ˜ D ( p ) k [ ω ( q ) ] k . (107)  \nRecall that these matrices are power series in ω , and the explicit upper limit of summation depends on the modified theory.  \nTo explain the further steps for solving these equations, let us focus on the polar perturbations as an example. For the polar perturbations, we impose the following conditions to break the homogeneity or linear-scaling invariance of Eq. (106),  \n[ v (0) ] n =0 ,ℓ = | m | k =1 = 1 , [ v (1) ] n =0 ,ℓ = | m | k =1 = 0 . (108)  \n̸  \nThese conditions are allowed by the homogeneity of Eq. (106). If [ v (1) ] n =0 ,ℓ = | m | k =1 = 0, we can always divide the whole v by [ v (0) + ζ v (1) ] n =0 ,ℓ = | m | k =1 to meet the above conditions.  \n̸  \nThe zeroth-order-inζ equation (105) can be solved using the Newton-Raphson algorithms discussed in [34]. After numerically obtaining v (0) and ω (0) , we can use them to solve for v (1) and ω (1) . Since we have set [ v (1) ] n =0 ,ℓ = | m | k =1 = 0, we are left with the following unknowns  \n̸  \n̸  \nx (1) = { [ v (1) ] n =0 ,ℓ = | m | k =1 , ω (1) } . (109)  \nWe can solve Eq.(106) for x (1) by computing the MoorePenrose inverse of the following Jacobian matrix  \n[ J ] ij = ∂ [ D (0) ( ω ) · v (1) ] ∂ [ x (1) ] j ∣ ∣ ∣ ∣ ∣ x = x ( n ) , (110)  \nsuch that  \nx (1) = -[ J ] -1 · [ ˜ D (1) ( ω (0) ) · v (0) ] . (111)  \nThis Jacobian matrix is exactly the same as that defined by Eq. (41) of [34], which was evaluated at the GR x (thus, we here denote them by the same symbol J ). In actual numerical computations, we only need to save the inverse of the Jacobian matrix and x at the last Newton-Raphson step when solving the linearized Einstein equations in GR, and then use them to compute x (1) by Eq. (111). Equation (111) is a key result of this work.  \nThe above procedure can also be similarly performed for the axial perturbations using the following initial guess,  \n[ v (0) ] n =0 ,ℓ = | m | k =5 = 1 , [ v (1) ] n =0 ,ℓ = | m | k =5 = 0 , (112)  \nand for scalar-led perturbations using the following initial guess,  \n[ v (0) ] n =0 ,ℓ = | m | k =7 = 1 , [ v (1) ] n =0 ,ℓ = | m | k =7 = 0 . (113)  \nSince the steps are identical to the polar case, we will not present them here.  \nFrom Eq. (111), we can immediately conclude that isospectrality may not persevere in modified gravity. The v (0) of the axial, polar, and scalar-led perturbations is different. Given a gravity theory, even though ˜ D (1) ( ω ) is the same,, ˜ D (1) ( ω (0) ) · v (0) can be different because v (0) of the axial, polar and scalar perturbations is different. Thus, in general, ω (1) of perturbations led by different sectors is different, which implies the departure from the isospectrality in GR.  \nThe above is not the only way to perturbatively solve for the QNM frequencies through a spectral-expansion method. Recall that the metric perturbation was earlier decomposed into the product of an asymptotic factor A k ( r ) and a correction factor u k ( r, χ ) [see Eq. (85)]. We could have therefore worked directly with the correction factor; in fact, the above procedure is equivalent to writing  \nu j = u (0) j + ζu (1) j , (114)  \nthen first solving  \n7 ∑ j =1 α + β ≤ 3 ∑ α,β =0 2 ∑ γ =0 d z ∑ δ =0 d χ ∑ σ =0 K k,η =0 ,γ,δ,σ,α,β,j ω γ z δ χ σ ∂ α z ∂ β χ u (0) j = 0 ,  \n(115)  \nfor u (0) j , and then using them to solve  \n7 ∑ j =1 α + β ≤ 3 ∑ α,β =0 2 ∑ γ =0 d z ∑ δ =0 d χ ∑ σ =0 K k,η =0 ,γ,δ,σ,α,β,j ω γ z δ χ σ ∂ α z ∂ β χ u (1) j = -7 ∑ j =1 ∑ α,β =0 ∑ γ =0 d z ∑ δ =0 d χ ∑ σ =0 K k,η =1 ,γ,δ,σ,α,β,j ζω γ z δ χ σ ∂ α z ∂ β χ u (0) j (116)  \nfor u (1) j . This, of course, is more complicated than the method we introduced above because would then have to solve a system of coupled partial differential equations for the correction factor. Having said that, the above procedure is not the same as writing  \nh j = h (0) j + ζh (1) j , (117)  \nand then first solving  \n7 ∑ j =1 α + β ≤ 3 ∑ α,β =0 2 ∑ γ =0 d z ∑ δ =0 d χ ∑ σ =0 G k,η =0 ,γ,δ,σ,α,β,j ω γ z δ χ σ ∂ α z ∂ β χ h (0) j = 0 , (118)  \nand then solving  \n7 ∑ j =1 α + β ≤ 3 ∑ α,β =0 2 ∑ γ =0 d z ∑ δ =0 d χ ∑ σ =0 K k,η =0 ,γ,δ,σ,α,β,j ω γ z δ χ σ ∂ α z ∂ β χ h (1) j = -7 ∑ j =1 ∑ α,β =0 ∑ γ =0 d z ∑ δ =0 d χ ∑ σ =0 K k,η =1 ,γ,δ,σ,α,β,j ζω γ z δ χ σ ∂ α z ∂ β χ h (0) j (119)  \nfor h (1) j . In the latter approach, since the asymptotic factor of h (0) j is a GR one, so is that of h (1) j , which means that the effects of modified gravity on the asymptotic factor (or behavior) of h j are ignored. On the other  \nhand, by solving u j perturbatively, the effects of modified gravity on the asymptotic factor are also included.  \nWe conclude this section by drawing a parallel between the eigenvalue perturbation theory of METRICS and the time-independent perturbation theory of quantum mechanics. Equation (111) is the black hole perturbation theory equivalent to the equations of the first-order eigenenergy shift of a quantum mechanical system with a slight time-independent modification to its potential. In particular, ω (0) plays the role of the unperturbed eigenfrequency, v (0) the role of the wave function of an unperturbed eigenstate, D (1) the role of the time-independent modifications to the potential of the quantum mechanical system, [ J ] -1 the role of normalization in quantum mechanics, and the dot products between [ J ] -1 , D (1) and v (0) the volume integral of the modulus square of the unperturbed wave function. In this work, we focus only on the leadingζ order because we work on effective field theories that are only valid to this order. But also as in quantum mechanics, one can, in principle, further develop the eigenvalue perturbation theory in METRICS to an arbitrary order in ζ , provided that the background BH spacetime in the gravity theory is also available up to the corresponding order in ζ .', 'VI. APPLICATION TO ROTATING BHS IN SCALAR-GAUSS-BONNET GRAVITY': 'To illustrate the application of METRICS to BHs in modified gravity, in this section, we apply the eigenvalue perturbation scheme developed above to BHs in sGB gravity, which belongs to a popular class of effective theory that has recently been studied extensively [74-76, 104108]. Specifically, we will apply METRICS to compute the QNM frequencies of the axial and polar perturbations of the 022, 033, and 021 modes. We focus on these l and m modes because there have been claims that suggest the presence of these QNMs in detected ringdown signals 5 [112, 113]. We focus on the fundamental modes ( n = 0) because they have the longest lifetimes. Finally, we focus more extensively on the axial and polar QNMs of the metric perturbation because these are the ones that can be detected with GW interferometers; sGB gravity only possesses the same two tensor polarizations as GR [98]. Nonetheless, for completeness, we do present the scalarmode frequencies of the 022 mode.', 'A. Numerical implementation': 'We first numerically solve Eq. (105) according to the procedure detailed in Sec. IV of [34]. Throughout our nu-  \nmerical computations, M is set to one (i.e. M = 1), thus defining a system of units. To speed up the convergence of our numerical results, we initialize the Newton-Raphson iterations with the known values of the Kerr QNM frequencies. This is a justified procedure, given that we have already showed that METRICS can accurately compute the (known) values of the QNM frequencies of Kerr BHs in GR [34]. The inverse of the Jacobian matrix is computed using the built-in PseudoInverse function of Mathematica to double precision.  \nFormally, H i ( r, χ ), ϑ ( r, χ ), Ω (1) and κ (1) are an infinite series in a . In practice, it is neither possible nor necessary to include H i ( r, χ ) and ϑ ( r, χ ) to infinite order for actual computations that are sought to a desired precision. Instead, one only needs to truncate these quantities at a given order in a , such that the inclusion of higher-order terms do not affect the calculation of our frequencies beyond a desired precision. Through this work, we will truncate H i ( r, χ ), ϑ ( r, χ ), their derivatives, Ω (1) , and κ (1) at the same order in a ,  \nH i ( r, χ ) ≈ N a ∑ k =0 N r ( k ) ∑ p =0 N χ ( k ) ∑ q =0 h i,k,p,q a k χ q r p , ϑ ( r, χ ) ≈ N a ∑ k =0 N r ( k ) ∑ p =0 N χ ( k ) ∑ q =0 ϑ i,k,p,q a k χ q r p , (120)  \nwhere N a is the truncation order in a (yet to be specified). We have checked that the H i ( r, χ ) and ϑ ( r, χ ) used here indeed satisfy the field equations (8) up to the N a th order of a by substituting directly the corrected metric into the field equations.', 'B. Validation Example: the 022 mode at a = 0 . 1': "We start by studying the 022 mode of sGB BHs with a = 0 . 1, whose results already reflect the general QNM features that we will observe for higher a BHs. The blue, red and green symbols in Fig. 2 represent the base-10 logarithm of the backward modulus difference of ω (1) of the axial (left panel) and polar perturbations (right panel),  \nB ( N ) = | ω (1) ( N ) -ω (1) ( N -1) | , (121)  \nas a function of the spectral order N , for N a = 2, 4 and 6, respectively. Observe that, for both parities, B ( N ) first decreases approximately exponentially, until reaching a certain spectral order, at which point B ( N ) fluctuates around an approximate constant.  \nThe saturation of the backward modulus difference is related to the fact that the background metric satisfies the field equations only up to a chosen N a -th order in a . At small N spectral order, the error in the corrected metric is not well resolved, since the error induced by the truncation of the spectral expansion dominates. But as the spectral order is gradually increased, the error due to  \n<!-- image -->  \nN  \nFIG. 2. The backward modulus difference B ( N ) = | ω (1) ( N ) -ω (1) ( N -1) | , where ω (1) is the leading-order modifications to the nlm =022-mode frequency (defined by Eq. (104)) of the axial (left) and polar (right) perturbations of a BH of a = 0 . 1 in shift-symmetric Einstein dilaton Gauss-Bonnet (sGB) gravity as a function of the spectral order ( N ). The background metric of the rotating BH includes modifications due to shift symmetric sGB gravity up to N a = 2, 4, and 6 orders in a . We observe that for all N a , B ( N ) first decreases exponentially and then fluctuates around a constant after reaching a specific spectral order. Nonetheless, the constant about which B ( N ) fluctuates can be significantly reduced if we include the metric corrections of sGB of a larger order in a . This is reasonable because a background metric with sGB corrections of a higher order of a satisfies the field equations better. The ω (1) of the scalar-led perturbations shows a similar tendency, but in the interest of space, we omit the scalar-led ω (1) from the figure.  \n<!-- image -->  \na  \nFIG. 3. The real (left) and imaginary (right) parts of ω (1) of the axial (blue dots) and polar (red triangles) metric perturbations of the nlm = 022, 033 and 021 modes as a function of dimensionless spin. We observe that ω (1) of the axial and polar perturbations does not coincide, which indicates that the isospectrality is broken in sGB gravity. Moreover, we observe that | ω (1) | is smaller for the axial perturbations compared to that of the polar perturbations, and for a ≲ 0 . 3, Re ω (1) and Im ω (1) for the polar perturbations. These two features are consistent with the previous studies regarding slowly rotating BHs in sGB gravity. Finally, we observe that, although we have included the numerical uncertainty (defined by Eq. (122)) of our frequency computations as error bars, the error bars for a ≤ 0 . 8 are too small to be seen in the figure and are visualized in Fig. 4. The small uncertainties indicate that the ω (1) we computed is accurate.  \n<!-- image -->  \nthe N a truncation becomes more important and begins to affect the computation. From Fig. 2, observe that B ( N ) decreases exponentially to a significantly smaller value and begins to fluctuate around a constant value at a larger spectral order as N a is increased.  \nObserve that, at relatively small spectral order (e.g.,  \nfor N ≤ 8 for the axial perturbations and N ≤ 11 for the polar perturbations), the backward modulus difference is approximately independent of N a . This is because, at these small spectral orders, the ω (1) computed with any N a is the same. This feature indicates that increasing N a can help improve the 'convergence' of ω (1) with respect  \n)  \nN  \n(  \nolar  \n-  \n-  \np  \nB  \n10  \nlog  \n-  \n-  \n-  \n0  \n1  \n2  \n3  \n4  \n5  \n5  \n10  \n15  \n20  \n25  \nN  \nN  \nN  \na  \na  \na  \n= 2  \n= 4  \n= 6  \nto the spin expansion.  \nThe results shown in Fig. 2 also offer a guideline to choose N a such that the resulting ω (1) is of the desired accuracy. As estimated in [114], the relative measurement uncertainty of the real and imaginary parts of the QNM frequencies combined across O (10 3 ) ringdown signals detected by the next-generation detectors is O (10 -4 ). Since ζ must be small [i.e. smaller than O (0 . 1)], we should ensure the uncertainty in ω (1) is smaller than O (10 -3 ). Figure 2 suggests that to achieve such an accuracy, we need to select N a such that a N a ≤ 10 -2 at least. This implies that if a = 0 . 1 then N a > 2, but if a = 0 . 3 then N a > 4, and if a = 0 . 7 then N a > 13. Finally, the exponential convergence of ω (1) also justifies the validity of the asymptotic behavior derived in Sec. III C and our choice of ρ ( k ) H and ρ ( k ) ∞ . If any of these details were not correct, ω (1) would not converge exponentially.", 'C. More rapidly rotating BHs': "We apply METRICS to compute the QNM frequencies of more rapidly rotating BHs in sGB gravity. For the purpose of the computations, we have recalculated the sGB metric corrections to 40th order in a , and this is the BH background we employ in our calculations. For a ≤ 0 . 7, we include metric modifications up to N a th order in a , such that N a is the least even integer that satisfies a N a < 10 -4 (we chose an even N a order because the H i ( r, χ ) functions are a power series in a 2 ). We will terminate our computations at a = 0 . 85 because (0 . 85) 40 ∼ 1 . 5 × 10 -3 , smaller than 10 -2 , a necessary criterion that we noted in Sec. VI B. If one wishes to consider BHs that are more rapidly spinning, one would need to use a background metric to higher than 40th order in a small-spin expansion (or an appropriate resummation or numerical background), which is not impossible but is computationally expensive. The QNM frequency for a > 0 . 85 will be extrapolated by constructing a fitting expression for the QNM frequency of a ≤ 0 . 85 (see Eq. (124)). Nevertheless, we expect our QNM frequencies are still useful for the analyses of the detected ringdown signals because the spin of many remnant BHs is < 0 . 85 [7, 11, 115].  \nAt a given a and N a , we compute ω (1) from N = 1 to 25. We stop at 25 because we find that the backward modulus difference is usually saturated or minimized from N = 20 to N = 25. Then, we select the optimal spectral order, defined as  \nN opt = arg min N B ( N ) , (122)  \nwhich allows us to define ω (1) ( N opt ) as the QNM frequency perturbation and  \nδ = B ( N opt ) (123)  \nas the numerical uncertainty in ω (1) .  \nFigure 3 shows the real and imaginary parts of ω (1) of the 022-, 033- and 021-mode frequency of the axial  \nFIG. 4. The base-10 logarithms of the numerical uncertainty, the minimum backward modulus difference of ω (1) for N ≤ 25. We observe that the uncertainty increases with a in general. This is reasonable because the background metric that we constructed only satisfies the field equations to a finite order of dimensionless spin. The error of the background increases with the dimensionless spin even we have adjusted the spintruncation order of the background metric which we keep for QNM frequency computations. Nonetheless, for a ≤ 0 . 7, we can keep the numerical uncertainty of all modes below ∼ 10 -4 . Thus, the ω (1) that we computed using METRICS is accurate enough (see Sec. VI B for the definition of 'desired' accuracy) to be applied to analyze astrophysical ringdown signals.  \n<!-- image -->  \nand polar perturbations, together with their uncertainties (which, for the most part, cannot be seen in the figure). Tables I, II and III of Appendix B in show the numerical value of real and imaginary parts of ω (1) of the 022-, 033-, and 021-mode frequency of both parities. From the figure, we immediately observe that isospectrality is broken, which is expected in modified gravity theories in general [70-76, 116]. Explicitly, we observe that the frequencies of the polar perturbations change more significantly compared to that of the axial perturbations. This is reasonable because the axial perturbations couple more weakly to sGB terms [77, 117, 118] (also see Sec. VI F). We also observe that, for slowly rotating BHs (e.g. a ≲ 0 . 3), the real and imaginary parts of the polar 022-mode frequency both decrease as sGB gravity comes into effect, i.e. Re( ω (1) ) < 0 > Im( ω (1) ) when a ≲ 0 . 3 for the polar 022 mode. Both of these features are consistent with previous studies of the QNM frequencies of slowly rotating BHs in sGB gravity [77, 117, 118] 6 . Apart from the broken isospectrality, observe that the numerical uncertainty (defined by the minimal backward modulus difference across a range of spectral orders) is included in Fig. 3 as error bars. Nonetheless, for a ≤ 0 . 8, the error bars are too small to be visualized while keeping the symbols a reasonable size. This indicates that the numerical uncertainty of ω (1) is tiny.  \nFigure 4 plots the numerical uncertainty as a function  \nof a , which clearly shows that the uncertainty, in general, increases with a . Observe also that the uncertainty of the axial and polar frequencies is different, which is because the ω (1) of the axial and polar modes spans different numerical ranges. The growth of the uncertainty with a is reasonable because the background metric that we constructed satisfies the field equations only to an N a th order in a , which means that the error of our background metric will at least be of ( N a + 1)th order. However, a N a +1 still grows with a even with N a adjusted such that a N a < 10 -4 . For example, at a = 0 . 005, N a = 4 and an order-of-magnitude estimate of the error is a N a +1 ∼ 10 -10 , whereas at a = 0 . 85, N a = 40, and a N a +1 ∼ 10 -3 . These levels of uncertainty in ω (1) are well within the desired accuracy defined in Sec. VI B for the analysis of astrophysical ringdown signals. Moreover, the numerical uncertainty for the relatively large a cases can always be reduced if the metric corrections are kept to a higher order in a .  \nUsing the METRICS frequencies we computed at different spins, we now fit their real and imaginary parts  \nof all modes except the 021 axial mode with a degree-8 polynomial,  \nω (1) = 8 ∑ j =0 w j a j , (124)  \nwhere in this fit we have set M = 1 and the w j are complex constants. We fit these modes as a degree-8 polynomial because we find this degree minimizes the fitting loss functions. As for the frequency of the 021 axial mode, we fit its real and imaginary parts as a degree 14 polynomial because we observed unphysical oscillations for small a if we fit the frequency with a polynomial of a smaller degree. The numerical value of w j is obtained using the builtin NonLinearModelFit function in Mathematica , with numerical uncertainty included as the error of the real and imaginary parts of ω (1) . The explicit numerical value and uncertainty of w j are given in Tables V and VI respectively. For the convenience of the reader, below we provide the fitting polynomials truncated at a 4 :  \nω 022 , A =(0 . 055241 + 0 . 00684399 i ) + (0 . 985294 + 0 . 0688128 i ) a +( -18 . 4902 + 0 . 207235 i ) a 2 +(157 . 802 -3 . 3872 i ) a 3 +( -682 . 224 + 14 . 8323 i ) a 4 + . . . , (125) ω 022 , P =( -0 . 215202 -0 . 0734094 i ) + ( -2 . 51816 -0 . 411031 i ) a +(48 . 4659 + 9 . 59898 i ) a 2 +( -431 . 428 -82 . 4765 i ) a 3 +(1935 . 01 + 378 . 832 i ) a 4 + . . . , (126) ω 033 , A =(0 . 109706 + 0 . 0061152 i ) + (0 . 685561 + 0 . 231546 i ) a +( -11 . 6733 -4 . 20702 i ) a 2 +(108 . 698 + 37 . 6761 i ) a 3 +( -510 . 61 -178 . 723 i ) a 4 + . . . , (127) ω 033 , P =( -0 . 872789 -0 . 113506 i ) + ( -1 . 20198 -0 . 339974 i ) a +(2 . 82484 + 9 . 09087 i ) a 2 +( -35 . 6024 -75 . 7708 i ) a 3 +(165 . 435 + 341 . 913 i ) a 4 + . . . , (128) ω 021 , A =(0 . 0665864 + 0 . 00342921 i ) + ( -1 . 59526 + 0 . 789588 i ) a +(38 . 1384 -16 . 7089 i ) a 2 +( -347 . 625 + 147 . 414 i ) a 3 +(1622 . 8 -663 . 824 i ) a 4 + . . . , (129) ω 021 , P =( -0 . 22612 -0 . 0754879 i ) + (0 . 0933968 + 0 . 0504285 i ) a +( -5 . 91521 -0 . 968127 i ) a 2 +(54 . 1801 + 10 . 6838 i ) a 3 +( -252 . 83 -51 . 1785 i ) a 4 + . . . , (130)  \nTo check that the fitting polynomials obtained can accurately estimate the ω (1) computed using METRICS, we compute ω (1) using the fitting polynomials and compare it with the ω (1) using METRICS. The top panels of Fig. 5 show the real (left) and imaginary (right) parts of ω (1) for the 022, 033 and 021 modes of both parities, computed using the degree-8 fitting polynomials and as a function of a ∈ [0 , 0 . 85]. By visual inspection, the fitting polynomials pass through the ω (1) data computed with METRICS perfectly, indicating that these polynomials can indeed accurately compute ω (1) . To further qualitatively gauge the accuracy of the fitting polynomials, the bottom panels of Fig. 5 show the relative fractional errors of the real (bottom, left) and imaginary (bottom, right) parts of ω (1) ,  \ncomputed using the fitting polynomials and the METRICS data. Observe that the relative fractional error is at most ∼ 10%, consistent with our observation that the fitting polynomials accurately estimate the METRICS ω (1) . The QNM frequencies could also be fitted using other functional forms, such as a rational function in a through a Pad'e function. However, we experimented with this idea and found that fitting the frequencies with a polynomial gives the most accurate estimate of the QNM frequencies at the computed a . This result makes sense since the ω (1) we presented in Fig. 3 clearly show no poles, and thus, a Pad'e resummation is not guaranteed to increase accuracy.  \n<!-- image -->  \nFIG. 5. Real (top left) and imaginary (top right) parts of ω (1) for the 022, 033 and 021 modes of both parities, computed using METRICS (symbols) and using the optimal degree-8 fitting polynomials (curves), as a function of a ∈ [0 , 0 . 85]. The relative fractional errors of the real (bottom left) and imaginary (bottom left) parts of ω (1) , computed using the fitting polynomials and the METRICS data, show the accuracy of the fit. Observe that the optimal fitting polynomials are excellent approximations to the METRICS ω (1) , with an error of at most ∼ 10%.  \n<!-- image -->", 'D. Comparisons with the previous results in slow-rotation expansions': 'Fitting the QNM frequencies as a polynomial in a also allows us to compare our results to those in [77]. In [77], f (Φ) = exp(Φ), but in the small-coupling limit, Φ ≪ 1, we can approximate the coupling function well as f (Φ) ≈ 1 + Φ, which reduces to the coupling function that we considered in this work, as we have explained in the Introduction. Hence, we can compare the w j of our fitting polynomials with the fitting coefficients of the leading-coupling terms in [77]. By reading Eqs. (40) and (47) and Tables I-III of [77], we find that the leading ζ -order modifications to the polar-mode frequencies in [77] are given by 7  \nω (1) , PG 022P = ( -0 . 22496 -0 . 0752 i ) +( -0 . 33536 + 0 . 00064 i ) a +( -0 . 32 + 0 . 30432 i ) a 2 , ω (1) , PG 033P = ( -0 . 87248 -0 . 11504 i ) +( -1 . 03488 -0 . 05664 i ) a (131) +( -0 . 96224 + 0 . 39792 i ) a 2 , ω (1) , PG 021P = ( -0 . 22496 -0 . 0752 i ) +( -0 . 16768 + 0 . 00032 i ) a +(0 . 08176 + 0 . 15408 i ) a 2 .  \nBy comparing the above expressions and Eq. (127), we find that the w j =0 , 1 , 2 of our fitting polynomials are close to the coefficient shown above. We identify several possible causes for the small differences, including but not limited to the following. First, the order in dimensionless spin and ζ that is used in the two studies is different. In this work, at a = 0 . 1, we use N a ≥ 4, whereas in [77] N a = 2 throughout their computations. Moreover, in this work we compute only the leading ζ modifications to the frequencies, whereas terms of higher degree in ζ are also included in the fitting expression of [77]. Second, the two studies used different numerical methods to compute the QNM frequencies. We compute the QNM frequencies using METRICS, whereas in [77] the shooting method is used, which is can be sensitive to the accuracy of the boundary conditions and the shooting point. Third, in this work f (Φ) = Φ, whereas in [77] f (Φ) = exp(Φ). When applying the shooting method, there might be regions where Φ is beyond the small-coupling limit, rendering the comparison not entirely fair. All these discrepancies can lead to differences between the QNM frequencies in the two studies.  \nIn spite of all these possible discrepancies, our fitting polynomials are actually consistent with the results of [77], allowing us therefore to estimate the spin at which the second-order-in-spin approximation becomes highly inaccurate. Let us then compare the frequencies computed using Eq. (131) and that computed using our fitting polynomials by computing the relative fractional differences between them,  \nrelative frac. diff. = ω (1) Re / Im ( a ) -ω (1) , PG Re / Im ( a ) ω (1) Re / Im ( a ) . (132)  \nFIG. 6. The relative fractional difference [defined by Eq. (132)] between fitting polynomials of the real and imaginary parts of the polar frequency of the 022 (solid), 033 (dashed) and 021 (dash-dotted) modes obtained using METRICS and that in [77] through the small-spin expansions.  \n<!-- image -->  \nFigure 6 shows the base-10 logarithms of the absolute value of the relative fractional differences of the real and imaginary parts of the polar frequency of the polar 022 (solid lines), 021 (dashed lines) and 033 (dash-dotted lines) modes. This comparison indicates that our polar frequencies agree well with those computed in [77] using the small-spin expansions for small spins. However, the imaginary part of the dominant mode (022) computed in [77] presents an error larger than 30 % relative to the METRICS frequencies above a = 0 . 3; for other modes [like the (021) mode], the small spin expansion becomes inaccurate for much small spins. As a increases, the frequencies computed using the small-spin expansions become increasingly inaccurate, and the METRICS frequencies should be used instead.', 'E. Parity content of metric perturbations': 'One great advantage of METRICS is its convenience in reconstructing metric perturbations, which requires only the reading and rearrangement of the elements of the eigenvector. This advantage allows us to examine the parity content of the metric perturbations, which can be put to good use to check the sanity of our calculations. Since we compute the modification to the QNM frequency of the axial and polar perturbations separately, we expect that the corresponding v (1) should also be solely axial and polar for consistency.  \nTo quantify the parity content of the metric perturbations, we define the parity dominance (PD) in the following way.  \nPD = Amp(P) Amp(A) , (133)  \nwhere  \n[Amp(P)] 2 = 4 ∑ k =1 N ∑ n =0 ∑ ℓ =even ( ℓ + m )! (2 ℓ +1)( ℓ -m )! | v (1) ,nℓ k ( N ) | 2  \n+ ∑ k =5 , 6 N ∑ n =0 ∑ ℓ =odd ( ℓ + m )! (2 ℓ +1)( ℓ m )! | v (1) ,nℓ k ( N ) | 2 ,  \n2  \n-(134)  \n[Amp(A)] 2 = ∑ k =5 , 6 N ∑ n =0 ∑ ℓ =even ( ℓ + m )! (2 ℓ +1)( ℓ -m )! | v (1) ,nℓ k ( N ) |  \n+ 4 ∑ k =1 N ∑ n =0 ∑ ℓ =even ( ℓ + m )! (2 ℓ +1)( ℓ -m )! | v (1) ,nℓ k ( N ) | 2 . (135)  \nHere ∑ ℓ =even and ∑ ℓ =odd , respectively, stand for the sum of over all ℓ even and odd integers between | m | and N + | m | . Heuristically, the PD gives a rough estimate of the ratio between the amplitude of the polar perturbations to that of the axial ones. Note that unlike the PD which we defined in [34], when estimating the amplitude of the polar perturbations (i.e. the numerator), we do not only take u i =1 ,..., 4 , but also the u i =5 , 6 , in accordance with the definition of the axial and the polar parity that we developed in this paper. Recall from Sec. II C that the polar perturbations satisfy ˆ Ph (P) µν = ( -1) l h (P) µν , where ˆ P is the parity operator. Since we are using associated Legendre polynomials as our angular spectral basis, which obey Eq. (31) upon parity transformation, h i =1 ,..., 4 and h i =5 , 6 will take turns being of polar and axial type. For example, for the l = 2 modes, when ℓ = 2, h i =1 ,..., 4 are polar and h i =5 , 6 are axial; however, when ℓ = 3, because of the parity of the associated Legendre polynomials (see Eq. (31)), h i =1 ,..., 4 are axial and h i =5 , 6 are polar.  \nWith this definition of PD, the PD of purely polar perturbations is ∞ , and that of axial perturbations is 0. However, due to numerical error, the PD of v (1) which we obtain will never be ∞ or 0. Thus, a PD ≫ 1 indicates that the perturbations are dominantly polar, while PD ≪ 1 indicates that the perturbations are dominantly axial. Figure 7 shows the PD of the axial and polar perturbations as a function of a . We see that for all a and the QNMs for which we compute the QNM frequencies, PD ⪆ 10 2 for the polar perturbations and PD ≲ 10 -3 otherwise, despite fluctuations. This indicates that v (1) of the polar perturbations remains mostly polar and that of the axial perturbations remains mostly axial for all a explored.  \nThe astute observer will notice that the parity dominance for some modes does not separate as cleanly as for other modes. For example, at a spin of a = 0 . 2, the 021 mode has a parity dominance of ∼ 10 4 for the polar mode and ∼ 10 -3 for the axial mode; however, the 033 mode has a dominance of ∼ 10 for the polar mode and ∼ 10 -4 for the axial mode. The reason that the polar mode has such a low parity dominance is a technical one. As explained in Sec. V, when solving the linearized field equations for the modified gravity corrections, one must first solve the linearized Einstein equations. In this paper, when doing so, we did not optimize the solution to the linearized Einstein equations to maximize its accuracy; instead, we were content with GR frequencies that are  \nFIG. 7. The base-10 logarithm of the parity dominance (see Eq. (134)) of the modification to the axial and polar perturbations as a function of dimensionless spin ( a ). PD gives a rough estimate of the ratio between the amplitude of the polar perturbations to that of the axial ones. PD ≫ 1 indicates that the perturbations are dominantly polar, and PD ≪ 1 indicates that the perturbations are dominantly axial. We observe that the PD is much larger than 1 (marked by the solid horizontal line in black) for the polar perturbations and much smaller than 1 for the axial perturbations. This feature indicates that the v (1) (defined by Eq. (104)) of the polar perturbations remain dominantly polar; and that of the axial perturbations remain mostly axial.  \n<!-- image -->  \naccurate to better than 10 5 relative to continuous fraction solutions to the Teukolsky equation. Skipping this optimization step (which can be remedied in the future, if desired) means that certain GR modes will be computed more accurately than others, with the 033 one being only accurate to 10 -5 while others are accurate to 10 -8 . This difference in accuracy in the calculation of the GR modes contaminates the accuracy of the sGB corrections, leading to a lower parity dominance for some modes. Having said that, we have ensured that the inaccuracies in the GR modes are not contaminating the sGB corrections to the complex frequencies beyond what was quoted in this paper.', 'F. Understanding the role of the scalar field ϑ': 'In Sec. VI C, we observed that | ω (1) | is larger for the polar perturbations, and we attributed this feature to the fact that axial perturbations couple more weakly to the sGB terms. One way to justify this claim is to examine the amplitude of ϑ , which can be directly read from the elements of v (1) . Inspecting Eqs. (8) and (12), we realize that the linearized field equations contain no terms that involve the product of G and h µν (and their derivatives) in the small and shift-symmetric coupling. Thus, h µν and their derivatives can only interact with G through ϑ , which can be obtained by solving □ ϑ + A ϑ = 0, which nothing but □ ϑ + G = 0 for the shift-symmetric case.  \nFIG. 8. The ϑ abundance (defined by Eq. (136)) of the axial and polar perturbations of different QNMs as a function of a . Heuristically, the ϑ abundance is the ratio between the amplitude of the scalar-field perturbations ( ϑ (1) ) and that of the axial or polar metric perturbations. We observed that the ϑ abundance for the polar perturbations is about 2 orders of magnitude stronger than that for the axial perturbations, indicating that the polar perturbations of different QNMs are always accompanied by a larger scalar field. In other words, the sGB terms exert stronger effects on the polar perturbations than on the axial perturbations, explaining why | ω (1) | of the polar perturbations is larger (see Fig. 3).  \n<!-- image -->  \nMoreover, G (1) is dominated by the polar parity terms 8 . Hence, the amplitude of ϑ indicates the strength of the interplay between the metric perturbations and the sGB terms in the field equations.  \nLet us quantify the strength of ϑ explicitly by defining the scalar-field abundance (i.e. the ϑ abundance) as the ratio between the L 2 norm of ϑ and the amplitude of the axial or polar perturbations. For the modes of even l + | m | , the ϑ abundance for the axial/polar perturbations is then defined as  \nϑ abundance(A/P) = Amp( ϑ ) Amp(A / P) , (136)  \nwhere  \nAmp( ϑ ) = N ∑ n =0 N + | m | ∑ ℓ = | m | ( ℓ + m )! (2 ℓ +1)( ℓ -m )! | v (1) ,nℓ 7 ( N ) | 2 . (137)  \nWe divide Amp( ϑ ) by Amp(A) or Amp(P) before comparison because the amplitude of the axial and polar perturbations may at first be different. Such a division can make the comparison more fair and informative.  \nFigure 8 shows the ϑ abundance of the axial and polar perturbations of different QNMs as a function of a . Observe that the abundance of the polar perturbations is at least about 2 orders of magnitude larger than that  \nof the axial perturbations, which indicates that the polar perturbations are accompanied by a larger ϑ . In other words, the polar perturbations couple to the terms that stem from G more strongly than the axial perturbations do.', 'VII. CONCLUDING REMARKS': 'In this work, we extended METRICS to the study of gravitational perturbations of rotating BHs in modified gravity, where we have focused on sGB gravity as an example. Our work yielded three important results. The first is the asymptotic behavior of gravitational perturbations at the event horizon of a general and axially symmetric BH in modified gravity. The second is the eigenvalue perturbation theory of METRICS, which allows us to compute the modifications to the gravitational QNM frequencies of BHs in modified gravity. The third is the leading-order modifications to the gravitational QNM frequency of the nlm = 022, 033 and 021 modes as a function of dimensionless spins in sGB gravity. The numerical uncertainty of the frequency of these modes is ≲ 10 -5 for a ≤ 0 . 6, ≲ 10 -4 for 0 . 6 < a ≤ 0 . 7, and ≲ 10 -3 for 0 . 7 < a ≤ 0 . 85.  \nOur work demonstrates that METRICS is an effective tool to study the perturbations of BHs in modified gravity. Even without decoupling and simplifying the linearized field equations, METRICS enables us to accurately compute the QNM frequencies in sGB gravity. To the best of our knowledge, our work is the first accurate computation of gravitational QNM frequencies of a rapidly rotating BH (of a ∼ 0 . 85) coupled to a scalar field in modified gravity theories. By simply reading off the eigenvectors of the METRICS solution, we can swiftly reconstruct the metric and scalar perturbations, which could be highly involved in other formalisms [58, 116]. These results offer insight into how different types of perturbations affect the QNM frequencies of BHs in modified gravity theories. The successful application of METRICS to sGB gravity shows that METRICS has the high potential to unlock the QNM spectra and metric perturbations of BHs in other modified gravity theories.  \nMETRICS also fundamentally alters the nature of black hole perturbation studies. As mentioned before, the conventional approach when computing QNM frequencies has traditionally focused on simplifying the linearized field equations into several master equations, which often requires special analytical transformations that have to be devised on a case-by-case basis. METRICS provides a unified framework to compute the QNM frequencies of BHs without these transformations in different modified gravity theories (provided that these theories reduce to GR continuously when the coupling constant vanishes). Our results show that, as long as sufficient computational resources are available, we can always compute the QNM spectra in a modified gravity theory to high accuracy. Although in this work, as a proof of principle, we have computed the QNM frequencies up to a = 0 . 85 in sGB  \ngravity, METRICS can, in principle, be applied to larger a BHs with similar accuracy as long as a sufficiently accurate background metric is available.  \nFuture work could focus on refining the METRICS approach even further. Thus far, METRICS has only been applied within the Regge-Wheeler gauge. Although the Regge-Wheeler gauge is adequate for many modified gravity theories, it will also be beneficial to extend METRICS to other gauges, because other formalisms and theories may need them. For example, the modified Teukolsky formalism of [47, 50, 53, 119] uses the in-going and outgoing radiation gauges. Being able to compute the QNM frequencies and reconstruct the perturbations in other gauges would facilitate cross-checks and comparisons between different formalisms. Also, in this work, we have developed an eigenvalue perturbation theory approach only to leading order in the coupling parameter of the modified theory. In principle, one can further develop this approach to arbitrary order in the coupling parameter. Although for tests of GR, knowing the leading-order modifications to the QNM frequencies should be sufficient, knowing the next-to-leading-order modifications can consolidate our understanding of the effects of modified gravity theories on the QNM response and also further inform us about the stability of BHs in modified gravity.  \nApart from the refinements of METRICS, this framework also opens up several other lines of research. First, the fitting functions for ω (1) (Eq. (124)) that we obtained could be used to analyze detected astrophysical ringdown signals and place the first (sGB) theory-specific constraints with ringdown data; the analyses are ongoing and results will be presented in a future paper. We also plan to extend the studies in this paper to other modified gravity theories, including sGB gravity with different potentials and coupling functions [e.g. V (Φ) = 1 2 µ 2 Φ 2 or f (Φ) = exp( β Φ), for some real µ and β ], and dynamical Chern-Simons gravity; this is also in progress and results will be reported elsewhere. Once we obtain the modifications to the QNM frequencies in these theories, we can construct fitting expressions and place these theories to the test with astrophysical ringdown data.  \nSecond, we can apply METRICS to model the waveforms of extreme-mass-ratio inspirals (EMRIs) around a Kerr BH. Modeling EMRI waveforms usually involves self-force calculations, which require solving the inhomogenous linearized Einstein equations with source terms [120]. In this paper, we have demonstrated that METRICS is efficient in solving inhomogenous linearized field equations. Hence, METRICS could also, in principle, be useful when solving the inhomogenous linearized Einstein equations required in self-force calculations. Moreover, reconstructing the waveform of an EMRI through METRICS only requires reading the eigenvector at the end of the calculation, which spares us from the effort of reconstructing the metric perturbations from the Weyl scalars, the current standard approach in EMRI waveform modeling.  \nFinally, the QNM frequencies presented in this paper  \ncan help establish an accurate waveform model for GWs emitted by binary BHs in sGB gravity, which can ultimately lead to even more tests of sGB gravity that include all coalescence (inspiral, merger and ringdown) stages. Recent breakthroughs in numerical relativity make simulating binary black hole coalescence in sGB gravity possible [62, 66]. However, accurately extracting the QNM frequencies from numerical simulations in sGB gravity requires high resolutions and can be computationally challenging. Our QNM frequencies can spare the need for these high-resolution simulations for accurate QNM frequency extraction.', 'ACKNOWLEDGEMENTS': 'The authors acknowledge the support from the Simons Foundation through Award No. 896696, the NSF through award PHY-2207650 and NASA through Grant No. 80NSSC22K0806. The authors would like to thank Emanuele Berti, Alessandra Buonanno, Vitor Cardoso, Gregorio Carullo, Gregory Gabadadze, Leonardo Gualtieri, Dongjun Li, and Shinji Tsujikawa for insightful discussion about this work. N.Y. would like to thank Takahiro Tanaka for insightful discussion about this work. The calculations and results reported in this paper were produced using the computational resources of the Illinois Campus Cluster, a computing resource that is operated by the Illinois Campus Cluster Program (ICCP) in conjunction with National Center for Supercomputing Applications (NCSA), and is supported by funds from the University of Illinois at Urbana-Champaign, and used Delta at NCSA through allocation PHY240142 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by National Science Foundation grants 2138259, 2138286, 2138307, 2137603, and 2138296. The author would like to specially thank the investors of the IlliniComputes initiatives and GravityTheory computational nodes for permitting the authors to execute runs related to this work using the relevant computational resources.', 'Appendix A: Symbols': "The calculations presented in this paper involved numerous symbols. For the convenience of the reader, we provide a list of the symbols and their definitions in this appendix.  \n- · α , when not suscripted nor superscripted, is the dimensional coupling constant of modify, first defined in Eq. (1).\n- · a is the dimensionless spin of the BH, first defined in Eq. (16).\n- · A k ( r ) is the asymptotic prefactor of the k -th perturbation variable, first defined in Eq. (86).  \n- · (A) is the superscript or supscript which denotes the quantity concerning the axial perturbations, first defined in Table. I.\n- · b = √ 1 -a 2 , first defined in Eq. (17).\n- · B ( N ) is the backward modulus difference of the QNM frequency, first defined in Eq. (121).\n- · d r is the degree of r of the coefficient of the partial derivative of the linearized field equations, first defined in Eq. (93).\n- · d χ is the degree of χ of the coefficient of the partial derivative of the linearized field equations, first defined in Eq. (93).\n- · d z is the degree of z of the coefficient of the partial derivative of the compactified linearized field equations, first defined in Eq. (94).\n- · D ( ω ) is the coefficient matrix of spectral expansion, from one particular basis to another, first defined in Eq. (98).\n- · δ is the numerical uncertainty of the METRICS frequencies, first defined in Eq. (123).\n- · ∆ = ( r -r + )( r -r -), first defined in Eq. (17).\n- · G k,η,γ,δ,σ,α,β,j is the coefficient of ω γ r δ χ σ ∂ α r ∂ β χ h j of the linearized field equations of h j , first defined in Eqs. (93).\n- · h k ( r, χ ) is the functions of metric perturbations, first defined in Eqs. (27) and (28).\n- · H k ( r, χ ) is corrections to the Kerr metric due to modified gravity, first defined in Eqs. (16).\n- · i = √ -1 is the imaginary unit.\n- · k in the subscript is the component of the metric perturbation functions and k = 1 , ..., 7, first defined in Eqs. (27) and (28).\n- · κ is the surface gravity on BH horizon, first defined in Eqs. (24).\n- · K k,α,β,γ,δ,σ,j is the coefficient of ω γ z δ χ σ ∂ α z ∂ β χ ( ... ) of the linearized field equations in z and χ , first defined in Eq. (94).\n- · l is the azimuthal mode number of the gravitational QNMs, first defined in Sec. I.\n- · ℓ is the degree of associate Legendre polynomial used in spectral expansion, first defined in Eq. (89).\n- · M is the BH mass, which is taken to be M = 1 throughout this work, first defined in Eq. (16).\n- · m is the azimuthal number of the metric perturbations, first defined in Eqs. (27) and (28).\n- · N is the number of the Chebyshev and associated Legendre polynomials used in the full spectral expansion, first defined below Eq. (92).\n- · N χ is the number of the associated Legendre polynomials included in the spectral expansion, first defined in Eq. (92).\n- · N z is the number of the Chebyshev polynomials included in the spectral expansion, first defined in Eq. (92).\n- · (P) is the superscript which denotes the quantity concerning the parity-led perturbations, first defined in , first defined in Table. I.\n- · PD is the parity dominance, which characterizes the parity content of metric perturbations, first defined in Eq. (133).\n- · r ± = M (1 ± √ 1 -a 2 ) is the radial coordinate of the position of the event horizon of the Kerr BH, first defined below Eq. (17).\n- · r ∗ is the tortoise coordinate, first defined in Eq. (34).\n- · ρ ( k ) ∞ and ρ ( k ) H are the parameters that characterize the boundary conditions of h k in future null infinity\n- [1] B. P. Abbott et al. (LIGO Scientific, Virgo), 'Observation of Gravitational Waves from a Binary Black Hole Merger,' Phys. Rev. Lett. 116 , 061102 (2016), arXiv:1602.03837 [gr-qc].\n- [2] B. P. Abbott et al. (LIGO Scientific, Virgo), 'GW151226: Observation of Gravitational Waves from a 22-SolarMass Binary Black Hole Coalescence,' Phys. Rev. Lett. 116 , 241103 (2016), arXiv:1606.04855 [gr-qc].\n- [3] B. P. Abbott et al. (LIGO Scientific, VIRGO), 'GW170104: Observation of a 50-Solar-Mass Binary Black Hole Coalescence at Redshift 0.2,' Phys. Rev. Lett. 118 , 221101 (2017), [Erratum: Phys.Rev.Lett. 121, 129901 (2018)], arXiv:1706.01812 [gr-qc].\n- [4] B. . P. . Abbott et al. (LIGO Scientific, Virgo), 'GW170608: Observation of a 19-solar-mass Binary Black Hole Coalescence,' Astrophys. J. Lett. 851 , L35 (2017), arXiv:1711.05578 [astro-ph.HE].\n- [5] 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 (2017), arXiv:1709.09660 [gr-qc].\n- [6] B. P. Abbott et al. (LIGO Scientific, Virgo), 'GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral,' Phys. Rev. Lett. 119 , 161101 (2017), arXiv:1710.05832 [gr-qc].\n- [7] B. P. Abbott et al. (LIGO Scientific, Virgo), 'Tests of general relativity with GW150914,' Phys. Rev. Lett. 116 , 221101 (2016), [Erratum: Phys.Rev.Lett. 121, 129902 (2018)], arXiv:1602.03841 [gr-qc].\n- [8] B. P. Abbott et al. (LIGO Scientific, Virgo), 'Binary Black Hole Mergers in the first Advanced LIGO Ob-  \n- and at the horizon, first defined in Eq. (41) and (40).\n- · Ω H is the angular velocity of BH horizon, first defined in Eq. (24).\n- · ϑ is the scaled scalar field in modified gravity theories, first defined above Eq. (11).\n- · z = 2 r + r -1 is the variable that maps r into a finite domain, first defined in Eq. (88).\n- · ζ = α 2 /M 4 is the dimensionless coupling parameter of modified gravity, first defined in Eq. (10).", 'Appendix B: Additional Tables': "The following tables present additional details related to the results described in the main body of this paper. In particular, Tables I, II, and III present the sGB corrections to the quasinormal frequencies for the 022, 033 and 021 modes respectively, for various choices of BH spin. Meanwhile, Table IV presents the sGB quasinormal frequencies of the scalar mode. Tables V and VI show the coefficients of the fitting polynomials, as well as their uncertainties, respectively.  \n- serving Run,' Phys. Rev. X 6 , 041015 (2016), [Erratum: Phys.Rev.X 8, 039903 (2018)], arXiv:1606.04856 [gr-qc].\n- [9] B. P. Abbott et al. (KAGRA, LIGO Scientific, Virgo, VIRGO), 'Prospects for observing and localizing gravitational-wave transients with Advanced LIGO, Advanced Virgo and KAGRA,' Living Rev. Rel. 21 , 3 (2018), arXiv:1304.0670 [gr-qc].\n- [10] 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 (2019), arXiv:1811.12907 [astro-ph.HE].\n- [11] 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 (2019), arXiv:1903.04467 [gr-qc].\n- [12] R. Abbott et al. (LIGO Scientific, Virgo), 'GW190412: Observation of a Binary-Black-Hole Coalescence with Asymmetric Masses,' Phys. Rev. D 102 , 043015 (2020), arXiv:2004.08342 [astro-ph.HE].\n- [13] R. Abbott et al. (LIGO Scientific, Virgo), 'GW190814: Gravitational Waves from the Coalescence of a 23 Solar Mass Black Hole with a 2.6 Solar Mass Compact Object,' Astrophys. J. Lett. 896 , L44 (2020), arXiv:2006.12611 [astro-ph.HE].\n- [14] C. M. Will, 'The confrontation between general relativity and experiment,' Living Reviews in Relativity 17 , 4 (2014).\n- [15] I. H. Stairs, 'Testing general relativity with pulsar timing,' Living Reviews in Relativity 6 , 5 (2003).  \nTABLE I. ω (1) of the nlm = 022-mode gravitational perturbations of rotating BHs in scalar Gauss Bonnet (sGB) gravity at different dimensionless spins a (first column). The superscripts (A) and (P) respectively stand for axial and polar perturbations. N opt is the optimal spectral order, the spectral order which minimizes the backward modulus difference of ω (1) . δ is the numerical uncertainty of the frequency calculations, which is the respective minimal backward modulus difference of the real and imaginary parts of the frequency. The numerical value of the real and imaginary parts of ω (1) is rounded to the nearest decade which is larger than the numerical uncertainty.  \nTABLE II. Identical to Table I, except that nlm = 033.  \n- [16] N. Wex and M. Kramer, 'Gravity Tests with Radio Pulsars,' Universe 6 , 156 (2020).\n- [17] N. Yunes and X. Siemens, 'Gravitational-Wave Tests of General Relativity with Ground-Based Detectors and Pulsar Timing-Arrays,' Living Rev. Rel. 16 , 9 (2013), arXiv:1304.3473 [gr-qc].\n- [18] K. Yagi and L. C. Stein, 'Black Hole Based Tests of General Relativity,' Class. Quant. Grav. 33 , 054001 (2016), arXiv:1602.02413 [gr-qc].\n- [19] E. Berti, K. Yagi, and N. Yunes, 'Extreme Gravity Tests with Gravitational Waves from Compact Binary Coalescences: (I) Inspiral-Merger,' Gen. Rel. Grav. 50 , 46 (2018), arXiv:1801.03208 [gr-qc].\n- [20] 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].\n- [21] E. Berti, K. Yagi, H. Yang, and N. Yunes, 'Extreme Gravity Tests with Gravitational Waves from Compact Binary Coalescences: (II) Ringdown,' Gen. Rel. Grav. 50 , 49 (2018), arXiv:1801.03587 [gr-qc].\n- [22] R. Abbott et al. (LIGO Scientific, VIRGO, KAGRA), 'Tests of General Relativity with GWTC-3,' (2021), arXiv:2112.06861 [gr-qc].\n- [23] S. Perkins and N. Yunes, 'Are Parametrized Tests of General Relativity with Gravitational Waves Robust to Unknown Higher Post-Newtonian Order Effects?' (2022), arXiv:2201.02542 [gr-qc].\n- [24] A. D. Sakharov, 'Violation of CP Invariance, C asymmetry, and baryon asymmetry of the universe,' Pisma Zh. Eksp. Teor. Fiz. 5 , 32-35 (1967).\n- [25] K. Petraki and R. R. Volkas, 'Review of asymmetric dark matter,' Int. J. Mod. Phys. A 28 , 1330028 (2013),", 'arXiv:1305.4939 [hep-ph].': "- [26] M. Gell-Mann and J. B. Hartle, 'Time symmetry and asymmetry in quantum mechanics and quantum cosmology,' in 4th International Conference on Ion Sources (1991) arXiv:gr-qc/9304023.\n- [27] S. H. S. Alexander, M. E. Peskin, and M. M. SheikhJabbari, 'Leptogenesis from gravity waves in models of inflation,' Phys. Rev. Lett. 96 , 081301 (2006).\n- [28] S. Nojiri and S. D. Odintsov, 'Introduction to modified gravity and gravitational alternative for dark energy,' eConf C0602061 , 06 (2006), arXiv:hep-th/0601213.\n- [29] S. Tsujikawa, 'Modified gravity models of dark energy,' Lect. Notes Phys. 800 , 99-145 (2010), arXiv:1101.0191 [gr-qc].\n- [30] S. Perlmutter et al. (Supernova Cosmology Project), 'Measurements of Ω and Λ from 42 high redshift supernovae,' Astrophys. J. 517 , 565-586 (1999), arXiv:astroph/9812133.\n- [31] A. G. Riess et al. (Supernova Search Team), 'Observational evidence from supernovae for an accelerating universe and a cosmological constant,' Astron. J. 116 , 1009-1038 (1998), arXiv:astro-ph/9805201.\n- [32] Y. Sofue and V. Rubin, 'Rotation curves of spiral galaxies,' Ann. Rev. Astron. Astrophys. 39 , 137-174 (2001), arXiv:astro-ph/0010594.\n- [33] G. Bertone and D. Hooper, 'History of dark matter,' Rev. Mod. Phys. 90 , 045002 (2018), arXiv:1605.04909 [astro-ph.CO].\n- [34] A. K.-W. Chung, P. Wagle, and N. Yunes, 'Spectral method for metric perturbations of black holes: Kerr background case in general relativity,' Phys. Rev. D 109 , 044072 (2024), arXiv:2312.08435 [gr-qc].\n- [35] A. K.-W. Chung, P. Wagle, and N. Yunes, 'Spectral method for the gravitational perturbations of black holes:  \nTABLE III. Identical to Table I, except that nlm = 021.  \nTABLE IV. Identical to Table I, except that the mode is the scalar mode.  \nSchwarzschild background case,' Phys. Rev. D 107 , 124032 (2023), arXiv:2302.11624 [gr-qc].  \n- [36] G. Carullo et al. , 'Empirical tests of the black hole nohair conjecture using gravitational-wave observations,' Phys. Rev. D 98 , 104020 (2018), arXiv:1805.04760 [grqc].\n- [37] G. Carullo, G. Riemenschneider, K. W. Tsang, A. Nagar, and W. Del Pozzo, 'GW150914 peak frequency: a novel consistency test of strong-field General Relativity,' Class. Quant. Grav. 36 , 105009 (2019), arXiv:1811.08744 [grqc].\n- [38] R. Brito, A. Buonanno, and V. Raymond, 'Blackhole Spectroscopy by Making Full Use of GravitationalWave Modeling,' Phys. Rev. D 98 , 084038 (2018), arXiv:1805.00293 [gr-qc].\n- [39] 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].\n- [40] 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].\n- [41] B. P. Abbott et al. (LIGO Scientific, Virgo), 'Tests of general relativity with GW150914,' Phys. Rev. Lett. 116 , 221101 (2016), [Erratum: Phys.Rev.Lett. 121, 129902 (2018)], arXiv:1602.03841 [gr-qc].\n- [42] 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 (2021), arXiv:2010.14529 [gr-qc].\n- [43] S. A. Teukolsky, 'Rotating black holes - separable wave equations for gravitational and electromagnetic perturbations,' Phys. Rev. Lett. 29 , 1114-1118 (1972).\n- [44] S. A. Teukolsky, 'Perturbations of a rotating black hole. 1. Fundamental equations for gravitational electromag-  \nnetic and neutrino field perturbations,' Astrophys. J. 185 , 635-647 (1973).  \n- [45] W. H. Press and S. A. Teukolsky, 'Perturbations of a Rotating Black Hole. II. Dynamical Stability of the Kerr Metric,' Astrophys. J. 185 , 649-674 (1973).\n- [46] S. A. Teukolsky and W. H. Press, 'Perturbations of a rotating black hole. III - Interaction of the hole with gravitational and electromagnet ic radiation,' Astrophys. J. 193 , 443-461 (1974).\n- [47] D. Li, P. Wagle, Y. Chen, and N. Yunes, 'Perturbations of spinning black holes beyond General Relativity: Modified Teukolsky equation,' (2022), arXiv:2206.10652 [gr-qc].\n- [48] A. Hussain and A. Zimmerman, 'An approach to computing spectral shifts for black holes beyond Kerr,' (2022), arXiv:2206.10653 [gr-qc].\n- [49] P. Wagle, D. Li, Y. Chen, and N. Yunes, 'Perturbations of spinning black holes in dynamical Chern-Simons gravity I. Slow rotation equations,' (2023), arXiv:2311.07706 [gr-qc].\n- [50] P. A. Cano, K. Fransen, T. Hertog, and S. Maenaut, 'Quasinormal modes of rotating black holes in higherderivative gravity,' (2023), arXiv:2307.07431 [gr-qc].\n- [51] 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].\n- [52] P. A. Cano, K. Fransen, T. Hertog, and S. Maenaut, 'Gravitational ringing of rotating black holes in higherderivative gravity,' Phys. Rev. D 105 , 024064 (2022), arXiv:2110.11378 [gr-qc].\n- [53] P. A. Cano, K. Fransen, T. Hertog, and S. Maenaut, 'Universal Teukolsky equations and black hole perturbations in higher-derivative gravity,' Phys. Rev. D 108 , 024040 (2023), arXiv:2304.02663 [gr-qc].\n- [54] P. L. Chrzanowski, 'Vector Potential and Metric Perturbations of a Rotating Black Hole,' Phys. Rev. D 11 ,  \nTABLE V. The coefficients w j of the fitting polynomial (c.f. Eq. (124)) to the axial and polar frequencies of the nlm = 022 , 033 and 021.  \nTABLE VI. The uncertainty of the real and imaginary parts of w j of the fitting polynomial (c.f. Eq. (124)) to the axial and polar frequencies of the nlm = 022 , 033 and 021. Note that the uncertainty of the fitting polynomial of the 021 axial mode is undefined because the frequency is actually overfitted to avoid unphysical oscillations of the frequency for small a (see discussion around Eq. (124)).", '2042-2062 (1975).': "- [55] V. Toomani, P. Zimmerman, A. Spiers, S. Hollands, A. Pound, and S. R. Green, 'New metric reconstruction scheme for gravitational self-force calculations,' Class. Quant. Grav. 39 , 015019 (2022), arXiv:2108.04273 [grqc].\n- [56] B. F. Whiting and L. R. Price, 'Metric reconstruction from weyl scalars,' Classical and Quantum Gravity 22 , S589-S604 (2005).\n- [57] C. O. Lousto, 'Reconstruction of black hole metric perturbations from Weyl curvature. II. The Regge-Wheeler gauge,' Class. Quant. Grav. 22 , S569-S588 (2005), arXiv:gr-qc/0501088.\n- [58] J. L. Ripley, N. Loutrel, E. Giorgi, and F. Pretorius, 'Numerical computation of second order vacuum perturbations of Kerr black holes,' Phys. Rev. D 103 , 104018 (2021), arXiv:2010.00162 [gr-qc].\n- [59] J. L. Ripley, 'Computing the quasinormal modes and eigenfunctions for the Teukolsky equation using horizon penetrating, hyperboloidally compactified coordinates,' Class. Quant. Grav. 39 , 145009 (2022), arXiv:2202.03837 [gr-qc].\n- [60] K.-i. Maeda, N. Ohta, and Y. Sasagawa, 'Black Hole Solutions in String Theory with Gauss-Bonnet Curvature Correction,' Phys. Rev. D 80 , 104032 (2009), arXiv:0908.4151 [hep-th].\n- [61] F. Moura and R. Schiappa, 'Higher-derivative corrected black holes: Perturbative stability and absorption crosssection in heterotic string theory,' Class. Quant. Grav.\n- 24 , 361-386 (2007), arXiv:hep-th/0605001 [hep-th].\n- [62] L. Arest'e Sal'o, K. Clough, and P. Figueras, 'WellPosedness of the Four-Derivative Scalar-Tensor Theory of Gravity in Singularity Avoiding Coordinates,' Phys. Rev. Lett. 129 , 261104 (2022), arXiv:2208.14470 [gr-qc].\n- [63] A. D. Kov'acs and H. S. Reall, 'Well-Posed Formulation of Scalar-Tensor Effective Field Theory,' Phys. Rev. Lett. 124 , 221101 (2020), arXiv:2003.04327 [gr-qc].\n- [64] W. E. East and J. L. Ripley, 'Evolution of Einsteinscalar-Gauss-Bonnet gravity using a modified harmonic formulation,' Phys. Rev. D 103 , 044040 (2021), arXiv:2011.03547 [gr-qc].\n- [65] L. Arest'e Sal'o, K. Clough, and P. Figueras, 'Puncture gauge formulation for Einstein-Gauss-Bonnet gravity and four-derivative scalar-tensor theories in d+1 spacetime dimensions,' Phys. Rev. D 108 , 084018 (2023), arXiv:2306.14966 [gr-qc].\n- [66] L. Arest'e Sal'o, S. E. Brady, K. Clough, D. Doneva, T. Evstafyeva, P. Figueras, T. Fran¸ca, L. Rossi, and S. Yao, 'GRFolres: A code for modified gravity simulations in strong gravity,' (2023), arXiv:2309.06225 [gr-qc].\n- [67] P. A. Cano and A. Ruip'erez, 'Leading higher-derivative corrections to kerr geometry,' Journal of High Energy Physics 2019 , 189 (2019), arXiv: 1901.01315.\n- [68] L. Manfredi, J. Mureika, and J. Moffat, 'Quasinormal Modes of Modified Gravity (MOG) Black Holes,' Phys. Lett. B 779 , 492-497 (2018), arXiv:1711.03199 [gr-qc].\n- [69] C.-Y. Chen, M. Bouhmadi-L'opez, and P. Chen, 'Lessons from black hole quasinormal modes in modified gravity,' Eur. Phys. J. Plus 136 , 253 (2021), arXiv:2103.01249 [gr-qc].\n- [70] P. Wagle, N. Yunes, and H. O. Silva, 'Quasinormal modes of slowly-rotating black holes in dynamical ChernSimons gravity,' Phys. Rev. D 105 , 124003 (2022), arXiv:2103.09913 [gr-qc].\n- [71] 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 [grqc].\n- [72] 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].\n- [73] 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].\n- [74] J. L. Bl'azquez-Salcedo, F. S. Khoo, and J. Kunz, 'Quasinormal modes of Einstein-Gauss-Bonnet-dilaton black holes,' Phys. Rev. D 96 , 064008 (2017), arXiv:1706.03262 [gr-qc].\n- [75] 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].\n- [76] L. Pierini and L. Gualtieri, 'Quasi-normal modes of rotating black holes in Einstein-dilaton GaussBonnet gravity: the second order in rotation,' (2022), arXiv:2207.11267 [gr-qc].\n- [77] L. Pierini and L. Gualtieri, 'Quasinormal modes of rotating black holes in Einstein-dilaton Gauss-Bonnet gravity: The second order in rotation,' Phys. Rev. D 106 , 104009 (2022), arXiv:2207.11267 [gr-qc].\n- [78] 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].\n- [79] M. Okounkova, L. C. Stein, J. Moxon, M. A. Scheel, and S. A. Teukolsky, 'Numerical relativity simulation of GW150914 beyond general relativity,' Phys. Rev. D 101 , 104016 (2020), arXiv:1911.02588 [gr-qc].\n- [80] M. Okounkova, 'Numerical relativity simulation of GW150914 in Einstein dilaton Gauss-Bonnet gravity,' Phys. Rev. D 102 , 084046 (2020), arXiv:2001.03571 [grqc].\n- [81] 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].\n- [82] 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].\n- [83] J. L. Ripley and F. Pretorius, 'Gravitational collapse in Einstein dilaton-Gauss-Bonnet gravity,' Class. Quant. Grav. 36 , 134001 (2019), arXiv:1903.07543 [gr-qc].\n- [84] S. Alexander and N. Yunes, 'Chern-Simons Modified General Relativity,' Phys. Rept. 480 , 1-55 (2009), arXiv:0907.2562 [hep-th].  \n- [85] R. Jackiw and S. Y. Pi, 'Chern-Simons modification of general relativity,' Phys. Rev. D 68 , 104012 (2003), arXiv:gr-qc/0308071.\n- [86] S. H. S. Alexander and S. J. Gates, Jr., 'Can the string scale be related to the cosmic baryon asymmetry?' JCAP 06 , 018 (2006), arXiv:hep-th/0409014.\n- [87] N. Yunes and F. Pretorius, 'Dynamical Chern-Simons Modified Gravity. I. Spinning Black Holes in the SlowRotation Approximation,' Phys. Rev. D 79 , 084043 (2009), arXiv:0902.4669 [gr-qc].\n- [88] K. Yagi, N. Yunes, and T. Tanaka, 'Slowly Rotating Black Holes in Dynamical Chern-Simons Gravity: Deformation Quadratic in the Spin,' Phys. Rev. D 86 , 044037 (2012), [Erratum: Phys.Rev.D 89, 049902 (2014)], arXiv:1206.6130 [gr-qc].\n- [89] 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].\n- [90] D. Ayzenberg and N. Yunes, 'Slowly-Rotating Black Holes in Einstein-Dilaton-Gauss-Bonnet Gravity: Quadratic Order in Spin Solutions,' Phys. Rev. D 90 , 044066 (2014), [Erratum: Phys.Rev.D 91, 069905 (2015)], arXiv:1405.2133 [gr-qc].\n- [91] A. Maselli, P. Pani, L. Gualtieri, and V. Ferrari, 'Rotating black holes in Einstein-Dilaton-Gauss-Bonnet gravity with finite coupling,' Phys. Rev. D 92 , 083014 (2015), arXiv:1507.00680 [gr-qc].\n- [92] R. McNees, L. C. Stein, and N. Yunes, 'Extremal black holes in dynamical Chern-Simons gravity,' Class. Quant. Grav. 33 , 235013 (2016), arXiv:1512.05453 [gr-qc].\n- [93] T. Regge and J. A. Wheeler, 'Stability of a schwarzschild singularity,' Phys. Rev. 108 , 1063-1069 (1957).\n- [94] E. Berti, 'Black hole perturbation therory: lectures notes of ICTS Summer School,' (2016).\n- [95] D. Li, A. Hussain, P. Wagle, Y. Chen, N. Yunes, and A. Zimmerman, 'Isospectrality breaking in the Teukolsky formalism,' (2023), arXiv:2310.06033 [gr-qc].\n- [96] D. A. Nichols, A. Zimmerman, Y. Chen, G. Lovelace, K. D. Matthews, R. Owen, F. Zhang, and K. S. Thorne, 'Visualizing Spacetime Curvature via Frame-Drag Vortexes and Tidal Tendexes III. Quasinormal Pulsations of Schwarzschild and Kerr Black Holes,' Phys. Rev. D 86 , 104028 (2012), arXiv:1208.3038 [gr-qc].\n- [97] C. B. Owen, N. Yunes, and H. Witek, 'Petrov type, principal null directions, and Killing tensors of slowly rotating black holes in quadratic gravity,' Phys. Rev. D 103 , 124057 (2021), arXiv:2103.15891 [gr-qc].\n- [98] P. Wagle, A. Saffer, and N. Yunes, 'Polarization modes of gravitational waves in Quadratic Gravity,' Phys. Rev. D 100 , 124007 (2019), arXiv:1910.04800 [gr-qc].\n- [99] E. Poisson, A Relativist's Toolkit: The Mathematics of Black-Hole Mechanics (Cambridge University Press, 2009).\n- [100] B. Carter, 'Global structure of the Kerr family of gravitational fields,' Phys. Rev. 174 , 1559-1571 (1968).\n- [101] A. T. Franzen and P. M. Gir˜ao, 'Regularity of a double null coordinate system for Kerr-Newman-de Sitter spacetimes,' (2020), arXiv:2008.13513 [gr-qc].\n- [102] D. Langlois, K. Noui, and H. Roussille, 'Asymptotics of linear differential systems and application to quasinormal modes of nonrotating black holes,' Phys. Rev. D 104 , 124043 (2021), arXiv:2103.14744 [gr-qc].\n- [103] A. Jansen, 'Overdamped modes in Schwarzschild-de Sitter and a Mathematica package for the numerical computation of quasinormal modes,' Eur. Phys. J. Plus 132 , 546 (2017), arXiv:1709.09178 [gr-qc].\n- [104] 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.\n- [105] K.-i. Maeda, N. Ohta, and Y. Sasagawa, 'Black Hole Solutions in String Theory with Gauss-Bonnet Curvature Correction,' Phys. Rev. D 80 , 104032 (2009), arXiv:0908.4151 [hep-th].\n- [106] T. P. Sotiriou and S.-Y. Zhou, 'Black hole hair in generalized scalar-tensor gravity,' Phys. Rev. Lett. 112 , 251102 (2014), arXiv:1312.3622 [gr-qc].\n- [107] S. Mignemi and N. R. Stewart, 'Charged black holes in effective string theory,' Phys. Rev. D 47 , 5259-5269 (1993), arXiv:hep-th/9212146.\n- [108] S. Mignemi and N. R. Stewart, 'Charged black holes in effective string theory,' Phys. Rev. D 47 , 5259-5269 (1993).\n- [109] 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].\n- [110] M. Isi and W. M. Farr, 'Revisiting the ringdown of GW150914,' (2022), arXiv:2202.02941 [gr-qc].\n- [111] G. Carullo, R. Cotesta, E. Berti, and V. Cardoso, 'Reply to Comment on 'Analysis of Ringdown Overtones in GW150914',' Phys. Rev. Lett. 131 , 169002 (2023), arXiv:2310.20625 [gr-qc].\n- [112] H. Siegel, M. Isi, and W. M. Farr, 'Ringdown of GW190521: Hints of multiple quasinormal modes with a  \n- precessional interpretation,' Phys. Rev. D 108 , 064008 (2023), arXiv:2307.11975 [gr-qc].\n- [113] V. Gennari, G. Carullo, and W. Del Pozzo, 'Searching for ringdown higher modes with a numerical relativityinformed post-merger model,' Eur. Phys. J. C 84 , 233 (2024), arXiv:2312.12515 [gr-qc].\n- [114] A. Maselli, S. Yi, L. Pierini, V. Vellucci, L. Reali, L. Gualtieri, and E. Berti, 'Black hole spectroscopy beyond Kerr: agnostic and theory-based tests with nextgeneration interferometers,' (2023), arXiv:2311.14803 [gr-qc].\n- [115] B. P. Abbott et al. (LIGO Scientific, Virgo), 'Properties of the Binary Black Hole Merger GW150914,' Phys. Rev. Lett. 116 , 241102 (2016), arXiv:1602.03840 [gr-qc].\n- [116] D. Li, A. Hussain, P. Wagle, Y. Chen, N. Yunes, and A. Zimmerman, 'Isospectrality breaking in the Teukolsky formalism,' (2023), arXiv:2310.06033 [gr-qc].\n- [117] J. L. Bl'azquez-Salcedo, C. F. B. Macedo, V. Cardoso, V. Ferrari, L. Gualtieri, F. S. Khoo, J. Kunz, and P. Pani, 'Perturbed black holes in Einsteindilaton-Gauss-Bonnet gravity: Stability, ringdown, and gravitational-wave emission,' Phys. Rev. D 94 , 104024 (2016), arXiv:1609.01286 [gr-qc].\n- [118] 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].\n- [119] A. Hussain and A. Zimmerman, 'An approach to computing spectral shifts for black holes beyond kerr,' (2022).\n- [120] A. Pound and B. Wardell, 'Black hole perturbation theory and gravitational self-force,' (2021), arXiv:2101.04592 [gr-qc]."}