††thanks: NHFP Einstein fellow

# Analyzing black-hole ringdowns

Maximiliano Isi<EMAIL_ADDRESS>LIGO Laboratory, Massachusetts Institute of
Technology, Cambridge, Massachusetts 02139, USA  Center for Computational
Astrophysics, Flatiron Institute, 162 5th Ave, New York, NY 10010  Will M.
Farr<EMAIL_ADDRESS>Center for Computational Astrophysics, Flatiron
Institute, 162 5th Ave, New York, NY 10010  Department of Physics and
Astronomy, Stony Brook University, Stony Brook NY 11794, USA

###### Abstract

A perturbed black hole rings down by emitting gravitational waves in tones
with specific frequencies and durations. Such tones encode prized information
about the geometry of the source spacetime and the fundamental nature of
gravity, making the measurement of black hole ringdowns a key goal of
gravitational wave astronomy. However, this task is plagued by technical
challenges that invalidate the naive application of standard data analysis
methods and complicate sensitivity projections. In this paper, we provide a
comprehensive account of the formalism required to properly carry out ringdown
analyses, examining in detail the foundations of recent observational results,
and providing a framework for future measurements. We build on those insights
to clarify the concepts of ringdown detectability and resolvability—touching
on the drawbacks of both Bayes factors and naive Fisher matrix approaches—and
find that overly pessimistic heuristics have led previous works to
underestimate the role of ringdown overtones for black hole spectroscopy. We
put our framework to work on the analysis of a variety of simulated signals in
colored noise, including analytic injections and a numerical relativity
simulation consistent with GW150914. We demonstrate that we can use tones of
the quadrupolar angular harmonic to test the no-hair theorem at current
sensitivity, with precision comparable to published constraints from real
data. Finally, we assess the role of modeling systematics, and project
measurements for future, louder signals. We release ringdown, a Python library
for analyzing black hole ringdowns using the the methods discussed in this
paper, under a permissive open-source license.

##  I Introduction

Black-hole (BH) ringdowns provide an exceptional observational handle on the
nature of gravity: a perturbed BH radiates gravitational waves (GWs) in the
form of damped sinusoids that cleanly encode information about the structure
of the BH spacetime [1, 2, 3, 4], and make it accessible to LIGO [5], Virgo
[6] and future GW detectors [7, 8, 9]. The research program known as “BH
spectroscopy” [10, 11, 12] aims to analyze such quasinormal modes (QNMs) in
order to test general relativity (GR) and the BH paradigm [13, 14], as
embodied by the Kerr metric for astrophysical BHs [15, 16]. This includes
basic tenets such as the no-hair theorem and the area law—respectively, the
statements that astrophysical BHs are fully characterized by their mass and
spin [17, 18, 19, 20, 21, 22], and that the total area of classical BH
horizons may not decrease over time [23, 24, 25].

Although the literature on BH spectroscopy is vast (e.g., [11, 12, 26, 27, 28,
29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47]),
sensitivity projections have often relied on simplistic data analysis
treatments that reduce the scope of their conclusions. For instance, analyses
based on Fisher matrices tend to assume that a test of the Kerr hypothesis
would require _independently_ analyzing two or more QNMs in order to obtain
separate estimates of the BH parameters to be compared for consistency, e.g.,
visually on a plot [11, 12, 26, 30, 40, 42]. Although useful for projections
in some regimes, such a setup does not reflect how an optimal measurement
would be carried out: in order to verify that a given ringdown signal is
consistent with a Kerr spectrum, we should self-consistently fit for multiple
modes _simultaneously_ under a Bayesian framework [27, 28, 29, 35, 36, 37, 38,
39]. The isolated-mode treatment, together with overly pessimistic heuristics,
has led to misconceptions—for example, that BH overtones (fast-decaying modes
with similar frequencies) should only be spectroscopically useful for
uncommonly loud detections [40, 42], in contradiction with observational
results obtained with modest signal amplitudes [39, 44].

The prevalence of oversimplifications in ringdown studies is hardly
surprising: formulating a robust analysis to extract QNMs from GW data is an
unexpectedly challenging task. Currently, our best chance for gaining access
to a BH ringdown is to target the late stage of GWs from binary BH
coalescences, which consists of emission by the perturbed remnant BHs born in
the merger. Operationally, the ringdown portion of such a signal is defined to
begin at some time $t_{0}$, after which the strain can be accurately described
as a superposition of damped sinusoids with complex frequencies corresponding
to the QNMs of the final BH [34].111In other words, we _define_ “ringdown” as
the portion of the signal that can be accurately described as a superposition
of damped sinusoids (we do not consider polynomial tails [48]). This is a
phenomenological statement that establishes our data analysis target, not
necessarily a statement about BH physics. Thus, studies aiming to isolate the
ringdown should be able to probe data at times $t\geq t_{0}$ irrespective of
what preceded it. The presence of this time-domain distinction introduces
subtle complications to ringdown studies in the presence of “colored” (i.e.
time-correlated) noise that invalidate the naive application of standard data
analysis strategies [29, 31, 37, 39, 47].

There are two broad avenues for addressing this challenge: (i) enhance regular
inspiral-merger-ringdown (IMR) waveform models with additional freedom in the
ringdown stage [35, 49, 46]; or (ii) define a model only for the late data,
and discard the rest [29, 36, 37, 39, 47]. The former is simpler to
accommodate within the standard infrastructure for LIGO-Virgo analyses,
especially when starting from a waveform whose ringdown model is explicitly
based on perturbation theory [35, 46]. However, it has the disadvantage of
necessarily coupling the ringdown measurement to the inspiral-merger regime,
thus becoming susceptible to systematics induced by the attachment procedure,
or shortcomings of the pre-ringdown model.

On the other hand, having an independent model for the post-merger data is
conceptually cleaner, and allows for a ringdown analysis fully agnostic about
the inspiral dynamics. However, this approach demands truncating the GW signal
at a specific time, which is difficult to handle with the usual LIGO-Virgo
analysis techniques [50]. Instead, it calls for special treatment in the time
domain [39, 37], or an equivalent nontrivial procedure in the frequency domain
[47]. There is also potential uncertainty about the optimal truncation point,
although this choice can be reasonably informed by numerical studies
[38].222In fact, this issue is also implicitly present in the enhanced-
waveform strategy mentioned above, in the choice of ringdown start time within
the IMR waveform model itself [35, 51].

The main goal of this paper is to lay down, and then empirically validate, the
formalism for properly carrying out this second type of analysis, examining in
detail the foundations of recent observational results [39, 52, 53, 44, 45].
This involves a thorough description of how to construct a suitable ringdown-
only template for BH spectroscopy, as well as of the Bayesian infrastructure
and data manipulations required to use it without bias. The end result is a
comprehensive, self-contained account of ringdown analyses that expands,
clarifies and, when necessary, corrects previous literature, under a uniform
notation and set of conventions. Consequently, this paper synthesizes the
conceptual and methodological backbone of past observations in [39] and others
like it, elucidating the role of different assumptions and data analysis
choices made in the literature. More importantly, it provides a flexible and
robust platform for future, richer measurements.

The second goal is to build on those insights to revisit the concepts of
ringdown detectability and resolvability. To that end, we formulate a
procedure for operationally establishing the number of ringdown modes required
by certain data, in a way that is not strongly sensitive to arbitrary choices
of prior bounds (in contrast with, e.g., approaches based exclusively on Bayes
factors). With this in hand, we examine what it means to “resolve” two
ringdown modes, and the conditions required to spectrally characterized them.
We show that Rayleigh-like spectral separation criteria can be easily
misapplied to obtain overly pessimistic projections, as we argue has been the
case in some previous studies of BH spectroscopy.

The third and final goal of the paper is to showcase our overall conceptual
framework, as well as our specific implementation, by applying it to the
question of BH spectroscopy with overtones. Studying GW150914-like ringdowns
in colored noise, we show how our method is successful at reconstructing the
spectral and polarization properties of the signal. We then evaluate the
detectability of the first overtone, and demonstrate its potential to improve
measurements of BH parameters and, more crucially, to constrain deviations
away from the Kerr spectrum. Through a numerical relativity injection, we
confirm once more that ringdowns are detectable and resolvable with current
detectors, and that overtones offer a valuable resource for tests of the no-
hair theorem.

This paper is intended as a comprehensive reference that can be read in its
entirety or by selecting specific sections. We begin in Sec. II by showing
that the most generic ringdown template suitable for data analysis consists of
a superposition of elliptically polarized damped sinusoids, and that nuances
like spherical-spheroidal mixing are overall inconsequential; we also examine
key properties of the Kerr spectrum, with a focus on tones of the (usually
dominant) quadrupolar harmonic. In Sec. III, we discuss the technical
challenges intrinsic to ringdown-only analyses and show how to circumvent them
by operating fully in the time domain; we show how to compute the matched-
filter signal-to-noise ratio (SNR) of a signal in noisy data. In Sec. IV, we
demonstrate the efficacy of the method and show how to use it to detect and
characterize damped sinusoids; importantly, we show modes need not be
separated in frequency to be resolvable. In Sec. V, we turn to more realistic
simulations of BH ringdowns, including several tones of the quadrupolar
angular harmonic of a GW150914-like system [54]; we demonstrate the recovery
of analytic Kerr and non-Kerr injections at different SNRs, and replicate some
of those results with a numerical-relativity injection. In Sec. VI, we offer
an itemized summary of the main conclusions in the paper, and outline future
work.

As an important companion to this paper, we simultaneously release a software
package, which we name ringdown [55]. This software implements our analysis
framework within a flexible platform that users can use to easily develop
their own ringdown analyses. The code is written in the Python language, with
a Bayesian sampling component based in Stan [56, 57].

##  II Ringdown model

###  II.1 Generic template

There exist two qualitatively distinct families of Kerr QNMs: prograde modes,
corotating with the hole, and retrograde modes, counterotating with the hole
[58, 59].333It is standard for theoretical QNM studies to disregard modes with
negative frequency, in which case negative (positive) $m$ becomes synonymous
with “retrograde” (“prograde”) modes; we do not adopt that terminology. Taking
both types into account and setting $t_{0}=0$, the most generic template based
on Kerr perturbation theory for the complex-valued ringdown strain, $h\equiv
h_{+}-ih_{\times}$, is

$h=\sum C_{[p]\ell mn}\,e^{-i\tilde{\omega}_{[p]\ell
mn}t}\,{}_{\scalebox{0.55}{$-$}2}S_{[p]\ell mn}(\iota,\varphi)\,,$ (1)

summing over $p=+1$ for prograde modes and $p=-1$ for retrograde modes,
azimuthal number $\ell\geq 2$, magnetic number $-\ell\leq m\leq\ell$, and
overtone number $0\leq n$. Each QNM carries a complex frequency
$\tilde{\omega}_{[p]\ell mn}\equiv\omega_{[p]\ell mn}-i\gamma_{[p]\ell mn}$,
with damping time $\tau_{[p]\ell mn}\equiv 1/\gamma_{[p]\ell mn}$. The
overtone number $n$ is defined to order modes of a given $p$, $\ell$ and $m$
by decreasing damping time, such that $\tau_{[p]\ell mn}>\tau_{[p]\ell
m(n+1)}$. The mode amplitudes and phases are encoded in the complex amplitudes
$C_{[p]\ell mn}$, which are set by the initial state of the perturbation and
are not known a priori [60, 26].

The angular structure of the strain is given by the spin-weighted spheroidal
harmonics ${}_{\scalebox{0.55}{$-$}2}S_{[p]\ell mn}$, as a function of polar
and azimuthal angles $\iota$ and $\phi$ relative to the BH spin direction,

${}_{\scalebox{0.55}{$-$}2}S_{[p]\ell mn}(\iota,\varphi)\equiv
e^{im\varphi}{}_{\scalebox{0.55}{$-$}2}S_{\ell
m}(\chi\,t_{M}\tilde{\omega}_{[p]\ell mn},\cos\iota)\,,$ (2)

with a nontrivial dependence on the product of the dimensionless spin
magnitude, $\chi$, and the dimensionless complex frequency,
$t_{M}\tilde{\omega}_{[p]\ell mn}\equiv(GM/c^{3})\,\tilde{\omega}_{[p]\ell
mn}$, for BH mass $M$ [3, 61, 58, 62, 63]. These functions are not orthogonal
over the sphere, causing the modes to mix: the projection of the strain of
mode $\\{[p]\ell mn\\}$ onto the corresponding spherical harmonic
${}_{\scalebox{0.55}{$-$}2}Y_{\ell m}(\iota,\varphi)$ will, in principle, pick
up contributions from an infinite set of modes
$\\{[p^{\prime}]\ell^{\prime}mn^{\prime}\\}$ sharing the same magnetic number
$m$ [3, 64]. This can be important in the study of numerical relativity
waveforms (e.g., [65, 38]), wherein the strain is extracted from the
simulation through a projection into the ${}_{\scalebox{0.55}{$-$}2}Y_{\ell
m}$’s (which are themselves complete and orthogonal over the sphere, unlike
the ${}_{\scalebox{0.55}{$-$}2}S_{\ell m}$’s) [66]. One can address this by
expanding the ${}_{\scalebox{0.55}{$-$}2}S_{\ell m}$ as a series of
${}_{\scalebox{0.55}{$-$}2}Y_{\ell m}$’s, most often finding that
${}_{\scalebox{0.55}{$-$}2}S_{\ell m}\approx{}_{\scalebox{0.55}{$-$}2}Y_{\ell
m}$ is a sufficiently good approximation [62]. Either way, as will become
apparent below, none of this is relevant for our purposes because we will not
be interested in measuring the $C_{[p]\ell mn}$’s.

By definition, the prograde versus retrograde distinction implies
$\mathrm{sgn}(m)=p~{}\mathrm{sgn}(\omega_{[p]\ell mn})$ for $m\neq 0$.444For
$m=0$, there are still two families of modes indexed by $p=\pm 1$, but the
prograde vs retrograde terminology loses meaning. Although both prograde
($p=+1$) and retrograde ($p=-1$) modes are equally fundamental in principle,
studies have suggested the latter to be suppressed for regular binaries in
which the BHs rotate in the same sense as the orbit [12, 60, 65, 67, 68]
(however, see [69, 70]). Based on this, it is standard to focus on prograde
modes and set $C_{[-1]\ell mn}=0$. We will do so here and drop the $p$ index
below, even though our formalism can trivially accommodate retrograde modes.

We can further narrow the scope of Eq. (1) by appealing to symmetry. Parity-
time (PT) symmetry of the baseline (unperturbed) metric implies $\omega_{\ell
mn}=-\omega_{\ell\scalebox{0.55}{$-$}mn}$; indeed, in the Kerr case,
$\tilde{\omega}^{\rm(GR)}_{\ell
mn}=-\tilde{\omega}^{\rm(GR)*}_{\ell\scalebox{0.55}{$-$}mn}$, where ∗
indicates complex conjugation [3]. Assuming this holds and implicitly setting
$p=1$ everywhere, we may rewrite the sum in Eq. (1) as555This double counts
$m=0$ modes, which can be solved by redefining $C_{\ell 0n}\to C_{\ell 0n}/2$.

$\displaystyle h=\sum_{\ell}\sum_{0\leq m\leq\ell}\sum_{n}$
$\displaystyle\left[C_{\ell mn}e^{-i\tilde{\omega}_{\ell
mn}t}{}_{\scalebox{0.55}{$-$}2}S_{\ell mn}(\iota,\varphi)\,+\right.$
$\displaystyle\phantom{[}\left.C_{\ell\scalebox{0.55}{$-$}mn}e^{i\tilde{\omega}_{\ell
mn}^{*}t}{}_{\scalebox{0.55}{$-$}2}S_{\ell\scalebox{0.55}{$-$}mn}(\iota,\varphi)\right],$
(3)

where
${}_{\scalebox{0.55}{$-$}2}S_{\ell\scalebox{0.55}{$-$}mn}(\iota,\varphi)=(-1)^{\ell}{}_{\scalebox{0.55}{$-$}2}S^{*}_{\ell
mn}(\pi-\iota,\varphi)$ by symmetry [61, 63].666This property is sometimes
presented without writing out the angular dependence (e.g., [12, 65]), which
could be incorrectly taken to mean that
${}_{\scalebox{0.55}{$-$}2}S_{\ell\scalebox{0.55}{$-$}mn}$ and
${}_{\scalebox{0.55}{$-$}2}S^{*}_{\ell mn}$ are directly interchangeable (App.
A). Following this reasoning, we will always consider the two $\pm m$ terms in
the summand of Eq. (II.1) as a set, denoting them jointly by $h_{\ell|m|n}$.
(See also App. D of [71].)

In the absence of a physical model for the $C_{\ell mn}$’s, we can simply
absorb the ${}_{\scalebox{0.55}{$-$}2}S_{\ell mn}$ factors into redefined
complex amplitudes $C_{\ell
mn}^{\prime}\equiv{}_{\scalebox{0.55}{$-$}2}S_{\ell mn}(\iota,\varphi)C_{\ell
mn}$, so that the summand in Eq. (II.1) becomes

$\displaystyle h_{\ell|m|n}$ $\displaystyle=C_{\ell
mn}^{\prime}e^{-i\tilde{\omega}_{\ell
mn}t}+C_{\ell\scalebox{0.55}{$-$}mn}^{\prime}e^{i\tilde{\omega}^{*}_{\ell
mn}t}$ $\displaystyle=\left[C_{\ell mn}^{\prime}e^{-i\omega_{\ell
mn}t}+C_{\ell\scalebox{0.55}{$-$}mn}^{\prime}e^{i\omega_{\ell
mn}t}\right]e^{-t/\tau_{\ell mn}}.$ (4)

This expression makes it clear that assuming $\tilde{\omega}_{\ell
mn}=-\tilde{\omega}^{*}_{\ell\scalebox{0.55}{$-$}mn}$ was equivalent to
requiring that each pair of prograde $\\{\ell,\pm m,n\\}$ modes add up to a
single elliptically polarized contribution, to which we refer as “the
$\\{\ell,|m|,n\\}$ mode.”

To make the elliptical character of $h_{\ell|m|n}$ explicit, factor the new
complex amplitudes $C^{\prime}$ as

$C_{\ell\pm|m|n}^{\prime}=\frac{1}{2}\left(1\pm\epsilon_{\ell|m|n}\right)A_{\ell|m|n}\,e^{i\phi_{\ell\pm|m|n}},$
(5)

for some arbitrary real amplitude $A_{\ell|m|n}$, an ellipticity
$-1\leq\epsilon_{\ell|m|n}\leq 1$, and two independent phases
$\phi_{\ell\pm|m|n}$. The inverse transformation is simply

$A_{\ell|m|n}=|C_{\ell+|m|n}^{\prime}|+|C_{\ell-|m|n}^{\prime}|\,,$ (6)
$\epsilon_{\ell|m|n}=\frac{|C_{\ell+|m|n}^{\prime}|-|C_{\ell-|m|n}^{\prime}|}{|C_{\ell+|m|n}^{\prime}|+|C_{\ell-|m|n}^{\prime}|},$
(7)

and $\phi_{\ell\pm|m|n}=\arg({C_{\ell\pm|m|n}^{\prime}})$, so there is no loss
of generality. To simplify the indices, define $j\equiv\\{\ell,|m|,n\\}$ and
$\phi_{\pm j}\equiv\phi_{\ell\pm|m|n}$, so that Eq. (II.1) becomes

$\displaystyle h_{j}=\frac{1}{2}A_{j}\,e^{-t/\tau_{j}}$
$\displaystyle\left[\left(1+\epsilon_{j}\right)e^{-i(\omega_{j}t-\phi_{+j})}~{}+\right.$
$\displaystyle\phantom{[}\left.\left(1-\epsilon_{j}\right)e^{i(\omega_{j}t+\phi_{-j})}\right].$
(8)

It is now clear that this is an elliptical GW: $\epsilon_{j}=1$ yields a
right-hand circularly polarized mode; $\epsilon_{j}=-1$ yields a left-hand
circularly polarized mode; and intermediate values yield generic, elliptically
polarized waves, including linear polarizations for $\epsilon_{j}=0$.

When analyzing GW detector data, it is convenient to work in the linear
polarization basis. The plus and cross polarizations corresponding to Eq.
(II.1), $h^{(+)}_{j}=\Re\left[h_{j}\right]$ and
$h^{(\times)}_{j}=-\Im\left[h_{j}\right]$, can be written

$h^{(+)}_{j}=h^{c}_{j}\,\cos\theta_{j}-\epsilon_{j}h^{s}_{j}\,\sin\theta_{j}\,,$
(9)
$h^{(\times)}_{j}=h^{c}_{j}\,\sin\theta_{j}+\epsilon_{j}h^{s}_{j}\,\cos\theta_{j}\,,$
(10)

for cosine and sine quadratures

$h^{c}_{j}\equiv A_{j}\,e^{-t/\tau_{j}}\cos(\omega_{j}t-\phi_{j})\,,$ (11)
$h^{s}_{j}\equiv A_{j}\,e^{-t/\tau_{j}}\sin(\omega_{j}t-\phi_{j})\,.$ (12)

Here, we have defined two new angles:
$\theta_{j}\equiv-(\phi_{+j}+\phi_{-j})/2$ and
$\phi_{j}\equiv(\phi_{+j}-\phi_{-j})/2$.

The role of these quantities is illustrated in Fig. 1. At time $t=0$, the
polarization state of $h_{j}$ is given by a phasor lying somewhere along an
ellipse in the $h^{(+)}_{j}$, $h^{(\times)}_{j}$ plane. The ellipse has
semimajor and semiminor axes $A_{j}$ and $\epsilon A_{j}$, respectively,
oriented such that the semimajor axis defines an angle $\theta_{j}$ with the
$h^{(+)}_{j}$ axis. Finally, $\phi_{j}$ is the angle between the initial state
and the semimajor axis within the ellipse. As time passes, the phasor circles
the ellipse with angular velocity $\omega_{j}$, while the size of the ellipse
itself shrinks at an exponential rate with $e$-folding time $\tau_{j}$.

Figure 1: Polarization ellipse. At any given time, the phasor of a given
$j\equiv\\{\ell,|m|,n\\}$ ringdown mode lies on an ellipse with amplitude
$A_{j}$ and ellipticity $\epsilon_{j}$, with semimajor axis tilted by an angle
$\theta_{j}$ with respect to the plus-polarization axis (abscissa); a second
angle, $\phi_{j}$, determines the initial location of the phasor within the
ellipse. The phasor circles the ellipse with angular frequency $\omega_{j}$,
while the ellipse itself shrinks with $e$-folding time $\tau_{j}$.

Equations (9–12) refer to a single mode $j=\\{\ell,|m|,n\\}$. To construct a
template for the signal $s_{I}$ recorded by GW detector $I$, we need only add
all modes under consideration and project each polarization onto the detector
through the corresponding antenna patterns $F_{+/\times}^{I}$,

$s_{I}=\sum_{j}\left[F_{+}^{I}h^{(+)}_{j}(t-\delta
t_{I})+F_{\times}^{I}h^{(\times)}_{j}(t-\delta t_{I})\right].$ (13)

The $F^{I}_{+/\times}$ encode the relative orientation of the detector with
respect to the GW strain tensor, and so depend on the source right ascension
$\alpha$, declination $\delta$ and polarization angle $\psi$.777These response
functions vary only on the timescale of a sidereal day, so we can ignore their
time dependence when analyzing ringdowns. However, for fixed $\alpha$ and
$\delta$, the effect of $\psi$ is fully degenerate with a rotation of the
polarization ellipse through $\theta_{j}$, so we can evaluate the
$F^{I}_{+/\times}$ at an arbitrary value of $\psi$ without loss of generality
(see, e.g., App. A in [72]). The $\delta t_{I}$ represent time-of-flight
delays accounting for the different signal arrival times at each detector, and
are usually defined with respect to the geocenter; these are also functions of
$\alpha$ and $\delta$, but not $\psi$.

From Eq. (13), it should be clear why we need not worry about mode mixing.
Even though we are using $j=\\{\ell,|m|,n\\}$ as a label, the polarization
functions $h_{j}^{(+/\times)}$ know nothing about the meaning of those
indices: that information has been absorbed by $A_{j}$, $\epsilon_{j}$,
$\theta_{j}$ and $\phi_{j}$, which for us will be nuisance parameters [see Eq.
(II.1)]. In fact, the functions

$h^{(+/\times)}_{j}\equiv
h_{+/\times}(t;\,\omega_{j},\tau_{j},A_{j},\epsilon_{j},\theta_{j},\phi_{j})$
(14)

also ignore that $\omega_{j}$ and $\tau_{j}$ originate in BH perturbation
theory: we can construct templates
$h_{+/\times}(t;\,\omega,\tau,A,\epsilon,\theta,\phi)$ with any generic set of
arguments. If we wanted to incorporate a physical model for the amplitudes
(e.g., in terms of the progenitor parameters in a binary coalescence, as in
[67]), we could always do so by restoring the spheroidal harmonics in Eq.
(II.1), as we discuss in App. B; in the absence of such a model, the
spheroidal harmonics just introduce redundant parameters to the model, and are
thus detrimental.

In spite of the several assumptions we have used to motivate it, the template
in Eq. (13) is extremely general. It can easily be made to accommodate
retrograde modes by passing the corresponding $(\omega_{[-1]\ell
mn},\,\tau_{[-1]\ell mn})$ as input to Eq. (14). Although this would be
inappropriate when perturbing around the Kerr solution, it can also circumvent
the assumption of PT symmetry by simultaneously including pairs of modes like

$h^{(+/\times)}_{\pm j}=h_{+/\times}(t;\,\omega_{\pm j},\tau_{\pm j},A_{\pm
j},\pm 1,\theta_{\pm j},\phi_{\pm j})\,,$ (15)

defining $\pm j=\\{\ell,\pm|m|,n\\}$, and allowing $\omega_{\pm j}$ and
$\tau_{\pm j}$ to take independent values unrestricted by the PT condition,
$\tilde{\omega}_{+j}=-\tilde{\omega}^{*}_{-j}$. In that general case, the two
circularly-polarized modes, $+|m|$ and $-|m|$, will not combine into a single
elliptical contribution.

Figure 2:  Frequency $f_{22n}$ (left), damping time $\tau_{22n}$ (middle) and
quality factor $Q_{22n}=\pi f_{22n}\tau_{22n}$ (right) for different
$\ell=m=2$ tones, as a function of of dimensionless BH spin $\chi$. Times are
measured in units of $t_{M}\equiv GM/c^{3}$ for BH mass $M$. The fundamental
mode ($n=0$) is defined to be the mode with the largest $\tau_{22n}$; the
overtones ($n\geq 1$) are ordered with decreasing $\tau_{22n}$ for increasing
$n$. For the fundamental mode (blue), $Q_{220}$ asymptotes to
$Q_{\mathrm{min}}\simeq 2.1$ as $\chi\to 0$, while it diverges to infinity as
$\chi\to 1$ (right, dashed lines).

Under most circumstances, we expect the signal to be dominated by modes with
$\ell=|m|=2$. If so, we may wish to construct a template including $N$
overtones of the $\ell=|m|=2$ angular mode in Eq. (13), i.e.

$s_{I}^{(22N)}=\sum_{n=0}^{N}\left[F_{+}^{I}h^{(+)}_{22n}+F_{\times}^{I}h^{(\times)}_{22n}\right],$
(16)

suppressing time dependence. This is a template with $N+1$ modes, for a total
of $6\left(N+1\right)$ damped-sinusoid parameters in the general case (fewer
if we restrict to a Kerr spectrum, as discussed below). The choice of $N$ is,
in principle, arbitrary. If $\ell>2$ modes are also expected, we may
incorporate them with the same or different number of respective overtones.

###  II.2 Kerr ringdowns

In the case of a Kerr BH, all QNM frequencies and damping rates are determined
fully by the hole’s mass $M$ and dimensionless spin $\chi$. We may thus write
$\tilde{\omega}_{\ell mn}=\tilde{\omega}_{\ell mn}^{\rm(GR)}(M,\chi)$, so that
the set of free parameters reduces to

$\left\\{M,\chi,A_{j},\epsilon_{j},\theta_{j},\phi_{j}\right\\},$ (17)

for a total of $2(2N+3)$ degrees of freedom for an $N$-overtone model
including only $\ell=|m|=2$. Throughout, we use the qnm Python package to
calculate Kerr ringdown frequencies and damping times [73].

The Kerr metric imposes some notable restrictions on the allowed QNM
frequencies $f_{\ell mn}\equiv\omega_{\ell mn}/(2\pi)$ and damping rates
$\gamma_{\ell mn}\equiv 1/\tau_{\ell mn}$. As we mentioned above, the
symmetries of the perturbation equations imply $\tilde{\omega}_{\ell
mn}^{\rm(GR)}=-\tilde{\omega}_{\ell\scalebox{0.55}{$-$}mn}^{\rm(GR)*}$, or,
equivalently,

$f_{\ell mn}^{\rm(GR)}=-f_{\ell\scalebox{0.55}{$-$}mn}^{\rm(GR)}\,,$ (18)
$\tau_{\ell mn}^{\rm(GR)}=\tau_{\ell\scalebox{0.55}{$-$}mn}^{\rm(GR)}\,.$ (19)

Further structure arises when we consider the functional dependence of these
quantities on the BH parameters. For example, Fig. 2 shows the frequency and
damping times for different prograde tones with $\ell=|m|=2$ as a function of
BH spin. The mass $M$ of the Kerr BH acts simply as an overall scale on the
frequency and damping times, so that dimensionless numbers like the product
$f\tau$ or ratio $f_{j}/f_{j^{\prime}}$ are functions of the spin parameter
$\chi$ only.

Although the frequencies and damping rates can be made to take a broad range
of values by varying $M$ and $\chi$, there are some restrictions. For the
fundamental $\ell=|m|=2$ mode ($n=0$), the quality factor (suppressing
indices)

$Q\equiv\pi f\tau$ (20)

is required to be greater than
$Q_{\mathrm{min}}=Q_{220}\left(\chi=0\right)\simeq 2.1$, so it is not possible
to match _every_ underdamped ($Q>1/2$) mode with some combination of $M$ and
$\chi$, nor is it possible to match _any_ critically damped or overdamped mode
(see Fig. 3).

Figure 3:  Quality factor versus spin magnitude for the fundamental
$\ell=|m|=2,\,n=0$ Kerr ringdown mode (cf., Fig. 2). A fractional deviation
$\delta f_{220}^{\prime}$ around the Kerr fundamental frequency results in an
equal perturbation to $Q_{220}$; a model that assumes GR ($\delta
f_{220}^{\prime}=0$) can account for this by inferring a BH spin ($\chi$)
different than the truth ($\chi^{\prime}$). The figure illustrates this for
$\delta f_{220}^{\prime}=0.9$ and true spin $\chi^{\prime}=0.6$, resulting in
an inferred value $\chi\approx 0.92$. The correspondence between $\delta
f_{220}^{\prime}$ and $\chi$ breaks down for negative values of $\delta
f_{220}^{\prime}$ that would yield $Q_{220}<Q_{\mathrm{min}}$ (horizontal
dashed line).

We have found that the spin-dependence of the $f$ and $\tau$ parameters for
low order Kerr modes in GR can be well approximated by a linear combination of
$\log(1-\chi)$, and powers from $\chi^{0}$ to $\chi^{4}$: for each mode the
dimensionless frequency $f\,t_{M}$ satisfies

$f\,t_{M}\simeq c_{l}\log\left(1-\chi\right)+\sum_{i=0}^{4}c_{i}\chi^{i},$
(21)

and similarly for the dimensionless damping rate $\gamma\,t_{M}$. The
approximating coefficients $c$ can be found by, e.g., least-squares fits to
the complex mode frequencies computed by the qnm package [73]. Our ringdown
package does this automatically for any requested $\\{\ell,|m|,n\\}$, but we
tabulate the values for a number of tones of the quadrupolar harmonic in App.
C as an example. Implementing the relation between the Kerr parameters and the
mode frequencies and damping times in this way is advantageous in contexts
where it is required to differentiate mode parameters with respect to Kerr
parameters (e.g., when using Hamiltonian Monte Carlo to sample from a
posterior density defined in the Stan language [57], as we will do below).

###  II.3 Deviations from Kerr

In order to verify that a given ringdown signal is consistent with a Kerr
spectrum, one could imagine applying a fully-general model by which all
frequencies and damping rates are allowed to vary freely in Eq. (14). However,
such a generic model provides no straightforward way of quantifying agreement
with the Kerr hypothesis: it can produce a $(2\times D)$-dimensional posterior
for the frequencies and damping times of $D$ modes, but cannot in itself
evaluate the degree of consistency with combinations allowed by the Kerr
model. Furthermore, parameter degeneracies render such a model impractical if
more than one or two modes are included, presenting challenges even for high
SNRs. It is also useful for practical reasons (specifying priors, sampling
efficiently, etc.) to be able to make direct connections between the observed
mode properties and physically relevant variables in GR, such as masses and
spins.

While a fully-general analysis can match any modes observed in the data, under
most circumstances we expect that the modes that do appear will be at least
close to the predictions for a Kerr BH in GR. Thus, instead of allowing all
frequencies and damping rates to vary freely, we may replace at least two of
those parameters with the Kerr values derived from some $M$ and $\chi$, while
allowing other modes to float around their corresponding Kerr values. In other
words, $M$ and $\chi$ are standing in for $f$ and $\tau$ of one of the
spectroscopic modes, upon which we therefore implicitly impose a minimum
allowed $Q$ (see Fig. 3); in a sufficiently modified theory these parameters
need not have, even approximately, the usual Kerr interpretation. In that way,
the set of free parameters becomes

$\left\\{M,\chi,\delta{f}_{j},\delta{\tau}_{j},A_{j},\epsilon_{j},\theta_{j},\phi_{j}\right\\},$
(22)

where $\delta{f}_{j}$ and $\delta{\tau}_{j}$ are fractional deviations away
from the Kerr prediction, i.e.,

$f_{j}=f_{j}^{\rm(GR)}(M,\chi)\left(1+\delta{f}_{j}\right),$ (23)
$\tau_{j}=\tau_{j}^{\rm(GR)}(M,\chi)\left(1+\delta{\tau}_{j}\right),$ (24)

for all but two of the included $j$’s; for those two designated modes, we set
$\delta f_{j}=\delta\tau_{j}=0$, so that $f_{j}=f_{j}^{\rm(GR)}(M,\chi)$ and
$\tau_{j}=\tau_{j}^{\rm(GR)}(M,\chi)$. Equivalently, we could chose to work
with fractional deviations on the damping rate $\gamma_{j}$, with
$\delta\gamma_{j}\approx-\delta\tau_{j}$ for small deviations; although we do
not do this here, taking $\delta\gamma_{j}$ as the primary quantity would have
the advantage of avoiding technical issues that arise when
$\delta\tau_{j}\to-1$ (see discussion toward the end of Sec. V.2.1). Finally,
another alternative is to work with deviation quantities that exhibit no
singularities for finite parameter values via

$\displaystyle f_{j}$
$\displaystyle=f_{j}^{\rm(GR)}\left(M,\chi\right)\exp\left(\delta
f_{j}\right)$ (25) $\displaystyle\tau_{j}$
$\displaystyle=\tau_{j}^{\rm(GR)}\left(M,\chi\right)\exp\left(\delta\tau_{j}\right).$
(26)

In all cases, the Kerr prediction is recovered when $\delta
f_{j}=\delta\tau_{j}=0$, and all parameterizations agree in the limit $\delta
f_{j},\delta\tau_{j}\ll 1$.

When the modes in question are overtones of a given $(\ell,|m|)$, we will
additionally require that the $\delta\tau_{j}$ values be bounded so as to
preserve tone ordering, i.e. $\tau_{22,n+1}<\tau_{22,n}$ (a similar condition
could be imposed on the frequencies when dealing with modes of different
$|m|$). For an $N$-overtone template, Eq. (16), this model has $6(N+1)$
degrees of freedom.

The fractional deviations $\delta f_{j}$ and $\delta\tau_{j}$ encapsulate all
relevant information about the agreement of a given signal with the Kerr
hypothesis. If the signal is well described by a Kerr spectrum for any valid
combination of $M$ and $\chi$, then we should find the posterior to be
consistent with $\delta f_{j}=\delta\tau_{j}=0$ for all modes under
consideration. Moreover, the characteristic width of the $\delta f_{j}$ or
$\delta\tau_{j}$ posteriors (say, the 90%-credible interval) quantifies the
degree to which we can establish agreement with the Kerr scenario, providing a
natural measure of the test’s precision.

In generic beyond-Kerr models, it is reasonable to expect the $\delta f_{j}$
and $\delta\tau_{j}$ to themselves be functions of the mass and spin, plus any
other parameters intrinsic to the model (e.g., the value of a scalar field).
Although this has little bearing on the analysis of individual signals, it
does complicate the pooling of results across catalogs of sources, which would
demand a hierarchical Bayesian treatment [74, 75]. Alternatively, it is
possible to write the $\delta f_{j}$ and $\delta\tau_{j}$ as a power series
expansion on the BH spin $\chi$, with source-independent coefficients to be
specified by any given theory [76, 77]. Under that framework, the $\delta
f_{j}$ and $\delta\tau_{j}$ would be replaced by the new set of coefficients
(whose number depends on the order of the expansion) as free parameters,
potentially facilitating the combination of observations888This is only true
if we assume that all sources analyzed belong to the same population;
hierarchical treatments would still be needed if we allowed for the
possibility of mixtures (e.g., a population composed of both Kerr BHs, and
some exotic BH mimicker). and the derivation of theoretical implications.
Although we do not adopt such a parameterization here, it would be trivial to
do so, without altering the qualitative nature of our conclusions.

####  II.3.1 Four-parameter models

We now turn to the question of how to best assign the deviation parameters in
Eqs. (23) and (24) to different QNMs. As implied above, the spectroscopic
analysis requires that at least two modes be detected: for a single-mode
measurement, $\delta f$ and $\delta\tau$ would be fully degenerate with the
mass and spin, modulo the $Q\geq Q_{\rm min}$ constraints imposed by the Kerr
spectrum (Sec. II.2). In a model with $D$ modes, we may consider up to $2D-2$
deviation parameters; in the simplest case, $D=2$ and we are left with four
spectroscopically relevant parameters: $M$ and $\chi$, plus a single $\delta
f$ and a single $\delta\tau$. As we show below, even when considering angular
harmonics other than $\ell=|m|=2$, it is generally preferable to assign $M$
and $\chi$ to the best measured mode.999Readers not interested in
parameterizations of the two-mode model may skip the rest of this section with
impunity; we continue the discussion of the formalism in Sec. III.

While keeping our focus on the $\ell=|m|=2$ angular harmonic, the most natural
two-mode model to consider is the one made up of the two longest lived tones,
i.e., the fundamental and the first overtone [39, 44]. In that case, the Kerr-
deviation parameters may be assigned to any two of
$\\{f_{220},\tau_{220},f_{221},\tau_{221}\\}$. Since, we expect the
fundamental mode to be more easily measured, it is reasonable to leave its
frequency and damping time unperturbed, and assign both deviation parameters
to the first overtone. The spectral model would thus become

$\displaystyle f_{220}$ $\displaystyle=f^{\rm(GR)}_{220}\left(M,\chi\right),$
(27a) $\displaystyle\tau_{220}$
$\displaystyle=\tau^{\rm(GR)}_{220}\left(M,\chi\right),$ (27b) $\displaystyle
f_{221}$ $\displaystyle=f^{\rm(GR)}_{221}\left(M,\chi\right)\left(1+\delta
f_{221}\right),$ (27c) $\displaystyle\tau_{221}$
$\displaystyle=\tau^{\rm(GR)}_{221}\left(M,\chi\right)\left(1+\delta\tau_{221}\right).$
(27d)

Whenever the fundamental mode dominates the measurement (which we expect to be
the case essentially always), it will pin down $M$ and $\chi$ in this model,
while the overtone is allowed to explore alternative values around Kerr
solution. Except for the restriction to $Q_{220}>Q_{\mathrm{min}}$ (Fig. 2),
this parameterization is fully general for the two-mode case, in that it can
fit (almost) any pair of QNMs. Equations (27) were put to use in [39] and,
subsequently, [44] to produce observational results.

We will adopt the above as our default parameterization. It is instructive to
see why this choice is a good one by exploring inferior alternatives. For
example, another option for the $N=1$ model would be to allow deviations in
the two frequencies by writing

$\displaystyle f_{220}$
$\displaystyle=f_{220,\mathrm{Kerr}}\left(M^{\prime},\chi^{\prime}\right)\left(1+\delta
f_{220}^{\prime}\right)$ (28a) $\displaystyle\tau_{220}$
$\displaystyle=\tau_{220,\mathrm{Kerr}}\left(M^{\prime},\chi^{\prime}\right)$
(28b) $\displaystyle f_{221}$
$\displaystyle=f_{221,\mathrm{Kerr}}\left(M^{\prime},\chi^{\prime}\right)\left(1+\delta
f_{221}^{\prime}\right)$ (28c) $\displaystyle\tau_{221}$
$\displaystyle=\tau_{221,\mathrm{Kerr}}\left(M^{\prime},\chi^{\prime}\right),$
(28d)

using primes to distinguish this from our default parameterization. While not
strictly degenerate, Eqs. (28) are _practically_ degenerate: although they
appear to involve four degrees of freedom, there are only three in practice.

We can see this by studying the relation between
$\\{M^{\prime},\,\chi^{\prime},\,\delta f_{220}^{\prime},\,\delta
f_{221}^{\prime}\\}$ and $\\{M,\,\chi,\,\delta f_{221},\,\delta\tau_{221}\\}$
through Eqs. (28) and Eqs. (27). To obtain a mapping between the two sets of
coordinates, first note that the quality factor $Q_{220}=\pi
f_{220}\tau_{220}$ can be written in both parameterizations as

$\displaystyle Q_{220}$
$\displaystyle=f_{220}^{\rm(GR)}\left(M,\chi\right)\tau_{220}^{\rm(GR)}\left(M,\chi\right)$
(29)
$\displaystyle=f_{220}^{\rm(GR)}\left(M^{\prime},\chi^{\prime}\right)\tau_{220}^{\rm(GR)}\left(M^{\prime},\chi^{\prime}\right)\left(1+\delta
f_{220}^{\prime}\right).$

Since $f_{220}^{\rm(GR)}\tau_{220}^{\rm(GR)}$ is independent of $M$, we can
numerically solve the above expression for $\chi$ as a function of the primed
quantities. In other words, as long as $\delta f_{220}^{\prime}$ is not so
small as to render $Q_{220}<Q_{\mathrm{min}}$, we can always find a value
$\chi$ that yields the same quality factor as a different value
$\chi^{\prime}$ plus the deviation $\delta f_{220}^{\prime}$ (Fig. 3). The
restriction to $Q_{220}>Q_{\mathrm{min}}$ is not worrisome because, since we
can assume that the $n=0$ mode dominates the signal, a ringdown with
$Q_{220}<Q_{\mathrm{min}}$ would be easily identified as anomalous by regular
parameter estimation of the full IMR signal within GR.

Figure 4:  Ratio of the first-overtone frequency (top) and damping time
(bottom) to the fundamental mode values, as a function of BH spin $\chi$,
assuming a Kerr BH. The labels on the right of each plot present the ratios as
fractional deviations of the value at $\chi=0.5$. The ratios only vary by a
few percent over the full range of $\chi$. (See also Fig. 2.) Figure 5:
Comparison between the distributions of $f_{220}$,
$\gamma_{220}=1/\tau_{220}$, $f_{221}$, and $\gamma_{221}=1/\tau_{221}$
induced by flat priors on $1/2<M,M^{\prime}<2$ (arbitrary units),
$0<\chi,\chi^{\prime}<1$, $-0.9<\delta f_{220}^{\prime}<0.9$, $-0.9<\delta
f_{221},\delta f_{221}^{\prime}<0.9$, and $-0.9<\delta\tau_{221}<0.9$ in the
Isi  _et al._ [39] parameterization of Eqs. (27) (blue ) and the
parameterization of Eqs. (28) (orange). The degeneracy in the latter
parameterization is apparent in that it produces $\gamma_{220}$–$\gamma_{221}$
pairs that lie (nearly) on a single plane, while the Isi  _et al._ [39] one
explores a finite volume in the four-dimensional frequency-damping rate space
(highlighted by red arrow).

We can proceed similarly to get an expression for $M$ as function of primed
quantities by solving

$f_{220}=f_{220}^{\mathrm{(GR)}}\left(M,\chi\right)=f_{220}^{\mathrm{(GR)}}\left(M^{\prime},\chi^{\prime}\right)\left(1+\delta
f_{220}^{\prime}\right)$ (30)

for $M$. With $M$ and $\chi$ in hand, we can find $\delta f_{221}$ from

$\displaystyle f_{221}$
$\displaystyle=f_{221}^{\mathrm{(GR)}}\left(M^{\prime},\chi^{\prime}\right)\left(1+\delta
f_{221}^{\prime}\right)$ (31)
$\displaystyle=f_{221}^{\mathrm{(GR)}}\left(M,\chi\right)\left(1+\delta
f_{221}\right),$

and, finally, $\delta\tau_{221}$ from

$\tau_{221}=\tau_{221}^{\mathrm{(GR)}}\left(M^{\prime},\chi^{\prime}\right)=\tau_{221}^{\mathrm{(GR)}}\left(M,\chi\right)\left(1+\delta\tau_{221}\right).$
(32)

The above equations suggest that, excepting the $Q_{220}\geq Q_{\rm min}$
restriction, there is a one-to-one mapping between the primed and unprimed
coordinates. However, the structure of the Kerr spectrum implies that Eq. (32)
can only produce a small range $\delta\tau_{221}$ values, and so the
coordinate transformation is effectively degenerate. This becomes apparent by
considering the fundamental-to-overtone damping time ratio as written in each
parameterization,

$\frac{\tau_{220}}{\tau_{221}}=\frac{\tau_{220}^{\mathrm{(GR)}}\left(M^{\prime},\chi^{\prime}\right)}{\tau_{221}^{\mathrm{(GR)}}\left(M^{\prime},\chi^{\prime}\right)}=\frac{\tau_{220}^{\mathrm{(GR)}}\left(M,\chi\right)}{\tau_{221}^{\mathrm{(GR)}}\left(M,\chi\right)}\frac{1}{1+\delta\tau_{221}}.$
(33)

Figure 4 shows this dimensionless ratio as a function of the spin parameter;
notably, the ratio varies by only $\sim\pm 1.2\%$ across the entire allowed
range of $0\leq\chi<1$. Thus, irrespective of the values of $\chi$ and
$\chi^{\prime}$, Eq. (33) requires $\left(1+\delta\tau_{221}\right)\simeq 1$,
or $\delta\tau_{221}\simeq 0$. Therefore, varying $M^{\prime}$,
$\chi^{\prime}$, $\delta f_{220}^{\prime}$, and $\delta f_{221}^{\prime}$ is
nearly identical to varying $M$, $\chi$, and $\delta f_{221}$ with
$\delta\tau_{221}=0$, and the seemingly four-parameter model of Eq. (28) all
but reduces to the three-parameter model $\\{M,\,\chi,\ \delta f_{221}\\}$.
The degeneracy can only be broken by a $\tau_{221}$ measurement with better
than ${\sim}1\%$ precision.

Figure 5 demonstrates this near-degeneracy empirically. We show the joint
distribution of $f_{220}$, $\gamma_{220}=1/\tau_{220}$, $f_{221}$, and
$\gamma_{221}=1/\tau_{221}$ induced by a uniform prior on $1/2<M,M^{\prime}<2$
(arbitrary units), $0<\chi,\chi^{\prime}<1$, $-0.9<\delta
f_{220}^{\prime}<0.9$, $-0.9<\delta f_{221},\delta f_{221}^{\prime}<0.9$, and
$-0.9<\delta\tau_{221}<0.9$ in the parameterization of Eq. (27) (un-primed)
and the parameterization of Eqs. (28) (primed). Using this latter
parameterization, the prior only permits ratios of $\gamma$ values that lie on
a (nearly) degenerate plane of (nearly) constant $\gamma_{220}/\gamma_{221}$;
using the former parameterization, the prior allows frequencies and damping
rates that explore a finite volume in the four-dimensional space. (See also
Sec. V.2.3.)

As for the damping times, Fig. 4 also shows the ratio of Kerr ringdown
frequencies is a small $\sim\pm 3.5\%$ over $0\leq\chi<1$; this suggests that
a parameterization with $\delta\tau_{220}$ and $\delta\tau_{221}$ would suffer
from a near-degeneracy analogous to Eq. (28).

####  II.3.2 Three-parameter models

The first overtone’s damping time, $\tau_{221}$, is often quite poorly-
measured [39]. Under the assumption that it is effectively unconstrained by
the data, the two-mode problem reduces to a three-dimensional parameter space.
This is equivalent to $\delta\tau_{221}$ in Eq. (27d) being unconstrained.
Therefore, the fully-general (up to the $Q_{220}>Q_{\mathrm{min}}$ constraint)
three-parameter likelihood can be parameterized by Eqs. (27) with
$\delta\tau_{221}=0$.

As above, it is useful to consider an alternative parameterization like

$\displaystyle f_{220}$
$\displaystyle=f_{220}^{\mathrm{(GR)}}\left(M^{\prime},\chi^{\prime}\right)\left(1+\delta
f_{220}^{\prime}\right)$ (34a) $\displaystyle\tau_{220}$
$\displaystyle=\tau_{220}^{\mathrm{(GR)}}\left(M^{\prime},\chi^{\prime}\right)$
(34b) $\displaystyle f_{221}$
$\displaystyle=f_{221}^{\mathrm{(GR)}}\left(M^{\prime},\chi^{\prime}\right),$
(34c)

where again primes denote quantities in the alternative parameterization. By
the same logic we applied to Eq. (29), it is clear that we can transform
between primed and unprimed coordinates in this case also, as long as $\delta
f_{220}^{\prime}$ does not push $Q_{220}$ below the Kerr constraint. Figure 6
shows the joint distribution of $\delta f_{220}^{\prime}$ and $\delta f_{221}$
under this one-to-one mapping when $1/2<M^{\prime}<2$ (arbitrary units),
$0<\chi^{\prime}<1$ and $-0.9<\delta f_{220}^{\prime}<0.9$ are drawn uniformly
within their respective ranges. In spite of the wide distribution of
$M^{\prime}$ and $\chi^{\prime}$ values, there is a very narrow one-to-one
correspondence between $\delta f_{220}^{\prime}$ and $\delta f_{221}$
(insensitive to the arbitrary $M$ range chosen).

Figure 6:  Joint distribution of $\delta f_{220}^{\prime}$ and $\delta
f_{221}$ under a uniform prior on $0.5<M^{\prime}<2$, $0<\chi^{\prime}<1$, and
$-0.9<\delta f_{220}^{\prime}<0.9$ under the one-to-one mapping between the
reduced three-dimensional parameterizations in the text. In spite of the range
of $M^{\prime}$ and $\chi^{\prime}$ values permitted at fixed $\delta
f_{220}^{\prime}$, the relation between $\delta f_{220}^{\prime}$ and $\delta
f_{221}$ is exceedingly tight.

The points in the $\delta f_{220}^{\prime}-\delta f_{221}$ plane closely
follow a hyperbola. We can see that this is the case by considering the ratio
of frequencies

$\displaystyle\frac{f_{220}}{f_{221}}$ $\displaystyle=\frac{f_{\rm
220,Kerr}(M^{\prime},\chi^{\prime})}{f_{\rm
221,Kerr}(M^{\prime},\chi^{\prime})}\left(1+\delta f_{220}\right)$
$\displaystyle=\frac{f_{\rm 220,Kerr}(M,\chi)}{f_{\rm
221,Kerr}(M,\chi)}\frac{1}{\left(1+\delta f_{221}\right)}\,.$ (35)

The ratios of Kerr frequencies are independent of mass and only weakly
dependent on the spin (Fig. 4), so they are approximately equal no matter the
values of $\chi$ and $\chi^{\prime}$. To good accuracy, this implies

$\left(1+\delta f_{220}^{\prime}\right)\left(1+\delta
f_{221}\right)\approx\text{constant}\,,$ (36)

which describes a hyperbola that asymptotes to $\delta f_{220}^{\prime}=-1$ as
$\delta f_{221}\rightarrow\infty$, and vice versa. Therefore, if we restrict
$-1<\delta f_{220}^{\prime}<1$, we will necessarily have $-0.5<\delta
f_{221}<\infty$, with the upper (lower) limit on $\delta f_{220}^{\prime}$
setting the lower (upper) limit on $\delta f_{221}$. This explains the range
of the $\delta f_{221}$ marginal distribution in Fig. 6 (lower right). The
full shape of that curve can be well approximated analytically by applying the
Jacobian implied by Eq. (36) to the $\delta f_{220}^{\prime}$ distribution.

Since there is a tight, nearly one-to-one mapping between $\delta
f_{220}^{\prime}$ and $\delta f_{221}$ even over a broad range of $M^{\prime}$
and $\chi^{\prime}$, using both parameterizations to fit real ringdown data
would not return independent measurements unless the overtone damping time can
be constrained to sub-percent precision.

##  III Inference framework

A signal template, such as discussed above, is only one of the many moving
pieces that make up the Bayesian infrastructure required to extract
information from ringdown data. In this section, we motivate and describe our
time domain domain framework for Bayesian inference with ringdown signals
(Sec. III.1), and show it is mathematically equivalent to a recently proposed
alternative (Sec. III.2); we also provide implementation details (Sec. III.3).

The reason for spending time carefully examining data analysis techniques is
that the framework standard to LIGO-Virgo transient analyses is not readily
suitable for ringdown-only studies. The overwhelming majority of such analyses
are constructed in the Fourier domain, the fundamental reason for this being
that nominal instrumental noise is well described as a stationary Gaussian
process. Assuming periodic boundary conditions, this means that the noise
covariance matrix diagonalizes in the frequency domain. The noise Fourier
amplitudes become independent random variables with variance described by some
one-sided power spectral density (PSD) as a function of frequency, $S(f)$ [78,
79]. This vastly simplifies Bayesian computations because, in that case, the
log-likelihood for a signal $\tilde{s}(f)$ is given by a simple noise-weighted
inner product,

$\ln P(\tilde{d}\mid\tilde{s})=-2{\Delta
f}\sum_{i=0}^{N-1}\frac{|\tilde{d}_{i}-\tilde{s}_{i}|^{2}}{S(f_{i})}+\mathrm{const.},$
(37)

for Fourier-domain data $\tilde{d}$ sampled at $N$ frequencies $f_{i}\geq 0$,
with $\Delta f=f_{i+1}-f_{i}$. Evaluating Eq. (37) takes $O(N)$
computations—an improvement over the $O(N^{3})$ required for an arbitrary
covariance matrix. The properties of the signal are secondary, as long as it
is smooth, contained within the frequency band, and sufficiently short for the
stationarity and periodic-boundary assumptions to apply.

Figure 7: The challenge of ringdown-only analyses. Targeting the ringdown
portion of an IMR signal requires discarding the data before some truncation
time $t_{0}$, here the peak of the signal (dotted gray line). A regular
Fourier domain analysis requires smoothly windowing the data segment—not doing
so is equivalent to applying a sharp rectangular window (top) that will lead
to spectral leakage. Avoiding spectral leakage requires a smooth roll-on
(e.g., a Tukey window), which is not compatible with a signal that has
significant content at the edge of the segment: a window that offers support
for $t<t_{0}$ will allow contamination from non-ringdown data (middle), while
one that does not will necessarily corrupt the start of the ringdown signal
(bottom). The issue can be circumvented by remaining in the time domain (Sec.
III.1), as we propose, or by applying nontrivial modifications to the
frequency-domain likelihood (Sec. III.2), as was proposed in [47, 80].

Unfortunately, analyses aimed at isolating the ringdown fail to meet these
basic requirements: it is not possible to enforce a cyclic boundary, and avoid
spectral leakage, without corrupting part of the signal. Regular CBC analyses
based on Eq. (37) rely on stretches of time-domain detector data chosen such
that the signal lies a safe distance away from the segment edges, which are
tapered smoothly to prevent leakage when Fourier transforming [81]. However,
that is not an option when targeting the ringdown in isolation: the signal in
which we are interested is necessarily located at the edge of the data
segment, and would be corrupted by any windowing procedure (Fig. 7). Retaining
data before the ringdown is not viable because the whitening filter for
colored noise will necessarily couple points across the cutoff, contaminating
the targeted signal [31]. Simply Fourier transforming an unwindowed data
segment starting at the beginning of the ringdown will both cause spectral
leakage and correlate the beginning of the segment to its end.

In spite of such fundamental obstacles, there have been several attempts to
formulate ringdown analyses fully or partially in the traditional Fourier
domain framework [82, 83, 31, 43]. The original LIGO-Virgo study in [82]
followed a mixed strategy, which consisted of inverse-Fourier transforming
overwhitened detector data to compute the likelihood in the time domain [83].
Although this method avoids explicitly Fourier-transforming a discontinuous
template, it is mathematically equivalent to Eq. (37): the overwhitening
filter couples data before and after the truncation point, contaminating the
ringdown measurement as illustrated in Fig. 10. The approach was successful
for GW150914 in [82] because the chosen truncation times lay far away from the
signal peak: with low SNR immediately before the truncation point, the
systematics caused by overwhitening were likely smaller than the statistical
noise; however, this would fail for times closer to the peak or for louder
signals. Other proposed strategies, like zeroing out data before the cutoff
[31], or replacing it with detector noise [43], fail to prevent the whitening
filter from coupling data across the boundary—which occurs irrespective of the
nature of the pre-cutoff data (whether the data before the cutoff are zeros or
noise is immaterial)—and can also be subject to spectral leakage induced by
sharp features in the data. None of these problems arise if we formulate our
analysis in the time domain without enforcing cyclic boundary conditions.

###  III.1 Time domain likelihood

symmetric
noncirculant
Toeplitz matrix

$\begin{pmatrix}\rho_{0}&\rho_{1}&\dots&\rho_{N-1}\\\
\rho_{1}&\rho_{0}&\dots&\rho_{N-2}\\\ \vdots&&\ddots&\vdots\\\
\rho_{N-1}&\dots&\rho_{1}&\rho_{0}\\\ \end{pmatrix}$ (38)

Figure 8: Symmetric Toeplitz matrices. _Left:_ illustration of the symmetric-
Toeplitz matrix structure of Eq. (42) in both circulant (main matrix) and
noncirculant (highlighted submatrix) forms; the shade of each square
represents the value of a matrix entry $C_{ij}$: the circulant case satisfies
$C_{ij}=C_{(N-i)j}$ for any entry $0\leq i,j<N$, while the noncirculant case
does not. _Right:_ schematic form of a noncirculant Toeplitz matrix, i.e., a
matrix that satisfies $C_{ij}=C_{i+1,j+1}$ and $C_{ij}=C_{ji}$, or
equivalently $C_{ij}=\rho(|i-j|)$ as in Eq. (42), but not $C_{ij}=C_{(N-i)j}$,
just as in the highlighted submatrix on the left.

We model instrumental noise as a discrete-time random process, represented by
a set of random variables $n_{i}\equiv n(t_{i})$ when sampled at arbitrary
times $t_{i}$. Assuming Gaussianity, the process will be fully characterized
by its mean $\mathrm{E}\left[n_{k}\right]$ and covariance matrix
$C_{ij}\equiv\mathrm{E}\left[n_{i}n_{j}-\mathrm{E}\left[n_{k}\right]\right],$
where $\mathrm{E}$ denotes expectation values. In general, the mean can be
trivially enforced to be zero $\mathrm{E}\left[n_{k}\right]=0$, while the only
restriction on $C$ is that it be positive semidefinite. For an $N$-vector
$n=\\{n_{0},\dots,n_{N-1}\\}$ drawn from such a process, $C$ will be an
$N\times N$ matrix,101010In this section, and only in this section, $N$ stands
for the number of samples analyzed, not number of overtones in a template. and
the log-probability of the draw will be given by

$\ln P(n)=-\frac{1}{2}\sum_{i,j}n_{i}\,C_{ij}^{-1}\,n_{j}+\mathrm{const.},$
(39)

where $C^{-1}$ is the inverse covariance matrix. Equation (39) requires
$O(N^{2})$ computations if the covariance matrix is pre-processed into a
convenient form. Such preprocessing costs $O(N^{3})$.

To reduce the computing cost per likelihood evaluation, we may for example
Cholesky-decompose the covariance matrix into a lower-triangular factor $L$
and its transpose, such that

$C_{ij}=\sum_{k}L_{ik}L_{kj}\,.$ (40)

The inverse of $L$ can be computed by backsubstitution efficiently, and can
act on an arbitrary time series $x$ to produce a new, uncorrelated
(“whitened”) time series,

$\bar{x}_{i}=\sum_{j}L_{ij}^{-1}x_{j}\,.$ (41)

If the covariance of $n_{i}$ is given by Eq. (40), then each component of
$\bar{n}$ will be independently drawn from a unit normal, i.e.,
$\bar{n}_{i}\sim\mathcal{N}(0,1)$. With that definition, the likelihood in Eq.
(39) simplifies to $\ln P(n\mid
C)=-\frac{1}{2}\sum_{i,j}\bar{n}_{i}\,\bar{n}_{j}+\mathrm{const.}$

Figure 9: Obtaining an aperiodic ACF: starting from a long stretch of time-
domain noise with duration $N_{\rho}\Delta t$ (left), we can obtain a PSD
through usual techniques like Welch estimation with a coherent integration
length of $N_{\rm FFT}$ (middle); we can then inverse-Fourier transform the
PSD to obtain a long, periodic estimate of the noise ACF (right, main panel),
which we can finally truncate to the desired analysis length $N$ (right,
inset) to form the covariance matrix through Eq. (42). Alternatively, we could
estimate the ACF directly by autocorrelating the time domain data, or we could
start the process from a PSD obtained by a different method [84, 80].

For a stationary random processes, the covariance takes a particularly simple
(symmetric Toeplitz) form,

$C_{ij}=\rho\left(\left|i-j\right|\right),$ (42)

where $\rho(k)$ is the autocovariance function (ACF). This can be estimated
empirically by autocorrelating a long stretch of noise-only data (i.e., a
segment of length $N_{\rho}\gg N$),

$\hat{\rho}(k)=\frac{1}{N_{\rho}}\sum_{i=0}^{N_{\rho}-1}n_{i}n_{i+k}\,,$ (43)

for $|k|<N_{\rho}$, and zero otherwise111111This is a biased estimator of the
true ACF; the bias at lag $k$ is a factor of $N_{\rho}/(N_{\rho}-|k|)$. We
favor Eq. (43), however, because it reduces the influence of large lag terms
where the variance in the ACF estimate is larger; it is thus more stable than
the unbiased estimator..

If, in addition to stationarity, we impose periodic boundary conditions, then
$\rho(k)=\rho(N-k)$ and $C$ will be circulant (Fig. 8). Circulant matrices are
diagonalized by the discrete Fourier transform [79], meaning the noise Fourier
amplitudes $\tilde{n}_{i}$ will be drawn from a Gaussian process with
covariance matrix

$\tilde{C}_{ij}\equiv
E[\tilde{n}_{i}\tilde{n}_{j}]=\frac{1}{2}\,T\,S(|f_{i}|)\,\delta_{ij}\,,$ (44)

where $T=N\Delta t$, $S(|f|)$ is the one-sided PSD, and $\delta_{ij}$ is the
Kronecker delta. The PSD is derived from the cyclic ACF through a discrete
Fourier transform,

$S(|f_{j}|)=2\Delta t\sum_{k=0}^{N-1}\rho(k)\,e^{-2\pi ijk/N}\,,$ (45)

for $f_{j}=j\Delta f=j/T$, and $0\leq j<N$. Naturally, because $\rho(k)$ is
real-valued, $S(f_{j})=S(f_{N-j})$.

In the limit of infinite observation time ($N\to\infty$), Eq. (44) is a
consequence of the Wiener-Khinchin theorem and the diagonalization it implies
is exact, no matter the specific structure of $\rho(k)$; for large but finite
$N$, the diagonalization is only approximate, unless $C$ is exactly circulant
[79, 85]. In either case, we can always choose to work with Fourier-domain
quantities without loss of generality. However, $\tilde{C}$ will only be
diagonal if and only if $C$ is circulant. If that is not the case, Eq. (37)
fails to apply and there is little incentive to switch to the Fourier
domain.121212There will always exist a coordinate transformation (known as
Karhunen-Loéve transformation) to diagonalize any covariance matrix [79], but
there will generally not exist an efficient algorithm (like the fast Fourier
transform) to effect it.

Conversely, the Fourier domain likelihood of Eq. (37) is predicated on Eq.
(44) and, therefore, assumes that the time-domain data are periodic. Starting
from a given PSD estimate $\hat{S}(|f|)$ sampled at $N$ frequencies $f_{j}$,
we may invert Eq. (45) to obtain the corresponding ACF estimate,

$\hat{\rho}_{S}(k)=\frac{1}{2T}\sum_{j=0}^{N-1}\hat{S}(|f_{j}|)\,e^{2\pi
ijk/N}\,.$ (46)

This estimator is normalized in the same way as Eq. (43), but will be cyclic
by construction. This additional symmetry is not implied by Eq. (42), or Eq.
(43), but is rather induced by Eq. (44). With $C_{ij}$ derived from
$\hat{\rho}_{S}(k)$, the time-domain expression in Eq. (39) is formally
equivalent to Eq. (37), with $n_{i}=d_{i}-s_{i}$.

The key to our approach is to do without the covariance structure imposed by
Eq. (44), enabling us to analyze a short segment of ringdown data without
corrupting it. Without tapering, our analysis is also free from spurious
Fourier bin covariances [86]. We can achieve this simply by using the time-
domain likelihood of Eq. (39) with a covariance matrix constructed from an
acyclic estimate of the ACF. This can be obtained directly from noise samples
in the time domain via Eq. (43). It can also be derived from a preexisting PSD
estimate using Eq. (46), as long as the PSD was estimated from data segments
of length $N_{\rm FFT}$ much longer than the analysis segment length $N$, with
$N_{\rm FFT}$ replacing $N$ in Eq. (46). A long segment allows us to truncate
the resulting $N_{\rm FFT}$-long $\hat{\rho}_{S}(k)$ to length $N$ before
constructing the covariance matrix (Fig. 9). After truncation,
$\hat{\rho}_{S}(k)\neq\hat{\rho}_{S}(N-k)$, thus breaking circularity.

Once we have constructed an acyclic covariance matrix, we can analyze a noisy
data stream $d(t)$ with a model defined only after some truncation time
$t_{\rm start}$. Assuming the data were presampled at some set of times
$t_{i}$, we will index them such that $i=0$ corresponds to the first sample at
or after the specified truncation time,

$t_{0}\equiv\min(\left\\{t_{i}\text{ such that }t_{\rm start}\leq
t_{i}\\}\right).$ (47)

With this convention, the log-likelihood for data after $t_{\rm start}$
containing a signal $s$ is nothing but

$\ln P\left(d\mid
s\right)=-\frac{1}{2}\sum_{i,j=0}^{N-1}\left(d_{i}-s_{i}\right)C^{-1}_{ij}\left(d_{j}-s_{j}\right),$
(48)

up to a constant, as implied by Eq. (39) and $d_{i}=s_{i}+n_{i}$ for $i\geq
0$. This likelihood is completely agnostic about times before $t_{0}$, and
does not impose a periodic boundary.

Note that it is essential to discard the data preceding the truncation point
before evaluating Eq. (48): it is not sufficient to set $s(t)=0$ for
$t<t_{0}$, as was done in [37]. Such a model predicts no signal until $t_{0}$,
at which point the template turns on sharply—this is decidedly not agnostic
about data for $t<t_{0}$. If one adopts such a model, whether the covariance
is cyclic or not becomes irrelevant: the beginning of the ringdown will be
unavoidably corrupted by the zeroes that preceed it, resulting in shortcomings
similar to those of the Fourier-domain strategy in [31].

The likelihood of Eq. (48) is the key element required to obtain a Bayesian
posterior on the parameters $\theta$ of a potential signal in the data,
$p(\theta\mid d)$. By Bayes’ theorem,

$p(\theta\mid d)\propto p(d\mid s_{\theta})\,p(\theta)\,,$ (49)

where $s_{\theta}\equiv s(t;\theta)$ is the signal corresponding to parameters
$\theta$, $p(d\mid s_{\theta})$ is the likelihood given by Eq. (48), and
$p(\theta)$ is the prior. If multiple, independent data streams are available
(e.g., from different GW detectors), the posterior and likelihood generalize
trivially, such that

$p(\theta\mid\\{d_{I}\\})\propto p(\theta)\prod p(d_{I}\mid
s_{\theta}^{I})\,,$ (50)

where $d_{I}$ represents each independent data stream, $s^{I}_{\theta}$ is the
signal expected at the $I$th stream given parameters $\theta$, and the product
is over all values of $I$.

Just as with the likelihood, it is convenient to have a notion of SNR that is
able to account only for data from some truncation point onward. Such a notion
is useful to quantify the amount of signal power in the ringdown for a given
signal in the data, and to make projections about ringdown detectability.

Equation (39) provides a natural definition for an inner product incorporating
data only after $t_{\rm start}$, for any two time series $x_{i}=x(t_{i})$ and
$y_{i}=y(t_{i})$:

$\langle x\mid
y\rangle_{t_{0}}\equiv\sum_{i,j=0}^{N-1}x_{i}\,C^{-1}_{ij}\,y_{j}=\sum_{i,j=0}^{N-1}\bar{x}_{i}\bar{y}_{j}\,,$
(51)

where we extend the sum $N-1$ steps up to our last data point at time
$t_{N-1}=t_{0}+(N-1)\Delta t=t_{0}+T-\Delta t$. We may use the notion of
distance induced by Eq. (51) to define the _time-bounded optimal SNR_ ,
$\mathrm{SNR}[t_{0}]$, of an arbitrary timeseries, $x_{i}$, as the norm of the
timeseries, namely

$\mathrm{SNR}[t_{0}]\equiv||x||_{t_{0}}\equiv\left\langle x\mid
x\right\rangle_{t_{0}}^{1/2}\,.$ (52)

By the same token, the _time-bounded matched-filter SNR_ , $\mathrm{SNR}_{\rm
mf}[t_{0}]$, of a signal, $h_{i}$, in some noisy data, $d_{i}=s_{i}+n_{i}$, is

$\mathrm{SNR}_{\rm mf}[t_{0}]\equiv{\left\langle s\mid
d\right\rangle_{t_{0}}}{\left\langle s\mid s\right\rangle^{-1/2}_{t_{0}}}\,.$
(53)

We will usually drop simplify notation by dropping the “$[t_{0}]$” specifier
when the chosen truncation time is clear from context. With these definitions
in place, Eq. (48) can be written succinctly as: $\ln P(d\mid
s)=-\frac{1}{2}||d-s||^{2}_{t_{0}}+\text{const.}$

For multiple detectors, Eq. (50) implies that the expressions for the SNR
generalize as usual to the network SNR,

$\mathrm{SNR}_{\rm
net}[t_{0}]\equiv\sqrt{\sum\left(\mathrm{SNR}_{I}[t_{0}+\delta
t_{I}]\right)^{2}}\,,$ (54)

where the sum is over each instrument $I$, and $\mathrm{SNR}_{I}[t_{0}+\delta
t_{I}]$ is the SNR at the $I$th instrument (matched filter or optimal),
evaluated at respective truncation times shifted by $\delta t_{I}$ with
respect to $t_{0}$. This global reference time can be chosen to be the
truncation time at one of the detectors (for which $\delta t_{I}$ would vanish
by definition), or some arbitrary reference like the geocenter.

If $C$ is derived from a circular ACF as in Eq. (45), then Eq. (51) is
equivalent to the usual Fourier-domain inner product implied by Eq. (37),
which weights data by the PSD (see, e.g., [87]). Accordingly, Eqs. (52) and
(53) reduce to the usual quantities in that case, assuming the analysis
segment is extended to encompass the full signal.

###  III.2 Modified frequency-domain likelihood

An alternative to our approach above is to circumvent truncation issues while
remaining in the frequency domain via nontrivial modifications to the
likelihood as proposed in [47, 80]. This approach destroys the diagonality of
the Fourier covariance matrix, thus reducing the incentive for switching to
the frequency domain in the first place, and generally increases the
computational cost. However, it is a valid strategy that is formally
equivalent to ours, so we review it here for completeness.

Consider a long stretch of data with enough samples before $t_{0}$ to permit
windowing so that the complete stretch can be treated as periodic, with a
cyclic, stationary noise covariance matrix $C$; let $M$ be the total length of
this segment, with $N$ samples corresponding to $t\geq t_{0}$ preceded by
$I\equiv M-N$ samples with $t<t_{0}$. The frequency domain noise covariance
for the $M$-long segment, $\tilde{C}$, is diagonal as in Eq. (44), but the
usual likelihood function of Eq. (37) is not appropriate because it depends on
data and signal values for $t<t_{0}$, as detailed in the introduction to this
section. However, we can eliminate this dependence by modifying the time-
domain noise covariance matrix, adding large terms along the diagonal for
entries $i<I$ corresponding to times before the desired analysis start time:

$C\to C+a\sum_{i<I}e_{i}e_{i}^{T},$ (55)

where $e_{i}$ is the $i$th time-domain basis vector, with zeros in every entry
except $i$ (i.e., with the $k$th component given by $\delta_{ik}$), and we
will take $a\to\infty$ at the end of our calculation to render the modified
likelihood completely insensitive to data and signal entries $i<I$. The effect
of this manipulation is to artificially impose infinite uncertainty about the
noise contribution to the data before $t_{0}$.

With the modified covariance matrix, the likelihood becomes

$\displaystyle\ln P\left(d\mid s\right)$
$\displaystyle=-\frac{1}{2}\left(d-s\right)^{T}\left(C+a\sum_{i<I}e_{i}e_{i}^{T}\right)^{-1}\left(d-s\right)$
$\displaystyle=-\frac{1}{2}\left(\tilde{d}-\tilde{s}\right)^{\dagger}\left(\tilde{C}+a\sum_{i<I}\tilde{e}_{i}\tilde{e}_{i}^{\dagger}\right)^{-1}\left(\tilde{d}-\tilde{s}\right),$
(56)

where $\tilde{e}_{i}$ is the Fourier transform of $e_{i}$ (not the $i$th basis
vector in the Fourier domain), $T$ denotes a transpose and $\dagger$ a
Hermitian conjugate. The additional terms in the frequency domain covariance
are not diagonal, so the modified likelihood cannot be computed in $O(M)$ time
in the frequency domain.

Consider the singular value decomposition (SVD) of the extra terms:

$a\sum_{i<I}\tilde{e}_{i}\tilde{e}_{i}^{\dagger}=U\Sigma V^{\dagger}$ (57)

where $\Sigma$ is diagonal $I\times I$ and $U$ and $V$ are unitary $M\times
I$. We have $\Sigma=O(a)$ so $\Sigma^{-1}\to 0$ as $a\to\infty$. The Woodbury
matrix identity states that

$\left(\tilde{C}+a\sum_{i<I}\tilde{e}_{i}\tilde{e}_{i}^{\dagger}\right)^{-1}\\\
=\tilde{C}^{-1}-\tilde{C}^{-1}U\left(\Sigma^{-1}+V^{\dagger}\tilde{C}^{-1}U\right)^{-1}V^{\dagger}\tilde{C}^{-1}$
(58)

Taking $a\to\infty$ introduces infinite uncertainty about the noise
contribution to times $t_{i}$ with $i<I$, eliminating any dependence in the
likelihood on the data or signal at these times. The effect is to eliminate
the $\Sigma^{-1}$ term:

$\lim_{a\to\infty}\left(\tilde{C}+a\sum_{i<I}\tilde{e}_{i}\tilde{e}_{i}^{\dagger}\right)^{-1}\\\
=\tilde{C}^{-1}-\tilde{C}^{-1}U\left(V^{\dagger}\tilde{C}^{-1}U\right)^{-1}V^{\dagger}\tilde{C}^{-1}.$
(59)

The modified likelihood thus becomes

$\ln P\left(d\mid
s\right)=-\frac{1}{2}\left(\tilde{d}-\tilde{s}\right)^{\dagger}\tilde{C}^{-1}F\left(\tilde{d}-\tilde{s}\right)$
(60)

where $F$ is a projection matrix ($F^{2}=F$) defined by [80]

$F=1-U\left(V^{\dagger}\tilde{C}^{-1}U\right)^{-1}V^{\dagger}\tilde{C}^{-1}.$
(61)

The effect of $F$ is to project out the parts of the data and signal that
would render the likelihood sensitive to samples before the desired start
time. If used with a covariance matrix $C_{ij}$ that agrees with the matrix
used in Eq. (48) for $i,j\geq I$, the likelihood functions will be identical.

The likelihood function in Eq. (60) uses the noise covariance in the frequency
domain $\tilde{C}$, much like standard parameter estimation analyses [50],
except the matrix is no longer diagonal. However, the frequency-domain method
requires that the data be circular; in practice, this means that there must be
enough padding before the start of the “analysis segment” and after the
effective end of the ringdown to implement a tapering of the data to enforce
the circularity assumption that renders $\tilde{C}$ diagonal. This means that
$I$ must be fairly large (perhaps even comparable to the total length of the
initial data segment, $M$, if the ringdown portion is short compared to the
necessary smooth taper).

The computational cost of the method at startup is $I^{3}$ to compute the SVD
and $\left(V^{\dagger}\tilde{C}^{-1}U\right)^{-1}$, $IM$ to compute the matrix
multiplications, and $M\log M$ to compute the PSD that determines
$\tilde{C}^{-1}$; at each iteration the cost is $IM+I^{2}$ to implement the
projection $F$. This is compared to $N^{3}$ at startup and $N^{2}$ for the
pure time-domain likelihood in Eq. (48). For $I\ll M$ (short taper, long
signal) the frequency-domain likelihood will be more efficient; for $I\sim M$
(long taper, short signal), usually the case for ringdown signals, the time-
domain likelihood is more efficient.

###  III.3 Implementation

We will use the time domain likelihood of Eq. (48), as we did in [39]. Besides
this function, we also need to specify priors for the signal parameters. In
the following, we apply uniform priors almost exclusively, but the quantities
for which the prior is uniform depends on the specific parameterization used,
as detailed in each respective section. For instance, when using the
elliptical QNM template of Eq. (14), we will apply flat priors on the
parameters in Eq. (17) or, if deviations from Kerr are allowed, Eq. (22). The
elliptical-mode parameterization has the advantage of allowing us to directly
define priors for the overall amplitude of each QNM, while remaining agnostic
about polarization content; this would not be straightforward had we chosen a
different parameterization, e.g., in terms of the right- and left-handed
polarization amplitudes $|C_{\ell\pm mn}^{\prime}|$ in Eq. (5), as in [52, 53,
44]: a uniform prior on those quantities would induce a triangular prior on
the overall QNM amplitude $A_{\ell|m|n}$, disfavoring $A_{\ell|m|n}=0$.

In order to estimate posterior densities, it is straightforward to implement
the time-domain likelihood of Eq. (48) to work with any standard Markov chain
Monte Carlo (MCMC) sampler. We choose to use the Stan platform [56, 57] in
order to take advantage of its efficient no-U-turn Hamiltonian Monte Carlo
(HMC) algorithm [88, 89]. Whatever the sampler used, it is often the case that
parameterizations of the likelihood that are most conceptually meaningful are
ill suited for sampling purposes. This is certainly true for our elliptical-
mode template, Eq. (14): although the
$\\{A_{j},\epsilon_{j},\theta_{j},\phi_{j}\\}$ quantities are a natural choice
due to their self-evident meaning, they are cumbersome to sample directly. For
example, the $\theta_{j}$ parameter becomes meaningless when $\epsilon=\pm 1$,
resulting in an ambiguity that carries over also to $\phi_{j}$ (cf., Fig. 1).
This can be circumvented by reparameterizing the likelihood in terms of more
convenient variables, and applying the requisite Jacobian. When sampling with
the elliptical-mode template, we have found that efficiency is optimized with
the reparameterization offered in App. D.

To study measurements in colored noise, we synthesize data from the Advanced
LIGO (aLIGO) design PSD [90]; when studying multiple-detector measurements, we
use that PSD for all instruments to avoid confounding factors arising from a
heterogeneous detector network. Unless otherwise stated, we simulate the
treatment of real data by first generating 4096 s of time-domain noise sampled
at 8192 Hz, then bandpassing between $[20,\,2048]\,\mathrm{Hz}$ (an
appropriate range for the signals we study), downsampling to a sample rate
4096 Hz and, finally, truncating at $t_{0}$; we always inject signals _before_
conditioning. As we would with real data, we estimate the ACF from a long
stretch of signal-free noise samples. We speed up the likelihood computation
through the Cholesky decomposition of $C$, Eq. (40), but one could also take
advantage of algorithms tailored to Toeplitz systems, like Levinson recursion
[91, 92].

##  IV Analyzing damped sinusoids

We first demonstrate some basic features of ringdown analyses under our
framework, starting from simple superpositions of damped sinusoids without the
complexities of the template introduced in Sec. II.1. For a single data stream
at times $t\geq t_{0}$, the simplest template that captures the features in
which we are interested is a generic superposition of some $D$ damped
sinusoids,

$s(t)=\sum_{n=0}^{D-1}A_{n}e^{-(t-t_{\rm ref})/\tau_{n}}\cos[2\pi
f_{n}(t-t_{\rm ref})+\phi_{n}]\,,$ (62)

with frequencies $f_{n}$ and damping times $\tau_{n}$, as well as amplitudes
$A_{n}$ and phases $\phi_{n}$ defined at some arbitrary reference time $t_{\rm
ref}$ (unless otherwise stated, we will set $t_{\rm ref}=t_{0}$ by default).
To prevent label switching degeneracies, it is helpful to impose an ordering
on $f_{n}$ or $\tau_{n}$; since BH overtones will be the main focus of our
analysis (Sec. V), in this paper we define $\tau_{n}>\tau_{n+1}$. (This can be
implemented through a coordinate transformation to an un-constrained space
[93, 56].) Even for enhanced models that account for multiple polarizations
and detector sensitivities as in Eq. (13), the goal will always be to infer
some set of $\\{f_{n},\tau_{n},A_{n},\phi_{n}\\}$ from the data using the
acyclic time-domain likelihood of Eq. (48).

Figure 10: Parameter recovery of our method (solid green) versus the our
implementation of the original LIGO-Virgo method in [82] (dashed red), on a
ringup-ringdown injection (bottom right diagram). From the peak onward (solid
trace in diagram), the injection is a simple damped sinusoid with frequency
(abscissa) and damping time (ordinate) indicated by the dotted lines. Contours
represent the 90%-credible levels obtained from each method for varying
injected postpeak SNR (line thickness and tone), Eq. (53); the top and right
panels show the corresponding marginals. The red distributions are biased by
the prepeak data, which is coupled into the analysis by the frequency-domain
analysis method as discussed in the introduction to Sec. III.)

###  IV.1 Method efficacy

By construction, our template is left undefined for times before $t_{0}$. This
means that we can apply it even if we do not expect data at those times to be
at all described by Eq. (62). We demonstrate this in Fig. 10, where we
synthesize noise from the aLIGO design PSD [90], and add it to a signal
conforming to Eq. (62) only after some $t_{0}$. We make use of a “ringup-
ringdown” injection,

$s_{\rm inj}(t)=A\,e^{-|t-t_{0}|/\tau}\cos[2\pi f(t-t_{0})]\,,$ (63)

for all times before and after some arbitrary $t_{0}$, with $f=250\,{\rm Hz}$
and $\tau=4\,{\rm ms}$ (consistent with the fundamental mode of the GW150914
remnant [82, 39]), and $A$ chosen to yield the post-peak SNRs indicated by the
legend. We condition the data as we normally would (including filtering and
downsampling), and discard times before $t_{0}$ to carry out a Bayesian
analysis based on Eq. (48), with a single damped sinusoid as a template. The
green distributions in Fig. 10 show that we are able to correctly infer the
properties of the injected signal, even though the injection departs from Eq.
(62) before $t_{0}$. On the other hand, a different method based on Eq. (37)
shows a bias induced by the pre-$t_{0}$ data (see App. E for details).

The example in Fig. 10 is designed to concretely illustrate the need for our
time-domain formalism, with the ringup-ringdown injection as a proxy for the
full inspiral-merger-ringdown signal. When other analysis strategies are
adopted, the nature and magnitude of the bias depends on the specific
properties of the noise and signal, as well as the details of the method
implementation (including conditioning). In general, we may expect the bias to
be more pronounced when targeting short-lived modes that are only visible near
the truncation time, where contamination from pre-$t_{0}$ data will tend to be
greatest. This is why we found the time-domain formalism to be indispensable
when searching for overtones in the GW150914 ringdown [39].

To demonstrate the robustness of our method, we must go beyond the concrete
example in Fig. 10 to show that we can recover arbitrary signals without
systematic bias, as long as they conform to Eq. (62) after $t_{0}$ and no
matter what came before it. We do this through a probability-probability (PP)
test [94, 50, 95, 96], also called “simulation based calibration.” We simulate
the measurement of a large number of signals drawn from our prior, with random
instantiations of noise drawn from our likelihood. If our method is working
properly, we expect that the resulting posteriors will be such that the true
values are recovered with X% credibility for X% of the simulations.
Importantly, we do not only simulate times $t\geq t_{0}$, but instead make use
of ringup-ringdown templates similar to Eq. (63), and put them through the
same conditioning process that we would apply to real data (Sec. III.3). This
ensures that we are sensitive to any potential degradation induced by the
bandpassing and downsampling filters. Therefore, the PP test allows us to
validate the entirety of our setup: conditioning, truncation, and likelihood
treatment.

Figure 11: PP plot. Deviation from uniformity in the recovered quantiles for
300 injections drawn from the prior, as described in the main text. Gray bands
mark the expected $1\sigma$, $2\sigma$ and $3\sigma$ variations under the null
hypothesis.

Since we will later focus mostly on two-mode models, we pick $D=2$ for this
test and generalize Eq. (63) so that our injections are now described by

$s_{\rm inj}(t)=\sum_{n}A_{n}\,e^{-|t-t_{0}|/\tilde{\tau}_{n}(t)}\cos[2\pi
f_{n}(t-t_{0})+\phi_{n}]\,,$ (64)

where we sum over two modes ($n=0,\,1$), and allow the pre-$t_{0}$ $e$-folding
time $\tau^{\prime}_{n}$ to differ from the regular post-$t_{0}$ value
$\tau_{n}$ by defining the piecewise auxiliary function
$\tilde{\tau}_{n}(t)\equiv\tau^{\prime}_{n}\Theta(t_{0}-t)+\tau_{n}\Theta(t-t_{0})$,
in terms of Heaviside steps $\Theta$. This choice is designed to provide a
variety of pre-$t_{0}$ morphologies, while preserving smoothness across the
$t_{0}$ boundary. For each simulation, we synthesize time-domain noise based
on the aLIGO design PSD, and draw injection parameters
$A_{n}\sim\mathcal{U}(0,3\times 10^{-21})$, $\phi_{n}\sim\mathcal{U}(0,2\pi)$,
$f_{n}/\mathrm{Hz}\sim\mathcal{U}(635,1285)$,
$\ln(\gamma_{n}/\mathrm{Hz})\sim\mathcal{U}(\ln 635,\ln 1285)$ and same for
$\ln(\gamma_{n}^{\prime}/\mathrm{Hz})$, for both $n=0,1$, and enforce
$\tau_{0}>\tau_{1}$ with no imposed ordering on $\tau_{0/1}^{\prime}$; these
distributions match our prior—except for the $\gamma_{n}^{\prime}\equiv
1/\tau_{n}^{\prime}$ parameters, which only control the pre-$t_{0}$ morphology
of the injection and, thus, are not part of our recovery template, Eq. (62).
The resulting posteriors conform to our statistical expectations, as shown by
the agreement between the measurements (colored lines) with our statistical
expectation (gray bands) in Fig. 11.

###  IV.2 Detecting modes

Figure 12: Amplitude posteriors from the analysis of two damped sinusoids
($D_{\rm true}=2$) in noisy data, using templates with a varying number of
modes in Eq. (62): $D=1$ (blue), $D=2$ (orange), $D=3$ (green). The diagonal
displays amplitude marginals, while the upper and lower corners respectively
show credible levels (90%, 50%, 10%) and individual samples; dotted lines mark
the true values (Table 1). Noise follows the aLIGO design PSD, resulting in an
injected matched filter SNR of 20, as defined in Eq. (53). The posterior for
$D=3$ allows some of the amplitudes to vanish, so we conclude $D_{\rm
best}=2$.

A template like Eq. (62) can accommodate an arbitrary number of damped
sinusoids $D$. Actual BBH merger ringdowns in GR excite a (countably) infinite
number of modes. Determining how many modes are _necessary_ to adequately
describe a given set of data is a model selection problem: to identify the
number of modes that are clearly detected in a fully Bayesian setting, we
could compare the relative evidences of models made up of different
combinations of modes. Formally, this procedure is encapsulated by the
Bayesian odds between models, or just their Bayes factors if we have no _a
priori_ preference between the different alternatives. However, Bayes factors
are only meaningful insofar as our priors actually represent our belief about
the expected distribution of parameter values. Assigning compelling priors to
the ringdown parameters can be quite challenging, especially when it comes to
the damped sinusoid amplitudes and relative phases. Therefore we do not
advocate for using relative evidences between models to choose the “correct”
number of modes in a signal. The procedure we outline below is considerably
less sensitive to prior choices, depending only on how those choices change
the _shape_ of the inferred posterior, unlike the model evidences that depend
on the prior in parameter regions without any posterior support.

We propose to determine the number of modes $D_{\rm best}$ that the data
require by demanding that:

  1. 1.

if $D=D_{\rm best}$, the posterior indicate $A_{n}>0$ with a reasonable degree
of certainty (we use 90% credibility in this work), for _all_ $0\leq n<D$;
while,

  2. 2.

if $D=D_{\rm best}+1$, this is no longer be the case.

In other words, we stop adding modes to the template once the only effect of
doing so is to introduce degeneracies. This procedure relies on the existence
of an intrinsic hierarchy of modes, so that there is a natural order in which
to progressively add them to the model. For example, modes with shorter
durations are generally harder to measure, so we may add modes in order of
decreasing damping time (suitable for a series of overtones). Alternatively,
we can add modes to the model in order of expected amplitude (suitable for a
series of angular harmonics).

We illustrate this in Fig. 12, where we inject two damped sinusoids ($D_{\rm
true}=2$) into Gaussian noise and carry out measurements with $D=1,2,3$ in Eq.
(62). The posterior for a single-mode template (blue distribution, top left)
confidently indicates that the data contain a signal, inconsistent with zero
amplitude. Including a second damped sinusoid in the template results in a
posterior (orange distribution) that is consistent with the injection, and
that offers support for neither $A_{0}=0$ nor $A_{1}=0$. This stops being true
with the addition of a third mode, in which case the posterior (green
distribution) allows any _one_ of the amplitudes to vanish, as long as the
other two do not. The fact that no two amplitudes may vanish simultaneously is
reflected in the characteristic “arching” of the 2D joint marginals for each
amplitude pair, excluding $A_{n}=A_{m}=0$ for any choice of $n\neq m$. At this
point, it is no longer helpful to continue adding modes to the template, so we
conclude $D_{\rm best}=2$ for this example.

Figure 13: Effect of phasing on mode detectability, illustrated through a pair
of $D=2$ signals sharing a spectrum $f_{n},\tau_{n}$, and mode-amplitude ratio
$A_{1}/A_{0}$, but with different mode phases $\phi_{n}$. Although both
signals integrate to the same SNR, the $n=1$ mode is down-weighted in example
A (blue) relative to example B (orange), potentially making it harder to
detect.

The number of modes required by the data $D_{\rm best}$ need not correspond to
the true number $D_{\rm true}$, and it certainly will not in the analysis of
true (as opposed to synthetic) BH ringdowns, for which $D_{\rm true}=\infty$
always. Since $D_{\rm best}$ is the number of _discernible_ modes, it will
depend on the overall SNR. It will also be contingent on the specific
characteristics of the signal (particularly the relative mode phases in the
true signal) and the specific noise realization. Any of those factors may
impact the features of the $D$-dimensional amplitude posterior, and therefore
our determination of $D_{\rm best}$. For instance, identifying a critically
damped ($Q=1/2$) or overdamped ($Q<1/2$) mode would require exceptionally high
SNR, since such modes do not exhibit any cycles. Moreover, even for
underdamped signals, the specific phasing determines the distribution of SNR
across modes, and can thus have a big impact on $D_{\rm best}$ (Fig. 13).

Figure 14: Reconstruction of two damped sinusoids injected into simulated
aLIGO noise, as recovered using a model with $D=2$ modes (orange distribution
in Fig. 12). The top panel shows the injected signal (blue line) and the
90%-credible reconstruction (blue shading), laid over the noisy data (black).
The bottom panel shows the two individual modes making up the injection
(lines), and their respective reconstructions (shading). The true parameters
are based on a GW150914-like remnant BH (Table 1); the shortest-lived mode
($n=1$) exhibits little over one full cycle.

Having identified certain number of damped sinusoids in the data, we may
reconstruct them (Fig. 14) and quantify their properties. Our definition of
$D_{\rm best}$ is designed such that this a meaningful thing to do: if the
data favored $A_{n}=0$ for some $n$, then the posterior on the corresponding
$f_{n}$ and $\tau_{n}$ would necessarily span the range of the prior, making
constraints harder to interpret (e.g., the 90%-credible interval would
necessarily vary with the prior bounds, which are usually arbitrary).

Table 1: GW150914-like parameters injected in Fig. 12. $n$ |  $f_{n}$ |  $\tau_{n}$ |  $A_{n}$ | $\phi_{n}$
---|---|---|---|---
0 | 250 Hz | 4.0 ms |  $3\times 10^{-21}$ |  $5.34$ rad
1 | 245 Hz | 1.4 ms |  $4\times 10^{-21}$ |  $1.79$ rad

###  IV.3 Characterizing modes

Figure 15: Joint frequency (abscissa) and damping rate (ordinate) posterior
for the longest- (blue) and shortest- (orange) lived damped sinusoids in the
$D=2$ analysis from Fig. 12. Contours enclose 90%, 50%, and 10% of the
probability mass; crosses mark the injected values (Table 1). For comparison
to Eqs. (65) and (66), we show the symmetric 68%-credible ($1\sigma$)
intervals for each marginal (dashed lines, top and right). Note that this
figure hides some of the information yielded by the Bayesian analysis, since
it does not show the correlations between the $n=0$ and $n=1$ parameters.

####  IV.3.1 Resolving modes

Having a preference for $A_{n}>0$ is also a sufficient condition for a
meaningful measurement of $f_{n}$ and $\tau_{n}$ under a specific model. This
is in spite of claims to the contrary in the ringdown literature. For
instance, Refs. [40, 42] presuppose that a frequency measurement in an $D=2$
model is only possible if the $f_{n}$ are “resolved,” by which they mean that
the frequencies satisfy

$|f_{0}-f_{1}|>\max(\sigma_{f_{0}},\sigma_{f_{1}})\,,$ (65)

where $\sigma_{f_{n}}$ is the standard deviation of the $f_{n}$ posterior
samples. By the same token, constraining the damping times would only be
possible if

$|\tau_{0}-\tau_{1}|>\max(\sigma_{\tau_{0}},\sigma_{\tau_{1}})\,.$ (66)

An insistence that modes be spectrally separated in this way led Refs. [40,
42] to conclude that a constraint of the overtone frequencies in a BH ringdown
would require $\mathrm{SNR}\sim 100$, disagreeing with the results in [39,
44].

Figure 16: Reconstruction of signal made up of two damped sinusoids (color)
with equal frequency $f_{0}=f_{1}=250\,\mathrm{Hz}$ (abscissa) but different
damping times $\tau_{0}=20\,\mathrm{ms}$ and $\tau_{1}=50\,\mathrm{ms}$
(ordinate, left panel), that is, damping rates $\gamma_{0}=50\,\mathrm{Hz}$
and $\gamma_{1}=200\,\mathrm{Hz}$ (ordinate, right panel). Contours enclose
90%, 50% and 10% of the probability mass, while vertical lines over the
marginals mark the 68% credible interval. From the left panel, it is clear
that $\tau_{1}$ is more sharply constrained than $\tau_{0}$, which is what we
expect if $\gamma_{1}$ is less constrained than $\gamma_{0}$ (right).

However, the possibility of measuring $f_{n}$ or $\tau_{n}$ does not hinge on
any such condition. In a Bayesian analysis, all damped sinusoids in the
template are modeled simultaneously and self consistently. A $(2\times
D)$-dimensional posterior captures all potential degeneracies between the
different $f_{n}$ and $\tau_{n}$. Two modes will be distinguishable whenever
the multidimensional posterior inferences about their frequencies and damping
times do not overlap. Non-overlapping frequency and damping time posteriors do
not necessarily produce non-overlapping one-dimensional marginal posteriors on
the frequency or damping time. The correct condition for resolvability is that
the separation between the frequency and damping time of two modes be large
compared to the posterior _covariance_ of these quantities. When the
posteriors do not have significant correlations between $f$ and $\tau$ the
distinguishability condition looks like

$r^{2}\equiv\frac{\left(f_{0}-f_{1}\right)^{2}}{\sigma_{f_{0}}^{2}+\sigma_{f_{1}}^{2}}+\frac{\left(\tau_{0}-\tau_{1}\right)^{2}}{\sigma_{\tau_{0}}^{2}+\sigma_{\tau_{1}}^{2}}\gtrsim
1.$ (67)

Whenever the posterior distributions for the two modes’ parameters are
distinguishable, each mode’s frequency and damping rate can be constrained. In
other words, we can meaningfully constrain the values of the $f_{n}$’s from
the data even if Eq. (65) is not satisfied, and same for $\tau_{n}$ and Eq.
(66).

As an example, consider the $D=2$ measurement in Fig. 12 (orange). We
summarize the associated frequency and damping time posteriors in Fig. 15, by
overlaying the joint distributions for $\\{f_{0},\tau_{0}\\}$, and for
$\\{f_{1},\tau_{1}\\}$ (blue and orange, respectively). The injected
frequencies are much closer than the characteristic width of the corresponding
marginals (top panel), and Eq. (65) is flagrantly violated. Yet we constrain
the frequency to be $f_{0}=252^{+8}_{-15}\,\mathrm{Hz}$ for the longest-lived
mode, and $f_{1}=229^{+36}_{-38}\,\mathrm{Hz}$ for the shortest-lived one (at
68% credibility, or $1\sigma$, to match the resolvability criteria), and the
two modes making up the signal can be reconstructed reasonably well (Fig. 14).

We chose the parameters injected in Fig. 15 to be consistent with the ringdown
of a GW150914-like remnant (Table 1). However, the observation that frequency
measurements are not predicated on Eq. (65) is general and applies to more
extreme cases. This includes hypothetical signals in which the true mode
frequencies are exactly identical. For example, in Fig. 16 we constructed an
injection with $f_{0}=f_{1}$; in spite of this, the two modes are correctly
identified thanks to their different decay rates, and we are able to measure
the two frequencies and damping times successfully. Of course, the point of
this example is not to argue that modes with similar parameters are _ideal_
for this kind of measurement, but simply that they are not _a priori_
unsuitable.

Correlations aside, the Rayleigh criteria can indeed be a useful guide in
forecasting $D_{\rm best}$ for hypothesized measurements, especially if high
SNRs can be expected [12, 97]. That is because, in the absence of cross-mode
correlations, failing _both_ Eq. (65) and Eq. (66) would imply that the data
are consistent with a shared frequency and damping time for the two modes.
Whatever their amplitudes and phases, the sum of two damped sinusoids with
identical $f$ and $\tau$ is just another damped sinusoid. Therefore, if the
data prefer $f_{0}=f_{1}$ and, simultaneously, $\tau_{0}=\tau_{1}$, then one
of the modes in the $D=2$ template can be zeroed out, and we would not have
found $D_{\rm best}=2$ in the first place, i.e., we would have stopped at
$D=1$ per the procedure outlined above (Sec. IV.2).

####  IV.3.2 Relative versus absolute precision

Figure 17: Fractional uncertainties corresponding to the measurement in Fig.
16, as proxied by the posterior on the quantites in Eq. (68). The origin
represents the location of the mean of the distributions in the left panel of
Fig. 16, $\hat{f}_{n}$ and $\hat{\tau}_{n}$, and crosses mark the relative
locations of the truth values of $f_{n}$ and $\tau_{n}$; contours enclose 90%
of the probability, while vertical lines over the marginals indicate the 68%
creidble interval. The fractional uncertainties are similar for $\tau_{0}$ and
$\tau_{1}$, but the fractional uncertainty in $f_{1}$ is larger than $f_{0}$.

Figure 16 serves to illustrate another important point: being able to separate
two modes in the $f,\tau$ plane does not imply we can measure those parameters
with high _relative_ precision. From the left panel of Fig. 16, it is clear
that the analysis can very successfully distinguish the two modes in the
signal through their damping times, and that the $\tau_{n}$ measurement is
more accurate for $n=1$ than $n=0$. This is a statement about the absolute
precision with which we can determine $\tau_{0}$ or $\tau_{1}$. Instead, we
may ask about the relative precision of these measurements by looking at the
quantities

$\delta\hat{\tau}_{n}\equiv\frac{\tau_{n}-\hat{\tau}_{n}}{\hat{\tau}_{n}}~{}~{},~{}~{}\delta\hat{f}_{n}\equiv\frac{f_{n}-\hat{f}_{n}}{\hat{f}_{n}}\,$
(68)

where $\hat{\tau}_{n}$ is the $\tau_{n}$ posterior mean as estimated from our
samples, and same for $f_{n}$. This definition gives us a proxy for fractional
deviations of the parameters around the true mean.

Figure 17 shows the posterior on these quantities as derived from Fig. 16.
Even though we were able to separate $\tau_{0}$ and $\tau_{1}$ cleanly in Fig.
16, both the $\delta\hat{\tau}_{0}$ and $\delta\hat{\tau}_{1}$ measurements
are quite broad percentage-wise: $\delta\hat{\tau}_{0}\approx\pm 20\%$ and
$\delta\hat{\tau}_{1}\approx\pm 30\%$, at 90% credibility. On the other hand,
the $\delta\hat{f}_{n}$ frequency deviations are much more tightly
constrained: $\delta\hat{f}_{0}\approx\pm 1\%$ and
$\delta\hat{f}_{1}\approx\pm 2\%$. So, even though the two modes have quite
different damping times, it is the $\delta\hat{f}$’s which are determined more
precisely. Therefore, it should come as no surprise that, when we attempt to
test the Kerr hypothesis using the first overtone of the $\ell=|m|=2$ mode, we
can better constrain fractional deviations in its frequency than its damping
time (as was found in [39, 44], and as we will show in Sec. V.2).

Incidentally, although $\tau_{n}$ is more sharply constrained for $n=1$ than
$n=0$ (Fig. 16, left), casting the measurement in terms of
$\gamma_{n}=1/\tau_{n}$ inverts the relationship between the modes:
$\gamma_{0}$ is better measured than $\gamma_{1}$ (Fig. 16, right). This
straightforward observation reveals the fact that $\omega_{n}$ and
$\gamma_{n}$ are the natural basis that diagonalizes this measurement. In any
case, this is of little consequence when it comes to the precision of
fractional deviation measurements.

This concludes our exploration of generic damped sinusoids based on the
simplified template of Eq. (62). In the following section, we turn to more
realistic examples of BH ringdowns, for which we will apply the full template
of Eq. (13). Our examples will focus on tones of the $\ell=|m|=2$ harmonic.

##  V Analyzing black-hole ringdowns

Figure 18: Simulated measurement of the least-damped $\ell=|m|=2$ mode for a
GW150914-like remnant BH. The left panel shows the joint posterior for $M$ and
$\chi$, quantities in which the measurement parameterized, with blue contours
enclosing 90%, 50% and 10% of the probability mass; the right panel shows the
implied $f_{0}$ and $\tau_{0}$ distribution. On the left (right) gray lines
mark contours of equal $f_{0}$ and $\tau_{0}$ ($M$ and $\chi$). A uniform
prior on $(M,\chi)$ induces a nonuniform prior on $(f_{0},\tau_{0})$, and
$(f_{0},\tau_{0})$ values below the $\chi=0$ contour are disallowed (gray
region on the right). Crosshairs mark the true values; parameters not shown
are $\phi_{0}=5.34\,\mathrm{rad}$, $\epsilon=-1$ and amplitude chosen to yield
ringdown $\mathrm{SNR}=14$ over the network.

The study of GW signals is not fundamentally different from the simplified
damped-sinusoid examples in the previous section. This is true even though the
existence of multiple polarizations and detectors complicates the
implementation in practice, as is reflected in the increased complexity of
both the template, Eq. (13), and the likelihood, Eq. (50).

In the general case, having two GW polarizations introduces an additional pair
of degrees of freedom per QNM, represented by $\epsilon_{j}$ and $\theta_{j}$
in Eq. (5). It also affects the structure of the template in Eq. (13), making
it always a linear combination of the two polarizations, with the antenna
patterns $F^{I}_{+/\times}$ serving as respective weights for each detector
$I$. However, since we can neglect the polarization angle $\psi$ (Sec. II.1),
the only material effect of the $F^{I}_{+/\times}$ factors consists of
regulating the expected amplitude _ratios_ of the signal as measured by
different detectors. Such ratios, as well as the relative time delays $\delta
t_{I}$, implicitly make the template depend on the source sky location. If
this was not known a priori, we would have to include the right ascension
$\alpha$ and declination $\delta$ as additional degrees of freedom in our
analysis. However, for signals seen by multiple detectors, we always have an
accurate measurement of both the $\delta t_{I}$’s and the amplitude ratios
from the full IMR analysis; for signals seen by a single detector, we
effectively only measure a single polarization, so the template reduces to Eq.
(62) and we can ignore the sky location altogether.

The IMR analysis also provides information about the overall time of arrival
of the signal, which is important in choosing a start time for the ringdown
analysis. Specifically, we can use the IMR results to reconstruct the time,
$t_{\rm peak}$, at which the peak of the complex strain envelope,
$|h|^{2}=h_{+}^{2}+h_{\times}^{2}$, arrived at the geocenter, with $t_{\rm
peak}^{I}\equiv t_{\rm peak}+\delta t_{I}$ the corresponding times at each
detector. We can use this reference point to define the truncation times,
$t_{0}^{I}$, at which to start the ringdown analysis [cf. Eq. (47)].
Generally, the peak time can be reconstructed with sufficient accuracy.
Otherwise, we could incorporate uncertainty in the arrival time and sky
location by sampling over $\alpha$ and $\delta$ without varying the truncation
time, as long as we ensured that the truncation time always lies on or after
whatever waveform feature we use to define the start of the ringdown, but we
do not implement such a model in this work.

Table 2: GW150914-like extrinsic parameters. $\alpha$ | $\delta$ | $\psi$ |  $t_{0}$ (GPS)
---|---|---|---
1.95 rad |  $-1.27$ rad | 0.82 rad | 1126259462.423 s

In the remainder of this paper, we will simulate GW signals as originating
from a sky location consistent with GW150914, as specified in Table 2.
Although $\psi$ plays no role in the ringdown analysis, an arbitrary choice
must be made in the injection in order for the amplitudes and phases at each
detector to be uniquely specified by
$\\{A_{j},\epsilon_{j},\theta_{j},\phi_{j}\\}$; our choice, shown in Table 2,
is the same as in [39]. We also specify a geocenter-based reference GPS time
$t_{0}$, whose only relevance is in the computation of the antenna patterns
and inter-detector delays.

After properly accounting for all this additional complexity, the lessons from
Sec. IV carry over directly to the GW case, including the discussions of mode
detectability and resolvability. We can use our infrastructure to measure the
properties of Kerr BHs, and to look for deviations away from the Kerr
spectrum. With recent observational results in mind [39, 44], we focus on
modes with $\ell=|m|=2$ but many of the observations are applicable to
arbitrary modes including other harmonics. Unless otherwise stated, we place
uniform priors on all quantities mentioned.

###  V.1 Kerr ringdowns

Within GR, an interesting application of BH spectroscopy is to measure the
mass and spin of a Kerr BH from ringdown data alone (e.g., [98, 99, 39, 52,
53, 44, 45, 47]). Under our Bayesian framework, the most straightforward way
of doing this consists of parameterizing the frequencies and damping times in
Eq. (13) as a function of $M$ and $\chi$ to directly obtain a posterior on
these quantities (Sec. II.2). The template can include as many QNMs as deemed
appropriate, e.g., by following the procedure outlined in Sec. IV.2.

Figure 19: As in Fig. 18, but for a template and injection that incorporate
the first overtone, in addition to the least-damped $\ell=|m|=2$ mode
(orange). For comparison, we also show the 90%-credible level from Fig. 18
(blue). Solid gray lines mark contours of constant $(f_{1},\,\tau_{1})$;
dotted, fainter lines mark $(f_{0},\,\tau_{0})$. Unlike Fig. 18, the right
panel here shows the implied posterior on the overtone, rather than
fundamental, frequency and damping rate; since the simulation in Fig. 18 only
included the fundamental, the blue contour on the right is fully determined by
the $n=0$ parameters, which are in one-to-one correspondence thanks to the
Kerr assumption. For this configuration, the overtone sharpens the $(M,\chi)$
recovery even with fixed SNR, but this behavior is not universal. Both here
and in Fig. 18, the ringdown network SNR is 14.

As a simple example, Fig. 18 shows a simulated mass and spin measurement based
on the $\ell=|m|=2,\,n=0$ mode of a GW150914-like remnant BH with
$M=68\,M_{\odot}$ and $\chi=0.69$, and other injected parameters as in Table
3. We simulate the signal as would be detected by the LIGO Hanford and
Livingston instruments, with noise corresponding to the aLIGO design
sensitivity (see Sec. III.3), and rescale the signal amplitude to yield a
network SNR of 14. The mass and spin are recovered correctly, and their joint
posterior shows a characteristic, elongated shape that extends along an equal-
frequency contour. In fact, although the priors are diagonal in $M$ and
$\chi$, the likelihood is effectively diagonal in $f_{0}$ and $\tau_{0}$
instead. The uniform $(M,\,\chi)$ prior translates into a highly nonuniform
$(f_{0},\,\tau_{0})$ prior, reflecting the coordinate geometry implied by the
Kerr spectrum per Fig. 2. In particular, this results in the exclusion of
points with $Q\lesssim 2.1$, which are unachievable for the fundamental
$\ell=|m|=2$ mode with any combination of mass and spin (gray region).

Table 3: Kerr $\ell=|m|=2$ tones injected in Figs. 18–21. Only the amplitude ratio $A_{1}/A_{0}$ is relevant here because we rescale the signal to a given SNR; $\theta_{n}$ is undefined for $|\epsilon|=1$. $M$ | $\chi$ | $A_{1}/A_{0}$ | $\epsilon_{0}$ | $\epsilon_{1}$ | $\phi_{0}$ | $\phi_{1}$
---|---|---|---|---|---|---
68 $M_{\odot}$ | 0.69 | 1.36 | $-1$ | $-1$ | 5.34 rad | 1.79 rad

It is possible for additional modes, if present, to reduce the uncertainty in
the $(M,\,\chi)$ measurement, even for fixed SNR. We show an example of this
in Fig. 19, where we add the first overtone ($n=1$) to the signal in Fig. 18,
while rescaling the amplitudes so as to keep the SNR unchanged. Improvements
of this kind are contingent on the additional modes actually appearing in the
data in a way that provides new information to narrow the viable $(M,\,\chi)$
space and that overcomes the increased uncertainty due to a weakened
fundamental mode (since we are keeping the SNR constant; e.g., see Fig. 13).
As underscored in Sec. IV, whether this is the case on any specific instance
will depend on the particular characteristics of the signal and noise
instantiation (including nuisance parameters, like $\phi_{n}$, $\epsilon_{n}$
or $\theta_{n}$). In Fig. 19, the quality of the $(M,\,\chi)$ measurement
depends nontrivially on the interplay between the loss of information due to a
weakened fundamental and the gain afforded by the overtone.

In realistic situations, we are interested in overtones because they can also
increase the available ringdown SNR when analyzing real GW data [38, 39],
which is not the scenario demonstrated in Fig. 19. Such an enhancement must
generally be balanced against the growth in uncertainty associated with an
increased number of degrees of freedom. Even when additional tones do
contribute, the longest-lived mode will usually be the best constrained and
the structure of the $(M,\,\chi)$ posterior will not differ qualitatively from
that in Fig. 18.

Figure 20: Effect of increasing the number of detectors on the $N=1$ analysis
in Fig. 19, as manifested on the recovered mass $M$, spin $\chi$, overtone-to-
fundamental amplitude ratio $A_{1}/A_{0}$, and fundamental ellipticity
$\epsilon_{0}$, with 90% and 10% credibility. Color corresponds to networks
containing interferometers at the locations of LIGO Hanford (H), LIGO
Livingston (L), Virgo (V) and Kagra (K), although we take all detectors to
have the design sensitivity of aLIGO. The blue (HL) curves correspond to the
orange result in Fig. 19, and we enforce a network SNR of 14. Although not
shown here, the result for $\epsilon_{1}$ is comparable to that for
$\epsilon_{0}$; we have also omitted $\theta_{0/1}$ because that parameter is
undefined for $|\epsilon_{0/1}|=1$, as is the case here (Table 3).

For a given network SNR, having multiple detectors does not fundamentally
change the above picture for $M$ and $\chi$ either. Nevertheless, a larger
network can allow for better inferences on polarization-related quantities,
like $\epsilon_{n}$ and $\theta_{n}$, which can in turn lead to more accurate
measurement of the spectrum by breaking degeneracies. In Fig. 20 we show the
effect of increasing the number of detectors for the $N=1$ analysis in Fig.
19, keeping the network SNR constant. In this case, the amplitude ratio
$A_{1}/A_{0}$ and the ellipticity $\epsilon_{0}$ are both better constrained
with a larger network, and $(M,\,\chi)$ are recovered with slightly better
accuracy. Of course, in a realistic situation, having more detectors usually
also brings a higher network SNR, which, all else being equal, always leads to
a better measurement. Both the SNR gain and the breaking of polarization
degeneracies become more important for templates with larger numbers of modes.
Unfortunately, the two LIGO detectors are nearly coaligned, so that their
ability to distinguishing GW polarizations is diminished.

Finally, there are circumstances in which it may be acceptable to use a
restricted version of the elliptical model in Eq. (13). As mentioned above,
this is the case for signals recorded by a single detector (or, to good
approximation, by just the two LIGO detectors), for which the polarization
parameters $\epsilon_{j}$ and $\theta_{j}$ are fully redundant. It is also
true more generally whenever the inclination is known and the system can be
assumed to be (reflection) symmetric over the equator; if so, it may be
reasonable to take all the modes in the signal to have a known ellipticity,
e.g., $\epsilon=+1\,(-1)$ for a face on (off) system (see App. B). With such a
model, the recovery of amplitudes and phases is aided by the reduced number of
parameters (Fig. 21). However, this does not necessarily improve the accuracy
of the $M$ and $\chi$ measurements.

Figure 21:  Effect of assuming a known ellipticity for all modes in the $N=1$
analysis of Fig. 19, as reflected on the BH mass $M$, spin $\chi$, mode
amplitude ratio $A_{1}/A_{0}$, and mode relative phase $\phi_{1}-\phi_{0}$,
with 90% and 10% credibility (contours). Blue distributions were obtained
assuming a circular-polarization model with $\epsilon_{0}=\epsilon_{1}=-1$, as
in the injected signal; the orange distributions come from the same analysis
as in Fig. 19.

###  V.2 Deviations from Kerr

Figure 22: Simulated constraints on Kerr deviations from a model including the
two longest-lived tones with $\ell=|m|=2$, using the parameterization of Eq.
(27) on the Kerr $N=1$ injection from Fig. 19. The corner plot shows posterior
densities for the BH mass $M$, spin $\chi$, overtone frequency deviation
$\delta f_{1}$, and damping-time deviation $\delta\tau_{1}$ (blue, “Kerr+”);
for comparison, we overlay the $(M,\,\chi)$ measurement obtained assuming
$\delta f_{1}=\delta\tau_{1}=0$ (orange, “Kerr;” same as in Fig. 19). In this
plot, we have marginalized over nuisance parameters
$\\{A_{n},\epsilon_{n},\theta_{n},\phi_{n}\\}$ for both tones in the template
($n=0,1$), contours enclose 90%, 50% and 10% of the probability, and vertical
lines over the marginals mark the 90% credible interval; dotted straight lines
indicate the truth.

Beyond mass and spin measurements within GR, we can look for signs of new
physics by allowing for deviations from the Kerr prediction. As discussed in
Sec. II.3, the most natural way to do this is to introduce Kerr deviation
parameters $\delta f$ and $\delta\tau$ (or $\delta f$ and $\delta\gamma$) for
the frequency and damping time of one or more QNMs (at most $2D-2$ new
parameters for a model with $D$ modes). We can use this parameterization to
establish whether any two or more damped sinusoids are consistent with a Kerr
spectrum, obtaining a quantitative answer to which we can assign a definite
credibility.

This kind of analysis can be carried out using however many modes are
confidently detected in the data (Sec. IV.2), whether they be tones of a given
angular harmonic (same $\ell$ and $|m|$, different $n$), a series of
fundamental modes with varying angular structure (different $\ell$ or $|m|$,
same $n$), or any other combination. Continuing our focus on the two-tone
$\ell=|m|=2,\,N=1$ model, we will study simulated measurements of the $\delta
f_{221}$ and $\delta\tau_{221}$ overtone parameters of Eqs. (27) (henceforth
just denoted $\delta f_{1}$ and $\delta\tau_{1}$ for simplicity). This serves
the dual purpose of demonstrating our procedure and validating our
infrastructure, as well as replicating some of the conclusions obtained with
real data in [39, 44], verifying they are in agreement with
expectation.131313However, note that neither the model nor the infrastructure
applied here correspond exactly to those used in [39] or [44]. The former used
a model like the one described in App. B, while the latter applied different
priors; both used different samplers.

In Sec. V.2.1, we demonstrate constraints based on signals that do follow a
Kerr spectrum; in Sec. V.2.2, we demonstrate measurements obtained from non-
Kerr signals.

####  V.2.1 Constraints from Kerr signals

We begin by reanalyzing the $\ell=|m|=2$, $N=1$ Kerr signal from the previous
section (Fig. 19, orange), this time extending our recovery model through the
two Kerr-deviation parameters $\delta f_{1}$ and $\delta\tau_{1}$. We show the
resulting posterior on these parameters, as well as $M$ and $\chi$, in Fig. 22
(blue). For this particular configuration, we obtain a reasonable constraint
on the overtone frequency deviation, measuring $\delta
f_{1}=-0.12^{+0.44}_{-0.47}$ with 90% credibility, or $\delta
f_{1}=-0.12^{+0.34}_{-0.28}$ with 68% credibility; on the other hand, the
damping rate is only poorly constrained to $\delta\tau_{1}>-0.34$ with 90%
credibility, with support extending to the upper edge of our prior,
$\delta\tau_{1}=1$. Unsurprisingly, the presence of these additional
parameters results in a broader range of allowed $M$ and $\chi$ values, as
seen in the comparison to the Kerr result (orange) in the upper left corner of
Fig. 22.

Figure 23 shows the corresponding waveform reconstructions under both the Kerr
and non-Kerr models: both are consistent with the true signal with comparable
uncertainty on the overall reconstruction, but the latter shows greater
uncertainty on the individual modes. This is a direct manifestation of the
flexibility introduced by $\delta f_{1}$ and $\delta\tau_{1}$, which opens up
more combinations of mode morphology that can add up to the observed signal.
For this configuration, the effect of the overtone on the overall waveform is
to lower the amplitude near $t_{0}$ by interfering destructively with the
fundamental mode. At late times, once the overtone has decayed, the signal is
fully given by the fundamental mode to great accuracy.

Figure 23: Waveform reconstructions from the two analyses in Fig. 22, assuming
a Kerr spectrum (top), or allowing deviations in the overtone (bottom). Solid
lines show the true waveforms for the injected fundamental mode (blue),
overtone (orange) and total signal (magenta), with colored envelopes
representing the reconstruction at the 90% credible level. We only show the
result for the Hanford (H1) observatory, since the Livingston reconstruction
is similar.

The result from this simulation is consistent with what was found in [39], in
particular with the expectation that, at moderate SNRs, $\delta\tau_{1}$
should be generally harder to pin down than $\delta f_{1}$ (also, e.g., [27]).
In this case, the reason has to do with the correlations between the overtone
deviations and $M$ and $\chi$, as displayed in Fig. 22: the model is able to
support high values of $\delta\tau_{1}$, as long as both $M$ and $\chi$ are
lowered with respect to their true values. This degeneracy is not fully
intrinsic to the $N=1$ model but arises from the way in which the two damped
sinusoids combine to produce the true injected signal in this example, as
determined by their relative phasing (here simply the difference
$\phi_{1}-\phi_{0}$, since $\epsilon=-1$ for both modes in the injection, so
the $\theta_{n}$ are meaningless).

Although the posterior structure will generally vary for different data, this
specific configuration is representative of what we expect from binary BH
numerical relativity simulations, so it is worth examining it further. For a
Kerr BH, smaller $M$ results in overall higher frequencies and lower damping
times for all modes. If the change in frequency is smaller than ${\sim}50\%$
(at this SNR), it can be counterbalanced by a decrease in $\chi$, which
however will also result in lower damping times (left and center panels in
Fig. 2). By decreasing both $M$ and $\chi$, one can thus obtain a spectrum of
modes with unaltered frequencies but faster decays. Meanwhile,
$\delta\tau_{1}$ frees the overtone from following this pattern, so that a
lower $M$ and $\chi$ can be used to decrease $\tau_{0}$ independently of
$\tau_{1}$, while keeping both $f_{0}$ and $f_{1}$ approximately constant. In
fact, these correlations can be used to decrease $\tau_{0}$ while
simultaneously _increasing_ $\tau_{1}$ through $\delta\tau_{1}$. If one also
alters the amplitudes so as to lower the ratio $A_{1}/A_{0}$, modifications to
the waveform induced by the changes to the damping times can be made to cancel
out, since one can compensate a shorter (longer) damping time with a larger
(smaller) initial amplitude. While $\delta\tau_{1}$ can vary over a broad
range, $\delta f_{1}$ is more constrained, since the frequency of the overtone
must remain close to the true value to properly interfere with the frequency
of the fundamental mode (which is itself well determined by data at late
times). All this is encoded in Figs. 23 and 24.

Figure 24: As in bottom panel of Fig. 23, but limiting $\delta\tau_{1}>+0.6$
(top) or $\delta\tau_{1}<-0.6$ (bottom). High values of $\delta\tau_{1}$ can
be accommodated with small changes to the overall signal (magenta) by reducing
$M$ and $\chi$ with respect to the true values, in such a way as to shorten
the fundamental mode (blue) while keeping its frequency roughly the same, and
by simultaneously tuning the amplitudes to decrease the overtone (orange) to
fundamental ratio, $A_{1}/A_{0}$. No such perfect tuning is possible for
overly negative values of $\delta\tau_{1}$, as highlighted in the inset where
there is an evident disagreement between the 90%-credible reconstruction and
the true signal.

Similar observations explain the asymmetry between positive and negative
values of $\delta\tau_{1}$ in Fig. 22. For sufficiently negative values of
$\delta\tau_{1}$, the overtone vanishes too quickly and the balance between
the two modes cannot be restored by changes to the other parameters. As
$\delta\tau_{1}\rightarrow-1$, the overtone approaches a Kronecker delta at
$t_{0}$ and its contribution to the total signal vanishes (except at the first
sample), even for high amplitudes $A_{1}$. For such values of
$\delta\tau_{1}$, then, the amplitude of the fundamental mode must be
decreased to match the signal near $t_{0}$, but this is inconsistent with
later times, during which the longest-lived mode dominates the injection.
Since the fundamental cannot be made to match both the early and late data in
this regime (see inset in Fig. 24), negative values of $\delta\tau_{1}$ are
disfavored.

Note that, had the presence of the overtone been unclear (e.g., if the
posterior had significant support for $A_{1}=0$), then letting
$\delta\tau_{1}\rightarrow-1$ would be a way for the sampler to get rid of
this fast-decaying mode for arbitrary values of the $n=1$ parameters; this can
cause an artificial preference for extremely negative values of
$\delta\tau_{1}$, and result in sampling problems. The LIGO-Virgo
collaboration has previously encountered this issue when analyzing low-SNR
ringdowns [44]. This issue could probably be eliminated by reparameterizing
the modifications in terms of $\gamma=1/\tau$ and $\delta\gamma$ because
$\delta\tau\to-1$ corresponds to $\delta\gamma\to\infty$; alternately, the
deviation parameters could be redefined non-linearly as $\tau\to\tau
e^{\delta\tau}$, with the same result that the problematic point
$\delta\tau\to-1$ is removed to infinity. We leave explorations of these
alternative parameterizations to future work.

Figure 25: Constraints on Kerr deviations for increasing ringdown SNR (color).
Same beyond-Kerr analysis as in Fig. 22, but for varying injected network SNR.
Contours enclose 90% of the probability mass, and dotted lines mark the truth.
As the SNR increases, the analysis converges on the correct values for $M$,
$\chi$, $\delta f_{1}$ and $\delta\tau_{1}$ (dotted black lines), with
progressively smaller uncertainties.

As with most other parameters, our ability to constrain $\delta f_{1}$ and
$\delta\tau_{1}$ improves with SNR. We demonstrate this in Fig. 25 with the
same injection as before but scaled in overall amplitude to produce
progressively higher network SNRs, with the injected matched-filter SNR as
defined in Eq. (54). The $\delta\tau_{1}$ measurement improves slowly (bottom
right panel), with support significantly restricted to $\delta\tau_{1}\leq 1$
only for SNR $\geq 50$. The tightening of $\delta f_{1}$ and $\delta\tau_{1}$
posteriors is accompanied by an improvement in the $M$ and $\chi$ measurement,
for which the posteriors converge to the right values. The individual mode
amplitudes and phases are also better determined for higher SNR (not shown).

In the examples above, the model used to produce the injection matches that of
the recovery. In realistic situations, however, higher SNRs will tend to
reveal additional signal features that are not explicitly captured by our QNM
model; it is thus important to understand at which point the induced
systematic error may mimic a deviation from Kerr. In the context of BH
overtones, for example, we may expect to see contributions from an increasing
number of tones with $n>1$ as we analyze data closer to the signal peak [38,
65, 33, 100, 70]. When fitting an $N=1$ model, such additional modes could
conceivably bias the $\delta f_{1}$ or $\delta\tau_{1}$ measurements. The
magnitude of this systematic bias was quantified as a function of SNR in [38],
both through a heuristic calculation based on template mismatches, and through
simulated Bayesian measurements of $M$ and $\chi$. Here, we expand upon those
results to show that, indeed, contamination due to $n>1$ will only bias
$\delta f_{1}$ or $\delta\tau_{1}$ for significantly high SNRs—loud enough for
the additional modes to be detected by the procedure in Sec. IV.2.

We simulate the same Kerr signal as above (Figs. 22–25), but add to it a third
mode corresponding to the second $\ell=m=2$ overtone ($n=2$). Our choice of
$M$ and $\chi$ (Table 3) determines the frequency and damping time of this new
mode, but we have freedom to choose its amplitude and phase parameters. In
order to produce a somewhat realistic example, we set those quantities based
on an $N=2$ least-squares fit of a GW150914-like numerical relativity waveform
(SXS:BBH:0305 [101, 102, 66]), modeled starting at peak strain; this serves as
a proxy for features in a real signal that would not be captured by the $N=1$
ringdown model. Based on that fit, we set $A_{2}/A_{0}=1.1$,
$\phi_{2}=0.83\,\mathrm{rad}$ and $\epsilon_{2}=-1$ in our injection, but we
do not modify the $n=1$ parameters. Assuming only LIGO Hanford and Livingston
are operational, we analyze this new $N=2$ Kerr signal with the same $N=1$
model we had been using so far, progressively increasing the injected SNR.

Figure 26: Effect of not accounting for the second overtone in a spectroscopic
measurement assuming $N=1$ for different total SNRs (color). This is the same
analysis on the Kerr signal in Fig. 25, except for the injection of the $n=2$
mode, which is not accounted in the recovery model. This causes a systematic
bias in the recovered $\delta f_{1}$ but, for this configuration, this only
becomes significant at 90% credibility when the injected network SNR reaches
100.

Figure 26 shows the results in the $\delta f_{1}$, $\delta\tau_{1}$ plane. For
this configuration, the posterior density does not begin to show an obvious
systematic bias until the total injected SNR reaches ${\sim}50$. The presence
of the $n=2$ mode in the data manifests in an overall shift of the posterior
towards $\delta f_{1}\geq 0$, but the damping rate parameter $\delta\tau_{1}$
remains completely unaffected even for the highest SNR we consider. Even
though $\delta f_{1}$ is clearly sensitive to the $n=2$ mode, it is not until
we reach $\mathrm{SNR}\gtrsim 100$ that the Kerr value ($\delta
f_{1}=\delta\tau_{1}=0$) falls outside the 90%-credible region.

However, this SNR is sufficiently high that the presence of the second
overtone would have been revealed by the mode-detection procedure outlined in
Sec. IV.2. Indeed, even before the SNR reaches ${\sim}70$, the joint
$A_{1}{-}A_{2}$ posterior stops displaying the characteristic “arching” that
we expect if only two modes are identified in the data (Fig. 27). This
suggests that it would be unlikely for us to be misled by the bias seen at
high SNRs in Fig. 26 when looking at a signal like this. Nevertheless, we do
expect the kind of bias to compound in the batch analysis of multiple signals
(as in [44]), especially if extreme care is not taken to identify signals
which present more visible modes than expected.

Although, the results Fig. 27 provide some idea of the SNRs at which we expect
systematics to overcome statistical uncertainty, the actual level of bias will
vary with the details of the signal and noise instantiation. In practice, this
means that potential systematics should be studied assiduously for each case,
especially if there are hints of a Kerr deviation. This task would be aided by
a dedicated study of numerical relativity simulations, which we leave for
future work.

Figure 27: Amplitude posteriors from the $N=2$ analysis of the Kerr $N=2$
signal in Fig. 26, carried out assuming a Kerr spectrum and for different
total injection SNRs (color). Contours represent the 90%-credible regions for
each of the three amplitudes, measured as a fraction of the respective true
injected value, so that unity represents the truth (dotted lines). For SNRs 14
and 20, the joint $A_{1}{-}A_{2}$ posterior displays the arching seen in Fig.
12; this is not the case for higher SNRs, revealing the presence of the third
mode ($n=2$) and preventing us from being misled by the spurious Kerr
deviation in in Fig. 26.

####  V.2.2 Measurement of non-Kerr signals

Figure 28: Simulated $\delta f_{1}$ and $\delta\tau_{1}$ measurements for
different injected values (color) and SNRs (14, thin trace; 50, thick trace).
The signal is as in Fig. 22, except for the nonzero values of the deviation