Nonlinear quasi-normal modes: uniform approximation

Recent works have suggested that nonlinear (quadratic) effects in black hole perturbation theory may be important for describing a black hole ringdown. We show that the technique of uniform approximations can be used to accurately compute 1) nonlinear amplitudes at large distances in terms of the linear ones, 2) linear (and nonlinear) quasi-normal mode frequencies, 3) the wavefunction for both linear and nonlinear modes. Our method can be seen as a generalization of the WKB approximation, with the advantages of not losing accuracy at large overtone number and not requiring matching conditions. To illustrate the effectiveness of this method we consider a simplified source for the second-order Zerilli equation, which we use to numerically compute the amplitude of nonlinear modes for a range of values of the angular momentum number.

1 Introduction Quasi-normal modes (QNM) of black holes (BH) are a highly valuable tool for accurately describing the last part of a gravitational wave (GW) signal, known as the ringdown phase [1].Far from the BH, the metric perturbation is expressed as a sum of exponentially decaying oscillations.One remarkable result of BH perturbation theory is that the frequency of these oscillations is quantized, meaning that it depends on a integer «overtone number» n and on the spherical harmonic numbers (ℓ, m) [2].While QNMs do not form a complete basis for describing the full late-time signal of a BH collision [3], they still fit very well the results from numerical relativity [4,5].
Various techniques have been developed to compute QNM frequencies with high precision.All of these techniques rely on studying the simplified problem of a monochromatic linearized perturbation on Schwarzschild or Kerr spacetime obeying boundary conditions infalling at the BH horizon and outgoing at infinity.To date, the most accurate of these is Leaver's continued-fraction method [6].Another well-known approximation technique is an application of WKB method to the QNM problem [7][8][9][10][11].This last method has the advantage of being analytic, thus providing insights into the mechanism of QNM generation and allowing generalizations to beyond-GR theories [12,13].Other techniques include the Bender-Wu approach [14], the monodromy method [15] (see e.g.[1] for a review), methods based on a spectral decomposition of the equations in hyperboloidal slices [16,17] and more recent developments related to supersymmetric and conformal field theories (see e.g.[18][19][20][21]).Finally, closely related to the WKB method, the method of uniform expansions was developed in [22], where it was applied together with Borel-Padé summation.
In this article we will also rely on the method of uniform expansions (see e.g.[23,24]), performing detailed computations of the linear frequencies for various n and ℓ and extending it to compute the amplitude of nonlinear QNMs.We will provide a detailed explanation of this technique in section 3. Interestingly, this method allows to circumvent the matching conditions traditionally used in the Schutz and Will WKB approach [9].Indeed, our approximate QNM wavefunction is accurate both near the maximum of the potential and at large distances from it.Similar to the WKB method, this approach does not strongly depend on the precise form of the potential in GR, making it potentially applicable to other modified-gravity settings.However, unlike the WKB method, we will demonstrate in section 4 that its accuracy remains high even as the overtone number n increases.We will achieve this without the need for high-order matching, resulting in more compact expressions compared to the WKB method.
Our primary motivation for introducing this new analytic technique is related to the study of nonlinearities in BH ringdown.Indeed, since the merger of two BHs is a highly nonlinear process, it is not surprising that the ringdown phase may not be entirely described by a linearized perturbation on the BH background.Several works have highlighted the significance of nonlinearities in BH ringdown [25][26][27][28][29][30][31].In particular, due to mode-coupling effects at quadratic order in perturbation theory, there exist a set of «nonlinear quadratic QNMs» whose frequencies are given by the sum or difference of linear QNM frequencies [31].These modes emerge at the second order in perturbation theory because the Regge-Wheeler and Zerilli equations involve a source term with a product of two linear QNMs [32][33][34].Notice, however, that the complete gravitational-wave signal after a BH merger does not only consist in a superposition of QNMs, as it also contains a flat-space and tail piece both at first and second-order [3,31,[35][36][37]; our work only concentrates on the QNM part of the Green's function.
Including these nonlinearities into ringdown models could prove to be advantageous for fitting the signals.Given the anticipated increase in precision of BH spectroscopy made possible by the space-based interferometer LISA [38], it is crucial to gain a quantitative understanding of the structure of these nonlinear QNMs, both analytically and numerically.In the following, we will loosely use the term "nonlinear" to explicitly refer to second-order perturbation theory of BH spacetimes, even if in full generality the nonlinearities also contain cubic and higher-order perturbations.
An open problem in the study of nonlinear QNMs, which this article aims to contribute to, is determining the amplitude of these modes.In ringdown models, the amplitude of linear QNMs is typically a free parameter dependent on initial conditions, which should be fitted against data [38].Being generated by nonlinear processes involving the multiplication of linear modes, we expect that the amplitude of nonlinear QNMs can be entirely determined by the amplitudes of the linear modes themselves.Such a statement has already been proved in numerous previous works [27-31, 36, 39, 40].By considering an idealized model of a pure nonlinear QNM that satisfies the same boundary conditions as the linear ones, we will once again demonstrate that this intuition is true.In other words, our work aims to quantitatively answer the question: «Given two linear QNMs with amplitudes A 1 and A 2 , what is the amplitude A NL of the associated nonlinear QNM?».Knowing how to accurately solve this open problem would greatly benefit ringdown modeling, as it would make it possible to account for nonlinearities without introducing redundant free parameters at the level of data analysis.
Our strategy concerning the study of nonlinear QNM is the following.First, we will introduce in section 2 the equation obeyed by second-order perturbations of the metric field around a Schwarzschild black hole, which turns out to be a Regge-Wheeler (RW)/Zerilli equation with a source term containing the nonlinearities [29,30,[32][33][34][41][42][43].When dealing with the nonlinear source term, particular care should be given to the choice of gauge as its asymptotic behavior crucially depends on that choice [42].For simplicity and to focus the discussion on the analytic method we developed to solve the RW/Zerilli equation, we will use a simplified expression for the source term featuring the correct asymptotic behavior.It's important to note that there are no obstruc-tions in using the full source term, as we show in a concrete example in appendix A. In section 5 we next employ the method of uniform expansions introduced in section 3 to obtain an analytical approximation to the solution of this second-order equation.This allows us to write in eq.(5.8) the amplitude of the nonlinear QNM at large distance, which we subsequently evaluate numerically in tables 2 and 3. We finally estimate the accuracy of the method of uniform expansion in section 6, finding that the approximate solutions should be correct at the percent level.
Except for section 3, we will denote all linear quantities with an overbar ( Ψ) in order to distinguish them from nonlinear ones (Ψ).Our convention for the Fourier transform is Newton's constant is G.All plots in this article are displayed in units where GM = 1, where M is the mass of the black hole.
Note added: As we were finalizing the writing of this article, we became aware of a related work [44], similar to ours in spirit.While using a different technique than the one presented in [44], we can confirm their main conclusions concerning the study of nonlinear QNMs.For example, our estimates for the amplitudes of nonlinear QNMs agree in the particular case illustrated in appendix A. However, in addition to presenting new results concerning also the study of linear QNM frequencies, we believe that our contribution offers a more quantitative estimate of the amplitudes of nonlinear QNMs, as we evaluate the accuracy of our approximation, see section 6.As an additional difference, we insist on requiring the appropriate asymptotic behavior for the term sourcing the nonlinearities (the 'source' term) in the second order RW/Zerilli equation (see section 2).This requirement is necessary to ensure that second-order QNMs are well defined.This point appears to be less emphasized in Ref. [44].
Furthermore, another recent work, Ref. [45], employs uniform approximations to study the Poschl-Teller potential, without focusing on the QNMs, but rather studying the accuracy of the method.

Nonlinear perturbations of a Schwarzschild black hole
It is well-known that linear perturbations around a Schwarzschild black hole of mass M obey the Regge-Wheeler [46] (RW) or Zerilli [47] equation, depending on the parity of the mode.Focusing on eigenfunctions at definite frequency ω and angular momentum ( l, m), and denoting by Ψ(r * ) the dimensionless amplitude of the propagating degree of freedom at linear order (either in the even or the odd sector), we have where r * is the tortoise coordinate r * = r + 2GM log r/(2GM ) − 1 (while r is the standard Schwarzschild radial coordinate), and V (r) is the Regge-Wheeler or Zerilli potential, given by (2.2) At second order in perturbation theory, it turns out that the second-order perturbations of the metric can also be encoded into one parity-odd and one parity-even function which we will denote by Ψ (also chosen to be dimensionless in the following) [29,30,[32][33][34][41][42][43].These second-order amplitudes satisfy the same RW/Zerilli equation, albeit with a source term that is proportional to the product of linear modes: 3) The explicit expression of S(r) depends on the choice of gauge and dynamical variable Ψ, see e.g.[29,33] for some specific examples.Generically, the source consists in the product of two firstorder modes Ψ1 and Ψ2 and their derivatives multiplied by functions of r, constructed so as to be quadratic in the linear amplitudes.Furthermore, even for a given parity of the second-order mode Ψ, the source can depend on both even and odd-parity first-order modes.Thus, we will rewrite S as S(r) = Ŝ(r) Ψ1 (r) Ψ2 (r) , where now Ŝ does not scale with the linear amplitudes 1 .It can further be shown that, with an appropriate gauge choice, one can impose the following asymptotic behavior of Ŝ [29,42]: These asymptotics are crucial for the study of nonlinear quasi-normal modes.This requirement is not specifically linked to our method; any study of nonlinear QNMs should employ dynamic variables Ψ1,2 defined so that eq.(2.5) are ensured, otherwise second-order QNMs would display a divergent power-law scaling, as discussed in Ref. [31].We will later observe that this translates in our formalism into the convergence of the integrals described in appendix B.
In order to keep the discussion as simple as possible, in the spirit of Refs.[31,36] we will use the following toy-model for Ŝ: which respects the boundary conditions (2.5).This choice is mainly motivated by two considerations: 1. Our focus in this article is demonstrating the effectiveness of an analytic method for computing quasi-normal mode amplitudes at the second order, rather than reconstructing a complete waveform template accounting for second-order QNMs.As a result, we do not (at this point) aim at a maximally realistic source term.When using our results to estimate the quantitative impact of second-order modes on ringdown waveforms, it will be important to go beyond the toy-model approximation (eq.(2.6)) and use the complete expression of S in a gauge adapted for gravitational radiation.We leave this issue to further work.
2. Expressions for the source S featuring the correct asymptotic behavior are rare in the literature.To our knowledge, ready-to-use expressions for S are only available for specific values of ℓ in Refs.[29,33,42].In addition, the only one that exhibits the asymptotics (2.5) (particularly near the horizon) is the source present in [29], which however focuses only on (ℓ = 2) × (ℓ = 2) → (ℓ = 4).We further show in appendix A that our toy-model (or simple modifications thereof) can quite accurately fit the function Ŝ provided in [29].
Having presented the necessary details of the dynamics we are interested in, we now turn to introduce our method for solving the RW/Zerilli equation.

Uniform expansions
In this section, we introduce the concept of uniform expansions, first set out in [48,49].We refer the reader to [23,24] for more details and applications of the formalism to quantum-mechanical problems.The basic idea is to compare the Zerilli/RW equation (2.3) with a simpler, auxiliary differential equation which is exactly solvable.While this step sounds similar to what we would do with a Poschl-Teller potential, we will then relate the exact solutions of the new auxiliary differential equation to approximate solutions of the original one.We will freely switch between r and r * inside the argument of functions, leaving the conversion r * (r) implicit.In this section we will not distinguish between Ψ and Ψ, since the discussion applies to both.
We choose to consider the following auxiliary differential equation for ϕ(σ), where σ(r) will be interpreted as a local rescaling of r where σ(r) and Λ(σ) are two functions we will solve for later on.We now want to choose the potential Θ to be "close" in some sense to ω 2 − V (r) in Eq. (2.3).We know that physically quasinormal modes are generated near the maximum of the potential present in Eq. ( 2.3), so that this region should matter the most.In order to exploit this, we choose Θ to be a quadratic function of σ where we have introduced a parameter ν and some convenient normalizations.The two solutions of the homogeneous part of equation (3.1) are where D ν is the parabolic cylinder function.We now relate the exact solutions ϕ A , ϕ B (σ) to approximate solutions of equation (2.3), thanks to a local rescaling that will transform r to σ.The function σ(r) will now be determined by simply substituing the following ansatz into Eq.(2.3).f (r) will later be chosen in a convenient way.Plugging into equation (2.3) and using the auxiliary differential equation (3.1) on ϕ ′′ we obtain where a prime denotes differentiation with respect to the tortoise coordinate r * .We now choose the following expressions for f and Λ: This choice makes everything outside the first parenthesis vanish.We are thus left with a differential equation on the only remaining undetermined variable σ, because ϕ factors out: The last term is called the schwarzian derivative (up to a normalization).If Θ(σ) has been chosen appropriately, σ will be a slowly varying function of r * and the schwarzian should be negligible.It is worth pausing to mention that the usual WKB method would have been recovered had we chosen Θ(σ) ≡ ±1 [23]; the schwarzian is the parameter controlling the validity of the approximation for that case also.This should reassure the reader that, despite our crude auxiliary potential Θ(σ), great accuracy can be achieved with this method2 .To achieve even better results, one could perturbatively include the effects of the schwarzian derivative and Borel-Padé resum the resulting series as suggested in [22].
Assuming for now that the schwarzian is negligible, we can approximate where the sign choice matches σ to r * preserving orientation.Notice that this is the only approximation we have made so far.It is possible to go beyond this approximation and include higher-order corrections, see e.g.[50][51][52], although we will not do it here.Eq. (3.8) completely determines σ up to an additive shift.
WKB approximation is invalid near to the turning points, that is the zeros of ω 2 − V .This is readily seen in this framework, because the schwarzian derivative could blow up.To ensure it remains finite in the whole domain of interest (i.e. the real r * line), turning points should be matched so that dσ dr * remains finite; the schwarzian derivative is then not necessarily divergent (necessary condition for the approximation to work).This will impose a relation between ω and ν.The (uniform) smallness of the schwarzian compared to ω 2 − V (sufficient condition) is to be assessed by hand (see section 6 and also [23] for a discussion).
Let us denote by σ − and σ + the two zeros of Θ, and by r − * and r + * the equivalent zeros of ω 2 −V i.e. its turning points.First, we fix the shift ambiguity of σ by setting σ(r + * ) = σ + .This choice has both the numerator and the denominator inside the square root of Eq. (3.8) vanish linearly, keeping the ratio finite.Then integrating equation (3.8) gives where we have used the explicit form of Θ in Eq. (3.2).This implicitly defines the coordinate σ.In addition, for the schwarzian in Eq. (3.7) to remain bounded, by integrating Eq. (3.8) between the turning points we get an "area law" we get a constraint for the integral of the potential near the turning points: This condition will be used to obtain the value of ω when solving for linear quasi-normal modes (where we will impose a quantization condition on ν).Alternatively, it will give the value of ν for nonlinear quasi-normal modes, where this time we will start from a given value of ω.
To recap, we now have found the generic solution of Eq. (2.3): where ϕ is a solution to the simpler equation (3.1), Θ is given in Eq. (3.2) and σ is determined from Eq (3.9).We highlight that, contrary to WKB approximate solutions, Eq. (3.12) exhibits no divergence close to the turning points we matched, making our expression ready to plot (see fig. 1).We now have to impose the correct boundary conditions to this solution.Let us discuss separately the cases of linear and nonlinear modes.

Linear modes
We study linear modes of angular number l which will be different from the nonlinear angular number ℓ.The linear mode amplitude Ψ obeys Eq. ( 2.3) where the source term S is set to zero.The auxiliary field φ obeys the homogeneous part of equation (3.1), whose solutions are given in Eq. (3.3).From Eq. (3.12) this means that where A and B are two amplitudes that are determined by initial conditions and boundary conditions at the horizon and infinity, and we have introduced an overbar on all quantities related to firstorder modes.Notice that, in order to make the field Ψ dimensionless, the amplitudes should scale as A, B ∝ ω1/2 .Additionally, the fundamental assumption that Ψ is a small perturbation of a background Schwarzschild spacetime translates into the fact that its amplitude is small, i.e.
This assumption will ensure that second-order modes are always smaller than first-order (linear) ones, as they should be proportional to the square of the linear amplitude, see eq. (2.4).Now, it remains to impose the correct boundary conditions.To ensure that waves leave the system at infinity and enter the black hole close to its horizon, we must impose the following QNM boundary conditions: where the choice of sign inside of the exponential is dictated from our convention for the Fourier transform (1.1), and Ψ∞ and ΨH are two constants.The following subsection works out the constraints that Eq. (4.2) imposes on A, B.

Quantization condition
We now determine the quantization condition that will give the quasi-normal modes spectrum.At the same time, we compute the asymptotic behavior of the wavefunction at large |r * | accurately to O(1) phases, because it will play a role in section 5.The hasty reader that is only interested in the quantization condition should focus on equations (4.4,4.5, 4.10-4.13)and equations (4.2, 4.3), jumping straight to eq. (4.16).
In order to determine the behavior of Ψ in Eq. (4.1) we have to relate the asymptotic behavior of σ and r * .By noticing that V vanishes and its integral is finite for both large positive and negative r * , we get from Eq. (3.9) the following asymptotic behavior: We will need in the following the expression of σ2 up to O(1).Evaluating the expansion of Eq. (3.9) to the next order, we get where C∞ and CH are two constants given by ) (4.9) Finally, there remains to obtain the behavior of the parabolic cylinder functions for large σ, which are σ ∼ e iπ ν/4 σν e −iσ 2 /4 , (4.10) for σ → ∞, and We see that, in order to impose the QNM boundary conditions eqs.(4.2) and ( 4.3) we have to choose since 1/Γ(−n) = 0.

Discussion
To recap, we have now found the (linear) solution to the homogeneous equation (2.3) which is D n e iπ/4 σ , (4.17   4.17) for l = 2, n = 0 versus a very accurate numerical solution ΨLeaver using Leaver's algorithm [6], as a function of the tortoise coordinate r * and in units in which GM = 1.The solutions are normalized so that Ψ ≃ e −iωr * for r * → ∞.The inset shows the fractional deviation ∆ = |( Ψ − ΨLeaver )/( Ψ + ΨLeaver )| of the uniform approximation solution versus Leaver's wavefunction.Our approximation is accurate both close to the maximum of V and far from it, although the accuracy degrades for large values of r * .
where A is a small but otherwise arbitrary amplitude that can be fixed e.g. by the initial conditions of a ringdown signal.This solution is plotted in fig. 1 for l = 2 and n = 0, showing that as advertised in section 3 the profile for Ψ is accurate both close to the minimum of V and at infinity.The quantization condition giving the value of ωn is obtained from eq. (3.11): This equation can be numerically solved for ω once given the RW/Zerilli potential in eq.(2.2).On the technical side, we have to ensure that in eq.(4.18) the integration path is chosen correctly.This is nontrivial due to the presence of branch cuts and many Riemann sheets.The prescription we use is to start from real values of ω, where turning points r± * are real and the correct path is obviously identified as a straight line.As the imaginary part of ω is increased, we track the turning points as they move away from the real axis.This prescription also correctly identifies which two of the many (complex) turning points of ω2 − V are to be used.Our integral in eq.(4.18) is then completely well defined.We plot the position of the turning points for l = 2 as a function of the overtone number in fig. 2.
Alternatively, one can also recover the usual WKB formula of Schutz and Will [9] by approximating the RW/Zerilli potential with a Taylor expansion up to second order near its maximum situated at r * : The QNM frequencies given by our approximation eq.(4.18), and 6-th order WKB results are compared to the most accurate numerical estimates in table 1.The comparison reveals that 6-th order WKB performs better than our method for small n, albeit at the price of significant complexity: this comes as no surprise since we are only performing a lowest-order computation, and our result could be improved by taking into account the schwarzian in eq.(3.7).Both approaches show improved accuracy as l increases.A remarkable feature of our approach is that, while WKB's agreement with the numerical solution worsens as n increases, our approximation remains accurate independently of n; this is clearly showcased by looking at the (very extreme) l = 2, n = 1000 mode.
We will now explain this behavior.Employing higher and higher order matching polynomials when using WKB has two positive, but distinct, effects 1.It improves the approximation of V that one uses very close to the top of the potential, which is crucial to obtain very precise results.

2.
V is well approximated in a larger and larger domain, meaning that even for turning points which are quite separated (this happens for large n), the approximate potential is accurate.
We already showed in eq.(4.19) that our method reconstructs a quadratic WKB matching.This is not surprising, since we used a quadratic auxiliary potential Θ(σ).The improvement that our method brings is that, thanks to the local deformations that σ(r) accounts for, V is now well approximated by Θ even far from the maximum.Then, in some heuristic sense, our approximation works at first order in the first effect, but at all orders in the second one.Since the second effect determines how accuracy behaves with n, while the first one has more to do with "absolute" accuracy, the results in table 1 can be understood.These observations also suggest that, if one were to increase the order of the polynomial Θ(σ) beyond σ 2 , one would consistently do better than WKB at the same order.A numerical fit of the relative errors of uniform approximations and 6-th order WKB using the data in table 1 indicates that uniform approximations should outperform 6-th order WKB roughly for n(ℓ) ≳ 4 + ℓ/3 .(4.20) Note that we have used the Zerilli potential for our numerical estimates, however since the Regge-Wheeler and Zerilli potential are isospectral [2] we could have as well used the RW potential.We preferred the Zerilli potential because for the RW potential in the complex r * plane, the relevant turning points of ω2 − V RW approach the same limiting r * point for large n, making the integral more susceptible to numerical instabilities.As can be seen from figure 2, this is not the case for the Zerilli potential.
We now briefly explain how we computed the first column of table 1.
Efficient computation of ω Inverting equation (4.18) to find ω for given n may seem like a daunting task.While (unless approached in special limits) the integral is not doable analytically, we mention how we obtained the frequencies without a painful scanning in ω ∈ C.
If we have a target ν ≡ n ∈ N and a rough guess ω0 of the related frequency, we can do the integral and compute its ν0 .Then the following holds (up to subleading corrections in δ ω) where all integrals are all convergent and between the same turning points r± *3 .We are then led to so that our guess ω0 can be improved to ω0 + δ ω if δν is small.To obtain a guess, we used the frequency of modes with slightly smaller n or l; we then run the above algorithm, observing very fast convergence (1-2 iterations are sufficient).For n = 1000 in table 1, we used the numerical result as seed.A prediction independent of this input would require a more careful study of the convergence properties of our algorithm.

Nonlinear modes
Having proven the validity of our method to compute the linear spectrum of QNM, we now use the technique to compute the amplitude of second-order modes using the nonlinear RW/Zerilli equation described in section 2. We first explain how to compute the nonlinear amplitude using uniform expansions before discussing the implications of our results.

Ratio of amplitudes
We now consider the differential equation (2.3) for a nonlinear mode, sourced by two linear modes Ψ1 and Ψ2 of frequencies ω1 , ω2 and mode numbers ( l1 , m1 , n 1 ), ( l2 , m2 , n 2 ).We also denote by A 1 and A 2 the amplitudes of the linear modes entering eq.(4.17).Since the source term is quadratic in the linear modes we know that the frequency ω of the nonlinear mode is given by [31] where the star denotes complex conjugation.In the following, we will focus on the first case ω = ω1 + ω2 , the treatment of the second one being very similar4 .On the other hand, a whole range of values of ℓ are allowed according to the usual rules of multiplication of spherical harmonics [31][32][33][34]: We want to solve eq.(2.3) using the same method of uniform approximation as before, by comparing it with the simpler eq.(3.1).Once ω and ℓ are given, we can compute the value of ν for nonlinear modes following eq.(3.11).We give in The quasi-normal frequencies of a Schwarzschild black hole in units where GM = 1, comparing our uniform approximation method, 6 th order WKB and numerical results.In parenthesis the relative error on the absolute value with respect to the numerical result | ωexact−ωapprox.ωexact |.Note that, while we present truncated values of ω, we used more accurate values to compute errors.The data in the second and third column of table 1 (for l ≤ 4) can be found for example in [11].We also computed the exact frequencies using the convenient implementation of Leaver's method that is found in the Black Hole Perturbation Toolkit [53].Similarly, [54] provides code to compute the WKB prediction.Lastly, numerical frequencies for l = 2, n ≤ 1000 can be found in [55].values of l1 , l2 and ℓ.Then, using the method of variation of constants we first find the solution to the differential equation on ϕ, eq.(3.1): where ϕ A/B are given in eq. ( 3.3), and c A/B are two constants that we will tune in order to ensure QNM boundary conditions.Notice that we have conveniently normalized ϕ by the product of amplitudes A 1 A 2 to which it should be proportional.In appendix B we show that the two integrals defining F A and F B are convergent for r * → ±∞ 5 .As we anticipated in section 2, this is intimately related to the asymptotic features of the source term Ŝ, and the integrals may have failed to converge had we worked in another gauge where Ŝ did not display these asymptotics.Thus, it is now trivial to read the asymptotic behavior of the nonlinear mode Ψ given in eq.(3.12), as it just parallels the case of linear modes: where C ∞ and C H can be obtained from the linear expressions in eqs.(4.8) and (4.9) just by replacing ω → ω, ν → ν, l → ℓ and r+ * → r + * .QNM boundary conditions impose Finally, we can get the amplitude of the nonlinear mode normalized by the product of linear modes at infinity, which is the important quantity for ringdown models in gravitational-wave observations: From this equation one can numerically compute this ratio of amplitudes.Indeed, notice that the source term Λ is related to the original source S present in the RW/Zerilli equation by eq.(3.6), and S itself depends on r which is related to σ via eq.(3.9).Moreover, notice that the √ ω factor is needed in order to make the ratio of amplitude dimensionless since Λ, A 1 and A 2 are dimensionful.

Discussion
In tables 2 and 3 we give the numerical value of the ratio in eq.(5.8) for the Zerilli equation (i.e. both Ψ, Ψ1 and Ψ2 are of even parity) with a range of values of ℓ, l1 and l2 .Some regularities emerge from the table; for example, notice that while the real part of ν is nearly constant, its imaginary part is approximately equal to l1 + l2 − ℓ.This behavior could be confirmed e.g. by looking at the eikonal limit of our computations, which we plan to investigate in a near future.A clear trend that we can deduce from table 3 is that the ratio (5.8) is maximized when the Clebsch-Gordan upper bound is saturated, i.e. ℓ = l1 + l2 6 .Since the damping times of the nonlinear modes, equal to the inverse of the imaginary part of ω, are all of the same order in table 2, this means that for fixed linear amplitudes the nonlinear modes with ℓ = l1 + l2 should be the dominant ones.This is a robust conclusion of our toy-model; of course, it remains to be elucidated whether this feature persists in a more realistic ringdown model of nonlinear modes, where one would have to take into account the exact expression of the source term (instead of the simple model eq.(2.6)) and Clebsch-Gordan coefficients.Moreover, the actual amplitude of the nonlinear modes will also depend on the product of amplitudes of linear modes A 1 × A 2 , which tends to be suppressed as one moves away from the dominant mode l1 , l2 = 2 produced in black-hole mergers, at least for approximately equal mass, quasi-circular binaries [56][57][58][59][60].
Nonetheless, our results show that it should be technically possible to deduce the amplitude of nonlinear modes in a ringdown signal solely from the measurement of the linear amplitudes A 1 , A 2 .As such, we can envision two interesting applications of our method, which we plan to explore in the near future: • Improvement of ringdown models: Typical ringdown models depend on the N complex amplitudes of the N linear modes that one includes in the model [61].Within the regime of perturbation theory, it should be possible to add to the signal, on top of the linear modes, the associated nonlinear modes without introducing any free parameter.Our approach, were the correct source S in eq. ( 2.3) known, would enable one to compute the amplitude of the nonlinear modes and to improve ringdown models without any additional cost from a data analysis perspective.
• Tests of General Relativity: Alternatively, one can view our results as hinting properties of the nonlinear-to-linear amplitude ratio in GR.For instance, it should be possible to include additional free amplitudes in a ringdown model at the frequencies of nonlinear modes, in the spirit of [27,28], and compare the value obtained from data to the GR prediction.This would constitute a new test probing GR further into its non-linear regime.

Accuracy of uniform expansions
All our previous results have been obtained using the method of uniform expansions which is itself an approximation.It is thus quite natural to ask what is the accuracy of this approximation.Going back to section 3, we see that the only approximation we have made so far is to neglect the second term on the right-hand side of eq.(3.7) in order to be able to compute σ ′ .This is valid when Once we have obtained the numerical solution for σ(r) from eq. (3.9), it is easy to check this inequality.In fig. 3 we show the magnitude of the ratio in eq.(6.1) both for the l = 2, n = 0 linear solution and the ℓ = 4, l1 = l2 = 2, n 1 = n 2 = 0 nonlinear solution, as a function of r * .In the linear case, the ratio (6.1) reaches a maximal value of ≃ 0.08.Notice, however, that it does not translate directly into a 8% error in the estimation of the QNM frequencies, as from table 1 the error of our method, compared to a more accurate numerical result, is only of 3% for the linear l = 2, n = 0 QNM frequency.In the nonlinear case, the approximation seems to work much better, as the maximum of the ratio (6.1) is only of 0.02.Since we cannot yet precisely compare our results to the numerical relativity ones in e.g.[27,28] 7 , it remains quite difficult to estimate the actual error in the ratio of nonlinear-to-linear amplitudes, eq.(5.8), but we can give an upper bound of 2% in this particular case.This accuracy could then be improved by computing next-to-leading order terms in the uniform expansion as is done in e.g.[50][51][52].Whether or not these higher-order corrections would be needed in ringdown models given the sensitivity of gravitational-wave detectors such as LISA [62] is an interesting question that we leave to further work.

Conclusions
In this work we have discussed the application of the so-called method of uniform approximation to black hole perturbation theory.As we have shown, this method provides reliable predictions of both the linear quasi-normal frequencies (table 1), and the amplitude of the nonlinear modes (tables 2 and 3).Concerning the linear spectrum, one important advantage of this technique is providing accurate formulas for the QNM frequencies even at large overtone number n, while the WKB   2, but for ℓ = 4, 5.We have highlighted in blue the values of the ratio whose order-of-magnitude is 10 −1 , which precisely correspond to the values for which ℓ = ℓ 1 + ℓ 2 .
At the moment we do not have any physical explanation on the sign difference of the blue ratios, which seems to be related to a slower decay of the integrand when ℓ = ℓ 1 + ℓ 2 .
approach loses accuracy as n increases.Other advantages of our formalism compared to WKB are that it does not require any matching between approximate solutions in different regions of the integration domain, and its relative simplicity and compact formulas.Contrary to what was done in [22] where V (r) was then expanded around the minimum, we preferred to numerically evaluate all integrals to maximize precision.
We then investigated the predictions of our method for the amplitude of second-order perturbations sourced by linear QNMs.For simplicity, we assumed a toy model source eq.(2.6) for the RW/Zerilli equation at second order in the main body of our work.As we show in appendix A for the case of ℓ = 4, l1 = l2 = 2, the same method can be straightforwardly applied to the full GR The ratio is always small, signalling that the error on QNMs frequencies and nonlinear amplitudes is at the percent level.
source, when the dynamical variables are chosen in such a way that this source displays the appropriate falloff conditions.In the context of this particular toy-model, our numerical results suggest that the amplitude of nonlinear QNMs is maximized when the right Clebsch-Gordan inequality is saturated, ℓ = l1 + l2 .A realistic source will have a more complex dependence on ℓ, ℓ 1 , ℓ 2 , making it difficult to establish at this point whether the same happens in GR or not.We evaluated the accuracy of our technique in section 6, finding an error of at most the percent level due to the uniform approximation, with respect to an exact solution of the quadratic problem with the toy model source.
The wavefunction can also be approximately computed, both numerically at finite r * and analytically in the limit of r * → ±∞.We believe one of the main strengths of the method to be the analytic control that it provides on the problem, a feature that could allow one to analytically study eq.(4.18) and eq.(5.8) for excited overtones n → ∞ or in the eikonal limit ℓ → ∞.
Our results indicate that the amplitude of second-order QNMs, properly normalized by the product of amplitudes of the linear modes which generate it, can be computed in General Relativity.Thus, it should be possible to use this fact in order to improve ringdown modelling or design tests of GR, as we discussed in section 5.2.For this we would need to derive the expression of the source term in the RW/Zerilli-type equation in a suitable gauge, and then relate the amplitude of the RW/Zerilli-type scalar to the asymptotic waveform, a task which we plan to undertake soon.Moving to more and more realistic scenarios, it would be interesting to employ our method in the case of Kerr spacetime.Even more ambitiously, it would be interesting to model the source terms in cases of departure from GR and standard BH geometry, e.g. using the formalism in [12] and leveraging as much as possible the properties of the near-light ring geometry.treatment of the limit r * → −∞ being similar.When sending r * to real infinity, the coordinate σ does instead go to a complex infinity following σ ≃ 2(ωr * ) 1/2 .Because of the imaginary part of ω, the integral defining F A/B is in fact exponentially divergent.Let us now show that this divergence is immaterial and that we could as well assume that σ → ∞.
Let us introduce a real coordinate σ R "close" to σ at large r * .We choose σ R = 2(|ω|r * ) 1/2 for definiteness, but we could have worked with other choices as well.We now split the integral defining F A/B in two pieces, σ σ + = σ R σ + + σ σ R (assuming analyticity of the integrand in the region of interest).Our aim is to first show that the second piece J A/B = σ σ R gives a subleading contribution to ϕ(σ), eq. ( 5.3), in the limit r * → ∞.Noting that both σ R and σ are large, we can take the large-r * limit in the integrand to find: where we have changed the variable of integration from σ to r * , used the asymptotic limit Ŝ ≃ r −2 * for r * → ∞ and from now on in this appendix we will not compute the exact value of the constants in front of the dominant term of the asymptotic expansions.The integral can be performed exactly.
The lower bound vanishes for large r * , while the upper bound still diverges exponentially because of the imaginary part of ω.Similarly, we find This shows that, at large distances, the contribution to Ψ of the second integral, from σ R to σ, is (from eq.(3.12)) O r −2 * × e −iωr * , i.e. an outgoing wave suppressed by the falloff of the source term.This is asymptotically negligible compared to the standard outgoing wave e −iωr * and does not change the ratio of amplitudes at large distances.The reasoning is exactly the same for the limit r * → −∞.Notice that the falloff conditions assumed by the source term are quite crucial in order for these terms not to contribute.
Let us now show that the first piece of the integral defining F A/B , eq. (5.4) with σ R as upper bound, is convergent for σ R → ±∞.To avoid clutter, we will from now on remove the R superscript on σ, keeping in mind that we are sending σ to real infinity.Let us first evaluate the convergence as σ → ∞.We have to relate the asymptotic behavior of the background coordinates σ1 and σ2 to σ; this can be easily obtained from eq. (4.6).Using eqs.(4.10) and (4.17 /2 e −σ 2 /(8GM ω) + Cst × σe −σ 2 /(8GM ω) , (B.9)

Figure 1 :
Figure1: Wavefunction in the uniform approximation: Plot of the real and imaginary parts of the wavefunction defined in eq.(4.17) for l = 2, n = 0 versus a very accurate numerical solution ΨLeaver using Leaver's algorithm[6], as a function of the tortoise coordinate r * and in units in

Figure 2 :
Figure 2: Turning points (in pairs) r±* for the Zerilli potential with l = 2 as a function of the overtone number n in the complex plane, in units where GM = 1.When n increases, the turning points approach the two poles of the Zerilli potential (cyan * ).This feature allows the integral in eq.(4.18) to blow up as n → ∞.

Figure 3 :
Figure3: Accuracy of the uniform approximation: We plot the ratio (6.1) both for the l = 2, n = 0 linear solution and the ℓ = 4, l1 = l2 = 2, n 1 = n 2 = 0 nonlinear solution, as a function of r * in units where GM = 1.The ratio is always small, signalling that the error on QNMs frequencies and nonlinear amplitudes is at the percent level.

2 )
This time the second integral completely vanishes at both ends.Plugging this result into eq.(5.3) defining ϕ, we getJ A ϕ B − J B ϕ A ≃ O r

Table 2 a
sample of values of ν for a range of

Table 3 :
Same as table