# High-overtone fits to numerical relativity ringdowns: Beyond the dismissed
$n=8$ special tone.

Xisco Jiménez Forteza1,2, Pierre Mourier1,2 1 Max Planck Institute for
Gravitational Physics (Albert Einstein Institute), Callinstraße 38, 30167
Hannover, Germany  2 Leibniz Universität Hannover, 30167 Hannover, Germany

###### Abstract

In general relativity, the remnant object originating from an uncharged black
hole merger is a Kerr black hole. This final state is reached through the
emission of a late train of radiation known as the black hole ringdown. In
linear perturbation theory around the final state, the ringdown morphology is
described by a countably infinite set of damped sinusoids — the quasinormal
modes — whose complex frequencies are solely determined by the final black
hole’s mass and spin. Recent results advocate that ringdown waveforms from
numerical relativity can be fully described from the peak of the strain
onwards if quasinormal mode models with $N_{\mathrm{max}}=7$ overtones (beyond
the fundamental mode) are used. In this work we extend this analysis to models
with $N_{\mathrm{max}}\geq 7$ up to $N_{\mathrm{max}}=16$ overtones by
exploring the parameter bias on the final mass and spin obtained by fitting
the nonprecessing binary black hole simulations from the SXS catalog. To this
aim, we have computed the spin weight $-2$ Kerr quasinormal mode frequencies
and angular separation constants for the $(l=m=2,n=8,9)$ co- and counter-
rotating overtones, which all approach a Schwarzschild algebraically special
mode at low spins. We provide tables of the values obtained for these modes,
which are in agreement with previous results. From the systematic
variable-$N_{\mathrm{max}}$ analysis, we find that $N_{\mathrm{max}}\sim 6$
overtones are on average sufficient to model the ringdown from the peak of the
strain, although about $21\%$ of the cases studied require at least
$N_{\mathrm{max}}\sim 12$ overtones to reach a comparable accuracy on the
final state parameters. Considering the waveforms from an earlier or later
point in time, we find that a very similar maximum accuracy can be reached in
each case, occurring at a different number of overtones $N_{\mathrm{max}}$. We
also provide new error estimates for the SXS waveforms based on the
extrapolation and the resolution uncertainties of the gravitational wave
strain, which dominate over the errors obtained from the quasilocal measures
of the final mass and spin. Finally, we observe substantial instabilities on
the best-fit amplitudes of the tones beyond the fundamental mode and the first
overtone, that, nevertheless, do not impact significantly the mass and spin
estimates.

##  I Introduction

The number of gravitational wave (GW) observations is increasing along with
the upgrades of GW interferometers. Up to date, the LIGO-Virgo collaboration
has reported a total of 48 binary black hole merger candidates Abbott  _et
al._ (2019, 2021a). Those observations are providing unprecedented constraints
on general relativity in its strong-field regime, with the merger-ringdown
phase in particular providing a promising channel for such studies.

A binary black hole merger is generally decomposed in three different regimes
that depict its orbital evolution: inspiral, merger and ringdown. The inspiral
regime represents the _slow_ far-field solution and it is well described by
post-Newtonian and effective-one-body theories. At the merger phase, that is,
when the two bodies get closer to each other, these analytic solutions break
down due to the strong general-relativistic effects and full numerical
relativity is needed. The final merger results in an initially perturbed
space-time that evolves towards the final Kerr solution by emitting an
ultimate tail of radiation better known as ringdown (RD). The strain
$h(t,\theta,\phi)$ of the RD waveform is predicted by linear perturbation
theory to decompose as a sum of damped sinusoids:

$\displaystyle
h(t,\theta,\phi)=\sum_{l,m,n}{\mathcal{A}}_{lmn}e^{-\iota\omega_{lmn}(t-t_{r})}\,{}_{-2}\mathcal{Y}_{lm}(\theta,\phi)\,,\,t\geq
t_{r}\,.$ (1)

Here, $l=2,3,\dots$ and $m=-l,-l+1,\dots,l-1,l$ account for the two angular
indices of the spheroidal decomposition, while $n=0,1,2,\dots$ labels the
tone; ${}_{-2}\mathcal{Y}_{lm}(\theta,\phi)$ are the spin-weighted spheroidal
harmonics of spin weight $s=-2$, as functions of the polar angle $\theta$ and
azimuthal angle $\phi$; ${\mathcal{A}}_{lmn}=A_{lmn}\,e^{\iota\varphi_{lmn}}$
is the tone complex amplitude; and $t_{r}$ is some undefined time beyond which
linear perturbation theory is expected to accurately describe the RD regime
Giesler  _et al._ (2019); London  _et al._ (2014); Bhagwat  _et al._ (2020).
In particular, and for non-charged black holes, the
$\omega_{lmn}=w_{lmn}-\iota/\tau_{lmn}$ defines a countably infinite set of
complex frequencies solely determined by the final black hole’s mass $M_{f}$
and spin $a_{f}$, where the values of $\omega_{lmn}$ correspond to poles of
the Green function to the inhomogeneous Teukolsky equation — the quasinormal
modes (QNMs) of the final black hole Leaver (1985); Detweiler (1980); Kokkotas
and Schmidt (1999). Here $\mathrm{Re}[\omega_{lmn}]=w_{lmn}$ and
$-\mathrm{Im}[\omega_{lmn}]=1/\tau_{lmn}$ take the role of the oscillation
frequency and the damping rate (inverse of the damping time) respectively. As
a rule of thumb, if one considers fixed the value of the $(l,m)$ indices, the
mass $M_{f}$, and for moderate spins $a_{f}$, the values of the damping times
$\tau_{lmn}$ decrease as the tone index $n$ increases. This sets the $n=0$
(fundamental) tone as the dominant tone while the $n\geq 1$ tones (overtones)
rank down continuously as $n$ increases. Moreover, one finds two branches of
solutions for $\omega_{lmn}$ also known as the corotating (dominant for $m>0$)
and counter-rotating (subdominant for $m>0$) modes Finch and Moore (2021);
Dhani (2020); Berti  _et al._ (2009); Jiménez Forteza  _et al._ (2020); Cook
(2020); Dhani and Sathyaprakash (2021) which both contribute111Such
contribution involves an additional sum over a binary index labeling the co-
and counter-rotating QNM frequencies and amplitudes which have been dropped
here to simplify the notation. In the following we shall rather explicitly
state, when needed, whether a given QNM is associated to a co- or counter-
rotating branch. Note that the counter-rotating modes excited in a binary
black hole merger are usually expected to have negligible amplitudes compared
to the corotating modes (for $m>0$ harmonics) Jiménez Forteza  _et al._
(2020); Finch and Moore (2021). to Eq. (1). Co- and counter-rotating modes are
distinguished by the sign of the real part $w_{lmn}$ of their frequencies.
While multiple conventions exist, in this work we follow the convention of,
_e.g._ , Berti  _et al._ (2009); Cook and Zalutskiy (2014), and we denote as
corotating the family of modes with positive $w_{lmn}$ regardless of the sign
of $m$ or of the sign attributed to the black hole’s spin $a_{f}$.

The black hole no-hair and uniqueness theorems in general relativity imply
that the final state of an uncharged black hole merger, and the associated QNM
spectrum, are uniquely determined by the values of the final mass and spin.
This has led to two main avenues to test such theorems. The first one consists
on performing an inspiral–merger–ringdown (IMR) consistency test, which relies
on independently estimating the final black hole mass and spin from both the
inspiral–merger and the ringdown phases Abbott  _et al._ (2021b). The second
approach is to perform black hole spectroscopy, which typically aims at
independently estimating the parameters of the fundamental tone of the
dominant angular mode, $(l=2,m=2,n=0)$, plus at least another mode, either i)
the first corresponding overtone, $(l=2,m=2,n=1)$, or ii) another angular
fundamental mode, either the $(l=m=3,n=0)$ or the $(l=2$, $m=1$, $n=0)$ mode
(in order of importance). So far and for unequal-mass-ratio binaries, the
higher angular mode remains the most promising approach to test the
implications of the black hole no-hair theorem Jiménez Forteza  _et al._
(2020). Successful independent evidence of the $(l=m=2,n=0)$ and the
$(l=m=3,n=0)$ modes in the ringdown phase of a GW event (in this case
GW190521) has been recently provided in Capano  _et al._ (2021). On the other
hand, channel i) becomes a promising possibility when dealing with near equal-
mass-ratio nonspinning binaries. For such events, the higher harmonic modes
are only weakly excited, while the overtones would still represent a valid
channel in the ringdown regime. A first attempt to observe overtones in GW
observational data has been performed in Isi  _et al._ (2019); Isi and Farr
(2021) on GW150914. However, the full spectroscopic analysis performed by
Capano  _et al._ (2021) on GW190521 could not find evidence of tones other
than the fundamental ones.

Current studies of the $(22n)$ ringdown modes rely on fits to numerical
relativity (NR) waveforms Bhagwat  _et al._ (2020); Jiménez Forteza  _et al._
(2020); Giesler  _et al._ (2019); London and Fauchon-Jones (2019), which are
shown to be consistent with current GW observations. In particular, using NR
waveforms has the following advantages: i) the underlying theory is well-
known; ii) the mass and the spin of the final black hole (BH) are accurately
estimated, hence accurately determining the QNM spectrum; and iii) numerical
errors in the simulated waveforms are typically smaller than current GW
detectors noise. In such studies, considering the $(l=m=2)$ spherical harmonic
of the strain222It is worth mentioning here that in NR codes, the strain
$h(t,\theta,\phi)$ is decomposed in terms of the spin-weighted spherical
harmonics basis instead of the spheroidal harmonics $\mathcal{Y}_{lm}$ used to
define QNMs, since it is a better adapted basis to the inspiral-merger regime.
This adds mode-mixing artifacts between both bases principally at modes other
than the $(22)$ mode Berti and Klein (2014); Cook (2020)., $h_{22}(t)$, from a
given time $t_{0}$ onwards, one fits for the successive complex amplitudes
$\mathcal{A}_{22n}$ of the $(22n)$ QNM tones for a running index
$n\in\\{0,\dots,N_{\mathrm{max}}\\}$, with various choices for the total
number $N_{\mathrm{max}}$ of overtones to be included in the model. The
$N_{\mathrm{max}}=7$ model has been shown to provide the best estimates of the
true final parameters (mass and spin) Giesler  _et al._ (2019); Finch and
Moore (2021), although no models beyond $N_{\mathrm{max}}=7$ have been studied
up to date. In this work we extend this analysis to $N_{\mathrm{max}}\geq 7$.
In particular, we have not found any publicly available catalog of Kerr QNM
data that provides a correct description for the $n=8$ tone nor for the
neighboring corotating tone that we label as $n=9$ here333The catalogs Berti
_et al._ (2006, 2009); Berti ; Cardoso  and Cook and Zalutskiy (2014); Cook
(2019) provide the QNM solutions up $n=7$. On the other hand, the solutions
provided by the qnm Python package Stein (2019) up to much larger $n$ values
are incorrect at ${n=8}$ (both for the co- and counter-rotating modes) due to
the erroneous estimate of the Schwarzschild limit, and are missing the
neighboring corotating branch that we label here as $n=9$. Solutions for these
$n=8,9$ modes have been previously obtained by Onozawa (1997) ($n=8$ modes
only) and Berti  _et al._ (2003); Cook and Zalutskiy (2014, 2016a, 2016b)
(with a different labeling of what we here call the $n=8,9$ corotating modes
as two corotating $n=8$ branches), where they are shown as
$\mathrm{Re}(\omega_{lmn})$ — $\mathrm{Im}(\omega_{lmn})$ frequency plots but
the data obtained for $n\geq 7$ modes were not made publicly available. . We
discuss these tones further in Sec. II below and show our results for their
frequencies in Sec. II.3, comparing to the results from Onozawa (1997); Berti
_et al._ (2003); Cook and Zalutskiy (2014, 2016a). In Sec. III we revisit the
definitions of the mass and spin in NR simulations and we provide two methods
to compute their uncertainties. Finally, in Sec. IV, we show the results on
the mass and the spin obtained from fitting models with up to
$N_{\mathrm{max}}=12$ or $16$ overtones to NR waveforms from the SXS and RIT
catalogs The SXS Collaboration (2016); Campanelli  _et al._ (2016) while we
further discuss the model instabilities.

##  II The ringdown QNM spectrum

###  II.1 The ringdown wave equation

The Teukolsky master equation Teukolsky and Press (1974) describes the
propagation of linear perturbations of fields of general spin weight $s$ in a
Kerr background Chandrasekhar (1985); Berti  _et al._ (2009); Kokkotas and
Schmidt (1999). The angular ($\mathcal{Y}_{lm}$) and radial ($R_{lm}$) sector
of this equation read, respectively,

$\displaystyle\frac{\partial}{\partial
u}\left[\left(1-u^{2}\right)\frac{\partial}{\partial
u}\,{}_{s}\mathcal{Y}_{lm}\right]+\left[(a_{f}\,\omega\,u)^{2}\phantom{\frac{\left(m+su\right)^{2}}{1-u^{2}}}\right.$
$\displaystyle{}-$ $\displaystyle
2\,a_{f}\,\omega\,s\,u+s+\left.{\mathscr{A}}-\frac{\left(m+su\right)^{2}}{1-u^{2}}\right]{}_{s}\mathcal{Y}_{lm}=0\,;$
$\displaystyle\Delta\,\partial^{2}_{r}R_{lm}+(s+1)(2r-2M_{f})\partial_{r}R_{lm}+VR_{lm}=0\,,$
(2)

where ${s=\pm 2}$ for gravitational perturbations. Here $u=\cos\theta$;
$\Delta=(r-r_{-})(r-r_{+})$, $r$ is the radial coordinate while $r_{+,-}$
stand for the coordinate radii of the outer and inner BH horizons
respectively; $l,m$ are the usual angular indices; $a_{f}$ is the black hole’s
spin and $M_{f}$ the black mass; $V=V(r,M_{f},a_{f},\omega,\mathscr{A},s,m)$
is the potential term for a Kerr BH (see Eq.(26) of Berti  _et al._ (2009));
$\omega$ is the complex frequency of the perturbation; and $\mathscr{A}$ is
the corresponding so-called angular separation constant.

For each value of the final spin $a_{f}$, each spin weight $s$ and each
angular mode $(l,m)$, the $(l,m,n)$– quasinormal modes are obtained by
imposing outgoing boundary conditions at spatial infinity and ingoing boundary
conditions at the black hole horizon. The QNMs and associated angular
separation constants form the only discrete set of (complex) values
$\left\\{\omega,\mathscr{A}\right\\}$ that are compatible with these boundary
conditions. Eq. (2) can be solved and its associated QNMs values obtained
following the algorithm proposed in Leaver (1985). The frequency and
separation constant solutions are then labeled by the integers $l$, $m$ and
$n$: $\omega\equiv\omega_{lmn}$, $\mathscr{A}\equiv\mathscr{A}_{lmn}$, where
$n=0,1,2,\ldots$ is the overtone index. The dependence on the spin weight $s$
usually remains implicit; we only consider gravitational perturbations here
and we set $s=-2$ throughout this work. See Leaver (1985) for a method for
numerically calculating the QNM spectrum, Berti ; Cardoso ; Berti  _et al._
(2006, 2009) and Cook (2019); Cook and Zalutskiy (2014) for a compilation of
the values in different situations and up to $n=7$ for the Kerr scenario, and
Stein (2019) for a Python package, qnm, to evaluate the QNM spectrum of Kerr
black holes for a variety of $(l,m,n)$ modes and spin weights $s$.

###  II.2 Computing the ($22n$) quasinormal modes

Most of the QNM frequency values used for this work (which all correspond to
the $(l=2,m=2)$ harmonic) were computed using the dedicated qnm Python package
Stein (2019). However, the method used (Leaver’s method Leaver (1985)) is
known to fail for the $(l=m=2,n=8)$ mode in the Schwarzschild limit as it
becomes an algebraically special mode Chandrasekhar and Detweiler (1975);
Berti  _et al._ (2003). Hence, this tone was flagged as unreliable in this
code since the spectrum computation relies on the Schwarzschild limit. Indeed,
the results from qnm at any spin for this mode appear to be inconsistent with
the neighboring modes, and are in disagreement with Onozawa (1997); Berti  _et
al._ (2003); Cook and Zalutskiy (2014), for both the co- and counter-rotating
branches. Hence, we rather computed the QNM frequencies $\omega_{lmn}$ — along
with the angular separation constants $\mathscr{A}_{lmn}$ — for this mode from
a modified version of the publicly available Mathematica code for Kerr QNMs
from Berti  _et al._ (2006, 2009) (available online at Berti ; Cardoso ).

This code is also based on Leaver’s continued fraction method Leaver (1985) —
but it may be used to directly compute the QNMs for any given spin, without
relying on the Schwarzschild limit. Following this method, estimates of
$\omega_{lmn}$ and $\mathscr{A}_{lmn}$ are found successively as roots of
infinite generalized continued fractions, which are approximated with a finite
numbers of fraction steps $n_{\mathrm{frac}}$. The coefficients involved in
the fraction used for each of both variables depend on the estimate of the
other variable, thus the alternated estimation of $\omega_{lmn}$ and
$\mathscr{A}_{lmn}$ is iterated until convergence is reached. We have modified
the continued fraction computation to include Leaver’s inversions (Eq. (14) in
Leaver (1985)) allowing for a more stable recovery of any given overtone
($n\geq 1$), and replaced the use of Mathematica’s time- and memory-consuming
root-finding algorithm by a direct implementation of the secant method. We
ensure the convergence both in terms of $n_{\mathrm{frac}}$ and of the
$(\omega_{lmn},\mathscr{A}_{lmn})$ loop by increasing $n_{\mathrm{frac}}$ by a
constant factor $c_{\mathrm{frac}}>1$ at each iteration, until a convergence
criterion is met. This criterion amounts to requiring that the (absolute)
variations of the estimates of both variables over three consecutive
iterations do not exceed a certain threshold (which we set at $3\cdot
10^{-11}$). This progressive increase of $n_{\mathrm{frac}}$ was necessary for
the computation of the $n=8$ modes (as well as the neighboring corotating
branch which we label here as $n=9$), as these modes typically require rather
large values of $n_{\mathrm{frac}}$ (further increasing as the spin gets
closer to zero or one) to reach such an accuracy; we set
$c_{\mathrm{frac}}=1.2$ for this computation. We have checked for consistency
that this algorithm provides the same $s=-2$ solutions as those for available
$(l=2,m=\pm 2,n\neq 8)$ tones from Berti  _et al._ (2006, 2009); Berti ;
Cardoso ; Stein (2019).

We make this modified Mathematica code available here: Mourier and Jimenez
Forteza ; Jimenez Forteza and Mourier , along with its translation into Python
and into Fortran. We also made use of this (much faster) Fortran equivalent,
to obtain the frequency solutions for these modes at a few spin values in the
regimes where convergence to the solution is particularly difficult to
achieve.

The algorithm does require an initial guess for $\omega_{lmn}$444More
precisely, since we modified the code to use the secant method rather than
Newton’s method to find the roots of the continued fractions, two initial
guesses on $\omega_{lmn}$ are required instead of one. These may simply be
chosen as two close yet distinct estimates such as bounds on the expected
solution or simply perturbations around a given estimate. This is also
required for initial guesses on $\mathscr{A}_{lmn}$. For this variable we used
systematic small deviations above and below the single guess value that was
originally used to initialize Newton’s method — that is, at each iteration
beyond the first, the $\mathscr{A}_{lmn}$ result from the previous iteration,
and at the first iteration, the Schwarzschild-limit solution
$\mathscr{A}_{lmn}=l(l+1)-s(s+1)$. , which simply needs to lie closer to the
desired mode than to any other tone of the same $(l,m)$ harmonic — in practice
for the modes discussed here, a $\sim 5\%$ accurate initial estimate is
typically sufficient. This allowed us to recover the $n=8$ co- and counter-
rotating QNMs over a wide spin range (see below) by simply using initial
guesses based on a few points of the $n=8$ curves in Fig. 4 of Onozawa (1997)
and interpolation and extrapolation between and beyond them. Using slightly
lower imaginary values for the initial guesses on the frequency with respect
to the above corotating solution, we also recovered the additional corotating
mode that also nears the imaginary axis at low spin identified by Berti  _et
al._ (2003); Cook and Zalutskiy (2014) and missing in Onozawa (1997) and in
the qnm package.

#### Remarks on the tone labeling convention

We label in this work this additional corotating branch as the $n=9$ mode and
the subsequent ones (with even smaller values of $\mathrm{Im}[\omega_{lmn}]$)
as $n=10,11,\dots$ The counter-rotating mode associated with each of the
$n\neq 8,9$ corotating tones, connecting to the same Schwarzschild limit up to
a $\mathrm{Re}(\omega)\mapsto-\mathrm{Re}(\omega)$ symmetry, is attributed the
same $n$ index. This overall leads to an unusual convention for $n\geq 8$ and
can be somewhat confusing: our $(l=2,m=2,n=n_{0})$ co- and counter-rotating
QNMs for each $n_{0}\geq 10$ are equivalent to the $(l=2,m=2,n=n_{0}-1)$
solutions for the qnm package Stein (2019) or from Cook and Zalutskiy (2016a)
for instance, and in the Schwarzschild limit, they match the Schwarzschild QNM
that is usually attributed the overtone index $n_{0}-1$ in the literature.
This offset is due to the presence of a single Schwarzschild ($l=2,m=2$) QNM,
traditionally labeled $n=8$, at $\omega=-2\,\iota$ — which coincides with an
algebraically special mode Maassen van den Brink (2000); Cook and Zalutskiy
(2014) — while two distinct Kerr QNMs are found near this value at low spins
Berti  _et al._ (2003); Cook and Zalutskiy (2014). This is reconciled in Berti
_et al._ (2003) by simply considering these two Kerr solutions as a double
$n=8$ branch. The same choice is made in Cook and Zalutskiy (2014, 2016a,
2016b) where, specifically, these Kerr modes that we here refer to as $n=8$
and $n=9$ overtones are labeled as the $n=8_{0}$ and $n=8_{1}$ modes,
respectively. On the other hand, the convention that we adopt here allows for
a consistent sorting of the found Kerr QNMs by decreasing imaginary frequency
(or decreasing damping time) for a given spin, and for the preservation of a
roughly equal spacing between $\mathrm{Im}[\omega_{l,m,n}]$ and
$\mathrm{Im}[\omega_{l,m,n+1}]$ for corotating modes for all values of $n$,
over most of the spin range.

This latter property could in fact be used to get reasonable initial guesses
for the frequency values of both the $n=8$ and $n=9$ corotating modes over
most of the spin range ($a_{f}\gtrsim 0.1$), without prior knowledge about
these values, by extrapolating the neighboring tones’ frequencies as a
function of $n$ for each $a_{f}$. In particular, we did generate initial
guesses in this way to systematically compute the $n=9$ mode frequencies for
spins $a_{f}\geq 0.1$. Initial guesses for lower spin values for this mode
were obtained by successive extrapolations of the results previously obtained
at larger spins.

###  II.3 The $\omega_{22n}$ curves

In Fig. 1 we show the range of $\omega_{lmn}$ values for the $(l,m)=(2,2)$ and
$n\in\left\\{6,\dots,10\right\\}$ corotating (solid curves) and counter-
rotating modes (dashed curves) for a unit-mass ($M_{f}=1$) Kerr black hole as
its dimensionless spin $a_{f}$ varies. These curves correspond to the full
range of spin values $a_{f}\in[0,1]$ for $n\not=8,9$, and to the ranges over
which we could obtain solutions for the $n=8,9$ tones:
$a_{f}\in\left[a_{f,n}^{\mathrm{start}},1\right]$ for the corotating modes,
with $a_{f,8}^{\mathrm{start}}\equiv{\color[rgb]{0,0,0}3.5}\cdot 10^{-3}$ for
$n=8$ and $a_{f,9}^{\mathrm{start}}\equiv{\color[rgb]{0,0,0}5.4}\cdot 10^{-3}$
for $n=9$, and $a_{f}\in{\color[rgb]{0,0,0}\left[0,0.997\right]}$ for the
associated counter-rotating mode. The counter-rotating mode curves shown have
actually been flipped around the imaginary axis (using the symmetry
transformation $\omega\mapsto-\omega^{*}$ where ∗ stands for the complex
conjugation) for better visualization. This choice can alternatively be
interpreted as a representation of the corresponding corotating solutions
either for the $(l=2,m=-2)$ mode (as in Berti  _et al._ (2003); Cook and
Zalutskiy (2014, 2016a, 2016b)) or for negative spins, from the following
symmetry relations Berti  _et al._ (2009); Cook and Zalutskiy (2014):

$\displaystyle\omega_{lmn}\left(a_{f}\right)=-\omega_{l(-m)n,c}^{*}\left(a_{f}\right)\,;$
(3)
$\displaystyle\omega_{lmn}\left(a_{f}<0\right)=-\omega_{lmn,c}^{*}\left(|a_{f}|\right)\,,$
(4)

where $\omega_{{l(\pm m)n},c}$ stand for the counter-rotating mode
frequencies. The negative-spin interpretation explains that most of these
curves continuously connect to the ($a_{f}\geq 0$) corotating branches. The
Schwarzschild limit for each $n\neq 8,9$ mode then appears at the transition
point between dashed and solid curves in the figure and is marked with a dot;
while the extremal-Kerr limit $a_{f}\rightarrow 1$ is found as
$\left\\{\mathrm{Re}[\omega_{lmn}],\mathrm{Im}[\omega_{lmn}]\right\\}\rightarrow\left\\{1,0\right\\}$
for all of the corotating modes shown here (_cf._ Chandrasekhar and Detweiler
(1975)).

For the corotating branch, $\mathrm{Re}[\omega_{lmn}]$ increases monotonically
with the spin of the final object $a_{f}$. Similarly,
$\mathrm{Im}[\omega_{lmn}]$ increases with $a_{f}$ over most of the spin range
for $n\geq 9$ and over the whole range for $n\leq 8$. Moreover, one can notice
that the values of the corotating modes span a broader domain (both in
$\mathrm{Re}[\omega_{lmn}]$ and $\mathrm{Im}[\omega_{lmn}]$) with respect to
the counter-rotating branch. The same features apply as well for all tones
other than those shown, except for $n=5$ which has a different high-spin
behavior Onozawa (1997). Regarding the counter-rotating branches, the
solutions decrease monotonically with $a_{f}$ both in real and imaginary part
up to $n=6$. For $n\geq 7$ the solutions become degenerate in
$\mathrm{Re}[\omega_{lmn}]$, where one can find multiple values of
$\mathrm{Im}[\omega_{lmn}]$ given one fixed $\mathrm{Re}[\omega_{lmn}]$.

Figure 1: Co- and counter-rotating (CR) QNM frequencies $\omega_{lmn}$ on the
complex plane for a final black hole mass set to $M_{f}=1$, for
$n\in\left\\{6,\dots,10\right\\}$ and for spin $a_{f}$ spanning $[0,1]$ or a
wide subset of this range (as discussed in the main text, Sec. II.3). The
solid lines account for the corotating solutions while the dashed ones
correspond to the counter-rotating tones. The latter are represented under the
$\omega\mapsto-\omega^{*}$ transformation, or equivalently, as per Eqs.
(3)–(4), as corotating modes with either $m=2\mapsto m=-2$, or, $a_{f}\mapsto-
a_{f}$. For $n\neq 8,9$, the Schwarzschild limit is recovered at the joining
point between the dashed and solid lines, and is materialized by a dot. At the
extremal-Kerr limit for all ($a_{f}>0$) corotating modes shown here,
$\omega_{lmn}\rightarrow 1$ as expected. The roughly equal spacing in
$\mathrm{Im}(\omega_{22n})$ between successive ($a_{f}>0$) corotating tones
can be noted on these curves for $\mathrm{Re}(\omega_{22n})\gtrsim 0.2$,
corresponding to spins $a_{f}\gtrsim 0.3$. This holds down to smaller spins
$a_{f}\gtrsim 0.1$ up to small shifts in $\mathrm{Re}(\omega_{22n})$ at fixed
$a_{f}$ between the tones shown here.

While the other modes were computed using the qnm Python package, as mentioned
above the $n=8$ (co- and counter-rotating) and $n=9$ curves that we present in
this figure have been obtained from the adapted version of the Mathematica
code555We also used this code to complete the curves for the other corotating
modes up to spins close to $1$, as the qnm package results become unreliable
—typically swapping different tones— at very high spins ($a_{f}\gtrsim$
0.995).  from Berti  _et al._ (2006, 2009); Berti ; Cardoso  described in Sec.
II.2. Our results for these three modes are in good qualitative agreement with
those shown in Onozawa (1997) (limited to the $n=8$ co- and counter-rotating
modes), Berti  _et al._ (2003), and Cook and Zalutskiy (2014). Their low-spin
behavior discussed in more detail below is well compatible with the higher-
accuracy investigation of Cook and Zalutskiy (2014) in particular — while the
low-spin range is too limited in Onozawa (1997); Berti  _et al._ (2003) to
unambiguously compare the trends.

We provide our complex frequencies and angular separation constants results
for these modes as three tables, each one listing a range of $|a_{f}|$ values
with the corresponding $\mathrm{Re}[\omega_{lmn}]$,
$\mathrm{Im}[\omega_{lmn}]$, $\mathrm{Re}[\mathscr{A}_{lmn}]$ and
$\mathrm{Im}[\mathscr{A}_{lmn}]$ for $M_{f}$ set to $1$. We make these tables
available here: Mourier and Jimenez Forteza ; Jimenez Forteza and Mourier ,
along with the codes that we used. Given the convergence criterion mentioned
above, we consider each of these values to be accurate up to the number of
digits provided, that is, to an absolute precision of $10^{-10}$ for each of
these quantities. These results complement the data of Berti  _et al._ (2006,
2009); Berti ; Cardoso  and they correct and complement the results of the qnm
package Stein (2019) for these three branches. We use these tables, along with
the qnm package for all other modes, to produce the fits described in Section
IV.

The corotating modes in these tables are provided for the ranges
$a_{f}\in[a_{f,n}^{\mathrm{start}},1]$ mentioned above, with a step on $a_{f}$
set to $\delta a_{f}=10^{-4}$ and further refined close to $a_{f}=1$ ($\delta
a_{f}=10^{-5}$ for $0.9990\leq a_{f}\leq 0.9999$ and $\delta a_{f}=10^{-6}$
for $a_{f}\geq 0.9999$). For either mode, convergence was extremely slow at
and in the vicinity of its respective $a_{f,n}^{\mathrm{start}}$ spin value,
preventing the investigation of a large number of $a_{f}$ values below this
point or a decrease of the spin step near this point. This difficult
convergence is likely a consequence the known failure of Leaver’s method in
the vicinity of the algebraically special Schwarzschild mode
$\omega=-2\,\iota$ Chandrasekhar and Detweiler (1975); Berti  _et al._ (2003),
where these modes lie at low spins. We could not achieve any convergence —
even at a much lower precision level — for the few
$a_{f}<a_{f,n}^{\mathrm{start}}$ values that we probed (_e.g._ , at
$a_{f}=a_{f,n}^{\mathrm{start}}-\delta a_{f}$ and
$a_{f}=a_{f,n}^{\mathrm{start}}-2\,\delta a_{f}$ in both cases with $\delta
a_{f}=10^{-4}$, at $a_{f}=3\cdot 10^{-3}$ for $n=8$, or at $a_{f}=5\cdot
10^{-3}$ for $n=9$). In both cases,
$\mathrm{Re}\big{[}\omega_{lmn}(a_{f}=a_{f,n}^{\mathrm{start}})\big{]}$ is
very close to zero, and extrapolating the frequency solutions to lower spins
would make them cross the imaginary axis at a finite spin value
$a_{f,8}^{\mathrm{start}}-2\,\delta
a_{f}<a_{f}<a_{f,8}^{\mathrm{start}}-\delta a_{f}$ for $n=8$ and at a finite
spin value $a_{f,9}^{\mathrm{start}}-\delta
a_{f}<a_{f}<a_{f,9}^{\mathrm{start}}$ for $n=9$, with the imaginary part of
the frequency remaining distinct from $-2$ at the crossing point in both
cases.

While the lack of convergence could simply be due to a complete failure of the
method in this range, these results — including the extrapolated values of
$a_{f}$ and $\mathrm{Im}(\omega_{22n})$ where the imaginary axis would be
crossed, although we find them with lower accuracy — are fully consistent with
the findings of Cook and Zalutskiy (2014). With an investigation extended even
closer to the imaginary axis, Cook and Zalutskiy (2014) indeed found both
branches to reach the axis at a finite spin and away from the algebraically
special Schwarzschild mode (with $\omega_{22n}\simeq-1.96384\,\iota$ at
$a_{f}\simeq 3.4826\cdot 10^{-3}$ and $\omega_{22n}\simeq-2.04223\,\iota$ at
$a_{f}\simeq 5.3279\cdot 10^{-3}$ for $n=8,9$, respectively), and to disappear
at lower spins. Cook and Zalutskiy (2016a) additionally showed that these QNMs
do not exist either _on_ the imaginary axis itself.

For the counter-rotating mode that we present in a third table (more
precisely, this table corresponds to the corotating $n=8$ mode for $m=2$,
$a_{f}<0$ or for $m=-2$, $a_{f}>0$, tied to the $m=2$, $a_{f}>0$, $n=8$
counter-rotating mode by the symmetry relations (3)–(4)), convergence was also
slower at spins very close to $0$ but could still be achieved down to
$|a_{f}|=10^{-6}$ — yet not at the Schwarzschild $a_{f}=0$ limit itself, as
expected. The values we obtain for $|a_{f}|>0$ are however compatible with the
Schwarzschild algebraically special limit $\omega_{22n}\rightarrow-2\,\iota$
(along with the $(m,n)$–independent Schwarzschild limit
$\mathscr{A}_{22n}\rightarrow 4$) for $a_{f}\rightarrow 0$, to within less
than $5\cdot 10^{-10}$ by extrapolation. This is in agreement with the
analytical prediction of this limit for the counter-rotating mode by Maassen
van den Brink (2000), while the solution also obtained down to
$|a_{f}|=10^{-6}$ in Cook and Zalutskiy (2014) was also compatible with it. We
accordingly assumed the validity of this limit and added it to the table at
$a_{f}=0$.

However, for this mode, convergence was much poorer at high spin values. We
accordingly provide the results for this mode over the spin range
$|a_{f}|\in{\color[rgb]{0,0,0}[0,0.997]}$ (the values at $a_{f}=0$ being
assumed as mentioned above), with a step $\delta a_{f}=10^{-6}$ at
$|a_{f}|\leq 2\cdot 10^{-5}$, $\delta a_{f}=10^{-5}$ at $2\cdot
10^{-5}\leq|a_{f}|\leq 2\cdot 10^{-4}$, $\delta a_{f}=10^{-4}$ at $2\cdot
10^{-4}\leq|a_{f}|\leq 0.97$, and $\delta a_{f}=10^{-3}$ at
$0.97\leq|a_{f}|\leq 0.997$. For both of the next spin values $|a_{f}|=0.998$
and $|a_{f}|=0.999$, convergence could not be reached, with the algorithm even
appearing to be divergent in the second case.

The disappearance of the $n=8,9$ corotating modes — or at least the failure of
the algorithm — at low spin values prevents the association of the counter-
rotating mode to either corotating branch. On the other hand, like the
previous studies, we have found only this single counter-rotating solution for
any given spin value in this region of the complex plane. Hence, as for the
$\omega=-2\,\iota$ Schwarzschild QNM, we associate this counter-rotating
branch to both corotating ones. Since we label these two corotating branches
as the $n=8$ and $n=9$ tones, the associated Schwarzschild mode and counter-
rotating branch may be considered as a degenerate $n=8$ and $n=9$ mode
simultaneously.

##  III The waveform catalog

(a)

(b)

Figure 2: Left panel: Distribution of the final mass $M_{f}$ and final spin
$a_{f}$ of the 610 SXS waveforms used in this work. Right panel: same
distribution but for the effective spin $\chi_{\mathrm{eff}}$ and the mass
ratio $q$ instead.

In this work, we have used 620 non-precessing waveforms from the SXS catalog
The SXS Collaboration (2016), and two from the RIT catalog Campanelli  _et
al._ (2016) for comparison. We excluded $\fpeval{620-610}$ out of the 620 SXS
cases, which we did not consider accurate enough for our analysis or which had
seemingly inconsistent final parameters (see Appendix A). We show in Fig. 2
the parameter space corresponding to the 610 SXS simulations analysed. The
left panel shows the $(M_{f},a_{f})$ distribution, with
$M_{f}\in[0.883,0.993]$ and $a_{f}\in[-0.344,0.997]$. Here and in the
following, we make the masses dimensionless by setting the total initial mass
of the two merging black holes $M=m_{1}+m_{2}$ to unity. On the right panel,
we alternatively show the distribution of the same SXS setups in terms of the
mass ratio $q=m_{1}/m_{2}\in[1,10]$ and the effective spin
$\chi_{\mathrm{eff}}=(\chi_{1}m_{1}+\chi_{2}m_{2})/(m_{1}+m_{2})\in\left[-0.97,0.9988\right]$.
The visible correlation between $M_{f}$ and $a_{f}$ is physical: at fixed mass
ratio, the relative energy radiated, $E_{f}=1-M_{f}$, increases with the value
of $\chi_{\mathrm{eff}}$ which is itself correlated with $a_{f}$ Rezzolla  _et
al._ (2008); Jiménez-Forteza  _et al._ (2017); Hofmann  _et al._ (2016). Among
all the cases analysed here, there are 18 cases with $a_{f}<0$, for which we
have used the symmetries given by Eq. (4).

###  III.1 Estimates of the final mass and spin

The final mass and final spin provided by NR catalogs are usually estimated
from their quasilocal definitions The SXS Collaboration (2016); Szilagyi  _et
al._ (2009); Ashtekar and Krishnan (2004); Iozzo  _et al._ (2021) on the
apparent horizon (AH). In the ringdown regime, the distorted black-hole
spacetime evolves quickly towards its stationary state. Then, the black hole
spin $S$ is obtained by computing the set of approximate Killing vectors
$\phi^{i}_{(k)}$ and the extrinsic curvature $K_{ij}$ at the AH and
integrating them over the induced AH area as,

$S_{\phi_{(k)}}=\frac{1}{8\pi}\int_{\mathrm{AH}}\phi_{(k)}^{i}\,s^{j}K_{ij}\,dA\,,$
(5)

where $s^{j}$ is the outgoing spacelike unit normal vector and
$S_{\phi_{(k)}}$ is the spin component estimate along the Killing vector
$\phi^{i}_{(k)}$. The spin magnitude is then evaluated as
$S\equiv\sqrt{S_{\phi_{(1)}}^{2}+S_{\phi_{(2)}}^{2}+S_{\phi_{(3)}}^{2}}$. The
final mass relies on the spin value $S$ and it is obtained by using the
Christodoulou formula for uncharged black holes Christodoulou (1970),

$\left(M_{f}^{\mathrm{l}}\right)^{2}=M^{2}_{\mathrm{irr}}+\frac{S^{2}}{4M_{\mathrm{irr}}^{2}}\,,$
(6)

which depends on the value of the final spin $S$ and the irreducible mass
$M_{\mathrm{irr}}$ (for further details see Sec.2.2 of The SXS Collaboration
(2016)). We work with the local dimensionless spin, namely,

$a_{f}^{\mathrm{l}}=\frac{S}{(M_{f}^{\mathrm{l}})^{2}}\,.$ (7)

The superscripts l stand for quasilocal mass and spin. Alternatively, the mass
and the spin can also be estimated from the energy and angular momenta
radiated away in the form of gravitational radiation. These _radiation_ -based
quantities (labeled with a superscript r) are obtained in terms of the Newman-
Penrose scalar $\psi_{4}$,

$\displaystyle\begin{split}M_{f}^{\mathrm{r}}&={\color[rgb]{0,0,0}M_{\mathrm{in}}}-\lim_{r\rightarrow\infty}\frac{r^{2}}{16\pi}\int\sum_{lm}\left|\frac{dh_{lm}}{dt}\right|^{2}\,dt^{\prime}\,;\\\
a_{f}^{\mathrm{r}}&=J_{\mathrm{in}}+\lim_{r\rightarrow\infty}\frac{r^{2}}{16\pi}\,\mathrm{Re}\\!\left[\int\sum_{lm}m\,h^{*}_{lm}\,\frac{d(h_{lm})}{dt}\,dt^{\prime}\right]\,,\end{split}$
(8)

where $h$ is the gravitational wave strain, $\psi_{4}=d^{2}h/dt^{2}$, and
$J_{\mathrm{in}}$ and $M_{\mathrm{in}}$ are the initial ADM dimensionless
angular momentum and initial ADM energy, respectively. To dissipate the local
gauge effects, the radiation quantities are evaluated at a distance
$\mathcal{O}(100M)$ away from the black holes apparent horizons (which have a
$\mathcal{O}(1M)$ radius) and extrapolated to null infinity. The integrals are
evaluated starting at a time $t^{\prime}=t_{\mathrm{in}}$ which is always
taken to be later than the emission of the junk radiation.

###  III.2 Resolution and extrapolation errors

Usually, the local grid the near-horizon zone of binary BH simulations is
better resolved that the radiation zone so that $a_{f}^{\mathrm{l}}$ and
$M_{f}^{\mathrm{l}}$ are estimated to larger accuracy than
$a_{f}^{{\mathrm{r}}}$ and $M_{f}^{{\mathrm{r}}}$. Moreover, the radiative
quantities are also affected by extrapolation errors when extrapolating from
$r=\mathcal{O}(100M)$ to null infinity, by conversion errors from $\psi_{4}$
to $h$ or by still non-zero residual gauge effects Jiménez-Forteza  _et al._
(2017); Keitel  _et al._ (2017); Hinder  _et al._ (2014); Iozzo  _et al._
(2021). Since the fit results presented below in Sec. IV.3 attempt to recover
the final mass and spin from radiative quantities (_i.e._ , from the strain
mode $h_{22}$) and not from their quasilocal definitions, the errors on these
estimates will be better described by the errors on the parameters computed
from Eq. (8) rather than the errors from Eqs. (6) and (7). Thus, we consider
two type of basic errors estimates666Since most of the simulations used in
this work are shown to be in the convergent regime The SXS Collaboration
(2016), the main source of errors is either the resolution or the
extrapolation to null infinity of the NR datasets. On the other hand, other
sources of errors such as the conversion from $\psi_{4}$ to $h$ Hinder  _et
al._ (2011) have not been considered here.: the local error
$\delta\epsilon_{\mathrm{l}}$ as in Finch and Moore (2021) and the radiation
error $\delta\epsilon_{\mathrm{r}}$ that are both defined from the following
estimates:

$\small\delta\epsilon_{\mathrm{l},\mathrm{r}}=\sqrt{\left(\frac{{}^{(N)\\!}M_{f}^{\mathrm{l},\mathrm{r}}-{}^{(N-1)\\!}M_{f}^{\mathrm{l},\mathrm{r}}}{M}\right)^{2}+({}^{(N)}a_{f}^{{\mathrm{l},\mathrm{r}}}-{}^{(N-1)}a_{f}^{\mathrm{l},\mathrm{r}})^{2}}\,.$
(9)

Here, $M=m_{1}+m_{2}=1$, and the superscripts (N) and (N-1) stand either for
consecutive resolution levels or for consecutive extrapolation orders777One
can also use the same formula to compute a mass and spin discrepancy
$\delta\epsilon_{\mathrm{l}}$ or $\delta\epsilon_{\mathrm{r}}$ between
simulation results from two different codes — rather than comparing different
resolution levels or extrapolation orders within a given code — as we will do
in Sec. IV.3. . Our final radiative error estimates are obtained by combining
the resolution and the extrapolation errors, while only resolution effects are
relevant for the local estimates (see Sec. III.2). We have restricted the sum
of Eq. (8) to $l=m=2$, in order to avoid the error contributions sourced by
the higher angular modes. Notice that this is only valid to get the error on
the fit estimates from the $(22)$ mode but not accurate enough to compute the
final mass and final spin. The same holds for the initial ADM mass and angular
momentum in Eq. (8): their exact value and time at which they are computed
become irrelevant for the error estimation, since these terms are suppressed
when computing differences between resolutions or extrapolation orders as per
Eq. (9).

Figure 3: Distributions of error values obtained for the NR local error
(dashed gray), the NR radiation error (dashed-dotted black) and its
extrapolation (solid green) and resolution (dotted blue) contributions, for
the SXS waveforms considered here and when multiple resolutions were
available. We observe that the extrapolation errors are typically larger than
the resolution ones. The median values obtained are $\left\\{1.7\times
10^{-4},7.6\times 10^{-4},8.1\times 10^{-4}\right\\}$ for resolution,
extrapolation and combined radiation errors respectively. The distribution of
(resolution-based) local errors is additionally shown (dashed gray line) and
they can be noticed to be substantially smaller than the radiative errors,
with a median value $\widetilde{\delta\epsilon_{\mathrm{l}}}=4.6\times
10^{-5}$.

More specifically, the local errors have been obtained from the differences on
the mass and spin between the highest ($N$) and the second-highest ($N-1$)
resolution datum per NR case, following Eq. (9). This type of error results
from the discreteness of the grid, thus, its value depends on the sampling of
the numerical domain. In the case of the local errors, the discreteness
affects the computation of the integrals on the AH in Eq. (7). Since the
resolution is usually finer at the black hole ‘near-horizon’ length scale
$\mathcal{O}(1M)$, the local errors are in general smaller than the radiative
ones. Moreover, the radiative errors account for all types of inaccuracies
that have been propagated to the strain $h_{22}(t)$. Here, we have estimated
them from i) the resolution errors $\delta\epsilon_{\mathrm{r},\mathrm{res}}$
and ii) the extrapolation errors $\delta\epsilon_{\mathrm{r},\mathrm{extr}}$.
Like for the local errors, the resolution errors are computed from the
difference between the highest ($N$) and second-highest ($N-1$) resolutions,
but now estimated on the strain888BH merger simulations are solved by
splitting the whole space-time in a set of subdomains that range from the
black hole scale to the waves scale, where the black hole scale is usually the
finest and the wave scale the coarsest. $h_{22}(t)$ . The extrapolation errors
arise from the extrapolation of the strain to null infinity. Such
extrapolation is performed by fitting with second- to fourth-order polynomials
the phase and the amplitude of the strain multiplied by radius $r\,h_{22}(t)$
extracted on a set of several finite $\mathcal{O}(100M)$ distances from the
black hole local domain The SXS Collaboration (2016); Keitel  _et al._ (2017);
Jiménez-Forteza  _et al._ (2017); Iozzo  _et al._ (2021). Here, we estimated
the associated error also from Eq. (9) by taking the differences between the
successive extrapolation orders $N=2$ and $N=3$ on the waveform $h_{22}(t)$ at
the highest resolution999A lower-order polynomial typically performs better at
extrapolating the ringdown regime (see Sec. 2.4.1 of The SXS Collaboration
(2016)), hence we did not consider the fourth order, and we always take the
second-order extrapolation level as the default $h_{22}(t)$. The final
radiative error is estimated as,

$\delta\epsilon_{\mathrm{r}}=\sqrt{(\delta\epsilon_{\mathrm{r},\mathrm{res}})^{2}+\delta\epsilon_{\mathrm{r},\mathrm{extr}}^{2}}\,.$
(10)

In Fig. 3 we show the distribution of the errors $\delta\epsilon_{\mathrm{l}}$
and $\delta\epsilon_{\mathrm{r}}$ as well as the separate contributions
$\delta\epsilon_{\mathrm{r},\mathrm{res}}$,
$\delta\epsilon_{\mathrm{r},\mathrm{extr}}$ for the SXS waveforms considered.
To compute $\delta\epsilon_{\mathrm{l}}$,
$\delta\epsilon_{\mathrm{r},\mathrm{res}}$ and $\delta\epsilon_{\mathrm{r}}$,
we have only used the 343 simulations that are provided with multiple
resolution data. Notice that the extrapolation error (solid green) becomes the
major contribution to the uncertainty of the SXS dataset used in this work,
being typically larger than the resolution one. This can also be seen from the
median values of both error estimates: as expected, this value is larger for
the radiative error, with $\widetilde{\delta\epsilon_{\mathrm{r}}}=8.1\times
10^{-4}$, than for the local error, with
$\tilde{\delta\epsilon_{\mathrm{l}}}=4.6\times 10^{-5}$.

##  IV Setup and fit results

In this section we show the results obtained from the fits of the $(22)$ mode
of the NR waveforms by a range of ringdown models. Each model corresponds to
the $(22)$ mode of Eq. (1) where we have simply set $t_{r}=0$, and restricted
to a total number $N_{\mathrm{max}}$ of QNM overtones with
$N_{\mathrm{max}}\in\left\\{0,\dots,12\right\\}$ — and occasionally up to
$N_{\mathrm{max}}=16$. In particular, we consider two main scenarios: i) the
final mass and final spin are fixed and set equal to the known NR values; or
ii) we seek for the mass and the spin that minimise the fit mismatch
$\mathcal{M}$ (see below).

In these two scenarios, we have neglected the counter-rotating modes of Eq.
(1) since they are expected to have negligible amplitudes compared to the
corotating ones Jiménez Forteza  _et al._ (2020); Finch and Moore (2021) and
to have a negligible impact on the recovery of the final mass and spin Dhani
and Sathyaprakash (2021). Similarly, mode-mixing effects are as well discarded
due to their small impact on the $(22)$ mode Finch and Moore (2021); Cook
(2020); Dhani and Sathyaprakash (2021).

###  IV.1 The fitting algorithm

In scenario i), we fit for the $2(N_{\mathrm{max}}+1)$ parameters
$\vec{\lambda}=\\{A_{22n},\varphi_{22n}\\}$ with
$n\in\left\\{0,\dots,N_{\mathrm{max}}\right\\}$, for known final mass and
spin; while in scenario ii) we fit as well for the mass $M_{f}$ and the spin
$a_{f}$, while in scenario ii) we fit for the same set of parameters
$\vec{\lambda}$ over a range of $(M_{f},a_{f})$ values and then optimize the
results over this mass and spin range, thus accounting for
$2(N_{\mathrm{max}}+2)$ parameters in total.

Notice that once the values of the mass and spin are fixed, the RD ansatz (1)
is linear in the complex amplitudes $\mathcal{A}_{lmn}$. Therefore, one may
use a linear least-squares algorithm to obtain the fit results Giesler  _et
al._ (2019); Bhagwat  _et al._ (2020); Cook (2020). That is, for a given value
of the $(M_{f},a_{f})$ pair, the complex amplitudes $\mathcal{A}_{lmn}$ are
obtained by minimising the $\chi^{2}$,

$\chi^{2}=\sum_{k}\left|\bar{h}_{22}\\!\left(\vec{\lambda}\right)(t_{k})-h_{22}(t_{k})\right|^{2},$
(11)

where the subscript $k$ labels the values of the time axis of the NR waveform,
$t_{k}\in\left[t_{0},t_{f}\right]$ for a certain fit starting time $t_{0}$ and
with $t_{f}=90M$; and $\bar{h}_{22}(\vec{\lambda})$ denotes the model
$(2,2)$-mode strain for the set of parameters $\vec{\lambda}$. By default in
the following, the starting time is set to $t_{0}=0$, which corresponds to the
peak of the $(22)$ mode of the strain $h_{22}(t)$. We however let this value
vary in Secs. IV.2 and IV.3.4 as specified there. The best-fit parameters per
RD model $\vec{\lambda}^{\mathrm{bf}}(N_{\mathrm{max}})$ are chosen as the
ones that minimise Eq. (11).

The above fully describes the fitting procedure in scenario i). In scenario
ii), the same process is iterated over a range of $(M_{f},a_{f})$ values to
find the optimal one. To this aim, we build a two-dimensional adaptive grid on
the final mass $M_{f}$ and the final spin $a_{f}$, with a grid minimum step
set to $3.2\cdot 10^{-6}$ in both variables. Every point of the grid is then
treated as a linear least-squares minimization problem on the parameters
$\vec{\lambda}=\\{A_{22n},\varphi_{22n}\\}$ as above Finch and Moore (2021);
Mourier and Jimenez Forteza ; Jimenez Forteza and Mourier . Closely related to
the $\chi^{2}$ and recurrently used in GW astronomy, we compute the mismatch
for each best-fit RD model, namely101010In this framework, one can easily show
that both $\mathcal{M}$ and the $\chi^{2}$ provide the same qualitative
behavior. In particular, for a model closely fitting the NR waveform,
${\chi^{2}\simeq 2\,\mathcal{M}\,(\sum_{k}|h_{22}(t_{k})|^{2})}$. Therefore, a
minimum on $\chi^{2}$ directly translates to a minimum in $\mathcal{M}$ and
_vice versa_.,

$\mathcal{M}=1-\frac{\left\langle
h_{22}\,\middle|\,\bar{h}_{22}\left(\vec{\lambda}^{\mathrm{bf}}\right)\right\rangle}{\sqrt{\Big{\langle}h_{22}\,\Big{|}\,h_{22}\Big{\rangle}\left\langle\bar{h}_{22}\\!\left(\vec{\lambda}^{\mathrm{bf}}\right)\,\middle|\,\bar{h}_{22}\\!\left(\vec{\lambda}^{\mathrm{bf}}\right)\right\rangle}}\,,$
(12)

where

$\langle f|g\rangle=\int_{t_{0}}^{t_{f}}f(t)\,g(t)^{*}\,dt\,.$ (13)

Finally, the best-fit mass and spin values
$(M_{f}=M_{f}^{\mathrm{fit}},a_{f}=a_{f}^{\mathrm{fit}})$ for the given
waveform and the given number of overtones $N_{\mathrm{max}}$ of the RD model
are selected as the grid point where $\mathcal{M}$ from Eq. (12) is minimal.

The associated value of the minimum $\mathcal{M}$ for each RD
$N_{\mathrm{max}}$ model is sufficient to assess the fit accuracy but
insufficient to determine whether the fitting parameters are physically
reliable. A decreasing value of the mismatch $\mathcal{M}$ between different
models is particularly sensitive to overfitting, especially if it is applied
to nested models such as the RD models we have considered in this work (the RD
model with $N_{\mathrm{max}}-1$ overtones corresponds to the subclass of the
RD model with $N_{\mathrm{max}}$ overtones with
$\mathcal{A}_{N_{\mathrm{max}}}$ set to $0$). To overcome this issue we use
the mass and spin bias $\epsilon$ defined in Eq. (4) of Giesler  _et al._
(2019),

$\epsilon=\sqrt{\left(\frac{\delta M_{f}}{M}\right)^{2}+\delta a_{f}^{2}}\;,$
(14)

where $\delta M_{f}=M_{f}^{\mathrm{fit}}-M_{f}^{\mathrm{true}}$ and $\delta
a_{f}=a_{f}^{\mathrm{fit}}-a_{f}^{\mathrm{true}}$. Thus, $\epsilon$ measures
the combined deviation of the final mass $M_{f}$ and the final spin $a_{f}$
with respect to the true parameters $M_{f}^{\mathrm{true}}$ and
$a_{f}^{\mathrm{true}}$ of the NR simulation, that are estimated from the mass
and spin quasilocal definitions The SXS Collaboration (2016), _i.e._ ,
following Eqs. (6) and (7).

###  IV.2 Fits with the mass and spin fixed to their true values

First, we show the results obtained by fitting the RD models (1) to the NR
waveform SXS:0305 following the same methodology described in Bhagwat  _et
al._ (2020); Jiménez Forteza  _et al._ (2020); Mourier  _et al._ (2021). With
$a_{f}^{\mathrm{true}}=0.692$ and $M_{f}^{\mathrm{true}}=0.952$, this waveform
is consistent with the first gravitational wave event GW150914, and it has
been recurrently used in several RD studies Giesler  _et al._ (2019); Bhagwat
_et al._ (2020). In Fig. 4 we show the mismatch curves for a set of models
with a number of overtones $N_{\mathrm{max}}$ spanning
$\left\\{0,\dots,10\right\\}$. In the RD models (1) used for these fits, the
final mass and spin are fixed (scenario i), which implies that the whole set
of QNM frequencies and damping times are fixed to their GR values. We analyse
the fit results as a function of the fitting starting time $t_{0}/M$. Notice
that the mismatch diminishes as the number of overtones $N_{\mathrm{max}}$
increases for all RD models and for any fit starting time $t_{0}$. For each RD
model, we find a local minimum in $\mathcal{M}$ as $t_{0}$ varies, after an
initial strong decrease and followed by a plateau of nearly-constant
$\mathcal{M}$. This minimum111111This first local minimum is the global
minimum in $\mathcal{M}$ for $N_{\mathrm{max}}\leq 6$. For larger
$N_{\mathrm{max}}$ values, the global minimum of $\mathcal{M}$ is different
and occurs at a later fit starting time $t_{0}\simeq 20M$, but it is still
only marginally smaller that the first local minimum.  occurs at increasingly
early starting times $t_{0}$ as the number of overtones $N_{\mathrm{max}}$
increases. In particular, it occurs at $t_{0}\simeq 0$ for the
$N_{\mathrm{max}}=7$ model as it has been observed in Giesler  _et al._
(2019); Bhagwat  _et al._ (2020); Finch and Moore (2021). For the new
$N_{\mathrm{max}}=8,9,10$ models, the local minimum in mismatch occurs at some
$t_{0}<0$. The same trend continues for all the subsequent models tested in
this work (that is, with ${10<N_{\mathrm{max}}\leq 16}$), which have not been
included here for the sake of the plot clarity.

At large $N_{\mathrm{max}}$, the decrease in the mismatch value (hence also of
the $\chi^{2}$) with increasing $N_{\mathrm{max}}$ may be mostly due to
overfitting of the data. Namely, increasing $N_{\mathrm{max}}$ increases the
number of free parameters in the model accordingly, which induces a decrease
in $\mathcal{M}$ and may become the main source of the observed decrease as
$N_{\mathrm{max}}$ gets large. In the next section we discuss how $\epsilon$
may be used as an approximate indicator to detect the overfitting in our RD
models.

Figure 4: We show the mismatch ${\mathcal{M}}$ at best-fit complex amplitudes
$\mathcal{A}_{22n}$ from Eq. (12) for a range of RD models with
$N_{\mathrm{max}}\in\left\\{0,\dots,10\right\\}$ as a function of the fit
starting time $t_{0}/M$. Notice that as the number of tones of the model
increases, the mismatch decreases progressively. In addition, for each
$N_{\mathrm{max}}$, one finds a local minimum on ${\mathcal{M}}$ as a function
of $t_{0}/M$, which occurs at progressively smaller fit starting times
$t_{0}/M$ as $N_{\mathrm{max}}$ increases.

###  IV.3 Varying the final mass and spin

####  IV.3.1 Results comparison for SXS/RIT waveform pairs

We now set $t_{0}=0$ and we study the behavior of $\epsilon$ for two pairs of
NR simulations, (SXS:0305, RIT:0062) and (SXS:0259, RIT:0118). Each pair
corresponds to merger simulations from two different catalogs (the SXS catalog
The SXS Collaboration (2016) and the RIT catalog Campanelli  _et al._ (2016))
with consistent values of every physical parameter. The true values of the
final mass and spin are $M_{f}=\left\\{0.952\pm 1.2\times 10^{-5},0.966\pm
3\times 10^{-5}\right\\}$ and $a_{f}=\left\\{0.692\pm 1.2\times
10^{-4},0.581\pm 2\times 10^{-5}\right\\}$ for the first and second pair of
simulations respectively, which correspond to merging binary BHs with mass
ratio $q=\left\\{1.22,2.5\right\\}$ and effective spin
$\chi_{\mathrm{eff}}=\left\\{-0.0165,0\right\\}$ respectively. The values of
uncertainty quoted on $(M_{f},a_{f})$ are computed from the differences on the
local final mass and spin as
$\left|(M_{f}^{\mathrm{l}})_{\mathrm{SXS}}-(M_{f}^{\mathrm{l}})_{\mathrm{RIT}}\right|$
and
$\left|(a_{f}^{\mathrm{l}})_{\mathrm{SXS}}-(a_{f}^{\mathrm{l}})_{\mathrm{RIT}}\right|$,
that translate to a local discrepancy
$\delta\epsilon^{\mathrm{l}}_{\mathrm{SXS}-\mathrm{RIT}}=\left\\{1.2\times
10^{-4},3.6\times 10^{-5}\right\\}$ respectively for the two simulation pairs.

In Fig. 5 we show on a log–log scale, the values obtained for ${\epsilon}$ and
${\mathcal{M}}$ for the two pairs of simulations and for a set of RD models
with $N_{\mathrm{max}}\in\left\\{0,\dots,16\right\\}$. On the top panel,
corresponding to the first simulation pair (SXS:0305, RIT:0062), we observe
that both ${\epsilon}$ and $\mathcal{M}$ decrease as $N_{\mathrm{max}}$
increases up to $N_{\mathrm{max}}=7$, where $\epsilon$ reaches a minimum at
$\epsilon\sim 3\cdot 10^{-4}$ as observed in Giesler  _et al._ (2019); Finch
and Moore (2021). This has been considered as one possible empirical evidence
that i) one needs $N_{\mathrm{max}}=7$ overtones to describe the post-peak
data Giesler  _et al._ (2019); Finch and Moore (2021) and ii) post-peak
nonlinearities are subdominant even at $t_{0}=0$. At $N_{\mathrm{max}}>7$ and
for both waveforms, the mismatch keeps decreasing at a reduced rate, while
$\epsilon$ increases. This trend reaches a saturation point at $\epsilon\sim
5\cdot 10^{-3}$ for $N_{\mathrm{max}}\sim 16$. Thus, in this case, we do not
improve any further the accuracy on the estimate of the mass and spin beyond
$N_{\mathrm{max}}=7$, for both SXS:0305 and RIT:0062. This could be taken as
the threshold point beyond which overfitting could be significantly affecting
the fits. However, this behavior is rather variable when studying other NR
cases (see bottom panel and its discussion below, and Sec. IV.3.2), where one
finds that the minimum $\epsilon$ point is case-dependent. On the other hand,
in Sec. IV.3.3 we show that large instabilities in the best-fit amplitudes
could be affecting the tones at $n\geq 2$, thus, any claim about the onset
point of overfitting shall be taken with caution.

The shaded gray area delimits the domain for which the mismatch is lower than
the mismatch between the two waveforms
$\mathcal{M}_{\mathrm{SXS}-\mathrm{RIT}}$, and where
$\epsilon\leq\delta\epsilon_{\mathrm{r},\mathrm{SXS}-\mathrm{RIT}}$.
Conversely, the smaller shaded orange area near the lower-left corner of each
plot stands for the radiation error on the SXS data alone (since only one
resolution level per case is provided for the RIT catalog); see Sec. III.2 for
further details on the computation of
$\delta\epsilon_{\mathrm{r},\mathrm{l}}$. In particular, the upper bound of
the SXS error on the mismatch axis — the mismatch horizontal orange line — is
estimated as the maximum mismatch that results from comparing both the two
highest resolution and the two best extrapolation levels ($N=2$ and $3$),
namely, ${\rm
max}\left(\mathcal{M}_{\mathrm{res}},\mathcal{M}_{\mathrm{extr}}\right)$.

In the bottom panel we show the $\mathcal{M}-\epsilon$ results for the higher
mass ratio pair (SXS:0259, RIT:0118). Notice that the trend on $\epsilon$
changes substantially compared to the previous case. As expected, the mismatch
always decreases, but flattens out at $N_{\mathrm{max}}=6$, especially for
SXS:0259. Here however, for SXS:0259, the value of $\epsilon$ decreases to
eventually reach its minimum only at $N_{\mathrm{max}}=13$ with $\epsilon\sim
3\cdot 10^{-4}$ — a similar minimal value as the one reached at a smaller
$N_{\mathrm{max}}$ in the case of SXS:0305. For $N_{\mathrm{max}}>13$,
$\epsilon$ increases again although it is not yet saturated at
$N_{\mathrm{max}}=16$. In the case of RIT:0118, $\epsilon$ decreases and hits
its minimum at $N_{\mathrm{max}}=5$ with a larger value $\epsilon\sim 3\cdot
10^{-3}$ to thereafter grow, saturate at $\epsilon\sim 10^{-2}$ for
$N_{\mathrm{max}}\sim 9$ and decrease again. These discrepancies between the
behaviors observed for both waveforms on the lower panel arise in the
$\mathcal{M}\leq\mathcal{M}_{\mathrm{SXS}-\mathrm{RIT}}$,
$\epsilon\leq\delta\epsilon_{\mathrm{r},\mathrm{SXS}-\mathrm{RIT}}$ domain,
thus they could be affected the by NR errors — or other systematics — of each
code.

It is moreover noteworthy that the minimal $\epsilon$ values — or the turning
points on the $\mathcal{M}-\epsilon$ plane — occur close the boundary
delimited by the SXS radiative error $\delta\epsilon_{\mathrm{r}}$ (orange
areas) for both SXS waveforms analysed here. This could indicate that the
change of trend for $N_{\mathrm{max}}\geq
N_{\mathrm{max}}(\epsilon_{\mathrm{min}})$ for the SXS models may be dominated
by the NR uncertainties.

(a)

(b)

Figure 5: We show the best-fit ${\epsilon-\mathcal{M}}$ plane for a range of
models with $N_{\mathrm{max}}\in\left\\{0,\dots,16\right\\}$ and for two NR
waveforms from the SXS catalog (green) and two waveforms from the RIT catalog
(orange), with initial parameters $(q=1.22,\chi_{\mathrm{eff}}=-0.0165)$ (top
panel) and $(q=2.5,\chi_{\mathrm{eff}}=0)$ (bottom panel). The gray lines and
shaded areas on each panel delimit the mismatch and $\epsilon$ values that are
respectively smaller than the mismatch
$\mathcal{M}_{\mathrm{SXS}-\mathrm{RIT}}$ and radiative discrepancy on the
mass and spin $\delta\epsilon_{\mathrm{r},\mathrm{SXS}-\mathrm{RIT}}$ between
the two waveforms considered in the panel. The smaller orange shaded areas at
the lower-left corners show the same in terms of the maximum
resolution/extrapolation mismatch ${\rm
max}\left(\mathcal{M}_{\mathrm{res}},\mathcal{M}_{\mathrm{extr}}\right)$, and
of the radiative error $\delta\epsilon_{\mathrm{r}}$, of each SXS waveform.
The two green and orange points furthest on the upper-right corner of each
panel correspond to the $N_{\mathrm{max}}=0$ model. In the top panel, notice
that as $N_{\mathrm{max}}$ increases the ${\epsilon-\mathcal{M}}$ points are
progressively shifted to the left bottom corner until $N_{\mathrm{max}}=7$,
where the minimum $\epsilon$ is achieved for both NR simulations. Beyond
$N_{\mathrm{max}}=7$, $\mathcal{M}$ keeps decreasing —albeit more marginally—
while $\epsilon$ increases. Conversely, for the case shown in the bottom
panel, we observe that the trend and the values at which
$\epsilon_{\mathrm{min}}$ is achieved are significantly different between the
two simulations. In particular, $\epsilon$ now reaches its minimum at
$N_{\mathrm{max}}=13$ for the SXS simulation while the much larger minimal
value of $\epsilon$ for the RIT waveform is reached at $N_{\mathrm{max}}=5$.
The difference on the true parameters between the waveforms from both codes is
$\delta\epsilon_{\mathrm{l},\mathrm{SXS}-\mathrm{RIT}}\lesssim 10^{-4}$, thus
much smaller than the radiative errors
$\delta\epsilon_{\mathrm{r},\mathrm{SXS}-\mathrm{RIT}}$ (gray vertical lines).

####  IV.3.2 Mass and spin recovery biases for the set of non-precessing SXS
simulations

Once studied individually the above two NR cases, we extend this analysis to
the set121212We excluded $\fpeval{620-610}$ out of the 620 such waveforms in
the catalog, which did not appear reliable enough for our analysis. We list
these cases in Table 1 along with the reasons of their exclusion.  of non-
precessing SXS binary-black hole waveforms The SXS Collaboration (2016). In
particular, we want to explore whether the consideration of
$N_{\mathrm{max}}>7$ models allows us to find which number of overtones is
statistically preferred over this set of NR waveforms. In Fig. 6 we show the
distributions of the values obtained for $\epsilon$ over the 610 waveforms
considered, for each of five RD models with
$N_{\mathrm{max}}\in\\{1,3,7,8,9\\}$. Consistently with the particular cases
shown on Fig. 5, we find that among these models the distributions for
$N_{\mathrm{max}}=1$ and $N_{\mathrm{max}}=3$ provide the largest values of
$\epsilon$, with the median values $\tilde{\epsilon}\simeq 10^{-1}$ and
$\tilde{\epsilon}\simeq 1.3\times 10^{-2}$ respectively, while the
distributions on $\epsilon$ for the $N_{\mathrm{max}}=7,8,9$ models are
shifted to significantly lower values. For instance, we have obtained a median
value $\tilde{\epsilon}=2.2\times 10^{-3}$ for $N_{\mathrm{max}}=7$,
consistent with Giesler  _et al._ (2019); Finch and Moore (2021). We do not
observe significant differences between the $N_{\mathrm{max}}=7,8,9$ models,
where all three distributions overlap within the $10$–$90$ percentiles. We
also show on this figure the distribution of the NR radiative
($\delta\epsilon_{\mathrm{r}}$) and local ($\delta\epsilon_{\mathrm{l}}$)
errors. For the radiative errors, we have taken into account the resolution
and extrapolation errors131313274 out of the 610 SXS waveforms discussed here
are only available at a single resolution. Thus, these cases have not been
accounted for in our NR error estimates.. Notice that the distribution on the
NR local errors does slightly overlap with the $N_{\mathrm{max}}=7,8,9$
distributions. On the other hand, the radiative errors broadly overlap with
the $N_{\mathrm{max}}=7,8,9$ distributions of $\epsilon$. As described in Sec.
III.2, we have obtained $\widetilde{\delta\epsilon_{\mathrm{r}}}\simeq
8.1\times 10^{-4}$ for the radiative error, thus a slightly smaller but
comparable value to $\tilde{\epsilon}(N_{\mathrm{max}}=7)$. In contrast, we
have obtained a much smaller median value,
$\widetilde{\delta\epsilon_{\mathrm{l}}}\simeq 4.6\times 10^{-5}$, for the
local error.

Figure 6: We show the ${\epsilon}$ distribution for five models with
$N_{\mathrm{max}}\in\\{1,3,7,8,9\\}$ computed from the set of non-precessing
NR waveforms from the SXS catalog, together with the distribution of NR error
estimates for the waveforms for which multiple resolutions are available. The
models with $N_{\mathrm{max}}=1,3$ show significantly larger values for the
$\epsilon$ distributions, with $\tilde{\epsilon}=10^{-1},1.3\cdot 10^{-2}$
respectively. On the other hand, we have obtained $\tilde{\epsilon}\sim
2\times 10^{-3}$ for the $N_{\mathrm{max}}=7,8,9$ models, thus not showing
significant differences among them.

In Fig 7 we show the $\epsilon$ estimates (shaded colored curves) for the 610
cases we have analysed as a function of the number of overtones
$N_{\mathrm{max}}$ of each RD model. The corresponding median values
$\tilde{\epsilon}$ (diamonds) and the $10$–$90$ percentiles of the
distribution (bars) are also shown for each $N_{\mathrm{max}}$. The shaded
curves are split in terms of the final spin value as follows: $a_{f}>0.7$ in
light gray, $0\leq a_{f}<0.7$ in light blue and the 18 cases with $a_{f}<0$ in
light red141414Since we restrict our models to the corotating modes, we did
use in our fits the prior knowledge of the sign of the true parameter
$a_{f}^{\mathrm{true}}$ to select those corotating modes adapted to this sign.
For $a_{f}^{\mathrm{true}}<0$, these modes may then be directly obtained from
the positive-spin _counter-rotating_ solutions _via_ the symmetry relation
(4). Note that in this case $N_{\mathrm{max}}$ is still to be understood as
the total number of distinct overtones in the model, so that the
$N_{\mathrm{max}}=10$ model for instance will be comprised of the
$n=0,1,\dots,7,8,10,11$ $a_{f}>0$ counter-rotating branches in our notations,
since there is no distinct $n=9$ such branch (see Sec. II.3). For the
consistency of the mode selection, we restrict the allowed spin range on the
$(M_{f},a_{f})$ grid to spins of the same sign as $a_{f}^{\mathrm{true}}$ —
either $a_{f}\in[0,0.99]$ or $a_{f}\in[-0.99,0]$. In a few cases where
$a_{f}^{\mathrm{true}}$ is close to $0$ (with either sign) and with small
$N_{\mathrm{max}}$ values (when $\epsilon$ is large), this may lead to an
underestimated $\epsilon$ when the minimum-$\mathcal{M}$ solution lies at the
$a_{f}=0$ boundary of the allowed spin range. . First, notice that all
$\epsilon$ curves with $a_{f}\geq 0$ (light gray and light blue) show a
similar trend as we vary the number of overtones $N_{\mathrm{max}}$. This
allow us to disregard possible artifacts originating from the mode mixing,
i.e., those possible artifacts that result from decomposing the NR strain in
terms of spherical rather than spheroidal harmonics, and that would be
expected to mostly arise at high spins Cook (2020); Finch and Moore (2021). On
the other hand, we do observe a higher concentration of the negative spins
(light red curves) at high $\epsilon$. For comparison, at $N_{\mathrm{max}}=8$
as an example, the median value of $\epsilon$ over the negative-spin cases
(light red curves) is $\tilde{\epsilon}_{a_{f}<0}\,\sim 2.0\cdot 10^{-2}$,
substantially larger than the medians of the intermediate positive-spin cases
($0<a_{f}<0.7$, corresponding to the light blue curves) and of the high-spin
cases ($a_{f}>0.7$, corresponding to the light gray curves),
$\tilde{\epsilon}_{0<a_{f}<0.7}\,\sim 2.4\cdot 10^{-3}$ and
$\tilde{\epsilon}_{a_{f}>0.7}\,\sim 1.2\cdot 10^{-3}$ respectively. We have
found that the radiative errors $\delta\epsilon_{\mathrm{r}}$ obtained for
$a_{f}<0$ are the largest among the NR setup, hence the NR uncertainties could
explain the high $\epsilon$ values obtained in those cases. The minimum-
mismatch value also deteriorates for $a_{f}<0$ with respect to $a_{f}\geq 0$
cases, for $N_{\mathrm{max}}\geq 5$. Beyond this $N_{\mathrm{max}}$, the
median of the best-fit mismatch values over all negative-spin cases
approximately plateaus at $\widetilde{\mathcal{M}}_{a_{f}<0}\,\sim 4\cdot
10^{-6}$, while for intermediate positive spins ($0<a_{f}<0.7$) and high spins
($a_{f}>0.7$) the corresponding medians only saturate at larger
$N_{\mathrm{max}}$ and reach lower values,
$\widetilde{\mathcal{M}}_{0<a_{f}<0.7}\,\sim 2\cdot 10^{-7}$ and
$\widetilde{\mathcal{M}}_{a_{f}>0.7}\,\sim 7\cdot 10^{-8}$ respectively for
these two positive-spin classes. The source of these discrepancies may lie on
the numerical setup of the NR simulations and its full exploration may require
some further investigation. On the other hand, we find that $\tilde{\epsilon}$
decreases before flattening out at ${N_{\mathrm{max}}\sim 5,6}$. Beyond this
point, the values of $\tilde{\epsilon}$ remain approximately stable at
$\tilde{\epsilon}\sim 2\times 10^{-3}$. The 18 cases with $a_{f}<0$ do not
increase substantially the value of $\tilde{\epsilon}$ since they still
represent a small fraction of the NR simulations studied. However, this may
need to be reviewed if more NR simulations with negative final spin are added
to the catalog. The orange shaded area in the lower half of the plot accounts
for the $10$–$90$ percentiles obtained from the radiative error distribution
while the dashed black line stands for its median value. Note that the $90$th
percentile (upper bound) of this error lies above the median values
$\tilde{\epsilon}$ of the $\epsilon$ distribution for $N_{\mathrm{max}}\geq 6$
(with the exception of $\tilde{\epsilon}(N_{\mathrm{max}}=12)$ lying slightly
above this line), which suggests that the waveform inaccuracies could be
affecting the estimates of $\tilde{\epsilon}$ at high $N_{\mathrm{max}}$.
Furthermore, it is noticeable that the median $\epsilon$ values for the
$N_{\mathrm{max}}=3-4$ models lie within the $10$–$90$ percentile bands of the
$N_{\mathrm{max}}\geq 5-6$ models. For comparison, the black dotted line shows
the median of the local error estimate
$\widetilde{\delta{\epsilon}_{\mathrm{l}}}$. This value is far below the
estimates obtained for $\epsilon$, but we recall here that the radiative error
$\widetilde{\delta{\epsilon}_{\mathrm{r}}}$ provides a conceptually more
appropriate measure of the error since it is computed directly from the (22)
mode of the strain, _i.e._ , from the data used to compute our fits.

Figure 7: We show the median value ${\tilde{\epsilon}}$ (diamonds) of the
distribution of $\epsilon$ values over the SXS waveforms studied, as a
function of the number of overtones $N_{\mathrm{max}}$ included in the model.
For each $N_{\mathrm{max}}\geq 3$, a second value is shown (circles) which has
been computed from the median $\tilde{\epsilon}_{\mathrm{min}}$ of all the
studied NR cases for which $\epsilon$ is minimum at $N_{\mathrm{max}}$ among
the $N_{\mathrm{max}}=0,\dots,12$ models shown here (there are no such cases
for $N_{\mathrm{max}}<3$, and a single such case for $N_{\mathrm{max}}=3$).
The error bars represent the $10$–$90$ percentiles of the $\epsilon$
distribution obtained for each of the respective $N_{\mathrm{max}}$ models.
The values obtained for $\epsilon$ as a function of $N_{\mathrm{max}}$ for
each of the 610 cases considered in this work are also shown individually as
shaded colored curves. The different colors of these curves represent the
cases belonging to different classes of final spin values (see details in the
main text, Sec. IV.3.2). We observe that the median values reach an
approximate plateau regime at ${N_{\mathrm{max}}\sim 5-6}$ and
${N_{\mathrm{max}}\sim 7}$ for $\tilde{\epsilon}$ and
$\tilde{\epsilon}_{\mathrm{min}}$ respectively. The shaded orange band shows
the $10$–$90$ percentiles of the radiation error distribution, with the dashed
black horizontal line denoting the median value
$\widetilde{\delta\epsilon_{\mathrm{r}}}$. The dotted black horizontal line
accounts for the median value for the local error,
$\widetilde{\delta\epsilon_{\mathrm{l}}}$.

In addition, we present on the same figure a second estimate
$\tilde{\epsilon}_{\mathrm{min}}$ for each value of $N_{\mathrm{max}}\geq 3$
(circles), computed as the median of $\epsilon$ over all the cases for which
this value of $N_{\mathrm{max}}$ minimizes $\epsilon$ among the
$N_{\mathrm{max}}\in\\{0,\dots,12\\}$ models considered in this analysis. We
have not found any case among the SXS waveforms considered for which
$\epsilon$ reaches its minimum at an $N_{\mathrm{max}}<3$. The value of
$\tilde{\epsilon}_{\mathrm{min}}$ for each $N_{\mathrm{max}}$ is smaller than
$\tilde{\epsilon}$, since all the cases for which $\epsilon$ is not at its
minimum at $N_{\mathrm{max}}$ have been excluded from the distribution in
computing $\tilde{\epsilon}_{\mathrm{min}}$. Similarly to $\tilde{\epsilon}$,
the values of $\tilde{\epsilon}_{\mathrm{min}}$ decrease with
$N_{\mathrm{max}}$ before approximately stabilizing for
${N_{\mathrm{max}}\gtrsim 7}$ at $\tilde{\epsilon}_{\mathrm{min}}\simeq
6\times 10^{-4}$, which is smaller than the value of $\tilde{\epsilon}$ for
those models. Nevertheless, these values of $\tilde{\epsilon}_{\mathrm{min}}$
lie within the ${\epsilon}$ $10$–$90$ bands and well within the radiative
error ${\delta\epsilon_{\mathrm{r}}}$ distribution.

Figure 8: Fraction of the NR cases studied at which $\epsilon$ reaches its
minimum value (among the $N_{\mathrm{max}}\leq 12$ models) at a number of
overtones $N_{\mathrm{max}}$, as a function of $N_{\mathrm{max}}$. The color
function provides the logarithm of the median values over the cases selected
in this way $\log_{10}(\tilde{\epsilon}_{\mathrm{min}})$ for each
$N_{\mathrm{max}}$ as in Fig. 7 (displayed as circles there). We have not
found any model reaching its minimum $\epsilon$ at any $N_{\mathrm{max}}<3$.

Finally, in Fig. 8 we show in histogram form, the fraction of cases for which
$\epsilon$ reaches its minimum at $N_{\mathrm{max}}$ — among the
$N_{\mathrm{max}}\in\\{0,\dots 12\\}$ models considered for this figure —, for
each $N_{\mathrm{max}}$. Here, the color map recalls in log scale the median
value of $\epsilon$ among the corresponding cases for each $N_{\mathrm{max}}$,
that is, the $\log_{10}$ of the same values $\tilde{\epsilon}_{\mathrm{min}}$
as indicated by circles on Fig. 7. The bulk of the distribution peaks at
$N_{\mathrm{max}}\sim 7$ although we observe that a large fraction (about
$21\%$) of the cases have a minimum $\epsilon$ at $N_{\mathrm{max}}=12$. These
include many cases which do not actually reach their minimum $\epsilon$ within
the range $N_{\mathrm{max}}\leq 12$ considered here and that would ideally
require to be fit with $N_{\mathrm{max}}>12$ overtones, such as SXS:0259
reaching its minimum $\epsilon$ at $N_{\mathrm{max}}=13$ as shown on the lower
panel of Fig. 5. On the other hand, the corresponding median value of
$\epsilon$ for these cases,
$\tilde{\epsilon}_{\mathrm{min}}(N_{\mathrm{max}}=12)$, is moderately higher
than the values obtained for cases that reach their minimum earlier,
$\tilde{\epsilon}_{\mathrm{min}}(N_{\mathrm{max}}=7,\dots,11)$, and is still
compatible with the NR error estimates given by $\delta\epsilon_{\mathrm{r}}$.
Therefore, and based on the trend observed for
$\tilde{\epsilon}_{\mathrm{min}}$ in Fig. 7 and on the examples of Fig. 5, we
do not expect these values to get significantly smaller at
$N_{\mathrm{max}}>12$.

####  IV.3.3 Stability of the fit amplitudes

We study the behavior of the recovered (best-fit) amplitudes of the first five
tones $A_{n}^{N_{\mathrm{max}}}$, $n=0,...,4$, as we increase the number of
overtones $N_{\mathrm{max}}$ of our models. We require as a criterion for a
stable recovery of a given tone $n$, that $A_{n}^{N_{\mathrm{max}}}$ remain
approximately constant as we modify the number of overtones
$N_{\mathrm{max}}$. In Fig. 9, we show in log scale the relative variation (in
percent) of the best-fit amplitudes of each tone $\delta
A_{n}(N_{\mathrm{max}})=\left|A_{n}^{N_{\mathrm{max}}}-A_{n}^{N_{\mathrm{max}}-1}\right|/A_{n}^{N_{\mathrm{max}}-1}$
between successive models as a function of $N_{\mathrm{max}}$. The shaded
colored curves stand for the estimates of $\delta A_{n}$ for each of the $610$
SXS simulations used in this work while the dots stand for the median value
$\delta\tilde{A}_{n}$ of $\delta A_{n}$ for each ringdown model with a number
$N_{\mathrm{max}}$ of overtones.

For the fundamental mode, we observe that $\delta\tilde{A}_{0}$ decreases
exponentially with $N_{\mathrm{max}}$ before stabilizing at
$N_{\mathrm{max}}\geq 6$. Remarkably, the flattening of the curve is very
similar to the one observed for $\tilde{\epsilon}$ on Fig. 7. Beyond
$N_{\mathrm{max}}=6$ the median relative variation of the amplitude $A_{0}$
remains nearly constant with $\delta\tilde{A}_{0}~{}\sim 0.2\%$ . This is
consistent with the results shown in Appendix B, where we have performed a
similar analysis (focusing on $\delta\tilde{A}_{0}$ and $\delta\tilde{A}_{1}$)
but varying the fit starting time $t_{0}/M$ at fixed $N_{\mathrm{max}}$.

On the other hand, while an initial exponentially decreasing trend with
$N_{\mathrm{max}}$ is also observed for all the overtone modes shown here, the
relative variation on the overtone amplitudes is larger than that of the
fundamental mode, and increases with $n$. For instance, notice that the $n=1$
amplitude typically varies by about $10\%$ at $N_{\mathrm{max}}\sim 5$ and
that this variation only achieves the $\sim 1\%$ level at
$N_{\mathrm{max}}\gtrsim 8$. The deviations become increasingly larger for
$N_{\mathrm{max}}=2,3,4$, where $\delta\tilde{A}_{2}\gtrsim 5\%$,
$\delta\tilde{A}_{3}\gtrsim 20\%$ and $\delta\tilde{A}_{4}\gtrsim 50\%$ for
all $N_{\mathrm{max}}$ values. This symptom of instability seems not to be
affecting the estimates of $\tilde{\epsilon}$ for $N_{\mathrm{max}}\geq 5$,
which surprisingly coincides with the flattening of the $\delta\tilde{A}_{0}$
curve.

Figure 9: In this figure we show the relative variation of the best-fit
amplitudes of each tone
$\delta{A}_{n}(N_{\mathrm{max}})=\left|{A}_{n}^{N_{\mathrm{max}}}-{A}_{n}^{N_{\mathrm{max}}-1}\right|/{A}_{n}^{N_{\mathrm{max}}-1}$
as a function of the number of overtones $N_{\mathrm{max}}$ included in the
model, for the $610$ SXS waveforms used in this work (shaded colored curves)
and for $n=0,\dots,4$. As we have done for $\tilde{\epsilon}$ in Fig. 7, for
each tone $n$, the dots represent the median values $\delta\tilde{A}_{n}$ per
$N_{\mathrm{max}}$ model. Each color corresponds to a given tone index $n$,
consistently for the median values and the individual curves. Since we
evaluate the differences of the best-fit amplitudes $A_{n}$ between the
consecutive $N_{\mathrm{max}}$–overtone and $(N_{\mathrm{max}}-1)$–overtone
models, the first point we can evaluate for each $\delta{A}_{n}$ curve
corresponds to the variation with respect to the $N_{\mathrm{max}}=n+1$ model,
$\delta{A}_{n}(N_{\mathrm{max}}=n+1)$.

####  IV.3.4 Analysis at different fit starting times $t_{0}$

So far, we have studied the NR waveforms using a single starting time
$t_{0}=0$ for all fits as in Finch and Moore (2021); Giesler  _et al._ (2019).
This specific time corresponds to the time of the peak of the $h_{22}(t)$
strain. However, since this particular time does not have any special physical
meaning151515One could reasonably choose $t_{0}$ as the time at which the
final common horizon is formed. However, this time is only well-defined
locally while it is causally disconnected from the events happening in the
radiation zone. See Mourier  _et al._ (2021) and references therein., we
extend here the analysis for a broader range of fit starting times, with
$t_{0}/M=\left\\{-5,0,5,10,15\right\\}$. Thus, in Fig. 10 we show the median
value $\tilde{\epsilon}$ as a function of the number of overtones
$N_{\mathrm{max}}$ and for the five values of $t_{0}/M$ selected above.

First, notice that all the curves studied here, show at first a progressively
decreasing value for $\tilde{\epsilon}$ as we incorporate more tones
$N_{\mathrm{max}}$ into the model. This decreasing trend reaches approximately
a minimum $\tilde{\epsilon}\sim(2-3)\cdot 10^{-3}$ (with similar minimum
values for all five $t_{0}/M$ choices shown) at
$N_{\mathrm{max}}\sim\left\\{12,6,3,2,1\right\\}$, for the
$t_{0}/M=\left\\{-5,0,5,10,15\right\\}$ curves respectively161616In the last
two cases, $t_{0}/M=10,15$, $\tilde{\epsilon}$ displays two local minima as a
function of $N_{\mathrm{max}}$ within the range considered, with the second
minimum being slightly lower than the first. In either case, both minima are
nevertheless very similar, so that a value of $\tilde{\epsilon}$ close to its
global minimum is already reached at the first local minimum. .

We observe that as the starting time $t_{0}$ increases, a lower number of
tones is required to get close to the minimum value of $\tilde{\epsilon}$. In
particular, the models with $N_{\mathrm{max}}=3,2,1$ are appropriate at
$t_{0}/M=5,10,15$ respectively. Furthermore, the fact noted above that the
minimum value reached by $\tilde{\epsilon}$ does not vary significantly with
the starting time $t_{0}$ — and remains consistent with the distribution of
the NR radiative error $\delta\epsilon_{\mathrm{r}}$ indicated by the orange
shaded area — is particularly intriguing, especially given that this still
applies to the curve obtained at $t_{0}/M=-5$. Notice that at negative times,
the amplitude of the strain $(2,2)$ mode is still increasing, so that the
morphology of the waveform at such an early stage still differs significantly
from a typical exponentially decaying ringdown wave (see for instance Figs.
16-22 of Bhagwat  _et al._ (2018)). On the other hand, the fact that
$\tilde{\epsilon}(N_{\mathrm{max}}=2)_{t_{0}/M=10}\simeq\tilde{\epsilon}(N_{\mathrm{max}}=12)_{t_{0}/M=-5}$
for instance, suggests that the lower tones $n\lesssim 2$ of the
$(N_{\mathrm{max}}=12)$ model may contribute significantly to the estimates of
$\epsilon$ obtained for $t_{0}/M=-5$ — with the addition of higher tones into
the model helping to better constrain the lower ones.

We justify this behavior relying on two main hypotheses: i) high overtones
could be fitting some fraction of the NR noise or ii) high-overtone amplitudes
and phases are flexible enough to fit well the early part of the waveform.
Hypothesis i) could only be better tested when more accurate waveforms are
added to the catalogs. On the other hand, hypothesis ii) would imply that the
high-tone degrees of freedom can capture the morphology of the waveform at
early times, while the values of $\epsilon$ result predominantly from
modelling increasingly better the lower tones $n=0,1..$. This could be the
case even irrespective of the NR errors. For instance, the short and similar
damping times of the high overtones could induce strong correlations between
the higher tones amplitudes combined to a very short time range where they are
still measurable, thus impeding one from getting physically reliable
information from them. The latter point is also reinforced by the results of
the amplitude stability analysis of Sec. IV.3.3.

This consideration is also relevant for deciding which starting time and model
should be preferred to estimate the final mass and the final spin of a given
GW event from its RD phase. As a rule of thumb, the statistical uncertainty on
a given parameter $\sigma_{\lambda}$ scales as the inverse of the RD signal-
to-noise ratio (SNR) $\rho\propto 1/\sigma_{\lambda}$. For an event consistent
with GW150914, about half of the total ringdown SNR is lost at $t_{0}/M=10$
and about $70\%$ at $t_{0}/M=15$ Bhagwat  _et al._ (2018, 2020); Jiménez
Forteza  _et al._ (2020). Therefore, given the minor variation of the minimum
$\tilde{\epsilon}$ with different values of $t_{0}$ observed here, we expect
that for real GW events, only the RD waves with a large number of overtones
$N_{\mathrm{max}}$, at early starting times and that are compatible with or
slightly dominated by the statistical error — i.e., with
$\sigma_{\lambda=M_{f},a_{f}}\gtrsim\epsilon(t_{0}/M,N_{\mathrm{max}})$ — will
be appropriate to place as accurate constraints on the mass and the spin as
possible.

Figure 10: We show the median value $\tilde{\epsilon}$ of $\epsilon$ as a
function of the number of overtones $N_{\mathrm{max}}$ in the RD model, for a
set of normalised fit starting times $t_{0}/M=\left\\{-5,0,5,10,15\right\\}$.
Notice that at $t_{0}/M=\left\\{5,10,15\right\\}$, the curves already nearly
hit their minimum at $N_{\mathrm{max}}\sim\left\\{3,2,1\right\\}$
respectively. We show again as a shaded orange band the 10–90 percentiles of
the distribution of the radiation error $\delta\epsilon_{\mathrm{r}}$ (cutting
part of the lowest values), as well as its median value as a black dashed
line.

##  V On the interpretation of the fit models

The magnitude $\epsilon$ measures systematic deviations on the recovered final
mass and final spin with respect to the true parameters of each NR simulation.
These deviations may be produced either by i) an insufficient number of tones
in our ringdown model, Eq. (1); ii) numerical errors propagated to the strain
$h_{22}(t)$; or iii) the possible presence of nonlinearities in the waveform
Bhagwat  _et al._ (2020); Okounkova (2020). We have observed in Figs. 6 and 7
that missing tones on Eq. (1) become the dominant source of deviations on
$(M_{f},a_{f})$ for models with $N_{\mathrm{max}}\lesssim 4-5$, and that such
deviations stabilize on average at $N_{\mathrm{max}}\gtrsim 6$. Furthermore,
there exists some nontrivial correlations between the values of
$\tilde{\epsilon}$ and the amplitudes of the tones. In particular, the
fundamental-mode best-fit amplitude $A_{0}$ stabilises at roughly the same
$N_{\mathrm{max}}$ as $\tilde{\epsilon}$ does. The amplitudes and phases of
the tones are not predicted from the solutions of black hole ringdown
perturbation theory but rather fixed by the initial conditions of each NR
simulation, _i.e._ , they do not hold any _a priori_ unique relation with the
final mass and final spin Berti  _et al._ (2006); London  _et al._ (2014);
Jiménez Forteza  _et al._ (2020). Therefore, the improvement on the
$(M_{f},a_{f})$ estimate at high $N_{\mathrm{max}}$ is achieved by adding
further information to the RD models through the complex frequencies
$\omega_{lmn}$, which depend uniquely on the mass and the spin. Concerning the
physical reliability of the amplitudes, we observe that the $n>1$ overtones
typically suffer from larger than $5\%$ variations on their values when
varying the number of overtones $N_{\mathrm{max}}$ of the model, in agreement
with Jiménez Forteza  _et al._ (2020).

As discussed in the main text, with $t_{0}=0$, the biases on $(M_{f},a_{f})$
decrease exponentially up to $N_{\mathrm{max}}\sim{5}$. The gain in accuracy
in the estimates of these parameters occurs irrespective of the stability
issues observed for the high-overtone amplitudes. Therefore, if one assumes GR
to be correct, and assuming a high enough SNR, those models shall provide
accurate estimates of the final mass and spin of real GW events, at least in a
majority of cases. In line with that argument, one could perform IMR
consistency tests, where any inconsistency observed on the final mass and spin
(as measured by $\epsilon$) between the inspiral and merger-ringdown regimes
could be hinting for new physics as it would suggest a violation of the no-
hair theorem. The results obtained in Giesler  _et al._ (2019); Finch and
Moore (2021) and complemented in this work, suggest that an QNM model with
${N_{\mathrm{max}}\sim 3-6}$ overtones would be able to constrain the final
parameters up to $\epsilon\sim 10^{-2}-10^{-3}$ in many cases (although some
signals would only reach such constraints with more overtones). This level of
accuracy is beyond the current LIGO-Virgo typical SNR-limited uncertainties on
the mass and the spin Abbott  _et al._ (2021b), which together result in
$\epsilon\sim 0.1$. However, this accuracy may be achieved and surpassed with
the third-generation detectors LISA, Einstein Telescope and Cosmic Explorer.

On the other hand, the stability issues observed on the amplitudes, together
with the variable $t_{0}$ analysis, become relevant in order to assess to
which extent these models and, in particular, high–overtone number models can
be used for performing black hole spectroscopy. This would imply estimating
_independently_ the frequencies and damping times of _each tone_ together with
phases and amplitudes. It is likely that the instabilities observed in the
tone amplitudes may become even larger when adding the frequencies and damping
times as extra free parameters. This would thus eventually induce systematic
errors that may have an important impact on the final estimate of the QNM
spectrum itself Bhagwat  _et al._ (2016); Capano  _et al._ (2021); Jiménez
Forteza  _et al._ (2021).

##  VI Conclusions

The aim of this work has been to study the behavior of high overtones for a
set of $610$ nonprecessing NR waveforms, extending the fit results obtained in
Giesler  _et al._ (2019); Finch and Moore (2021) to higher than
$N_{\mathrm{max}}=7$ overtones. To this end, we have computed the quasinormal
mode frequencies for the overtones with indices $n=8$ and $n=9$ both for the
corotating and the counter-rotating branches.

The $n=8$, $n=9$ corotating modes have been computed for a range of spins
$a_{f}\in\left[a_{f,n}^{\mathrm{start}},1\right]$ with
$a_{f,8}^{\mathrm{start}}=3.6\cdot 10^{-3}$ and
$a_{f,9}^{\mathrm{start}}=5.4\cdot 10^{-3}$, while the one counter-rotating
mode associated to these tones is provided for spins ranging from
$a_{f}=10^{-6}$ (extended by assumption to $a_{f}=0$ with the Schwarzschild
$\omega_{228}=-2\,\iota$ solution) to $0.997$. Our results are consistent with
Onozawa (1997); Berti  _et al._ (2003); Cook and Zalutskiy (2014, 2016a,
2016b) for these three branches, and are made available here: Mourier and
Jimenez Forteza , completing the results provided by Berti  _et al._ (2006,
2009); Berti ; Cardoso ; Cook and Zalutskiy (2014); Cook (2019); Stein (2019).

First, we have used these results to extend the RD fits to the SXS waveform
SXS:0305 to $N_{\mathrm{max}}>7$, as a function of the fit starting time
$t_{0}$, with the final mass and spin fixed to the simulation’s true final
parameters. We observe that the mismatch $\mathcal{M}$ keeps marginally
decreasing with $N_{\mathrm{max}}$ for the models with $N_{\mathrm{max}}=8,9$
and beyond. On this respect, we have found that the first local minimum in
mismatch $\mathcal{M}$ as a function of $t_{0}$ occurs at negative starting
times $t_{0}<0$ for $N_{\mathrm{max}}\geq 8$, which is possibly due to data
overfitting.

Second, we have estimated the value of the final mass – final spin recovery
bias $\epsilon$ for QNM models with
$N_{\mathrm{max}}\in\left\\{0,\dots,16\right\\}$ overtones, starting at the
peak of the $(2,2)$ strain component, $t_{0}=0$, and for two pairs of SXS and
RIT waveforms with identical parameters, $(a_{f},M_{f})=(0.692,0.952)$ for the
first pair {SXS:0305, RIT:0062} and $(a_{f},M_{f})=(0.581,0.966)$ for the
second pair {SXS:0259, RIT:0118}. We have found that the trend on $\epsilon$
can be significantly different between the two simulations and that it is, in
general, case-dependent. We estimate a very similar minimum $\epsilon$ of
$\epsilon_{\mathrm{min}}\sim 3\times 10^{-4}$ for the two SXS waveforms, at
$N_{\mathrm{max}}=7$ for SXS:0305 and at $N_{\mathrm{max}}=13$ for SXS:0259.
For the RIT waveforms with the same parameters, we have obtained
$\epsilon_{\mathrm{min}}\sim 4\times 10^{-4}$ at $N_{\mathrm{max}}=7$ and
$\epsilon_{\mathrm{min}}\sim 3\times 10^{-3}$ at $N_{\mathrm{max}}=6$, for
RIT:0062 and RIT:0118 respectively.

Then, we have applied the fitting algorithm described in Sec. IV.1 to 610 out
of the 620 SXS non-precessing binary black hole simulations and for
$N_{\mathrm{max}}\in\left\\{0,\dots,12\right\\}$, still with $t_{0}=0$.
Specifically, our results for $\epsilon$ are consistent with the
$N_{\mathrm{max}}=3,7$ models shown in Finch and Moore (2021). We observe that
the median value $\tilde{\epsilon}$ of the distribution of $\epsilon$ over
these simulations decreases exponentially with $N_{\mathrm{max}}$ to about
$\tilde{\epsilon}\simeq 2\times 10^{-3}$ at $N_{\mathrm{max}}\sim 5-6$ —
although there is a significant overlap between $N_{\mathrm{max}}=3,4,5,6$ in
the distribution of $\epsilon$ values between the multiple SXS cases.
Moreover, $\tilde{\epsilon}$ does not change significantly beyond
$N_{\mathrm{max}}\sim 6$, which also applies to our new $N_{\mathrm{max}}=8$,
$N_{\mathrm{max}}=9$ and $N_{\mathrm{max}}>9$ models. We noted nevertheless
that for about $21\%$ of the cases, models with $N_{\mathrm{max}}\geq 12$ were
required to hit the (similar) minimum value on $\epsilon$. The value of
$\tilde{\epsilon}$ appears to always be bounded by the NR errors. We provide
optimistic and pessimistic estimates of the NR errors which we have here
referred to respectively as i) local errors $\delta\epsilon_{\mathrm{l}}$ and
ii) radiative errors $\delta\epsilon_{\mathrm{r}}$. The latter should be a
more accurate representation of the NR on the strain since it is directly
derived from it. We notice in particular that our $90$th percentiles of the
radiation error distribution are above the plateau value $\tilde{\epsilon}\sim
2\times 10^{-3}$, hence they could be affecting or even dominating the
$\tilde{\epsilon}$ values at $N_{\mathrm{max}}\geq 5$.

Furthermore, we have studied the stability of the best-fit amplitude values
$A_{n}$ for a range of tones $n\in\left\\{0,...,4\right\\}$. For the
fundamental $n=0$ mode and for the first overtone $n=1$, we have found that
the median relative amplitudes variations $\delta\tilde{A}_{0}$ and
$\delta\tilde{A}_{1}$ between successive $N_{\mathrm{max}}$-overtone and
$(N_{\mathrm{max}}-1)$-overtone models are below the $1\%$ level for RD models
with $N_{\mathrm{max}}\geq 4$ and $N_{\mathrm{max}}\geq 8$, respectively. On
the other hand, we observe that the amplitudes of the overtones with $n>2$ are
unstable. We have also found a significant correlation between the typical
variation $\delta\tilde{A}_{0}$ of $A_{0}$ as a function of
$N_{\mathrm{max}}$, and the value of $\tilde{\epsilon}$. These elements could
indicate that, for a majority of the studied cases, the improvement on
$\epsilon$ is predominantly achieved by increasingly improving the constraints
on $A_{0}$ — and possibly on the first few overtones to a lesser extent.

Finally, we have repeated the $N_{\mathrm{max}}\in\\{0,\dots,12$} RD models
study over the 610 SXS cases considered over a few different values of the fit
starting time $t_{0}$, with $t_{0}/M\in\\{-5,0,5,10,15\\}$. We have found that
the minimum value reached by $\tilde{\epsilon}$ as a function of
$N_{\mathrm{max}}$ does not vary significantly if we vary the fit starting
time $t_{0}$; a model with $N_{\mathrm{max}}=7$ overtones at $t_{0}/M=0$
provides, on average, a similar (slightly higher) accuracy on $\epsilon$ as a
model with $N_{\mathrm{max}}=1$ at $t_{0}/M=15$. This is relevant since the
effects of the overtones $n\geq 1$ are expected to be small at $t_{0}/M\gtrsim
10$ due to their short damping times. This further supports the hypothesis
that the constraints on $\epsilon$ may be predominantly induced by an
improvement in modelling the low tones $n=0,1,\dots$, regardless of the weak —
or unstable — constraints one obtains for the higher tone amplitudes. In this
regard, the ansätze and, in particular, the higher tones, appear to be
sufficiently flexible to accurately fit the strain at times around the peak
even when those high overtones take amplitude values that are likely unstable,
or even nonphysical.

We note that our ringdown models Jimenez Forteza and Mourier  may be suitable
for performing IMR consistency tests for current and next-generation GW
observatories, where the gain in sensitivity may allow us to hit and surpass
the accuracy levels observed here for numerical data. On the other hand, based
on the results of this paper, we are more skeptical about using overtone
models of a given $(l,m)$ mode to robustly perform black hole spectroscopy due
to the instabilities observed on the amplitudes of the overtones with $n>1$ —
which are likely to propagate to frequencies and damping times when those are
added as free parameters. For the first overtone ($n=1$), such amplitude
instabilities are reduced to the $1\%$ level for models with
$N_{\mathrm{max}}\gtrsim 8$. Thus, a two-tone $n=0,1$ spectroscopy may remain
possible provided that one considers a large number of additional tones in the
model, at the expense of adding a large number of free parameters.

###### Acknowledgements.

We acknowledge the Max Planck Gesellschaft for support and we are grateful to
the Atlas cluster computing team at AEI Hannover for their help. The authors
are also thankful to Swetha Bhagwat, Collin Capano, Sumit Kumar, Alex Nitz and
Paolo Pani for useful discussions and comments on this paper. We also thank
the anonymous referee for insightful comments and suggestions that led to
further improvements to this document.

## References

  * Abbott  _et al._ (2019) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. X 9, 031040 (2019), arXiv:1811.12907 [astro-ph.HE] .
  * Abbott  _et al._ (2021a) R. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. X 11, 021053 (2021a), arXiv:2010.14527 [gr-qc] .
  * Giesler  _et al._ (2019) M. Giesler, M. Isi, M. A. Scheel, and S. A. Teukolsky, Phys. Rev. X 9, 041060 (2019), arXiv:1903.08284 [gr-qc] .
  * London  _et al._ (2014) L. London, D. Shoemaker, and J. Healy, Phys. Rev. D90, 124032 (2014), [Erratum: Phys. Rev.D94,no.6,069902(2016)], arXiv:1404.3197 [gr-qc] .
  * Bhagwat  _et al._ (2020) S. Bhagwat, X. J. Forteza, P. Pani, and V. Ferrari, Phys. Rev. D101, 044033 (2020), arXiv:1910.08708 [gr-qc] .
  * Leaver (1985) E. Leaver, Proc. Roy. Soc. Lond. A A402, 285 (1985).
  * Detweiler (1980) S. L. Detweiler, Astrophys. J. 239, 292 (1980).
  * Kokkotas and Schmidt (1999) K. D. Kokkotas and B. G. Schmidt, Living Rev. Rel. 2, 2 (1999), arXiv:gr-qc/9909058 [gr-qc] .
  * Finch and Moore (2021) E. Finch and C. J. Moore, Phys. Rev. D 103, 084048 (2021), arXiv:2102.07794 [gr-qc] .
  * Dhani (2020) A. Dhani, (2020), arXiv:2010.08602 [gr-qc] .
  * Berti  _et al._ (2009) E. Berti, V. Cardoso, and A. O. Starinets, Class. Quant. Grav. 26, 163001 (2009), arXiv:0905.2975 [gr-qc] .
  * Jiménez Forteza  _et al._ (2020) X. Jiménez Forteza, S. Bhagwat, P. Pani, and V. Ferrari, Phys. Rev. D 102, 044053 (2020), arXiv:2005.03260 [gr-qc] .
  * Cook (2020) G. B. Cook, Phys. Rev. D102, 024027 (2020), arXiv:2004.08347 [gr-qc] .
  * Dhani and Sathyaprakash (2021) A. Dhani and B. S. Sathyaprakash, (2021), arXiv:2107.14195 [gr-qc] .
  * Cook and Zalutskiy (2014) G. B. Cook and M. Zalutskiy, Phys. Rev. D 90, 124021 (2014), arXiv:1410.7698 [gr-qc] .
  * Abbott  _et al._ (2021b) R. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. D 103, 122002 (2021b), arXiv:2010.14529 [gr-qc] .
  * Capano  _et al._ (2021) C. D. Capano, M. Cabero, J. Westerweck, J. Abedi, S. Kastha, A. H. Nitz, A. B. Nielsen, and B. Krishnan, (2021), arXiv:2105.05238 [gr-qc] .
  * Isi  _et al._ (2019) M. Isi, M. Giesler, W. M. Farr, M. A. Scheel, and S. A. Teukolsky, Phys. Rev. Lett. 123, 111102 (2019), arXiv:1905.00869 [gr-qc] .
  * Isi and Farr (2021) M. Isi and W. M. Farr, (2021), arXiv:2107.05609 [gr-qc] .
  * London and Fauchon-Jones (2019) L. London and E. Fauchon-Jones, Class. Quant. Grav. 36, 235015 (2019), arXiv:1810.03550 [gr-qc] .
  * Berti and Klein (2014) E. Berti and A. Klein, Phys. Rev. D90, 064012 (2014), arXiv:1408.1860 [gr-qc] .
  * Berti  _et al._ (2006) E. Berti, V. Cardoso, and C. M. Will, Phys. Rev. D73, 064030 (2006), arXiv:gr-qc/0512160 [gr-qc] .
  * (23) E. Berti, https://pages.jh.edu/~eberti2/ringdown/.
  * (24) V. Cardoso, https://centra.tecnico.ulisboa.pt/network/grit/files/ringdown/.
  * Cook (2019) G. B. Cook, “Kerr quasinormal modes: s=-2, n=0–7 (version 1),” Zenodo, https://doi.org/10.5281/zenodo.2650358 (2019).
  * Stein (2019) L. C. Stein, J. Open Source Softw. 4, 1683 (2019), arXiv:1908.10377 [gr-qc] .
  * Onozawa (1997) H. Onozawa, Physical Review D 55, 3593–3602 (1997).
  * Berti  _et al._ (2003) E. Berti, V. Cardoso, K. D. Kokkotas, and H. Onozawa, Phys. Rev. D 68, 124018 (2003), arXiv:hep-th/0307013 [hep-th] .
  * Cook and Zalutskiy (2016a) G. B. Cook and M. Zalutskiy, Phys. Rev. D 94, 104074 (2016a), arXiv:1607.07406 [gr-qc] .
  * Cook and Zalutskiy (2016b) G. B. Cook and M. Zalutskiy, Classical and Quantum Gravity 33, 245008 (2016b), arXiv:1603.09710 [gr-qc] .
  * The SXS Collaboration (2016) The SXS Collaboration, “SXS Gravitational Waveform Database,”  (2016).
  * Campanelli  _et al._ (2016) M. Campanelli, J. Healy, C. Lousto, and Y. Zlochower, “CCRG@RIT Catalog of Numerical Simulations,”  (2016).
  * Teukolsky and Press (1974) S. A. Teukolsky and W. H. Press, Astrophys. J. 193, 443 (1974).
  * Chandrasekhar (1985) S. Chandrasekhar,  _The mathematical theory of black holes_ (Oxford Classic Texts in the Physical Sciences, 1985).
  * Chandrasekhar and Detweiler (1975) S. Chandrasekhar and S. L. Detweiler, Proc. Roy. Soc. Lond. A344, 441 (1975).
  * (36) P. Mourier and X. Jimenez Forteza, https://codeberg.org/GW_Ringdown.
  * (37) X. Jimenez Forteza and P. Mourier, https://github.com/frcojimenez/GW_Rdown.
  * Maassen van den Brink (2000) A. Maassen van den Brink, Phys. Rev. D 62, 064009 (2000), arXiv:gr-qc/0001032 [gr-qc] .
  * Rezzolla  _et al._ (2008) L. Rezzolla, E. Barausse, E. N. Dorband, D. Pollney, C. Reisswig, J. Seiler, and S. Husa, Phys. Rev. D78, 044002 (2008), arXiv:0712.3541 [gr-qc] .
  * Jiménez-Forteza  _et al._ (2017) X. Jiménez-Forteza, D. Keitel, S. Husa, M. Hannam, S. Khan, and M. Pürrer, Phys. Rev. D95, 064024 (2017), arXiv:1611.00332 [gr-qc] .
  * Hofmann  _et al._ (2016) F. Hofmann, E. Barausse, and L. Rezzolla, Astrophys. J. 825, L19 (2016), arXiv:1605.01938 [gr-qc] .
  * Szilagyi  _et al._ (2009) B. Szilagyi, L. Lindblom, and M. A. Scheel, Phys. Rev. D80, 124010 (2009), arXiv:0909.3557 [gr-qc] .
  * Ashtekar and Krishnan (2004) A. Ashtekar and B. Krishnan, Living Rev. Rel. 7, 10 (2004), arXiv:gr-qc/0407042 [gr-qc] .
  * Iozzo  _et al._ (2021) D. A. Iozzo, N. Khera, L. C. Stein, K. Mitman, M. Boyle, N. Deppe, F. Hébert, L. E. Kidder, J. Moxon, H. P. Pfeiffer, and et al., Physical Review D 103 (2021), 10.1103/physrevd.103.124029.
  * Christodoulou (1970) D. Christodoulou, Phys. Rev. Lett. 25, 1596 (1970).
  * Keitel  _et al._ (2017) D. Keitel, X. Jimenez-Forteza, _et al._ , Phys. Rev. D96, 024006 (2017), arXiv:1612.09566 [gr-qc] .
  * Hinder  _et al._ (2014) I. Hinder _et al._ , Class. Quant. Grav. 31, 025012 (2014), arXiv:1307.5307 [gr-qc] .
  * Hinder  _et al._ (2011) I. Hinder, B. Wardell, and E. Bentivegna, Phys. Rev. D84, 024036 (2011), arXiv:1105.0781 [gr-qc] .
  * Mourier  _et al._ (2021) P. Mourier, X. Jiménez Forteza, D. Pook-Kolb, B. Krishnan, and E. Schnetter, Phys. Rev. D 103, 044054 (2021), arXiv:2010.15186 [gr-qc] .
  * Bhagwat  _et al._ (2018) S. Bhagwat, M. Okounkova, S. W. Ballmer, D. A. Brown, M. Giesler, M. A. Scheel, and S. A. Teukolsky, Phys. Rev. D97, 104065 (2018), arXiv:1711.00926 [gr-qc] .
  * Okounkova (2020) M. Okounkova, (2020), arXiv:2004.00671 [gr-qc] .
  * Bhagwat  _et al._ (2016) S. Bhagwat, D. A. Brown, and S. W. Ballmer, Phys. Rev. D94, 084024 (2016), [Erratum: Phys. Rev.D95,no.6,069906(2017)], arXiv:1607.07845 [gr-qc] .
  * Jiménez Forteza  _et al._ (2021) X. Jiménez Forteza _et al._ , (2021), in preparation.

##  Appendix A Outliers

Index | Issue
---|---
0002 |  Large extrapolation error
0084 | Large extrapolation error
0090 | Large extrapolation error
0091 | Large extrapolation error
0158 | Large extrapolation error
0170 | Reported $M_{f}>1$
0171 | Reported $M_{f}>1$
0218 | Large extrapolation error
1110 | Large extrapolation error
1134 | Reported $M_{f}>1$
Table 1: List of the $\fpeval{620-610}$ out of 620 non-precessing SXS binary-
black hole waveforms (labeled under the form SXS:BBH:_index_) that we do not
consider in our analysis. These waveforms are excluded either due to an
unphysical value being reported for $M_{f}$ ($M_{f}>1$) in their respective
metadata files, or due to a large extrapolation error (as measured by a
mismatch value $\mathcal{M}\geq 10^{-3}$ between the waveforms provided at
extrapolation orders $N=2$ and $N=3$), as per indicated in the second column.

(a)

(b)

Figure 11: Relative variation of the best-fit amplitudes
$\delta\tilde{A}_{0}(N_{\mathrm{max}})$ and
$\delta\tilde{A}_{1}(N_{\mathrm{max}})$ at $N_{\mathrm{max}}=5$ (left panel)
and at $N_{\mathrm{max}}=12$ (right panel) as a function of the starting time
$t_{0}/M$. The solid curves represent the median values of the 610 cases
studied in this work while the shadowed areas stand for the 10–90 percentiles
of the corresponding distributions.

We list in Table 1 the $\fpeval{620-610}$ non-precessing SXS waveforms that we
have removed from our analysis. We used the mismatch as a landmark to identify
the cases with larger uncertainty. In Fig. 3 we have seen that the dominant
contribution to the numerical uncertainty comes from the extrapolation of the
waveform to null infinity. Accordingly, for each of the non-precessing SXS
cases, we have computed the mismatch $\mathcal{M}$ as given by Eq. (12)
between the waveforms with successive extrapolation orders $N=2$ and $N=3$. To
this end, we aligned beforehand the two numerical waveforms
$h_{22}^{N=2,3}(t)$ in time and in phase so that the peak of the strain is
located at time $t=0$ for both, with the same initial phase. We then excluded
the cases for which we found an extrapolation mismatch $\mathcal{M}\geq
10^{-3}$. We moreover excluded another three cases with a a seemingly
incorrect reported final mass $M_{f}>1$ (with in fact even $M_{f}>2$ in each
of these cases).

##  Appendix B Stability of the fit amplitudes at varying starting times

In this section we study the respective fractional variations $\delta A_{0}$,
$\delta A_{1}$ of the best-fit amplitudes of the fundamental mode and of the
first overtone between successive models with $N_{\mathrm{max}}$ and
$N_{\mathrm{max}}-1$ overtones — as in Fig. 9 — at different starting times
$t_{0}/M$ and for two separate values of $N_{\mathrm{max}}$. The first value
chosen is $N_{\mathrm{max}}=5$, _i.e._ , a value at which
$\delta\tilde{A}_{0}$ has not yet plateaued at $t_{0}/M=0$. The second value
is $N_{\mathrm{max}}=12$ for which we expect the values of
$\delta\tilde{A}_{0}$ to have plateaued for all the starting times considered
in this work, and $\delta\tilde{A}_{1}$ to have also reached a plateau at
$t_{0}/M=0$. In Fig. 11 we show the results for the fractional variations (in
percent), with $N_{\mathrm{max}}=5$ on the left panel and
$N_{\mathrm{max}}=12$ on the right panel. The solid curves represent the
median values $\delta\tilde{A}_{n}$ as shown in Fig. 9 while the shadowed
areas stand for the the 10–90 percentiles of the distributions over the 610
SXS cases considered in this work. The resolution on $t_{0}/M$ in this figure
is $\Delta(t_{0}/M)=2.5$.

For $N_{\mathrm{max}}=5$, we obtain that $\delta\tilde{A}_{0}$ decreases with
$t_{0}$ until $t_{0}/M\simeq 2.5$, beyond which it approximately stabilizes at
$\delta\tilde{A}_{0}\sim 0.3\%$ (slightly increasing again towards later
times). The larger values at earlier times $t_{0}/M<2.5$ are due to the need
for a larger number of tones to accurately fit the data — while the
insufficient number of tones induces the further instabilities observed on
$\delta\tilde{A}_{0}$ (_cf._ Fig. 9 at $t_{0}/M=0$). On the other hand, the
values for $\delta\tilde{A}_{1}(N_{\mathrm{max}}=5)$ decay until reaching a
minimum of about $2\%$ at $t_{0}/M\sim 7.5$ to thereafter grow back at late
times. As for the case of the $n=0$ amplitude, the early decreasing of
$\delta\tilde{A}_{1}$ is sourced by the lack of higher overtones in the RD
model, while the late increase may be due to the effects of the NR noise
together with the suppression of the first overtone amplitude at late times
(see the further discussion of this tone below).

For $N_{\mathrm{max}}=12$, the values of $\delta\tilde{A}_{0}$ remain
approximately constant (or very slightly increasing), around $0.1\%$ to
$0.2\%$, for all the starting times studied, where the little variations on
their values are consistent with the values of the dispersion of the data
given by the 10–90 percentiles. We recover approximately the same values as in
the large–$N_{\mathrm{max}}$ plateau regime as a function of
$N_{\mathrm{max}}$ shown for $t_{0}/M=0$ in Fig. 9. Like for
$N_{\mathrm{max}}=5$, we observe an increasing trend for
$\delta\tilde{A}_{1}(N_{\mathrm{max}}=12)$, at $t_{0}/M>0$. Since here the
number of tones is much larger, this indicates that the $n=1$ tone is less
stable compared to $n=0$ at late times regardless of the number of tones
included in our models. This lack of stability at late times may be originated
by the suppression of the first overtone at these times, that is sourced by
its high damping factor $\exp(-t/\tau_{1})$. Recall that for instance, for a
GW150914–like event, $\tau_{1}/M\sim 3.8$ and the first overtone amplitude is
reduced by $\sim{e^{-4}}$ between $t/M=0$ and $t/M\sim 15$, thus leaving the
amplitude values of this tone more exposed to the effects of the NR noise.