An analytic approach to quasinormal modes for coupled linear systems

Quasinormal modes describe the ringdown of compact objects deformed by small perturbations. In generic theories of gravity that extend General Relativity, the linearized dynamics of these perturbations is described by a system of coupled linear differential equations of second order. We first show, under general assumptions, that such a system can be brought to a Schr\"odinger-like form. We then devise an analytic approximation scheme to compute the spectrum of quasinormal modes. We validate our approach using a toy model with a controllable mixing parameter $\varepsilon$ and showing that the analytic approximation for the fundamental mode agrees with the numerical computation when the approximation is justified. The accuracy of the analytic approximation is at the (sub-) percent level for the real part and at the level of a few percent for the imaginary part, even when $\varepsilon$ is of order one. Our approximation scheme can be seen as an extension of the approach of Schutz and Will to the case of coupled systems of equations, although our approach is not phrased in terms of a WKB analysis, and offers a new viewpoint even in the case of a single equation.


Introduction
Quasinormal modes (QNMs) describe damped vibrational modes of dissipative systems [2], formally associated with a non-hermitian boundary value problem.They are essential tools in the study of the dynamics of astrophysical compact objects [3] and they encode hydrodynamical properties of strongly coupled systems through the gauge/gravity duality [4].The study of quasinormal modes of black holes was initiated long ago in the early 1970s with the works of Vishveshwara [5], Press [6], Teukolsky [7], Chandrasekhar and Detweiler [8] (see for instance [9] for an introduction), and efficient computational methods are available for the well-known black hole solutions of pure General Relativity (GR) [10][11][12].Renewed interest has been spawned by experimental detection [13] which allows the observational test of theories that extend or modify GR.This includes the study of QNMs in selected theories beyond GR [14][15][16][17][18][19], the development of effective field theories to parametrize and test extensions of GR [20][21][22][23][24][25][26][27][28] and the proposal of new non-perturbative methods that relate the perturbation theory of explicit black hole solutions to Seiberg-Witten gauge theories [29][30][31][32].Recent developments in the understanding of the validity of perturbation theory further raise the hope of precision measurements in the future [33][34][35][36][37].
In this work we present a new analytic approximation scheme for the computation of the spectrum of QNMs of a system of (in general N ) coupled linear differential equations in Schrödinger-like form.Our approach, based on a generalized adiabatic theorem, will also offer a new viewpoint on the regime of validity of the approximation of Schutz and Will [1].Moreover we shall prove that, when an action formulation is available, it is in general possible to put the system of linear equations for the perturbations of a compact object in Schrödinger-like form and we shall describe an algorithm to do so, starting from the quadratic action in frequency space, performing an angular decomposition and some field redefinitions.The problem of finding the QNMs for a system of coupled linear differential equations has received significant attention in the recent literature.Alternative approaches to their study have been devised, including the phenomenological approaches of Refs.[38][39][40], the eikonal approach of Refs.[41][42][43] and the first-order formulation of Refs.[44][45][46].
The results of our analysis can be conveniently summarized in a simple procedure for the computation of the quasinormal modes of a coupled linear system in Schrödinger-like form.We state it in the case of two coupled equations for simplicity.Under the conditions detailed in the paper, the procedure is as follows.Given a 2 × 2 matrix potential V(τ ), where τ is the radial tortoise coordinate, one has to: • find the eigenvalues v i (τ ) of V(τ ) and identify the maxima of v i (τ ), determined by the condition v ′ i (τ i ) = 0.Each physical maximum will be associated with a tower of quasinormal modes.
• for each value of τi , go back to the original formulation in terms of V(τ ) and consider the Taylor expansion of V(τ ) around τi .
• perform a τ -independent change of basis that diagonalizes the constant term in the expansion of V(τ ).The potential takes now the form (4.2), up to row exchange, where x = τ − τi .
• in terms of these coefficients a i , b i , . . ., the leading order quantization condition is given by eq.(4.15), the quadratic mixing correction by eq.(4.25) and the one-loop anharmonic corrections are captured by eq.(4.32).
These analytic formulas for the computation of the quasinormal frequencies, together with the conceptual formulation of the approximation scheme, are the main results of this work.
They can be thought of as a generalization of the classic results of Schutz and Will, which apply for a non-matrix (i.e.1×1) potential [1].It is also worth stressing that our formalism applies also to the case of an N × N potential V. The steps outlined above are applicable to the full set of eigenvalues of V and the only non-trivial question is how to generalize the computation of the mixing corrections.At the order we consider in perturbation theory, it is enough to consider separately the mixing corrections from each off diagonal term (i, j) (with i ̸ = j), which can be extracted by considering the 2 × 2 submatrix and are given by the analogue of eq.(4.25).The corrections to ω 2 from all the mixing terms with j ̸ = i should then be summed.
Even though the leading order quantization condition (4.15) may look formally analogous to that derived for the N = 1 case in Ref. [1], the result of our procedure, already at leading order, is inequivalent to that obtained by just neglecting the mixing terms from the start.Indeed our analysis is necessary to correctly identify the point around which to expand the potential.Ignoring the mixing terms from the beginning will lead to an incorrect identification of the appropriate maximum.In addition, the constant off-diagonal terms are reabsorbed in the third step of the previous procedure: the mixing terms (not necessarily small) are partly taken into account already by the leading order result.Finally, the subleading mixing and anharmonic corrections we provide a systematic way to go beyond the leading order quantization condition.
The article is organized as follows.In Section 2 we introduce the Schrödinger-like equations governing the linear dynamics of perturbations around black holes and review the usual formulation of the boundary value problem associated with QNMs.In Section 3 we then reformulate the problem, first in terms of a classical Hamiltonian dynamical system and then in terms of a quantum mechanical system evolving with a non-hermitian time-dependent evolution operator.We use this formulation to outline our approximation scheme, based on a generalized adiabatic theorem.We then derive analytically the (approximate) quantization conditions for the QNMs for a general system of two coupled equations, in Section 4, including anharmonic corrections up to one-loop (corresponding to a second order WKB analysis) and mixing correction up to quadratic order in the mixing parameter.We validate our approach by studying an explicit example in Section 5, and then draw our conclusions in Section 6.The detailed proofs of some technical results of relevance for the paper are provided in the appendices.In particular, Appendix A contains a proof that the system of equations for the perturbations can be always put in Schrödingerlike form whenever an action formulation is available, and Appendix B contains a proof of the generalized adiabatic theorem for non-hermitian evolution operators acting on a finite dimensional vector space.

Perturbations around black holes
In Einstein gravity the equations governing the linear dynamics of massless spin s perturbations ψ s (r ⋆ ) (s = 0, 1, 2) around black holes can be cast in the form of a Schrödinger-like equation of the form where r ⋆ corresponds to the radial tortoise coordinate, the time dependence of the perturbation has been taken as e −iωt and the angular dependence has been separated through the use of the appropriate angular functions -spherical harmonics for Schwarzschild black holes and spin-weighted spheroidal harmonics for Kerr black holes.The potential function V (r ⋆ ) -which depends on the background, the spin s of the field and the angular momentum ℓ of the perturbation -can be chosen to be real and in the case of gravitational (s = 2) perturbations on the Schwarzschild background is given by the Regge-Wheeler (odd sector) and Zerilli (even sector) potentials; with these choices, in the rotating case, it also depends on the frequency ω.More generally, in the effective field theory approach to black hole perturbations in scalar-tensor theories it is always possible to obtain a Schrödinger-like equation for the perturbations around spherically symmetric [22] or slowly rotating [23] hairy black holes, backgrounds with a non trivial static scalar profile.We show this explicitly and in full generality, under a stability requirement, in Appendix A. However, typically, one has to deal with a system of coupled equations rather than a single equation.This is the case, in particular, for the even sector of a spherically symmetric background, where the scalar degree of freedom is coupled to the even gravitational perturbation.
Our goal is to develop an analytic approximation scheme to compute the quasinormal frequencies of a Schrödinger-like system of N coupled equations in terms of a finite number of coefficients, generalizing the approach first introduced in Ref. [1].

The quasinormal modes
The quasinormal modes are defined by the following boundary value problem: given a smooth N × N matrix function V(r ⋆ ) (the potential) defined for −∞ < r ⋆ < ∞ and satisfying the fall-off condition we want to find the values of ω such that the Schrödinger-like system of equations admits a solution satisfying the asymptotic conditions Physically, these conditions correspond to an outgoing (right-moving) wave at spatial infinity r ⋆ → +∞ and to an ingoing (left-moving) wave at the black hole horizon r ⋆ → −∞.
The special values of ω that admit such a solution are known as the quasinormal frequencies ω qnm and the corresponding functions are the quasinormal modes.Due to the boundary conditions, this boundary value problem is not self-adjoint, so that the quasinormal frequencies can be (and usually are) complex numbers and the corresponding functions are usually not normalizable.We shall assume that the eigenvalues of the potential V(r ⋆ ) are real.In some instances the potential can be a function of the frequency ω itself; in such a case we shall assume that the eigenvalues are real for real ω.In what follows we shall work under the simplifying assumption of a potential independent of the frequency and then comment on how to generalize our methods to the case of frequency-dependent potentials.In Sec. 5 we treat an example in which the eigenvalues of V are complex in a range of r ⋆ and discuss the regime of validity of our approach.Notice that, as we show in Appendix A, it is always possible to obtain a Schrödinger-like system with hermitian potential for the perturbations of a stationary spacetime background if an action formulation is available.In those cases the eigenvalues of the potential will be real (for real ω), and our approach will be thus applicable.Moreover, as it will become clear in the following sections, the validity of our approach (and of that of Schutz and Will [1] as well) is related to the smoothness of the potential.Examples with perturbations, such as those of Ref. [47], can violate this assumption, since, depending on the shape of the perturbation (e.g. a bump), the derivative of the evolution operator can be large even for perturbations of arbitrarily small amplitude.However, as discussed in [48,49], over short time scales the ringdown is described by the QNMs of the unperturbed potential and the use of our approach is justified.

The Hamiltonian formulation and a generalized adiabatic theorem
The system of differential equations (2.3) admits an Hamiltonian formulation.In order to make contact with the standard notation of Hamilton's mechanics, in this section we shall map the variable r ⋆ to a fictitious time coordinate τ and the fields ψ i to the generalized coordinates q i , and introduce the canonical conjugate momenta p i .The time dependent Hamiltonian H(q, p, τ ) = 1 2 reproduces the system of equations (2.3), as can be easily checked by using Hamilton's equations In this language it becomes transparent that it is not possible in general to bring the system of differential equations (2.3) in diagonal form through a point transformation Q i = f i (q j , τ ).Indeed, such a change of variable corresponds to a canonical transformation with generating function F 2 (q i , P i ) = P i f i (q j , τ ), so that the generalized coordinates and the Hamiltonian transform as (3. 3) The term H(q(Q), p(P, Q)) in the new Hamiltonian does not include terms linear in the generalized momenta P i , whereas the term ∂F 2  ∂τ is linear in the momenta, being equal to P i ∂f i (q j ,τ ) ∂τ . Therefore, the point transformation always introduces a mixing term between the new coordinates and momenta, unless f i (q j ) is time independent.But a time independent coordinate transformation can diagonalize the potential term q i V ij (τ )q j for all times τ only if the latter has the simple form V ij (τ ) = v(τ ) Vij , with constant Vij , which is not the case in general.
Let us go back to the boundary problem and try to gain further insights in the Hamiltonian language.Hamilton's equations written in matrix form are d dτ which, when multiplied by the imaginary unit i, is formally the time-dependent Schrödinger equation with a non-hermitian evolution operator that we denote with K. Since |τ V ij (τ )| → 0 for τ → ±∞, it is natural to perform a (complex) canonical transformation that diagonalizes the time evolution in the asymptotic region. 1 Defining the evolution equation becomes i d dτ while the boundary conditions are simply In this basis it becomes transparent that for ω real in the asymptotic region |τ | → ∞ the evolution is unitary and governed by the asymptotic eigenvalues ±ω.The boundary conditions correspond to the requirement that the system undergoes a transition from an eigenstate with eigenvalue +ω at early times to one with eigenvalue −ω at late times.
In the case of a slowly varying Hermitian evolution operator, the adiabatic theorem of quantum mechanics implies that transitions occur efficiently only when two eigenvalues become (almost) degenerate, as in the well-known Landau-Zener effect (i.e.close to a level crossing, usually lifted by perturbations in a quantum mechanical system).
We shall establish a mild generalization of the adiabatic theorem for evolution operators that are diagonalizable, have a real spectrum, and act on a finite dimensional vector space.

An adiabatic theorem
The evolution operator K can be regarded as an operator acting on the 2N -dimensional (complex) vector space with elements (q, p) T .Let us assume that the potential V(τ ), acting on the N -dimensional q subspace, is a diagonalizable operator with real eigenvalues v i , i = 1, . . ., N .2Then the evolution operator is diagonalizable for every time τ , with eigenvalues determined by the zeros of the characteristic polynomial: where in the second equality we used that the identity matrix commutes with every operator.From eq. (3.9) it follows that the eigenvalues occur always in pairs and are such that their square is an eigenvalue of (ω 2 1 N − V).In particular: Since the operator K is not Hermitian its eigenvectors |k a ⟩ (with a = 1, . . ., 2N ) do not form an orthonormal basis with respect to the standard inner product; however, they still form a basis.Moreover, since the Hilbert space is finite dimensional we can consider the dual basis ⟨ ka |, with a = 1, . . ., 2N , defined by which exists and is unique.The dual basis vectors act as projectors on the desired eigenspaces.By the uniqueness of the dual basis, it follows that given an eigenvector |k a ⟩ of K with eigenvalue λ a , the dual vector ⟨ ka | is a left eigenvector with eigenvalue λ a : As a last ingredient for the proof of the adiabatic theorem we notice that for any fixed value of ω 2 > 0, at early times τ → −∞ the eigenvalues λ i,± are real, for every i.This follows simply from the eigenvalue equation (3.10), the assumption that ω and v i (τ ) are real and the fall-off condition (2.2), which implies v i (τ → −∞) → 0. By continuity, the reality of the eigenvalues can be violated only after a level crossing: as v i (τ ) → ω 2 the two eigenvalues λ i,± → 0 and become degenerate.
As we discuss in detail in Appendix B, these preliminary considerations ensure that the standard proof of the adiabatic theorem of quantum mechanics goes through for this system, even though the evolution operator K(τ ) is not Hermitian.As a result, away from level crossing and for smooth potentials, level transitions are suppressed.
In our argument this is the only step where the assumption of real ω comes into play: the generalized adiabatic theorem, valid for ω real, allows us to single out level crossing as special points to focus on.

A heuristic argument on level crossing and quasinormal modes
The boundary conditions for the quasinormal modes (3.7) require a level transition from an eigenstate with eigenvalue +ω to one with with eigenvalue −ω.From the generalized adiabatic theorem of the previous section it follows that for well-behaved potentials, with slow variation and real eigenvalues, such a transition occurs efficiently only if a pair of eigenvalues λ i,± (τ ) (one positive and one negative) becomes degenerate, by touching or crossing 0 at a given time τi .Thinking about a fixed value of i for definiteness, the condition for this to happen is: Starting from λ 2 i,± (τ ) > 0 at τ → −∞ the eigenvalues approach 0 as τ → τi .Assuming the v i (τ ) is a smooth function of τ this can happen either with To gain some intuition we can consider the one-dimensional case N = 1.In this case v i (τ ) ≡ V (τ ) and the condition V ′ (τ ) = 0 corresponds to the well-know condition that the real part of the quasinormal frequencies ω 2 is equal to the maximum of the potential, familiar from the approximation scheme of [1].Heuristically this can be justified by noticing that the quasinormal boundary conditions can be formulated requiring the outgoing waves to have equal amplitudes both at the horizon and at spatial infinity.In the analysis of [1], based on matching the ingoing and outgoing WKB asymptotics, this condition can not be accommodated if the frequency ω 2 (or its real part) is smaller than the maximum of the potential.Indeed in this scenario there would be a region with exponential damping, making it impossible to satisfy the boundary conditions.We shall make the assumption that also in the multidimensional case, N > 1, only the crossing points with v ′ i (τ i ) = 0 are relevant for the computation of the quasinormal frequencies.
For well-behaved physical potentials, the total number of real maxima in the physical coordinate region is expected to be N , the same as the dimension of the system of linear equations.In this case, every maximum v i (τ i ) with i = 1, . . ., N will give rise to a tower of quasinormal modes.
Up to now we have been assuming that ω is real and that the eigenvalues v i (τ ) are real.The former should be understood as an approximation, under the condition that the imaginary part of ω is much smaller than the real part.The latter condition can be relaxed in physically realistic examples, again on a heuristic basis.

An approximation scheme for computing quasinormal frequencies
Having determined the interesting points τi as the values of τ that correspond to stationary points of the eigenvalues of the potential matrix, defined by the condition v ′ i (τ i ) = 0, we shall focus on the dynamics of the system by going back to the original formulation (2.3) and expanding the potential in a power series around every maximum v i (τ i ).Moreover, for notational convenience we change independent variable from τ to x = τ − τi , for every fixed value of i. 3To illustrate the procedure let us focus on the case N = 2 and choose i = 1 for concreteness.Consider the Taylor expansion of the matrix V(x) around x = 0, corresponding to a maximum of v 1 (x): As a first step we bring the constant matrix V (0) to diagonal form by performing a constant change of variables. 4In this basis, the potential takes now the form where the absence of the term a 1 follows from the fact that we are expanding around a local maximum of v 1 (x), while b 0 and c 0 have been reabsorbed by the constant diagonalization.The coefficients d 1 , b 2 , c 2 , . . .are in general all non-zero, but will give subleading corrections to the quasinormal frequencies, as we shall explain in the following.We shall denote the first diagonal component of V(x) in this particular basis as V 11 (x).
In a situation where the mixing terms b k , c k (with k ≥ 1) are small we can adopt a perturbative approach and introduce a formal counting parameter ε.This is the case, for instance, when they arise from a weakly coupled extension of GR, in which they are proportional to the coupling constants of the underlying theory.Notice, however, that the effect of b 0 , c 0 has been fully captured in our approach by diagonalizing the constant potential matrix at the maximum; such an effect is not necessarily small, so that our approach can allow to compute the quasinormal modes also in scenarios in which the mixing terms are of order one, but the residual coupling after diagonalization of the constant matrix are small as we shall see in the next section.
The leading order approximation for the tower of quasinormal frequencies associated with v 1 (τ 1 ), is obtained by neglecting the mixing terms and considering only the dynamics of q 1 (x), as governed by V 11 (x).At lowest order, the behaviour of the quasinormal modes can then be captured by the analysis of [1], considering only the quadratic approximation for the potential V 11 (x) ≃ a 0 − a 2 x 2 .Anharmonic corrections to V 11 (x) can be included by either performing a higher order WKB analysis of the asymptotics [11,51,52] or by mapping the question to a bound state problem [12,[53][54][55][56].In order to include also the corrections deriving from the mixing terms it proves useful to adopt the latter approach, which allows both for a conceptually clear formulation of the computation and a straightforward method of calculation.We shall adopt a perturbative approach and compute separately corrections of two classes: mixing terms among q 1 and q 2 , and anharmonic terms in the potential.

Mixing terms
Let us consider without loss of generality the system of equations obtained by expanding around the maximum of v 1 (x), and keep only terms that are formally of order x 2 or lower with respect to V 11 (x).That is: if we formally assign units to x, we only keep terms whose contribution to the tower of quasinormal frequencies associated with V 11 (as defined in the limit ε → 0) is of the same order as the contribution generated by a 2 .As it will become clear (see eq. (4.25)), mixing contributions of this order are generated only by b 1 , c 1 , d 0 .
Introducing explicitly the counting parameter ε, we have the system with boundary conditions such that q i is left-moving at the horizon and right-moving at infinity, as in (3.7).Since we are studying the set of quasinormal modes associated with v 1 (τ 1 ), we are interested in the eigenmodes that reduce to the q 1 eigenfunction in the limit ε → 0. To understand the asymptotic behaviour of the system we can use the ordinary WKB approximation for a system of coupled differential equations [57], considering the limit x → ±∞. 5 At leading order in the WKB expansion we look for solutions of the form with A = B = const and S ′′ (x) ≪ (S ′ (x)) 2 .Neglecting terms proportional to S ′′ (x), our system of coupled differential equations becomes an algebraic system in u ≡ S ′ (x): The system admits non-trivial solutions when it becomes degenerate; the vanishing of the determinant gives the conditions (4.6) In the limit x → ±∞ the two WKB solutions then become: From the asymptotic analysis we learn that the eigenfunctions that satisfy the quasinormal mode boundary conditions and reduce to the q 1 eigenfunction in the limit ε → 0 are those associated with u + : Moreover, the asymptotic behaviour of the solutions suggests performing the analytic continuation x → e +i π/4 z, so that the eigenvalue problem (4.3) is mapped to a new eigenvalue problem with regular boundary conditions: ⃗ q(z) → 0 z → ±∞.(4.10) We shall work under the assumption that the analytic continuation can be performed without encountering singularities in the complex plane, so that the procedure allows to recover the quasinormal modes from the eigenvalues of the associated bound state problem.
A similar (but distinct) procedure has been discussed in the case of a single differential equation in Schrödinger-like form in Refs.[53][54][55][56].
Multiplying by −i/2 and defining we obtain the system Working perturbatively in ε we first of all recover the leading order formula of Schutz and Will [1].Indeed, working at zeroth order in ε, the eigenstate with q (0) 2 = 0 satisfies the Schrödinger equation and the boundary conditions for the bound states of a quantum mechanical harmonic oscillator. 6The eigenfunctions and eigenvalues are then easily obtained: where H n (•) are Hermite polynomials and we normalize the eigenfunction to have unit L 2 norm, with Expressing the eigenvalues in terms of the original variables we recover the quantization condition for the quasinormal modes of [1], without performing a WKB matching analysis: In order to compute the leading ε-corrections to the quantization condition we use the well-known Rayleigh-Schrödinger perturbation theory approach, familiar from the timeindependent eigenvalue problems of ordinary quantum mechanics.Since the perturbation is off-diagonal, the first non-trivial correction arises only at second order in perturbation theory.
We denote by δH ε the order ε perturbation in equation (4.12): Going through the usual derivation, the order ε 2 correction to the eigenvalue ν n can be expressed as n (z).(4.17) In the ordinary treatment given in quantum mechanics textbooks it is customary to express ⃗ q (1) n (z) as a linear superposition of zeroth-order eigenfunctions ⃗ q (0) m (z).In our case it will be convenient to solve directly the differential equation for the first-order eigenfunction ⃗ q (1) n (z).Since q (0) 2,n (z) = 0 and the perturbation δH ε is off-diagonal we shall only need q (1) 2,n (z).Working at first order in ε, the differential equation for q where we have defined α n so that its real part is positive and We can treat the term on the right hand side as a source term and solve the differential equation using the Green's function: (4.20)The Green's function for the differential operator d 2 dz 2 − α 2 n satisfying the appropriate boundary conditions is given by: from which it follows that Plugging back in eq.(4.17) and using eq.(4.13), after straightforward manipulations we obtain where P n (σ) is a polynomial of degree (2n+2) with only even terms, defined by the integral where the function F n (y) defined by is known analytically and tabulated for the first few values of n in Table 1 together with the polynomials P n (σ).An asymptotic approximation for F n (y) can be obtained in the limit |y| → ∞, as detailed in Appendix C.Moreover, it is possible to check that |F 0 (y)| < 1 for Re(y) > 0, so that the mixing correction to the fundamental mode is always bounded by |ε 2 b 1 c 1 /a 2 | and the corresponding perturbative expansion is under control as long as this parameter is small.In particular, the approximation scheme holds also for cases where ε ∼ O(1), provided that |b 1 c 1 /a 2 | is small, which depends on the potential.

Anharmonic corrections
Higher order corrections to the quadratic potential can be included systematically and the effect of these terms can be captured straightforwardly in a perturbative approach.In this spirit, in this section we shall consider only the anharmonic corrections to the leading result (4.15), neglecting the mixing terms analyzed in the previous section.Following the procedure outlined at the beginning of Section 4 we focus on the single differential equation We consider the Taylor expansion of V 11 (x) around x = 0, using for convenience the following notation, consistent with the previous sections: Since we are expanding around a maximum, the coefficient a 2 is positive in our sign conventions.
The asymptotic WKB analysis of the previous section and the analytic continuation we performed relied on the quadratic form of the potential, so the latter is no longer warranted in this case.However, as recently pointed out by Hatsuda [12], it proves convenient to perform a different analytic continuation which allows us to include efficiently anharmonic corrections up to very high orders.Following [12] we introduce a deformation parameter ℏ such that: For the value ℏ = i we recover our original problem, while ℏ has formally the same role as Planck's constant in the quantum mechanical Schrödinger equation for a particle in the inverted potential −V 11 (x).Defining x = x/ √ ℏ, multiplying by 1/2 and rearranging the terms, the equation becomes which makes it manifest that ℏ acts as a loop counting parameter for anharmonic interactions.The quasinormal modes are obtained as the analytic continuation at ℏ = i of the eigenvalues for the bound state problem (4.30). 7y dimensional analysis, the anharmonic corrections to the leading order eigenvalues will be a function of the dimensionless coupling constants We shall assume here that the dimensionless coefficients are small so that the perturbative expansion is under control, which is a valid assumption in the explicit examples we consider.
A straightforward method to compute corrections to the eigenvalues from anharmonic terms would be to use Rayleigh-Schrödinger perturbation theory.A more efficient approach was devised long ago by Bender and Wu [59], as recently emphasized in this context in Ref. [12].This amounts to a systematic expansion in ℏ -or equivalently in the number of loops -and reproduces the quantization conditions obtained through the approach of [1,51], which is also an expansion in ℏ (although with a different formulation).The general formulation for a bound state problem of the form (4.30) and the resulting recurrence relations have been discussed in detail in Ref. [60].The 1-loop quantization condition, corresponding to a second order WKB analysis, is given by [12,51] Comparing with the leading order quantization (4.13) we see that the relative correction ν n depends on the expansion parameters (a 3 , a 4 ) as expected.We also notice that the relative correction grows linearly with n in the limit of large n, so that the approximation is under control only for n small, in a range controlled by the dimensionless couplings (a 3 , a 4 ).Corrections from quintic or higher interactions enter only at higher loop (or WKB) order.

Explicit example
In general, the linearized dynamics of the perturbations around a given black hole solution, described by the system (2.3), cannot be recast in the form of a set of fully decoupled Schrödinger-like equations.This is the case, for instance, of the perturbations of charged Kerr-Newman black holes (see for instance [9], section 111), of massive and partially massless spinning fields on Schwarzschild or Kerr spacetimes [61][62][63], 8 and of perturbations of hairy black holes beyond GR, such as in scalar-tensor theories where the scalar has a nontrivial background profile (see, e.g., Refs.[15,22,23,45,68,69]). 9 Additional examples of coupled systems of equations arise in the study of the hydrodynamics of strongly-coupled gauge theories through the gauge/gravity duality, see for instance [72,73].
In the following, we shall consider an explicit example for the potential matrix V in (2.3), and we shall use our analytical approach to compute explicitly the quasinormal modes, and compare them against the values obtained with numerical methods.We stress that the choice of example is non-generic and rather peculiar, and is selected to push the approximation to its extreme condition.It does not correspond to the typical situation, but should be instead interpreted as a worst case scenario.This explicit system is simple enough that it allows analytic control, yet, as we shall discuss, it reduces to well-known results in the ε = 1 limit, allowing us to explore the regime of large ε and to cross-check our computations with previous results in the literature. 10In addition, it is an instructive example which allows us to discuss possible issues associated with the method in the case in which the potential is not in a hermitian form.
The toy model we shall focus on is where the radial coordinate r is related to the tortoise variable r ⋆ in (2.3) through the relation dr ⋆ /dr = 1/f (r).It can be checked that this potential satisfies the adiabatic condition of Appendix B, away from level crossing.ε is a control parameter that we shall vary within the interval [0, 1].This will allow us to explore different cases, from a regime of perturbative mixing (ε ≪ 1) to a regime where the coupling is order one at r ∼ r s (ε ∼ 1).Note that, when ε = 1, eq. ( 5.1) recovers the potential of parity-odd partially massless spin-2 fields on a Schwarzschild-de Sitter spacetime in the limit of zero cosmological constant 8 Partially massless fields are special irreducible representations of massive spinning particles on Einstein spacetimes [64][65][66] (the case of more general spacetimes has been discussed e.g. in [67]).They exist for particles with spin greater than 1 when their masses take very particular values.Their main property is that they carry one degree of freedom less than generic massive fields with the same spin, thanks to a residual gauge symmetry that is responsible for removing the helicity-zero component. 9If the scalar background is constant or, in particular, vanishing, it is always possible to decouple the equations for the perturbations [21].Another example where it is possible to diagonalize the potential matrix V in (2.3) through a coordinate independent linear transformation is given by charged dilaton black holes [70,71]. 10To test our method on explicit examples of scalar-tensor theories, like that of a scalar-Gauss-Bonnet model, we would need the explicit solution of the background.However, this is known analytically only perturbatively in the coupling; as a consequence, it would allow us to test only the small ε regime.2: Comparison between the analytic estimate and the numerical computation for the first set of QNMs for (n = 0, ℓ = 2) perturbations, in the toy model described by the potential in eq.(5.1).We include both the mixing corrections of Sec.4.1 and the anharmonic (one-loop) corrections of Sec.4.2.This set of quasinormal modes corresponds to the local maximum located in the region r/r s > 3/2, where the eigenvalues of V(r) are real and the approximation scheme we use is justified.

Re(ωr
(see, e.g., eq.(4.13) of Ref. [63]).In this limit, the quasinormal modes are expected to recover the odd spectra of massless spin-1 and spin-2 fields, as we shall check explicitly for various values of ℓ.The toy example (5.1) will also give us the opportunity to discuss possible subtleties and limitations of our approach.
Our results are summarized in Tables 2 and 3, where we compare the quasinormal modes obtained from our analytic method, including both mixing and anharmonic corrections, with the ones computed numerically. 11ome comments are in order here.Let us start by looking at Table 2, where we report the first set of frequencies with n = 0 and ℓ = 2.Note that, as anticipated above, when ε = 1 we recover the quasinormal modes of spin-1 fields on a Schwarzschild spacetime [61,63,75].The precision of the analytic estimates is quite remarkable, with an error less than a few percent in most cases and ≲ O(0.1%) for the real part of ω for small values of ε.The percent error with respect to the numerical values can be further reduced by including higher orders corrections.The case of the second set of modes, reported in Table 3, is instead slightly different.As it can be seen from the table, the error with respect to the numerical values becomes large as ε increases, signaling a breakdown of the analytic approximation.The reason for this is that, in the particular example under consideration,  3: Comparison between the analytic estimate and the numerical computation for the second set of QNMs for (n = 0, ℓ = 2) perturbations, in the toy model described by the potential in eq.(5.1).We include both the mixing corrections of Sec.4.1 and the anharmonic (one-loop) corrections of Sec.4.2.This set of quasinormal modes corresponds to the local maximum located in the region 1 < r/r s < 3/2, where the eigenvalues of V(r) are complex close to level crossing or for large ε.The approximation scheme is not justified and the results for the imaginary part of the quasinormal frequencies are off for sizeable mixing parameter ε, especially for the imaginary part.
the eigenvalues of the potential (5.1) acquire a nonzero imaginary part on a finite range of r.For small values of ε, this range is small and the local maxima are located in the region where the potential eigenvalues are real.The approximation is therefore well defined and provides an accurate estimate of both the two sets of quasinormal frequencies for small ε.However, as ε increases, the maximum of one of the two eigenvalues of V(r) falls close to the boundary of the region where the eigenvalue is complex, approaching the edge until it is lost.The approximation scheme applied to this eigenvalue is no longer justified for sizeable values of ε and cannot be used to estimate the corresponding set of quasinormal normal modes. 12This explains why the frequencies in Table 3 are significantly off for sizeable mixing parameter ε, especially for the imaginary part.We should stress again that what happens to the modes in Table 3 is not the typical situation.In fact, as we show explicitly in Appendix A, when an action formulation is available it is in general possible to recast the N coupled equations for the perturbations in a Schrödinger-like form (2.3) with hermitian potential.As a result, the eigenvalues of the new V(r) are real and the breakdown shown in Table 3 does not happen.In particular, if N = 2 and the potential is frequency independent, the required transformation can be always found analytically and the hermitian potential can be written in closed form.[75].
explicitly for the particular case of the system (5.1) with ε = 1.
However, we find the specific example (5.1) still particularly instructive.Even if not generic, it is useful to show what the limitations of our analytical approximation scheme might be when a closed-form expression for the hermitian potential is not available, or when an action formulation is not available and our proof in Appendix A does not apply.For a generic two by two potential matrix it is easy to show that the eigenvalues are always real when the off-diagonal elements have the same sign, and can develop an imaginary part only when −b • c > (a − d) 2 /8.This explains the properties previously emphasized for the QNMs of our toy example (5.1), since the quantity b • c changes sign for r/r s < 3/2 and its size is controlled by ε, elucidating the results of Tables 2 and 3.
In Table 4 we compare, as a function of ℓ and for ε = 1, the analytic computation for the n = 0 quasinormal mode associated with the local maximum located in the region r/r s > 3/2, with the numerical value of the n = 0 quasinormal frequency of a spin-1 field on a Schwarzschild spacetime [75].As discussed previously, in this limit the system (5.1)recovers the odd QNM spectrum of massless spin-1 and spin-2 fields on a Schwarzschild background and the spin-1 field is the one associated with the set of modes for which our approximation scheme is valid.The accuracy for the real part is remarkable and reaches the sub-percent level in the eikonal limit, whereas for the imaginary part it asymptotes to around 10%, signaling a systematic shift probably induced by the approximate treatment of the off-diagonal contributions, which are sizeable for ε = 1.Note that this can well be a feature of the field basis chosen to write the potential.In some cases convergence can be improved by writing the potential in a hermitian form, using the general procedure outlined in Appendix A below.An explicit example of this is precisely given by partially massless spin-2 fields on Schwarzschild-de Sitter spacetime: in [63], it is shown that, in the basis where the potential is hermitian, the off-diagonal mixing terms go as ∼ 1/ℓ in the large-ℓ limit.We thus expect in this case a better precision of the estimate of the imaginary part of the frequencies in the eikonal limit.
We conclude mentioning that a different approach to estimate the QNMs in the eikonal limit for a system of coupled linear equations has been discussed in Refs.[41][42][43].We stress though that our method is more general as it applies, under the assumptions discussed above, to any value of ℓ.

Discussion
In this work we have introduced a new analytic approximation scheme to estimate the QNM frequencies for systems of coupled linear differential equations of second order in Schrödinger-like form.Examples where coupled equations arise include massive and partially massless spinning fields on Schwarzschild-(anti-)de Sitter spacetimes, and perturbations around black holes beyond GR, e.g. in the context of scalar-tensor theories.Moreover, as proved in Appendix A, whenever an action formulation is available the system of equations can be put in Schrödinger-like form with a hermitian potential.
The approximation scheme is always under control when the mixing corrections are perturbative and small.However, we stress that our method can be applied equally well to cases where the mixing is large and provides up to order-one corrections close to the black hole, as long as an emergent perturbativity condition for the expansion of the potential around a critical point is satisfied, as discussed at the end of Section 4.1.Another main advantage of this approach with respect, for instance, to more phenomenological parameterizations of the ringdown in theories beyond GR (see e.g.Ref. [39]) is that it allows to extract the quasinormal frequencies from the shape of the potential matrix around the maxima of its eigenvalues.As a consequence, with a finite experimental accuracy, the knowledge of the full shape of the potential is not needed to make accurate predictions. 14n addition, the potential does not need to be expressed as a series of (inverse) power corrections to the GR potential, but can have an arbitrary functional form.The only assumptions are that V(r ⋆ ) in (2.3) is a smooth function of the radial tortoise coordinate r ⋆ , satisfying the fall-off condition (2.2) and having real eigenvalues.In fact, our method is even more general as it applies even to cases where the potential depends also generically on the frequency, V ≡ V(r ⋆ , ω), and has real eigenvalues for real values of ωsee Refs.[23,63] for some explicit examples.In such cases, one can formally go through the same procedure described in Section 4, regarding the various quantities -such as for instance the potential eigenvalues -as functions of the frequency, and at the end solve self-consistently for the quasinormal modes.
The main results of our analysis are the leading order quantization condition (4.15), the quadratic mixing correction (4.25) and the one-loop anharmonic correction (4.32).The analysis of anharmonic corrections could be easily extended at higher loop orders using the methods of [59,60].It would be interesting to extend these methods to the analysis of a system of equations and investigate the Borel summability of the perturbative series along the lines of [12].
Our method is well justified when the eigenvalues of the potential matrix are real, and when the imaginary part of the frequency ω is small compared to its real part.The first condition can be always accommodated (for real values of the frequency), as long as the potential matrix is chosen to be hermitian, which can be always done when an action formulation is available, as previously discussed.In order to validate our approach and show its limitations we considered an example modeled on the case of partially massless spin-2 fields on a Schwarzschild-de Sitter spacetime in the limit of zero cosmological constant studied in Ref. [63], using a non-hermitian potential.The accuracy, when the approximation is justified, is remarkable and it is expected to improve when using the hermitian formulation.
The fact that our approximation scheme does not depend on the full shape of the potential away from the maxima of the eigenvalues makes it particularly suitable when combined with an effective field theory approach [22,23].It would be interesting to generalize the light-ring expansion introduced in [22] for a single equation to coupled systems, and understand more systematically how to connect possible deviations in the observed ringdown frequencies to the effective couplings using our approximation scheme.We leave this and related aspects for future work.

B Proof of the generalized adiabatic theorem
In this appendix we present the detailed proof of the generalized adiabatic theorem described in section 3.1.The proof follows the standard approach of the quantum mechanical adiabatic theorem (see for instance [58]), but does not rely on the hermiticity of the evolution operator.
Theorem.Consider a finite dimensional system described by a time-dependent Schrödinger equation with (possibly non-Hermitian) evolution operator K(τ ).Assume moreover that the spectrum of K is real and non-degenerate for a range of τ ∈ [τ i , τ f ], with instantaneous eigenvalues λ a (τ ) in increasing order for a = 1, . . ., 2N and corresponding eigenstates |k a (τ )⟩.Then, in the limit ∂ τ K(τ ) → 0 and away from level crossing, if the system at time τ i is in the instantaneous eigenstate |k a (τ i )⟩ of K(τ i ) for some fixed a, at time τ f it will be in the eigenstate |k a (τ f )⟩, up to a scalar factor.
Proof.Consider the basis of instantaneous eigenstates of K(τ ): The term e iθ b (τ ) is an adiabatic factor that in the quantum mechanical case, corresponding to Hermitian K(τ ), can be shown to be a pure phase and gives rise to Berry's phase when a parameter of the system is varied adiabatically and cyclically.

C Asymptotic approximation for the mixing correction
The function appearing in the mixing correction to the quasinormal modes in equation (4.25) admits an asymptotic expansion in the limit |2y| → ∞ under the assumption that |arg(2y)| < π 2 .Indeed, expanding the function P n (σ)e −σ 2 in powers of σ around σ = 0, P n (σ)e −σ 2 = ∞ k=0 a n,k σ k , and applying Watson's lemma (see e.g.Theorem 15.2.8 of [77]) it follows that in the limit |2y| → ∞ (for |arg(2y)| < π 2 ): The polynomial P n (σ) has only even powers of σ and has degree 2n + 2: P n (σ) = P n,0 + P n,2 σ 2 + • • • + P n,2n+2 σ 2n+2 , (C.3) therefore the coefficients a n,k with odd k are all identically zero and the first few coefficients a n,k can be expressed in term of the P n,k as follows: a n,0 = P n,0 , a n,2 = P n,2 − P n,0 , a n,4 = P n,4 − P n,2 + P n,0 /2, . . .(C.4) The leading term of the asymptotic expansion is proportional to P n,0 , which can be computed explicitly from the definition of P n (σ) to be from which it follows that F n (y) ∼ n + 1 2 1 y 2 + . . ., (C.6) so that in this limit the mixing correction to the quasinormal modes is approximately: Whether this asymptotic expansion provides an accurate approximation or not depends on the value of the parameters characterizing the potential and should be checked on a case by case basis.

Table 1 :
Explicit form of the polynomials P n and functions F n defined respectively in eqs.(4.24) and (4.26), relevant for the computation of the mixing correction (4.25). Table s )

Table 4 :
Comparison between the analytic estimate and the numerical computation for the fundamental QNM n = 0 of a spin-1 field on a Schwarzschild spacetime computed from the potential in eq.(5.1) with ε = 1, as a function of ℓ.We include both the mixing corrections of Sec.4.1 and the anharmonic (one-loop) corrections of Sec.4.2.The numerical values are from Ref.