Non-linear Black Hole Ringdowns: an Analytical Approach

Due to the nature of gravity, non-linear effects are left imprinted in the quasi-normal modes generated in the ringdown phase of the merger of two black holes. We offer an analytical treatment of the quasi-normal modes at second-order in black hole perturbation theory which takes advantage from the fact that the non-linear sources are peaked around the light ring. As a byproduct, we describe why the amplitude of the second-order mode relative to the square of the first-order amplitude depends only weakly on the initial condition of the problem.


Introduction
During the final stage, called ringdown, of the merger of two Black Holes (BHs), Quasi Normal Modes (QNMs) describe the BH response to any kind of disturbance and their frequencies only depends on the BH mass, spin and charge of the final BH [1].As such, QNMs provide a privileged tool to understand the properties of BHs, which is one of the major goals of gravitational wave astronomy [2].
The Gravitational Wave (GW) strain is generically described using first-order BH perturbation theory.However, non-linearities are an intrinsic property of general relativity.As a consequence, second-order effects turn out to be relevant in describing ringdowns from BH mergers and have received much attention recently [3,4,5,6,7,8,9] (see Refs. [10,11,12,13,14,15,16,17] for some older results).For example, from numerical results like those of [3,4], it can be seen that non-linearities could account for a relevant fraction of the amplitude signal, being approximately proportional to 10% or more of the linear amplitude squared.
The GWs produced during the ringdown and sufficiently far from the BH horizon are the sum of exponentially damped QNMs, characterized by two angular harmonic numbers (l, m) and an overtone number n.Their amplitude is denoted A (l,m,n) and their oscillation frequency and decay timescale are given by the real and imaginary parts of ω (l,m,n) , respectively, Email addresses: davide.perrone@unige.ch(A.Riotto), thomas.barreira@polytechnique.edu (A.Riotto), kehagias@central.ntua.gr(A.Riotto), antonio.riotto@unige.ch(A.Riotto) r h l, m (t − r) = n≥0 A l,m,n e −iω l,m,n (t−r) . ( The non-linearities are generically caught by the ratios where the amplitudes A l 1 ,m 1 ,n 1 ×l 2 ,m 2 ,n 2 of the second-order nonlinear GW strain is constructed from the linear GWs with amplitudes A l 1 ,m 1 ,n 1 and A l 2 ,m 2 ,n 2 with frequency (1.3) Some generic outputs have been pointed out in the literature, independently from the properties of the BH or from the type of the BH merger, be it originating from the quasi-circular merger of two BHs giving rise to a Kerr BH or from head-on spinless BH collisions resulting in a spinless BH: 1. the amplitude A l 1 ,m 1 ,n 1 ×l 2 ,m 2 ,n 2 scales like A l 1 ,m 1 ,n 1 A l 2 ,m 2 ,n 2 , such that in the ratio NL the amplitudes cancel out; 2. the relative amplitude of second-order modes depends mildly on the initial conditions [18].
The goal of this paper is to offer an analytical treatment of the second-order QNMs and an explanation of these two general results.We will do so by offering some generic considerations about BH perturbation theory at second-order and showing that, upon mild assumptions and based on the nature of the non-linear sources, properties 1 and 2 can be recovered.We will also offer more quantitative results by using the WKB approximation.
The paper is organized as follows.In section 2 we offer generic considerations and focus on the first-order solution.Section 3 contains the description of the second-order evolution.In section 4 we devote our attention on the dependence on the initial conditions, while the WKB approximation is used in section 5. Section 6 provides various examples, and section 7 contains our conclusions.Finally, the paper is supplemented by various appendices.We will use units of G = ℏ = 1

Setting the stage
We start by writing the BH equations for the linear and quadratic perturbations, ψ 1 and χ 2 respectively, in their generic and somewhat symbolic form [2] where the subscript indicates the order in perturbation theory and where the tortoise coordinate is defined as For practical purpouses we will use in the following the Zerilli potential V(r(x)) where l is the angular momentum number, but our considerations are not limited by this choice.We have introduced a generic source term S , which is a quadratic function of the firstorder solution ψ 1 (x, t) and its derivatives in time and in space.
To fully solve the problem we need to also specify two initial conditions for the first-and second-order solutions ) and a set of boundary conditions, which will depend on the problem we would like to solve.In our case we are interested in QNMs, and we will take purely ingoing boundary conditions at the horizon r = 2M and purely outgoing at infinity r → ∞.Schematically, using the tortoise coordinate and going in frequency space ψ(x, ω) → e iωx , x → ∞, (2.8) (2.9) We will now proceed by the Laplace technique whose main properties we will briefly summarise for the reader's sake.

Laplace transform
We define a Laplace transform operator as and its inverse as where the integral is performed in the complex s plane and ϵ → 0 + .Usually to perform the integral one closes it with a semi-circular path sent to infinity, but one can also deform it to take into account the presence of branch cuts on the negative real axis.Furthermore, one can think of the s variable as s ∼ −iω, to connect with the usual Fourier analysis.In this case the causality properties are linked to absence of poles in the half-plane Re(s) > 0. A very useful property of Laplace transform is which will allow us to introduce the initial conditions in the differential equation as an extra source term.

The first-order solution
We start by focusing on the first-order problem.Most of the material in this section is standard and the expert reader can skip it.
We first perform the Laplace transform in the variable t, to get where the initial conditions enter in the equation as a source term because of Eq. (2.12) [19].We define the source term for the linear problem as to write the differential equation in the familiar form.
It is also straightforward to write purely ingoing and outgoing boundary conditions in Laplace space, ) (2.17) We can now look for the corresponding Green's function To find it in Laplace space we solve first the homogeneous problem where ϕ + is the solution of the homogeneous equation with only the boundary condition at x → +∞ and ϕ − is the one with the boundary condition at x → −∞.This means that the asymptotic behaviour of the solutions is The Green's function can be written as where and we have introduced the Wronskian of the solution, as a function of s Notice that in this type of differential equations the Wronskian is independent of the specific x chosen, and it will have a zero in s whenever the two solutions ϕ ± coincide, satisfying both boundary conditions at the same time.We proceed by performing a convolution of the Green's function with the source term, containing our initial conditions.This will give us the solution of the differential equation in Laplace space or, explicitly, (2.24) Using the inverse Laplace transform we have the solution at first-order

.25)
As we integrate in the complex s plane we will pick up poles of the function ψ 1 (x, s) and eventual contour deformations if branch cuts are present.If we assume that the numerator has no poles in s, we are left with the zeros of the Wronskian W(s) in the denominator, which give us the damped oscillatory behaviour of the QNMs.The branch cuts that could arise from the numerator are linked to the subleading polynomial behaviour which is observed after the ringdown.We will ignore the effects of such polynomial behaviours because our goal is to confront the analytic calculation with a specific ratio of amplitudes extracted from the dominant oscillating behaviour in numerical experiments.
Picking poles of the Wronskian forces ϕ + = ϕ − , and we can write (2.26) and the distinction between x > and x < can be dropped, choosing one of the two ϕ's as function of x and the other as function of x ′ .If we are interested only in the exponential behaviour, we can define ) and build up the solution in a series expansion as (2.28) In the following sections we will further approximate this result by picking only the first zero of the Wronskian s 0 , the fundamental mode, because it is the least damped at infinity (2.29)

Source term
In full generality, the source is the sum of terms of the type (2.30) where D[t, x] is a differential operator applied to the first-order solution, and R i (r) is a generic ratio of polynomial functions in the r variable.The differential operators D[t, x] can always be put in the form with ( j, k) generic non-negative integers, because we can always exchange ∂ 2 t with (∂ 2 x − V(x)), being ψ 1 a solution for the first-order differential equation.Inserting the explicit form of ψ 1 (x, t) found previously and approximating it further by keeping only the dominant mode Re(s) Im(s) we get where It is interesting to notice that in general we have a factorisation of the constant c 0L which has all the information about the initial conditions and we also have the expected oscillation with frequency 2s 0 ∼ 2ω lm0 .Furthermore we know that the source is peaked around a certain x 0 , which is where the light-ring is placed and approximately at the peak of the potential V(x).

The evolution at second-order
Having set the stage, we are now ready to solve the secondorder problem We start again with the Laplace transform in t, which gives where now the source J 2 is and the Laplace transform of the source is simply Since the Green's function is the same as at first-order, the solution at second-order is just which explicitly reads Performing the inverse Laplace transform we can identify two types of contribution, depending on the pole structure where the first is the term obtained from the initial conditions and the second comes from the source (3.9) Let us look at them separately, starting with the contribution from initial conditions.It has exactly the same analytic structure in s as the one obtained at linear order, so it will contribute to the solution with frequency s 0 , which will of course depend upon the angular momentum we are considering.Exactly as before we can integrate over the poles of 1/W, sending ϕ + = ϕ − , with the same analyticity assumptions on the functions we made for the linear case.Consequently, (3.12) As for the source term, we have two contributions for the poles, the zeros of the Wronskian and the pole of the source.Separating them we get and At this stage, it is crucial that the spectrum s n coming from the Wronskian have no pole corresponding to 2s 0 , otherwise we shall end up with a double pole.This situation does not happen for the QNM spectrum, at least for the dominant frequencies.
For the Wronskian poles part we can define and we can pick again only the dominant oscillatory contribution to write We stress again that in general the s 0 found at second-order could be different from the s 0 found at first-order, because it could be characterized by an angular momentum l different from the one at first-order.This is a crucial point and we will come back to it in the next section.More importantly, the solution at second-order has a contribution which oscillates as 2s 0 , being s 0 the linear frequency.and m = 0, with head-on initial conditions, as in Ref. [11].For the plot we used units of M = 1.
The last step is to take now the limit x → ∞ for the source term at 2s 0 .Here comes our basic assumption.If the source is localised in a region x ′ L < x ′ < x ′ R , we can assume that we are always in the case x > x ′ .This assumption is justified because we know that, for the QNMs, the light ring region dominates in the source at the initial time (aftwerwards, the source propagates), see in Fig. 2, and the observers are far away from it.
Hence we can write in this limit (an alternative derivation is offered in Appendix C) which, defining can be rewritten as Recalling the shape of ϕ − , it might seem that this integral is divergent as it grows exponentially at infinity, a typical behaviour when studying QNMs only on the x variable, as also pointed out in [10].The convergent behaviour is ensured by causality, i.e. reintroducing the time variable t and a theta function θ(t − x).
We have not introduced it for the sake of clarity of our equations, and because in our estimate we will take only the peak of the source, assuming its convergence.We get a form which resembles the linear case one.It is this property which will allow us to considerably simplify the non-linearities of the QNMs.

The dependence on the initial conditions
We are now in the position to offer some considerations about the dependence on the initial conditions.As we have seen, the second-order solution has various contributions.
1.The first one χ I 2 comes from the solution of the homogeneous equation.Depending upon the angular momentum we are interested in, this can renormalize or not the firstorder solution, i.e. in the specific non-linearity ratio, the linear frequency part at the denominator will receive contributions also from the non-linear part.For instance, one can consider an l = 4 source term generated from the linear solution at l = m = 2.This was the case considered in Refs.[4,3].In such a case, the Zerilli potentials are different at first-and second-order and therefore the solutions of the homogeneous equations are different and the linear solution is not renormalized.2. The second contribution χ S 2 W poles comes from the source term picking up only the poles of the Wronskian and again the corresponding frequency can be equal or different from the linear one, thus it will renormalize or not the linear solution.

Finally, the third contribution χ S
Let us consider the case in which the source is constructed from an angular momentum which is different from the linear solution, like in Refs.[4,3].If so, extracting from the second-order solution only the third contribution and comparing to the firstorder solution at null infinity, we get where we used the asymptotic relation which shows that the ratio does not depend on the initial condition of the problem as c S scales like c 2 0L , see Eq. (3.18).This explains the mild dependence on the initial conditions found in Ref. [18].Notice also that at null infinity ϕ + (x, 2s 0 ) ∼ ϕ 2 − (x, s 0 ) ∼ e −2s 0 x and the dependence on x correctly cancels out.
Let us summarize here the basic results which might a bit obscured by the many expressions we have written so far: 1.If the second-order source is well localized and the homogeneous solution of the second-order mode has a frequency which is different from the frequency of the linear mode, then defining NL by extracting from the secondorder result the piece with such double frequency, the ratio NL depends mildly on the initial conditions of the problem at null infinity.Furthermore, the amplitude of the secondorder contribution scales indeed as the square of the amplitude of the linear one.2. Still under the assumption of a well localized source, if the homogeneous solution of the second-order mode has a frequency which is equal to the frequency of the linear mode, then defining NL by extracting from the second-order result the piece with such double frequency, the ratio NL depends on the initial conditions of the problem.Such a dependence becomes mild when the renormalization of the linear piece from the second-order homogeneous solution is negligible.3. The amount of non-linearity will depend on some intrinsic properties of the final black hole state, such the mass and the spin.They will show up in the second-order source and the non-linearity will automatically depend on it.

Evaluating the non-linearity from the WKB approximation
In order to estimate the level of second-order non-linearities, we can use the WKB approximation which we summarize in Appendix D. Our work differs from Ref. [5] mostly because we are not using a toy model source, but a realistic one to give the estimate.Furthermore we are interested in explicitly evaluate the nonlinearity predicted in Ref. [3].We start from the crucial term (3.18), which we rewrite as where we have introduced the angular momentum l to keep track of the different solutions and sources From the WKB analysis of the linear solution we can provide an estimate for h l (x).The source is peaked close the potential peak, and is therefore a well localized function in the region I 2 of the WKB, defined in Appendix D, where ϕ − can be written as an exponential, as in Eq. (D.40).We focus on the case for which the linear solution is not renormalized.The procedure can be also followed in the case in which the linear solution is renormalized, and the result have a similar form.We get (assuming a mild dependence of x 0 on the multipoles) where k l 1 is defined in Eq. (D.13) and A − (l) is evaluated for the specific l and for s (l, 0) , which in this case is s 0 , found at linear level.The region I 2 in the WKB approximation is identified as the one between the two inversion points, close to the potential peak.Those regions I 1,2,3 are fully defined in the Appendix D. We can estimate the integral using the method of the steepest descent close to the point x = x 0 .Being ϕ (l 1 ) − an exponential, we can factor out the function and write h (l) (x) as where we defined for simplicity the rational function and P n i (x) and P m i (x) are polynomials of degree n i and m i .Substituting in the integral we get and the steepest descent method can be immediately evaluated by taking x ′ = x 0 inside the polynomial terms and computing the resulting phase The non-linearity parameter NL is then extracted from (5.9) Using the approximate form for the Wronskian in WKB and the solution ϕ − extended to +∞, defined in Appendix D we get where ν satisfies [20] We can notice that most of the coefficients connecting the various pieces of the solution simplify in the end and we are left with the ratio between the matching of the regions I 2 and I 3 , A + (l 1 ) at linear level and the same matching coefficient but at nonlinear level A + (l).They are all defined in Appendix D.

Examples
We can provide now a few example.We will start with a source term for the second-order with l = 4, generated from a linear solution with l 1 = l 2 = 2.We are therefore in the case in which no renormalization of the linear coefficient takes place and the non-linearity does not depend on the initial conditions.In this case we can immediately apply the formulas derived above, in particular Eq.(5.10), after having evaluated the source term at the peak of the potential.
We use the source obtained in Ref. [13], reported in Appendix A, Eq. (A2), valid for a Schwarzschild BH.We notice that the only dependence on m that we obtain is in the source term, which is generated with a linear term l 1 = l 2 = 2 and If we instead calculate the non-linearity of the GW strain h, accounting for the corresponding normalization ψ = 2 ḧ and accounting for the two polarization, we get, in units of the BH mass M NL h = ω 2 220 NL ≃ 0.069, ω 220 ≃ 0.37, (6.2) which fits well the result in Ref. [18] (see their Fig. 3).
In the case in which the non-linearity is evaluated for l = 2 and l 1 = 2, we follow a similar procedure with respect to the one presented more generally in section 5.The main differences are the solutions ϕ ± , which now are the same as the one found at linear level.So A ± are the same for the two orders and k l 1 = k l 2 = k l .For such a case, we consider the source reported in Appendix A, Eq. (A1), derived in Ref. [11] for the head-on collision of two Schwarzschild BHs resulting in a spinless BH.The two BHs start very close to each other at a distance L and with momentum P, as described in the series of works [21,22,23,11].This behaviour is encoded in the initial conditions, which is Misner's wormhole solution [24], where the two BH solution is constructed as a wormhole solution at fixed time.The corresponding starting conditions for the linear solution are and have a more complicated expression for the non-linear solution [11].Such initial conditions are the evolved using the usual BH theory through the Zerilli potential.The form of the non-linearity in this case can be written as Inserting the various parameters in the formula above we obtain a non-linearity equal to which is much smaller than the one generated at l = 4 found above.We have numerically checked that the dependence on the initial conditions is rather mild, being the renormalization of the linear solution basically absent.

Conclusions
Non-linearities in QNMs are a relevant subject as it will help to improve our understanding of the BH ringdown and probe the non-linear nature of gravity.In this paper we have described a perturbative second-order method to estimate analytically the non-linearity for the QNMs generated at the ringdown of a BH.
The method relies only on the knowledge of the specific source and on the fact that it is well peaked, typically around the light ring, where the potential has a maximum.The amplitude of the second-order contribution scales like the square of the amplitude of the linear contribution and we have described under which conditions the non-linearity depends only mildly on the initial conditions.Our perturbative approach has certainly limitations due to its inability to account for the back-reaction of the GWs onto the BH geometry.
Furthermore, we have analyzed the propagation of the linear and quadratic BH perturbations in the eikonal limit and, by using the WKB approximation, we have obtained solutions to the linear and quadratic perturbation equations.However, the implemented WKB approximation works better in the large multipole limit, while the GW signals are dominated by lowfrequency waves.It could be possible to extend this work further by trying to apply the same techniques to a spinning BH, which could complement the analysis and allow us to compare also with other results present in the literature.

Appendix B. Comments on the non-linear overtones in the source
In this work we have used heavily the approximation that at linear level only the first overtone n = 0 matters, However in the source we have a quadratic dependence on the first-order solution, so the sum squared will contain several overtones.We show here that they do not count in the case of a non-linearity evaluated as in Ref. [4].This is because of the pole structure of the Laplace transformed source.Indeed, we have seen that in general the source will have the form 2) If we plug in the full linear solution ψ 1 we get For the sake of clarity, we relabel and the source term becomes c n c m e t(s n +s m ) h nm (x).(B.5) As we Laplace transform we get simply so that the pole structure will enter in the "renormalisation" of the linear order, but only the selected pole 1/(s − 2s 0 ) will contribute to the dominant frequency non-linearity asymptotically.

Appendix C. Another derivation of Eq. (3.19)
Consider the case in which the first-order solution is constructed from an l = 2 mode and the second-order from an l = 4 mode.The solution at second-order χ l=4 2 (x, t) satisfies where, we have use Eq.(2.33) with S 0 (t) = c 2 0L h(x).It is then easy to verify that the general solution to Eq. (C.1) is where ϕ ±,l=4 (x, s 4,0 ) are the two in independent solutions of the equation and Ψ l=4 2 (x, 2s 2,0 ) is a particular integral of the equation Therefore, the only unknown at second-order is the particular integral Ψ l=4 2 (x, 2s 2,0 ) of Eq. (C.4).To determine the latter, we use the method of the variation of parameters for solving nonhomogeneous differential equations where the particular integral is expressed as a linear combination of the homogeneous problem, but with x-depended coefficients.We write 2 (x, 2s 2,0 ) = u + (x)ϕ +,l=4 (x, 2s 2,0 ) + u − (x)ϕ −,l=4 (x, 2s 2,0 ), (C.5)where ϕ ±,l=4 (x, 2s 2,0 ) solves the homogeneous equation Then imposing the condition we find from Eq. (C.4) Solving Eqs.(C.7) and (C.8) we find that A simple integration then gives where we have used the fact that h(x) is a localized source and we're interested in values of x → ∞.Therefore, we get that the partial integral is given by = cS ϕ +,l=4 (x, 2s 2,0 ), (C.12) which reproduces Eq. (3.19).

Appendix D. The WKB approximation
In this Appendix we follow Ref. [20] using WKB techniques to find the Green's function.We start again from the equation The QNMs can be understood as solutions where the transmission and reflection of the wave have approximately the same amplitude.So we expect to be able to find those solutions by considering "energies" of order of the maximum of the potential.This means that we require, around the maximum x 0 of V(x) and, expanding V(x) we need that This expansion in a quantum mechanical point of view means that we are considering "energies" just below the peak of the potential, so we can identify two inversion points, x 1 , x 2 , located at s 2 + V(x) = 0, (D.5) which will result in two inversion points, located at x 1 and x 2 .We can then identify three regions to rewrite the problem in a WKB fashion In each of the regions certain approximations hold, so we have • region 1: I 1 = (−∞, x 1 ), here the solution is standard WKB up to the inversion point, namely

Transmission and matching
Let us be more precise about the functional forms of ϕ ± and review the matching procedure.The main property that we will use is that the expansion of the parabolic cylinder function solution behaves at ±∞ differently lim x→+∞ ϕ (2) and we are interested only in the term exp(−i , cancelling the other oscillatory behaviour.Let us start from the matching procedure to find ϕ ± in the central region and then we'll match to find the whole solution.At x → −∞ we are in the I 1 region, so we have and to match it with ϕ (2) we need to impose both the conditions to get the correct behaviour.This gives us the definition of the matching coefficient for ϕ (D.26) Instead to match the solution from +∞ (i.e.ϕ + ) is sufficient to impose only B = 0, (D.27) and ν could be generic, hence leading to Now there is a little caveat that should be addressed.From the equations above we have that which means that the inversion points Q(s, x = 0) are complex.Furthermore we need to impose the continuity in a point outside x 1 and x 2 .A possible solution is to take x 1 and x 2 on the real axis where we have the correct oscillatory behaviour from the WKB solution.This is reasonable because it is here that the parabolic cylinder D function has the correct behaviour so it could be matched.The behaviour that we want to match, like in Ref. [20], requires that the quadratic part is bigger than the imaginary part, i.e.
but we also need to stay within the parabolic approximation range, because otherwise we end up in the asymptotic region which evolves like a free wave.An estimate can be done by considering the third derivative and where this becomes important in the expansion, i.e. when which identifies the last point where this matching can be performed, Hence we can choose a point in the range with ϵ ≪ 1.We define x 1 and x 2 as those points respecting this matching property.Now we go on with the matching starting from −∞ and we match ϕ − between the region I 2 and the region I 3 , to get the full solution in the whole space.In I 3 we have where T is the transmission coefficient.We need to match it to

Appendix D.2. The Wronskian
For our evaluation another ingredient is needed, i.e. the Wronskian, built from the functions ϕ ± .As we've seen those are the solutions of the homogeneous equation, but with only one of the two boundary condition satisfied.With those available we can find an approximate form for the Wronskian This function is independent of the specific point picked, but we would like to extract as much information as possible by picking a specific easy point.In particular we would like to choose x 0 , but to be able to do this we need to match the two WKB solutions to the central parabolic region.
We have already seen that the ϕ − solution connected to the region within the inversion points is fully constrained, This is the form of the Wronskian we used to evaluate the nonlinearities in Section 5. We could also expand for ν around n, to show that at linear order we get a zero around s − s 0 .The Wronskian around the zero becomes W(s) = 2πi 2 n+1/2 e iπ/4 Q 1/4 (x 1 , s 0 )Q 1/4 (x 2 , s 0 ) D n (t(x 1 ))D n (t(x 2 )) × × π 2 s 0 Γ(n) (2Q ′′ (x 0 )) 1/4 (s − s 0 ) + O (s − s 0 ) 2 (D.48) where we have defined s 0 as the solution for a specific n and we evaluated ν = n in the A + coefficient.Notice that the same Wronskian but at nonlinear level does not have a zero in s = 2s 0 and in general the ν dependence will lead to both a real and imaginary part.

Figure 2 :
Figure2: Numerical evaluation of the initial source in the specific case of l = 2 and m = 0, with head-on initial conditions, as in Ref.[11].For the plot we used units of M = 1.

∂ 2 xFigure D. 3 :
Figure D.3: Q(x)and its quadratic approximation, The various regions are identified by the solutions of Q(x) = 0, as done in Ref.[20].The plot has units M = 1.

Figure D. 4 :
Figure D.4: Real part of the WKB solution at linear level.Highlighted in red are the two inversion points.Notice that the function seems to diverge at ±∞, but when we reintroduce the time dependence and the causality conditions we end up with an asymptotically vanishing solution.The plot is in units M = 1.