# Simulation-based inference of black hole ringdowns in the time domain

Costantino Pacilio  0000-0002-8140-4992<EMAIL_ADDRESS>Dipartimento di Fisica “G. Occhialini”, Universitá degli Studi di Milano-
Bicocca, Piazza della Scienza 3, 20126 Milano, Italy  INFN, Sezione di Milano-
Bicocca, Piazza della Scienza 3, 20126 Milano, Italy  Swetha Bhagwat
0000-0003-4700-5274 School of Physics and Astronomy & Institute for
Gravitational Wave Astronomy, University of Birmingham,
Birmingham, B15 2TT, United Kingdom  Roberto Cotesta  William H. Miller III
Department of Physics and Astronomy, Johns Hopkins University, 3400 N. Charles
Street, Baltimore, Maryland, 21218, USA

###### Abstract

Gravitational waves emitted by a ringing black hole allow us to perform
precision tests of General Relativity in the strong field regime. With
improvements to our current gravitational wave detectors and upcoming next-
generation detectors, developing likelihood-free parameter inference
infrastructure is critical as we will face complications like non-standard
noise properties, partial data and incomplete signal modeling that may not
allow for an analytically tractable likelihood function. In this work, we
present a proof-of-concept strategy to perform likelihood-free Bayesian
inference on ringdown gravitational waves using simulation based inference.
Specifically, our method is based on truncated sequential neural posterior
estimation, which trains a neural density estimator of the posterior for a
specific observed data segment. We setup the ringdown parameter estimation
directly in the time domain. We show that the parameter estimation results
obtained using our trained networks are in agreement with well-established
Markov-Chain methods for simulated injections as well as analysis on real
detector data corresponding to GW150914. Additionally, to assess our
approach’s internal consistency, we show that the density estimators pass a
Bayesian coverage test.

_Introduction_ — The detection of gravitational waves (GWs) from binary black-
hole (BBH) mergers [1] enables the possibility to perform novel tests of
General Relativity (GR) in the strong-field regime [2, 3, 4]. GR predicts that
stationary black holes (BHs) are described by a strikingly simple two-
parameter family metric known as the Kerr solution [5], a result which is
commonly known as the “no-hair theorem” [6, 7]. This is particularly relevant
to the final phase of a BBH evolution, where the two black holes have merged
to form a perturbed remnant that settles down to a stationary Kerr BH [8].

Perturbed BHs shed away multiple moments by emitting GWs [9] with
characteristic complex oscillation frequencies called quasi-normal modes
(QNMs) [10, 11]. In GR, the oscillation frequency spectrum can be predicted
from the mass and spin of the remnant BH [12, 11] in accordance with the no-
hair theorem. Detecting multiple QNMs from the observed signal, also known as
black-hole spectroscopy [13, 14], allows us to test deviations from the BH
geometry and in the perturbative equation of motion in GR as well as the
boundary conditions imposed to solve the perturbative equation of motion [15,
16].

Detecting multiple QNMs unambiguously with present GW detectors has been
challenging thus far due to the low signal-to-noise ratio (SNR) in the
ringdown portion of the signal [3, 4, 14] (but see [17, 18] for recent hints
of multiple QNMs detections from GW data). On the other hand, it is expected
that multiple ringdowns at high SNRs will be observed with next-generation GW
detectors such as LISA [19], Einstein Telescope [20], and Cosmic Explorer
[21].

Next-generation detectors are expected to routinely observe ringdown SNRs as
high as $\mathcal{O}(100)$, with exceptionally loud events having SNRs greater
than $1000$ [22, 23, 24, 25, 26]. Such observations will enable to detect
several superimposed QNMs and measure deviations from GR at sub-percent level
[25, 26]. These will require more powerful Bayesian techniques to isolate the
signal from the background noise, to handle non-standard noise properties like
non-stationarities and non-Gaussianities [27, 28, 29], and gaps in data as
well as systematics in modelling the signal [30]. Hence, exploring alternative
inference methods is of fundamental interest.

The aim of this paper is to propose a likelihood-free method for the parameter
estimation of BH ringdowns, based on simulation-based inference (SBI) [31].
SBI methods have been previously applied to GW inference of a full inspiral-
merger-ringdown signals, starting with the pioneering works of [32, 33, 34,
35] and followed by manifold improvements [36, 37, 38, 39]. All these studies
show that SBI reaches an accuracy comparable with traditional Markov-Chain
methods for a full inspiral-merger-ringdown waveform [40, 41, 42, 43] and also
allow for a considerable speed-up in the analysis.

In addition, SBI is highly efficient in estimating automatically marginalized
posterior densities [39]. This is crucial in the parameter estimation of BHs
in the ringdown phase when the observable signal includes a very large number
of QNMs but, in practice, one is only interested is a small subset of well-
resolved modes for testing GR purposes [44, 45, 46, 47].

In this work we introduce a time domain implementation of likelihood-free
inference of BH ringdowns using SBI. In contrast to frequency domain
approaches presented in a recent study [48], this technique facilitates noise
simulation and requires fewer training samples. Likelihood-free ringdown
inference in the time domain was first explored in [49] based on the technique
of conditional variational auto-encoders [50, 51]. Here, we demonstrate for
the first time that SBI matches the performance of well-established Markov-
Chain methods for several simulated signal injections [52, 53]. We also
provide a real-data application of estimating the posterior mass and spin of
the final BH remnant from GW150914 [1], which agrees well with previous
Bayesian analyses [54].

_Sequential neural posterior estimation_ — We use sequential neural posterior
estimation (SNPE) to estimate posterior densities. SNPE was first introduced
in [55] and further expanded in [56, 57]. The aim of SNPE is to train an
approximate density estimator
$q_{\boldsymbol{\phi}}(\boldsymbol{\theta}|\boldsymbol{x})$ of the true
posterior density $p(\boldsymbol{\theta}|\boldsymbol{x})$. Here,
$\boldsymbol{x}$ is the data segment and $\boldsymbol{\theta}$ is the set of
model parameters to be estimated; in particular, $\boldsymbol{x}$ is the sum
of a deterministic component $\boldsymbol{h}(\boldsymbol{\theta})$ (hereafter
identified with the strain of the GW projected on the detector) and a
stochastic component $\boldsymbol{n}$ (the noise of the detector),
$\boldsymbol{x}(\boldsymbol{\theta})=\boldsymbol{h}(\boldsymbol{\theta})+\boldsymbol{n}$.
It is assumed that one can generate $\boldsymbol{h}(\boldsymbol{\theta})$
though a numerical process, so as to prepare a training set of ordered couples
$\\{\boldsymbol{\theta}_{i},\boldsymbol{x}_{i}\\}$ to train the density
estimator.

The training process maximises the likelihood
$\prod_{i}q_{\boldsymbol{\phi}}(\boldsymbol{\theta}_{i}|\boldsymbol{x}_{i})$
w.r.t. the model parameters $\boldsymbol{\phi}$. As shown in [55], the
estimator converges to the true posterior density asymptotically with the
dimension of the training set. A particular class of estimators called
normalizing flows [58] has become increasingly popular for SNPE. In the
following, we model the estimator as a neural spline flow (NSF), a spline-
based neural normalizing flow [59, 35]. We use the implementation from the sbi
package [60] — see also Appendix A for implementational details.

Once trained, the model $q_{\boldsymbol{\phi}}$ can be sampled from in a
fraction of a second, returning fast posterior samples
$\\{\boldsymbol{\theta}_{n}\\}$. Moreover, differently from Markov-Chain
methods, the model can also be evaluated and returns fast accurate estimates
of $p(\boldsymbol{\theta}|\boldsymbol{x})$.

The word sequential refers to the fact that the approximate density estimator
$q_{\boldsymbol{\phi}}(\boldsymbol{\theta}|\boldsymbol{x})$ targets a
particular observation $\boldsymbol{x}_{o}$, as opposed to amortized NPE which
targets the entire prior volume of the training set. Specifically, sequential
training proceeds through a series of adaptive steps, or rounds: at each
round, new training samples are drawn from a proposal distribution
$\tilde{p}(\boldsymbol{\theta})$, which is initially set to the prior
$p(\boldsymbol{\theta})$ and then it shrinks adaptively at each round to
prioritize samples with high posterior density
$q_{\boldsymbol{\phi}}(\boldsymbol{\theta}|\boldsymbol{x})$ conditioned to
$\boldsymbol{x}=\boldsymbol{x}_{o}$. While the original implementation of SNPE
in [55] proposes to draw the new samples directly from the posterior of the
previous round, we opt for a refined version called truncated SNPE (TSNPE)
[57], as we observed it to give more stable results. We briefly summarize it:

  1. 1.

During the first round $r=1$, the training samples
$\\{\boldsymbol{\theta}_{i},\boldsymbol{x}_{i}\\}_{r=1}$ are drawn from the
prior $p(\boldsymbol{\theta})$ and the model is optimized w.r.t.
$\\{\boldsymbol{\theta}_{i},\boldsymbol{x}_{i}\\}_{r=1}$.

  2. 2.

During the $k$-th round ($k>1$), new training samples are still drawn from the
prior $p(\boldsymbol{\theta})$, but they are rejected if they fall outside the
$1-\epsilon$ highest posterior density (HPD) region of
$q_{\boldsymbol{\phi}}(\boldsymbol{\theta}|\boldsymbol{x})$ from the previous
round, evaluated at $\boldsymbol{x}=\boldsymbol{x}_{o}$. The model from the
previous round $r=k-1$ is then further optimized w.r.t. to the collection of
samples from all rounds
$\\{\boldsymbol{\theta}_{i},\boldsymbol{x}_{i}\\}_{r=1}\cup\dots\cup\\{\boldsymbol{\theta}_{i},\boldsymbol{x}_{i}\\}_{r=k}$.

  3. 3.

The training proceeds until a stopping criterion is met.

Here, $\epsilon$ is an hyper-parameter which is fixed in advance and plays a
crucial role: a large $\epsilon$ leads to under-covered posteriors due to
large truncation of the original prior volume, while a small $\epsilon$ does
not adapt the prior enough to advance the training. Ref. [57] reports that the
range $10^{-5}<\epsilon<10^{-2}$ is working for a variety of benchmarking
problems. In this work, we fix $\epsilon=10^{-4}$ as we found it to give the
best results.

To avoid overfitting, we stop the training when the truncated volume at the
current round encompasses more than $80\%$ of the truncated volume at the
previous round. Moreover, at rounds $r>1$, we also discard simulations from
the first round since they are less informative.

More details about the data preparation and the density estimator are provided
in Appendix A.

_Ringdown waveform_ — The GW strain that we use for the inference has the form

$\boldsymbol{h}(\boldsymbol{\theta})=F_{+}h_{+}+F_{\times}h_{\times}$ (1)

where $F_{+}$ and $F_{\times}$ are the pattern functions of the detector [61],
depending on the sky position and relative orientation of the source, of the
polarization of the waveform and of the starting time of the ringdown. The
plus- and cross- waveform components are expressed as superpositions of damped
sinusoids as

$\displaystyle h_{+}=\sum_{l,m,n}\mathcal{A}_{lmn}e^{-\frac{(t-t_{\rm
start})}{\tau_{lmn}}}\cos(\Phi_{lmn})Y_{lm}^{+}(\iota)$ (2a) $\displaystyle
h_{\times}=\sum_{l,m,n}\mathcal{A}_{lmn}e^{-\frac{(t-t_{\rm
start})}{\tau_{lmn}}}\sin(\Phi_{lmn})Y_{lm}^{\times}(\iota)$ (2b)

where $l\geq 2$, $|m|\leq l$ and $n\geq 0$, and

$\Phi_{lmn}\equiv 2\pi f_{lmn}(t-t_{\rm start})+\phi_{lmn}\,.$ (3)

The discrete indices $(l,m,n)$ label the (complex) QNMs
$\tilde{\omega}_{lmn}=2\pi f_{lmn}+i/\tau_{lmn}$ of the remnant black hole,
with $f_{lmn}$ and $\tau_{lmn}$ being the $(l,m,n)$ frequency and damping
time, respectively. While the QNM indices $(l,m,n)$ span a countably infinite
set, numerical simulations show that only a finite subset is significantly
excited in the aftermath of a binary black-hole merger [13, 14, 44]. The
amplitudes $\mathcal{A}_{lmn}$ quantify the extent to which different QNMs are
excited [14, 22, 62, 44, 63, 64, 65].

The expressions (2a) assume non-precessing progenitors, implying the
equatorial symmetry

$\mathcal{A}_{l-mn}e^{i\phi_{l-mn}}=(-1)^{l}\mathcal{A}_{lmn}e^{-i\phi_{lmn}}\,.$
(4)

Together with the symmetry $\tilde{\omega}_{l-mn}=-\tilde{\omega}_{lmn}^{*}$
valid for the QNMs of Kerr BHs, it allows to reabsorb negative-$m$ modes into
positive-$m$ ones and to restrict the sums over $m>0$ only.

The plus- and cross- spherical harmonics are functions of the inclination
angle $\iota$ defined as

$Y_{lm}^{+,\times}(\iota)=\tensor{[}_{-2}]{Y}{{}_{lm}}(\iota,0)\pm(-1)^{l}\tensor{[}_{-2}]{Y}{{}_{l-m}}(\iota,0)$
(5)

in terms of the spin-weighted spherical harmonics
$\tensor{[}_{-2}]{Y}{{}_{lm}}$. Note that we set the azimuthal angle to zero
because it is degenerate with the QNM phases. We also take the the spin-
weighted spherical harmonics $\tensor{[}_{-2}]{Y}{{}_{lm}}$ to approximate the
spin-weighted spheroidal harmonics [66],
$\tensor{[}_{-2}]{Y}{{}_{lm}}\approx\tensor{[}_{-2}]{S}{{}_{lmn}}$.

For the parameter estimation, we follow [67, 68, 17] and keep sky position,
polarization, starting time and inclination as fixed parameters.

As a consequence, the space of independent model parameters is spanned by
$\boldsymbol{\theta}=\\{M,\chi\\}\cup\\{\mathcal{A}_{lmn},\phi_{lmn}\\}_{lmn}$
and has $2(1+N_{\text{modes}})$ dimensions. We assume that GR is valid and we
map $\\{M_{f},\chi_{f}\\}$ into the Kerr QNMs using the fits in [45].

_Simulated injections_ — We perform simulated injections into zero noise [69,
70] to benchmark the ability to recover a known set of injected parameters. We
simulate three systems:

  1. 1.

$\text{Kerr}_{220}$: A system containing only the $(2,2,0)$ mode. This system
is designed to benchmark our inference strategy on the simplest problem of a
single excited mode.

  2. 2.

$\text{Kerr}_{221}$: A system containing the $(2,2,0)$ mode and the $(2,2,1)$
mode, checking the ability to recover the fundamental tone and its first
overtone.

  3. 3.

$\text{Kerr}_{330}$: A system containing the $(2,2,0)$ and the $(3,3,0)$ mode,
checking the ability to recover the fundamental tone and a higher angular
mode.

Figure 1: Marginalized posteriors for $(M_{f},\chi_{f})$ for the three
injected systems. The systems $\text{Kerr}_{220}$ and $\text{Kerr}_{221}$ have
$\text{SNR}\approx 14$ while the system $\text{Kerr}_{330}$ has
$\text{SNR}\approx 53$, which justifies its narrower posteriors. 2-D contours
denote $68\%$ and $90\%$ credible regions. Black dashed lines denote the
injected parameters.

In all cases, the mass and spin of the remnant BH are fixed to
$M_{f}=67~{}M_{\odot}$ and $\chi_{f}=0.67$. The other parameters are listed in
Table 1. For the $\text{Kerr}_{220}$ and the $\text{Kerr}_{221}$ systems, we
fix the inclination angle to the GW150914-like value $\iota=\pi$ [67, 68].
Since the spherical harmonics $Y_{330}^{+,\times}$ vanish at $\iota=\pi$, we
fix $\iota=\pi/4$ for the $\text{Kerr}_{330}$ system, in order to excite the
higher angular mode.

Following [67], we fix the right ascension $\alpha$, declination $\delta$ and
polarization $\psi$ to the GW150914-like values
$\\{\alpha,\delta,\psi\\}=\\{1.95,-1.27,0.82\\}$. Similarly, we fix the
starting time of the ringdown to the median peak time estimated from the
posterior samples of GW150914 [71, 68]; specifically, we fix the GPS starting
time at H1 to be $t_{\rm H1}=1126259462.42323$, while the starting time at L1
is derived from the time delay from H1 using the sky location fixed above.

The SNRs are computed from the noise power spectral density (PSD), estimated
from the data segments in the vicinity of the event GW150914 [53]. The
amplitudes of $\text{Kerr}_{220}$ and the $\text{Kerr}_{221}$ systems are
adjusted to give a signal-to-noise ratio $\text{SNR}\approx 14$, comparable to
the SNRs of ringdown events observed so far, e.g., GW150914. On the other
hand, $\text{SNR}\approx 14$ is too low to allow for a detection of the
subdominant $(3,3,0)$ mode: following [14, 23], we expect that the higher mode
$(3,3,0)$ is resolved from the $(2,2,0)$ starting from a ringdown SNR of
$\mathcal{O}(50)$, for non-spinning progenitors with mass ratio $q\approx
1.5$. Therefore, we fix the corresponding amplitudes to
$\\{\mathcal{A}_{220},\mathcal{A}_{330}\\}=\\{30\times 10^{-21},3\times
10^{-21}\\}$, which results into $\text{SNR}\approx 53$.

$\text{Kerr}_{220}:~{}(2,2,0)$
---
$M_{f}$ | $\chi_{f}$ | $\mathcal{A}_{220}$ |  | $\phi_{220}$ |  | $\iota$ | SNR
$67~{}M_{\odot}$ | $0.67$ | $5$ |  | $1.047$ |  | $\pi$ | 14
$\text{Kerr}_{221}:~{}(2,2,0)+(2,2,1)$
$M_{f}$ | $\chi_{f}$ | $\mathcal{A}_{220}$ | $\mathcal{A}_{221}$ | $\phi_{220}$ | $\phi_{221}$ | $\iota$ | SNR
$67~{}M_{\odot}$ | $0.67$ | $8.92$ | $9.81$ | $1.047$ | $4.19$ | $\pi$ | 14
$\text{Kerr}_{330}:~{}(2,2,0)+(3,3,0)$
$M_{f}$ | $\chi_{f}$ | $\mathcal{A}_{220}$ | $\mathcal{A}_{330}$ | $\phi_{220}$ | $\phi_{330}$ | $\iota$ | SNR
$67~{}M_{\odot}$ | $0.67$ | $30$ | $3$ | $1.047$ | $5.014$ | $\pi/4$ | 53
Table 1: Parameters of the simulated zero-noise injections. Amplitudes are
expressed in units of $10^{-21}$. Figure 2: Violin plots of the recovered
amplitudes $\mathcal{A}_{lmn}$ for each injected model. We subtract the
injected values $\mathcal{A}_{lmn}^{\rm(true)}$ for better visualization. All
amplitudes are expressed in units of $10^{-21}$.

We infer the posteriors from box-uniform priors in the model parameters,
within the prior ranges listed in Table 2. We benchmark our results against
posterior samples obtained with the time-domain inference software pyRing
(v2.3.0) [72, 52, 67, 3]. We sample with the nested sampler cpnest [73], using
4096 live points and 4094 maximum Markov-Chain steps, which typically results
in $\sim 20$k posterior samples. For consistency, we extract 20k posterior
samples from the trained estimator
$q_{\boldsymbol{\phi}}(\boldsymbol{\theta}|\boldsymbol{x}_{o})$ when doing the
comparisons.

parameter | prior range
---|---
$M_{f}$ | $[20,300]~{}M_{\odot}$
$\chi_{f}$ | $[0,0.99]$
$\mathcal{A}_{220}$ | $[0.1,50]\times 10^{-21}$
$\mathcal{A}_{221}$ | $[0,50]\times 10^{-21}$
$\mathcal{A}_{330}$ | $[0,50]\times 10^{-21}$
$\phi_{lmn}$ | $[0,2\pi]$
Table 2: Prior ranges for inferring the model parameters. All priors are box-
uniform distributions within the specified ranges.

We train the posterior density estimators on a single GPU A100. Individual
training processes take $\\{291,324,249\\}$ epochs and last $\\{81,70,52\\}$
mins for the $\\{\text{Kerr}_{220},\text{Kerr}_{221},\text{Kerr}_{330}\\}$
models respectively. Figure 1 displays the marginalized posteriors for
$(M_{f},\chi_{f})$ obtained with pyRing and with our SBI approach, showing an
excellent agreement between the two methods. Figure 2 shows that amplitudes,
containing the physical information on the degree of excitation of the QNMs,
are recovered consistently with Markov-Chain methods for all injected systems.
We also display full corner plots for the $\text{Kerr}_{221}$ and
$\text{Kerr}_{330}$ in Appendix C.

As an internal diagnostics, we perform coverage tests [74, 75, 76]. In
particular we check that, over a set of simulated injections from the
restricted priors $\tilde{p}(\boldsymbol{\theta})$, the Bayesian credible
intervals can be used similar to frequentist confidence regions; equivalently,
that the true injected values fall within the $\gamma\%$ credible regions in a
fraction $\gamma$ of the injections. This is a necessary condition for the
density estimators to provide a consistent parameter estimation. Further
details are provided in Appendix B. Figure 3 displays the cumulative
distribution ${\rm c.d.f.}(\gamma)$ of the Bayesian credible intervals for the
model ${\rm Kerr}_{220}$: in order to check consistency with the identity, we
evaluate $\text{c.d.f.}(\gamma)$ for 100 catalogs, each consisting of 100
events. We see that ${\rm c.d.f.}(\gamma)$ has support well within the
expected $90\%$ credible interval around the diagonal line, represented by the
gray region [4]111Note that such tests are difficult with ordinary Markov-
Chain methods, in which typically the densities are sampled from but not
evaluated exactly, often restricting the scope of coverage diagnostics to the
marginalized 1-d distributions [77].. We obtain similar consistency results
for the ${\rm Kerr}_{221}$ and ${\rm Kerr}_{330}$ models, as displayed in
Appendix B.

Figure 3: Cumulative distribution of the coverage
$\gamma(\boldsymbol{\theta}_{*},\boldsymbol{x}_{*})$ for the injected model
${\rm Kerr}_{220}$. Each of the 100 blue lines corresponds to the c.d.f. of
$\gamma$ from $N_{s}=100$ draws of $\boldsymbol{\theta}_{*}$. The shaded grey
area denote the $90\%$ uncertainty over $\text{c.d.f.}(\gamma)$. We also quote
the median and $90\%$ confidence bounds of the KS test p-values across all
draws.

Our results show that the time-domain SBI inference for ringdown returns
posterior parameter predictions that are internally consistent and it
reproduces the results obtained by the established Markov-Chain methods.

_Inference of real data_ — We perform a Bayesian parameter estimation on the
real event GW150914 [1]. A full inspiral-merger-ringdown (IMR) analysis
identified it as a binary BH with detector-frame total mass
$M=72^{+4}_{-3}~{}M_{\odot}$ and binary mass ratio $q<1.42$ at $90\%$
credibility, while being compatible with vanishing progenitor spins [71]. In
order to reconstruct the remnant properties, we analyze the event at $t_{\rm
H1}+3~{}{\rm ms}$, corresponding to $\sim 10M$ after the merger ($t_{\rm
start}\sim t_{\rm peak}+10M$ was found to mark the onset of the ringdown stage
from the analysis of nearly-equal mass NR simulations [44, 78]). This choice
results in a ringdown-only ${\rm SNR}\sim 8.5$ [54]. As in [67], we fix
$\\{\alpha,\delta,\phi\\}=\\{1.95,-1.27,0.82\\}$ and $\iota=\pi$. We assume a
model containing only the fundamental mode $(2,2,0)$, which dominates the late
ringdown stage for nearly equal-mass non-spinning binaries [79, 62, 22, 44].
The training takes 282 epochs and lasts 91 minutes. Figure 4 shows the $90\%$
credible regions of the reconstructed final mass and final spin using our SBI
implementation, compared to the estimates from the pyRing software. The two
results are in good agreement with each other. In the plot, we also display
the final mass and spin deduced from the IMR posterior samples [71] using the
phenomenological fits in [80, 81].

Figure 4: Remnant properties of GW150914, obtained by analysing the publicly
released data from LVK [82]. Solid lines denote the $90\%$ credible regions
reconstructed with our SBI implementation (blue) and with the pyRing analysis
(orange). We also display the final mass and spin obtained from the full IMR
posterior samples [71] through the phenomenological remnant fits in [80, 81].

_Conclusions and outlook_ — Our study shows that simulation-based GW inference
of BBH ringdowns gives results that are both internally consistent and
compatible with Markov-Chain methods. After benchmarking SBI on a control set
of injections at zero noise, we applied our methods to the analysis of real GW
data from LIGO detectors. We see that inferences done using our model on the
event GW150914 agrees with the inference of the traditional setup. It is
important to note that zero-noise injections include several simplifications,
but real data present noise and a morphological bias between the template and
the underlying signal. The fact that our results are consistent in both cases
constitutes a robust demonstration of the feasibility of our approach.

We envisage a number of future studies to further expand our results. First, a
detailed investigation of all the promising ringdown events from the GW
observations is ongoing. While the significance for detecting a second mode in
present ringdown data is still debated [83, 3, 4, 68, 84, 48], there is
tentative evidence that GW150914 and GW190521 [85] allow us to measure the
presence of a sub-dominant overtone [67, 86] and of a higher angular harmonic
[17, 18], respectively. We seek to confirm or refute these findings with our
independent, likelihood-free method.

A promising application leverages on the ability of SBI to automatically
marginalize over nuisance parameters. In particular, when inferring on
ringdown signals, one is only interested in the small subset of modes which
are better resolved. The presence of multiple modes, which are challenging to
resolve individually, is a potential source of confusion noise and therefore
it must be modelled. Moreover, an ideal model would allow for a variable
number of such modes, in order to fit the residuals with greater flexibility.
This problem is especially suited for SBI, as a variable number of background
modes can be regarded as a collection of nuisances, and we will present a
dedicated exploration in a follow-up work.

Lastly, we use sequential (i.e., non-amortised) NPE because it is convenient
for analyzing individual observations, which is the scope of our work. On the
other hand, amortised NPE might prove useful when assessing the expected
performances of a detector on large population catalogs. Such studies are
well-timed, given the recent adoption of the LISA mission [19] by the European
Space Agency, and the strong efforts to build next-generation ground-based GW
detectors [20, 21]. Relevant examples include [87] exploring a battery of
diverse science cases, and [23, 25, 26, 88, 89] with a focus on BH
spectroscopy. Since performing rigorous Bayesian analysis on tens of thousands
of sources is computationally challenging, all these works approximate
individual posterior densities with Fisher matrices, which are known to
introduce systematic biases in the shape of the posterior recovery [90].
Training an amortised density estimator on the expected noise spectral
densities would allow to drop such approximations and to draw much more robust
conclusions.

_Acknowledgments_ — We are especially grateful to Enrico Barausse and Marco
Crisostomi for fruitful exchanges during an earlier stage of this work. CP
thanks Elliot Finch for useful discussions. C.P. is supported by ERC Starting
Grant No. 945155–GWmining, Cariplo Foundation Grant No. 2021-0555, MUR PRIN
Grant No. 2022-Z9X4XS, and the ICSC National Research Centre funded by
NextGenerationEU. CP acknowledges support from the European Union’s H2020 ERC
Consolidator Grant “GRavity from Astrophysical to Microscopic Scales” (Grant
No. GRAMS-815673) during a visit to SISSA. S.B. is supported by UKRI Stephen
Hawking Fellowship No. EP/W005727. R.C. is supported by NSF Grants No.
AST-2006538, PHY-2207502, PHY-090003 and PHY-20043, and by NASA Grants No.
20-LPS20-0011 and 21-ATP21-0010. Computational work was performed on the
LEONARDO HPC facility hosted at CINECA, with allocations through INFN and
Bicocca.

#### Software.

We used the following public software: scipy [91] and pycbc [41] to estimate
the PSD and the correlation matrix from noisy time series; sbi [60] for
simulation based inference; pyRing [72] for Markov-Chain ringdown inferences;
numpy [92], matplotlib [93] and seaborn [94] for chain visualizations.

##  Appendix A Technical details

In this Appendix, we provide technical details on the process of data
preparation, and on the numerical implementation and training of the neural
density estimator.

###  A.1 Data preparation

When performing simulated injections, we sample the strains at $2048$ Hz and
truncate them at a duration of $0.1$ s, thus resulting in a data segment of
$204$ bins per detector. When analyzing real data, we use the publicly
available data from GWOSC [82] at $4096$ Hz and down-sample them to $2048$ Hz.

We whiten the data directly in the time domain following [53],

$h_{\rm
white}^{a}(\boldsymbol{\theta})=\sum_{b}\tensor{(L^{-1})}{{}^{a}_{b}}h^{b}(\boldsymbol{\theta})$
(6)

where the matrix $L$ is the (upper-triangular) Cholesky factor of the
covariance matrix. Specifically, we first compute the auto-correlation
function (ACF) from the inverse-Fourier transform of the PSD, then, we compute
the covariance matrix as the Toeplitz matrix of the ACF (see Eq.s (39)-(44) of
[53]). Noise realizations can be added to the whitened ringdown strain as
simply as

$\boldsymbol{x}_{\rm white}(\boldsymbol{\theta})=\boldsymbol{h}_{\rm
white}(\boldsymbol{\theta})+\mathcal{N}(0,1)\,.$ (7)

###  A.2 Density estimator

We model the density estimator as a neural spline flow (NSF) [58, 35] and use
the implementation from the sbi package [60]. The relevant hyper-parameters
are listed in Table 3. In particular, the NSF is a flow of 5 transforms and
each hidden layer contains 150 units.

The NSF does not take $\boldsymbol{x}_{\rm white}$ in input directly, but we
apply an embedding network for dimensional reduction: the raw input is a
concatenation of $\boldsymbol{x}_{\rm white}$ from two detectors (LIGO-Hanford
and LIGO-Livingstone), thus resulting into an input dimension of $204+204$
bins; they are then mapped into 128 bins by a fully connected neural network
with two hidden layers consisting of 150 units. The output of the embedding
network is fed into the NSF.

We train and validate the density estimator in batches of 512 samples. During
the first round we generate 50k training samples, while at later rounds we
generate 100k additional samples. We use the Adam optimizer with learning rate
0.001.

Since neural networks learn better from standardized and/or normalized data
[95], we linearly rescale $\boldsymbol{x}_{\rm white}$ to have zero mean and
unit variance along the feature dimension [96], and we normalize the
parameters $\boldsymbol{\theta}$ between 0 and 1.

Finally, we vary the noise realizations in (7) at each training epoch to make
the inference resilient to specific noise realizations222To avoid confusion,
we clarify the distinction between training rounds and training epochs. A
training round is the training update of the density estimator at each
truncation of the prior volume. A training epoch is a single feed-forward and
back-propagation step of the training set into the network. Each training
round consists of several-to-many training epochs, necessary to optimize the
model at the current round.. Note that the original implementation of sbi does
not allow to vary noise realizations between training epochs. We leverage on
the flexibility of object-oriented programming to wrap the original SNPE
implementation of sbi and we redefine its loss method, so as to resample the
noise at each evaluation of the loss333Such a data-augmentation strategy is
feasible, because sampling a standard normal in (7) is a fast numerical
operation and it does not impact on the training time. In cases where noise
generation is computationally costly, one would opt for alternative data
augmentations, e.g., preparing the training set offline with multiple copies
of the same raw strains $\boldsymbol{h}(\boldsymbol{\theta})$ but different
noise realizations. See [48] for an example of the last strategy..

Neural spline flow
---
num_blocks | 2
hidden_features | 150
num_transforms | 5
num_bins | 10
batch_norm | True
Embedding FC network
input_dim | 408
num_hidden_layers | 2
hidden_dim | 150
output_dim | 128
Training hyper-parameters
num_simulations | [50k,100k,…]
batch_size | 512
learning_rate | 0.001
validation_fraction | 0.1
trunc_quantile $\epsilon$ | $10^{-4}$
stopping_ratio | 0.8
varying_noise | True
Table 3: Hyper-parameters for the architecture and training of the the density
estimator.

##  Appendix B Coverage Test

Given model parameters $\boldsymbol{\theta}_{*}$ and a corresponding model
realization
$\boldsymbol{x}_{*}=\boldsymbol{h}(\boldsymbol{\theta}_{*})+\boldsymbol{n}$,
we define the coverage $\gamma(\boldsymbol{\theta}_{*},\boldsymbol{x}_{*})$ as
the approximate posterior probability contained within the highest posterior
density (HPD) region which has $\boldsymbol{\theta}_{*}$ at its boundary,

$\gamma(\boldsymbol{\theta}_{*},\boldsymbol{x}_{*})=\int
d\boldsymbol{\theta}~{}q_{\boldsymbol{\phi}}(\boldsymbol{\theta}|\boldsymbol{x}_{*})\mathbbm{1}\left[q_{\boldsymbol{\phi}}(\boldsymbol{\theta}|\boldsymbol{x}_{*})>q_{\boldsymbol{\phi}}(\boldsymbol{\theta_{*}}|\boldsymbol{x}_{*})\right]$
(8)

where $q_{\boldsymbol{\phi}}$ is the density estimator.

The integral (8) can be evaluated efficiently via importance sampling as the
expectation value

$\gamma(\boldsymbol{\theta}_{*},\boldsymbol{x}_{*})\approx\mathbb{E}_{\boldsymbol{\theta}\sim
q_{\boldsymbol{\phi}(\boldsymbol{\theta}|\boldsymbol{x}_{*})}}\big{(}\mathbbm{1}\left[q_{\boldsymbol{\phi}}(\boldsymbol{\theta}|\boldsymbol{x}_{*})>q_{\boldsymbol{\phi}}(\boldsymbol{\theta_{*}}|\boldsymbol{x}_{*})\right]\big{)}$
(9)

thanks to the fact that the neural density estimator is fast to sample and to
evaluate. It can be shown [75, 76] that, if the density estimator
$q_{\boldsymbol{\phi}}(\boldsymbol{\theta}|\boldsymbol{x})$ approximates the
true posterior $p(\boldsymbol{\theta}|\boldsymbol{x})$, then the coverage
across the prior $p(\boldsymbol{\theta})$ is distributed uniformly within
$[0,1]$. The latter statement is equivalent to
$\text{c.d.f.}(\gamma)\equiv\int_{0}^{\gamma}d\gamma^{\prime}p(\gamma^{\prime})=\gamma$
for each $\gamma\equiv\gamma(\boldsymbol{\theta}_{*},\boldsymbol{x}_{*})$ and
$\boldsymbol{\theta}_{*}\sim p(\boldsymbol{\theta})$. This is a necessary
condition for the density estimator to be a valid inference and we check for
it.

Since we are performing truncated SNPE, we cannot sample the injected
parameters from the original prior, $\boldsymbol{\theta}_{*}\sim
p(\boldsymbol{\theta})$. Indeed, sequential training progressively truncates
the prior volume as the training proceeds and the density estimator is only
guaranteed to work on the final version of the truncated prior
$\tilde{p}(\boldsymbol{\theta})$. Therefore, in the case of TNSPE, the prior
for the Bayesian coverage test must correspond to the final truncated prior,
$\boldsymbol{\theta}_{*}\sim\tilde{p}(\boldsymbol{\theta})$.

Figure 5 displays coverage diagnostics for the the injected models
corresponding to ${\rm Kerr}_{221}$ and ${\rm Kerr}_{330}$, complementing the
diagnostics for the ${\rm Kerr}_{220}$ presented in Figure 3 in the main text.

Figure 5: Cumulative distributions of the coverage
$\gamma(\boldsymbol{\theta}_{*},\boldsymbol{x}_{*})$ for the injected models
${\rm Kerr}_{221}$ and ${\rm Kerr}_{330}$. Each blue line corresponds to the
c.d.f. of $\gamma$ from $N_{s}=100$ draws of $\boldsymbol{\theta}_{*}$. Each
panel plots the cumulatives from $100$ experiments. Shaded grey areas denote
the $90\%$ uncertainties over $\text{c.d.f.}(\gamma)$. For each panel, we also
quote the median and $90\%$ confidence bounds of the KS test p-values across
all draws.

##  Appendix C Additional corner plots

Figure 6: Full corner plot of the posterior recovered from the injected system
$\text{Kerr}_{221}$. 2-D contours denote the $68\%$ and $90\%$ credible
regions. Black dashed lines denote the injected parameters from Table 1.
Figure 7: Full corner plot of the posterior recovered from the injected system
$\text{Kerr}_{330}$. 2-D contours denote the $68\%$ and $90\%$ credible
regions. Black dashed lines denote the injected parameters from Table 1.

## References

  * Abbott  _et al._ [2016a] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 061102 (2016a), arXiv:1602.03837 [gr-qc] .
  * Yunes  _et al._ [2016] N. Yunes, K. Yagi, and F. Pretorius, Phys. Rev. D 94, 084002 (2016), arXiv:1603.08955 [gr-qc] .
  * Abbott  _et al._ [2021a] R. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. D 103, 122002 (2021a), arXiv:2010.14529 [gr-qc] .
  * Abbott  _et al._ [2021b] R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), (2021b), arXiv:2112.06861 [gr-qc] .
  * Kerr [1963] R. P. Kerr, Phys. Rev. Lett. 11, 237 (1963).
  * Israel [1968] W. Israel, Commun. Math. Phys. 8, 245 (1968).
  * Carter [1971] B. Carter, Phys. Rev. Lett. 26, 331 (1971).
  * Teukolsky [2015] S. A. Teukolsky, Class. Quant. Grav. 32, 124006 (2015), arXiv:1410.2130 [gr-qc] .
  * Vishveshwara [1970] C. V. Vishveshwara, Nature 227, 936 (1970).
  * Kokkotas and Schmidt [1999] K. D. Kokkotas and B. G. Schmidt, Living Rev. Rel. 2, 2 (1999), arXiv:gr-qc/9909058 .
  * Berti  _et al._ [2009] E. Berti, V. Cardoso, and A. O. Starinets, Class. Quant. Grav. 26, 163001 (2009), arXiv:0905.2975 [gr-qc] .
  * Teukolsky [1972] S. A. Teukolsky, Phys. Rev. Lett. 29, 1114 (1972).
  * Berti  _et al._ [2006] E. Berti, V. Cardoso, and C. M. Will, Phys. Rev. D 73, 064030 (2006), arXiv:gr-qc/0512160 .
  * Berti  _et al._ [2007] E. Berti, J. Cardoso, V. Cardoso, and M. Cavaglia, Phys. Rev. D 76, 104044 (2007), arXiv:0707.1202 [gr-qc] .
  * Dreyer  _et al._ [2004] O. Dreyer, B. J. Kelly, B. Krishnan, L. S. Finn, D. Garrison, and R. Lopez-Aleman, Class. Quant. Grav. 21, 787 (2004), arXiv:gr-qc/0309007 .
  * Barausse and Sotiriou [2008] E. Barausse and T. P. Sotiriou, Phys. Rev. Lett. 101, 099001 (2008), arXiv:0803.3433 [gr-qc] .
  * Capano  _et al._ [2023] C. D. Capano, M. Cabero, J. Westerweck, J. Abedi, S. Kastha, A. H. Nitz, Y.-F. Wang, A. B. Nielsen, and B. Krishnan, Phys. Rev. Lett. 131, 221402 (2023), arXiv:2105.05238 [gr-qc] .
  * Siegel  _et al._ [2023] H. Siegel, M. Isi, and W. M. Farr, Phys. Rev. D 108, 064008 (2023), arXiv:2307.11975 [gr-qc] .
  * Amaro-Seoane  _et al._ [2017] P. Amaro-Seoane, H. Audley, S. Babak, J. Baker, E. Barausse, P. Bender, E. Berti, P. Binetruy, M. Born, D. Bortoluzzi, J. Camp, C. Caprini, V. Cardoso, M. Colpi, J. Conklin, N. Cornish, C. Cutler, K. Danzmann, R. Dolesi, L. Ferraioli, V. Ferroni, E. Fitzsimons, J. Gair, L. Gesa Bote, D. Giardini, F. Gibert, C. Grimani, H. Halloin, G. Heinzel, T. Hertog, M. Hewitson, K. Holley-Bockelmann, D. Hollington, M. Hueller, H. Inchauspe, P. Jetzer, N. Karnesis, C. Killow, A. Klein, B. Klipstein, N. Korsakova, S. L. Larson, J. Livas, I. Lloro, N. Man, D. Mance, J. Martino, I. Mateos, K. McKenzie, S. T. McWilliams, C. Miller, G. Mueller, G. Nardini, G. Nelemans, M. Nofrarias, A. Petiteau, P. Pivato, E. Plagnol, E. Porter, J. Reiche, D. Robertson, N. Robertson, E. Rossi, G. Russano, B. Schutz, A. Sesana, D. Shoemaker, J. Slutsky, C. F. Sopuerta, T. Sumner, N. Tamanini, I. Thorpe, M. Troebs, M. Vallisneri, A. Vecchio, D. Vetrugno, S. Vitale, M. Volonteri, G. Wanner, H. Ward, P. Wass, W. Weber, J. Ziemer, and P. Zweifel, arXiv e-prints , arXiv:1702.00786 (2017), arXiv:1702.00786 [astro-ph.IM] .
  * Punturo  _et al._ [2010] M. Punturo _et al._ , Class. Quant. Grav. 27, 194002 (2010).
  * Abbott  _et al._ [2017] B. P. Abbott _et al._ (LIGO Scientific), Class. Quant. Grav. 34, 044001 (2017), arXiv:1607.08697 [astro-ph.IM] .
  * Gossan  _et al._ [2012] S. Gossan, J. Veitch, and B. S. Sathyaprakash, Phys. Rev. D 85, 124056 (2012), arXiv:1111.5819 [gr-qc] .
  * Berti  _et al._ [2016] E. Berti, A. Sesana, E. Barausse, V. Cardoso, and K. Belczynski, Phys. Rev. Lett. 117, 101102 (2016), arXiv:1605.09286 [gr-qc] .
  * Cabero  _et al._ [2020] M. Cabero, J. Westerweck, C. D. Capano, S. Kumar, A. B. Nielsen, and B. Krishnan, Phys. Rev. D 101, 064044 (2020), arXiv:1911.01361 [gr-qc] .
  * Bhagwat  _et al._ [2022] S. Bhagwat, C. Pacilio, E. Barausse, and P. Pani, Phys. Rev. D 105, 124063 (2022), arXiv:2201.00023 [gr-qc] .
  * Bhagwat  _et al._ [2023] S. Bhagwat, C. Pacilio, P. Pani, and M. Mapelli, Phys. Rev. D 108, 043019 (2023), arXiv:2304.02283 [gr-qc] .
  * Kumar  _et al._ [2022] S. Kumar, A. H. Nitz, and X. J. Forteza, (2022), arXiv:2202.12762 [astro-ph.IM] .
  * Baghi  _et al._ [2022] Q. Baghi, N. Korsakova, J. Slutsky, E. Castelli, N. Karnesis, and J.-B. Bayle, Phys. Rev. D 105, 042002 (2022), arXiv:2112.07490 [gr-qc] .
  * Spadaro  _et al._ [2023] A. Spadaro, R. Buscicchio, D. Vetrugno, A. Klein, D. Gerosa, S. Vitale, R. Dolesi, W. J. Weber, and M. Colpi, Phys. Rev. D 108, 123029 (2023), arXiv:2306.03923 [gr-qc] .
  * Toubiana  _et al._ [2023] A. Toubiana, L. Pompili, A. Buonanno, J. R. Gair, and M. L. Katz, (2023), arXiv:2307.15086 [gr-qc] .
  * Cranmer  _et al._ [2020] K. Cranmer, J. Brehmer, and G. Louppe, Proceedings of the National Academy of Science 117, 30055 (2020), arXiv:1911.01429 [stat.ML] .
  * Chua and Vallisneri [2020] A. J. K. Chua and M. Vallisneri, Phys. Rev. Lett. 124, 041102 (2020), arXiv:1909.05966 [gr-qc] .
  * Gabbard  _et al._ [2022] H. Gabbard, C. Messenger, I. S. Heng, F. Tonolini, and R. Murray-Smith, Nature Phys. 18, 112 (2022), arXiv:1909.06296 [astro-ph.IM] .
  * Green  _et al._ [2020] S. R. Green, C. Simpson, and J. Gair, Phys. Rev. D 102, 104057 (2020), arXiv:2002.07656 [astro-ph.IM] .
  * Green and Gair [2021] S. R. Green and J. Gair, Mach. Learn. Sci. Tech. 2, 03LT01 (2021), arXiv:2008.03312 [astro-ph.IM] .
  * Dax  _et al._ [2021] M. Dax, S. R. Green, J. Gair, J. H. Macke, A. Buonanno, and B. Schölkopf, Phys. Rev. Lett. 127, 241103 (2021), arXiv:2106.12594 [gr-qc] .
  * Dax  _et al._ [2023] M. Dax, S. R. Green, J. Gair, M. Pürrer, J. Wildberger, J. H. Macke, A. Buonanno, and B. Schölkopf, Phys. Rev. Lett. 130, 171403 (2023), arXiv:2210.05686 [gr-qc] .
  * Wildberger  _et al._ [2023] J. Wildberger, M. Dax, S. R. Green, J. Gair, M. Pürrer, J. H. Macke, A. Buonanno, and B. Schölkopf, Phys. Rev. D 107, 084046 (2023), arXiv:2211.08801 [gr-qc] .
  * Bhardwaj  _et al._ [2023] U. Bhardwaj, J. Alvey, B. K. Miller, S. Nissanke, and C. Weniger, Phys. Rev. D 108, 042004 (2023), arXiv:2304.02035 [gr-qc] .
  * Veitch  _et al._ [2015] J. Veitch _et al._ , Phys. Rev. D 91, 042003 (2015), arXiv:1409.7215 [gr-qc] .
  * Biwer  _et al._ [2019] C. M. Biwer, C. D. Capano, S. De, M. Cabero, D. A. Brown, A. H. Nitz, and V. Raymond, Publ. Astron. Soc. Pac. 131, 024503 (2019), arXiv:1807.10312 [astro-ph.IM] .
  * Ashton and Talbot [2021] G. Ashton and C. Talbot, Mon. Not. Roy. Astron. Soc. 507, 2037 (2021), arXiv:2106.08730 [gr-qc] .
  * Ashton  _et al._ [2022] G. Ashton _et al._ , Nature 2, 10.1038/s43586-022-00121-x (2022), arXiv:2205.15570 [stat.CO] .
  * London  _et al._ [2014] L. London, D. Shoemaker, and J. Healy, Phys. Rev. D 90, 124032 (2014), [Erratum: Phys.Rev.D 94, 069902 (2016)], arXiv:1404.3197 [gr-qc] .
  * London and Fauchon-Jones [2019] L. London and E. Fauchon-Jones, Class. Quant. Grav. 36, 235015 (2019), arXiv:1810.03550 [gr-qc] .
  * Baibhav  _et al._ [2023] V. Baibhav, M. H.-Y. Cheung, E. Berti, V. Cardoso, G. Carullo, R. Cotesta, W. Del Pozzo, and F. Duque, Phys. Rev. D 108, 104020 (2023), arXiv:2302.03050 [gr-qc] .
  * Cheung  _et al._ [2024] M. H.-Y. Cheung, E. Berti, V. Baibhav, and R. Cotesta, Phys. Rev. D 109, 044069 (2024), arXiv:2310.04489 [gr-qc] .
  * Crisostomi  _et al._ [2023] M. Crisostomi, K. Dey, E. Barausse, and R. Trotta, Phys. Rev. D 108, 044029 (2023), arXiv:2305.18528 [gr-qc] .
  * Bhagwat and Pacilio [2021] S. Bhagwat and C. Pacilio, Phys. Rev. D 104, 024030 (2021), arXiv:2101.07817 [gr-qc] .
  * Kingma and Welling [2013] D. P. Kingma and M. Welling, arXiv e-prints , arXiv:1312.6114 (2013), arXiv:1312.6114 [stat.ML] .
  * Doersch [2016] C. Doersch, arXiv e-prints , arXiv:1606.05908 (2016), arXiv:1606.05908 [stat.ML] .
  * Carullo  _et al._ [2019] G. Carullo, W. Del Pozzo, and J. Veitch, Phys. Rev. D 99, 123029 (2019), [Erratum: Phys.Rev.D 100, 089903 (2019)], arXiv:1902.07527 [gr-qc] .
  * Isi and Farr [2021] M. Isi and W. M. Farr, (2021), arXiv:2107.05609 [gr-qc] .
  * Abbott  _et al._ [2016b] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 221101 (2016b), [Erratum: Phys.Rev.Lett. 121, 129902 (2018)], arXiv:1602.03841 [gr-qc] .
  * Papamakarios and Murray [2016] G. Papamakarios and I. Murray, Advances in neural information processing systems 29 (2016), arXiv:1605.06376 [stat.ML] .
  * Greenberg  _et al._ [2019] D. Greenberg, M. Nonnenmacher, and J. Macke, in  _International Conference on Machine Learning_ (PMLR, 2019) pp. 2404–2414, arXiv:1905.07488 [cs.LG] .
  * Deistler  _et al._ [2022] M. Deistler, P. J. Goncalves, and J. H. Macke, Advances in Neural Information Processing Systems 35, 23135 (2022), arXiv:2210.04815 [stat.ML] .
  * Papamakarios  _et al._ [2021] G. Papamakarios, E. Nalisnick, D. J. Rezende, S. Mohamed, and B. Lakshminarayanan, Normalizing flows for probabilistic modeling and inference (2021), arXiv:1912.02762 [stat.ML] .
  * Durkan  _et al._ [2019] C. Durkan, A. Bekasov, I. Murray, and G. Papamakarios, arXiv e-prints , arXiv:1906.04032 (2019), arXiv:1906.04032 [stat.ML] .
  * Tejero-Cantero  _et al._ [2020] A. Tejero-Cantero, J. Boelts, M. Deistler, J.-M. Lueckmann, C. Durkan, P. J. Gonçalves, D. S. Greenberg, and J. H. Macke, Journal of Open Source Software 5, 2505 (2020).
  * Jaranowski  _et al._ [1998] P. Jaranowski, A. Krolak, and B. F. Schutz, Phys. Rev. D 58, 063001 (1998), arXiv:gr-qc/9804014 .
  * Kamaretsos  _et al._ [2012a] I. Kamaretsos, M. Hannam, and B. Sathyaprakash, Phys. Rev. Lett. 109, 141102 (2012a), arXiv:1207.0399 [gr-qc] .
  * London [2020] L. T. London, Phys. Rev. D 102, 084052 (2020), arXiv:1801.08208 [gr-qc] .
  * Jiménez Forteza  _et al._ [2020] X. Jiménez Forteza, S. Bhagwat, P. Pani, and V. Ferrari, Phys. Rev. D 102, 044053 (2020), arXiv:2005.03260 [gr-qc] .
  * Forteza  _et al._ [2023] X. J. Forteza, S. Bhagwat, S. Kumar, and P. Pani, Phys. Rev. Lett. 130, 021001 (2023), arXiv:2205.14910 [gr-qc] .
  * Berti and Klein [2014] E. Berti and A. Klein, Phys. Rev. D 90, 064012 (2014), arXiv:1408.1860 [gr-qc] .
  * Isi  _et al._ [2019] M. Isi, M. Giesler, W. M. Farr, M. A. Scheel, and S. A. Teukolsky, Phys. Rev. Lett. 123, 111102 (2019), arXiv:1905.00869 [gr-qc] .
  * Cotesta  _et al._ [2022] R. Cotesta, G. Carullo, E. Berti, and V. Cardoso, Phys. Rev. Lett. 129, 111102 (2022), arXiv:2201.00822 [gr-qc] .
  * Nissanke  _et al._ [2010] S. Nissanke, D. E. Holz, S. A. Hughes, N. Dalal, and J. L. Sievers, Astrophys. J. 725, 496 (2010), arXiv:0904.1017 [astro-ph.CO] .
  * Rodriguez  _et al._ [2014] C. L. Rodriguez, B. Farr, V. Raymond, W. M. Farr, T. B. Littenberg, D. Fazi, and V. Kalogera, Astrophys. J. 784, 119 (2014), arXiv:1309.3273 [astro-ph.HE] .
  * gwt [2022] https://dcc.ligo.org/ligo-p1800370/public (2022).
  * Carullo  _et al._ [2023] G. Carullo, W. Del Pozzo, and J. Veitch, pyRing: a time-domain ringdown analysis python package, git.ligo.org/lscsoft/pyring (2023).
  * [73] https://johnveitch.github.io/cpnest/.
  * Cook  _et al._ [2006] S. R. Cook, A. Gelman, and D. B. Rubin, Journal of Computational and Graphical Statistics 15, 675 (2006).
  * Talts  _et al._ [2018] S. Talts, M. Betancourt, D. Simpson, A. Vehtari, and A. Gelman, arXiv e-prints , arXiv:1804.06788 (2018), arXiv:1804.06788 [stat.ME] .
  * Karchev  _et al._ [2022] K. Karchev, R. Trotta, and C. Weniger 10.1093/mnras/stac3785 (2022), arXiv:2209.06733 [astro-ph.CO] .
  * Romero-Shaw  _et al._ [2020] I. M. Romero-Shaw _et al._ , Mon. Not. Roy. Astron. Soc. 499, 3295 (2020), arXiv:2006.00714 [astro-ph.IM] .
  * Bhagwat  _et al._ [2018] S. Bhagwat, M. Okounkova, S. W. Ballmer, D. A. Brown, M. Giesler, M. A. Scheel, and S. A. Teukolsky, Phys. Rev. D 97, 104065 (2018), arXiv:1711.00926 [gr-qc] .
  * Kamaretsos  _et al._ [2012b] I. Kamaretsos, M. Hannam, S. Husa, and B. S. Sathyaprakash, Phys. Rev. D 85, 024018 (2012b), arXiv:1107.0854 [gr-qc] .
  * Barausse  _et al._ [2012] E. Barausse, V. Morozova, and L. Rezzolla, Astrophys. J. 758, 63 (2012), [Erratum: Astrophys.J. 786, 76 (2014)], arXiv:1206.3803 [gr-qc] .
  * Hofmann  _et al._ [2016] F. Hofmann, E. Barausse, and L. Rezzolla, Astrophys. J. Lett. 825, L19 (2016), arXiv:1605.01938 [gr-qc] .
  * gwt [2020] https://gwosc.org/gwtc-1/ (2020).
  * Abbott  _et al._ [2020a] R. Abbott _et al._ (LIGO Scientific, Virgo), Astrophys. J. Lett. 900, L13 (2020a), arXiv:2009.01190 [astro-ph.HE] .
  * Isi and Farr [2022] M. Isi and W. M. Farr, (2022), arXiv:2202.02941 [gr-qc] .
  * Abbott  _et al._ [2020b] R. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 125, 101102 (2020b), arXiv:2009.01075 [gr-qc] .
  * Finch and Moore [2022] E. Finch and C. J. Moore, Phys. Rev. D 106, 043005 (2022), arXiv:2205.07809 [gr-qc] .
  * Branchesi  _et al._ [2023] M. Branchesi _et al._ , JCAP 07, 068, arXiv:2303.15923 [gr-qc] .
  * Maselli  _et al._ [2024] A. Maselli, S. Yi, L. Pierini, V. Vellucci, L. Reali, L. Gualtieri, and E. Berti, Phys. Rev. D 109, 064060 (2024), arXiv:2311.14803 [gr-qc] .
  * Yi  _et al._ [2024] S. Yi, A. Kuntz, E. Barausse, E. Berti, M. H.-Y. Cheung, K. Kritos, and A. Maselli, (2024), arXiv:2403.09767 [gr-qc] .
  * Vallisneri [2008] M. Vallisneri, Phys. Rev. D 77, 042001 (2008), arXiv:gr-qc/0703086 .
  * Virtanen  _et al._ [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, Nature Methods 17, 261 (2020).
  * Harris  _et al._ [2020] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, Nature 585, 357 (2020).
  * Hunter [2007] J. D. Hunter, Computing in Science & Engineering 9, 90 (2007).
  * Waskom [2021] M. L. Waskom, Journal of Open Source Software 6, 3021 (2021).
  * LeCun  _et al._ [1998] Y. LeCun, L. Bottou, G. B. Orr, and K. R. Müller, Efficient backprop, in  _Neural Networks: Tricks of the Trade_, edited by G. B. Orr and K.-R. Müller (Springer Berlin Heidelberg, Berlin, Heidelberg, 1998) pp. 9–50.
  * George and Huerta [2018] D. George and E. A. Huerta, Phys. Rev. D 97, 044039 (2018), arXiv:1701.00008 [astro-ph.IM] .