Nonlinear effects in the black hole ringdown: absorption-induced mode excitation

Gravitational-wave observations of black hole ringdowns are commonly used to characterize binary merger remnants and to test general relativity. These analyses assume linear black hole perturbation theory, in particular that the ringdown can be described in terms of quasinormal modes even for times approaching the merger. Here we investigate a nonlinear effect during the ringdown, namely how a mode excited at early times can excite additional modes as it is absorbed by the black hole. This is a third-order secular effect: the change in the black-hole mass causes a shift in the mode spectrum, so that the original mode is projected onto the new ones. Using nonlinear simulations, we study the ringdown of a spherically-symmetric scalar field around an asymptotically anti-de Sitter black hole, and we find that this"absorption-induced mode excitation"(AIME) is the dominant nonlinear effect. We show that this effect takes place well within the nonadiabatic regime, so we can analytically estimate it using a sudden mass-change approximation. Adapting our estimation technique to asymptotically-flat Schwarzschild black holes, we expect AIME to play a role in the analysis and interpretation of current and future gravitational wave observations.


I. INTRODUCTION
Gravitational waves provide an unprecedented ability to witness the approach to equilibrium of strongly-deformed black holes. Immediately following a binary merger, the remnant is in a highly nonequilibrium state, which emits gravitational radiation via characteristic ringing. Observations of these waves by LIGO and Virgo can be used to spectroscopically characterize the final black hole, to test Einstein's theory of gravity in an extremely dynamical context, and to establish consistency with no-hair theorems [1][2][3] and the final state conjecture [4].
To leading order in the deviation from a stationary black hole, the ringdown is mostly 1 characterized by a collection of quasinormal modes (QNMs) [5][6][7]. These have complex frequency, corresponding to oscillations and decay. The QNM spectrum is uniquely determined by the mass and spin of the black hole, which makes ringdown analysis an especially exciting prospect for gravitationalwave astronomy. Indeed, consistency between measured quasinormal frequencies (QNFs) and predicted final states based on inspiral measurements give rise to stringent tests of general relativity [8,9]. For particularly loud events, there is some observational evidence for modes beyond the fundamental (l, m) = (2, 2), consistency of which would test the no-hair theorem [10][11][12].
Although it is not yet possible to routinely observe multiple ringdown modes for individual events, this will change with the enhanced sensitivity of upgraded and future detectors [13][14][15]. On the theoretical side, highaccuracy signal modeling will be critical for controlling * laura.sberna@aei.mpg.de 1 The ringdown can also contain an initial direct component, as well as a late-time power-law tail. We do not consider these components in this work.
systematic errors in parameter estimates and for confidently constraining deviations from GR. To that end, it will be necessary to understand the ringdown beyond linear order-either to incorporate nonlinear corrections into models or to show they are negligible. General relativity is a fundamentally nonlinear theory, and the ringdown is no exception: QNMs can interact. This can be analyzed by comparing linear theory against numerical relativity simulations. Some examples include second-order effects in the ringdown following a head-on black-hole collision [16], in the black-hole response to strong incident waves [17,18], mode-doubling during the ringdown [19,20], and even gravitational-wave turbulence in asymptotically anti-de Sitter (AdS) spacetimes [21,22] (and possibly asymptotically flat spacetimes [23]).
In contrast, for typical mergers, sums consisting purely of (linear) QNMs have been shown to provide a good fit to numerical relativity from surprisingly early times 2 [24][25][26][27][28][29][30] (although this does not automatically imply that the solution itself is linear). It therefore remains to be understood to what degree nonlinear mode interactions are present in astrophysical mergers: it has been proposed that the formation of a common horizon "hides" nonlinearities from asymptotic observers [31], or that the black-hole potential "filters" nonlinear features out of the signal. If true, this hints at a remarkable further property of black holes: not only do they hide singularities, but also nonlinearities! Here we study nonlinear interactions among QNMs and uncover a key effect that causes new modes to be excited. We begin with a toy-model system consisting of an asymptotically AdS black hole and a perturbing scalar field in spherical symmetry. This system allows for precise numerical studies with initial data comprising a single 2 Models based on higher QNFs have been less successful in inferring remnant properties [24,25].
QNM-enabling us to track the excitation of new modes over time. We observe clear QNM excitation during the ringdown that scales nonlinearly with the initial data.
The dominant effect that we uncover is as follows: (1) as the seed mode decays, its absorption causes the blackhole mass to grow; (2) this change to the background spacetime results in a modified QNM spectrum; and finally (3) the remaining seed mode (no longer a QNM of the final black hole) is projected onto the new spectrum, resulting in excitation of additional modes. We call this effect "absorption-induced mode excitation" (AIME). For all cases studied, the QNM decay-and hence the change in the background spacetime-is fast compared to the typical oscillation time scale. AIME can therefore be estimated by taking the approximation of a sudden (nonadiabatic) perturbation to the background, in analogy with quantum mechanics [32]. We show that such analytic estimates are in good agreement with simulations.
Guided by the intuition gained from the toy model, we estimate the magnitude of AIME for gravitational perturbations of (asymptotically flat) Schwarzschild black holes. We estimate the change in mass due to absorption for a typical QNM excited in a merger, and then we project this mode onto the spectrum of the late-time black hole. We find that for realistic amplitudes, this gives percent-level corrections that could be relevant for gravitational-wave observations. In addition to AIME, QNMs can be excited via direct mode coupling [20], but the nature of our simulations isolates the AIME effect. At the practical level, our work provides guidance for future ringdown models and illuminates the transition to the linear regime. It also clarifies the role of nonlinearities and cautions against over-interpreting the success of fits that use a set of linear modes.
Although conceptually clear, the calculations involved in carrying out the mode projection are somewhat nontrivial. In particular, they have required us to generalize (in the Schwarzschild case) the work of Teukolsky and Press [33] to complex frequencies. (To our knowledge this is the first time this has appeared.) We also make use of a recently-presented finite bilinear form on Teukolsky perturbations [34,35]: for QNMs this gives orthogonality, but since the initial perturbation is not a QNM of the final black hole, it can also be used to perform the projection. We provide full details of these calculations. This work is organized as follows. Sections II-IV are devoted to the toy model. In Sec. II we review the theory of scalar QNMs in Schwarzschild-AdS. In Sec. III we describe AIME, and in Sec. IV we describe our numerical studies to elucidate the effect. In Sec. V we analyze AIME for asymptotically-flat spacetimes, making use of the estimates we validate earlier for Schwarzshild-AdS. We summarize our work and discuss implications for ringdown analysis and modeling in Sec. VI. We use units with G = c = 1 throughout.

II. BACKGROUND: QUASINORMAL MODES IN SCHWARZSCHILD-ANTI-DE SITTER
We study Einstein gravity with a negative cosmological constant coupled to a massless complex scalar field. The equations of motion are where G ab is the Einstein tensor, L is the AdS length scale, and T ab is the scalar field stress-energy tensor, We consider perturbations of the Schwarzschild-AdS (SAdS) metric, where We set L = 1 throughout.

A. Scalar quasinormal modes
At leading order, the scalar field evolves in a fixed SAdS background. The response to a small perturbation is characterized by the scalar QNMs. With the ansatz φ = e −iωt R(r)Y m l (θ, ϕ), with Y m l a spherical harmonic, we obtain the radial equation, (6) By redefining the radial function R(r) = r −1 X(r) and using the tortoise coordinate dr * dr = f −1 SAdS (r), the equation is cast into a Schrödinger-like form, where The potential V SAdS asymptotes to zero at the horizon and at spatial infinity. The horizon condition identifying scalar QNMs is the standard ingoing condition, Near spatial infinity, the radial function has two solutions, φ ∼ constant and φ ∼ r −3 . We impose reflective boundary condition at the AdS boundary and therefore choose the solution φ ∼ r −3 . Last, we compute the scalar QNM frequencies and QNM radial profiles numerically using Leaver's method [36,37] as detailed in Ref. [38]. Spherically-symmetric QNFs are simply labeled by the overtone number n.

B. Excitation coefficients
The amplitudes and phases of QNMs are determined from the initial data through the excitation coefficients [39][40][41]. These coefficients naturally arise in the initial data problem for perturbations of a black hole spacetime, in a theoretical framework pioneered in Refs. [36,42,43]. Here we adapt the formalism to the case of a scalar perturbation of SAdS.
Re-introducing the scalar-field time dependence, we want to solve in full generality the equation We define the Laplace transform of the field variable X, and the inverse transform The transform satisfies where the function I(ω, r) is related to the value and derivative of the field at 3 t = 0, φ(0, r) = r −1 X (0) (r) anḋ φ(0, r) = r −1Ẋ (0) (r). By the properties of the derivative of the transform, The general solution to Eq. (13) can be found via the Green's function of the homogeneous equation, defined in terms of two independent solutionsX r+ ,X ∞ , As we see below, it is not necessary to specify the behavior ofX +∞ at the horizon. We can then write the solution to Eq. (13) aŝ 3 The equations are invariant under time translation and we are therefore free to set the initial time to t = 0. ω Figure 1. Schematic integration contour used to evaluate the transform of (17). The integral is originally defined with a contour running above the real line from left to right. The crossed circles mark zeros of the Wronskian (i.e., the QN frequencies), while the semicircle at infinity carries a non-zero contribution associated with the direct propagation of the perturbation to the observer along null geodesics.
The Wronskian of two independent solutions is independent of r, so we can choose to compute it at spatial infinity, To go back to the time domain, we need to plug this solution into the inverse transform (12) and integrate over ω. The integrand has poles in the complex-ω plane at the zeros of the Wronskian, B ∞ (ω) = 0, as depicted in Fig. 1. As seen by the asymptotic behavior in Eq. (15), these are precisely the QNFs ω n = ω n,R + iω n,I , with ω I < 0 and n the overtone number. Thanks to the asymptotic behavior of the potential in the SAdS background, the integrand does not exhibit any branch cuts (present in asymptotically flat spacetimes).
The inverse transform, whose integral runs along and above the real line, is equivalent to the integral along the closed contour drawn in Fig. 1 minus the integral along the semi-circle at infinity. In this work, we focus on the QNM contribution to the response (i.e., the poles), and neglect the contribution of the semicircle at large |ω| (corresponding to the direct propagation of the perturbation to the observer). The latter is absent when the initial data extends to null infinity, as is the case in our simulations. We are left with the integral on the closed contour, which amounts to the sum of the residuals associated with all the poles. Asymptotically, at the AdS boundary, where the QNM excitation coefficients are defined as and where we defined the QN radial profileX ωn = X ∞ (ω n , r * ) =X r+ (ω n , r * )/A ∞ , so that r 2X ωn = 1 at the AdS boundary.
The excitation coefficients are proportional to the overlap integral of the initial condition with the relevant QNM profile. The prefactors, which are independent of the initial condition and only depend on the properties of the background geometry and the type of perturbation considered, are known as the excitation factors, In this work, we are interested in how QNMs interact nonlinearly while being absorbed by the black hole. This is best explored using QNM initial data. For a QNM initial perturbation of frequencyω, the excitation coefficients read The overlap integral is absolutely convergent at the horizon only if the source has compact support. For a source extending to the horizon, as is the case for QNM initial data, a regularization procedure is needed to handle the divergence of the absolute value of the QN profile at the bifurcation surface. To this end, we implement a subtraction procedure (as discussed, e.g., in [34,35] for Kerr), independent of the regularization point R * → −∞.
Thanks to the regularization, we verified that QNMs with different overtone number are orthogonal, C(ω n , ω n ) = 0 for ω n = ω n (see again [34,35] for Kerr). We show the convergence of the overlap integral and the orthogonality of a few modes in Appendix A.

III. ABSORPTION-INDUCED MODE EXCITATIONS
Several studies have identified nonlinear behavior, in particular turbulence, in asymptotically AdS black holes even under weak perturbations [21,22,44]. These studies demonstrated that, for sufficiently large AdS black holes, there can be transfer of energy from short-to long-wavelength metric perturbations ("inverse cascade"), which is not captured by the standard QNM expansion. This phenomenology might have an asymptotially-flat analog, for Kerr black holes near extremality [23].
For generic black holes, the structure of the perturbative equations suggests that first-order modes should drive, at second order, modes with a characteristic combination of harmonics. This was identified using numerical simulations in the gravity sector of asymptotically flat [17] and AdS black holes (in 5 dimensions) [19]. More recently, this second order effect was observed with a perturbative numerical scheme in Kerr black holes [20,45].
At second order in perturbation theory, QNMs can also backreact on the background spacetime, changing the black hole mass and spin. This mediates a third order couplings between the modes. In this work, we describe the absorption of QNMs by the black hole and the third order mode coupling. In this section and the next, we focus on scalar QNMs of AdS black holes. We consider gravitational perturbations of asymptotically flat black holes in Sec. V.

A. Absorption of quasinormal modes
Quasinormal modes have ingoing energy flux through the black hole horizon, which causes the black hole area to increase over time. Waves can be absorbed by the black hole horizon at all stages of a binary merger. Incoming waves, however, need to overcome the potential barrier outside the horizon. During the inspiral, most waves are reflected by the barrier, and the effect of absorption is suppressed [46,47]. Indeed, absorption appears at high orders in post-Newtonian (PN) expansions in the orbital velocity v orb (M f ) 1/3 (M is the total mass, f the gravitational-wave frequency): 4PN in Schwarzschild [46], 2.5PN in Kerr [48]. During the ringdown, however, M f can be larger and the effect of absorption non-negligible. As we will see, black-hole absorption could be a significant effect in the early ringdown, appearing at third order in black hole perturbation theory.
To see this, we first compute the flux of QNMs across the black hole horizon. We write a generic spherically symmetric black hole metric in ingoing Eddington-Finkelstein coordinates, where v = t + r * . Close to the horizon, a scalar QNM can be approximated as where A r+ is the complex horizon amplitude. The timelike Killing vector t a = ∂ a v −f ∂ a r defines a conserved scalar field energy current, where the scalar-field stress-energy is given in Eq. (3). The integrated energy flux across the event horizon is where n a = −t a is the inward normal to the horizon, 4 and M is the black hole mass (technically, the Misner-Sharp mass at the horizon [50,51]). Using the asymptotic limit of the scalar field (24), one finds that for a single mode [52], for any (l, m). For a perturbation containing modes up to overtone N , this generalizes to (28) Perturbations with different (l, m) do not mix at this order in the flux formula, thanks to the orthogonality of spherical harmonics. 5 Equation (27) confirms that QNMs have a positive energy flux through the horizon. Energy falls into the black hole at an exponential (decaying) rate, with a time scale given by half of the mode decay time, This result is valid for scalar QNMs of any sphericallysymmetric static black hole. The analogous calculation of the flux of angular momentum can be found in Appendix D. The energy flux is also proportional to the square of the scalar field amplitude at the horizon, which is itself proportional to the amplitude observed far away from the black hole, although the proportionality constant can be nontrivial. To see how the black hole horizon grows over time as a result of the QNM flux, we write the black hole area A in terms of the irreducible mass, In SAdS, the black-hole mass and the irreducible mass are related by M = M irr (1 + 4M 2 irr L −2 ). The area evolution is then given by where dM irr /dM > 0. This implies where A(∞) > A(0). The total area change from v = 0 is given by The expressions above can also be written in terms of t, the time measured by an observer on the same null slice at the AdS boundary (see Fig. 2).

B. Heuristics of mode excitation via absorption
The evolution of the black-hole background can affect the ringdown in two ways-by altering the QNM spectrum and the mode content of the ringdown (the excitation of modes). We can understand these effects in two extremes: 1. the adiabatic regime, when the background evolves very slowly compared to the typical dynamical time of QNMs already excited, τ back. τ osc. . We associated the dynamical timescale of a QNM with its real frequency, i.e., τ osc. = 2π/ω R ; 2. the nonadiabatic regime (also known as "diabatic"), when the background evolves rapidly compared to the modes, τ back. τ osc. .
This distinction was previously explored in AdS black holes in Ref. [53]. The authors found universal behavior in the black hole response to arbitrarily fast perturbing quenches, i.e., impulsive interactions with a scattering packet of energy, as well as adiabatic behavior in the opposite regime.
In the adiabatic regime, the perturbation remains in the same QN state, and the mode content does not change. This is analogous to adiabatic processes in quantum mechanics [32]. If the initial data is composed of a single seed QNM, the mode evolves into the corresponding mode in the spectrum of the final black hole (with the same angular and overtone numbers). In the end, no other modes of the late-time spectrum are excited.
In the nonadiabatic regime, modes feel a nearinstantaneous shift in the background. In the new background, the seed mode with n =n is no longer a pure QNM, and therefore consists of a collection of modes of the new black hole. This results in the excitation of additional modes of the new spectrum (with n =n), a process we call "absorption-induced mode excitation", or AIME.
To re-expand the old mode in terms of the new ones, we compute the excitation coefficients (21) of the new modes with ω = ω n (M final ), with initial data consisting of the original moden withω = ωn(M init ), The excitation coefficients should approximate the AIME amplitudes. This re-expansion sources all overtones n of the new spectrum, as QNMs of different black holes are not orthogonal. In a realistic scenario, we will have something in between: modes will be excited with a fraction of the nearinstantaneous excitation coefficient (21), and their frequencies will adjust to the background evolution. As we will see in Section IV C, the near-instantaneous excitation coefficients are a good approximation even when the system is not in the fully nonadiabatic regime. AIME is obviously nonlinear in nature: the amplitude of the excitations is cubic in the initial perturbation. We can understand this as follows. In the fully nonadiabatic regime, we need to project the initial perturbation (then mode of a SAdS spacetime with M irr = M irr,init ) onto all the modes of the final spacetime, having mass M irr = M irr,final . From Eq. (32), we expect that at leading order the change in mass will be proportional to the square of the mode amplitude at infinity: The amplitude of each excited mode n will be given by its excitation coefficient (21). If the difference between the initial and final mass is equal to zero, the excitation coefficient of any mode n =n vanishes due to their orthogonality. We therefore expect the excitation coefficient corresponding to n =n to be proportional to ∆M , to leading order in the perturbation amplitude. Taking into account that the excitation coefficient is also proportional to the amplitude of the initial data, we find forn = n. We will confirm this expectation with fully nonlinear calculations in the next section.

IV. NUMERICAL RESULTS
In this section, we numerically evolve a scalar field coupled to gravity to investigate nonlinear effects in the ringdown of a spherically symmetric black hole in asymptotically AdS with reflecting boundary conditions. We do so in two steps: 1. In Sec. IV B, we numerically evolve the full system in spherical symmetry. We choose initial data containing a single scalar QNM withn > 0, which is ideally suited to study nonlinear behavior. This is because, after some time, nonlinearly excited modes with n <n will come to dominate the ringdown, simply because of their longer decay times (despite their initially small amplitudes). We observe nonlinearities in all of our simulations, for a variety of perturbing overtone numbers and amplitudes. These simulations confirm that the horizon area grows as a result of a QNM flux as predicted by Eq. (32), and suggest that the dominant nonlinear process is AIME, as defined in Sec. III B.
2. In Sec. IV C, we further show that these nonlinear excitations are independent of how we generalize QNM initial data to a nonlinear setting. Moreover, we pin down the mechanism behind the nonlinear excitations. We do so by using linear simulations of a scalar field in an evolving background designed to mimic the backreacted black hole spacetime. We show that these simulations, which single out the effect of absorption, are an excellent approximation to the fully nonlinear runs. Linear simulations also allow us to explore AIME in the adiabatic limit.
A. Numerical and fitting setup

Numerical Method
We solve the full system of equations (1)-(2) in spherical symmetry using the numerical code developed in Ref. [38,54], which is based on a characteristic approach [55]. We restrict to spherical symmetry, using the following ansatz for the metric in ingoing Eddington-Finkelstein coordinates, The code uses a finite difference scheme with mixed second and fourth order radial derivative operators satisfying summation by parts [56,57], and fourth order Runge-Kutta for the time evolution. We compactify the spatial domain by using the inverse radial coordinate ρ = 1/r. This allows us to include the AdS boundary ρ = 0, so the computational domain is 0 ≤ ρ ≤ 1/r 0 , with the inner radius r 0 lying several grid points within the apparent horizon, and thereby excising the singularity (Fig. 2). Equations are discretized using a uniform grid throughout the compactified spatial domain. For the results below, we typically use 4801 points to cover the domain in ρ, and set dv/dρ = 0.4 for numerical stability. (See Appendix B for truncation error estimates.) Boundary data consist of the ADM mass M = 0.104, and the initial data is fully specified by the scalar field profile on the initial null slice φ(v = 0).
The system is solved by integrating radially, starting from the AdS boundary, along v = constant null curves to obtain the remaining field values and their derivatives; the scalar field is then integrated one step dv forward in time, and the procedure is iterated.
In this work, we choose a specific QNM with overtone numbern as initial data for the scalar field. The radial QNM profile is obtained by specifying the target overtone frequency and solving the radial Klein-Gordon equation on SAdS, following the methods detailed in Ref. [38]. The radial profiles φn(r) are normalized using the maximum norm with An setting the overall amplitude of the field. For large amplitudes An, a significant fraction of the ADM mass M will be contained within the scalar field itself. As a consequence, the black hole will have a mass significantly smaller than M and the actual black hole background will not precisely match the SAdS background for which the initial data were constructed. With this in mind, we can construct two types of initial data: • For "global-mass" initial data, we obtain the radial scalar field profiles for the target overtonen on a SAdS background of mass M . Then, this is used as initial data for the scalar field without updating the ADM mass M of the spacetime. Thus, the perturbation frequency and profile matches the final black hole mass (after the perturbation has been completely absorbed) but has a small discrepancy with the black hole mass at the beginning of the evolution.
• For "matched" initial data, we follow the same procedure outlined above and extract the irreducible mass of the black hole on the first slice v = 0. Next, we use this extracted mass to compute the new background mass, update the target overtone frequency ωn, and finally obtain the updated radial profiles. This process can be iterated, however we observe that with a single iteration we achieve very good agreement between the mass used for the initial data calculations and the mass computed on the first slice of the evolution. These initial data should therefore better represent a single QNM of the initial black hole.
The asymptotic expansion of the scalar field about infinity, with reflecting boundary conditions, takes the form [38,54] The quantity φ ∞ (v) is an output of the simulation, and contains the information about the mode content of the solution. We extract it at the AdS boundary as a time series φ ∞ (v). Moreover, the time coordinates v and t coincide at the AdS boundary; hence, unless otherwise stated, we will refer to this time series φ ∞ (v) as φ(t) for simplicity. Another useful output of the simulation that we extract is the apparent horizon area time series A H (v), which we denote A(v).

Fitting Method
We analyze the QNM content of the scalar field ringdown by fitting the time series with the model containing N + and N − positive and negative frequency modes, respectively. The latter are sometimes referred to as mirror modes [13], and satisfy ω −n = −ω * n = −ω n,R + iω n,I , with ω n,R > 0.
We determine the value of the amplitudes A ±n and phases ϕ ±n using a least squares fit to the boundary data, while holding the frequencies fixed. The QNFs used for the fit are computed based on the late-time mass of the black hole, using Leaver's method with the code presented in Ref. [38] using the. The fits are performed between an early time t 0 and a late time T . We assess the goodness of the fits by plotting the fit residuals and computing the mismatch between the model and the data, as described in Appendix B.

B. Nonlinear evolution
The typical nonlinear response of the scalar field to QNM initial data is exemplified by the ringdown shown in Fig. 3. The response begins with the perturbing mode ringing (in this case the first overtone,n = 1), with no sign of an "immediate response" phase. This can be explained by the fact that the initial data is a QNM defined on the entire ingoing null slice, from the AdS boundary up to the black hole horizon, as opposed to compactly supported initial data. At intermediate times, we observe a transition where eventually the decay of the first overtone causes it to drop below the fundamental,  which takes over for the rest of the evolution. At late times, as expected for asymptotically AdS black holes, we do not detect a power-law tail. We therefore use a pure QNM model (39) to fit the data from t 0 = 0 and T = 300M .
Eventually, after the transition, the dynamics are dominated by a mix of the positive and negative frequency fundamental modes, as seen in the oscillating decay from t = 110M onwards in Fig. 3. These modes are longerlived and smaller in amplitude compared to the perturbing first overtone (n = 1). This simple three-mode model suggested by visual inspection is confirmed by a least squares fit; see Fig. 12 and Table III in Appendix B. The residual with just the n = 1 included in the fit is well above the numerical error, but significantly decreases with the inclusion of the n = +0 and n = −0 modes in the model. The asymmetric excitation of the positive and negative fundamental modes should not come as a surprise; although the spectrum and the radial profiles are symmetric under ω → −ω * , X → X * (see Eqs. (7), (8)), so that X ω = X * −ω * , our choice of positive frequency initial conditions picks out a preferred sign.
The evidence for overtones above the perturbing mode, n > 1, or including n = −1, is very weak. The residuals do not improve significantly when adding these modes. In particular, residuals never improve significantly at early times, perhaps due to the fact that our fit model does not allow for evolution of the quasinormal frequencies and amplitudes. Furthermore, the amplitudes of the first three modes remain consistent with the minimal model as overtones are added, but the amplitudes of the higher tones are unstable; for example, the n = 2 mode amplitude doubles when adding the third overtone. Similarly, the  mismatch decreases significantly when reaching the minimal model, but fluctuates around the minimum value as overtones are added. Generically, nonlinear theory does not prohibit the excitation of higher overtones, hence they are likely present in the nonlinear data at late times. These higher overtones are hard to identify due to a combination of perturbatively small amplitudes and short decay times. For this reason, we limit the scope to the minimal fit models that include n < n pert overtones in the remainder of the section. We justify this choice further in Appendix B.
We verify that modes are excited nonlinearly by studying the excitations as a function of the QNM perturbation amplitude. We find that mode amplitudes behave cubically, as predicted in Eq. (35), with O(0.1 − 10) prefactors. We show this in Fig. 4 forn = 1 perturbations and in Appendix C forn = 2 perturbations, but we confirmed this pattern 6 up ton = 10. We verified that modes are excited cubically also in The backreaction timescale of our prototypical n = 1 mode is shorter than the oscillation timescales of the modes up to at least n = 3. This is also true for the other modes used as initial data in this work, n = 2, 3. Based on this estimate, their nonlinear excitations should behave non-adiabatically.
larger black holes, with r + L. We omit these results for brevity, and because for larger black holes the late time data is more susceptible to numerical error.
As discussed in Sec. III A, we expect QNMs to backreact on the black hole, increasing its horizon area over time. Assuming that the perturbing mode is the dominant source of flux across the horizon, we can use Eq. (32) to predict the change in area. We simply read off the perturbing mode horizon amplitude (appearing in Eq. (32)) from its numerically determined profile, as well as the black hole area on the initial slice. The result, shown in Fig. 5 forn = 1 andn = 3 perturbations, is in excellent agreement with the numerical evolution of the horizon area, and shows an area variation of 2-4%.
The agreement between the numerical and analytic area evolution confirms that its timescale is set by half the decay time of the dominant QNM in the ringdown. We can then use this timescale to determine whether the system will behave adiabatically or not. The comparison between the area evolution timescales and the mode oscillation frequencies (all shown in Table I) tells us that all our simulations fall under the nonadiabatic regime, as defined in Sec. III B.

C. Linear evolution with time-dependent background
The nonlinear simulations contain the full dynamics and mode excitations produced by the scalar field backreaction, with the expected cubic relation. To further isolate this effect we carry out a special linear simulation, where we evolve the scalar field on a time-dependent background. The spacetime is prescribed analytically to have a evolving mass corresponding to the dominant QNM flux through the horizon, as in Eq. (32). These linear simulations single out the effect of the absorption-induced mode excitation, allowing us to verify where the cubic excitation is actually stemming from.
In Fig. 6, we show the ringdown of the scalar field in nonlinear simulations with two choices of scalar initial data: the matched and the global-mass initial data defined in Sec. IV A. Both simulations contain modes excited nonlinearly and to a similar degree 7 . This tells us that the excitations are mostly independent of the preparation of QNM initial data. The mismatch between the background assumed for the initial QNM profile and the actual background on the initial slice could have been responsible for most of the nonlinear excitations-but this is not the case.
We compare these nonlinear calculations in Fig. 6 to the the result of a linear simulation where we evolve the scalar field using the same initial data and the same parameters as used in the (matched ID) nonlinear evolution, but on a SAdS background with a time dependent mass set by the imaginary frequency of the initial data. We control the time dependence of the background black hole area as follows, where τ back. = 2|ω I,n | with ω I,n the imaginary frequency of the initial data used, and ∆A = A fin − A init . The three parameters ω I,n , A init , and A fin are inputs of the simulation and are chosen to match the nonlinear simulations. Specifically, A fin = 0.104, and A init is the mass used to compute the ID, which will be different depending on whether we use the global-mass or matched ID. In linear simulations, it is the evolution of the background alone that produces the additional modes seen at 7 Matched initial data produces slightly stronger nonlinear excitations. This is because this choice of initial data for the scalar field induces a smaller black hole area on the initial slice, resulting in a larger area excursion to reach the final value (set by the ADM mass, the same in all simulations). The frequency of the initial perturbation, impacting the area evolution, is also slightly different in the two choices of initial data. late times in Fig. 6. This shows good agreement between linear and nonlinear simulations with the same choice of scalar initial data. This indicates that, among the possible sources nonlinearities in our system, AIME is indeed the dominant one. We also confirm that the excitations seen in the linear simulation are related to the initial data through a cubic scaling (assuming the black-hole mass change is dictated by the square of the perturbing field; see Fig. 4). This shows that absorption excites modes in the QN spectrum cubically, as argued in Section III B.
The linear setup, moreover, gives insight on the adiabaticity of AIME in our system. We would like to understand if the system is behaving adiabatically or nonadiabatically, and whether our criterion based on the mode frequencies in Tab. I is reliable. To investigate this we evolve then = 1 mode in a background with a modified timescale for the black hole area change given by Eq. (41) with τ back. = α × τ back. true = α × (2|ω I,n |).
In the nonadiabatic limit, α → 0, the background instantaneously jumps to its final configuration, producing the largest excitations. In the opposite limit, α → ∞, the background is frozen and no modes are excited that were not present in the initial data; we observe the presence of then = 1 overtone well into intermediate times. For α = 1, we recover the true evolution induced by the perturbing mode. Figure 7 shows the resulting evolution of the scalar field for a range of values of α, spanning the two extremes. Modes are excited with smaller amplitude when the background evolution is slower (large α). In those cases, we observe the mode's imaginary frequency shifting slowly in time. The clearest example in Fig. 7 is for α = 100, the slowest evolving background, where we observe the slope of the scalar field amplitude settling into the final frequency, at late times. This slow evolution of the scalar field amplitude is due to the time-dependence of imaginary frequency induced by the slowly evolving background.
The right panel of Fig. 7 shows how the excitation amplitudes converge to the nonadiabatic limit as α → 0. The true background evolution (α = 1) is itself close to the nonadiabatic limit, as its excitation amplitudes are only ∼ 30% and ∼ 70% smaller than the ones obtained with α = 0.01 for n = +0 and n = −0, respectively. We conclude that, when the background and mode timescale respect the condition τ back. τ osc. , the nonadiabatic approximation can give a good order-of-magnitude estimate of the modes excited through this nonlinear process. Figure 7 also shows an analytic estimate of the nonadiabatic excitation amplitude for n = +0. This estimate is obtained by computing the excitation coefficient of the n = +0 mode of the late-time black hole, with mass M . The initial data is set by the moden = 1 of the initial black hole. We regularize the overlap integral between the modes with the subtraction procedure described in Sec. II B. The result is in good agreement (within a factor of 2) with the amplitude measured in the simulations in the small α limit. We were not able to obtain a reliable estimate of the excitation coefficient of the n = −0 mode, as our regularization method has a precision comparable to the mode amplitude (∼ 10 −3 ).
The adiabatic limit, on the other hand, is harder to explore numerically. When artificially slowed down, most of the evolution of the background takes place when the amplitude of the perturbing mode has already decayed by orders of magnitude. This is shown by comparing the area evolution in the inset and the main plot in the left panel of Fig. 7. In this regime, it is hard to disentangle the effect of an adiabatic background evolution from the decay of the perturbation, as they both contribute to a reduced excitation amplitude; when AIME becomes relevant, the perturbation has already decayed a few orders of magnitude in amplitude. Indeed, the knee in the excitation coefficients as a function of the background timescale is consistent with the drop in the perturbing amplitude (right panel of Fig. 7). At intermediate timescales, the decayed perturbation amplitude can be used to obtain a better analytic estimate of the nonlinear excitation, as shown in Fig. 7. At very long timescales, on the other hand, the excitation amplitude does not fall as steeply as the perturbing mode. In this limit, the excitation could be explained by the partial area increase taking place before the complete decay of the perturbing mode.

V. TOWARDS ASTROPHYSICAL BLACK HOLES
In this section, we connect our results to gravitationalwave observations of ringdowns. Although so far we have considered scalar perturbations of asymptotically AdS black holes, the generic effect of overtones being excited due to an evolving black-hole background will occur just as well for gravitational perturbations of asymptotically flat black holes, e.g., in the ringdown of a merger remnant. Although a full analysis of this case is beyond the scope of this paper, here, building on the success of the analytic predictions in reproducing the nonlinear excitations seen in SAdS, we make similar predictions for gravitational perturbations of Schwarzschild black holes. We estimate the effect of AIME using the same methodology as described in Sec. III A. After summarizing the theoretical framework, we evaluate the excitation coefficients numerically to estimate the significance of the backreaction phenomenon for gravitational wave observations. Though here we restrict our analysis to non-spinning black holes for simplicity, we also outline the steps required to extend our results to Kerr.
At second order in perturbation theory, vacuum gravitational perturbations of Kerr were computed numerically for the first time in Ref. [20]. At this order, perturbations are excited through mode mixing, so that a linear perturbation of azimuthal number m = ±m can only excite modes with m = ±2m and 0. Our focus is on nonlinear excitations mediated by the black hole area evolution, which first appears at third order in perturbation theory. Naïvely, its perturbative order suggests that AIME might be subdominant. This might not be the case, as we discuss in Section VI.
We work again in the horizon penetrating, ingoing Eddington-Finkelstein coordinates, Eq. (23) with f (r) = 1 − 2M/r. The tortoise coordinate is defined as dr * /dr = 1/f (r) and the horizon is located at r + = 2M .

A. QNM horizon flux and area evolution
We begin by deriving the evolution of the Schwarzschild black hole area induced by the QNM flux. We follow closely the derivation of Teukolsky and Press [33], adapting it to complex frequencies.
The area of a Schwarzschild black hole is given by Thus, the evolution of the mass and the area are related by The energy flux can be expressed in terms of the Weyl scalar ψ 0 , giving Taking the horizon limit, we find This is more conveniently expressed in terms of the amplitude of the Weyl scalar ψ 4 , The Teukolsky-Starobinski identities relate the amplitudes of ψ 4 and ψ 0 via the Starobinsky constant C, The Starobinsky constant for complex frequencies is obtained via the analytic continuation where 8 E = E |s| lm is the eigenvalue of the spinweighted spherical harmonics s S lm , normalized as dθ| s S lm (θ)| 2 sin θ = 1. The amplitude of the Weyl scalar ψ 4,r+ at the horizon can be related to the amplitude at infinity ψ 4,∞ through the global solution, which we evaluate numerically. Finally, ψ 4,∞ is related to the two gravitational wave polarizations at infinity ψ 4 (r → ∞) = 1 2 (ḧ + − iḧ × ) through the asymptotic limit, Dots denote derivatives with the respect to the Boyer-Lindquist time t. 8 Note that the eigenvalue E is independent of the sign of s in the absence of the black hole spin.

B. Nonlinear excitations due to absorption
The nonlinear excitation of modes due to absorption depends, as we argued in the SAdS case, on how fast backreaction takes place compared to the mode oscillation. Figure 8 shows that for the l = m = 2 mode, which dominates most observations of ringdowns, the fundamental and first overtone can both give rise to nonadiabatic background evolution for spins a 0.9. In particular, the AIME excitation of the n = 3 mode driven by the absorption of the fundamental mode should be very similar to the SAdS example discussed in the previous section. Thus, from the point of view of modes in Schwarzschild, backreaction happens almost instantaneously. This justifies using an instantaneous projection to estimate nonlinear excitation amplitudes, as done in Sec. IV C. In the following, we summarize the theory of the Schwarzschild excitation coefficients, specializing to the case of QNM initial data.
Given initial data in Schwarzschild, the excitation coefficient of a QNM can be computed in the Laplace-transform formalism [41]. The Laplace transform of the radial Teukolsky equation reads where the radial function is defined via the separation of the Teukolsky perturbation variable and where with λ = E − s(s + 1). We assume that the initial data is a QNM with frequencyω = ωn ,l,m , so that R s lm (t = 0) = R(ω) anḋ R s lm (t = 0) = −iωR(ω), where we suppressed the QNM indexes. Then, the source in Eq. (51) reads [58] The excitation coefficient for a mode of frequency ω is defined as where we defined the overlap integral The integral is defined on a complex deformation of the real line, C, so as to be convergent in the two limits r * → ±∞. Based on the asymptotic form of the QN radial functions at infinity and at the horizon, the contour is given by [34,35] r * → Λ+ρe iθ , θ ∈ [−π−arg(ω+ω), − arg(ω+ω)], (58) for r * > Λ > 0, and for r * < λ < 0, with ρ > 0 in the two regions. The resulting integral, which is the sum of three contributions (r * < λ, r * > Λ, and the intermediate integral on the real line) is independent of Λ, λ, and the choice of θ within its limits. We adopt the following definition for the excitation factors [34,35] B which guarantees the overall normalization of the excitation coefficients. The excitation factors of a Schwarzschild black hole can also be computed from the asymptotic expansion of the solutions of the Teukolsky equation, see for instance Refs. [59,60].

C. Numerical estimates
As an example, we consider an initial perturbation composed of a single fundamental (n = 0) mode with 9 l = 2, m = 0, and estimate the nonlinear excitation of overtone modes with the same angular numbers under the approximation that the mass evolves instantaneously. The timescales for these modes are shown in Fig. 8 (in Schwarzschild, mode frequencies are independent of m).
We find the QNFs and the radial functions for s = −2 using the Black Hole Perturbation Toolkit [61] for a given initial black hole mass M init . For initial data consisting of a single (the fundamental) QNM, we then estimate the final mass of the black hole due to the inward flux using Eqs. (46)- (50). We then re-expand the initial mode withω = ω 0 (M init ) in terms of the overtones of the final black hole, ω = ω n (M final ). We verified that the overlap integral has the expected scaling with the mass variation, O ∼ ∆M for ∆M M init , by evaluating it numerically for a range of mass increments; see Appendix A.
The excitation coefficients of the first three overtones are given in Tab. II, for two values of the initial perturbation amplitude. For simplicity, we only consider 9 We choose a mode with vanishing m because it does not backreact on the black hole spin, and we can consistently study the problem in a Schwarzschild background.
To visualize the effect of these nonlinearities, we show the corresponding total gravitational wave at infinity in Fig. 9. For simplicity, we assume that all modes are sourced with the same (zero) phase. The impact of AIME on the gravitational wave signal is significant: it is above the percent level for both choices of the perturbation amplitude. The nonlinear excitation in the fundamental mode has the strongest impact, both at early and late times. The effect of the nonlinearly excited overtones is also significant at early times.
Following a black hole binary inspiral, the merger itself may seed overtones. In a realistic scenario, therefore, the nonlinearly generated modes will add to those sourced directly. Furthermore, we emphasize that in producing Fig. 9, for simplicity, we did not take into account any delay between the linear and the nonlinear excitation. The precise onset of nonlinear excitations requires further investigation, as discussed for SAdS.

D. The spinning case: expectations
In this work, we studied nonlinearities in the ringdown of nonspinning (Schwarzschild and Schwarzschild-AdS) black holes. The generic merger of an astrophysical black hole binary will, however, result in a spinning (Kerr) black hole. It is therefore essential to extend the study of nonlinear effects to the spinning case to make better contact with gravitational wave observations. A detailed analysis of the Kerr ringdown beyond linear order is left for future work. Here we summarize the main conceptual differences compared to the nonspinning case.
Backreaction onto a spinning black hole will change its angular momentum in addition to its mass, according to for a single mode; see Appendix D. As a consequence, a linear perturbation will need to be re-expanded not only in terms of modes with different black hole mass, but also with different spin. Since the angular solution (a spinweighted spheroidal harmonic) now depends on the spin, different l modes would generically be excited through the nonlinear process. Modes with different m will remain orthogonal. Note that the excitation amplitude will depend on the perturbing mode number m through the change in the black hole spin, as predicted by Eq. (61).

VI. DISCUSSION AND CONCLUSIONS
The final phase of a black hole binary merger, the ringdown, is an ideal testing ground for general relativity. Some time following the merger, the signal should be well described by linear perturbations-in particular, QNMsof an isolated black hole with fixed mass and spin. The goal of the black hole spectroscopy program is to identify these modes and accurately measure their frequencies in gravitational wave data.
Our goal in this work was to explore the dynamics of perturbed black holes beyond linear order. With numerical simulations of a toy model (a scalar field coupled to gravity in SAdS), we elucidated an important nonlinear effect, entering at higher order in perturbation theory and with secular characteristics. We verified that the flux of QNMs through the horizon can change the black hole mass during the early ringdown. This is a second order effect on the metric, which affects the QNFs (which depend on the background). The evolution of the black hole can also excite QNMs at third order in perturbation theory. We call this effect AIME. Although for gravitational waves (in contrast to a scalar field), nonlinear effects start at second order (e.g., mode doubling), the mechanism described here secularly accumulates as the initial modes decay, and can become comparatively significant.
Through detailed analysis in asymptotically AdS blackhole contexts, we find that AIME can be captured with fully nonlinear calculations, as well as through linear simulations on a time-dependent background. We also provide a simple method to estimate the amplitude of the excited modes, based on QN excitation coefficients. This method hinges on the fact that the background evolves non-adiabatically compared to the mode oscillations. We expect that the intuition developed for the toy model carries over to asymptotically flat spacetimes. We estimate the amplitude of AIME in Schwarzschild and show that this effect could make a percent-level contribution to the ringdown following a merger.
Our analysis has implications for ringdown studies of generic black hole merger gravitational wave events: • Overtones are excited generically and dynamically in the ringdown at higher orders in perturbation theory.
• The mass (and spin) of the remnant black hole necessarily evolve over time in the aftermath of the merger. Ringdown analysis looking for QNMs near the peak of the signal should take into account that their frequencies can evolve with time. In particular, inference based on signals detected early in the ringdown, while showing a higher number of modes, will bias the extracted value towards a smaller black hole mass. Later on, longer-lasting lower modes will allow for a more accurate measurement of the final mass. For this reason, higher overtones could actually bias rather than improve the inference of the (final) remnant properties, as seen in some numerical simulations [24,25].
• The amplitude of the QNMs should evolve also over time in the early ringdown as a result of AIME. The evolution should happen on the same timescale as the black hole mass and spin.
• Overall, our results suggest that ringdown analysisespecially when starting early after the mergercould benefit from being made more flexible, allowing for the evolution of mode frequencies and amplitudes.
Beyond these practical observations, the larger question of nonlinear behavior in the post-merger regime is interesting in its own right. At second order in perturbation theory, a known effect is the energy transfer between modes with different azimuthal number m, explored in Ref. [20]. Despite appearing at higher order in perturbation theory, AIME could be as significant to the ringdown as second order mixing. Although a direct comparison is not straightforward, the excitations of the Weyl scalar measured at null infinity in [20] are of the same order (10 −4 − 10 −5 for a 10 −2 perturbation) as the analytic estimates reported in our Table II. Indeed, AIME can be thought of as a secular effect, due to the cumulative change in the black hole mass, in contrast to mode mixing, whose source oscillates and strictly decreases as the linear mode decays. This may upend the order counting argument. Second order calculations might also be affected by the secular mass evolution discussed in this work, depending on the stage at which they are carried out. Finally, our mechanism can excite all the overtones of the initial perturbation of azimuthal numberm, and can therefore be complementary to mode-mixing excitations of different m modes. Overtones of the l = |m| = 2 mode are already being used in analyzing the ringdown of gravitational wave events to extract the mass and spin of the remnant [10], while higher angular modes in the ringdown have proven more elusive [11] (see however, [12]).
While our work elucidates nonlinear mode excitations in perturbed black holes in general settings through lessons drawn in asymptotically AdS scenarios, further studies in asymptotically flat spacetimes would be valuable. In particular, linear simulations on a time-dependent background, as applied in this work to asymptotically AdS black holes, could be easily used to further explore AIME in asymptotically flat black holes. The problem becomes richer in Kerr, where modes with different angular number l could be excited as well. Particularly intriguing would be exploring our mechanism for mode excitation in the extremal limit of Kerr, where QN frequencies seem to suggest that the adiabatic limit might apply (Fig. 8) so AIME would be suppressed. AIME should also apply in the context of superradiance, where modes, although not absorbed, nevertheless modify the background spacetime. Finally, modes should be similarly excited when black holes accrete matter.
Finally, our work shows that black holes do not act like censors of nonlinearities after merger. On the contrary, by efficiently absorbing energy and growing, they mediate the nonlinear AIME. Overlap integral between scalar QNMs. The background is SAdS with black hole mass M = 0.104. As the regulator R is taken closer to the horizon, the overlap between modes with different overtone number goes to zero. Black lines are power-law fits, ∼ ax b ; we find b > 0 for all combinations of modes. fit residuals as we add more modes to the model in Fig. 12. In all these diagnostic tools, the evidence for overtones above the perturbing mode, n > 1 and including n = −1, is very weak: the residuals do not improve significantly when adding these modes, and while the amplitudes of the first three modes remain consistent with the minimal model as overtones are added, the amplitudes of the higher tones are unstable: the n = 2 mode amplitude doubles when adding the third overtone. Likewise, the mismatch decreases significantly when reaching the minimal model, Fit residuals for models with increasing mode content, for a scalar field ringdown (black, solid) withn = 1 and An = 0.1 initial data. The minimal model providing a good fit to the data contains two fundamental modes and the first overtone (green). Fit residuals do not decrease significantly when higher overtones n > 1 are added (pink), which we interpret as weak evidence for their presence in the data. The numerical noise is shown for reference (black, dashed).  Table III. Mismatch between the data and the fitand the amplitudes of the overtones for the fits in Fig. 12.

Appendix C: Higher overtone perturbations
We confirmed that the amplitude of scalar nonlinear excitations in SAdS are cubic in the perturbation amplitude, for perturbation of any overtone number. Figure 13 shows the case ofn = 2 perturbations. At highern, the number of parameters to fit increases. Although the analysis becomes more challenging, the cubic relation is still distinguishable at high perturbation amplitudes.

Appendix D: Angular momentum flux
We compute the flux of angular momentum of a complex scalar field in a spherically symmetric black hole spacetime along the same lines of the energy flux calcula-tion in Section III A, see also Ref. [52]. The Killing vector associated with the ϕ symmetry is given by ϕ a = ∂ ϕ and defines a conserved angular momentum current, J a ang = T ab ϕ a .
The integrated angular momentum flux across the event horizon reads dJ dv = horizon dΩ n a J a ang .
Substituting the asymptotic form of the scalar field at the horizon, Eq. (24), we find Putting together this result and the energy flux computed in Eq. (27), we obtain the generalization of the standard angular momentum-energy relation to perturbations with complex frequencies, Under the assumption of real frequency, the relation reduces to the well known dJ = m ω dM .