Introduction. The birth of gravitational-wave (GW) astronomy [1] marks a new era in the exploration of strong-field gravity [2,3]. As the simplest macroscopic objects cloaking curvature singularities, black holes (BHs) play a special role as astrophysical laboratories to test gravity and to search for new physics [4][5][6][7][8][9]. The structure and dynamics of BHs in our Universe is well described by the two parameters (mass M and angular momentum J) characterizing the Kerr metric. In general relativity, the perturbed BHs formed in a binary merger approach a stationary state by emitting GWs in a discrete set of characteristic quasinormal modes (QNMs) with complex frequencies determined only by M and J. The "black hole spectroscopy" program consists in observing these "ringdown" waves, measuring the QNM frequencies, using them to estimate mass and spin [10], and (if more than one mode can be observed) test that the remnant is indeed consistent with a Kerr BH [11][12][13][14][15]. The observability of QNMs depends crucially on their excitation in the merger process. Even within linear perturbation theory, where one only considers linear metric perturbations to Einstein's equations in the Kerr background, determining which modes are excited is a formidable problem [16][17][18][19][20][21][22][23][24][25].
General relativity is an intrinsically nonlinear theory. The merger of two comparable-mass BHs leading to a perturbed Kerr BH is one of the most violent processes in the Universe, where these nonlinearities should play an important role. It is therefore surprising that merger simulations in numerical relativity result in a very smooth transition from inspiral to merger and ringdown [26,27]. Where are the nonlinearities of general relativity? This state of affairs has led many (including some of us) to conjecture that nonlinear effects may be hidden behind the horizon, suppressed by the presence of a photonsphere, or even absent altogether(see e.g. [28][29][30][31][32][33][34][35][36][37] and references therein). In this paper we show that merger simulations of BH binaries of comparable masses in quasicircular orbits (as well as high-energy, head-on BH collisions) do, in fact, excite nonlinear modes in the ringdown stage. Second-order quasinormal modes. In BH perturbation theory, the GW strain and the Newman-Penrose scalar Ψ 4 produced by a BH merger at late times can be approximated by a linear combination of damped sinusoids (in addition to a subdominant power-law tail as well as retrograde QNMs, which we disregard here) [14][15][16][17],
rh (1)
where r is the (luminosity) distance from the source. The spin-2 spin-weighted spheroidal harmonics S m = S m (θ, φ, χω n m ) depend on the angular variables (θ, φ), on the complex QNM frequencies ω n m , and on the dimensionless spin χ = J/M 2 of the remnant BH [38]. This expression, found by solving the Teukolsky equation [39], is valid when the GW amplitude is small enough that one can linearize Einstein's equations in the Kerr background. At second order in the GW amplitude one finds similar equations for the second-order perturbations h (2) , now sourced by first-order quantities [40][41][42][43][44][45][46][47]. Let k be a generic mode, which can be either a first-order mode (k = k i = i m i n i ) or a higher-order mode. We will denote a second-order mode sourced by the first-order
From a waveform modeling point of view, the second-order modes are just damped sinusoids, like the first-order modes. Spin-weighted spherical harmonics, rather than spheroidal harmonics, are commonly used for waveform extraction in numerical relativity [48]. In our analysis, the index k will belong to a set
ing all the indices of the QNMs present in the m spin-2 spin-weighted spherical harmonic component. Then a ringdown waveform including both first-and second-order modes can be schematically written as rh (2)
where A k, m and φ k, m are the amplitude and phase of the kth (linear or nonlinear) mode found in the m spin-weighted spherical harmonic component. Note that A 1 mn, 2 m could be nonzero even if 1 = 2 because the spheroidal harmonic S 1 m is not necessarily orthogonal to the spherical harmonic Y 2 m , even if 1 = 2 [48].
Because second-order QNM frequencies are sourced by first-order modes, their frequencies, amplitudes and phases are expected to obey the relationships [42][43][44][45][46][47]
Second-order modes are a robust prediction of the perturbative expansion in general relativity. Other nonlinearities in the ringdown, such as the memory effect [49] or absorption-induced mode excitation [50], have previously been observed in simulations. However nonlinear QNMs have never been confidently identified until recently [51], with the exception of pioneering work by London et al. [45] using greedy fitting algorithms. Second-order modes in merger simulations. We have looked for the second-order modes in two sets of binary BH merger simulations. The first set consists of ultrarelativistic head-on collisions of equal-mass, nonspinning BHs with different boosts γ, similar to the sequences considered in Refs. [52,53]. In this one-parameter family of solutions the amplitude of the linear mode increases with the boost parameter γ, so the amplitude of the second-order modes is also a monotonic function of γ. Axial symmetry allows us to simulate this problem in two dimensions with GRChombo [54,55] by applying dimensional reduction [56][57][58], thus saving computational time and allowing for better accuracy relative to previous work [59]. As the quadratic modes are sourced by a product of two first-order modes, and quadratic contributions (proportional to Y i mi Y j mj ) overlap with Y i + j mi+mj , we will look for the i m i n i × j m j n j mode in the i + j m i + m j ringdown waveform [43,45,47]. For head-on collisions, we will be fitting rΨ 4 = r ḧ instead of h, and all of the reported amplitudes refer to rΨ 4 .
In this case the 200 mode dominates the ringdown of the nonspinning remnant [52], so we focus mainly on the 200 × 200 mode in the m = 40 waveform.
The second set of simulations consists of quasicircular mergers of binary BHs with different mass ratios from the publicly available SXS waveform catalog, simulated in (3 + 1)-dimensions with the spectral code SpEC [60]. Recent waveforms produced using Cauchy characteristic extraction [61][62][63] may improve the quality of our fits, but the relatively small set of publicly available waveforms does not adequately cover the relevant parameter space for our study. For quasicircular mergers the 220 and 330 modes are typically dominant (with their amplitudes depending on the mass ratio and spins of the binary), and we focus our search on (i) the 220 × 220 mode in the m = 44 waveform, and (ii) the 220 × 330 mode in the m = 55 waveform.
Identifying QNMs in a waveform can be challenging, partly because of their rapid decay; for nonspinning BHs, their quality factor is of order ∼ 3 [14]. The search for subdominant modes, which decay faster, requires some care. Even if their inclusion yields smaller fit residuals, consistency checks are crucial to avoid overfitting. In this work we fit the waveforms by a linear combination of damped sinusoids, as in Eq. ( 2), using a least-squares fitting algorithm. The amplitude and phase of each mode are always free fitting parameters, while the complex QNM frequencies are either free or fixed depending on the mode, as shown in Fig. 1 and explained below.
We first try to find the second-order modes without assuming knowledge of their QNM frequencies, as follows. We consider the QNM frequencies as free fitting parameters, and we fit the waveform with a different number of QNMs as we vary the starting time of the fit t start . If a fitted QNM returns a frequency that is consistent with a linear mode expected to exist in the waveform over a wide range of t start , we assume that the mode is there. We then fix the frequency of that QNM (as calculated in BH perturbation theory) in our fit and we add more QNMs with free frequencies to search for additional modes. We iterate until we do not see returned frequencies that are consistent with any linear modes. For head-on high-energy mergers, we find a combination of the modes 200, 400 . . . 10 00 in the m = 40 waveform due to numerical contamination between modes. For the SXS waveforms, we only confidently identify the 440 mode in the m = 44 multipole, and the 550 mode in the m = 55 multipole (out of all possible linear modes). With these first-order modes identified, we use a fitting model that consists of all such modes (with fixed frequencies) to search for additional higher-order modes by adding one more damped exponential with free frequency. As shown in the top-row panels of Fig. 1, when we vary t start relative to a reference time t 0 (defined to be the time of peak luminosity of the dominant m = 22 multipole), the free mode hovers around the expected second-order mode frequency (from left to right: ω 200×200 , ω 220×220 , or ω 220×330 , respectively). We do not expect the free mode frequency to (tstart -t0)/M (tstart -t0)/M We highlight in brighter colors the results for γ = 1.5 (left), q = 1.22 (the "SXS:BBH:0305" simulation, center) and q = 1.88 ("SXS:BBH:0403" simulation, right), while we plot the results for all other simulations in gray. Top row: search for the second-order mode frequency. We use a mode with a variable complex frequency in our fitting model to search for the expected second-order modes, and we use modes with fixed frequencies (black solid circles) to remove the contribution from linear modes when they are present. The color scale (top bar) represents different starting times of the fit. For quasicircular mergers, the labeled modes correspond to those of the remnant BH in the highlighted simulation. The location of the target mode for other simulations may be slightly different, because it depends on the remnant spin. Second row: fractional deviation |δω| of the fitted complex frequency with respect to the expected second-order mode. Third row: amplitude of the second-order mode when tstart is varied across a window of length T0 centered around the value of minimum |δω|, topt, and the second-order mode frequency is fixed to its expected value in the fitting model. Bottom row: same as the third row, but for the phase of the second-order mode.
converge exactly to the expected frequency due to numerical noise and contamination from other effects (such as additional nonlinearities) in the waveform, especially for modes that decay significantly faster than the dominant mode. In the Supplemental Material, we show through a controlled experiment that a free frequency hovering near the target mode is the expected behavior in the presence of (small) unaccounted additional modes. We also searched for the 200 × 400 mode in the head-on simulations. The results (which are not as clean as those for the 200 × 200 mode, because 200 × 400 is subdominant) are shown in the Supplemental Material. Having established the presence of nonlinear modes in the simulations, we now perform further checks to verify their physical nature.
We cannot exclude a priori that the new mode we found is in accidental agreement with the expected second-order QNM frequency. A nontrivial consistency test requires that, in addition to the frequency, the amplitude of the second-order modes should be consistent across different fitting ranges. To check this, we first look for the "optimal starting time," t opt , for which the fractional deviation between the fitted and expected complex frequencies, i.e.
, has a minimum. In the three cases of interest, = ω 200×200 , ω 220×220 or ω 220×330 , respectively. Then we assume that the mode exists, we fix the frequency to the expected value in our fitting model, and we check the consistency of the fitted amplitude. More explicitly, we check whether the recovered amplitude has an error smaller than 10% when t start varies within a window of length T 0 centered around t opt , where T 0 is the period of oscillation of the fundamental mode across all m multipoles (T 0 = T 200 for head-on mergers, and T 0 = T 220 for inspirals). We choose this value of T 0 because it is at least two times larger than the period of the second-order mode that we are searching for. This threshold is further justified in the Supplemental Material by studying the impact of the numerical noise in the simulations on the quantities of interest. Later times are excluded because the second-order mode falls below the numerical noise floor.
We find that all the waveforms we considered satisfy this requirement on the amplitude. We also checked that the amplitudes obtained from a model with free frequency are consistent with those where the frequency is fixed, albeit with larger fluctuations, as expected. Independently of the chosen t start , we use the convention A k, m ≡ A k, m (t peak ). In other words, we take into account the known exponential time decay by extrapolating the fitted amplitudes back towards the peak of the dominant multipole. Second-order amplitude dependence. As a more stringent check, we can verify whether the recovered second-order mode amplitudes follow the dependence predicted in Eq. (3b) across different simulations. For each simulation, we extract the second-order mode amplitudes by taking the mean of the amplitude within the T 0 starting time window mentioned above. We extract the first-order mode amplitudes after t start -t 0 = 25M , when nonlinearities and overtones have died out. We estimate the errors on the amplitudes as detailed in the Supplemental Material.
In Fig. 2 we plot the second-order mode amplitudes versus their first-order counterparts on a log-log plot. The data are consistent with a power-law dependence when the errors are taken into account. The slope of the fitted line for A 200×200 vs A 200 (in the head-on waveforms) and A 220×220 vs. A 220 (in the SXS waveforms) is found to be consistent with 2 within 1σ, as expected. Similarly, the slopes of the fitted lines for A 200×400 vs A 200 A 400 (for head-on waveforms) and A 220×330 vs A 220 A 330 (for SXS waveforms) are consistent with 1. Unsurprisingly, the 200 × 400 mode search results in head-on mergers are not as clean as the 200 × 200 results (see the Supplemental Material).
Because of numerical errors in the simulations, we can confidently identify the 220 × 220 mode only for SXS waveforms with mass ratio q ≤ 1.5. Since q varies over a small range, the amplitudes of the 220 mode inferred from different simulations are similar to each other, and the amplitude of the 220 × 220 mode does not vary much across different simulations. For this reason the data points are relatively close to each other, and the error on the slope is larger than in the other cases we considered. Phase consistency. Similar to the amplitude tests, we can check the consistency of our fits with the fitted phases of the second-order modes. As shown in the bottom row of Fig. 1, the fitted phases of the second-order modes vary by less than 10% × 2π within the T 0 window.
Moreover, as the second-order modes are sourced by two linear QNMs, the relationship in Eq. (3c) between the phases of the modes should hold, modulo (possibly) a constant phase difference that can only be computed by a Green's function calculation. In the right column of Fig. 2 we show that the phases extracted from the SXS simulations follow the expected relationship. In the Supplemental Material we show similar plots for head-on mergers. The error bars are larger, but the results are still consistent with expectations. Conclusions. We have shown that nonlinear QNMs are excited in simulations of comparable-mass BH binary mergers in quasicircular orbits, as well as in high-energy head-on BH collisions. The detectability of nonlinear QNMs may require next-generation detectors, and it will be addressed in future work. In any case, the presence of nonlinear modes demonstrates that nonlinearities must be taken into account in the modeling of GWs from binary BH mergers, and it suggests that they may play an important role during the violent merger phase. This has far-reaching consequences for our understanding of strong-field BH dynamics and for the observational BH spectroscopy program. Note added. While preparing this Letter, we learned that Mitman et al. conducted a similar study, whose results agree with ours [63]. Acknowledgments. We thank Yanbei Chen, Macarena Lagos, Sizheng Ma, Keefe Mitman and Leo Stein for discussions. The QNM frequencies used in this paper were computed using the qnm package [64]. We also made use of the sxs package [60] to access the SXS waveform catalog, and of the matplotlib [65], numpy
for the quasicircular mergers. Note that, in general, the coefficients may depend on the remnant BH spin. For quasicircular mergers, our fit of the 220 × 220 mode coefficient in the amplitude dependence equation is in excellent agreement with Ref. [63]. The coefficients of the 220 × 220 and 220 × 330 mode amplitudes are also in good order-ofmagnitude agreement with the estimates in Ref. [45]. Head-on merger simulations. We use the GR-Chombo codebase to simulate the head-on merger of two BHs with Lorentz boost factor γ (see [59] for details). The axial symmetry along the collision axis allows us to use dimensional reduction (see [56][57][58]) and simulate the system using only two dimensions, significantly decreasing the computational cost and improving the precision of the simulations. We use initial data similar to Ref. [53], consisting of a superposition of two boosted BHs. The constraint violations in the initial data are mitigated by using a large initial distance R init between the two BHs (typically R init = 1600 -2000 M 0 , where M 0 is the rest mass of the BHs). Selection of SXS waveforms. We use a conservative criterion to select the SXS simulations included in our study: (i) we require the free-mode search to yield a fitted frequency that is closer to the nonlinear mode than any other possible linear mode for a duration of T 0 ; (ii) when we include a second-order mode with fixed frequency in our fitting model, we require the fitted amplitude to be consistent within 10% in a fitting window of length T 0 around the time where the fitted frequency in the free mode search is closest to the expected one. We find that only nonprecessing systems with individual progenitor dimensionless spins χ 1 , χ 2 0.5 meet these criteria. For the 220 × 220 mode search, we only use simulations with mass ratio q ≤ 1.50, while for the 220 × 330 mode search we only use those with 1.25 ≤ q ≤ 2.00, because simulations outside of these ranges do not meet the above criteria. We include all simulations that meet the criteria, excluding those with almost overlapping simulation parameters. We find the nonlinear modes consistently for most of these simulations, with the exception of earlier simulations that do not give accurate higher m multipoles (e.g., SXS:BBH:0194), or simulations where the higher m multipoles do not settle down to a clean QNM decay before t -t 0 ∼ 35M (e.g., SXS:BBH:0605).
The simulation numbers used and their corresponding parameters are shown in Table I. We stress again that the selection criteria are conservative, and if we allow for a larger variation of the amplitude we could include additional simulations that are still consistent with the quadratic amplitude dependence, albeit with larger error bars. We expect the quality of the fits to improve for waveforms using Cauchy Characteristic Extraction [61][62][63], which are not yet publicly available. Optimal starting times of the fits. In the main text, we checked the consistency of our fits around an optimal starting time of the fits, defined as the starting time t opt, m that gives a fitted frequency closest to the expected quadratic mode in the m harmonic. For reproducibility, in Table II we list t opt, m /M for each of the head-on and quasicircular merger simulations we considered. Variation of the starting time window. The estimated total error on the amplitudes often has a significant contribution from the amplitude fluctuation error δA fluc . This is computed as follows: we first find the fit starting time that gives a frequency closest to the expected one in a free mode search, then we compute the standard deviation of the amplitude within a time window of length T 0 around this time. In Fig. 3 we show that while varying the length of the window affects the errors in the inferred slope of the amplitude dependence, the results are still consistent with prediction within 1σ unless the time window is very narrow, i.e. T 0 0.5T 200/220 .
If a quadratic mode is present in the waveform, a model containing such a mode should fit the waveform significantly better than one that only includes linear modes. The "goodness of fit" can be quantified by considering the mismatch between the best fit waveform ψ and the actual waveform Ψ,
where the scalar product is defined as
the asterisk denotes complex conjugation, and t end is the ending time of the fitting window. In Fig. 5 we show that the mismatch always decreases significantly when the model includes a quadratic mode, rather than only the other (linear) modes.
We estimate the errors on the amplitudes, which include contributions from the fluctuation of the amplitude within the T 0 window, the error between simulation resolutions, and the waveform extraction errors, as follows.
For the SXS simulations, we take into account the fluctuation (standard deviation) δA fluc of the amplitude within the T 0 window, the difference δA res of the amplitude extracted from the highest and second-highest resolutions of the SXS simulations, and the difference δA extrap of the amplitude extracted from the waveforms with extrapolation orders of N = 2 and N = 3. The total error is estimated as δA tot = δA 2 fluc + δA 2 res + δA 2 extrap . For the head-on merger simulations, we estimate the error as δA tot = δA 2 fluc + δA 2 res + δA 2 extract + δA 2 mass , where δA res is computed by Richardson extrapolation, δA extract is the variation of the amplitude when the waveform is extracted at the radii 900, 1000, . . . 1500M 0 , where M 0 is the total rest mass of the binary, and δA mass is propagated from the error of the estimated remnant mass.
We do not include δA mass for the SXS waveforms because this contribution is negligible.
In the right column of Fig. 2 of the main text we show that the extracted phases of the SXS simulations follow the expected relationship, Eq. (3c) of the main text. Figure 4 shows analogous plots for the 200 × 200 and 200 × 400 modes of the head-on mergers. The results are again consistent with the expected relationship, although the extracted phases are less accurate when compared to the SXS results. This is likely because the head-on simulation waveforms are strictly real (while the quasicircular SXS waveforms are complex), which makes it harder to extract the phase evolution of the waveform. The 200×400 mode of head-on mergers. As discussed in the main text, we find strong evidence for second-order modes. In particular the m = 55 multipole contains the 220 × 330 mode, which is sourced by two different first-order modes. By the same argument, we also expect to observe the 200×400 mode in the m = 60 multipole of head-on merger waveforms. As shown in Fig. 6, the free mode search for the 200 × 400 mode does not converge as well as that for the 200 × 200 mode, while the real parts of the fitted frequencies are arguably consistent with the expected mode. The extracted amplitudes and phases are also not as flat as those of the 200 × 200 mode, even varying by more than 10% within the T 0 window for some values of γ. Note that similar to linear modes like the 200 and 400 modes, the 200 × 200 mode is mixed into the m = 60 harmonic due to numerical contamination, and we have to include it in our fitting model along with the other linear modes (shown as black circles in Fig. 6). While higher-resolution simulations might be required to cleanly extract the 200 × 400 mode, if we assume that the 200×400 mode exists and extract its amplitude and phase, we find that the values are still consistent with Eqs. (3b) and (3c) of the main text, as shown in the bottom left panel of Fig. 3b and in the right panel of Fig. 4.
A good way to test our mode searching procedure is by constructing toy model waveforms with simulated noise that we can control. We consider three simulated scenarios: the first (toy model A) is a contamination-free idealized scenario; the second (toy model B) is the realistic scenario we expect to encounter in our simulations; the third (toy model C) serves to show the distinguishability of nonlinear modes from linear modes with similar frequencies.
We first test the ability of our procedure to identify modes in waveforms with noise. We construct toy model A with two injected modes plus numerical noise to mimic the m = 44 harmonic waveform of the SXS:BBH:0305 simulation:
+ A 220×220 e -i(ω220×220t+φ220×220) + n 0305,44 (6) The numerical noise in the m = 44 SXS:BBH:0305 waveform, n 0305,44 , is conservatively (over)estimated by computing the difference between the waveform with the highest resolution (Lev 6) with extrapolation order N = 2, and the one with second-highest resolution (Lev 5) with N = 3. Before taking the difference, we take care of the cumulative phase error of the waveforms by aligning them such that the mismatch in the ringdown phase is minimum. The amplitudes and phases of the injected 440 and 220 × 220 modes are taken from the best-fit values at the optimal starting time (i.e., when |δω| is a minimum in the center column of Fig. 1), so as to mimic SXS:BBH:0305 as close as possible. As seen in the left column of Fig. 7, when we fit the waveform with one mode fixed at the 440 mode frequency and another mode with free frequency, the free mode converges to the 220 × 220 frequency, as expected. Moreover, when we also fix the frequency of the 220 × 220 mode, the fitted amplitude and phase are consistent within 10% in the T 0 starting time window around the optimal fitting time.
Of course, the waveform is not only contaminated by noise, but also by other modes and by nonlinearities close to the peak of the waveform. To determine how this contamination would affect our results, we construct toy model B by adding the 441 overtone: Ψ B (t) = Ψ A (t) + A 441 e -i(ω441t+φ441) , (7) where we take A 441 = 5A 440 and φ 441 = φ 440 + π. We use this value for the phase because the first overtone has been found (at least in the m = 22 waveform) to be almost in anti-phase with the fundamental mode, but changing the phase does not affect our results to a significant degree.
As shown in the middle column of Fig. 7, the free mode hovers in the vicinity of the 220 × 220 frequency, while the amplitude of the fixed-mode fit varies by ∼ 10% across the T 0 time window. This shows that if other contaminations are present, even if they decay faster than the 220 × 220 mode, we should not expect the free mode search to converge perfectly, nor the fitted amplitude to be perfectly consistent. Incidentally, we could not confidently resolve the 441 overtone in our waveform, so we expect that there should be more contaminations (other than the 441 mode) that are masking it. Therefore the fits in Fig. 1 of the main text should perform at most as well (if not worse) than the results shown here for toy model B.
Finally, suppose that the 220 × 220 mode does not exist in the signal. Would our procedure return false positive results? This is a particularly important test for the following reason: it is well known that it is easier to accurately identify the real part of the frequency of a mode in a ringdown waveform than its imaginary part. By coincidence, the real parts of ω 441 and ω 220×220 are similar, so we need to make sure that what we are seeing is not the 441 mode, rather than the 220 × 220 mode. We do this by considering toy model C, where the waveform only consists of the 440 and 441 modes: Ψ C (t) = A 440 e -i(ω440t+φ440) + A 441 e -i(ω441t+φ441) + n 0305,44 . (8) The amplitude and phase of the 441 mode are the same as in toy model B. As shown in the right column of Fig. 7, the fitted frequencies for the free-mode search scatter around ω 441 , significantly further away from ω 220×220 than the results in Fig. 1. The fitted amplitude also varies by a factor of ∼ 1, much more than ∼ 10%. However, note that while the 220 × 220 mode is not present, the fitted phase varies by less than ∼ 10% × 2π, because the 220 × 220 and 441 modes have similar real frequencies. This shows that it is important to consider the consistency of both the amplitude and the phase when searching for the secondorder modes, which often have a real frequency close to other linear modes.
In conclusion, toy models A, B, and C confirm that our mode-searching procedure is robust against numerical noise and against contamination by faster-decaying modes. The criterion of an amplitude consistent within ∼ 10% in a time window of T 0 is also adequate to rule out false positives, while accounting for the variation caused by contamination.