# Black Hole Spectroscopy for Precessing Binary Black Hole Coalescences

Hengrui Zhu  0000-0001-9027-4184<EMAIL_ADDRESS>Department of
Physics, Princeton University, Jadwin Hall, Washington Road, New Jersey,
08544, USA  Princeton Gravity Initiative, Princeton University, Princeton, New
Jersey, 08544, USA  Harrison Siegel  0000-0001-5161-4617<EMAIL_ADDRESS>Department of Physics, Columbia University, 704 Pupin Hall, 538 West 120th
Street, New York, New York 10027, USA  Center for Computational Astrophysics,
Flatiron Institute, New York NY 10010, USA  Keefe Mitman  0000-0003-0276-3856
<EMAIL_ADDRESS>Theoretical Astrophysics, Walter Burke Institute for
Theoretical Physics, California Institute of Technology, Pasadena, California
91125, USA  Maximiliano Isi  0000-0001-8830-8672 Center for Computational
Astrophysics, Flatiron Institute, New York NY 10010, USA  Will M. Farr
0000-0003-1540-8562 Center for Computational Astrophysics, Flatiron Institute,
New York NY 10010, USA  Department of Physics and Astronomy, Stony Brook
University, Stony Brook NY 11794, USA
Michael Boyle  0000-0002-5075-5116 Cornell Center for Astrophysics and
Planetary Science, Cornell University, Ithaca, New York 14853, USA  Nils Deppe
0000-0003-4557-4115 Cornell Center for Astrophysics and Planetary Science,
Cornell University, Ithaca, New York 14853, USA  Department of Physics,
Cornell University, Ithaca, NY, 14853, USA  Lawrence E. Kidder
0000-0001-5392-7342 Cornell Center for Astrophysics and Planetary Science,
Cornell University, Ithaca, New York 14853, USA  Sizheng Ma
0000-0002-4645-453X Perimeter Institute for Theoretical Physics, 31 Caroline
Streeth North, Waterloo, Onatrio NSL 2Y5, Canada  Jordan Moxon
0000-0001-9009-6955 Theoretical Astrophysics, Walter Burke Institute for
Theoretical Physics, California Institute of Technology, Pasadena, California
91125, USA  Kyle C. Nelli  0000-0003-2426-8768 Theoretical Astrophysics,
Walter Burke Institute for Theoretical Physics, California Institute of
Technology, Pasadena, California 91125, USA  Harald P. Pfeiffer
0000-0001-9288-519X Max Planck Institute for Gravitational Physics (Albert
Einstein Institute), Am Mühlenberg 1, D-14476 Potsdam, Germany  Mark A. Scheel
0000-0001-6656-9134 Theoretical Astrophysics, Walter Burke Institute for
Theoretical Physics, California Institute of Technology, Pasadena, California
91125, USA  Saul A. Teukolsky  0000-0001-9765-4526 Theoretical Astrophysics,
Walter Burke Institute for Theoretical Physics, California Institute of
Technology, Pasadena, California 91125, USA  Cornell Center for Astrophysics
and Planetary Science, Cornell University, Ithaca, New York 14853, USA
William Throwe  0000-0001-5059-4378 Cornell Center for Astrophysics and
Planetary Science, Cornell University, Ithaca, New York 14853, USA  Vijay
Varma  0000-0002-9994-1761 Department of Mathematics, Center for Scientific
Computing and Data Science Research, University of Massachusetts, Darmouth, MA
02747, USA  Nils L. Vu  0000-0002-5767-3949 Theoretical Astrophysics, Walter
Burke Institute for Theoretical Physics, California Institute of Technology,
Pasadena, California 91125, USA

###### Abstract

To accurately perform black hole spectroscopy, it is essential to know which
quasinormal modes dominate astrophysical ringdown signals. In this Letter, we
present a phenomenological description of the quasinormal modes that are
excited in the ringdowns of comparable mass, quasi-circular precessing binary
black hole coalescences. By analyzing an exhaustive catalog of numerical
relativity simulations, we confirm that the relative fundamental quasinormal
mode amplitudes of precessing systems are related to those of non-precessing
systems by a simple rotation, and that additional structure in the spectrum is
connected to the system’s kick velocity and other asymmetries in the orbital
dynamics. We find that the ringdowns of precessing systems need not be
dominated by the ${(\ell,m)=(2,\pm 2)}$ quasinormal modes. These results build
upon previous works on waveform modeling, and are consistent with a recent
ringdown analysis of the LIGO-Virgo gravitational wave signal GW190521.

_Introduction.—_ According to general relativity, the coalescence of two black
holes results in the formation of a perturbed remnant, which equilibrates to a
Kerr black hole by emitting gravitational waves in a process called the
_ringdown_ [1, 2, 3]. The ringdown emission consists of a transient burst that
is succeeded by a spectrum of quasinormal modes (QNMs), which then themselves
decay to reveal a power-law tail [4, 5]. The QNMs are exponentially damped
sinusoids with complex frequencies that are uniquely determined by the mass
and spin of the remnant black hole [6, 7, 1, 8, 9] and complex amplitudes that
are nontrivially related to the initial conditions of the system [10].

Much remains to be understood about how the progenitor properties of
coalescing black holes map to the QNMs radiated by the remnant black hole.
Numerical relativity studies of the ringdowns of comparable mass quasi-
circular binary black hole coalescences have largely been confined to non-
precessing systems, where the black hole spins are perpendicular to the
orbital plane. These studies have shown empirically that some progenitor
information, such as the mass ratio and spin magnitudes, is encoded in the QNM
amplitudes [11, 12, 13, 14, 15, 16, 17].

For precessing systems [18], where the progenitor black holes possess spin
components parallel to the orbital plane, the picture becomes more
complicated. In the extreme mass-ratio case, the relative amplitude ratios
between QNMs are dictated by the spin magnitude of the heavier black hole and
a small number of orbital parameters [19, 20, 21, 22]. For comparable mass
systems, waveform models exist which map non-precessing waveforms to
precessing ones through a non-inertial frame rotation to the co-precessing
frame [23, 24, 25], and these mappings have been extended to ringdown waveform
modeling in, e.g., Refs. [24, 26, 27]. Studies of the mismatch between NR
waveforms and explicit QNM models, building on Ref. [28], have also been
performed for precessing systems [29].

In this work, we focus on fitting the ringdowns of precessing systems from the
perspective of observational black hole spectroscopy and perturbation theory,
synthesizing conceptual aspects of previous results and extending them to
cover a larger parameter space of mass ratios, binary spin configurations, and
remnant spin misalignments. We fit the amplitude of each individual QNM in the
canonical frame for perturbation theory, i.e., the superrest frame of the
remnant [30, 31, 32]. This means that any large amplitudes we recover will
correspond to large contributions of the associated modes’ frequencies in
LIGO-Virgo-KAGRA [33, 34, 35] data. We analyze the ringdown within 250
numerical relativity simulations of both precessing and non-precessing binary
black hole coalescences [36], and we highlight links between the amplitudes of
the QNMs in precessing systems and two mechanisms: (1) the misalignment of the
total angular momentum flux $\mathrm{d}\mathbf{J}/\mathrm{d}t$ near the time
of merger with the final spin $\bm{{\chi}}_{f}$ of the remnant black hole, and
(2) asymmetries in the inspiral radiation on either side of the orbital plane.

The misalignment angle between $\mathrm{d}\mathbf{J}/\mathrm{d}t$ and
$\bm{{\chi}}_{f}$ approximately relates the QNM amplitudes of many precessing
and non-precessing systems by a simple rotation, and additional structure in
the QNM amplitude ratios is correlated with hemispheric asymmetries, thus
holding information about the direction of the system’s kick velocity and
other orbital dynamics. We show that in some precessing systems, the dominant
QNM will not be an ${(\ell,m)=(2,\pm 2)}$ QNM. We also find that the
precession-induced asymmetry in the gravitational wave emission can result in
$\mathcal{O}(10)$ amplitude ratios between $\pm m$ modes which share the same
frequency. Our findings have broad implications for the interpretation of
ringdown signals in data from LIGO-Virgo-KAGRA and future detectors, and may
also inform investigations into the nonlinear aspects of the ringdown.

Our plotting and fitting code, as well as an example simulation output, can be
found in Ref. [37].

_QNM Conventions.—_ Throughout this Letter, we will refer to individual QNMs
by four indices $(p,\ell,m,n)$, following the notation in [38]. The angular
structure of the QNMs is described by spheroidal harmonics with angular
indices $\ell$ and $m$. The radial structure of the QNMs is denoted by the
index $n$, and is also tied to the lifetime of the QNMs: the longest-lived
${n=0}$ QNMs are referred to as fundamental QNMs, and faster-decaying ${n>0}$
QNMs are called overtones. For a given set of $(\ell,m,n)$, when ${m\neq 0}$
there are two distinct QNMs which are labeled by an index
$p\equiv\mathrm{sgn}(m\,\mathrm{Re}\left[\omega\right])$ and whose phase
fronts are either co-rotating (${p=+}$) or counter-rotating (${p=-}$) with the
black hole; solutions with ${m=0}$ are azimuthally symmetric, so that there is
no notion of co- vs counter-rotating fronts and the two possible signs of
$\mathrm{Re}\left[\omega\right]$ directly encode the polarization state. We
define $A_{(p,\ell,m,n)}$ as the norm of the complex amplitude for a given
$(p,\ell,m,n)$ QNM.

_Methods.—_ We used simulations of 226 precessing and 24 non-precessing binary
black hole coalescences, with mass ratios ${q\leq 8}$ and component spin
magnitudes ${\chi\leq 0.8}$, created with the Spectral Einstein Code (SpEC)
[36]. We computed the strain $h$ and the five Weyl scalars at future null
infinity using the SpECTRE code’s Cauchy-characteristic evolution (CCE) module
[39, 40, 41]. With this asymptotic data, we mapped the simulated systems to
the superrest frame of their remnant black hole $250M$ past the peak of the
strain’s luminosity (where $M$ is the total Christodoulou mass of the binary),
using the BMS frame fixing procedure outlined in Ref. [32] and the code scri
[42, 43, 44, 45]. Consequently, each simulation has been transformed such that
its remnant black hole is at rest at the origin, has its spin in the
positive-$\hat{z}$ direction, and has no supermomentum [31, 32, 30]. This
post-processing is necessary in order to robustly extract physical QNM
amplitudes.

For each simulation, we simultaneously fit the strain using a QNM model with
all fundamental modes $(p,\ell,m,0)$ for $p\in\\{+,-\\}$, $\ell\in\\{2,3\\}$,
and $m\in\\{-\ell,-\ell+1,...,\ell\\}$. This model is fit over the two-sphere
using a linear least-squares routine [46, 14, 47] and the following procedure:

  * •

Fit the waveform from $t_{0}$ to $100M$ past $t_{\mathrm{peak}}$, where
$t_{0}$ is the start time of the fit and is evenly varied from $20M$ to $80M$
past $t_{\mathrm{peak}}$—note that $t_{\mathrm{peak}}$ is the time at which
the strain’s luminosity reaches its maximum value (see Eq. (2.8) of Ref.
[48]);

  * •

Over a series of windows in $t_{0}$ with a length of $20M$, compute the
fractional uncertainty for each of the QNM amplitudes;

  * •

Compute the mean of these fractional uncertainties over the fitted QNMs;

  * •

Check to see which $t_{0}$ window has the minimum mean fractional uncertainty
and use this window to extract the amplitude (as a mean over the window) of
each of the QNMs in the ringdown model.

This routine ensures that the QNMs are extracted over the window of start
times in which their recovered amplitudes are all, on average, the most
stable. We do this rather than simply fitting at a fixed time because the
timescales of the QNMs vary from system to system, depending on the magnitude
of the remnant’s spin. Over our entire catalog, this routine results in fits
that have window start times ${t-t_{\mathrm{peak}}\approx 33M\pm 7M}$ and
relative errors over the two-sphere of $\sim 0.8\%$.

Figure 1: Left: Amplitude ratios of the co-rotating $({\ell,m)=(2,\pm 1)}$ and
$({\ell,m)=(2,\pm 2)}$ fundamental QNMs for every simulation in our catalog,
where $A_{(p,\ell,\pm m,0)}^{2}=A_{(p,\ell,+m,0)}^{2}+A_{(p,\ell,-m,0)}^{2}$.
The misalignment angle $\theta$ between the final remnant spin and the angular
momentum flux at the time of peak luminosity (horizontal axis, Eq. (1)) is
correlated with the excitation of the co-rotating $({\ell,m)=(2,\pm 1)}$
fundamental QNMs (color), for a large range of mass ratios $q=m_{1}/m_{2}\geq
1$ (vertical axis). Higher mass-ratio systems can have larger values of
$\theta$, since the spin of the more massive black hole can contribute a
larger fraction of the remnant’s final spin in these systems. The monotonic
increase of $A_{(+,2,\pm 1,0)}/A_{(+,2,\pm 2,0)}$ does not indicate that
$A_{(+,2,\pm 1,0)}$ are the most dominant QNMs for large $\theta$: both
$A_{(+,2,\pm 1,0)}$ and $A_{(+,2,\pm 2,0)}$ decrease for $\theta>\pi/2$ and
are dominated by counter-rotating modes, e.g., $A_{(-,2,\pm 2,0)}$. Right:
Amplitude ratios of summed co-rotating and counter-rotating $({\ell,m)=(2,\pm
1)}$ and $({\ell,m)=(2,\pm 2)}$ fundamental QNMs, with $A_{(\pm,\ell,\pm
m,0)}$ defined in Eq. (5). One can roughly predict QNM amplitude ratios for
comparable mass precessing binaries by rotating an $({\ell,m)=(2,\pm 2)}$
spherical harmonic by $\theta$ off the $z$ axis and taking the resulting ratio
of different $(\ell,m)$ harmonics in the rotated frame. This prediction for
$q=1$ systems is shown as a black curve (Eq. (4) in the Supplemental
Material). For higher mass-ratio systems, the dominant ${\theta=0}$
perturbation to be rotated will include significant $\ell\not=m$ content,
dependent on the mass ratio and progenitor spins [11, 12]. An example for the
rotation of a $q=8$, $\chi_{\mathrm{eff}}=-0.6$ perturbation is shown as a
blue curve; see Sec. 6 in the Supplemental Material.

In addition to the QNM amplitudes, we also compute two other quantities in the
superrest frame of the remnant: the angle $\theta$ between the total angular
momentum flux $\mathrm{d}\mathbf{J}/\mathrm{d}t$ at $t_{\mathrm{peak}}$ and
the final remnant spin $\bm{{\chi}}_{f}$, and the angle $\phi$ of the remnant
kick velocity $\bm{v}$ with respect to $\bm{{\chi}}_{f}$. We formally define
the remnant spin misalignment angle $\theta$ as

$\displaystyle\theta=\mathrm{cos}^{-1}\left[\frac{\mathrm{d}\mathbf{J}/\mathrm{d}t\left(t=t_{\mathrm{peak}}\right)}{||\mathrm{d}\mathbf{J}/\mathrm{d}t\left(t=t_{\mathrm{peak}}\right)||}\cdot\frac{\bm{{\chi}}_{f}}{||\bm{{\chi}}_{f}||}\right],$
(1)

where $\mathrm{d}\mathbf{J}/\mathrm{d}t$ is the total angular momentum flux
(see Eq. (2.24) of Ref. [48]), $\bm{{\chi}}_{f}$ is the spin of the remnant
black hole measured at future null infinity (see Eq. (15) of Ref. [49]), and
$||\cdot||$ is the $L^{2}$ norm. This angle shares similarities with the
precession angles and optimal emission direction as identified in, e.g., Refs.
[27, 23, 24, 25], but will be used here only to characterize theoretical
predictions of the remnant perturbation structure. The kick angle $\phi$ is
given by

$\displaystyle\phi$
$\displaystyle=\cos^{-1}\left[\frac{\bm{v}_{\mathrm{CoM}}^{\mathrm{remnant}}-\bm{v}_{\mathrm{CoM}}^{\mathrm{binary}}}{||\bm{v}_{\mathrm{CoM}}^{\mathrm{remnant}}-\bm{v}_{\mathrm{CoM}}^{\mathrm{binary}}||}\cdot\frac{\bm{{\chi}}_{f}}{||\bm{{\chi}}_{f}||}\right]$
$\displaystyle=\cos^{-1}\left[-\frac{\bm{v}_{\mathrm{CoM}}^{\mathrm{binary}}}{||\bm{v}_{\mathrm{CoM}}^{\mathrm{binary}}||}\cdot\frac{\bm{{\chi}}_{f}}{||\bm{{\chi}}_{f}||}\right],$
(2)

where $\bm{v}_{\mathrm{CoM}}^{\mathrm{remnant}}=0$ in the superrest frame of
the remnant, and $\bm{v}_{\mathrm{CoM}}^{\mathrm{binary}}$ is the pre-merger
center-of-mass velocity (see Eq. (12) of Ref. [49]) of the binary in the same
frame. We extract the latter over a $500M$ window starting at
$t=t_{\mathrm{peak}}-1000M$ using the procedure outlined in Ref. [32].

We make use of the Wigner-D matrices
$\mathfrak{D}_{m^{\prime},m}^{\ell}(\mathbf{R})$ that govern the mixing of
spin-weighted spherical harmonics under rotations through the transformation
[43, 50, 51, 52]

$\phantom{}_{s}Y_{(\ell,m^{\prime})}^{\prime}=\sum_{m}\mathfrak{D}_{m^{\prime},m}^{\ell}\left(\mathbf{R}\right)\phantom{}_{s}Y_{(\ell,m)}\,,$
(3)

where $\mathbf{R}$ is some quaternion, $\phantom{}{}_{s}Y_{(\ell,m)}$ is a
spin-weight $s$ spherical harmonic, and primes indicate a rotated frame. The
Wigner-D matrices mix amplitudes of harmonics with the same $\ell$ but
different $m$. Because we will be strictly interested in static rotations off
the $z$ axis by the remnant spin misalignment angle $\theta$, we simplify our
notation to

$\displaystyle\mathfrak{D}_{m^{\prime},m}^{\ell}(\mathbf{R})\rightarrow\mathfrak{D}_{m^{\prime},m}^{\ell}(\theta),$
(4)

where now $\mathfrak{D}_{m^{\prime},m}^{\ell}(\theta)$ represents a rotation
by an angle $\theta$ off the $z$ axis, but not about it. If we consider co-
rotating and counter-rotating perturbations, it can be shown that a fully co-
rotating $+m$ perturbation, when rotated by ${\theta=\pi}$ with respect to the
remnant spin, would instead excite fully counter-rotating $-m$ modes [53].

Figure 2: Top: Amplitude ratio of the co-rotating $({\ell,m)=(2,1)}$ and
$(2,2)$ fundamental QNMs, versus remnant spin misalignment angle $\theta$; see
Eq. (1). Color shows the the kick velocity angle $\phi$ with respect to the
remnant spin direction; see Eq. (Black Hole Spectroscopy for Precessing Binary
Black Hole Coalescences). The black curve is the prediction for this amplitude
ratio based on the rotation of spin-weighted spherical harmonics; see Eq. (3).
We demonstrate that the kick direction corresponds to the spread of QNM
amplitude ratios; the kick direction is strongly tied to the asymmetry of the
most dominant modes in the inspiral, and these asymmetries are in turn passed
on to the QNM amplitudes. Bottom: Now showing the amplitude ratio of the
counter-rotating $(-,2,-2,0)$ and co-rotating $(+,2,2,0)$ fundamental QNMs.
Note for $\theta\gtrsim\pi/2$ the counter-rotating QNMs dominate [27].
Fractional uncertainties here are larger, since the absolute values of the
counter-rotating amplitudes are close to zero for small values of $\theta$.

_Results and Discussion.—_ As shown in the left panel of Fig. 1, there is a
clear correlation between the remnant spin misalignment angle $\theta$ and the
excitation of the $(+,2,\pm 1,0)$ QNMs over our entire catalog of both
precessing and non-precessing simulations. This correlation has been
previously theorized [24], but to our knowledge has not been rigorously
explored. Below we detail a conceptual framework which explains this
phenomenon by characterizing the structure of the remnant perturbations. This
framework elucidates the underlying perturbation theory that motivates
existing waveform modeling approaches to the ringdowns of precessing systems
[24, 26, 27, 22].

Assume that the angular content of the remnant perturbation is related to that
of the binary and is oriented along the orbital angular momentum $\mathbf{L}$
at a time near merger, e.g., the time of peak luminosity [11, 12, 15].
Although the notion of an orbit breaks down near merger, it is still useful
conceptually for tracking the angular content of the perturbations. Because we
do not have direct access to $\mathbf{L}$ in our simulations, in practice we
use the total angular momentum flux $\mathrm{d}\mathbf{J}/\mathrm{d}t$ which
is parallel to $\mathbf{L}$ up to 1 post-Newtonian (PN) order [54, 55]. In
linear perturbation theory, the QNMs are sourced by perturbations to a fixed
Kerr background with the $z$ axis aligned with the black hole spin. Provided
that this is a valid description of the ringdown at some time, we must study
the structure of the perturbation in a frame aligned with the final remnant
spin $\bm{{\chi}}_{f}$. In precessing systems, $\bm{{\chi}}_{f}$ will not
necessarily be parallel to $\mathbf{L}$. This means that, in the remnant
frame, the perturbations need not have $z$ as a symmetry axis. Lastly, because
the dominant angular content of the inspiral emission in the co-orbital frame
for precessing systems is similar to that of non-precessing systems (up to
asymmetries over the orbital plane) [52, 56, 25, 23, 44], it is natural to
expect that the geometry of the remnant perturbation is roughly the same in
both cases. Therefore, the QNMs of a precessing system should be approximately
sourced by a rotated version of the perturbation that acts on the remnant of a
non-precessing system, and the QNM amplitudes of precessing and non-precessing
systems should thus be related by a simple rotation.

Given the above conceptual picture, we can obtain an approximation for the
amplitude ratios of QNMs in precessing systems by rotating the angular content
of the dominant perturbation expected in a non-precessing $\theta=0$ system.
For simplicity we assume, to good approximation for equal mass binaries, that
the dominant ${\theta=0}$ perturbation is comprised solely of equal-amplitude
co-rotating ${(\ell,m)=(2,\pm 2)}$ harmonics. With this assumption, we predict
that counter-rotating QNM amplitudes are large for ${\theta>\pi/2}$, and
${\ell\neq m}$ QNM amplitudes are large for ${\theta\simeq\pi/2}$. We can also
account for larger mass ratios by additionally rotating ${l\neq m}$
perturbations [11, 12].

In the right panel of Fig. 1, we show that this rotational prediction tracks
the QNM amplitude ratio $A_{(\pm,2,\pm 1,0)}/A_{(\pm,2,\pm 2,0)}$ well as a
function of $\theta$, where we have defined

$\begin{split}A_{(\pm,\ell,\pm
m,0)}^{2}&=A_{(+,\ell,m,0)}^{2}+A_{(-,\ell,m,0)}^{2}\\\
&\phantom{=.}+A_{(+,\ell,-m,0)}^{2}+A_{(-,\ell,-m,0)}^{2},\end{split}$ (5)

i.e., the norm of a multi-dimensional amplitude vector. Notably we observe
that for ${\theta\simeq\pi/2}$, the most dominant QNMs in comparable mass
precessing systems will not be the co-rotating or counter-rotating
$({\ell,m)=(2,\pm 2)}$ QNMs, as is often assumed in spectroscopic tests, but
rather the $(2,\pm 1)$ or the $(2,0)$ QNMs (see Fig. 2 of the Supplemental
Material).

Additionally, precessing systems exhibit asymmetries in their inspiral
emission above and below the orbital plane; the asymmetries of the binary can
be passed on to the perturbations of the remnant to create differences in the
amplitudes of the $\pm m$ QNM polarizations [22, 17]. This is partly
responsible for the spread around our naive rotational prediction in Fig. 1.
Asymmetric emission is associated with nonzero remnant kick velocities [57,
58, 59, 45, 60, 58]. Here we show how the kick direction leaves imprints in
the ringdown which neatly correspond to deviations from our rotational
predictions.

In the top panel of Fig. 2 we compare the kick direction with spreads around
our rotational prediction for $A_{(+,2,1,0)}/A_{(+,2,2,0)}$; looking at only
one sign of $m$ should make us more sensitive to hemispheric asymmetries. For
$\theta\lesssim\pi/4$, where the $({\ell,m)=(2,\pm 2)}$ modes dominate, we
find that the spread is indeed correlated with the kick velocity’s direction:
a kick in the positive-$\hat{z}$ direction corresponds to less radiated power
in the $(2,2)$ mode, while a kick in the negative-$\hat{z}$ direction
corresponds to more radiated power in the $(2,2)$ mode. However, when
$\pi/4\lesssim\theta\lesssim 3\pi/4$, then ${A_{(+,2,1,0)}\simeq
A_{(+,2,2,0)}}$, and the kick velocity is no longer as strongly correlated
with only one $(\ell,\pm m)$ QNM.

In the bottom panel of Fig. 2, we show a similar result when comparing the
counter-rotating $(-,2,-2,0)$ and co-rotating $(+,2,2,0)$ QNM amplitudes;
these amplitudes should be related as per the discussion in our _Methods_
section. When ${\theta>\pi/2}$, the counter-rotating QNMs dominate over the
co-rotating QNMs [27]. Again, the spread of the amplitude ratio around the
rotational prediction is in part governed by the direction of the kick
velocity with respect to the remnant spin.

The spread of systems near ${\theta=0}$ in the right panel of Fig. 1, even for
nearly equal mass systems, can also be partly attributed to asymmetric
emission. Many of these systems are precessing despite having zero remnant
spin misalignment, as discussed in our Supplemental Material. Most have
progenitors whose spins are almost equal and opposite, meaning that the
remnant spin direction is dominated by the direction of $\mathbf{L}$. These
configurations can result in large kicks and correspondingly large asymmetries
in QNM emission.

Lastly, the phases and magnitudes of the precession-induced asymmetries in the
inspiral are not all the same across $\ell$ and $m$ indices in the remnant
frame, and consequently the QNM amplitude ratio
$A_{(p,\ell,m,n)}/A_{(p,\ell,-m,n)}$ need not always have the same value for
all $\ell$ and $m$ in a given system. PN expansions suggest that the magnitude
of asymmetry in the dominant modes during the inspiral can become
$\mathcal{O}(1)$ close to merger [44]; we show in Fig. 3 that the degree of
asymmetry in certain QNMs can be $\mathcal{O}(10)$ in the most extreme
systems, and can be distributed in a variety of ways, presumably based on the
orbital phase and the particular spin configuration of each system. These QNM
asymmetries can be especially prominent for large mass ratio systems, since
the smaller progenitor can predominantly merge into the larger progenitor’s
northern or southern hemispheres, thus cleanly exciting mostly one of the
${\pm m}$ polarizations in a way that may not be possible in nearly equal mass
systems [20]. As a result of emission asymmetries, precessing systems can
produce radically different observable QNM amplitude spectra depending on the
viewing angle of any terrestrial detectors.

Figure 3: Amplitude ratio of co-rotating $\pm m$ QNMs, using a small subset of
simulations in our catalog. In precessing systems there are asymmetries in the
emission power above and below the orbital plane, which in turn lead to
asymmetries in the QNM amplitudes over the remnant equatorial plane. The
asymmetries for given $\ell$ and $|m|$ harmonics can be large,
$\mathcal{O}(10)$ in the extreme cases shown here, and distributed in diverse
ways for different systems. Errors on these ratios are mostly too small to be
seen.

_Conclusion.—_ By fitting the first-order quasinormal mode amplitudes in the
ringdowns of 250 numerical relativity simulations of both non-precessing and
precessing binary black hole coalescences, we comprehensively demonstrate in
this Letter that the remnant perturbations in precessing systems may be
thought of as being approximately related to those of non-precessing systems
by a simple rotation. We characterize this rotation by the misalignment of two
quantities: the total angular momentum flux $\mathrm{d}\mathbf{J}/\mathrm{d}t$
near the time of merger and the final spin $\bm{{\chi}}_{f}$ of the remnant
black hole.

Notably, we show that in some precessing systems the dominant quasinormal mode
will not be an ${(\ell,m)=(2,\pm 2)}$ quasinormal mode, but rather an
${(\ell,m)=(2,\pm 1)}$ or an ${(\ell,m)=(2,0)}$ quasinormal mode, contrary to
what is commonly assumed in ringdown analyses of LIGO-Virgo-KAGRA data [61,
62, 63, 64, 65, 66]. Additionally, precession-induced emission asymmetries,
which are correlated with remnant kicks and other orbital dynamics, are
encoded in the quasinormal mode spectra. These asymmetries can be large, such
that amplitude ratios of individual $\pm m$ quasinormal modes which share the
same frequency can be $\mathcal{O}(10)$ in some precessing systems. This
implies that for real measurements, the viewing angle can dramatically alter
the observed quasinormal mode spectrum of a precessing system.

We emphasize that the excitation of ${\ell\neq|m|}$ modes does not arise from
operating in a non-inertial frame, as our fitting is performed in the
superrest frame of the remnant, i.e., the inertial frame in which quasinormal
modes are defined by perturbation theory. Our fits return the amplitude
attached to the specific frequency of a physical $(p,\ell,m,n)$ quasinormal
mode. This means that any large amplitudes we recover will correspond to large
contributions of those associated frequencies in observed LIGO-Virgo-KAGRA
data, which must be taken into account in order to accurately perform black
hole spectroscopy.

Our results are consistent with a recent analysis by Siegel et al. [65] of the
LIGO-Virgo gravitational wave signal GW190521_030229, wherein ringdown models
with a large ${(\ell,|m|)=(2,1)}$ fundamental quasinormal mode amplitude were
found to produce remnant mass and spin estimates in agreement with those of
most inspiral-merger-ringdown waveforms. Furthermore, our results are similar
to those of Hughes et al. [21], where it was found that the excitation of
quasinormal modes in extreme mass-ratio systems is correlated with the polar
angle of plunge of the smaller mass. It is possible that the framework we have
explored herein may be valid over all possible binary mass ratios. Our work is
also broadly consistent with numerous pivotal advances in gravitational
waveform modeling [52, 56, 44, 27, 23, 24, 25, 22, 26, 67], and seeks to
illuminate the fundamental physics underlying these models.

In future work, we will explore nonlinear effects in the ringdowns of
precessing systems. One might think that at very early times in the ringdown
of a precessing system, it could be necessary to consider the remnant spin
direction to be dynamically changing with time. Our Letter suggests that at
least to linear order, information from the point of peak luminosity imprints
itself on the ringdown of a precessing system in a straightforward way.

##  I Acknowledgements

We thank Katerina Chatziioannou, Scott Hughes, Yuri Levin, Frans Pretorius,
and Aaron Zimmerman for valuable conversations. We are also grateful to
Eleanor Hamilton, Mark Hannam, Sascha Husa, Lionel London, and Geraint Pratten
for their insight on this work in the context of gravitational waveform
modeling. Computations for this work were performed with the Wheeler cluster
at Caltech. This work was supported in part by the Sherman Fairchild
Foundation and NSF Grants No. PHY-2011968, PHY-2011961, PHY-2309211,
PHY-2309231, OAC-2209656 at Caltech., as well as NSF Grants No. PHY-2207342
and OAC-2209655 at Cornell. V.V. acknowledges support from NSF Grant No.
PHY-2309301. The Flatiron Institute is a division of the Simons Foundation.
H.S.’s research is supported by Yuri Levin’s Simons Investigator Award 827103.
Software: seaborn [68], matplotlib [69], jupyter [70], numpy [71], Python3
[72].

## Supplemental material

Figure 4: Relative error of QNM fit vs. the remnant spin, colored by the
misalignment angle $\theta$. We find that the relative error decreases with
increasing spin, as the QNMs are longer-lived. Precessing binaries in our
catalog in general yield remnant spins smaller than the non-precessing
systems, thereby resulting in higher error, but the relative error remains on
the order of 1%. Figure 5: Parity breaking of $\pm m$ QNMs over the remnant
equator as computed via Eq. (8) (vertical axis), versus the misalignment angle
$\theta$ (horizontal axis). Non-precessing systems maintain reflection
symmetry over the remnant equator in the ringdown, whereas precessing systems
do not: here the events with small asymmetries (below the vertical axis break)
are non-precessing. Note that many precessing systems have ${\theta\approx 0}$
or $\pi$, i.e. remnant spin misalignment is not necessitated by precession.
Especially at ${\theta\approx 0}$, many precessing systems are characterized
by large kick rapidities: many of these systems have roughly equal component
spins fully anti-aligned with each other [59], and thus the total angular
momentum of the binary and subsequently the spin of the remnant is dominated
by $\mathbf{L}$.

_1\. Error of QNM fits.—_ To characterize the quality of the QNM fits, we
define the following fitting error:

$E(t_{0})\equiv\sum_{\ell\in\\{2,3\\}}\sum_{|m|\leq\ell}\frac{\int_{t_{0}}^{100M}\left|h_{(\ell,m)}-h_{(\ell,m)}^{\mathrm{QNMmodel}}\right|^{2}dt}{\int_{t_{0}}^{100M}\left|h_{(\ell,m)}\right|^{2}dt}$
(6)

In Fig. 4, we plot the mean of $E(t_{0})$ over the most stable window of each
simulation versus the magnitude of the remnant’s spin, colored by the
misalignment angle $\theta$. As can be seen, the fitting error is on the order
of 1% for every simulation. Furthermore, we find that fitting error decreases
with the magnitude of the remnant’s spin. This is likely because the QNMs
decay faster as the spin magnitude decreases, which implies that more of the
fit to the numerical waveform is influenced by numerical noise. As for the
connection to the misalignment angle $\theta$, precessing black hole binaries
require less fine-tuning of their spin configurations to produce slowly-
spinning remnants, since the spin of the progenitor black holes need not be
aligned with the angular orbital momentum. In the future, it could be useful
to perform fits which have a varying end time based on the magnitude of the
remnant’s spin, but we save this for future work.

_2\. Parity Breaking Over the Remnant Equator.—_ One way to quantify the
asymmetry of the $\pm m$ QNMs is by looking at their behavior under the
$z$-parity operation, which is defined to act on the strain such that

$P_{z}(h_{(\ell,m)})=(-1)^{\ell}~{}\overline{h}_{(\ell,-m)},$ (7)

where the overbar denotes a complex conjugate. We can quantify the overall
parity violation in the superrest frame by computing

$\displaystyle
P=\sqrt{\frac{\sum_{\ell,m}|h_{(\ell,m)}-(-1)^{\ell}~{}\overline{h}_{(\ell,-m)}|^{2}}{4\sum_{\ell,m}|h_{(\ell,m)}|^{2}}},$
(8)

as described in [44]. For non-precessing simulations, $P$ will be essentially
zero. However, precessing systems can exhibit significant parity breaking in
the $\pm m$ QNMs.

In Fig. 5 we show the parity violation of all the simulations in our catalog.
The non-precessing systems are clearly separated from the precessing systems,
as denoted by the vertical axis break. Note that there are many precessing
systems which have a remnant spin misalignment angle ${\theta\approx 0}$ or
$\pi$. Remnant spin misalignment is not strictly required in precessing
systems: fine-tuning of the spins, such as is found in superkick
configurations [59] for example, can allow for precessing systems with
${\theta\approx 0}$, and large mass ratio systems can have essentially
arbitrary $\theta$.

Figure 6: Amplitude ratio between co-rotating $(3,2)$, $(3,1)$, $(3,0)$, and
the $(3,3)$ fundamental QNMs. The rotational predictions (black curves) are
obtained by rotating a $(3,\pm 3)$ spherical harmonic by $\theta$, i.e. by
implicitly assuming that the dominant perturbation for ${\theta=0}$ consists
solely of co-rotating $(3,\pm 3)$ content. Because our catalog mostly contains
$q\neq 1$ precessing simulations, this assumption is reasonable. The kick
velocity is not as straightforwardly correlated with ${\ell=3}$ QNM amplitudes
as it is with ${\ell=2}$ QNM amplitudes.

_3\. Excitation of $\ell=3$ QNMs.—_ Like the $\ell=2$ QNMs, the excitation
amplitudes of $\ell=3$ QNMs should also be dictated by rotations. For non-
equal mass non-precessing systems, the $(3,3)$ mode typically dominates over
the $(3,2)$ mode [11]; thus we will assume that we can predict the $\ell=3$
QNM amplitudes of precessing systems by rotating a $(3,\pm 3)$ co-rotating
perturbation by $\theta$. In Fig. 6, we show that the relative amplitudes
between the co-rotating $\ell=3$, $m<3$ fundamental QNMs and the $(+,3,3,0)$
QNM follow this prediction. However, the spread around the prediction is no
longer as clearly explained by the kick velocity angle $\phi$.

Figure 7: Amplitude ratio of the co-rotating $(2,0)$ and $(2,2)$ fundamental
QNMs versus the system’s misalignment angle. Color shows the kick velocity
angle $\phi$ with respect to the remnant’s spin axis. The black curve is the
prediction for this amplitude ratio based on the rotation of spin-weighted
spherical harmonics via the Wigner-D matrices.

_4\. Excitation of $(\ell,m)=(2,0)$ QNMs.—_ In Fig. 7, we show that the
relative amplitude between the $(+,2,0,0)$ and $(+,2,2,0)$ QNMs also follows
the rotational prediction. Notably, we find that the $(+,2,0,0)$ QNM can
dominate over both the $(+,2,1,0)$ and $(+,2,2,0)$ QNMs when the misalignment
angle is close to $\pi/2$.

Figure 8: Amplitude ratio of the co-rotating $(2,-2)$ and $(2,2)$ fundamental
QNMs versus the system’s misalignment angle. Color shows the kick velocity
angle $\phi$ with respect to the remnant’s spin axis. In non-precessing
systems, the two modes have the same amplitude due to presence of a reflection
symmetry about the orbital plane. For precessing systems, however, this
symmetry is broken: when the remnant is kicked in the direction of its spin,
the $(2,-2)$ mode is larger than the $(2,2)$ mode and vice versa for a kick in
the opposite direction.

_5\. Asymmetries and Kicks.—_ In Fig. 8, we provide more evidence of how the
system’s kick velocity correlates with the spread of QNM amplitude ratios
around the rotational prediction. The ratio of the co-rotating $(2,-2)$ and
$(2,2)$ QNMs as a function of the misalignment angle $\theta$ is cleanly
spread around 1 by the direction of the system’s kick velocity. In particular,
if the kick is in the same (opposite) direction as the remnant’s spin axis,
then the ${m=-2}$ (${m=2}$) QNM is larger in amplitude. This is expected
because a nonzero kick in the direction of the remnant’s spin axis implies
that there is asymmetric emission of gravitational waves across the orbital
plane, which means that the perturbation that sources the ringdown must also
be asymmetric in the mirror modes. When one ${l=2}$ QNM is dominant, which is
generally the case for ${\theta\lesssim\pi/2}$, then the asymmetries in the
$\pm m$ polarizations of that mode dictate the kick direction. However, when
multiple ${l=2}$ QNMs have comparable magnitudes as is often the case for
${\theta\approx\pi/2}$, the kick is not as strongly correlated with any one
QNM.

_6. ${\mathbf{q\not=1}}$ rotational predictions.—_ As outlined in the main
text, the excitation of QNMs in many precessing systems is correlated with a
rotated version of the QNMs in a non-precessing system. For ${q=1}$ systems
with ${\chi_{\rm eff}=0}$, this relationship is simple in that one only needs
to rotate a ${(\ell,m)=(2,2)}$ perturbation by $\theta$. For systems with
${q\not=1}$, however, our conceptual picture implies that one must rotate a
${(\ell,m)=(2,1)}$ perturbation in addition to an ${(\ell,m)=(2,2)}$
perturbation. This means that the rotational prediction will be

$\frac{\mathfrak{D}_{1,2}^{\ell,\pm}(\theta,q,\chi_{\mathrm{eff}})}{\mathfrak{D}_{2,2}^{\ell,\pm}(\theta,q,\chi_{\mathrm{eff}})}~{},$
(9)

where

$\displaystyle\begin{split}\left[\mathfrak{D}_{m_{1},m_{2}}^{\ell,\pm}(\theta,q,\chi_{\mathrm{eff}})\right]^{2}&=\left[\frac{A_{(\ell,m_{1})}}{A_{(\ell,m_{2})}}(q,\chi_{\mathrm{eff}})\mathfrak{D}_{m_{1},m_{1}}^{\ell,\pm}(\theta)\right]^{2}\\\
&\phantom{=.}+\left[\mathfrak{D}_{m_{1},m_{2}}^{\ell,\pm}(\theta)\right]^{2}\end{split}$
(10)

and

$\left[\mathfrak{D}_{m_{1},m_{2}}^{\ell,\pm}(\theta)\right]^{2}=\left[\mathfrak{D}_{m_{1},m_{2}}^{\ell}(\theta)\right]^{2}+\left[\mathfrak{D}_{-m_{1},m_{2}}^{\ell}(\theta)\right]^{2}.$
(11)

For $A_{(2,1)}/A_{(2,2)}$, we use the result of Ref [11], i.e.,

$A_{(2,1)}/A_{(2,2)}\approx|0.43(\sqrt{1-4\nu}-\chi_{\mathrm{eff}})|,$ (12)

where

$\nu=\frac{q}{(1+q)^{2}}.$ (13)

Eq. (9) is the equation we use to produce the rotation predictions in the
right panel of Fig. 1 in the main text. Note that one can imagine that a
rotated $(2,0)$ perturbation should also contribute to the spectrum across
$\theta$. This effect, however, should be subdominant compared to those of the
rotated $(2,2)$ and $(2,1)$ perturbations so we do not consider it in this
work.

## References

  * Teukolsky [1973] S. A. Teukolsky, Astrophys. J. 185, 635 (1973).
  * Buonanno  _et al._ [2007] A. Buonanno, G. B. Cook, and F. Pretorius, Phys. Rev. D 75, 124018 (2007), arXiv:gr-qc/0610122 .
  * Owen [2009] R. Owen, Phys. Rev. D 80, 084012 (2009), arXiv:0907.0280 [gr-qc] .
  * Andersson [1999] N. Andersson, Journal of Astrophysics and Astronomy 20, 269 (1999).
  * Frolov and Novikov [1998] V. P. Frolov and I. D. Novikov, eds.,  _Black hole physics: Basic concepts and new developments_ (1998).
  * Vishveshwara [1970] C. V. Vishveshwara, Phys. Rev. D 1, 2870 (1970).
  * Press [1971] W. H. Press, Astrophys. J. Lett. 170, L105 (1971).
  * Chandrasekhar and Detweiler [1975] S. Chandrasekhar and S. L. Detweiler, Proc. Roy. Soc. Lond. A 344, 441 (1975).
  * Leaver [1985] E. W. Leaver, Proc. Roy. Soc. Lond. A 402, 285 (1985).
  * Berti and Cardoso [2006] E. Berti and V. Cardoso, Phys. Rev. D 74, 104020 (2006).
  * Kamaretsos  _et al._ [2012a] I. Kamaretsos, M. Hannam, and B. Sathyaprakash, Phys. Rev. Lett. 109, 141102 (2012a), arXiv:1207.0399 [gr-qc] .
  * Kamaretsos  _et al._ [2012b] I. Kamaretsos, M. Hannam, S. Husa, and B. S. Sathyaprakash, Phys. Rev. D 85, 024018 (2012b).
  * London [2020] L. London, Phys. Rev. D 102, 084052 (2020).
  * Li  _et al._ [2022] X. Li, L. Sun, R. K. L. Lo, E. Payne, and Y. Chen, Phys. Rev. D 105, 024016 (2022).
  * Borhanian  _et al._ [2020] S. Borhanian, K. G. Arun, H. P. Pfeiffer, and B. S. Sathyaprakash, Class. Quant. Grav. 37, 065006 (2020), arXiv:1901.08516 [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).
  * 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] .
  * Gangardt  _et al._ [2021] D. Gangardt, N. Steinle, M. Kesden, D. Gerosa, and E. Stoikos, Phys. Rev. D 103, 124026 (2021).
  * Apte and Hughes [2019] A. Apte and S. A. Hughes, Phys. Rev. D 100, 084031 (2019).
  * Lim  _et al._ [2019] H. Lim, G. Khanna, A. Apte, and S. A. Hughes, Phys. Rev. D 100, 084032 (2019), arXiv:1901.05902 [gr-qc] .
  * Hughes  _et al._ [2019] S. A. Hughes, A. Apte, G. Khanna, and H. Lim, Phys. Rev. Lett. 123, 161101 (2019), arXiv:1901.05900 [gr-qc] .
  * Ghosh  _et al._ [2023] S. Ghosh, P. Kolitsidou, and M. Hannam, arXiv eprints (2023), arXiv:2310.16980 [gr-qc] .
  * Schmidt  _et al._ [2012] P. Schmidt, M. Hannam, and S. Husa, Phys. Rev. D 86, 104063 (2012), arXiv:1207.3088 [gr-qc] .
  * O’Shaughnessy  _et al._ [2013] R. O’Shaughnessy, L. London, J. Healy, and D. Shoemaker, Phys. Rev. D 87, 044038 (2013), arXiv:1209.3712 [gr-qc] .
  * Boyle  _et al._ [2011] M. Boyle, R. Owen, and H. P. Pfeiffer, Phys. Rev. D 84, 124011 (2011).
  * Pratten  _et al._ [2021] G. Pratten _et al._ , Phys. Rev. D 103, 104056 (2021), arXiv:2004.06503 [gr-qc] .
  * Hamilton  _et al._ [2023] E. Hamilton, L. London, and M. Hannam, arXiv (2023), arXiv:2301.06558 [gr-qc] .
  * Giesler  _et al._ [2019] M. Giesler, M. Isi, M. A. Scheel, and S. Teukolsky, Phys. Rev. X 9, 041060 (2019), arXiv:1903.08284 [gr-qc] .
  * Finch and Moore [2021] E. Finch and C. J. Moore, Phys. Rev. D 103, 084048 (2021), arXiv:2102.07794 [gr-qc] .
  * Magaña Zertuche  _et al._ [2022] L. Magaña Zertuche _et al._ , Phys. Rev. D 105, 104015 (2022), arXiv:2110.15922 [gr-qc] .
  * Mitman  _et al._ [2021] K. Mitman _et al._ , Phys. Rev. D 104, 024051 (2021), arXiv:2105.02300 [gr-qc] .
  * Mitman  _et al._ [2022] K. Mitman _et al._ , Phys. Rev. D 106, 084029 (2022), arXiv:2208.04356 [gr-qc] .
  * Aasi  _et al._ [2015] J. Aasi _et al._ (LIGO Scientific), Class. Quant. Grav. 32, 074001 (2015), arXiv:1411.4547 [gr-qc] .
  * Acernese  _et al._ [2015] F. Acernese _et al._ (VIRGO), Class. Quant. Grav. 32, 024001 (2015), arXiv:1408.3978 [gr-qc] .
  * Akutsu  _et al._ [2021] T. Akutsu _et al._ (KAGRA), PTEP 2021, 05A101 (2021), arXiv:2005.05574 [physics.ins-det] .
  * Boyle  _et al._ [2019] M. Boyle _et al._ , Class. Quant. Grav. 36, 195006 (2019), arXiv:1904.04831 [gr-qc] .
  * [37] https://github.com/HengruiPrinceton/precession_ringdown.
  * Isi and Farr [2021] M. Isi and W. M. Farr, arXiv eprints (2021), arXiv:2107.05609 [gr-qc] .
  * Deppe  _et al._ [2023] N. Deppe, W. Throwe, L. E. Kidder, N. L. Vu, K. C. Nelli, C. Armaza, M. S. Bonilla, F. Hébert, Y. Kim, P. Kumar, G. Lovelace, A. Macedo, J. Moxon, E. O’Shea, H. P. Pfeiffer, M. A. Scheel, S. A. Teukolsky, N. A. Wittek, _et al._ , SpECTRE v2023.10.11, 10.5281/zenodo.8431874 (2023).
  * Moxon  _et al._ [2020] J. Moxon, M. A. Scheel, and S. A. Teukolsky, Phys. Rev. D 102, 044052 (2020), arXiv:2007.01339 [gr-qc] .
  * Moxon  _et al._ [2023] J. Moxon, M. A. Scheel, S. A. Teukolsky, N. Deppe, N. Fischer, F. Hébert, L. E. Kidder, and W. Throwe, Phys. Rev. D 107, 064013 (2023), arXiv:2110.08635 [gr-qc] .
  * Boyle  _et al._ [2020] M. Boyle, D. Iozzo, and L. C. Stein, moble/scri: v1.2 (2020).
  * Boyle [2013] M. Boyle, Phys. Rev. D 87, 104006 (2013).
  * Boyle  _et al._ [2014] M. Boyle, L. E. Kidder, S. Ossokine, and H. P. Pfeiffer, Gravitational-wave modes from precessing black-hole binaries (2014), arXiv:1409.4431 [gr-qc] .
  * Boyle [2016a] M. Boyle, Phys. Rev. D 93, 084031 (2016a), arXiv:1509.00862 [gr-qc] .
  * Berti and Klein [2014] E. Berti and A. Klein, Phys. Rev. D 90, 064012 (2014), arXiv:1408.1860 [gr-qc] .
  * Zhu  _et al._ [2023] H. Zhu, J. L. Ripley, A. Cárdenas-Avendaño, and F. Pretorius, arXiv eprints (2023), arXiv:2309.13204 [gr-qc] .
  * Ruiz  _et al._ [2008] M. Ruiz, R. Takahashi, M. Alcubierre, and D. Nunez, Gen. Rel. Grav. 40, 2467 (2008), arXiv:0707.4654 [gr-qc] .
  * Iozzo  _et al._ [2021] D. A. B. Iozzo _et al._ , Phys. Rev. D 103, 124029 (2021), arXiv:2104.07052 [gr-qc] .
  * Boyle [2016b] M. Boyle, Journal of Mathematical Physics 57, 092504 (2016b), https://pubs.aip.org/aip/jmp/article-pdf/doi/10.1063/1.4962723/13685223/092504_1_online.pdf .
  * Wigner [2012] E. Wigner, Group theory: And its application to the quantum mechanics of atomic spectra (2012).
  * Schmidt  _et al._ [2011] P. Schmidt, M. Hannam, S. Husa, and P. Ajith, Phys. Rev. D 84, 024046 (2011).
  * Press and Teukolsky [1973] W. H. Press and S. A. Teukolsky, Astrophys. J. 185, 649 (1973).
  * Apostolatos  _et al._ [1994] T. A. Apostolatos, C. Cutler, G. J. Sussman, and K. S. Thorne, Phys. Rev. D 49, 6274 (1994).
  * Kidder [1995] L. E. Kidder, Phys. Rev. D 52, 821 (1995).
  * O’Shaughnessy  _et al._ [2011] R. O’Shaughnessy, B. Vaishnav, J. Healy, Z. Meeks, and D. Shoemaker, Phys. Rev. D 84, 124002 (2011), arXiv:1109.5224 [gr-qc] .
  * Bonnor and Rotenberg [1961] W. B. Bonnor and M. A. Rotenberg, Proceedings of the Royal Society of London Series A 265, 109 (1961).
  * Brügmann  _et al._ [2008] B. Brügmann, J. A. González, M. Hannam, S. Husa, and U. Sperhake, Phys. Rev. D 77, 124047 (2008).
  * Ma  _et al._ [2021] S. Ma, M. Giesler, V. Varma, M. A. Scheel, and Y. Chen, Phys. Rev. D 104, 084003 (2021), arXiv:2107.04890 [gr-qc] .
  * Alessio and Arzano [2019] F. Alessio and M. Arzano, Journal of High Energy Physics 2019, 28 (2019).
  * 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  _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] .
  * Abbott  _et al._ [2020] R. Abbott _et al._ (LIGO Scientific, Virgo), Astrophys. J. Lett. 900, L13 (2020), arXiv:2009.01190 [astro-ph.HE] .
  * 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] .
  * Abbott  _et al._ [2021] R. Abbott _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. D 103, 122002 (2021).
  * Ossokine  _et al._ [2020] S. Ossokine, A. Buonanno, S. Marsat, R. Cotesta, S. Babak, T. Dietrich, R. Haas, I. Hinder, H. P. Pfeiffer, M. Pürrer, C. J. Woodford, M. Boyle, L. E. Kidder, M. A. Scheel, and B. Szilágyi, Phys. Rev. D 102, 044055 (2020).
  * Waskom [2021] M. L. Waskom, Journal of Open Source Software 6, 3021 (2021).
  * Hunter [2007] J. D. Hunter, Computing in Science & Engineering 9, 90 (2007).
  * Kluyver  _et al._ [2016] T. Kluyver _et al._ , Jupyter notebooks – a publishing format for reproducible computational workflows (2016).
  * Harris  _et al._ [2020] C. R. Harris _et al._ , Nature 585, 357 (2020).
  * Van Rossum and Drake [2009] G. Van Rossum and F. L. Drake,  _Python 3 Reference Manual_ (CreateSpace, Scotts Valley, CA, 2009).


# Ultrafast generation and detection of propagating coherent acoustic phonon
wave packets in ultra-thin iron pnictide films

D. Cheng, B. Song  Department of Physics and Astronomy, Iowa State University,
and Ames Laboratory, Ames, IA 50011 USA  J. H. Kang  C. Sundahl  Department of
Materials Science and Engineering, University of Wisconsin-Madison, Madison,
WI 53706, USA  L. Luo, J-M. Park  Department of Physics and Astronomy, Iowa
State University, and Ames Laboratory, Ames, IA 50011 USA  Y. G. Collantes  E.
E. Hellstrom  Applied Superconductivity Center, National High Magnetic Field
Laboratory, Florida State University, Tallahassee, FL 32310, USA  M. Mootz  I.
E. Perakis  Department of Physics, University of Alabama at Birmingham,
Birmingham, AL 35294-1170, USA  C. B. Eom  Department of Materials Science and
Engineering, University of Wisconsin-Madison, Madison, WI 53706, USA  J. Wang
<EMAIL_ADDRESS><EMAIL_ADDRESS>Department of Physics and Astronomy,
Iowa State University, and Ames Laboratory, Ames, IA 50011 USA

###### Abstract

We observe pronounced oscillations in differential reflectivity of 9 nm and 60
nm BaFe2As2 (Ba-122) thin films using ultrafast optical spectroscopy. Our
studies show that the oscillations result from propagating longitudinal
acoustic (LA) phonon wave packets with strong thickness and temperature
dependence. Particularly, the experimentally measured oscillation frequency
approaches to 50 GHz for the ultra-thin film. Our calculations show that
Young’s modulus of 9 nm thin film is nearly four times as large as that of 60
nm thin film, consistent with the experiment. The increase in Young’s modulus
as thickness decrease was attributed to the decrease in parent Ba-122
tetragonality $c/a$ near the film-substrate interface due to material-
substrate mismatch effect. The temperature-dependent change in LA phonon
frequency was attributed to the change in parent Ba-122 othorhombicity
$(a-b)/(a+b)$.

††preprint: AIP/123-QED

The high-frequency phonons carry key information about the materials’ elastic
properties and were commonly detected by superconducting tunnel junctions and
heat pulsesreview . Both measurements were performed at low temperatures
(usually in superconductivity regimes). Thanks to the development of ultrafast
optical techniques, the modern approach is based on ultrafast lasers, which
can generate and detect acoustic phonon wave packets, as seen in the
pioneering work in 1984 by Thomsen et alfirstexp . After the early attempts,
enduring efforts have been devoted to generate and study phonon modes using
ultrafast lasersiron ; ref1 ; ref2 ; ref4 ; ref5 . The intense ultrafast
optical pulse (pump) disturbs the thin film and generates acoustic phonons,
which will bounce back and forth between film surfaces and cause a periodic
modulation of optical properties. The modulation is subsequently measured by a
time-delayed light pulse (probe) at sample’s surface. The detected phonon mode
is usually longitudinal acoustic (LA) waves since the transverse acoustic (TA)
waves, or shear waves, could only be generated through anisotropic thermal
expansion in thermally anisotropic materialsshearmode . Unlike superconducting
tunnel junctions and heat pulses techniques, the ultrafast laser method can
also be used at much higher temperaturesthermal ; a ; b ; c ; d ; e .

The ultrafast optical technique expands our knowledge beyond bulks, with the
unprecedented sensitivity, time resolution and applicability in ultra-thin
films even down to $\sim$10 nm. Specifically, we can examine the elastic
properties, like Young’s modulus, of ultra-thin films compared with its bulk
counterpart. The properties might have been significantly altered due to the
strain effect and lower mobility of dislocation. Such studies and pertinent
discussions Young ; cc  have been rare, especially for ultra-thin film down to
10nm and temperature dependence of Young’s modulus, which has motivated the
present study on BaFe2As2 (Ba-122) films.

The Ba-122 phase of iron pnictide materials hosts a series of layered
intermetallics HC ; 122_1 ; 122_2 ; 122_3 ; 122_4 ; 122_5  and their
derivatives bqs , subject to a tetragonal-orthohombic transition and a phase
diagram involving magnetism and superconductivity David . It has been shown
that superconductivity of ultra-thin film Ba-122 could be enhanced by tuning
the othorohmbicity $(a-b)/(a+b)$ of the unit cellsamplegrower . Change of
sample’s othorohmbicity would yield the change of elastic properties like
intrinsic strain, indicative of the fact that the Young’s modulus of ultra-
thin film used for the study should differ from that of the bulk. This makes
the system a good candidate to study the anomalous elastic properties in thin
films.

In this Letter, we demonstrate a powerful and versatile method to measure
Young’s modulus of Ba-122 ultra-thin films through coherent LA phonon
spectroscopy. This is implemented by ultrafast two-color differential
reflectivity spectroscopy utilizing the pump (800 nm, 1.55eV) and probe (400
nm, 3.1 eV) laser pulses. Two thin films of 9 nm and 60 nm are measured and
compared. We also performed temperature-dependent phonon frequency measurement
to study its connections to the structural phase transition. We found that
Young’s modulus of Ba-122 film increase with the decrease in film thickness,
which was mainly attributed to the material-substrate mismatch effect, i.e.,
tetragonality $c/a$ change near the film-substrate interface. The temperature-
dependent change in LA phonon frequency was attributed to the change in parent
Ba-122 othorhombicity $(a-b)/(a+b)$.

Parent Ba-122 thin films were grown on single-crystal cubic substrates by
pulsed laser deposition with a KrF (248 nm) ultraviolet excimer laser at
650-740 °C. High-quality epitaxial films were synthesized on lithium fluoride
(LiF). The 800 nm laser is directly from a Ti:Sa amplifier with 40 fs pulse
duration and 1 kHz repetition rate, which is used as pump to photoexcite the
samples out of equilibrium. The 400 nm laser is frequency-doubled from the 800
nm beam by using a Beta Barium Borate (BBO) crystal, which is used as a probe
to detect the transient reflectivity change of the samples. The pump-probe
time delay is controlled by an optical delay stage. Both pump and probe beams
are at near normal incidence of the samples.

Figures 1a and 1b illustrate the generation and detection methods of LA in
Ba-122 thin films and the crystal structure of bulk Ba-122 and LiF. Figure 1c
shows the schematics of $a$-$c$ over layers for Ba-122 thin film. The lattice
constants are different at the top and the bottom (in contact with substrate)
due to the intrinsic stress, more specifically the material-substrate
mismatch. The lattice constants near the top are $a=3.9612$ Å, and $c=13.0061$
Å. The lattice constant a of Ba-122 single crystal is smaller than that of LiF
substrate (4.027 Å), which causes larger lattice constant $a^{\prime}$ and
smaller lattice constant $c^{\prime}$ (thus smaller tetragonality
$c^{\prime}$/$a^{\prime}$) near the bottom. Generally, the larger the number
of layers, the less the material-substrate mismatch effect.

Figure 1: (a) Crystal structures of bulk Ba-122 and LiFcrystalstructure . (b)
Generation and detection of coherent acoustic phonons in Ba-122 thin films.
The coherent phonons generated by the 800 nm pump pulse can be detected by the
time-delayed 400 nm probe pulse through the photoelastic effect. (c)
Schematics of $a$-$c$ over layers for Ba-122 thin film.

Figures 2a and 2b show the transient reflectivity $\Delta R/R$ of 9 nm and 60
nm Ba-122 thin films, respectively, at $T$=5 K. Besides exponential decay of
$\Delta R/R$ signal, periodic oscillations in $\Delta R/R$ were observed in
both samples. The acoustic pulse generated at the surface propagates through
the oscillations of thin film’s sensitivity function and leads to the observed
oscillations. We did theoretical fitting of the $\Delta R/R$ data and found
that they cannot be fitted with a single exponential decay and that a
biexponential decay works pretty well, which can be described as $\Delta
R/R=R\textsubscript{0}+A\textsubscript{1}\text{exp}^{-(t-t_{0})/\tau\textsubscript{1}}+A\textsubscript{2}\text{exp}^{-(t-t_{0})/\tau\textsubscript{2}}$,
where $R\textsubscript{0}$ is a time independent constant,
$A\textsubscript{1}$ and $A\textsubscript{2}$ are the amplitudes of each
exponential term, $t$ is the time delay, and $\tau\textsubscript{1}$ and
$\tau\textsubscript{2}$ are the decay time for the first and second
exponential terms, respectively. According to our fitting,
$\tau\textsubscript{1}$ is 25 ps (17 ps) and $\tau\textsubscript{2}$ is 1 ns
(2 ns) for 9 nm (60 nm) thin film, respectively. After subtracting the fitted
biexponential decay background from $\Delta R/R$ signal, we performed Fourier
transform of the residual oscillation signal, which are shown in Fig. 2c and
2d, respectively. For the 9 nm film, there is one peak at $\sim$47.5 GHz (Fig.
2c) and for the 60 nm film, the peak is at $\sim$25 GHz (Fig. 2d). We
attribute the 47.5 GHz and 25 GHz resonances to the coherent LA phonon modes
generated by the pump beam since TA can only be generated in materials with
anisotropic thermal expansion coefficients.

The acoustic phonon mode resonance f is related to the sound velocity and
probe wavelength through $f=2nu\text{cos}(\theta)/\lambda$ (Ref. 10), with
$\theta$ being the angle of incidence of the probe beam (in our experiment the
$\theta$ is close to zero) and $\lambda$ being the probe wavelength (400 nm).
The refractive index $n$ for Ba-122 at 400 nm (25000 cm-1) is 3.875. From f
which has been measured in Fig. 2c and 2d, we can calculate the sound velocity
u, which is 2456.77 m/s for the 9 nm film and is 1290.32 m/s for the 60 nm
film. The sound velocity of the 60 nm film is nearly half of that of the 9 nm
film and it is related to Young’s modulus through the equation
$u=\sqrt{Y/\rho}$, with $Y$ being the Young’s modulus of the sample and $\rho$
the density of the sample. The density of the two thin films should be very
similar. Therefore, the main difference in calculated sound velocity must be
from the difference in Young’s modulus of the two samples. At fixed density
$\rho$, $Y\propto u^{2}$, thus Young’s modulus of 9 nm Ba-122 thin film is
nearly four times as large as that of 60 nm thin film. The cause will be
discussed later.

We also measured the LA phonon mode resonance as a function of temperature.
Figure 2e shows the transient reflectivity $\Delta R/R$ of the 9 nm film at
various temperatures and Fig. 2f plots the LA phonon mode resonance as a
function of temperature, extracted from Figure 2e. As seen from Fig. 2f, the
LA phonon mode resonance is nearly unchanged from 5 K to 50 K. At 80 K, the
resonance starts to drop. This trend is perfectly correlated with temperature-
dependent othorhombicity of the 9 nm Ba-122 samplesamplegrower . Therefore, we
attribute the decrease in the LA mode resonance to the decrease in Young’s
modulus originating from the decrease of orthorhombicity $(a-b)/(a+b)$.

Figure 2:  Transient reflectivity change ${\Delta}R/R$ with corresponding
biexponential fitting and residuals of (a) 9 nm and (b) 60 nm Ba-122 thin film
at 5 K. Data are offset for clarity. Fourier spectra for the residuals of (c)
9 nm and (d) 60 nm Ba-122 thin film transient reflectivity. (e) Temperature-
dependent transient reflectivity ${\Delta}R/R$ change (with offset) and (f)
acoustic phonon frequency change for 9 nm Ba-122 thin film.

To understand $u$ at different film thicknesses, we performed density
functional theory (DFT) calculations Perdew  with a set of tested parameters
CalPara . Note that the sound propagates perpendicular to the film surface and
thus mainly depends on the elastic properties along $c$-axis. On the other
hand, the strain caused by the mismatch of substrate LiFe ($a=4.027$ Å,) and
Ba-122 ($a=3.9612$ Å)Song  mainly affects the $a$-$b$ plane, i.e., the
tetragonal unit cell tends to be more expanded as it gets closer to the
substrate, which is subject to intense distortion. The 9 nm film corresponds
to 6 layers of the unit cell (Fig. 1b); while the 60 nm film corresponds to 40
layers, which is close to a bulk state. With this insight, we model the
influence of film thickness by the extension of distortion in the $a$-$b$
plane (with tetragonal symmetry). In other words, we try to address the
drastic change in $u$ by an elastic effect. Our purpose is to examine whether
the elastic distortion can give rise to the significant changes observed in
experiments.

The lattice parameter $a$ of substrate LiF sets an upper limit for distortion.
That is, we examine $a$ in a range from 3.9612 Å (bulk limit, 60 nm film) to
4.027 Å (distortion limit, 9 nm film). Figure 3a illustrates the lattice
parameter $c$ decreases with parameter $a$. Meanwhile, the cell volume $V$ is
not constant, but will increase with $a$ (Fig. 3b). That suggests the bottom
layer generally has a larger unit cell compared with the top (bulk). Elastic
property, like Young’s modulus, can be examined, which is expressed as Young

$Y=\frac{1}{V}\frac{{\partial}^{2}{E}}{\partial\varepsilon^{2}}.$ (1)

Volume $V$ is for the unit cell of the given distorted lattice $a$; $E$ is the
total energy, $\varepsilon$ is the strain, which specially refers to
$\varepsilon_{cc}$ in this scenario (the diagonal term along $c$-axis). With
DFT and unit cell method, we can evaluate
${\partial}^{2}{E}/{\partial\varepsilon^{2}}$ by a fixed $a$ (maintain
symmetry) and a slightly varied $c$ around the equilibrium $c_{eq}$ at given
$a$.

$\varepsilon_{cc}={\Delta}c/d,~{}{\Delta}c=c-c_{eq}(a)$ (2)

where $d$ is the thickness of thin film. The results are shown in Figure 3c
and 3d. Notably, with increased $a$ and $V$, Young’s modulus undergoes a four
times increase (Fig. 3c). Combined with

$u=\sqrt{Y/{\rho}}=\sqrt{YV/m},$ (3)

we find that for bulk sample, the sound speed is 1240 m/s and for thin limit
sample, the sound speed is 2416 m/s (Fig. 3e). Our theoretical calculations
are in good agreement with the experimental results.

Note that $\varepsilon$ is evaluated for the whole film, rather than for
different locations within a film. Thus, the evaluated $Y$ physically means
the average value for the whole thin film. It is important to realize one
fitting curve (parabolic) in Fig. 3d will only correspond to one point in Fig.
3c. For example, in Fig. 3d, the red fitting curve corresponds to the bulk
limit (red dashed line in Fig. 3c), and the points represent a series of small
deviations around bulk; while the blue curve corresponds to the highly
distorted limit (blue dashed lines in Fig. 3c).

Figure 3: (a) Lattice parameter $c$ as a function of lattice parameter $a$.
(b) Cell volumes as a function of lattice $a$. (c) Young’s modulus dependence
of $a$. (d) Elastic energy function E for 9 nm thin film (red dots, solid blue
line is fitting) and bulk (red dots, solid red line is fitting). (e) Sound
speeds as a function of lattice parameter a. (f) Energy as a function of
strain $\epsilon$ when $a$ = 3.99 Å(red dots, solid black line is fitting).

Noteworthy, discontinuity is observed when the lattice parameters are varied
(Fig. 3a-c, e). This discontinuity corresponds to a phase transition, which is
relevant to the double (or multiple) energy valleysPbond , which have been
observed in layer-stacking structures, e.g., Ba-122/Fe based 1144-phaseHC . It
occurs when $c$ parameter decreases, either for an applied pressure, or for an
extended $a$-$b$ plane as done in this case. The $c$ lattice parameter will
undergo a drastic change at a critical value, and it is called the collapse
phase HC . Consequently, from Fig. 3f, the strain energy function near the
discontinuity is obviously deviant from parabolic line. That means $Y$ will
exhibit quite different values at the two sides. It is interesting to check
such abnormality using nano-scales imaging JW1  and other THz spectroscopy
techniques XYang ; JW2 ; JW3 . In addition, although the mechanical properties
are not the most intriguing aspects for intermetalic materials like Ba-122,
the demonstrated ultrafast technique allows the sensitive measurement of $Y$
and $u$, which could provide a new avenue to obtain key information about the
rich phase transition, electron-lattice coupling, etc.

To summarize, we discovered propagating coherent acoustic phonon wavepackets
in Ba-122 thin films using ultrafast transient reflectivity measurement. We
reveal clear difference in frequency of the excited LA mode between 9 nm and
60 nm thin film sample. For the 9 nm thin film the LA frequency is nearly
twice as large as that of the 60 nm samples. Temperature dependent LA mode
frequency follows the structural phase transition of the sample. We attribute
the change in mode frequency to the change in othorhombicity, which leads to
the change in Young’s modulus, consistent with our calculations. Our results
show that ultrafast transient reflectivity technique represents a powerful and
sensitive tool to study the mechanical properties and the structural phase
transition of ultrathin films with thickness even less than 10 nm. Ultrast
acoustic phonon detection may be extended to diverse complex thin films
including superconducting S_1 ; S_2 , magnetic M_1  and photovoltaic materials
P_1 .

###### Acknowledgements.

This work was supported by National Science Foundation 1905981 (spectroscopy
and analysis). L.L. was supported by the Ames Laboratory, the US Department of
Energy, Office of Science, Basic Energy Sciences, Materials Science and
Engineering Division under contract No. DEAC0207CH11358 (technical
assistance). The work at UW-Madison (synthesis and characterizations of
epitaxial thin films) was supported by the US Department of Energy (DOE),
Office of Science, Office of Basic Energy Sciences (BES), under award number
DE-FG02-06ER46327 (C.B.E.). Theory work at the University of Alabama,
Birmingham (modeling) was supported by the US Department of Energy under
contract # DE-SC0019137 (M.M and I.E.P) and was made possible in part by a
grant for high performance computing resources and technical support from the
Alabama Supercomputer Authority.

## References

  * (1) M.N. Wybourne, J.K. Wigmore, and N. Perrin, J. Phys. Condens. Matter. 32, 5347 (1989).
  * (2) C. Thomsen, J. Strait, Z. Vardeny, H. J. Maris, J. Tauc, and J. J. Hauser, Phys. Rev. Lett. 53, 989 (1984).
  * (3) S. Kumar, L. Harnagea, S. Wurmehl, B. Buchner, and A. K. Sood, EPL 100, 57007 (2012).
  * (4) J. Wang, Y. Hashimoto, J. Kono, A. Oiwa, H. Munekata, G. D. Sanders, and C. J. Stanton, Phys. Rev. B 72, 153311 (2005).
  * (5) P. Ruello, T. Pezeril, S. Avanesyan, G. Vaudel, V. Gusev, I. C. Infante, and B. Dkhil, Appl. Phys. Lett. 100, 212906 (2012).
  * (6) B. Mansart, D. Boschetto, A. Savoia, F. Rullier-Albenque, A. Forget, D. Colson, A. Rousse, and M. Marsi, Phys. Rev. B 80, 085210 (2009)
  * (7) H. Takahashi, Y. Kamihara, H. Koguchi, T. Atou, H. Hosono, T. Katayama, J. Takeda, M. Kitajima, and K. G. Nakamura J. Phys. Soc. Jpn. 80 013707 (2011).
  * (8) O. Matsuda, O. B. Wright, D. H. Hurley, V. E. Gusev, and K. Shimizu, Phys. Rev. Lett. 93, 095501 (2004).
  * (9) Vaswani C et. al., Phys. Rev. X 10, 021013 (2020).
  * (10) Luo L et. al., Nat. Mater. 20, 329 (2021).
  * (11) Yang X et. al., npj Quantum Mater. 5, 13 (2020)
  * (12) GD Sanders, CJ Stanton, J Wang, J Kono, A Oiwa, H Munekata, Physical Review B  72, 245302 (2005).
  * (13) J Wang, GA Khodaparast, J Kono, A Oiwa, H Munekata, Journal of Modern Optics 51, 2771 (2004).
  * (14) Li T, et. al., Nature 496, 69 (2013).
  * (15) J. Wang, and F. G. Shi, Scr. mater. 42, 687-694 (2000).
  * (16) H. Ogi, M. Fuji, N. Nakamura, T. Yasui, and M. Hirao, Phys. Rev. Lett. 98, 195503 (2007).
  * (17) B. Q. Song, Mingyu Xu, Vladislav Borisov, Olena Palasyuk, C. Z. Wang, Roser Valenti, Paul C Canfield, K. M. Ho, Phys. Rev. Materials 5, 094802 (2021)
  * (18) Patz A et. al., Nat. Commun., 5, 3229 (2014).
  * (19) Patz A et. al., Phys. Rev. B., 95, 165122 (2017).
  * (20) Yang X et. al., Phys. Rev. Lett., 121, 267001 (2018).
  * (21) Yang X et. al., Phys. Rev. Lett., 121, 267001 (2018).
  * (22) C. Vaswani et. al., Nat. Commun., 12, 258 (2021).
  * (23) B. Q. Song, Manh Cuong Nguyen, C. Z. Wang, P. C. Canfield, and K. M. Ho, Phys. Rev. Materials 2, 104802 (2018).
  * (24) David C. Johnston, Adv. Phys. 59, 803 (2010).
  * (25) J. H. Kang, P. J. Ryan, J.-W. Kim, J. Schad, J. Podkaminer, N. Campbell, J. Suttle, T. H. Kim, L. Luo, D. Cheng, J. Wang, R. McDermott, M. S. Rzchowski, C. B. Eom (to be published)
  * (26) K. Momma and F. Izumi, J. Appl. Crystallogr., 44, 1272-1276 (2011).
  * (27) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * (28) All calculations of energy are based on density functional theory and Perdew Burke Ernzerhof (PBE) correlation functional, implemented by the Vienna ab initio simulation package (VASP). The projected augmented wave (PAW) pseudopotential method is employed. The plane wave energy truncation is 400 eV. The k-sampling has also been tested to ensure the convergence of energy. The magnetism is stripe AFM for Fe atoms. See Ref. Song ; bqs  for a detailed description in various 122-phase compounds. Such magnetic structures have enlarged the unit cell, we have adopted the minimum sufficient $\sqrt{2}{\times}\sqrt{2}{\times}1$ 122-unit cell that contains 20 atoms.
  * (29) B. Q. Song, Manh Cuong Nguyen, C. Z. Wang, and K. M. Ho Phys. Rev. B, 97, 094105 (2018).
  * (30) R. Hoffmann and C. Zheng, J. Phys. Chem. 89, 4175 (1985).
  * (31) Kim R H et. al., ACS Photonics 8, 1873 (2021).
  * (32) X. Yang, B. Song, C. Vaswani, L. Luo, C. Sundahl, M. Mootz, J-H. Kang, Y. Yao, K-M Ho, I. E. Perakis, C. B. Eom and J. Wang, ArXiv 2008.08732 (2020).
  * (33) Luo L et. al., Nat. Commun. 10, 607 (2019).
  * (34) Yang X et. al., Phys. Rev. B, 99, 094504 (2018).
  * (35) Vaswani C et al., Phys. Rev. Lett. 124, 207003 (2020)
  * (36) X. Yang et al., Phys. Rev. B. 99, 094504 (2019)
  * (37) A Patz et al., Physical Review B. 91, 155108 (2015)
  * (38) Liu Z et al., Phys. Rev. Lett. 124, 157401 (2020)


# Continual Learning for Infinite Hierarchical
Change-Point Detection

###### Abstract

Change-point detection (CPD) aims to locate abrupt transitions in the
generative model of a sequence of observations. When Bayesian methods are
considered, the standard practice is to infer the posterior distribution of
the change-point locations. However, for complex models (high-dimensional or
heterogeneous), it is not possible to perform reliable detection. To
circumvent this problem, we propose to use a hierarchical model, which yields
observations that belong to a lower-dimensional manifold. Concretely, we
consider a latent-class model with an unbounded number of categories, which is
based on the chinese-restaurant process (CRP). For this model we derive a
continual learning mechanism that is based on the sequential construction of
the CRP and the expectation-maximization (EM) algorithm with a stochastic
maximization step. Our results show that the proposed method is able to
recursively infer the number of underlying latent classes and perform CPD in a
reliable manner.

Index Terms—  Bayesian inference, continual learning, change-point detection
(CPD), chinese-restaurant process (CRP), expectation-maximization (EM)
algorithm.

##  1 Introduction

Change-point detection (CPD), which consists of locating abrupt transitions in
the generative model of the observations, is a problem with a plethora of
applications. For instance, CPD is widely used in finance [1, 2], the analysis
of social networks [3, 4], or cognitive radio [5, 6]. The main focus of CPD
methods has been traditionally on batch settings, where the entire sequence of
observations is available and has to be segmented. However, CPD is most useful
in online scenarios, where change points must be detected as new incoming
samples are observed. Online CPD methods have two intertwined tasks to solve:
i) segmentation of sequential data into partitions (or segments) and ii)
estimation of the generative model parameters for the given partitions.

Since each partition has a different generative distribution, the
identifiability of change points is related to the difference between such
distributions. In this context, Bayesian inference is useful for inferring the
distributions given a prior distribution in a reliable manner. The Bayesian
online change-point detection (BOCPD) approach [7] used this idea for
recursively performing density estimation, which yields a more robust
detection process as the propagation of uncertainty is considered. However, it
can be observed that, for complex likelihood models, which have a number of
parameters much higher than the number of observations between two consecutive
change points, reliable CPD becomes unfeasible. This can be the case of,
although is not restricted to, high-dimensional and/or heterogeneous
observations (mixture of continuous and discrete variables), which usually
have a prohibitive number of parameters.

To address the aforementioned issue, in [8] we presented a hierarchical
probabilistic model based on latent classes, i.e., a mixture model. The CPD
problem can be carried out directly on the lower-dimensional manifold, where
the discrete latent variables lie. Hence, this method requires less evidence
than the observational counterpart since the number of parameters is reduced,
which yields faster and more reliable detections. However, [8] requires that
the number of classes is fixed _a priori_.

The main contribution of this paper is to introduce a novel approach, based on
continual learning [9, 10, 11], to recursively infer the underlying sequence
of latent classes, its distributions, and the change points. The key idea of
the proposed model is to allow for an unbounded order on the latent model,
that is, the number of classes is not fixed and could even become infinite. In
particular, we use the Chinese-restaurant process (CRP) [12], which is a well-
known Bayesian non-parametrics method, to model the latent variables with an
unbounded number of classes. That is, the CRP may increase the number of
classes as new observations come in. Moreover, as with any mixture model, the
expectation-maximization (EM) algorithm [13] is used, but in this work the
maximization step (M-step) is substituted by a stochastic M-step [14].
Finally, the experimental results on real data show how both the latent-class
inference process and the change-point detection perform reliably.

##  2 Bayesian Online Change-Point Detection

We start by considering a time series
$\bm{x}_{1:t}=\\{x_{1},x_{2},\ldots,x_{t}\\}$, which is divided into non-
overlapping partitions, denoted by $\rho_{i},i=1,2,\ldots$ Each partition is
separated from its neighbors by change points (cp). Based on [7], we assume
that the data within each partition $\rho_{i}$ is independent and identically
distributed (i.i.d.) according to some generative probability distribution
$p(x_{t}|\bm{\theta}_{\rho_{i}})$, where the parameter vector,
$\bm{\theta}_{\rho_{i}}$, is unknown. Under this assumption, change points are
determined by changes in the parameters:

$\bm{\theta}_{t}=\begin{cases}\bm{\theta}_{\rho_{1}},&t<\textsc{cp}_{1},\\\
\bm{\theta}_{\rho_{2}},&\textsc{cp}_{1}\leq t\leq\textsc{cp}_{2},\\\
\bm{\theta}_{\rho_{3}},&\textsc{cp}_{2}\leq t\leq\textsc{cp}_{3},\\\
&\vdots\end{cases}$ (1)

The main idea in [7] is the run-length, $r_{t}$, which is defined as a
discrete random variable that counts the number of time-steps since the last
cp, that is,

$r_{t}=\begin{cases}0,&\textsc{cp}\text{ at time }t\\\
r_{t}+1,&\text{otherwise},\\\ \end{cases}$ (2)

and may be seen as a proxy for change points. The objective of the BOCPD
technique is to compute the posterior distribution $p(r_{t}|\bm{x}_{1:t})$
recursively, from which we will identify a cp if the probability mass
accumulates near $r_{t}=0$.

The posterior distribution $p(r_{t}|\bm{x}_{1:t})$ is obtained by
marginalizing the joint distribution $p(r_{t},\bm{x}_{1:t})$ over all the
$r_{t}$ values seen so far, which, in turn, is computed by marginalizing the
model parameters, $\bm{\theta}_{t}$. The learning of $\bm{\theta}_{t}$ given
the partition, required for the computation of $p(r_{t},\bm{x}_{1:t})$, is
carried out using a multiple thread inference mechanism induced by the run-
length. For instance, to learn $\bm{\theta}_{3}$ given $r_{3}=2$, only the
observations $\\{x_{2},x_{3}\\}$ are required. This parallel inference scheme
is depicted in Figure 1, where we illustrate the aforementioned example using
the notation $\bm{\theta}_{3}|\\{x_{2},x_{3}\\}$.

The inference of $p(r_{t}|\bm{x}_{1:t})$ in [7] may become unfeasible when the
complexity of the generative model increases, for instance, for high-
dimensional and/or heterogenous observations. That is, if the likelihood
$p(x_{t}|\bm{\theta}_{\rho_{i}})$ for the partition $\rho_{i}$ depends on an
extremely large number of parameters, it would not be possible to obtain
sufficient statistical evidence to detect change points. This problem may
yield the BOCPD method unusable in some problems.

##  3 CPD on Hierarchical Models

Fig. 1: Illustration of the parallel inference threads for the estimation of
$\bm{\theta}_{t}$ conditioned on the run-length $r_{t}$ given $\bm{x}_{1:t}$.

The aforementioned problem of the BOCPD for complex generative models can be
overcome by introducing hierarchical models. We propose to use latent classes
to obtain such hierarchical model. These latent classes, $z_{t}$, yield
observations, $x_{t}$, that belong to a lower-dimensional manifold, and allow
us to write the generative distribution of $x_{t}$ as

$p(x_{t}|\bm{\theta}_{t})=\sum_{z_{t}=1}^{K}p(x_{t}|z_{t})p(z_{t}|\bm{\theta}_{t}),$

where $z_{t}$ is a categorical random variable, with $K$ being the maximum
number of classes or categories, and $\bm{\theta}_{t}$ is the vector of
parameters, i.e., the probability of each class. This form of latent-class
model can be seen as a mixture model.

Even assuming a hierarchical model, we are still interested in
$p(r_{t}|\bm{x}_{1:t})$, which would require the marginalization over
$\bm{z}_{1:t}$ as follows

$p(r_{t},\bm{x}_{1:t})=\sum_{\bm{z}_{1:t}}p(r_{t},\bm{z}_{1:t},\bm{x}_{1:t}).$
(3)

However, for large values of $t$ and $K$, the marginalization in (3) is
computationally unfeasible due to the combinatorial sums. In [8], to avoid the
marginalization, we assumed that we observe $\bm{z}_{1:t}$, instead of
marginalizing them, by directly plugging in the values of the _maximum a
posteriori_ (MAP) estimates, which are given by

$z_{t}^{\star}=\arg\max_{z_{t}}p(z_{t}|x_{t}).$ (4)

Now, using the MAP estimates as observations and assuming that the joint
distribution on the right hand side (r.h.s.) of (3) factorizes as

$p(r_{t},\bm{x}_{1:t},\bm{z}_{1:t}^{\star})=p(\bm{x}_{1:t}|\bm{z}_{1:t}^{\star})p(r_{t},\bm{z}_{1:t}^{\star}),$

with

$p(r_{t},\bm{z}_{1:t}^{\star})=\int
p(r_{t},\bm{z}_{1:t}^{\star},\bm{\theta}_{t})d\bm{\theta}_{t},$

we are effectively considering that the change points occurred on the sequence
of latent classes. Using the extended recursion of [8], which is given by

$p(r_{t},\bm{z}_{1:t}^{\star})=\sum_{r_{t-1}}p(r_{t}|r_{t-1})\Psi^{(r)}_{t}p(r_{t-1},\bm{z}_{1:t-1}^{\star}),$
(5)

where $p(r_{t}|r_{t-1})$ is the conditional prior and

$\displaystyle\Psi^{(r)}_{t}$
$\displaystyle=p(z_{t}^{\star}|r_{t-1},\bm{z}_{1:t-1}^{\star})$
$\displaystyle=\int
p(z_{t}^{\star}|\bm{\theta}_{t})p(\bm{\theta}_{t}|r_{t-1},\bm{z}_{1:t-1}^{\star})d\bm{\theta}_{t},$
(6)

is the predictive distribution of the present latent variable conditioned on
previous data and the run-length, we have all the ingredients to compute

$p(r_{t}|\bm{z}_{1:t}^{\star})=\frac{p(r_{t},\bm{z}_{1:t}^{\star})}{\sum_{r_{t}}p(r_{t},\bm{z}_{1:t}^{\star})},$
(7)

which determines the location of the change points.

###  3.1 Infinite-dimensional Hierarchical BOCPD

The problem of the hierarchical BOCPD algorithm presented above is that the
number of classes, $K$, must be known and fixed _a priori_. That is, $K$ is
not allowed to vary over time, which can be a stringent condition in some
scenarios. In this section, we consider the more interesting case that $K$ is
unknown and can be time-varying, i.e., new classes may appear as
$t\rightarrow\infty$. Then, we cannot select the order of the latent-class
model in advance. A naive idea would be to fix an upper bound on $K$ and
proceed as in the previous section. However, this upper bound could not be
available and, even if it is, the performance can be poor, as we will see in
Section 5. In the following, we will present a method for unbounded and time-
varying $K$, that is, $K$ is incremented when an unseen type of observations
appears, which translates into a hierarchical BOCPD with unbounded $K$.

Using an unbounded number of classes results in the following problem when
integrating over $\bm{\theta}_{t}$ to compute $\Psi^{(r)}_{t}$. Assuming a
Dirichlet distribution for $\bm{\theta}_{t}$, which is the conjugate prior for
categorical distributions and therefore yields a tractable integral in (6),
the evidence $p(z_{t})\rightarrow 0$ as $K$ grows. To overcome this issue, we
can consider an exchangeable distribution of the form
$p([\bm{z}_{t}])=\sum_{\bm{z}_{1:t}\in[\bm{z}_{1:t}]}p(\bm{z}_{1:t})$, where
$[\bm{z}_{1:t}]$ is a given division of classes, which is independent of the
temporal assignments, i.e., $\bm{z}_{1:3}=\\{1,2,2\\}$ corresponds to the same
division of objects as $\bm{z}_{1:3}=\\{2,1,1\\}$. This is often known as the
_exchangeability_ property [15, 12] and is a safe assumption in our setup as
we are interested in changes in the probabilities of ${\mathbf{z}}_{t}$, not
in the particular sequences $\bm{z}_{1:t}$.

The latent-class model model with an unbounded dimension can be addressed
using the Chinese-restaurant process (CRP) [12], which is a Bayesian non-
parametrics method [16]. The CRP is based on the metaphor where clients
(observations $x_{t}$) are assigned to different tables (latent classes
$z_{t}$) in a sequential manner. The assignment of classes to objects in the
CRP is determined by the predictive posterior distribution, which is given by

$p(z_{t}=k|z_{1},\ldots,z_{t-1})=\begin{cases}\frac{m_{k,t-1}}{t-1+\alpha},&k\leq
K_{t-1},\\\ \frac{\alpha}{t-1+\alpha},&k=K_{t-1}+1,\\\ \end{cases}$ (8)

where $m_{k,t-1}$ counts the number of assignments to class $k$ up to time
$t-1$, $K_{t-1}$ is the number of classes associated with $m_{k,t-1}>1$ and
$\alpha$ is a hyperparameter, which corresponds to the natural parameter of a
symmetric Dirichlet prior distribution, and controls how likely is the
appearance of a new class.

Exploiting the aforementioned CRP construction, the computation of
$\Psi^{(r)}_{t}$ in (5) is straightforward, and is given by

$\Psi^{(r)}_{t}=p(z^{\star}_{t}=k|r_{t-1},\bm{z}^{\star}_{1:t-1}),$ (9)

where we now count the number of MAP estimates, $z_{t}^{\star}$, equal to $k$
up to time $t-1$. Notice that this expression is analogous to (8) for a given
run-length, i.e., for each parallel thread in Fig. 1. Then, we may proceed to
compute the posterior $p(r_{t}|\bm{z}_{1:t}^{\star})$.

One final comment is in order. So far, we have derived a tractable recursive
way to introduce latent-class models into Bayesian CPD methods with an
unbounded number of classes. However, nothing has been said on how to compute
the MAP estimates in a continual learning fashion, which are required in (7).
This task is explored in Section 4.

##  4 Continual Learning of the CRP

Algorithm 1 Infinite-dimensional Hierarchical BOCPD

1: Input: Observe $x_{t}$ and initialize $\hat{\bm{\varphi}}_{K_{t-1}}$.

2: Sample $z_{t}\sim p(z_{t}|\bm{z}_{1:t-1}^{\star})$

3: if $z_{t}=K_{t-1}+1$ then

4: Initialize $\hat{\bm{\varphi}}_{K_{t-1}+1}$

5: end if

6: Compute $p(z_{t}=k|\bm{z}_{1:t-1}^{\star}),\forall k\leq K_{t-1}+1$

7: Compute
$\mathbb{E}[\mathbb{I}\\{z_{t}=k\\}|\bm{z}_{1:t-1}^{\star},x_{t},\hat{\bm{\varphi}}_{k}^{(t-1)}],\forall
k\leq K_{t-1}+1$

8: Update parameters $\\{\hat{\bm{\varphi}}_{k}\\}^{K_{t-1}+1}_{k=1}$ using
(11)

9: Calculate
$z^{\star}_{t}=\arg\max(p(z_{t}|\bm{z}_{1:t-1}^{\star},x_{t},\\{\hat{\bm{\varphi}}_{k}^{(t)}\\}_{k=1}^{K_{t-1}+1})$

10: if $z^{\star}_{t}=K_{t-1}+1$ then

11: $K_{t}=K_{t-1}+1$

12: end if

13: for $r_{t}=1$ to $t$ do

14: Evaluate $\Psi^{(r)}_{t}$ using (9)

15: Calculate $p(r_{t},\bm{z}_{1:t}^{\star})$

16: Obtain $p(\bm{z}_{1:t}^{\star})=\sum_{r_{t}}p(r_{t},\bm{z}_{1:t}^{\star})$

17: Compute $p(r_{t}|\bm{z}_{1:t}^{\star})$

18: Update $m^{(r)}_{k,t}\leftarrow
m^{(r)}_{k,t-1}+\mathbb{I}\\{z_{t}^{\star}=k\\}$

19: end for

20: Return: $r^{\star}_{t}=\arg\max p(r_{t}|\bm{z}_{1:t}^{\star})$

In this section, we compute the MAP estimates of $z_{t}$ in an online and
recursive fashion. This task also involves the estimation of
$\\{\bm{\varphi}_{k}\\}_{k=1}^{K_{t}}$, which are the parameters of the
mapping between observations and latent variables, that is,
$p(x_{t}|z_{t}=k)=p(x_{t}|z_{t}=k,\bm{\varphi}_{k})$. Here, the number of
classes $K_{t}$ increases if when sampling from the CRP predictive
distribution the result is $K_{t-1}+1$. That is, at the beginning of each
iteration we create a new class with an emission probability given by (8),
which is only kept if the MAP estimate is $z_{t}^{\star}=K_{t-1}+1$.

Mixture models do not usually have closed-form solutions for the estimates of
the parameters and the class assignments. Therefore, it is necessary to resort
to the expectation-maximization (EM) algorithm [13], for which we need the
log-likelihood of the complete data, which is given by

$\mathcal{L}_{\bm{\varphi}}(\bm{x}_{1:t},\bm{z}_{1:t})=\log
p(\bm{x}_{1:t},\bm{z}_{1:t}|\\{\bm{\varphi}_{k}\\}^{K_{t}}_{k=1})=\\\ \log
p(\bm{z}_{1:t})+\sum_{\tau=1}^{t}\log
p(x_{\tau}|z_{\tau},\\{\bm{\varphi}_{k}\\}^{K_{t}}_{k=1}),$ (10)

where the prior distribution $p(\bm{z}_{1:t})$ factorizes as

$p(\bm{z}_{1:t})=p(z_{t}|\bm{z}_{1:t-1})p(z_{t-1}|\bm{z}_{1:t-2})\cdots
p(z_{1}).$

This factorization is possible due to the chain-rule and the CRP construction
described in Section 3.1. Once the complete data log-likelihood is available,
we may apply the expectation step (E-step) and the maximization step (M-step)
of the EM algorithm. In this work, we have slightly modified the M-step to
accept the proposed continual learning framework. concretely, the estimation
of the parameter at each step is simply performed by taking one iterate of a
steepest descent method, yielding a stochastic M-step [14]. The E-step amounts
to

$\mathbb{E}[\mathbb{I}\\{z_{t}=k\\}|\bm{z}_{1:t-1}^{\star},x_{t},\hat{\bm{\varphi}}_{k}^{(t-1)}]=p(z_{t}=k|\bm{z}_{1:t-1}^{\star},x_{t},\hat{\bm{\varphi}}_{k}^{(t-1)})\\\
\propto
p(x_{t}|z_{t}=k,\hat{\bm{\varphi}}_{k}^{(t-1)})p(z_{t}=k|\bm{z}_{1:t-1}^{\star}),$

where $\mathbb{E}[\cdot]$ is the expectation operator,
$\hat{\bm{\varphi}}^{(t)}_{k}$ is the estimate of $\bm{\varphi}_{k}$ at time
$t$, and we have exploited (8). In the M-step, the estimate of the parameters
$\\{\bm{\varphi}_{k}\\}^{K_{t}}_{k=1}$ is updated based on the gradient:

$\hat{\bm{\varphi}}^{(t)}_{k}\leftarrow\hat{\bm{\varphi}}^{(t-1)}_{k}+\eta_{k,t}\nabla_{\bm{\varphi}_{k}}\mathbb{E}[\mathcal{L}_{\bm{\varphi}}(\bm{x}_{1:t},\bm{z}_{1:t})],$
(11)

where $\eta_{k,t}$ is the (adaptive) learning rate for the $k$th class at time
$t$. In this expression, we have assumed that the same initial learning rate
is chosen for the parameters of a given class, but it is possible to select
multiple learning rates per class. Once we have the E- and M-steps, we can
compute the posterior of $z_{t}$ and maximize it to obtain $z_{t}^{\star}$ as
in (4). Finally, Algorithm 1 presents all the necessary computations of the
proposed recursive method at each time instant $t$ and the Python
implementation can be found in https://github.com/pmorenoz/continual_ihcpd for
reproducibility purposes.

##  5 Experiments

Fig. 2: Upper row plots show the well-drilling univariate signal for the
unbounded latent variable model (left) and the hierarchical CPD method (right)
with fixed $K$. The colors represent latent-class asignments. Bottom row plots
show the MAP estimates of the run-length.

In this section we evaluate the performance of the proposed method. We apply
the infinite-dimensional hierarchical BOCPD algorithm to real-world data, and
in particular, to a sequence of raw nuclear magnetic response measurements
taken during a well-drilling process. This data consists of $4500$ real-valued
univariate observations taken at a fixed sampling frequency. In the following,
we assume that the time steps are ordered and discrete for simplicity.

To apply the proposed model, we choose $p(x_{t}|z_{t}=k,\bm{\varphi}_{k})$ to
be Gaussian distributed with unknown mean and variance, that is,
$\bm{\varphi}_{k}=\\{\mu_{k},\sigma_{k}^{2}\\}$. Moreover, the model has two
hyperparameters that we need to select. The first one, which is related to the
CPD method, is the parameter $\lambda$ of the hazard function that is used as
the conditional prior, $p(r_{t}|r_{t-1})$. In the experiments, we have
selected $\lambda=10^{6}$. The second one is the parameter $\alpha$, which is
involved in the CRP construction, and controls how likely is the appearance of
a new unseen class. We set it to $\alpha=1.0$. For the stochastic M-step, we
use two different adaptive learning rates for the mean and variance whose
initial values are given by $\eta_{\mu}=1.0$ and $\eta_{\sigma}=0.02$.
Importantly, we made both learning rates decrease at a rate of $2\%$ per time-
step if $z_{t}=k$ was selected as the most likely latent class. This choice
avoids adapting very old parameters with new incoming data.

Figure 2 shows the results obtained for $t=4500$ iterations.111A video
demonstrating the complete simulation of the algorithms is available at
https://www.youtube.com/watch?v=ymZPNURhtIc. The unbounded model is compared
with the hierarchical CPD approach with an upper bound on the number of
classes $K=10$. In the upper figures we can see the well-drilling signals, as
well as the latent-class assignments in different colors for both approaches.
As can be observed, the final number of classes inferred by the CRP was
$K_{4500}=7$. In the bottom figures we show the MAP estimates of the run-
length, $r_{t}^{\star}$. These figures show that the MAP estimation of the
run-length aligns quite well with the signal transitions. Furthermore, it
should be noted that the proposed method is more robust to outliers as can be
seen for $t\approx 200$ and $t\approx 600$, where the outlier is captured by
the latent class assignment but a CP is not declared. In fact, the MAP
estimate of the run-length is noisier for the method with a fixed number of
classes than for the unbounded model.

Finally, it is important to note that, since the unbounded model creates new
classes as they become necessary, its computational complexity is smaller than
that of hierarchical CPD approach, which needs to estimate the parameters of
$K=10$ classes at every time step.

##  6 Discussion and Future Work

This work has extended the Bayesian online change-point detection method to
more complex scenarios by considering a hierarchical model, which is based on
latent-class variables. To prevent the limitation of fixing the order of the
hierarchical model _a priori_ , we allow for an unbounded number of classes
using the chinese restaurant process. Moreover, the inference of the class
assignments is done with an expectation-maximization algorithm, where the
M-step is carried out stochastically, that is, only one iteration of a
steepest descent method is taken. Finally, the performance of the proposed
method is validated empirically over real-world data. We show its robustness
and utility for the aforementioned purposes. In future work, it would be
interesting to extend it to multi-channel settings with multivariate
generative models.

## References

  * [1] E. Andersson, D. Bock, and M. Frisén, “Some statistical aspects of methods for detection of turning points in business cycles,” _Journal of Applied Statistics_ , vol. 33, no. 3, pp. 257–278, 2006.
  * [2] I. Berkes, E. Gombay, L. Horváth, and P. Kokoszka, “Sequential change-point detection in GARCH (p, q) models,” _Econometric Theory_ , vol. 20, no. 6, pp. 1140–1167, 2004.
  * [3] M. Raginsky, R. M. Willett, C. Horn, J. Silva, and R. F. Marcia, “Sequential anomaly detection in the presence of noise and limited feedback,” _IEEE Transactions on Information Theory_ , vol. 58, no. 8, pp. 5544–5562, 2012.
  * [4] V. Krishnamurthy, “Quickest detection POMDPs with social learning: Interaction of local and global decision makers,” _IEEE Transactions on Information Theory_ , vol. 58, no. 8, pp. 5563–5587, 2012.
  * [5] M. Arts, A. Bollig, and R. Mathar, “Quickest eigenvalue-based spectrum sensing using random matrix theory,” _arXiv preprint arXiv:1504.01628v1_ , 2015.
  * [6] L. Du, C.-H. Liu, M. Laghate, and D. Cabric, “Sequential detection of number of primary users in cognitive radio networks,” in _Asilomar Conf. Signals, Systems and Computers_ , 2015.
  * [7] R. P. Adams and D. J. C. MacKay, “Bayesian online changepoint detection,” _arXiv preprint arXiv:0710.3742_ , 2007.
  * [8] P. Moreno-Muñoz, D. Ramírez, and A. Artés-Rodríguez, “Change-point detection on hierarchical circadian models,” _arXiv preprint arXiv:1809.04197_ , 2018.
  * [9] M. B. Ring, “Continual learning in reinforcement environments,” Ph.D. dissertation, University of Texas at Austin, 1994.
  * [10] J. Schmidhuber, “Powerplay: Training an increasingly general problem solver by continually searching for the simplest still unsolvable problem,” _Frontiers in Psychology_ , vol. 4, p. 313, 2013.
  * [11] C. V. Nguyen, Y. Li, T. D. Bui, and R. E. Turner, “Variational continual learning,” _International Conference on Learning Representations (ICLR)_ , 2018.
  * [12] J. Pitman, “Combinatorial stochastic processes,” Dept. Statistics, UC Berkeley., Tech. Rep., 2002.
  * [13] A. P. Dempster, N. M. Laird, and D. B. Rubin, “Maximum likelihood from incomplete data via the EM algorithm,” _Journal of the Royal Statistical Society. Series B (Methodological)_ , vol. 39, no. 1, pp. 1–38, 1977\.
  * [14] O. Cappé and E. Moulines, “On-line expectation–maximization algorithm for latent data models,” _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , vol. 71, no. 3, pp. 593–613, 2009.
  * [15] J. F. C. Kingman, “The coalescent,” _Stochastic Processes and their Applications_ , vol. 13, no. 3, pp. 235–248, 1982.
  * [16] P. Orbanz and Y. W. Teh, “Bayesian nonparametric models,” _Encyclopedia of Machine Learning_ , pp. 81–89, 2010.