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

Linear perturbations of spherically symmetric spacetimes in general relativity are described by radial wave equations, with potentials that depend on the spin of the perturbing field. In previous work we studied the quasinormal mode spectrum of spacetimes for which the radial potentials are slightly modified from their general relativistic form, writing generic small modifications as a power-series expansion in the radial coordinate. We assumed that the perturbations in the quasinormal frequencies are linear in some perturbative parameter, and that there is no coupling between the perturbation equations. In general, matter fields and modifications to the gravitational field equations lead to coupled wave equations. Here we extend our previous analysis in two important ways: we study second-order corrections in the perturbative parameter, and we address the more complex (and realistic) case of coupled wave equations. We highlight the special nature of coupling-induced corrections when two of the wave equations have degenerate spectra, and we provide a ready-to-use recipe to compute quasinormal modes. We illustrate the power of our parametrization by applying it to various examples, including dynamical Chern-Simons gravity, Horndeski gravity and an effective field theory-inspired model.


I. INTRODUCTION
There are experimental and conceptual reasons to expect that general relativity (GR) and the standard model of particle physics should be modified at some level. Most modifications of GR involve additional gravitational degrees of freedom and higher-order curvature corrections [2,3]. Black holes (BHs) are a promising experimental playground to reveal or constrain these modifications. For example, light bosonic fields can trigger nonperturbative effects in astrophysical BHs via superradiance, affecting the spin distribution of astrophysical BHs and leading to potentially observable gravitational-wave and electromagnetic signatures [4][5][6][7][8][9][10][11].
In the absence of large, smoking-gun effects, we must rely on precision measurements. Astrophysical BHs in GR are remarkably simple, being characterized only by their mass and spin by virtue of the so-called "no-hair theorems" [12][13][14][15][16][17]. As such, they are ideal laboratories for precision measurements: any deviation from this simplicity is a potential hint of new physics. In particular, the relaxation of BH spacetimes to their equilibrium configuration in GR is very simple. Consider for example two BHs merging to form a single spinning BH, a process of particular relevance for gravitational-wave astronomy. The merger can be highly dynamical and violent. However, according to GR, at late times the remnant must * rmcmanu3@jhu.edu † berti@jhu.edu ‡ caiomacedo@ufpa.br § mkimura@rikkyo.ac.jp ¶ andrea.maselli@roma1.infn.it * * vitor.cardoso@ist.utl.pt be a slightly perturbed Kerr solution, described by only two parameters: its mass and spin. The relaxation to a Kerr remnant is well described by linear perturbation theory. This is known as the "ringdown" stage, where the gravitational-wave amplitude consists of a superposition of exponentially damped sinusoids or "quasinormal modes" (QNMs) with characteristic frequencies and damping times [18,19]. A superposition of QNMs accurately describes the merger waveform even before the peak of the gravitational-wave emission [20][21][22][23][24][25][26][27]. Therefore, one possible test of GR consists of testing the consistency between vacuum, linearized GR predictions and the observed QNM spectrum. This idea is now commonly called "black hole spectroscopy" [28][29][30][31].
The general procedure to test a given modification of GR is to find BH solutions, compute their QNM spectrum, and finally constrain deviations of the spectrum from the GR predictions through gravitational-wave observations. The differential equations that must be solved to determine the QNM spectrum have a relatively simple, "universal" structure. A general parametrization of Schwarzschild perturbations induced by scalar, vector and tensor fields shows that linearized perturbations always lead to wave-like equations [32][33][34][35]. However, in general these wave-like equations are coupled. Verifying which theories lead to coupled perturbation equations is a laborious task, but some known examples in the literature include the low-energy limit of string-motivated theories, such as Einstein-dilaton-Gauss-Bonnet [36][37][38][39] and dynamical Chern-Simons (dCS) gravity [40][41][42][43][44]. Coupling also occurs in effective field theory (EFT) modifications of GR [45][46][47][48].
We have recently computed QNM frequencies for scalar, vector and tensor perturbations of a spherically symmetric spacetime which can be described as small deviations from the corresponding GR perturbation equations (see Ref. [1], henceforth Paper I). We wrote deviations in the corresponding radial potentials as a power series in the (inverse) radial coordinate. We found that corrections to the QNM frequencies are (to leading order) linear in these perturbations, and we computed the coefficients that determine these corrections. Here we extend these results by calculating quadratic corrections in the perturbative potentials, as well as the corrections that arise from coupling power-law perturbative corrections between the scalar, vector, polar (even-parity) and axial (odd-parity) gravitational perturbation equations in GR. We still work under the assumption that the background solution is nonspinning (although, as shown in Paper I, the formalism can be applied to spinning black holes in the slow-rotation limit [43,59]) and that the perturbation equations are separable.

A. Executive summary
Our starting point is a generalized, matrix-valued master equation for the coupled radial perturbations induced by N Here f = 1 − r H /r, r H = 2M is the horizon radius, ω is the complex QNM frequency, and V(r) = V i j (r) is a N × N matrix of radial potentials. 1 The factor of f ensures that the effective potential terms vanish at the horizon. We assume the background spacetime to be asymptotically flat, and we parametrize V as a sum of the GR potentials V GR i j and small power-law series corrections δV i j : Here V GR i j denotes the potentials describing massless spin s = 0, 1, 2 perturbations in GR. Usually, V GR i j = 0 for i j and V GR ii 0, since the fields decouple in GR (but see [49,58] for counterexamples). The coefficients α (k) i j are independent of r, but they may be functions of ω [1,42]. For k ≥ 1 the potential matrix vanishes at spatial infinity, i.e. V i j (r) → 0 as r → ∞. Perturbations with k = 0 tend to a constant value, δV i j → α (0) i j /r 2 H . For simplicity, we neglect off-diagonal 1 In principle the coupled perturbation equations may also involve "frictionlike" terms of the form f Z∂ r Φ, i.e. terms of first order in radial derivatives [34]. However the matrix Z can be reabsorbed into V through field redefinitions, as shown in Appendix A.
contributions that fall off slower than the GR potentials [cf.
Eqs. (6)-(8) below]: α (0) i j = α (1) i j = 0 for i j. The terms δV i j in the master equation (1) will, in general, modify the GR QNM frequencies ω 0 by a correction that is perturbatively small in α (k) i j . The key result of this work is that the corrected QNM frequencies read where e i j pq (ks) = e pqi j (sk) , i, j, p, q = 1, . . . , N, k, s = 0, . . . , ∞, and we use the Einstein summation convention. A prime denotes a derivative with respect to ω, with all α (s) pq and α (s) pq evaluated at ω 0 . The derivation of Eq. (4) is presented in Appendix B.
The values of ω 0 for the tensor, vector and scalar perturbations [19,30] and the coefficients d ii (k) [1] are available online [60]. The values of d i j (k) and e i j pq (ks) for these same perturbations were first computed in this work, and they are also available online [60]. We stress that d i j (k) and e i j pq (ks) depend on the unperturbed potentials V GR i j and on the unperturbed QNM frequency ω 0 .
Equation (4) allows for the efficient calculation of the QNM frequencies to quadratic order through simple multiplications and additions. While this expression may look complex, its use is trivial: the prefactors α (s) pq and their derivatives are, in principle, all independent, with their values prescribed by the specific theory in question (cf. Section VI for examples).

B. Plan of the paper
The plan of the paper is as follows. In Section II we present the eigenvalue problem for the QNM frequencies, and we briefly review the numerical method to find them. In Section III we compute quadratic corrections for uncoupled fields. In Section IV we show that coupling fields whose spectra are nondegenerate leads to quadratic corrections in the QNM frequencies. In Section V we show that coupling fields whose spectra are degenerate leads to corrections in the QNM frequencies which are linear in the perturbation parameter. Finally, in Section VI we apply the formalism to some specific examples: dCS gravity [42,44], Horndeski gravity [34], and an EFT-inspired model [46]. In Section VII we discuss some limitations of our analysis and directions for future work.
To improve readability, we relegate several technical results to the appendices. As already mentioned, in Appendix A we demonstrate that friction-like terms (containing first derivatives of the fields) can be reabsorbed in the definition of the potentials, and in Appendix B we derive Eq. (4). In Appendix C we look at the case of three fields. There we show that (i) QNM frequency corrections arising from the coupling of two fields are independent of the total number of coupled fields up to quadratic order, and (ii) when two uncoupled fields have nondegenerate spectra, a linear coupling between the fields gives rise to quadratic QNM frequency corrections. Finally, in Appendix D we prove that when two uncoupled fields have degenerate spectra, a linear coupling between the fields gives rise to linear QNM frequency corrections.

A. Quasinormal modes in general relativity
Gravitational perturbations of the Schwarzschild geometry in GR can be classified by their behavior under parity.2 It is common to classify the metric perturbations as odd (or axial, or Regge-Wheeler) and even (or polar, or Zerilli). These are governed by master variables Φ ± which obey the master equations [62,63] f d dr The effective potential for odd perturbations reads while the effective potential for even perturbations is where λ = 2 + − 2, and is an angular harmonic index labeling the tensorial spherical harmonics used to separate the angular dependence of the perturbations. These two potentials are, quite remarkably, isospectral [64,65], and maintaining isospectrality under generic perturbations of the potentials requires fine tuning [1]. Perturbations of a Schwarzschild background induced by scalar and vector fields are of interest not only in modified gravity theories (that in general introduce additional degrees of freedom), but also in phenomena that involve coupling between gravitational and nongravitational fields. In GR, Schwarzschild perturbations induced by fields φ s of spin s = 0, 1 are also described by master equations similar to Eq. (5), with potentials [19] Note that V 2 = V − , i.e. the s = 2 potential in Eq. (8) corresponds to odd gravitational perturbations.

B. Calculation of the quasinormal frequencies
There are many techniques to compute QNM frequencies [19,43,66]. Since we are striving for generality, here we follow a direct integration approach [43]. The idea is to 2 It is possible to construct definite parity perturbations even in the Kerr background: see e.g. Appendix C of [61].
integrate the radial wave equations from the horizon to infinity given an initial guess of the QNM frequency, and to vary the frequency until Eq. (1) and the relevant boundary conditions are satisfied. The values of the fields at the horizon must be specified to perform the integration. The horizon fields form an N-dimensional vector Φ (i) H and we can perform the integration for any basis of N such Φ (i) H = Φ (i) j , with the final result independent of the choice of basis. For simplicity, in our integrations we consider a basis such that Φ (i) j = δ i j , and we then construct an N × N matrix S from the integration of the N fields under these N initial conditions. The eigenvalues of Eq. (1) are then the complex roots of which can be found numerically.
As we show below, by expanding Eq. (9) with respect to the small parameters α (k) i j we get semianalytic expressions for the coefficients d

III. DECOUPLED FIELDS: QUADRATIC CORRECTIONS
Let us start for simplicity from the case of a single field, and therefore a single nonzero α (k) i j = α in the expansion (3). The QNM frequencies are the roots of Eq. (9), where the matrix S (now a scalar) is a function of both ω and α. As α and α ≡ ∂ ω α are in principle independent, we can vary α while holding α constant. By expanding S up to second order in α we get Let us restrict the expansion of S(ω, α) to those points such that Eq. (9) is satisfied. These points will describe a curve ω(α) starting from ω = ω 0 , α = 0. Replacing the total derivative with respect to α with the directional derivative along ω(α), one finds Each term in the expansion must vanish identically along the curve ω(α). Near α = 0, this curve is approximated by the expansion ω ≈ ω 0 + αd + α 2 e: cf. Eq. (4). By inserting this expansion into the linear term of Eq. (11), we find We can now evaluate this expression at (ω 0 , 0) and solve for d:  Following the same steps for the quadratic term in Eq. (11) we find an expression for e: By construction, both of these expressions depend on the structure of S at α = 0, but not on the value of α in the expansion (4), as long as α is small. Therefore these results encode deviations from the GR spectrum in a theory-independent manner.
We now reinsert the k index labeling specific power-law corrections to the potentials. For linear corrections we get while quadratic corrections yield  (4). The values for e (ks) with k, s = 0, . . . , 10 and ≤ 5 for scalar, vector, axial gravitational and polar gravitational perturbations are available online [60].
Note that Eq. (14) and Eq. (16) agree when k = s, as a result of the definition of e (ks) in the expansion (4).
The linear corrections d (k) where found to five significant figures in Paper I. Here we present numerical results for the coefficients e (ks) for a single field perturbed by a power law potential.
In Table I we list r H e (kk) for axial and scalar gravitational perturbations with = 2 and selected values of k. The values for = 2, . . . , 5 and k, s = 0, . . . , 10 for the scalar, vector, axial gravitational and polar gravitational cases are available online [60]. The large-k behavior of e (kk) for axial perturbations is shown in Fig. 1. In Paper I we found that the linear coefficients d (k) in the large-k limit are well approximated by where (κ, β, γ, ζ) are numerical coefficients. Assuming the same functional form for e (kk) for axial perturbations with = 2, the best-fit parameters are β ≈ 1.7 and γ ≈ 2.4, to be compared with β ≈ 0.66 and γ ≈ 1.5 for d (k) . At quadratic order, cross-term corrections e (ks) also contribute. Representative correction coefficients e (ks) for axial gravitational perturbations with k ≥ s and = 2 are shown in Fig. 2.
These quadratic corrections are necessary when α > 10 −5 . The real and imaginary parts of e (kk) and d (k) are of the same order of magnitude (compare Table I of Paper I with Table I in this paper). If α (kk) e (kk) < 10 −5 d (k) the quadratic correction would be smaller than the numerical error in the leading term, which is currently available to five significant figures. By this argument, we expect the n-th corrections to be needed when α 10 −5/n .

IV. NONDEGENERATE COUPLED FIELDS
Let us now consider the coupling between any two of the scalar, vector, axial gravitational and polar gravitational perturbations (excluding for the moment couplings between axial and polar gravitational perturbations, which will be the subject of Section V below). We will show that d i j (k) is zero for i j (i.e., that leading-order corrections induced by the couplings are quadratic) and that the number of coupled fields does not change the values of d ii (k) = d (k) or e i j pq (ks) . We have produced an extensive list of the coefficients e i j pq (ks) for = 2, . . . , 5 and k, s = 0, . . . , 10, which is available online [60].
The unperturbed QNM spectrum is the union of the spectra for each unperturbed field. In this section we assume that these unperturbed spectra are nondegenerate. Corrections around the tensor QNM spectrum will be called tensor-led in the following. Similarly, corrections around the scalar (vector) QNM spectra are scalar-(vector-)led, respectively. Spectra with multiple branches, corresponding to different fields, have been observed in extreme-mass ratio simulations and nonlinear BH mergers in Einstein-Maxwell theory [67][68][69][70], Chern-Simons theory [71], and scalar Gauss-Bonnet gravity [39,72].
The argument used in the derivation of Eqs. (13) and (16) and e i j pq While these equations look more complex than Eqs. (13) and (16), they are entirely equivalent but expressed in full generality. Moreover, they reduce to the single-field case when We now discuss some important subtleties in problems involving coupled fields.

A. Effect of N coupled fields
One may worry that if we consider the case of three fields and only couple two of them, the resulting d i j (k) may differ from the case of only two fields. This is because the roots of Eq. (9) will remain the same when the number of fields increases, but the functional dependence of the determinant around the roots will not. Therefore, the derivatives in Eqs. (18) and (19) will be different depending on the number of fields we consider.
Fortunately, the structure of Eq. (9) together with the α i j (k) expansion implies that we only need to consider the corrections due to the roots of the 1×1 and 2×2 minors of S which contain the field corresponding to the spectra we perturb about. Hence for each ω 0 , one 1×1 minor and (N −1) 2×2 minors contribute corrections at quadratic order. This is shown in Appendix C. As a result, the value of N plays no role in the calculation of d The 1×1 minor of the matrix S corresponds to the case of an uncoupled field. Moreover, the 1 × 1 minor captures all effects of the additional potentials to quadratic order. Hence, we can make the identification d ii (k) = d (k) and e iiii (ks) = e (ks) , where d (k) and e (ks) are found from the uncoupled case in Section III.
The roots of the 2 × 2 minors will give corrections due to field couplings up to quadratic order. We will now only consider the case of 2 coupled fields and examine the roots of these 2 × 2 minors.

B. Coupling between nondegenerate spectra generates quadratic corrections
Diagonal perturbation terms δV ii (r) generate linear and quadratic corrections to the QNM frequencies. On the contrary, perturbative couplings δV i j (r) only generate quadratic corrections, as long as the spectra of the unperturbed fields are nondegenerate. This is shown in Appendix C by examining the structure of S.
While corrections to the QNM frequencies due to the coupling between two fields with nondegenerate spectra are  (20)] for scalar-axial gravitational and scalar-polar gravitational couplings. We show tensor-led = 2 corrections for a few selected values of k. The coefficients for k, s = 0, . . . , 10 are available online [60]. quadratic in α, the magnitude of these coupling-induced corrections when the values of α are specified need not be smaller than a given linear correction: for example, in dCS gravity α (k) 11 = O( ) and α (k) 12 = O( 1/2 ) (see e.g. [44] and Sec. VI), so the two corrections are formally of the same order in .
Since there are no linear corrections, Eq. (19) simplifies to e i j pq Some representative coefficients for scalar-axial gravitational and scalar-polar gravitational couplings are listed in Table II. We stress again that these results only hold when the two unperturbed fields are not isospectral. For example, in the derivation of Eqs. (18) and (19) we have made the assumption that ∂ ω S| (ω 0 ,0) 0, which does not hold when the union of the unperturbed spectra is degenerate at ω 0 . This assumption rules out the important case of a coupling between the axial and polar gravitational perturbations, which are known to be isospectral. We now turn to the effect of couplings between fields with degenerate spectra.

V. DEGENERATE COUPLED FIELDS
So far we made the assumption that the spectra of the coupled system in Eq. (1) are nondegenerate in the unperturbed case, i.e. when α (k) i j = 0 for all i, j, k. This was used to obtain Eqs. (18) and (19), and also in Appendix C. Unfortunately this assumption is not valid whenever there is a coupling between the axial and polar gravitational perturbations, because in GR the corresponding QNM frequencies are well known to be isospectral [64,65].
In this section we show that couplings between degenerate spectra yield linear corrections to the QNM frequencies. Furthermore, the expansion (4) does not apply in this case. If all elements of δV are nonzero, the total first-order correction to the spectra is not the sum of the corrections found from each individual element of δV, but it is rather a nonlinear combination of these corrections. To illustrate the nature of the problem, let us start again with a simple example. Consider the system defined by for some potential V 0 and coupling Z. The QNM spectra of φ 1 and φ 2 are trivially degenerate for α = 0. To see that corrections enter at linear order, diagonalize the system with the transformation The resulting equations are then It is clear that the corrections to the spectra will enter at linear order in α despite the initial perturbations being coupled, whereas the results of the previous section would imply that the correction should be quadratic in α.
In Appendix D we consider a general potential matrix of the form (3). Assuming a QNM frequency expansion of the form we show that The coefficients δV +− (for example) have the expansion Each of the factors δV (k) ±± and δV (k) ±∓ is related to the expectation value of the kth term in a power-series expansion of the potential perturbations: see Eqs. (D7) and (D8) for their definitions.  (34); bullets were computed through a direct integration method. The inset shows the relative difference between the two calculations for the real (solid black lines) and imaginary (dashed red lines) part of the modes.
One might have hoped that the corrections to the QNM spectra from the coupling of the two degenerate fields would allow for an expansion analogous to Eq. (4). Indeed, in the absence of coupling, the argument of the square root in Eq. (28) becomes a square, and we do recover a linear sum over singlefield expectation values. However, in general, the presence of couplings makes this relation nonlinear. Therefore we provide the values of the quantities δV (k) −+ δV (s) +− for k, s = 0, ...10 and = 2 . . . 5 online [60]. In Table III we also list a small sample of values of δV (k) −+ δV (k) +− for = 2. In order to find linear corrections to the QNM frequencies, these quantities must be plugged into Eqs. (28) and (29). Note that δV (k) ++ and δV (k) −− are just the coefficients d (k) for uncoupled (axial or polar) gravitational perturbations.

VI. EXAMPLES
For illustration, we now apply the formalism to compute QNM spectra for some classes of modified theories of gravity that are known to lead to coupled perturbation equations. Specifically, we consider two models where the coupling is between scalar and tensor modes (dCS gravity [44] and Horndeski gravity [34]) and a model where the coupling is between axial and polar gravitational perturbations (the EFT inspired model [46] not considered in detail in Paper I), so that the background QNM spectra are degenerate.

A. Dynamical Chern-Simons gravity
In dCS gravity, an effective low-energy theory with an additional scalar degree of freedom [41], nonspinning BHs are described by the Schwarzschild metric. The polar sector of gravitational perturbations is the same as in GR, whereas axial gravitational perturbations and scalar perturbations lead to a coupled system of the form (1) [42,71] with the following potentials, in the notation of Eqs. (3), (6) and (8): The parameter β appearing in the dCS action has dimensions [L] −4 and it sets the strength of the coupling, playing a role similar to the Brans-Dicke parameter ω BD . It is useful to introduce a small dimensionless coupling parameterγ ≡ β −1/2 r −2 H such that the equations decouple in the GR limitγ → 0. We first study how the parameterγ modifies the tensor-led mode. Using Eq. (4) and reading off the relevant coefficients from the potentials (30)-(32) we find Proceeding similarly for the scalar-led mode, we find ω = ω 0 + 2d (8) 144π ( + 1)γ 2 + e (88) 144π ( + 1)γ 2 2 + e 1221 (55) 12γ π These expressions illustrate the importance of specifying the form of the coupling: since the off-diagonal terms V 12 and V 21 are proportional to β −1/2 , coupling-induced corrections end up being of the same order as the corrections due to δV 22 . Tensor-led and scalar-led dCS QNM frequencies have been previosuly computed using various methods [71,73]. In Fig. 3 we compare Eqs. (33) and (34) (solid lines) with a numerical QNM calculation based on the direct integration of the perturbed field equations (bullets) for the fundamental = 2 mode. The left panel refers to tensor-led modes, while the right panel refers to scalar-led modes. In the inset we show the relative difference between the quadratic expansions of Eqs. (33) and (34) and the numerical calculation for the real (solid) and imaginary (dashed) parts of the QNM frequencies.

B. Horndeski gravity
The perturbations of scalar and tensor fields on a Schwarzschild background in Horndeski gravity were studied in [34] . The even-parity perturbation equations can be separated through field redefinitions, but there is coupling between even gravitational and scalar perturbations [32].
The master equation for the scalar-led modes takes the form The change to the spectrum is determined by an "effective mass" parameter µ and by a second parameter Γ, built out of the background values (denoted here by overbars) of the free functions appearing in the Horndeski action [74,75]: Let ω 0 belong to the spectrum of V s=0 . Then the scalar-led frequencies can be approximated as As noted in [34], for the class of Horndeski theories in which the GW speed propagation satisfies c T = 1, the Γ factor vanishes identically. In this case Eq. (35) depends on a single free parameter, given by the effective mass µ, and the scalar-led frequencies simply reduce to We compute the values of ω for the = 2 scalar mode in Horndeski gravity using a direct integration method as a function of the effective mass µ. In Fig. 4 we plot numerical results (bullets) against the results obtained from Eq. (39) (solid lines). The two approaches are in excellent agreement, with relative deviations ∆ω = (ω dir − ω fit )/ω fit 10 −4 for both the real and imaginary parts within the range of masses we consider.

C. An effective field theory model
One of the EFT models considered in [46] couples the axial and polar gravitational perturbations. The wave equations map onto Eq. (1) with where V − and V + were defined in Eqs. (6) and (7), V(r) can be found in Appendix A of [46], and the small dimensionless parameter is inversely related to the UV cutoff scale of the EFT. The spectra of V + and V − are degenerate when = 0, so the coupling should induce linear corrections, as discussed in Section V. This expectation is verified in Figure 5. There we consider perturbations about the fundamental QNM frequency with = 2, and we show how the two branches of QNM frequencies -corresponding to the two possible sign choices in Eq. (28) -change with . Note that the coupled QNM frequencies computed by direct integration in Fig. 4 are novel results that were not presented in [46].  Fitting the two branches of QNM frequencies to a seventhorder polynomial in gives a linear coefficient (−0.1479 − 0.2729i)r H for the solid black branch, and (0.1480 + 0.2719i)r H for the dashed red branch. The linear corrections are nonzero, as they should, because the uncoupled axial and polar spectra are degenerate.

VII. CONCLUSIONS AND A COMPUTATIONAL RECIPE
We have extended the formalism of Paper I to compute QNM frequencies of coupled fields as long as the perturbations they induce are small deviations from the perturbation equations for the Schwarzschild geometry in GR. Our main result is a convenient, ready-to-use recipe to compute QNM frequencies at quadratic order in the perturbations. We crucially allow for the possibility of coupling between the master equations. First-order (friction-like) terms in the field derivatives can be accommodated through field redefinitions (cf. Appendix A).
We have found the expansion of the QNM frequencies for uncoupled wave equations to quadratic order. Perhaps our most interesting findings concern the QNM spectra of coupled fields. When the coupling occurs between fields with nondegenerate spectra at zero order in the perturbations, linear-order corrections to the QNM frequencies vanish. However, when the coupling occurs between fields with degenerate spectra, the QNM frequency corrections are linear in the perturbations.
Our results significantly simplify the task of computing QNM frequencies in any modified theory of gravity, or any theory allowing for additional fields. The general recipe for this calculation can be summarized as follows: 1) Derive the master equations for the perturbation variables in the given theory; 2) Eliminate first-order (friction-like) terms in the field derivatives through field redefinitions, as described in Appendix A; 3) Identify the relevant coefficients α (k) i j in the perturbed potentials δV i j [Eq. 4b) If any two unperturbed spectra are degenerate, compute corrections to the QNM frequencies using Eq. (28) and the tabulated values of δV (k) ±± and δV (k) ±∓ defined in Eqs. (D7) and (D8), which are also available online [60].
In Sec. VI we illustrate this procedure for three classes of modified theories of gravity leading to coupled perturbation equations: two models coupling the scalar and tensor modes (dCS gravity [44] and Horndeski gravity [34]) and an EFT model coupling the axial and polar gravitational perturbations [46], where the background QNM spectra are degenerate.
While our expansion is theory-agnostic, we make assumptions about the effect of the modified gravity theory: the background should be perturbatively close to the Schwarzschild metric, and the corrections to the "ordinary" potentials in the GR master equations should be amenable to a power-series expansion in inverse powers of the radial variable.
This results by construction in small corrections to the GR QNM spectra (4). In general, as discussed in Paper I, new nonperturbative frequencies (e.g., quasibound states emerging from zero frequency for massive scalars) may appear in the spectrum, and these are not captured by our formalism.
The assumption that the background is only perturbatively different from the Schwarzschild solution is slightly less restrictive than one might think. For example, in Paper I we showed that slowly rotating Kerr BHs can be accommodated within the formalism. Recent work on higher-derivative corrections to the Kerr geometry [76] and on QNM frequencies of rotating solutions for small coupling [53,54] may allow us to make progress on the calculation of QNMs in modified gravity for rotating BH remnants, such as those observed by the LIGO/Virgo collaboration [77]. Related attempts at parametrizing deviations from the Kerr QNM spectrum [78,79] made use of the connection between the stability of null geodesics and QNMs [80][81][82]. Our formalism may help to clarify the conditions under which this "geodesic correspondence" applies [83,84].

Single-field case
For simplicity, let us first consider adding a friction-like term to the master equation: for some function Z(r * ) and small parameter (note that in this appendix, and here only, we redefine f V → V to simplify the notation). A field redefinition is sufficient to remove the term proportional to the first derivative [85]. The resulting field equation is where the last term in parentheses (proportional to 2 Z 2 ) can be ignored in our perturbative framework. Therefore we can use the formalism of the main text by a suitable redefinition of the radial wave function and of the radial potential. We will now extend this idea to the case of multiple, coupled fields.

Coupled case
Let us now add a matrix of friction-like terms to a coupled set of master equations: for some matrix of functions Z(r * ) and small parameter . For simplicity, here and below we redefine ω 2 − V 1,2 → V 1,2 , we use primes for derivatives with respect to r * , and we consider the case of only two fields. We make a field redefinition and we multiply the resulting equations on the left by This yields where W is a matrix with elements As a result, the field equations no longer depend on the first derivative of any field. When V and Z are diagonal we recover Eq. (A3) at linear order in , and when = 0 we recover the unperturbed QNM spectrum. In general, the removal of friction-like terms introduces new potentials. For these new potentials to fit into our perturbative formalism, all new contributions must vanish at the horizon and (at most) tend to a constant at infinity. One parametrization that would satisfy these requirements is for small parameters β (s) i j , which can be mapped to the parameters α (k) i j discussed in the body of the paper. The summation must start at s = 2, so that the integrals in Eqs. (A8)-(A11) do not diverge as r * → ∞.

Appendix B: Expansion of ω
The coupling parameters α (k) i j appearing in Eq. (3) can, in general, depend on the frequency ω. The computation of the QNM frequency itself depends on the couplings, so a nontrivial ω-dependence of the couplings will affect the QNMs. Consider the implicit expansion of ω about the GR value ω 0 : One may also Taylor expand α k i j (ω) about ω 0 , We assume that α (k) i j and all of its derivatives evaluated at ω 0 are small, so that we can consider them to be of the same order. By substituting (B2) into (B1) we find Finally, substitute (B3) into itself to get the result quoted in Eq. (4): Appendix C: The case of three fields Let us schematically write the coupled master equations for three fields as for some linear operators L 1,2,3 , perturbative potentials δV i j = α i jVi j and small parameters α i j (in this appendix we slightly change the notation to minimize clutter). We assume that the spectrum of L 1 is nondegenerate with the spectra of both L 2 and L 3 . Expanding the fields in powers of α i j to linear order we find φ = φ 0 + α 11 φ 11 + α 12 φ 12 + α 13 φ 13 , (C2) ψ = ψ 0 + α 21 ψ 21 + α 22 ψ 22 + α 23 ψ 23 , (C3) χ = χ 0 + α 31 χ 31 + α 32 χ 32 + α 33 χ 33 . (C4) Note that all of the first-order fields are only functions of the uncoupled zeroth-order fields. Let us write down explicitly some of the field equations: ω 2 φ 12 = L 1 φ 12 +V 12 ψ 0 , (C6) ω 2 ψ 23 = L 1 ψ 23 +V 23 χ 0 .
Recall from Section II B that we numerically integrate the system of equations order by order with respect to a diagonal basis of initial values at the horizon. We also need to specify the value of each perturbation at the horizon, which we set to zero for all but the zeroth-order component of each field, so when α i j = 0 the perturbations are not excited. We will denote by N ψ (φ 12 ), for example, the numerically integrated solution to Eq. (C6). In this notation, the superscript denotes which field is excited at the horizon during the integration.
The purpose of this appendix is to show that the perturbed spectrum of L 1 is calculated using Eqs. (18) and (19), where the determinant S refers to the 1 × 1 or 2 × 2 minors of the 3 × 3 matrix which contain N φ (φ). The proof goes as follows.
In conclusion, the matrix S reduces to and its determinant reads The α i j coefficients of each term can be varied independently, so each term in this sum must vanish. Note also that N ψ (ψ 0 ) 0 and N χ ( χ 0 ) 0, as we have assumed that L 1 is not degenerate with L 2 and L 3 . The first three lines in this sum give the corrections to linear and quadratic order, and are explicitly those formed from the 1 × 1 and the 2 × 2 minors of S containing the (1, 1) element, because we perturb around the spectrum of L 1 . The last two lines vanish identically, because -recalling that N ψ (ψ 0 ) 0 and N χ ( χ 0 ) 0 -the first line implies N φ (φ 11 ) = 0.
In conclusion: adding a third field does not affect our expansion, and the roots of the 1 × 1 or 2 × 2 minors of the matrix S [Eq. (C8)] which contain N φ (φ) can be used to compute the appropriate coefficients. This argument can be extended by induction to the case of N > 3 fields.
An immediate and important corollary of Eq. (C10) is that there are no corrections to S which are linear in α i j with i j. Moreover, by looking at the second and third lines of Eq. (C10) we conclude that couplings give quadratic corrections only if α i j 0 and α ji 0 for i j. In conclusion, a perturbative coupling between nondegenerate operators L i and L j gives at most quadratic corrections to the QNM frequencies.