Extracting linear and nonlinear quasinormal modes from black hole merger simulations

I. INTRODUCTION

The pioneering work on the linear perturbations of the Schwarzschild spacetime by Regge and Wheeler [2], Zerilli [3,4] and Vishveshwara [5] revealed that the response of a black hole (BH) to an incoming pulse of radiation is characterized by a superposition of damped exponentials with discrete frequencies and damping times, now commonly known as the "ringdown" by analogy with the dying tones of a vibrating bell. The damping occurs because BH spacetimes absorb gravitational waves (GWs) at the horizon and emit radiation at spatial infinity: the system is dissipative, hence the name "quasinormal" modes (QNMs), as opposed to the "normal" modes of self-adjoint physical systems [6][7][8][9][10].

The correspondence between BH spectra and atomic spectra was soon clear to early researchers who computed QNM spectra and clarified their physical interpretation [11][12][13][14][15][16][17][18][19]. The separation of angular variables results in multipolar components characterized by a discrete set of angular indices ðl; mÞ. For fixed value of these indices, the spectrum contains a tower of discrete QNM frequencies ω l;m;n that can be sorted by the magnitude of their imaginary part (the so-called "overtones"). Typically n ¼ 0 denotes the fundamental mode, while increasing values of n correspond to larger imaginary parts and shorter damping times. Teukolsky proved that the perturbation equations for rotating (Kerr) BHs are also separable [20][21][22][23]. Detweiler pointed out that the fundamental QNM frequency of a Kerr BH (the one with the smallest imaginary part and longest damping time) depends only on its mass and spin [24], so-at least conceptually-the relation can be inverted to identify the Kerr BH parameters from a knowledge of the frequency and damping time. In fact, all complex frequencies are fixed once we know the mass and spin of the remnant, so the observation of additional resonant frequencies can be used to identify Kerr BHs. This is the GW equivalent of atomic spectroscopy [25].

Within linearized BH perturbation theory, Leaver [19] proved that each multipolar component of the waveform at intermediate times-after the "prompt response", and before the onset of power-law tails-is described by a superposition of QNMs, h lm ðtÞ ≡ X n A l;m;n e -i½ω l;m;n ðt-t start Þþϕ l;m;n : ð1Þ

Here, t start is an arbitrary starting time. In linearized GR, the complex Kerr QNM frequencies ω l;m;n depend only on the remnant BH mass M rem and dimensionless spin χ rem , but not on the nature of the perturbation, and are known to very high accuracy [26,27]. On the contrary, the QNM amplitudes A l;m;n and phases ϕ l;m;n depend on the astrophysical process causing the perturbation.

The merger of two comparable-mass BHs was identified early on as one of the most promising LIGO-Virgo sources [28], but predicting which QNMs are most excited as a result of the merger was essentially a matter of guesswork before the first numerical BH merger simulations. This is important because practical attempts to implement the spectroscopy program in GW data analysis rely on the measurability of the frequency and damping time of the fundamental mode, and possibly of other modes [26,[29][30][31]), but only the frequencies of modes with sufficiently large amplitude are measurable. Green's function techniques imply that the QNM amplitudes A l;m;n can be factorized as a product of complex "excitation factors" B l;m;n that depend only on the remnant's mass and spin and complex-valued, initial-data dependent integrals I l;m;n [19,[32][33][34][35][36][37][38]. Before the first numerical BH merger simulations, heuristic arguments suggested that comparable-mass BH mergers may have ringdown signal-to-noise ratio (SNR) roughly comparable to the inspiral SNR [39], while other astrophysical processes would be much less efficient at exciting QNMs [40]. However, early work trying to quantify the measurability of different multipoles and overtones had to rely on educated guesses [26].

Our understanding of ringdown excitation improved after the 2005 numerical relativity (NR) breakthrough [41][42][43]. Fits of NR simulations showed that the radiation from a BBH merger is dominated by the l; jmj ¼ 2, 2 spherical harmonic multipole, while higher multipoles are subdominant [44,45]. Unbiased mass and spin estimates require the combination of different multipolar components [46] or the inclusion of overtones [47,48]. Since multiple modes will always be excited to some extent, we must first understand which combination of modes will dominate the signal [26]. Are we going to observe a combination of low-n modes for different multipoles, or are higher overtones of the l; m ¼ 2, 2 component dominant? Can we measure the frequencies and damping times accurately enough to resolve the modes? The answers to these questions depend on (i) the properties of the merger remnant progenitors, and (ii) the sensitivity of the detectors [45][46][47][49][50][51][52][53][54][55][56][57][58][59][60][61][62][63]. Initial attempts at fitting numerical simulations with a superposition of QNMs argued that it is possible to fit the numerical waveforms with a superposition of several overtones within linear perturbation theory [48,51,64], possibly including counterrotating modes and higher multipoles [59][60][61]63,[65][66][67][68][69][70][71][72][73][74][75]. More recent work (confirming early findings by London et al. [51]) showed that nonlinearities are important: quadratic modes are present in NR merger simulations of quasicircular and head-on mergers [76][77][78][79]. This caused a recent surge of activity in the analytical modeling of nonlinear QNMs [38,[80][81][82][83][84][85] (see e.g., [86][87][88][89][90][91][92][93][94] for earlier work, and [63,95] for a discussion of other nonlinearities that could affect the ringdown).

In this paper we continue the "agnostic" BH spectroscopy program of Ref. [96]. We develop a systematic fitting algorithm to extract the QNM content of NR waveforms. Our goal is, essentially, to build the most accurate ringdown model that avoids overfitting. We first identify which modes can be reliably found in NR data. Then we produce a catalog of their amplitudes, phases and (possibly) starting times. We obtain polynomial "hyperfits" of the resulting mode amplitudes as a function of the parameters of the merger progenitors, and we look for universal relations in the amplitude ratios of different modes.

The plan of the paper is as follows. In Sec. II we review the classification of all possible QNMs that can be found from fits of a NR simulation. In Sec. III we describe the algorithm underlying our QNM extraction procedure and show illustrative examples of its application to specific waveforms in the SXS catalog. This rather lengthy and technical section is the core of this paper. It is crucial to understand the assumptions underlying our QNM amplitude estimates, but it can be skipped by readers who want to focus on the physical implications of our work. In Sec. IV we apply the algorithm to the dominant multipolar components of simulations from the whole SXS catalog, and we use the results of this large-scale fitting campaign to produce "hyperfits" of the QNM ampitudes for generic BH binary systems. In Sec. V we study the relation between the amplitudes of various modes, including (i) the ratio between nonlinear modes and their "parent" linear modes, (ii) the ratio between the amplitudes of the first overtones and the fundamental modes, and (iii) the ratios between retrograde and prograde modes. In Sec. VI we present some considerations on the thorny issue of defining a starting time for the ringdown for data analysis purposes. In Sec. VII we present some conclusions and directions for future research.

To improve readability, consistency checks and more technical material are relegated to several appendixes. In Appendix A we clarify our definition of prograde and retrograde QNMs. In Appendix B we complement the discussion of Sec. III with some technical details about our fitting routine. In Appendix C we develop a toy model to test the accuracy of the fitting algorithm. In Appendixes D and E we "stress test" our fitting routines using this toy model, and show that they do indeed give consistent results. In Appendix F we study spherical-spheroidal mode mixing, and show that the fitting routines give results consistent with theoretical expectations. The current SXS catalog does not, in general, make use of Cauchy characteristic extraction (CCE), but some simulations using CCE are publicly available. In Appendix G we apply our fitting algorithm to the limited set of public CCE simulations, and show that the use of better wave extraction techniques gets rid of "spurious" QNMs. In Appendix H we investigate the effect of using different stability criteria to determine whether a mode is robustly present in the ringdown. Finally, in Appendix I we discuss some outliers in the large catalog of extracted QNM amplitudes and their possible origin.

II. QUASINORMAL MODE CLASSIFICATION

General relativity is a nonlinear theory of gravity. The evolution of the full nonlinear equations for binary black holes (BBHs) poses theoretical challenges and is computationally expensive. However, if our goal is to study the latetime postmerger behavior, the spacetime asymptotically approaches the Kerr solution, and we can linearize the Einstein equations on a Kerr metric background.

At late enough times, neglecting contributions from the back-scattering of GWs, the postmerger waveform can be modeled as a linear combination of QNMs with complex frequencies ωk

where A k , ϕ k , ω r;k and ω i;k are real. In general, ω i;k ≤ 0. Each mode is labeled by a mode multi-index k, which has the form k ¼ l; m; n for linear modes, and the form k ¼ l 1 ; m 1 ; n 1 × l 2 ; m 2 ; n 2 for quadratic modes (we will neglect contributions beyond quadratic order). The set of multi-indices K includes all the modes k that could be present in the waveform. We will also define the mode k ¼ constant, where ω r;constant ¼ ω i;constant ¼ 0, denoting a "pseudo-QNM" contribution that acts as a complex constant in time, A constant e iϕ constant . This is useful to model waveforms containing a constant component.

In this work, we are interested in recovering the list of all the modes K that are present in NR waveforms, and the corresponding values of A k and ϕ k . To do this, we will employ two types of fitting models: (i) a "frequencyagnostic" model with N f free QNMs, where A k , ϕ k , ω r;k and ω i;k of each mode are left as free parameters

or (ii) a model with only A k and ϕ k free, but with ω r;k and ω i;k fixed at the values corresponding to specific modes

The free parameters are underlined in Eqs. ( 3) and (4). In model (4), the fixed values of ω r;k and ω i;k for a particular mode k can be computed through standard methods in BH perturbation theory given the remnant mass M rem and dimensionless spin χ rem . We will call these fixed frequencies the Kerr QNMs frequencies. While the QNM contribution to the strain decays exponentially over time, the values of A k and ϕ k , as defined in Eqs. ( 3) and ( 4), are constant. Unless otherwise specified, we will always measure the time t starting from t peak , the time corresponding to the peak strain for the l; m ¼ 2, 2 multipole of any given simulation. This is equivalent to defining t peak ≡ 0 when using Eqs. ( 3) and ( 4), or to making a time shift by substituting t → tt peak . By doing this, all of the A k and ϕ k best-fit values are those "extrapolated back to t peak " regardless of the starting time t 0 (≥ t peak ) of our fitting window. In other words, we define the amplitude to be A k instead of A k e ω i;k ðt 0 -t peak Þ , and the phase to be

The perturbation of a Kerr BH spacetime can be decomposed in terms of spin-weighted spheroidal harmonics S lm , which are labeled by polar and azimuthal numbers l and m, similar to spherical harmonics Y lm [21,97]. For a nonrotating (Schwarzschild) BH, the spin-weighted spheroidal harmonics S lm reduce to the spin-weighted spherical harmonics Y lm , but the two are different in general. In particular, S l 1 m 1 is not orthogonal to Y l 2 m 2 as long as

In NR, GWs are usually extracted using a spin-weighted spherical harmonic decomposition [98]. While we still expect to see the l; m; n Kerr modes in the l; m harmonic because S l 1 m 1 overlaps the most with Y l 2 m 2 if l 1 ; m 1 ¼ l 2 ; m 2 , Kerr modes belonging to a different set of l 0 ; m 0 values can get mixed into the l; m harmonic by various means.

In a given l; m spherical harmonic, we expect to see a rich spectrum of different QNMs belonging to the following categories [96]:

(i) "Natural" linear modes: the l; m; 0 fundamental mode and the l; m; n overtones with n ≥ 1. These modes appear naturally in the l; m harmonic, in the sense that their l; m number are the same as the underlying harmonic. (ii) Spherical-spheroidal mixing modes: modes with indices l 0 ; m; n, l 0 ≠ l, mixed into the l; m harmonic because S lm is not orthogonal to Y lm [44,76,96,97,99,100]. (iii) Quadratic modes: These modes arise due to nonlinear coupling between two linear modes, and are labeled as l 1 ; m 1 ; n 1 × l 2 ; m 2 ; n 2 . They are a direct consequence of the nonlinearity of Einstein's field equations, and their frequencies can be computed with second-order BH perturbation theory [51,89,90],

The amplitude of a quadratic mode in a certain l; m harmonic depends on the Clebsch-Gordan coefficients [38]. The l 1 ; m 1 ; n 1 × l 2 ; m [38,89,90,92]. (iv) Retrograde modes: We define these to be the modes that counterrotate in space compared to the orbital angular momentum of the BBH. This is different from the more common definition in the literature (see e.g., Refs. [63,76]), which defines them to be those modes that are counterrotating with respect to the spin of the remnant BH. We clarify the rationale behind this definition in Fig. 12 and Appendix A.

Each corotating mode has a retrograde counterpart that we will label with the prefix "r": for example, the r2, 2, 0 mode is the retrograde counterpart of the 2,2,0 mode. As they have a frequency ω r of opposite sign compared to their prograde counterpart, they are counterrotating phasors in the complex plane [70,71]. Note that modes with m < 0 might have a negative frequency ω r , but they might not be retrograde in the sense of being counterrotating in space compared to the BBH orbital angular momentum. To identify the rotational direction, we multiply e -iω r t by e imϕ , where ϕ is the azimuthal coordinate. The mode is only retrograde in our definition if sgnðmÞ ≠ sgnðω r Þ (see Ref. [63] for a discussion of this point). (v) Recoil modes: These are either linear modes with l 0 ≠ l or quadratic modes with l 1 þ l 2 ≠ l. They exist in the "wrong" l; m harmonic because of imperfections in the choice of the Bondi-van der Burg-Metzner-Sachs (BMS) frame in NR simulations [76,101]. If the modes were extracted by CCE and the BMS frame was fixed appropriately, these modes would be (at least in principle) irrelevant [63]. Note that the above categories are not mutually exclusive: for example, it is possible to find quadratic recoil modes.

III. QUASINORMAL MODE EXTRACTION PROCEDURE

The detection of QNMs is a nontrivial task, due to their exponentially decaying nature. The difficulty arises both when detecting QNMs in real GW data and when trying to fit NR simulations, which are noise-free (modulo numerical noise). On one hand, this has stirred some controversy on whether we have already detected multiple overtones in a GW event [102][103][104][105][106][107][108]. On the other hand, it led to a better understanding of the physical excitation of higher overtones and of our (in)ability to extract them from NR simulations [48,96,109,110]. In this work we contribute to the ongoing discussion by pursuing the following goals:

(1) Mode finding. We identify the largest set of Kerr modes that are likely to be present in the postmerger waveform.

(2) Amplitude/phase extraction. We extract the amplitude and phase of each of the Kerr QNMs thus found. (3) Ringdown starting time determination. We determine the time at which each Kerr QNM starts to become stable. For such a procedure to work, it is clear that we must set well-defined (and inevitably, somewhat arbitrary) criteria to decide whether a Kerr mode is present in the waveform, or to determine when the mode is "stable". To do so in a quantifiable manner, we must specify some tolerance values for our QNM extraction procedure.

The three goals stated above are not mutually independent. For example, missing a certain Kerr mode in the mode-finding stage can bias the amplitudes and phases of the other modes, because the missed Kerr mode is not included in our fitting model.

A schematic representation of the full procedure is illustrated in Fig. 1 for the SXS:BBH:0305 simulation previously considered in Refs. [48,96]. The method is based on the following steps, which are explained in detail in the following subsections:

(A) Prepare postmerger waveform. We first select a gravitational waveform corresponding to a specific numerical simulation (e.g., SXS:BBH:0305) and a specific multipole (e.g., l; m ¼ 2, 2). We disregard the waveform data before t peak , defined to be the time at which jh 2;2 j attains its maximum value. This early segment of the waveform corresponds to the premerger signal.

(B) Agnostic frequency recovery. We identify potential modes in the waveform in a frequency-agnostic manner:

(B.1) We fit the waveform with N f free QNMs, with all four real parameters A; ϕ; ω r ; ω i left as free parameters to be fitted, as in model (4). For each N f we perform fits with a varying starting time t 0 to estimate the evolution of the QNM frequency content of the waveform with time.

(B.2) We identify potential modes by checking whether the fitted values of ω r and ω i are consistent with the expected Kerr frequencies, as computed in BH perturbation theory.

(C) Amplitude/phase stability test. We assume that all the potential modes are present in the waveform and test their stability in time:

(C.1) We fit the waveform with a fitting model consisting of all the potential modes fixed at their corresponding ω r and ω i , but with A and ϕ as free parameters to be fitted.

(C.2) We check whether the fitted amplitudes and phases of the modes are self-consistent (stable) across a specified time window. We iteratively remove the modes that fail the test.

(C.3) We repeat the check until all the modes pass the test.

(D) Select the best instance of stable modes. We repeat steps (B) and (C) for a range of N f and select the instance where we extract the largest number of stable modes. Now we have an exhaustive list of the modes that are present in the postmerger waveform and have stable complex amplitudes.

(E) Amplitude/phase extraction. We extract the amplitude and phase of each mode within a time window over which the mode is present and its amplitude is stable.

In some steps of the procedure we will define and use certain tolerance criteria to test whether a mode is present in the waveform. For the reader's convenience, the quantities defining these criteria are listed in Table I.

A. Prepare postmerger waveform

In this work, we will consider publicly available waveforms from two catalogs; the SXS waveform catalog [98] and the Ext-CCE catalog. The SXS catalog contains NR simulations of BBH mergers extracted at null infinity through extrapolation, while the Ext-CCE catalog benefits from improvements such as center-of-mass corrections, improved Cauchy-characteristic evolution, the inclusion of gravitational memory using BMS balance laws, and improved wave extraction [111][112][113][114]. The SXS waveforms are labeled as SXS:BBH:XXXX, and the Ext-CCE waveforms are labeled as SXS:BBH_ExtCCE:XXXX, where XXXX is a four-digit number. Each number corresponds to a specific numerical simulation of a BBH merger with different sets of intrinsic parameters (but, in general, simulations with the same label XXXX do not refer to simulations with the same intrinsic parameters in the SXS and Ext-CCE catalogs). The waveforms available are decomposed in spin-weighted spherical harmonics, so each number corresponds to a set of waveforms decomposed in different multipolar components labeled by l and m.

For a particular simulation, we will separately fit each multipolar component. Once we have chosen the simulation and the multipolar indices, we determine the peak merger time t peak where the absolute value of the strain of the dominant l; m ¼ 2, 2 multipole has a maximum. In what follows, we will be fitting the waveform within a time window that starts no earlier than t peak , so the waveform at t < t peak will not be relevant for our analysis. Certain harmonics in some simulations become noisy at late times, so we always trim the waveform at tt peak ≳ 120M. We find that the late time cutoff of the fitting window does not significantly affect the results of our analysis, because the modes decay exponentially in time.

The SXS catalog lists waveform data for many multipolar components l; m (up to l ∼ 8). However, in our analysis, we will only fit the multipolar component waveforms that have a sufficiently large peak strain: jh peak;l;m j > γjh peak;2;2 j, where γ ¼ 0.02 or 0.001. This is because waveforms with a low strain are more affected by contaminations due to mode-mixing and numerical noise, and they are less relevant for data analysis anyway. The quality of the higher l; m components is typically better for more recent simulations. Therefore, we use γ ¼ γ low ¼ 0.001 if the simulation number is SXS:BBH:XXXX > 0305 and l ¼ m ≤ 6. Otherwise, we use γ ¼ γ hi ¼ 0.02. We reduce the tolerance only for the l ¼ m ≤ 6 components because they are often more important, and we do not want to include too many other l; m components due to computational costs. For some of the next steps of the mode extraction procedure, we will have to compute the theoretically predicted QNM frequencies for specific values of l; m; n. These are easily computed [8,27] from the mass M rem and spin χ rem of the remnant BH, which are listed in the publicly available metadata file of each simulation. We will use the qnm package [115] to compute the frequencies within our procedure.

Agnostic frequency fit

We begin by fitting the waveform using N f free QNMs and all four real parameters A; ϕ; ω r , and ω i as free parameters to be fitted, as given by model (3). To capture the evolution of the free QNM frequencies over time, we vary the starting time of the fitting window from ðt 0t peak Þ=M ¼ 0 to 50 with a time step of Δt ¼ 0.1M, while keeping the ending time of the window fixed at t endt peak ¼ 120M. We consider up to N f ¼ 10 free modes, which corresponds to 40 real fitting parameters in model (3). To reduce the computational cost of the high-dimensional fit, we employ the JaxFit package [116], which implements a trust region method for nonlinear least-squares fitting. The implementation is similar to that in the curvefit method of the optimize subpackage in SciPy, but with auto-differentiation and hardware acceleration implemented with Jax [117]. The technical details of the fitting algorithm and the relevant hyperparameters are described in Appendix B. In Appendix D, we test the fitting algorithm and show that (with some caveats) it works well enough for our purposes.

In Fig. 2 we show the results of the frequency-agnostic fit for N f ∈ 5; 6; …; 10 using the l; m ¼ 2, 2 multipole of the SXS:BBH:0305 waveform. By visually inspecting the evolution of the frequencies, we can already identify some of the Kerr modes that have been captured by the agnostic fit. Notably, the 2,2,0, 2,2,1, 3,2,0, and constant modes are consistently approached by at least one free QNM for any value of N f we consider, strongly hinting at the presence of these Kerr modes in the waveform. In the following subsection, we will employ more precise criteria based on the results of the agnostic fit to determine the presence of a Kerr mode.

By using as many as N f ¼ 10 free QNMs we could be overfitting the waveform, in the sense that some of the QNMs might be fitting away signal components unrelated to the ringdown (such as the prompt response, GW memory, backscattering tails, or numerical noise). Indeed, as shown in Appendix D, with N f ≳ 7 the fit mismatch can be orders of magnitude below the estimated numerical noise floor. At this step of the mode extraction procedure, our main goal is to exhaustively identify all of the potential QNMs that may be present in the waveform; we do our best not to miss any significant QNMs that have relatively higher amplitudes and are longer-lived, at the cost of possibly including spurious ones, so we choose to overfit.

As shown in the N f ¼ 10 panel of Fig. 2, there could be as many as nine potential QNMs in the waveform, at least for the l; m ¼ 2, 2 multipole of the SXS:BBH:0305 waveform. Therefore using a value of N f as large as 10 might be necessary, in general, if we want to do our best at extracting QNMs. By inspecting the different panels, the N f ¼ 10 results are arguably the cleanest, with the free QNMs converging better to the known Kerr frequencies. Also, while other non-QNM contaminations could be present in the waveform, this does not necessarily mean that they would prevent the extraction of QNMs, even if their amplitudes are comparable to the QNMs. For example, Gaussian noise will not fatally impede a least-squares fitting algorithm, and the memory effect might be approximately mitigated by fitting away a constant at late times. Moreover, if a free QNM is fitting away a non-QNM component, it should not linger around the same frequency in the complex plane over time; conversely, if a free QNM lingers around the same Kerr QNM frequency over time, the component of the waveform it is fitting is then (by definition) a Kerr QNM.

Even if this stage of the fitting procedure were to pick up any spurious Kerr QNMs, at a later step we will use a more stringent criterion to rule them out.

Potential mode identification

We now examine the evolution of free QNMs and check whether the results indicate the presence of certain Kerr QNMs in the waveform.

We first need to decide, a priori, what Kerr modes are relevant to our analysis, or what modes we expect to find. For a strictly agnostic analysis, we could in principle deem all of the linear and quadratic modes with all possible combinations of indices (l, m and n) to be relevant, regardless of the l; m harmonic of the waveform we are considering. However, we expect that most of these modes will be negligible, especially when they are recoil modes. Therefore, while we deem all natural modes and sphericalspheroidal mixing modes relevant (including the retrograde modes within these categories), we only include recoil modes such that their natural-host component l; m has a peak strain jh peak;l;m j > γjh peak;2;2 j. This is the same requirement used for deciding which waveform components to fit, as explained in Sec. III A. In this way we ensure that we include the recoil modes that could potentially be present in the multipolar component in question because they could have a significant amplitude in their natural-host component.

Given a list of relevant Kerr QNM frequencies, we can determine their proximity to the free QNMs to check whether some of them are consistent with the agnostic fit. We do this by drawing an ellipse on the complex plane around each relevant mode frequency and checking whether any free QNM stays within the ellipse consistently FIG. 2. Agnostic fit results of the complex QNM frequencies for the l; m ¼ 2, 2 multipole of the SXS:BBH:0305 waveform. Each panel corresponds to a different number N f of free QNMs (labeled in the top-left corner) used in the fitting model (3). Dots of a particular color correspond to the evolution of a particular free mode over time, i.e. when we vary the starting time of the fitting window from ðt 0t peak Þ=M ¼ 0 to 50. Relevant Kerr modes are shown and labeled as gray circles, with those consistent with the agnostic fits shown in black. The gray ellipses around each circle are the tolerance zones of each mode. The gray vertical-dashed line divides the complex plane into a positive-frequency region (ω r > 0) and a negative-frequency one (ω r < 0). We show the region ω i > 0 in gray as it corresponds to exponentially growing modes, which should be absent.

(i.e., at least p ag % of the time, where "ag" stands for "agnostic") in any continuous time window t 0 ∈ ðt; t þ τ ag Þ. For example, if we choose τ ag ¼ 10M and p ag % ¼ 95%, and a free QNM stays within the ellipse around the complex frequency ω2;2;0 for at least 95% of the time between, say, t 0t peak ∈ ð16.7; 26.7ÞM, then we say that the 2,2,0 mode has passed the frequency-agnostic test. We choose the 95% criterion to avoid missing certain modes because of short fluctuations in the fitting results.

When we fit across different values of t 0 , we use the results of the previous time step (t 0 -Δt 0 ) as an initial guess for the fitting algorithm. In Fig. 2, we label the corresponding frequencies (a "stream" of free QNM frequencies connected by iterative initial guesses) across different t 0 with the same color. However, at a specific fit starting time t 0 , as seen in model ( 3), all free QNMs have the same form and contribute equally, meaning that we could very well shuffle the colors of different points corresponding to the same t 0 : when analyzing the results, at a specific t 0 we do not care about the color of each free mode (nor about the colors of the modes in a neighboring time step). Therefore, in the proximity test for a relevant Kerr mode we do not require that the mode has to be approached by a free mode of the same color. In other words, we do not care whether the mode "switched color" when it is within the ellipse. In fact, this "color switching" behavior is frequently observed for overtones, as shown in Fig. 2, where the 2,2,1 Kerr mode is often approached by modes with different colors.

We now have a list of potential Kerr modes, which will be subject to further screening in the next steps.

C. Amplitude/phase stability test

The frequency-agnostic test presented above checks whether a Kerr QNM frequency ω ¼ ω r þ iω i is consistently present, but a QNM also consists of two other parameters, A and ϕ. A QNM is defined to be given by Eq. ( 2) with A; ϕ; ω r , and ω i all being constant in time. A non-QNM component could, in principle, mimic a QNM by giving the same fitted frequency, but with time-evolving values of A and ϕ. Thus, an additional test is required to check whether A and ϕ are stable (i.e., constant) in time.

At this point, there are two routes we can take. The first is to examine the fitted A and ϕ of the free QNMs in the frequency-agnostic fits and check for stability, and the second is to do a frequency-fixed fit-as in model ( 4)-with frequencies chosen to match the potential Kerr modes identified by the frequency-agnostic test. We will undertake the second route for two reasons.

The first reason is that the frequency-agnostic fits could contain more free modes than the potential Kerr modes. The values of A and ϕ for a free QNM that picked up a Kerr QNM might be affected by other free QNMs that were "greedily" overfitting non-QNM components. One might wonder whether ω r and ω i of a free QNM approaching a Kerr mode could also be affected by the "greedy" QNMs. In fact, as seen in Fig. 2, increasing N f does not cause ω r nor ω i to deviate from a potential Kerr mode: the 2; 2; 0; 2; 2; 1; 3; 2; 0; 2; -2, 0 and constant modes are always approached by a free QNM, even if some other free QNM is not picking up any Kerr mode (see e.g., the brown mode in the N f ¼ 10 panel). This holds rather generally for all of the waveforms we examined.

The second reason is that in model ( 4) we fix ω r and ω i , reducing the number of free parameters, and therefore making it easier to examine the behavior of A and ϕ of a Kerr mode.

We fit the waveform with model ( 4), fixing ω r;k ; ω i;k with k ∈ K, where K is the set of potential Kerr modes, again over a wide range of starting times t 0 . We then test whether A and ϕ are both stable (i.e., they fluctuate less than a certain tolerance) within a time window

We quantify the fluctuation by a weighted quadrature sum of the fluctuations δA and δϕ,

where δA k and δϕ k are defined to be the p stable -percentile range of A k and ϕ k within the time window, and

< ϵ stable for some t, then we conclude that the potential mode is stable, and it is likely to be present in the waveform. If a subset of modes k ⊆ K is unstable, then we eliminate the mode with the highest min t Δ k;stable ðtÞ and redo the frequency-fixed fits with the remaining modes. The A 0 k term encourages the algorithm to pick up subdominant modes with low amplitudes A k ≲ A tol . Effectively, it penalizes them less for a fluctuating relative amplitude, while still screening them based on the fluctuation of their phase.

We choose to eliminate only one mode at a time because a spurious Kerr mode could interfere with a robust Kerr mode and destabilize it, especially if the two modes have similar frequencies (e.g., 3,2,0 and 3,3,0). By eliminating the most unstable mode first, which is likely the spurious mode, the robust mode will likely become stable in the next iteration and it will not be eliminated. As the natural fundamental mode of any given multipolar component must be present in the waveform, we will never eliminate it with this test, even if its amplitude fluctuates more than the tolerance. We iterate the procedure until all of the modes are stable. At this point we have a list of modes that are likely to be present in the waveform.

D. Finalizing stable modes

Steps (B) and (C) (agnostic frequency recovery and amplitude/phase stability tests) of the procedure are done in parallel for different values of N f . We choose to do multiple runs with different values of N f (instead of doing a single run with large N f , say N f ¼ 10) to make sure that we will not miss certain modes because of overfitting the waveform with too many free QNMs. Intuitively, there should be an optimal N f for a given waveform such that including fewer modes will cause underfitting, and including more causes overfitting. Therefore, across all the results with different N f , we select the instance that returns the longest list of potential modes to be our final result. In practice, we find that we rarely find fewer modes when we increase N f (at least for 5 ≤ N f ≤ 10), meaning that overfitting the waveform seldom causes us to miss a mode. Nevertheless, we still include the lower N f results to be safe, and because they are less costly in terms of computational resources.

There are alternatives to the above strategy. For example, one could have chosen instead the union of all the stable modes returned by the stability test for all N f . However, the correct extraction of a certain mode could depend also on the extraction of the other modes. For example, if for a certain N f we identify the Kerr modes f2; 2; 0j3; 3; 0g, and for a larger N f we find the set f2; 2; 0j2; 2; 1j3; 2; 0g, the union of the two will be f2; 2; 0j2; 2; 1j3; 2; 0j3; 3; 0g, but it might very well be the case that the lower N f run was mistaking the 3,2,0 mode as its neighbor 3,3,0 because it missed the 2,2,1 mode. Our strategy is more conservative and avoids this problem. Another alternative would be to do the frequency-agnostic fits for all different N f , identify the potential modes in each run, take the union between all runs, and do one single stability test, iteratively removing unstable modes from the joint list of QNMs. This raises concerns similar to those highlighted for the previous alternative strategy, so we do not use this method.

This step concludes the "mode finding" portion of the procedure. In principle, the ringdown waveform could contain infinitely many QNMs with arbitrary mode numbers. The main goal of our algorithm is to identify as many modes as possible, trying our best not to miss significant modes (i.e., those with a higher amplitude or those that are longer living) without including spurious modes that should have been undetectable. Indeed, our algorithm could miss certain modes, especially those with a high overtone number n; see Appendix E.

E. Amplitude/phase extraction

Given a set of Kerr modes K that are likely to be present in the waveform, we would like to extract A k ; ϕ k for all k ∈ K, so that we can specify completely the QNM content in the waveform. This is a nontrivial task because A k and ϕ k depend on the chosen fitting window. As explained previously, the ending time t end of the fitting window does not sensitively affect the fit. However, as can be seen in Fig. 3, the starting time t 0 does, so we should write the extracted amplitude and phase as A k ðt 0 Þ and ϕ k ðt 0 Þ. If the fit starts too early (t 0 is too small), the results will be systematically biased because the GW response has not stabilized into constant-amplitude/phase QNMs due to a number of reasons, including the time-evolving background, delayed build-up of QNM amplitudes, prompt radiation, and nonlinear effects [96,110]. If the fit starts too late, fast-decaying or small-amplitude modes can be missed because they dive beneath the noise floor.

Given these caveats, a simple way ahead is to extract A k and ϕ k at times such that the mode k is most stable. We can find, for each mode, the time t ¼ t best that minimizes the FIG. 3. Amplitude (top row) and phase (bottom row) stability test for the l; m ¼ 2, 2 multipole of the SXS:BBH:0305 waveform. In the leftmost panel we fit the waveform with all of the potential modes identified in the frequency-agnostic fit (see legend at the bottom). Moving from the left to the right panels, we iteratively remove the Kerr mode whose amplitude and phase fluctuate the most over time, until (in the rightmost panel) all modes are stable within the tolerance specified. quantity Δ k;stable ðt; τ stable Þ defined in Eq. ( 6). Then, we can extract A k and ϕ k as the median values within the time window t 0 ∈ ðt best ; t best þ τ stable Þ. We can also quantify the fluctuation of the mode, e.g., by the center 95th percentile of A k and ϕ k within the window (we do not use the maximum and minimum values within the window to quantify the fluctuation because we want to leave some margin for large sudden fluctuations in A k and ϕ k due to nonconvergence, which rarely happens). We show an example of this method in Fig. 4. Note that while we extract the amplitude and phases at times later than t peak , we always define their values to be those extrapolated back to t peak .

We now know which Kerr modes are likely to be present in the waveform, and with what amplitude and phase. However, as we have shown explicitly, QNMs are only a valid description of the postmerger waveform (modulo memory effects, noise, and backscattering) at late times. To fully specify our model, we must identify an optimal starting time. We can define the starting time of a particular mode to be the time at which it first becomes stable, i.e., the earliest time t ¼ t k;start for which Δ k;stable < ϵ stable . This can be interpreted as the time at which we confidently believe the mode to be present in the waveform. It can also be extracted from the frequency-fixed fits, which are labeled as filled circles in Fig. 4. Note in particular that t k;start ≠ t k;best ; the former is the time at which the mode first becomes stable (marked by filled circles in Fig. 4), while the latter is the starting time of the fitting window that gives the best stability (the starting time of the bold lines in Fig. 4). We will return to a more detailed discussion of the definition of the starting time of the mode in Sec. VI below.

We stress once again that the results for a particular Kerr mode depend on which other modes are included in the fit. Therefore, the amplitude/phases and starting times extracted from the fit refer to those obtained by fitting the whole waveform with a specific set of modes. For example, while the starting time of the 2,2,1 mode in Fig. 4 is t 0t peak ∼ 10M, we should expect it to start even later if we were to use a fitting model including (say) only the 2,2,0 and 2,2,1 modes, as often done when analyzing real GW data, because missing other modes will hinder the stability of the fit.

IV. HYPERFITS OF THE AMPLITUDES

The simplicity of BH spectroscopy is largely due to the no-hair theorem. As the remnant BH can be characterized fully by its final mass and dimensionless spin M rem and χ rem , its QNM frequencies do not depend on the initial conditions of the BBH merger, i.e., on the properties of the merger remnant progenitors. However, the degree of excitation of different QNMs varies depending on the initial conditions. The amplitude and phase of each QNM are functions of the parameters of the progenitors, including their spins χ 1 , χ 2 and the binary-mass ratio q ¼ m 1 =m 2 .

Here we choose the convention that m 1 ≥ m 2 (so that q ≥ 1).

Whether or not one cares about A and ϕ of a QNM depends on their goal. If one plans merely to test the no-hair theorem in a restrictive sense, then checking the consistency between the different QNM frequencies ω may suffice. However, given a set of initial BBH parameters, A and ϕ can be predicted in a well-defined manner within GR. More robust spectroscopy tests should also check the consistency of A and ϕ with the predictions of GR [75].

To do this we must build a model of A and ϕ as a function of the parameters of the binary progenitors. This is the goal of this section. We will restrict our attention to nonprecessing waveforms, i.e., simulations where both progenitor BHs have spins either parallel or antiparallel to the direction of the orbital angular momentum. These FIG. 4. Amplitude (top) and phase (bottom) extraction for SXS: BBH:0305, l; m ¼ 2, 2. A k and ϕ k are shown as bold lines within the τ stable time window when Δ k;stable is minimal, which is where the amplitudes and phases are extracted. We also label the time t k;start with a solid circle.

simulations are parametrized by q; χ 1 and χ 2 , where χ 1 ≡ χ 1;z ¼ AEjχ 1 j, where a positive (negative) sign means that χ 1 is parallel (antiparallel) to the orbital angular momentum of the binary, and similarly for χ 2 .

It will also be useful to use parameters that appear often in a post-Newtonian (PN) expansion of BBH dynamics,

Note that χ þ is also called χ eff in the literature.

For each simulation, we run the full mode-extraction algorithm presented in Sec. III. This yields data points for the amplitudes of different modes over the whole parameter space ðη; χ þ ; χ -Þ. In Fig. 5 we visualize the amplitude of the 2,2,0 mode in this three-dimensional parameter space. Interactive versions of the plots for this and other QNMs can be found on the Jaxqualin webpage [1].

Given the results of the fit at different points in parameter space, we can construct a phenomenological model of the amplitude as a function of the initial parameters of a BBH merger. We will only construct these models for modes that reside in their natural l; m multipole (e.g., for the 2,2,0 mode in the 2,2 multipole), but not for the mixing or recoil modes. This is because the mixing and recoil modes are subdominant, making their extraction noisier. In principle, one could map the amplitudes of a mode from their natural multipole to their mixing contribution in a different multipole by considering the spherical-spheroidal mixing coefficients (see Appendix F).

By inspection of Fig. 5, we observe that the amplitude of the 2,2,0 mode varies rather smoothly throughout the whole parameter domain. This is true in general for all the modes that we will consider in this section. In the literature, hyperfits of mode amplitudes are often performed by assuming a polynomial model in some combinations of the mass ratio (η or q) and of the individual progenitor spins (χ 1 , χ 2 , χ þ , or χ -) [47,51,54]. We will follow this approach, and we will adopt a strategy that tries to minimize overfitting due to the inclusion of too many terms in the polynomial.

When performing the hyperfits, we will consider the SXS simulations SXS:BBH:0209-0305 and SXS:BBH:1419-1509, for a total of 188 simulations. The first set of simulations was performed to validate and improve waveform models [118,119] (we also included SXS: BBH:0305, the most studied waveform in BH spectroscopy studies), while the second set was used to construct the NRHybSur3dq8 [120] and NRHybSur3dq8_CCE [121] NR surrogate models. These are high-quality simulations that cover the parameter space in question adequately. For completeness, we will also include the data points or fitting loss of the other nonprecessing simulations in the catalog, but those are not used to train nor validate the fitting model. All of the fitting and validation are performed with the SCIKIT-LEARN package [122].

B. Fixing the degree of the polynomial model

We now proceed to perform the hyperfits on the adjusted amplitudes Ãl;m;n . We choose to fit the amplitudes with a multivariate polynomial in three variables (η, χ þ and χ -), but we must choose the degree of the polynomial. To avoid underfitting or overfitting, we split the full set of data points randomly into a training set (consisting of 80% of all simulations) and a validation set (the remaining 20%). We fit the training set by a polynomial of three variables with a weighted least-squares loss for a varying degree N , including all cross terms. We use the inverse of the fluctuation ranges δA and δϕ as weights for the loss to penalize simulations that return less stable amplitude and phases. Then, we compute the loss of the validation set when using the fitted polynomial model, and identify the value of N that minimizes the validation loss. In Fig. 6, we can see that the polynomial degree that gives the lowest validation loss is N ¼ 4; this is the optimal degree to use without overfitting the data.

For modes other than the 2,2,0 mode shown in Fig. 6, we will often also find N ∼ 4. With N ¼ 4 and three variables, this polynomial is specified by 35 coefficients. In the next subsection, we will show that the number of coefficients in the polynomial can be reduced further. FIG. 5. The amplitudes of the 2,2,0 mode in the l; m ¼ 2, 2 harmonic waveform of all fitted nonprecessing SXS simulations. Top: the amplitudes visualized in the three-dimensional parameter space of η, χ þ and χ -. Each ball represents an SXS simulation at the corresponding location in parameter space, the color of the ball corresponds to the magnitude of A 2;2;0 extracted from that simulation. The size of the ball is proportional to the inverse of the fluctuation Δ k of the extracted amplitude; in other words, the larger balls represent amplitudes we measured more confidently, hence they will be weighted heavier when performing the hyperfits. Bottom row: same as the top plot, but projected onto 2D parameter spaces, each panel being the space spanned by a selected pair of the three parameters η, χ þ and χ -. Interactive versions of this plot and similar plots for other modes can be found on the Jaxqualin webpage [1].

C. Fixing the number of terms in the polynomial

A multivariate polynomial contains a significant number of cross terms even if the polynomial degree is moderate. In our model, while a degree of N might be necessary to capture the dependence of the amplitude, some of the terms in the polynomial might not be required for the model to work well. For this reason, we use an iterative scheme to remove unnecessary terms in the polynomial.

We begin with a polynomial of degree N containing all the i possible terms. For example, for N ¼ 2, there will be a set of i ¼ 10 terms, including 1; η;

We disregard one of the terms and fit the data with a polynomial that only includes the remaining i -1 terms. The training loss of the fit will be increased when we do so. We repeat this step for each of the i different terms, disregarding each term and keeping the other i -1 terms, keeping track of the induced increase in the training loss. In this way we can identify the term that leads to the smallest increase in the loss when it is neglected. This can be considered the term that contributes the least to the goodness of fit. We remove this least-useful term from the set. At this point we are left with i -1 terms, and we have completed the first iteration. We use the same procedure to iteratively remove the least-useful term, until we are only left with one term. This iterative procedure allows us to rank the importance of each term.

Given this ranking, we can try to determine the number of terms we have to include in order to adequately fit the data, i.e. the number j ≤ i for which the j most important terms give a good enough model. We do this, again, by inspecting the loss of the validation set. As shown in Fig. 7, for the 2,2,0 mode with N ¼ 4 and i ¼ 35, while the training loss decreases monotonically as we include more terms, the validation loss reaches a minimum and plateaus at around j ∼ 20. As a lower value of j could already give a validation loss that is almost as good as the minimum value, we pick the best value j best by searching for the lowest value of j such that the loss is less than 1.5 times the minimum validation loss. For most modes, we find that j best ≲ 15. In this way we can compress the dimensions of the hyperparameter space of the model by a factor of ≳2.

We can now repeat the same procedure for the amplitudes and phases of different modes. For some modes, we had to specify by hand a polynomial degree N that is larger than the "best" one found with the procedure in Sec. IV B, because the variation of the data contains maxima or minima within the fitting domain. Explicitly, we used N ¼ 5 for Ã2;1;0 ; Ã3;2;0 and Ã4;4;0 , and N ¼ 6 for Ã5;5;0 ; φ4;4;0 and φ5;5;0 . Overfitting is avoided by the steps explained in this sub-section to remove unnecessary terms from the polynomial.

For the modes considered, we found the following bestfit functions: FIG. 7. Fitting loss as a function of the number of terms we keep in the multivariate polynomial given the maximum degree N of the polynomial, for the case of Ã2;2;0 . In this case, we have already fixed N ¼ 4 (see Fig. 6). In practice, we start by including all the terms of the multivariate polynomial of degree N , and we remove terms iteratively, progressing from right to left in the figure. The best number of terms j best , determined by the validation set, is marked by a black square.

Ãr2;2;0 ¼ 0.001657

φr2;2;0 ¼ 13.14

From the hyperfits, the contribution of a particular QNM to the waveform strain at any time after t peak can be obtained as

η Ãl;m;n e -i½ω l;m;n ðt-t peak Þþð φl;m;n þmϕ 2;2;0 Þ=2 : ð13Þ

Here M is the total mass of the BBH progenitors, d L is the luminosity distance to the merger, and η ¼ η (η ¼ ηδ) for even (odd) values of m, respectively. The contribution of the QNM in question to the detector strain can be obtained by projecting Eq. ( 13) onto the antenna response function of the detector. If needed, one can map the progenitor parameters to the remnant mass and spin with, e.g., the surfinBH package [125]. All of the hyperfit functions are implemented in the Jaxqualin package, with examples shown on the code webpage. We omit modes that have a negligible amplitude in most of the simulations, as well as the overtones and quadratic modes whose amplitudes are too noisy. In fact, as we will discuss below, it is easier to model the quadratic modes as a function of the amplitude of the linear modes that source them. We have compared our amplitude hyperfits with Ref. [54] and confirmed that they follow the same trends, except at close-to-extremal values of χ þ or χ -, where there are fewer simulations available to train the hyperfit.

The hyperfit functions listed above should be used with caution. First of all, one should proceed with care when using the fits outside of the range of the parameter space covered by the data used for training. When doing the hyperfits for a given mode, we only used the simulations for which we identified the mode and extracted the amplitude and phase. Therefore, in a region of parameter space where we did not find the mode, the hyperfit is an extrapolation. Secondly, a high-dimensional multivariate polynomial might not be able to give a good fit for arbitrary data. In fact, for some of the modes, the fits should only be taken as an order-of-magnitude estimate, and they can have large relative errors in regions where the mode amplitudes are small. The error of the fits can be visualized on the Jaxqualin webpage [1], and it should be carefully considered before using the fits.

If the fitting error is too large for the user's purposes, we recommend falling back to an interpolation scheme. In the Jaxqualin package, we provide an interpolation function along with the hyperfit functions. The interpolator uses linear barycentric interpolation (as implemented in SciPy [126]) to interpolate between data points, avoiding the use of any fitting model. Using an interpolator bypasses the errors that inevitably arise from fitting, but it introduces fluctuations due to the noise of the data points. The interpolator does not return a value if the requested point is outside of the convex hull of the data points. This is somewhat restrictive, but it avoids returning spurious extrapolated results.

A. Quadratic modes

The hyperfitting models presented above are phenomenological. For quadratic modes, it has been shown, both theoretically and by examining NR simulations, that their amplitudes and phases depend quadratically on the linear modes that source them,

where all the amplitudes are assumed to be extracted from the "natural" multipolar components (l; m ¼ l 1 þ l 2 ; m 1 þ m 2 for the quadratic mode). The coefficients c l 1 ;m 1 ;n 1 ×l 2 ;m 2 ;n 2 and d l 1 ;m 1 ;n 1 ×l 2 ;m 2 ;n 2 can in principle be computed from the Teukolsky equation by plugging quadratic combinations of linear QNMs as source terms into the equation (although this has only been attempted in the nonrotating case and under several approximations; see e.g., [82,84,85,89,90]). In general, these coefficients are expected to depend on χ rem and (to some extent) on the initial conditions [83], and therefore on η; χ þ , and χ -.

With the fitting procedure in Sec. III, we have extracted both linear and quadratic modes. Therefore, we can compare the amplitudes and phases between these modes for each simulation to check whether the quadratic relationship holds. In Fig. 8 we show how the 2; 2; 0 × 2; 2; 0 mode amplitude and phase depend on the corresponding quantities for the "parent" 2,2,0 mode. We also illustrate the same relations for the 2; 2; 0 × 3; 3; 0 mode as a function of the properties of the 2,2,0 and 3,3,0 modes. These plots generalize and extend similar results presented in Ref. [78]. Interestingly, the dependence seems to be rather universal, in the sense that the constants c and d are only mildly (if at all) dependent on χ rem .

In Fig. 9 we plot c l 1 ;m 1 ;n 1 ×l 2 ;m 2 ;n 2 against χ rem . For both the 2; 2; 0 × 2; 2; 0 and 2; 2; 0 × 3; 3; 0 modes, the amplitude ratio follows an almost-constant linear trend with respect to χ rem for the majority of simulations. Only a small fraction of points deviate from this trend. At least some of these outliers could result from biases, either because the mode extraction algorithm misses some modes, or because it picks up some spurious ones. Note that A 3;3;0 → 0 as η → 1=4, which could explain the bias in the ratio of the 2; 2; 0 × 3; 3; 0 mode in the comparable-mass limit.

Assuming that c l 1 ;m 1 ;n 1 ×l 2 ;m 2 ;n 2 and d l 1 ;m 1 ;n 1 ×l 2 ;m 2 ;n 2 do not depend on χ rem , their best-fit values are c2;2;0×2;2;0 ¼ 0.144; ð15aÞ d2;2;0×2;2;0 ¼ 0.077; ð15bÞ c2;2;0×3;3;0 ¼ 0.448; ð15cÞ d2;2;0×3;3;0 ¼ 0.172;

where the overbar means that we assume these coefficients to be constants. If we allow the coefficient c 2;2;0×2;2;0 to be dependent on χ rem and assume that the dependence is linear, the best fit is c 2;2;0×2;2;0

For the other coefficients, we observe no significant dependence on χ rem .

B. Overtones

There has been some discussion on whether the 2,2,1 overtone has been observed in the first merger event GW150914 [104][105][106]127]. In fact, it is questionable whether the first overtone is useful at all for testing GR, even in the limit of infinite signal-to-noise ratio [96,109,110]. We contribute to this discussion by asking if there is some universal relation between the amplitude of the 2,2,1 overtone and the amplitude of the 2,2,0 fundamental mode.

In principle we could derive such a relation "directly" from the fitting Eqs. (12a) and (12e). By looking directly at the amplitudes, we find that their ratios is rather universal. As shown in the left panel of Fig. 10, the amplitude ratio A 2;2;1 =A 2;2;0 has a strong dependence on χ rem , but only a very mild dependence on the (symmetric) mass ratio. The ratio A 3;3;1 =A 3;3;0 follows a similar trend, although there is a larger spread in the data. An exponential in a power of χ rem seems to fit quite well most of the data. The best fit curves shown in the figure are given by A 2;2;1 =A 2;2;0 ¼ 6.519 ðχ 3.366 rem þ0.500Þ ;

Note that this relation might not hold in the extrapolation region, i.e., for χ rem ≲ 0.1 or χ rem ≳ 0.95.

In principle, it should be possible to compute the ratios in Eqs. (17a) and (17b) given the initial conditions for the perturbations of the remnant Kerr BH. In practice this is a difficult problem, as there is a large amount of arbitrariness in the choice of the initial time slice, gauges, etc. The ratios we extract from the numerical data are equivalent to the ratio between the so-called "excitation coefficients" C l;m;n of the two modes, which can be factorized as

Here the "excitation factors" B l;m;n ≡ B l;m;n ðχ rem Þ are functions only of the remnant BH spin, but independent of the initial conditions, while the integral I l;m;n depends on both the initial conditions and on the source exciting the perturbations [19,34,36,37]. Determining appropriate initial conditions for a comparable-mass BBH merger is highly nontrivial, and (to the best of our knowledge) the coefficients C l;m;n have never been computed "from first principles" in this context. However, recently there are claims that the ratio between overtone amplitudes is closely approximated by the ratio between the corresponding excitation factors B l;m;n . This would imply that the integral I l;m;n does not strongly depend on n for fixed valus of the multipolar indices l and m. In Fig. 10 we check the accuracy of this approximation by plotting the absolute value of the ratios between the excitation factors jB l;m;1 =B l;m;0 j (the numerical values were computed in Refs. [34,36] and are publicly available online [27]). The A l;m;n and B l;m;n ratios follow similar trends for low values of the dimensionless remnant spin χ rem , but there are significant deviations between the two at higher values of χ rem . This implies that the ratios A l;m;1 =A l;m;0 will, in general, depend on the initial data and on the source exciting the perturbations, although the dependence may be mild in certain favorable situations. This is an interesting problem that deserves further study.

A. An overoptimistic definition

Our goal so far has been to extract the amplitude and phase of a mode to the best accuracy from NR simulations. The natural procedure that we employed is to extract them within a time window where they are the most stable (i.e., where their amplitude varies the least).

However, when analyzing real data it is important to maximize the signal-to-noise ratio (SNR), so (ideally) we would like to extract the modes as early as possible without biasing the results. As the numerical noise in the NR simulations is negligible when compared to the instrumental noise of GW detectors, we can treat the mode extraction FIG. 10. The ratios A 2;2;1 =A 2;2;0 (left) and A 3;3;1 =A 3;3;0 (right) as functions of χ rem . The color of the circles corresponds to the value of η for that particular simulation. The black dashed lines are the best-fit curves given in Eqs. (17a) and (17b). The purple dotted lines are the excitation factor ratios jB 2;2;1 =B 2;2;0 j (left) and jB 3;3;1 =B 3;3;0 j (right). FIG. 9. The ratios A 2;2;0×2;2;0 =A 2 2;2;0 (left) and A 2;2;0×3;3;0 =A 2;2;0 A 3;3;0 (right) as a function of χ rem . The ratio shown on the left is equal to the coefficient c l 1 ;m 1 ;n 1 ×l 2 ;m 2 ;n 2 defined in Eq. (14a). The black dashed lines correspond to the lines shown in the left panels of Fig. 8, where c l 1 ;m 1 ;n 1 ×l 2 ;m 2 ;n 2 is assumed to be a constant (independent of χ rem ). The gray dotted lines are the best-fit lines if we instead assume a linear relationship between the ratios and χ rem .

results from NR as "ground truths" for finding the optimal extraction time in real detections.

For ringdown signals detected in real data, we can define the optimal starting time for extracting each QNM k to be the time t k;start at which the variation of the mode amplitude and phase (Δ NR A k and Δ NR ϕ k , as extracted from NR) is smaller than their statistical measurement errors (Δ stat A k and Δ stat ϕ k , as extracted from real data, e.g., at the 2σ level). Assuming that the statistical error is sensible (Δ stat A k < A k and Δ stat ϕ k < 2π at the very least), this can be seen as the earliest time at which the mode extraction from real data will not be significantly biased compared to the statistical error. The quantities Δ stat A k and Δ stat ϕ k are well-defined, and they can be read off from the posterior of the analysis. The quantities Δ NR A k and Δ NR ϕ k depend on the length of the time window τ NR over which we measure the spread of A and ϕ (compare the τ stable window used in Sec. III).

Clearly, this definition of t k;start depends on the level of instrumental noise in a GW detector. However, we can still develop some intuition on t k;start by imposing that Δ NR A k and Δ NR ϕ k should vary less than a sensible range. For example, we can reuse our criteria in Sec. III C, i.e., the quantities ϵ stable ; β A ; β ϕ ; A tol ; τ stable , and p stable %, using the "normal" stability condition listed in Table I, and define t k;start to be the earliest time t 0 at which Δ k;stable ðt 0 ; τ stable Þ < ϵ stable . For example, t k;start is labeled as a solid circles for each mode in Fig. 4. While there is some arbitrariness in the choice of the parameters, a value of ϵ stable of order Oð0.1Þ is in line with the sensitivity of current GW detectors, while τ stable ∼ 10M is comparable to half a period of the dominant 2,2,0 mode (the natural time scale of the ringdown).

Such a definition of t k;start is not only somewhat arbitrary, but it also has a crucial shortcoming. Because we define a different t k;start for each mode k, there is no guarantee that at the t k;start appropriate for a particular mode, the other modes will also be stable. This can be seen from Fig. 4. If we focus, for example, on t 2;2;0; start , then the amplitudes and phases of all other modes are unstable. As discussed in Ref. [96], by including multiple QNMs in the fitting model we can overfit the merger part of the signal, and thus drive the amplitude of the dominant QNMs (in this case, the 2,2,0 mode) to their "correct" late-time value at earlier and earlier times. In other words, while from Fig. 4 we may be tempted to conclude that t 2;2;0; startt peak ∼ 0M, this could be an artifact of overfitting the merger with all the other modes, as can be deduced from their unstable amplitudes. This implies that the t k;start defined in this way could be overoptimistic (i.e., it tends to be "too early").

Such overoptimism would only be exacerbated by using a model containing fewer modes than those used for extracting t k;start in this work (e.g., less than the 6 modes listed in Fig. 4). The goal of the algorithm specified in Sec. III is to find all identifiable modes within NR, and to include them all in the fitting model when extracting the amplitudes and phases of the modes. At current detector sensitivity, we could only hope to measure a few (∼2) QNMs in real data, so it would not be sensible to always include so many modes in the fitting model. When that is the case, we should use a larger t k;start than those shown in this work, because the negligence of subdominant modes could affect the stability of the dominant ones.

C. A recommended strategy

In practice, the "true" optimal starting time should be somewhere in between these overoptimistic and overconservative definitions. A recommended strategy is to first decide which modes should be included in the model, fit only with those modes, and then use the overconservative strategy to decide the optimal starting time for the full model (including all QNMs, not just one particular QNM).

A sketch of such a procedure is as follows:

(1) Prepare the waveform to be fitted.

(2) Decide the set of modes K ¼ fkg to use and the tolerance for stability, e.g., the values of ϵ stable and τ stable . The tolerance parameters could be determined from the statistical error obtained when performing data analysis. (3) Perform a fit with all modes k ∈ K, and find t k;start for each of them. (4) The optimal starting time t start is the largest t k;start among all modes k ∈ K. For example, if one is using the l; m ¼ 2, 2 multipole of a BBH simulation and chooses K ¼ f2; 2; 0g or f2; 2; 0j2; 2; 1g, one would find t start ≳ 10M in both cases.

If we want to quantify the optimal time for performing data analysis on real signals in detector noise, it is important to make sure that the waveform prepared in Step 1 should be as close as possible to the one expected in the detector. For example, the actual detected waveform should contain all of the multipolar components of the waveform projected onto the detectors via the antenna pattern functions, which depend on the sky position of the source and on the orientation of the detectors. This is significantly more complicated than dealing with the individual multipolar components considered in this work.

We stress once again that the t k;start values shown in our figures uses the overoptimistic definition of Sec. VI A and, they should not be applied to analyze real GW data. This is because the merger is generally overfitted by the subdominant modes, because the number of modes used in this work will usually be different from the number of modes included in real data analysis, and because we only considered individual multipoles rather than the (sensibly more complicated) sum of multipoles that is relevant for real detections. The steps laid out in this subsection suggest a better strategy for determining t start in real data.

Another possibility for determining t k;start is to take the amplitude and phases extracted for different modes in Sec. III (or the hyperfits in Sec. IV) as the ground truth, and to perform an array of data analysis runs for different values of the starting time t 0 to find the earliest time for which all the QNM frequencies, amplitude and phases are consistent with the "true" values. This could be a viable strategy, but it requires performing multiple data analysis runs (as in Refs. [105,106,131]) or marginalizing over the ringdown starting time (as in Ref. [104]).

VII. CONCLUSIONS

In this work, we have presented a strategy for extracting QNMs from BBH merger simulations. The core procedure is sketched in Fig. 1 and detailed in Sec. III, with tolerance criteria summarized in Table I. The strategy is implemented in the publicly available Jaxqualin package [1], and the tolerance parameters can be tuned for the user's purposes.

By fitting an exhaustive set of nonprecessingquasicircular BBH waveforms in the SXS catalog, we have constructed hyperfitting models for the amplitudes and phases of different QNMs. We have verified the expected scaling between quadratic modes and the linear modes that source them, confirming and extending previous results [51,77,78]. We have also found quasiuniversal relations between different modes, including the amplitude ratio between the first overtone and the fundamental mode shown in Fig. 10, and the amplitude ratio between the retrograde and prograde modes shown in Fig. 11. To the best of our knowledge, these results are new. These quasiuniversal relations warrant further investigation, both in waveform modeling and in BH perturbation theory. The apparent discrepancy between the remnant spin dependence of nonlinear modes extracted from BBH simulations (see Fig. 9) and the dependence found for Gaussian scattering at second order in BH perturbation theory [83] deserves further investigation.

Finally, we have presented some practical proposals to define a convenient ringdown starting time for data analysis purposes. These proposals should be tested and refined by analyzing ringdown signals in real data.

It is trivial to apply the Jaxqualin code to other NR waveforms or simulation catalogs, such as those using CCE, or precessing BBH simulations. We expect that higher-quality simulations including center-of-mass corrections, improved Cauchy-characteristic evolution, the inclusion of gravitational memory using BMS balance laws, and improved wave extraction methods would be very useful to remove spurious QNMs, to improve the accuracy of our hyperfitting models, and to shed further light on the relation between the amplitudes and phases of different modes. providing resources that have contributed to the research results reported within this paper [132]. We made use of the following PYTHON packages to produce the results of this paper: adjustText [133], Jax [117], JaxFit [116], Matplotlib [134], NumPy [135], PANDAS [136,137], QNM [115], SCIKIT-LEARN [122], SciPy [126], SCRI [138] and SXS [139].

APPENDIX A: DEFINITION OF PROGRADE AND RETROGRADE MODES

Given a BH remnant with spin χ rem , there are two solutions for the QNM frequencies with given values of l and m; one is corotating with the BH remnant, and the other is counterrotating. In the literature, these are usually called "prograde" and "retrograde" modes, respectively.

Here we employ a different definition that is more useful for our work. We define the modes corotating (counterrotating) with the BBH orbit to be prograde (retrograde), respectively. This definition is clarified in Fig. 12. We implicity assume that the waveform is extracted in a frame aligned with the BBH angular momentum, so prograde modes with positive frequency are corotating with the BBH orbit. We also assume that χ rem < 0 corresponds to the remnant BH spin being antialigned with the orbital angular momentum.

We use this unconventional definition because in BBH simulations, the orbital direction is easy to control, while the remnant spin direction depends nontrivially on the initial conditions. By using the present definition, when we vary the initial conditions (e.g., q; χ 1 ; χ 2 ) smoothly, the amplitude and phase of the mode will also vary smoothly, even when the spin direction of the remnant flips. A similar behavior can be observed for the QNM frequency of a mode: for example, both the prograde mode (2,2,0, right panel) and the retrograde mode (r2, 2, 0, left panel) shown in Fig. 12 evolve smoothly on the same line from negative to positive χ rem . If we were to use the conventional definition, the QNM frequencies, amplitude and phase would undergo nonsmooth transitions when the remnant spin flips.

To summarize, our definition follows the usual convention for remnants with spins aligned with the BBH orbital direction (which is the case for the large majority of the nonprecessing simulations in the SXS catalog). However, when the remnant spin is flipped, our definition is the opposite of the usual convention.

FIG. 12. Our definition of prograde and retrograde modes. For illustration, we show the two QNM frequency solutions for l; m ¼ 2, 2. We show the negative-ω r region of the complex plane in the left panel, and the positive-ω r region in the right panel. The blue (orange) branches correspond to the case when the remnant BH spin is aligned (antialigned) with the BBH orbit. Each branch goes from jχ rem j ¼ 0 up to 0.95, with the points at jχ rem j ¼ 0.9 labeled by circles. In our definition, we call both branches in the right panel the prograde mode (2,2,0), while those in the left panel are the retrograde mode (r2, 2, 0). This definition generalizes to other l; m multipoles, but for m < 0 the plot should be mirrored (i.e., ω r → -ω r ). With this definition, we have implicitly assumed that the waveform is extracted in a frame such that the z-axis is coaligned with orbital angular momentum of the BBH.

The second class of toy waveforms corresponds to the second hypothesis, for which the waveform contains clean overtones as early as the peak with little contamination. Explicitly, these toy waveforms contain QNMs with strictly constant amplitudes and phases,

For simplicity, we include only the modes K ¼ f2; 2; 0j 2; 2; 1j2; 2; 2j2; 2; 3j2; 2; 4j3; 2; 0g in the waveform. The amplitude and phase of the 2,2,0 mode are fixed at 1 and 0, respectively. The amplitude of the 2,2,1 mode is drawn uniformly from (3,5), 2,2,2 from (8,12), 2,2,3 from (12,24), 2,2,4 from (20,30), and 3,2,0 from ð1 × 10 -3 ; 1 × 10 -2 Þ. Similar to the previous toy waveform, we inject Gaussian noise nðtÞ with standard deviation σ n ¼ 5 × 10 -7 .

APPENDIX D: CONSISTENCY TEST OF THE FITTING ALGORITHM

In the main text we fit the postmerger waveform with up to N f ¼ 10 free QNMs using model (3), corresponding to 40 free parameters in total. In general it will be difficult for the algorithm to locate the global minimum of the leastsquares cost function in such a high-dimensional space. Nonetheless, we will see that we do not have to insist on finding the global minimum best fit, because a local minimum could already serve our purposes. As long as we have fitted all of the physical (QNM) content of the waveform by diving below the noise floor, we do not have to find a better minimum, because by doing so we are only fitting the noise better. Moreover, even if we only found a local minimum, we will show that a free QNM will approach a Kerr mode as long as it stays long enough above the noise floor.

We first test the behavior of the fitting algorithm with the two waveforms specified in Appendix C. Recall that, for the toy waveforms of the first class, the waveform approaches a combination of constant amplitude/phase QNMs at late times (modulo Gaussian noise), while those of the second class are always pure QNMs. To test the ability of the algorithm to find an adequate least-squares minimum, we perform the same fit with different initial guesses. For each free QNM k in the fitting model (3), we draw the initial guess of A k from a log-uniform distribution in ð0.1; 10Þh peak , where h peak is the amplitude of the waveform at t peak . We also draw the initial guesses for ϕ k from ð0; 2πÞ, ω r from ð-2; 2Þ, and ω i from ð-1; 0Þ. Note that we do not put bounds on ϕ in the trust region algorithm for nonlinear least squares fitting, but we take ϕ k mod 2π as the extracted phase value. Nonetheless, an initial guess drawn between in ð0; 2πÞ covers adequately the relevant space of values that ϕ k can take.

In Fig. 13 we plot the fit mismatch versus t 0 for 10 runs with different initial guesses on the same waveform. We do this separately for a realization of the toy waveform of the first class (top panel) and of the second class (bottom panel). Recall that for a particular time step t 0 , we use the best-fit value from the previous time step t 0 -Δt 0 as the initial guess, so the initial guess that we specify is only used for the first time step (t 0t peak ¼ 0M).

For the toy waveform of the first class, no matter which initial guess we choose for t 0t peak ¼ 0M, the fits always converge to the same mismatch at times t 0t peak ≳ 40, when the waveform approaches a combination of clean QNMs. For the toy waveform of the second class the fit mismatch is always the same, regardless of the fit starting time t 0t peak or the initial guess realization. This test FIG. 13. Fit mismatch versus fit starting time for different initial guesses on a realization of the toy waveform of the first kind (top panel) and of the second kind (bottom panel). Each orange curve corresponds to a series of fits where the initial guess of the fitting algorithm at t 0t peak ¼ 0M is drawn randomly (a different realization for a different curve), and subsequent time steps make use iteratively of results from the previous time step as an initial guess.

suggests that if a waveform were a combination of clean QNMs (as long as we fit it with a number N f of free QNMs greater than its actual QNM content), the goodness of fit would be minimally affected by the initial guess.

The disagreement between the fits with different initial guesses when the waveform is not a clean combination of QNMs might look concerning. However, we know that at early times the waveform contains non-QNM contributions anyway, and this is the reason for the disagreement in the first place. At early times, if we were to insist in looking for the global minimum, we could have used a large number of initial guesses that exhaustively cover the relevant parameter space and select the instance that gives the lowest mismatch. Such an attempt would not help us extract more of the QNM content of the waveform as long as the original (local) minimum already has a mismatch lower than that induced by non-QNM components, such as numerical noise. Nonetheless, as mentioned in Appendix B, we use 10 initial guesses for each waveform and select the best instance. In our experience this strikes a balance between the computational resources required and the ability of the fit to find an adequate minimum. In fact, as seen in the top panel of Fig. 13, the fit mismatch is always below 10 -7 -well below the noise level, regardless of the initial guess we use.

Next, we perform the same test on an NR waveform, namely the l; m ¼ 2, 2 multipole of SXS:BBH:0305. The results are shown in Fig. 14. For N f ∈ f5; 6; …; 10g, a majority of the 10 initial guesses give similar mismatch curves. More importantly, they either converge to the same mismatch (e.g., for N f ¼ 5 and 6, close to t 0t peak ∼ 0M), or they dive beneath the resolution-error induced mismatch for N f ≥ 7.

Although the mismatch converges for different initial guesses when the data matches the model, there is no guarantee that the recovered ω k;r and ω k;i are also minimally affected by the initial guess. As shown in Fig. 15, for SXS:BBH:0305 with l; m ¼ 2, 2, the values of ω k;r and ω k;i obtained with different initial guesses could differ initially at t 0t peak ¼ 0M, but most of them converge to the values corresponding to the initial guess that gives the lowest mismatch at late times. An exception are the points around the 2,2,1 mode, which converge at around t 0t peak ¼ 30M only, the time when the 2,2,1 mode can be cleanly identified above the noise without a significantly varying amplitude or phase. The two toy waveforms FIG. 14. Fit mismatch for different initial guess realizations for SXS:BBH:0305, l; m ¼ 2, 2. We use model ( 3) for these fits. Different panels correspond to fits using a different number of free QNMs N f . Each blue curve corresponds to a series of fits where the initial guess of the fit algorithm at t 0t peak ¼ 0M is drawn randomly, and subsequent time steps make use of the previous time step results as initial guesses. The dashed line indicates the mismatch induced by the resolution error of the simulation, which is estimated by the (phase-minimized) difference between the waveforms extracted from the best and second-best resolution simulations. At any given t 0 in any panel, either most of the blue curves are under the resolution-induced mismatch, or all of them converge to the same values. examined earlier, as well as other l; m multipoles of SXS: BBH:0305, show the same qualitative behavior as Fig. 15: the consistency between free QNMs of different initial guesses is better around a Kerr mode when it is stably above the noise. However, if certain Kerr modes exist only briefly above the noise floor (e.g., a mode with very low amplitude or one with an overtone number n ≥ 2), they might not be picked up by the free QNMs at all, especially if N f is smaller than the number of Kerr modes actually present in FIG. 15. Fitted QNM frequencies at different t 0 with different initial guesses, for the l; m ¼ 2, 2 multipole of SXS:BBH:0305. Different panels correspond to different t 0 , and different colored circles in the same panel correspond to the recovered f ω k ¼ ω k;r þ iω k;i for different initial guesses. We make the circles transparent to better visualize the locations where different points overlap, which appears with a deep gray color. For each initial guess, we fit the waveform with N f ¼ 10 free QNMs. Notice how at late times (t 0t peak ≳ 40), different initial guesses give mostly overlapping points. The same qualitative behavior can be seen in other waveforms in general; as long as a Kerr mode is stable above the noise long enough at a given time, free QNMs will converge around the mode, regardless of the initial guess. the waveform. Some free QNMs might try to pick up these modes, but the precise value of the free QNM could depend on the initial guess, i.e., at a particular t 0 different initialguess runs would give different frequencies around the brief Kerr mode, forming a small spread. In our procedure, we introduce a margin for error by requiring only that the free QNM approaches a Kerr mode within an ellipse in the complex plane. This allows for a small bias due to an initial guess that does not return the global minimum. If the bias is larger than this margin of error, the Kerr mode that the free QNM is trying to pick up has too low an amplitude (or decays too quickly) anyways, so it is not too harmful to lose the mode. On the other hand, if N f is larger than the number of Kerr modes that are consistently above the noise floor, some free QNMs will be (over)fitting the non-QNM components of the waveform, and their frequency evolution could depend significantly on the initial guess. As explained in the main text, this is not a big concern, because even if the overfitting causes us to pick up spurious modes, they will be removed in a later step of the mode-extraction procedure (i.e., they will fail the mode stability tests).

With the above caveats in mind, the results in this appendix indicate that while there is no guarantee that we can find the global least-squares minimum, as long as a Kerr mode stably exists above the noise floor for a long enough time, the mode will most likely be picked up, no matter the initial guess. Using ∼10 different realizations of the initial guess and picking the best instance (the one that gives the lowest mismatch) is adequate for extracting all of the QNM components of a waveform.

When n ≥ 2 overtones are present

Next, we test whether our procedure is capable of extracting overtones with mode number n ≥ 2. Similar to the previous subsection, we construct toy waveforms by including different modes with randomly drawn amplitudes and phases. However, in contrast to the previous case, we keep the mode amplitudes and phases constant close to the merger peak, meaning that the QNMs are "clean" throughout the whole ringdown waveform. We do this because it has been hypothesized in the literature that overtones exist cleanly close to the merger peak, and we want to test whether our procedure will work well if that is the case. Recent work has found that, if such a hypothesis were true, the overtone amplitudes will also respect a hierarchy, where the amplitude increase with n for n ≲ 5 [48]. Therefore, we impose such a hierarchy when drawing the amplitudes of the overtones, as specified at the end of Appendix C.

The top row of Fig. 17 shows the mode extraction results for the case when we include the n ¼ 2, 3, 4 overtones, in addition to the 10 modes included in the previous section. We find that our procedure always prefers to pick up modes that are less quickly decaying, instead of the n ≥ 2 overtones [110]. To test whether it is inherently impossible to find these overtones with our procedure, we do another test, this time including only the 3,2,0 mode in addition to the 2; 2; n fundamental mode and overtones, with n ≤ 4. In this case, we can find the n ¼ 2 overtone for nine out of ten realizations, but we still cannot identify the n ≥ 3 overtones. Note that the detectability of the modes depends on their amplitude, so if the modes do not follow the hierarchy that we prescribed, the results could be different.

The results in this subsection indicate that our procedure is not suitable for extracting overtones with n ≥ 3, while it can in principle extract the n ¼ 2 mode if few other modes are present. This is expected, because higher overtones decay faster. Other methods might be necessary for extracting these overtones. However it is hard to imagine that these higher overtones may be useful in GW data analysis when they cannot even be extracted with a leastsquares fit in a minimal amount of noise, as we show here.

3. When χ rem ≲ 0.1

In the limit χ rem → 0, modes labeled by the same l and n but by a different azimuthal index m will converge to the same (Schwarzschild) frequency, making it difficult to distinguish between them when fitting the waveform. This is a problem because of the presence of recoil modes, which can have the same l but a different m.

In Fig. 18, we show the mode extraction results for χ rem ¼ 0.1 and 0.01. When χ rem ¼ 0.1, at least one of the recoil modes is missed in most of the cases, while for χ rem ¼ 0.01 we always miss either some recoil modes or spherical-spheroidal mixing modes. Comparing with Fig. 16, we can also see that the amplitudes and phases of the other modes are not recovered as accurately when the remnant spin is low. For this reason, we do not include the simulations that give χ rem ≤ 0.1 when constructing the hyperfitting model in Sec. IV of the main text.

APPENDIX F: VERIFICATION OF SPHERICAL-SPHEROIDAL MIXING MODE AMPLITUDES

Spherical-spheroidal mode mixing implies that within a waveform of multipolar component l; m we can find QNMs with the same m but different l [44,76,96,97,99,100]. The degree of mixing can be characterized by the sphericalspheroidal mixing coefficients μ m;l;l 0 ;n , defined by

where S l 0 ;m 0 ;n is the spheroidal harmonic eigenfunction for the l; m multipole and nth overtone with spin weight s ¼ -2. Thus, the complex amplitude A l;m;njl;m of the l; m; n QNM in its natural multipole l; m is related to the amplitude A l;m;njl 0 ;m in the l 0 ; m multipole by

In our mode extraction procedure, we do not assume any relation between the amplitudes of a QNM in its "natural" multipole and its amplitude in other multipoles where it may appear because of spherical-spheroidal mixing. Therefore, we can use the known values for μ m;l;l 0 ;n to verify whether the amplitudes of the mixing modes we extract are indeed consistent with the theory. FIG. 16. Mode extraction results for a constructed toy ringdown waveform containing the modes K ¼ f2; 2; 0j2; 2; 1j3; 2; 0j4; 2; 0j 3; 3; 0j2; 1; 0j4; 4; 0jr2; 2; 0jr3; 2; 0jconstantg, for which the remnant BH spin is set to be χ rem ¼ 0.7. We use ten realizations of the waveform, each constructed by drawing the amplitudes and phases of different modes randomly, while loosely respecting a hierarchy as prescribed in Appendix. C. Each "column" of points represents the relative extraction error of the amplitude (left panel) and phase (right panel), defined as δA=A and δϕ=2π, for different realizations of the waveform. Gray circles are the points corresponding to realizations for which all of the injected modes are identified by our procedure without any spurious modes, while colored crosses correspond to realizations where some modes are missed by the algorithm. Crosses of the same color correspond to the same realization, and the mode (or modes) that are missed for each realization are labeled in the legend.

Given a QNM extracted in its natural multipole with amplitude A l;m;njl;m , we can compute by Eq. (F1) its theoretical amplitude A t l;m;njl 0 ;m when mixed into the l; m 0 mode, and compare with the same quantity extracted through our procedure, A e l;m;njl 0 ;m , by defining

This quantifies the relative difference between the theoretical mixing amplitudes and those extracted from the fit. We define the relative error δA=A by choosing the denominator to be A ¼ jA t l;m;njl;m j. We perform the test for the 2,2,0 mode mixing into the l; m ¼ 3, 2 multipole, and for the 3,2,0 mode mixing into the l; m ¼ 2, 2 multipole. We use these mixing modes as examples because they are often dominant. The results are shown in Fig. 19 for the main set of 188 SXS simulations. The difference between the theoretical and extracted mixing amplitudes shows broad consistency between the two. The exception are modes corresponding to remnants with small spin χ rem , because then the recoil modes can interfere with the spherical-spheroidal mixing modes, as explained in Appendix E.

APPENDIX G: COMPARING THE RESULTS WITH AND WITHOUT CAUCHY CHARACTERISTIC EXTRACTION

As discussed in the main text, the postmerger signal of a given multipole is contaminated by modes from other FIG. 17. Same as Fig. 16, but for toy waveforms including 4 overtones f2; 2; 1j2; 2; 2j2; 2; 3j2; 2; 4g and keeping all the other modes used in Fig. 16 (top), and including only the 2,2,0 and 3,2,0 mode in addition to the overtones (bottom). We label realizations where modes are missed as crosses and those with both missed modes and spurious modes as diamonds. Points shown as gray circles correspond to realizations with no missed modes or spurious modes. multipoles due to a combination of spherical-spheroidal mixing and gravitational recoil. By applying CCE and by appropriately fixing the BMS frame, the recoil modes can in principle be removed [63,76]. This can produce a cleaner ringdown signal and allow for better extraction of the modes, especially for slower spinning remnants. The publicly available CCE waveform catalog currently contains 13 waveforms, and they do not cover the progenitor parameter space as extensively as the original SXS waveform catalog. Nonetheless, we can compare the results between the SXS and CCE waveforms for selected cases to understand the improvements that may be possible by using the CCE waveforms. For the CCE waveforms, we use the scri code [138,[140][141][142] to map the waveform into the correct BMS frame for QNM extraction. In principle, this removes the recoil modes and the constant shift in the waveform [63,76].

While the waveforms we compare come from effectively the same progenitor settings, the sampling rate of the CCE waveforms is around an order of magnitude lower than the sampling rate of the SXS waveforms. Each data point for CCE is Δt ≈ 1M apart, while that for SXS is Δt ≈ 0.1M. This limits the scope of the comparisons done in this appendix.

As an example, we compare the CCE:0001 and SXS: BBH:0389 waveforms, which are both equal-mass simulations with nonspinning progenitors (the CCE waveforms are usually labeled as SXS:BBH_ExtCCE:XXXX; for simplicity, we will call them CCE:XXXX in this appendix). We first test the frequency-agnostic fitting procedure by fitting both waveforms with N f ¼ f5; 6; …; 10g, and selecting the case in which we find the highest number of potential modes (cf. Sec. III B 2). In Fig. 20 we compare the results between the best instance among N f ∈ f5; 6; …; 10g for both CCE:0001 and SXS:BBH:0389 for the ðl; mÞ ∈ fð2; 2Þ; ð3; 2Þ; ð4; 4Þg multipoles. In the SXS results we observe traces of the recoil modes, including for example the 4,4,0 and 2; -2, 0 modes in the 2,2 multipole. As expected, in the CCE waveforms, the recoil modes are suppressed and we do not observe them. An exception is the 2,2,0 mode in the 4,4 multipole, which is still identified as a potential mode in the CCE waveform. However, this mode did not pass the amplitude flatness test (unlike the corresponding mode found in its SXS counterpart).

As the CCE waveforms are cleaner, in some cases we can identify more modes, especially overtones and quadratic modes. For example, in the l; m ¼ 3, 2 multipole of the CCE waveform we observe the 2,2,1 and 3,2,1 overtones and even the retrograde r2; 0; 0 × 2; 2; 0 quadratic mode, with all of them passing the flatness test. By contrast, in the SXS waveform only the 2,2,1 and 3,2,1 overtones are listed as potential modes, and both of them fail the flatness test.

Next, we compare the results of the fixed-frequency fitting analysis. For both waveforms, following Sec. III C, we perform the stability test on the list of potential modes to arrive at a list of modes that we can consider "robust." Then we can do a final fixed-frequency fit to extract the amplitude of different modes. In Fig. 21 we compare the amplitudes extracted from the CCE and SXS waveforms for the same three multipoles examined above. While the list of modes found is different, the amplitudes of the modes that were found in both waveforms are similar. The amplitudes of the most dominant modes-in particular, modes with multipolar indices l and m that match the NR waveform multipole under consideration-are nearly equal across the two waveforms, unless some important overtones are missed. For example, as the 2,2,1 and 3,2,1 overtones failed the stability test and were discarded for the l; m ¼ 3, 2 multipole of the SXS waveform, the extracted amplitude of the 3,2,0 mode is different from the amplitude extracted from the CCE waveform at the percent level. Note that CCE alters the spurious constant in the waveforms, so a difference in the amplitude of the constant mode is expected.

The tests performed in this appendix show that using CCE waveforms is particularly important for robustly extracting QNMs. However, the results obtained by fitting non-CCE waveforms are still useful, as long as we identify all the important modes present in the waveform. A bias could be introduced if some modes are missed or misincluded, but the bias is typically small when compared to, e.g., the fluctuation of the mode amplitude with respect to t 0 , or to the error of the hyperfits presented in Sec. IV.

APPENDIX H: ADJUSTING THE STRENGTH OF THE STABILITY CONDITION

In the main text we use certain threshold parameters, such as p ag and ϵ stable , for determining whether a mode is robustly present in the ringdown. While the mode identification and stability conditions are qualitatively sound, these thresholds are largely arbitrary, and they were determined "experimentally" by trial and error. In our experience, the threshold parameter that most sensitively affects the mode extraction procedure is ϵ stable , the tolerance in the fluctuation of the amplitude and phase (Δ k;stable ). Changing the value of ϵ stable affects the number of modes that we find. In this appendix we show how this, in turn, affects the amplitudes of the dominant modes that we would like to extract. We will call our standard parameter choice in the main text (ϵ stable ¼ 0.2 and A tol ¼ 10 -3 ) the "normal" condition. Here we compare the "normal" results against "strong" and "weak" thresholds, such that ðϵ stable ; A tol Þ ¼ ð0.1; 10 -7 Þ and ð0.4; 10 -3 Þ, respectively. We call the former "strong" because potential modes found by the frequency-agnostic fit are more easily discarded by the stability test. The opposite is true for the "weak" condition. While ϵ stable affects mode extraction the most, we also change A tol to 10 -7 for the strong condition to verify that this parameter does not affect the extraction of the amplitudes and phases the dominant modes. As expected, using the strong (weak) condition we find fewer (more) modes. In Fig. 22 we show how this affects the extraction of amplitudes of different modes within their natural multipole. We define the relative error due to using a different condition to be

and we show the relative error δA cond =A by choosing A ¼ A weaker . For example, in the top row of Fig. 22 we use "strong" as the stronger condition and "normal" as the weaker condition. The results include all the simulations used in the training/validation set of Sec. IV. The relative errors δA and δϕ are mostly at the percent level, but they can be as large as Oð10%Þ for relatively subdominant modes. This is not too concerning: in real GW data, we would not expect the amplitudes of those modes to be measured better than that for current detectors.

APPENDIX I: OUTLIERS IN THE EXTRACTED AMPLITUDES

The mode extraction algorithm outlined in the main text uses tolerance criteria to check whether the modes are present in the waveform. Choosing the best criteria is inherently difficult; if the criteria are too stringent we will lose some modes and if they are too lax, we will identify spurious modes that should not be present. In principle there is an infinite number of modes in the waveform, but most of them will be negligible. Nonetheless, our goal is to extract all of the modes that we could identify. When we say that a mode is "present" in the waveform we mean that it is "not negligible", i.e., that the mode has a stable amplitude and phase, and missing it would cause a significant bias in the extraction of the other modes.

In fact, as we are using the same criteria for all modes across all simulations, there is no guarantee that there exists a set of criteria that would identify all modes present in all simulations without finding spurious modes. Therefore, as FIG. 22. Violin plots of the relative difference in the extracted amplitude (left column) and phase (right column) for different modes. In the top row we compare the "strong" and "normal" conditions, while in the bottom row we compare the "normal" and "weak" conditions. The violins denote the distribution of the errors for all of the simulations considered. The numbers above the plotting panel show the percentage of instances for which the amplitude extracted for that particular mode has a nonzero change across the different conditions. The change in the condition only affects the number of other modes included when we are extracting the mode in question. Thus, if the two conditions we are comparing find the same list of modes, the amplitude of the extracted mode will be exactly the same, because the fit includes the same modes across the two different conditions. the misidentification of other modes (either missing modes or finding spurious modes) could affect the extraction of a given mode, it is not surprising to see the extracted mode amplitude and phases being biased for a small fraction of simulations. For example, the data points lying above the trend in the top left panel of Fig. 8 (or in the left panel of Fig. 9) could be biased for this reason. In fact, by inspecting the mode extraction results for these simulations, we find that they often miss certain recoil modes that are identified in other simulations.

In some figures, we have removed one or two outliers that are significantly biased (e.g., one or more orders of magnitude away from the observed trends), or those with large fluctuations in the extracted amplitude and phase. These simulations are listed in Table II.

TABLE II. Outlier data points removed in some figures in the main text. These simulations were removed either because they are significantly far from the trends observed in other simulations, or because their amplitudes or phases fluctuate significantly (i.e., their error bars are one or more orders of magnitude larger than the norm).