I. INTRODUCTION

According to the theory of general relativity, binary black holes undergo an inspiral and merging process, generating ripples in spacetime known as gravitational waves. These waves carry valuable information about strong-field gravitational dynamics and are now routinely detected by the LIGO-Virgo-KAGRA Collaboration [1][2][3]. Subsequent to the coalescence of the two distinct black hole horizons into one, a phase known as the ringdown occurs, signifying the conclusion of the merger. In this phase, the remnant black hole settles down to an equilibrium state, emitting gravitational waves with distinctive frequencies that provide insights into the black hole's properties [4][5][6].

At some sufficiently late time in the merger process, the ringdown waveform ought to be well described by linear black hole perturbation theory [7][8][9]. In particular, the Teukolsky equation, governing dynamics of linear perturbations of Kerr black holes, admit eigensolutions known as quasinormal modes (QNMs), which manifest as damped sinusoids with frequencies uniquely determined by the mass and spin of the background Kerr geometry. One can then in principle independently measure the remnant mass and spin by fitting the ringdown signal with these QNMs [10][11][12]. However, there has been much debate regarding how soon after merger one can start the fit using a linear model of ringdown and how many overtones can faithfully be measured (whether from simulations or actual data) .

The debate is spurred by several complications that arise in the analysis. First, QNMs are not a complete description of black hole perturbations even at the linear level (see, e.g., the discussions in [34,35]). In contrast with, for example, certain quantum mechanics problems, where the selfadjoint Hamiltonian operator yields a set of eigenmodes that form a complete and orthonormal basis, the QNM solutions of a perturbed black hole are not complete, due in part to the dissipative nature of the theory. 1 Moreover, the transient, typically representing the high frequency non-QNM content of the perturbation, and the power-law tail, arising from backscattering of the gravitational radiation off spacetime curvature, are present in the ringdown waveform for generic initial data [38][39][40][41][42]. Their presence may bias the fitting using a model that only takes the QNMs into account [43][44][45].

Furthermore, a plethora of nonlinear effects may dominate the ringdown waveform at early time. At second order, linear quasinormal modes couple to give rise to quadratic quasinormal modes [46][47][48][49][50][51][52][53][54][55][56][57][58], which can be louder than their linear counterparts in certain angular harmonics following binary mergers [59][60][61][62][63].

Another nonlinear effect that has gained recent attention is the change in the black hole mass and spin due to radiation backreaction [64][65][66][67][68]. Given that in comparable mass mergers ∼5% of the energy of the binary is radiated in gravitational waves in the coalescence and early ringdown phase [69], one could expect a significant fraction of that would register as a change in the remnant mass and spin over a few light crossing times. After that, the remnant black hole's properties stabilize exponentially.

At the perturbative level, changes in the mass and spin can be introduced at the linear level through corresponding zero frequency modes. However, since these modes are zero frequency, they do not represent changes to the background geometry due to dynamical perturbations (including QNMs and transients). The latter only cause backreaction at second order in perturbation theory and hence will only produce an effective change in properties of linear QNMs at third order.

Since even second-order perturbations of Kerr black holes are not yet well understood (see, e.g., [70][71][72][73][74]), works attempting to address imprints of mass and spin change on QNMs have been heuristic in nature. Reference [75] studied the frequency modulation in the ringdown signal of a massless scalar field due to backreaction. Reference [64] explored what they call absorption-induced mode excitation (AIME) for black holes in asymptotically anti-de Sitter spacetime, where they showed that the overtones (fast decaying QNMs) can lead to excitation of the fundamental mode (slowest decaying QNM) due to nonlinear backreaction on the background spacetime; they further calculated the change in fundamental mode amplitude in asymptotically flat Schwarzschild due to AIME using the "sudden masschange approximation." Reference [66] developed a theoretical framework to model the backreaction of a scalar field on the background spacetime, which yields excellent agreement with numerical simulations. Reference [67] found evidence for time-dependent frequencies in Vaidya spacetimes, which describe accretion of null dust onto a black hole, yielding a prescribed mass function. Concurrent to our work, Ref. [76] is investigating AIME for asymptotically flat black holes using a perturbative approach. They find that the prograde QNM can mix into the retrograde QNM, potentially yielding observables of direct astrophysical relevance.

In this paper, we work without the aforementioned assumptions or approximations and directly investigate the imprints of nonlinear absorption during black hole ringdown through numerical experiments. To distinguish linear and nonlinear contributions, we adopt a close-limit type framework [77][78][79][80][81][82], in which we evolve a single perturbed black hole by numerically integrating both the linearized and nonlinear Einstein equations. With a consistent gauge choice and initial data to linear order, we attribute differences in the extracted gravitational waveform to nonlinear effects.

A similar numerical approach has been adapted in Refs. [81,82] for axially symmetric perturbations of a Schwarzschild black hole; here we focus on nonaxisymmetric perturbations of a Kerr black hole with dimensionless spin of 0.7, the expected spin for the remnant formed from the quasicircular inspiral of two nonspinning equal mass black holes [69,83]. Furthermore, in the earlier works the nonlinear waveforms were normalized by the initial Arnowitt-Deser-Misner mass (i.e., including all the radiation in the initial data), which results in a mismatch of the time axis compared to the linear waveforms. Here, we normalize the time axis by the initial mass of the background black hole to give a consistent comparison between linear and nonlinear evolution. We also present a detailed ringdown analysis in terms of the wave's quasinormal mode content, quantify different nonlinear effects, and discuss their relevance to observations. We focus on the dominant quadrupolar l ¼ m ¼ 2 harmonics of the waveform, which are traditionally thought to be dominated by linear contributions. Here we show that nonlinear effects can lead to relative changes of order unity in the amplitude and phase of certain QNMs for the largest, astrophysically reasonable perturbations we consider. Despite these large nonlinear distortions, we find that the resulting waveform can still be well modeled by a linear superposition of QNMs. At a first glance this seems counterintuitive, however ringdown modeling does not make use of theoretical estimates of the amplitude and phase of excited QNMs, and by the time the waveform is in a regime where the fitting algorithm is stable, the instantaneous frequencies of the QNMs have largely settled down to those of the final remnant black hole. Moreover, the time window over which modes can stably be extracted is essentially the same for the linear versus nonlinear results, suggesting that transients are the main culprit preventing straightforward interpretation of the coalescence/early ringdown phase of a merger in terms of black hole perturbation theory, linear or otherwise.

A. QNM conventions

In this paper, individual QNMs are denoted by four indices ðp; l; m; nÞ, following conventions in Ref. [84]. QNM angular structure is described by spheroidal harmonics with indices l and m, while their radial structure is indicated by the index n: n ¼ 0 QNMs are fundamental modes, while n > 0 are overtones, ranked by increasing decay rate. For a given ðl; m; nÞ, there are two QNMs labeled by p, signifying modes with wavefronts corotating (prograde, p ¼ þ) or counterrotating (retrograde, p ¼ -) with the black hole. We define A ðp;l;m;nÞ as the complex amplitude for each QNM, and ω ðp;l;m;nÞ ðM; χÞ its complex frequency, which is a function of the black hole mass M and dimensionless spin χ. To distinguish amplitudes of QNMs from linear and nonlinear waveform, we use superscripts, e.g., A linear ðp;l;m;nÞ and A nonlinear ðp;l;m;nÞ .

II. NUMERICAL SETUP

To obtain initial data for a single perturbed black hole, as in Ref. [56], hereby referred to as paper I, in which we investigated the spin and initial data dependence of quadratic quasinormal modes, we use the conformal thin sandwich approach developed in [85]. First, we superpose the background Kerr metric in spherical Kerr-Schild coordinates g 0 ij [86] with the following initial guess h 0 ij for the metric perturbation:

where A, r 0 , and w are the amplitude, radial location, and width of the initial perturbation, respectively, and Y lm ij is the pure-spin tensor spherical harmonics following conventions in, e.g., Ref. [87]. Here (in contrast to paper I) we are choosing time-symmetric initial data. 2 This trivially solves the momentum constraints, but more importantly for our purposes means the background shift vector is not modified, making for more straightforward comparison between the linear and nonlinear evolution. The Hamiltonian constraint is then solved for a conformal factor ψ rescaling the initial guess for the spatial metric, after which the perturbation can be read off as

This procedure yields consistent initial data for the nonlinear and linearized evolutions.

As mentioned in the Introduction, we restrict the initial guess for the perturbation to only have l ¼ jmj ¼ 2 harmonics in the angular sector and the background black hole to have a dimensionless spin of 0.7 initially. We further fix r 0 ¼ 30M and w ¼ 2M, though we have checked that our results are largely insensitive to such choices.

We then numerically evolve this family of initial data, varying A from the perturbative to the nonperturbative regime, using the Spectral Einstein Code (SpEC) for both the linear and nonlinear evolution. The nonlinear evolution employs the generalized harmonic formulation of the Einstein equations [88][89][90]. For the linearized evolution, we use the implementation by Refs. [91,92], which is a consistent linearization of the generalized harmonic equations about the given black hole background. In both cases we use a freezing gauge condition, where the gauge source functions are set to be those of the Kerr-Schild background. The first-order gauge source function for the linearized evolution is set to zero for consistency.

III. RESULTS

We performed linear and nonlinear evolutions of perturbed Kerr initial data with four amplitudes evenly spaced on a logarithmic scale, labeled with amp 0-3 from low to high. Amp 0 results in a relative mass and spin change of less than 10 -4 , while amp 3 causes a change in the mass of 3.4% and dimensionless spin of 5.6%, where we use the Christodoulou mass and dimensionless spin calculated through approximate Killing fields [93,94]. In Fig. 1 we show the amplitude of the linear and nonlinear l ¼ m ¼ 2 spin-weighted spherical harmonic component of the waveform extracted at a coordinate radius of 100M with the smallest (amp 0) and largest (amp 3) amplitude initial data. 3For the largest amp 3 run, the l ¼ m ¼ 2 waveform receives a distortion in the waveform amplitude of ∼30%. At a first glance this is surprising in that the amplitude distortion is a third-order effect (see Fig. 2). Though what this indicates is that a perturbation resulting in a 3.4% change in the black hole mass is well into the nonlinear regime, consistent with the expectation that the analogous "perturbation" (as measured by mass change) caused by a black hole merger should be considered a rather nonlinear excitation of the remnant.

Before describing the nonlinear changes to the phase evolution of the waveforms, it would be useful to understand their QNM content. Using the methods described in paper I, we fit the l ¼ m ¼ 2 component of the linear and nonlinear waveforms with superpositions of QNMs including up to the second overtone, using the final mass and spin to define the corresponding QNM frequencies. We then compute the mismatch as a function of fitting start time between the best-fit QNM model ψ QNM 22 and the numerical relativity waveform ψ NR 22 , defined as

:

In the bottom row of Fig. 1, we plot this mismatch as a function of t 0 for both the linear and nonlinear waveforms. At early time the mismatch for both the linear and nonlinear waveform is high, which is as expected as here the signal is dominated by a transient. From a time close to the peak amplitude in the waveform onward, the mismatch drops significantly. Surprisingly, the improvement in mismatch, both in magnitude and when it occurs, is largely insensitive to nonlinearities. This suggests that the time when nonlinearities are most relevant in modifying the properties of the QNM excitation is also when the transient is still a significant component of the signal, and hence it is the latter that is ultimately the limiting factor in how early on in a signal a QNM analysis can be used to interpret the ringdown. Figure 2 illustrates the nonlinear distortion to the complex QNM amplitudes extracted from the linear and nonlinear waveforms for three dominant modes: the fundamental prograde and retrograde modes ðAE; 2; 2; 0Þ and the first overtone for the prograde mode ðþ; 2; 2; 1Þ. The QNM amplitudes are extrapolated to the peak of the waveform, with uncertainties estimated from stability of fit [56,96]. What is shown there are the relative differences FIG. 1. Top: comparison between waveforms from linear and nonlinear evolution of a single perturbed black hole. The time axis is normalized by the initial black hole mass M i , and t ¼ 0 at the peak of the linear waveform. Waveforms from nonlinear evolution are shown as red solid lines and those from linear evolution as orange dashed lines. We further plot the difference between the waveform amplitudes, i.e., rMjjψ linear 4 jjψ nonlinear 4 jj, using thick blue lines. We note that the amplitude difference does not take into account the accumulation of a phase difference, which we show in Fig. 3. We find good agreement between the nonlinear and linear waveforms at the lowest amplitude (left); however, near order unity relative disagreement (more clearly evident in Fig. 2) is found at the highest amplitude (right). The latter causes a change of 3.4% (5.6%) in the mass (dimensionless spin) of the black hole in the nonlinear evolution. Bottom: mismatch between the best-fit QNM model (including up to the second overtones) and the waveform varying fitting start time, where we used the late-time mass and spin extracted at the horizon for the nonlinear waveforms. For the highest amplitude case, the mismatch for the nonlinear evolution is only marginally larger than its linear counterpart at late time, suggesting that a linear superposition of QNMs could still be a good model. In all panels, numerical error is estimated using the difference between the highest and second highest resolution nonlinear evolution waveforms; linear evolution shows similar error characteristics.

between nonlinear versus linear evolution as a function of amplitude (see the Appendixes for plots of the absolute amplitudes and the stability of the QNM fitting with time). The scaling of the difference with amplitude is as expected from perturbation theory.

The remaining QNMs, namely n ≥ 2 prograde overtones and any overtones for the retrograde mode, are not stable for a meaningful time window (∼10M) in either the linear or nonlinear waveform, suggesting the linear transient is obstructing the extraction of these modes.

Finally, we return to the drift in waveform phase due to nonlinear effects. We illustrate this in Fig. 3 by showing the phase difference (blue line) between the highest amplitude amp 3 nonlinear (red line) and linear (dashed orange line) waves versus time. At early times the two cases are in phase, but then around the time corresponding to the peak amplitude (t ¼ 0) there is a rapid shift in frequency difference between the two. This results in a phase difference that afterward is dominated by a linear drift in time, the slope of which is consistent with a change of the prograde mode frequency due to the change in mass and spin of the black hole; i.e., the slope of the solid black line is ℜfω ðþ;2;2;0Þ ðM i ; χ i Þω ðþ;2;2;0Þ ðM f ; χ f Þg, with M i and χ i (M f and χ f ) the initial (final) mass and spin. Subdominant to the linear phase drift is a small oscillation that eventually decays away. The oscillation appears to come from beating of the prograde ðþ; 2; 2; 0Þ mode with the frequency and amplitude-shifted retrograde mode ð-; 2; 2; 0Þ. This oscillation eventually decays away as the retrograde mode decays faster than the prograde mode. The dashed black line in the figure combines this beating effect with the linear drift, giving an excellent match to the late-time full nonlinear waveform.

IV. DISCUSSION AND CONCLUSION

We investigated nonlinear effects in the gravitational wave perturbation of a black hole by comparing evolution of the same initial perturbation solving the full Einstein equations versus the Einstein equations linearized about the given black hole background. Our main goal for doing this was to investigate how perturbations that result in FIG. 2. Shift of the complex QNM amplitudes due to nonlinear effects. On the vertical axis we plot the fractional amplitude difference between the linear and nonlinear waveform for the dominant QNMs jA linear ðp;l;n;mÞ -A nonlinear ðp;l;n;mÞ j=jA linear ðp;l;n;mÞ j. On the horizontal axis we used the linear A linear ðþ;2;2;0Þ as a proxy for linear perturbation amplitude. We find that the QNM amplitudes drift from the linear approximation at third order in the amplitude (hence the fractional amplitudes drift at second order), in accordance with formal perturbation theory. Despite being a higherorder effect, for the largest amplitude perturbation we consider, relevant to astrophysical mergers, the relative amplitude changes are close to ∼50%. We stress that this drift in complex amplitude comes from both the absolute amplitude and the phase. FIG. 3. Phase difference between the highest amplitude (amp 3) linear and nonlinear waveform. The phase of the nonlinear waveform is the solid red line, the linear waveform the dashed orange line, and their difference the solid blue line. Following the peak of the waveform (t ¼ 0), the phase difference can be characterized by a linear drift overlayed with an oscillation that decays away at late time. The linear drift is caused by the change in the ðþ; 2; 2; 0Þ frequency due to the mass and spin change; the solid black line is a "prediction" of this effect using the corresponding QNM frequencies. The oscillation is caused by the drifted retrograde ð-; 2; 2; 0Þ mode beating against the prograde mode, which decays with time since the retrograde mode decays faster than the prograde mode. We include this effect in the prediction of the difference as the dashed black line.

astrophysically relevant changes in the mass and spin of the black hole affect the QNM excitation and their extraction from the full signal using models based on linear QNMs. For such large amplitude perturbations we find an order unity relative change in the amplitude of the dominant quadrupolar QNM excited. However, despite such a large difference, we find almost no adverse effect in extracting this mode, its retrograde counterpart and its first overtone, using a fitting procedure based on the linear problem. In other words, the time in the fitting procedure where the linear waveform becomes stable is almost identical to the time when it becomes stable for the nonlinear waveform, and the mismatch in either waveform compared to the relevant sum of linear QNMs extracted becomes similarly small. The reason for this appears to be that the time when the nonlinear interaction is most relevant in affecting the amplitude of the QNMs, and indirectly their frequencies through the change in the black hole's mass and spin, is also when the transient part of the signal obscures the interpretation of the waveform as a sum of QNMs.

Our findings would apply to more generic initial conditions despite the use of a simple initial data configuration, namely a single tensor spherical harmonic. When multiple angular modes are present, the only added complication is quadratic coupling of quasinormal modes. However, the resulting quadratic quasinormal modes can be easily incorporated into the fitting model by adding their corresponding frequencies, and the challenge remains the transient, which already presents at the linear level where the angular modes decouple.

Though the particular form of the transient is a function of our initial data, it is difficult to imagine that an effective transient in the analogous problem of a remnant formed via a black hole merger will be any less significant. The "good news" in that regard is our work implies that linear QNM fitting for the ringdown in black hole mergers, whether from numerical predictions or detected events, does not need to model a changing mass and spin. The "bad news" is that even were a model of the changing mass and spin developed, it might not help in extending the time over which a stable extraction of QNMs can be obtained due to obscuration by the transient, which is already present at the linear level. 4 Note that our results certainly do not imply that nonlinear effects cannot be observed in the ringdown. For example, quadratic QNMs are present, though they could similarly be incorporated into a linear QNM model by including the appropriate frequency-doubled modes. A significant caveat on the above broader reaching implications of these results is that we have only studied an initial χ ¼ 0.7 spin black hole and a limited family of initial conditions.

Our framework can be easily extended to investigate other effects during ringdown, such as the quadratic QNMs just mentioned or precession. Regarding the latter, one could simply rotate the pulse relative to the black hole to study precession of the remnant spin due to radiation backreaction and test the toy model proposed in Ref. [96]. We leave a study of these extensions and a more thorough exploration of parameter space to future endeavors. period of fitting start time (∼10M). And the instability in the fitting is similar between the linear and nonlinear waveform, suggesting that the linear transient is the dominant obstruction for measuring these faster decaying modes. The exact cutoff in the extractability of overtones is dependent on initial data, but we expect that the transients always dominate over the nonlinear effects for obstructing the measurement of overtones for generic initial data.

In Fig. 6, we plot the horizon mass as a function of time on top of the linear and nonlinear waveform for the highest amplitude (amp 3) run, where the horizon time is shifted by 100M to roughly match the waveform time. We note that this matching is very heuristic as the horizon is by definition causally disjoint from the exterior. We find that indeed most of the mass change occurs before the peak of the waveform, during which the QNM fitting model fails even for the linear waveform [44]. FIG. 4. Amplitudes of QNMs as a function of fitting start time for the linear (dashed) and nonlinear (solid) waveforms. The QNM is stably extracted if its amplitude remains constant for a sufficiently long period of fit starting time. We note that the nonlinear effects of changing mass and spin does not seem to affect the stability of quasinormal mode fitting, even for the first overtone at the highest amplitude (blue lines in the rightmost plot). FIG. 5. Same as Fig. 4 but for overtones that are not stable. For our initial data, the prograde overtones with n ≥ 2 and any of the retrograde overtones cannot be stably extracted in both the linear and nonlinear waveforms. each spherical shell is summarized in Table I. We show convergence of the linear and nonlinear amp 3 waveform in Fig. 7 and convergence of the L 2 norm for the generalized harmonic constraint (see Ref. [88] for its definition) in Fig. 8. FIG. 7. Convergence of the linear (left) and nonlinear (right) waveform for the amp 3 runs. The resolutions are labeled by Lev0 (lowest) through Lev3 (highest); see Table I for details. We find the residual between different resolutions decreases as expected. FIG. 8. L 2 norm for the generalized harmonic constraint as a function of simulation time, from the largest amplitude (amp 3) run. The largest constraint violation is observed when the pulse hits the horizon, which converges away with resolution.

TABLE I. Number of collocation points/spectral elements for each resolution level (Lev). Lev N r N L 0 10 14 1 12 16 2 14 18 3 16 20