# Detecting gravitational-wave bursts from black hole binaries in the Galactic
Center with LISA

Alan M. Knee Department of Physics & Astronomy, University of British
Columbia, Vancouver, BC V6T 1Z1, Canada  Jess McIver Department of Physics &
Astronomy, University of British Columbia, Vancouver, BC V6T 1Z1, Canada
Smadar Naoz Department of Physics & Astronomy, University of California, Los
Angeles, CA 90095, USA  Mani L. Bhaumik Institute for Theoretical Physics,
Department of Physics & Astronomy, UCLA, Los Angeles, CA 90095  Isobel M.
Romero-Shaw Department of Applied Mathematics and Theoretical Physics,
Cambridge CB3 0WA, United Kingdom  Kavli Institute for Cosmology Cambridge,
Madingley Road Cambridge CB3 0HA, United Kingdom  Bao-Minh Hoang Department of
Physics & Astronomy, University of California, Los Angeles, CA 90095, USA
Mani L. Bhaumik Institute for Theoretical Physics, Department of Physics &
Astronomy, UCLA, Los Angeles, CA 90095

###### Abstract

Stellar-mass black hole binaries (BHBs) in galactic nuclei are gravitationally
perturbed by the central supermassive black hole (SMBH) of the host galaxy,
potentially inducing strong eccentricity oscillations through the eccentric
Kozai-Lidov (EKL) mechanism. These highly eccentric binaries emit a train of
gravitational-wave (GW) bursts detectable by the Laser Interferometer Space
Antenna (LISA)—a planned space-based GW detector—with signal-to-noise ratios
(SNRs) up to ${\sim}100$ per burst. In this work, we study the GW signature of
BHBs orbiting our galaxy’s SMBH, Sgr A∗, which are consequently driven to very
high eccentricities. We demonstrate that an unmodeled approach using a wavelet
decomposition of the data effectively yields the time-frequency properties of
each burst, provided that the GW frequency peaks between
$10^{-3}\,\,\mathrm{Hz}$–$10^{-1}\,\,\mathrm{Hz}$. The wavelet parameters may
be used to infer the eccentricity of the binary, measuring $\log_{10}(1-e)$
within an error of $20\%$. Our proposed search method can thus constrain the
parameter space to be sampled by complementary Bayesian inference methods,
which use waveform templates or orthogonal wavelets to reconstruct and
subtract the signal from LISA data.

††software: GWpy (Macleod et al., 2021), SciPy (Virtanen et al., 2020), NumPy
(Harris et al., 2020), Matplotlib (Hunter, 2007), LISA GW Response (Bayle et
al., 2022), PyTDI (Staab et al., 2022).

##  1 Introduction

The Laser Interferometer Space Antenna (LISA) is a planned space-borne
gravitational-wave (GW) observatory (Amaro-Seoane et al., 2017; Colpi et al.,
2024), currently set to launch in the late 2030s. With sensitivity in the
$10^{-4}\,\,\mathrm{Hz}$ to $10^{-1}\,\,\mathrm{Hz}$ frequency range, LISA
will unlock the source-rich mHz band of the GW spectrum, potentially detecting
GWs from stellar-mass black hole binary (BHB) inspirals, massive BHB mergers,
extreme mass-ratio inspirals, and millions of ultra-compact galactic binaries
composed of white dwarfs and neutron stars (Amaro-Seoane et al., 2023). The
ability of LISA to observe compact binaries long before they merge will
present unique opportunities to study the dynamics of hierarchical multi-body
systems (Amaro-Seoane et al., 2023), where the influence of the perturbing
body can have a measurable effect on the gravitational waveform, thus
revealing key insights into compact binary formation channels in dense stellar
environments (Sigurdsson & Hernquist, 1993; Portegies Zwart & McMillan, 2000;
McMillan & Portegies Zwart, 2003; Rodriguez et al., 2016; Belczynski &
Banerjee, 2020; Gerosa & Fishbach, 2021; Bartos et al., 2017; Tagawa et al.,
2020).

Galactic nuclei are expected to contain an abundance of stellar-mass black
holes (BHs), many of which may exist in binaries. Two-body relaxation
naturally causes heavier masses to migrate inwards, resulting in a dense
compact object core where BHBs are readily assembled and hardened through
repeated BH-BH encounters (Freitag et al., 2006; O’Leary et al., 2009, 2016;
Antonini & Rasio, 2016; Rodriguez et al., 2018a, b; Samsing, 2018). The Milky
Way is no exception, potentially hosting ${\sim}2\times 10^{4}$ BHs within the
inner parsec of the Galactic Center (Morris, 1993; Miralda-Escudé & Gould,
2000; Freitag et al., 2006; O’Leary et al., 2009; Naoz et al., 2018; Rose et
al., 2022). If the galaxy hosts a central supermassive black hole (SMBH), as
is believed to be true for most large galaxies (Kormendy & Ho, 2013) including
our own (Ghez et al., 2005; Gillessen et al., 2009), stellar-mass BHBs may
become bound to the SMBH on a relatively wide outer orbit, forming a
hierarchical triple (Antonini & Perets, 2012; Hoang et al., 2018; Rose et al.,
2020). We illustrate this special orbital configuration in Fig. 1.
Gravitational torques exerted by the SMBH can subsequently induce secular
eccentricity and inclination oscillations through the eccentric Kozai-Lidov
(EKL) mechanism, driving highly inclined BHBs towards eccentricities of nearly
unity (Kozai, 1962; Lidov, 1962; Naoz, 2016) and resulting in a GW source
relevant for LISA. These novel sources possess a unique GW signature, as their
extreme eccentricity concentrates the GWs into a train of bursts occurring
once every periastron passage (O’Leary et al., 2009; Gould, 2011; Kocsis &
Levin, 2012; Xuan et al., 2023b). Changes in the morphology and timing of the
bursts reflect the secular evolution of the BHB due to the SMBH, allowing one
to infer the orbital parameters of both the binary and its perturber (Romero-
Shaw et al., 2023). Eccentricity notably enhances GW emission and accelerates
the inspiral timescale (Peters & Mathews, 1963), possibly contributing to the
compact binary merger rate observed by ground-based GW detectors (Hoang et
al., 2018; Stephan et al., 2019; Lim & Rodriguez, 2020; Martinez et al., 2020;
Abbott et al., 2023).

Figure 1: Schematic diagram of a hierarchical triple system. The inner binary
(blue ellipse), consisting of two masses $m_{1}$ and $m_{2}$, orbits around a
tertiary mass $M$ on a wide outer orbit (grey ellipse). The angle $\iota$ is
the relative inclination between the inner and outer orbital planes. In this
work, we focus on triples where the inner binary is a stellar-mass BHB and the
tertiary body is a SMBH.

Previous works have established that SMBH-induced BHB eccentricity
oscillations are detectable by LISA out to a few Mpc (e.g., Randall & Xianyu,
2018; Hoang et al., 2019; Wang et al., 2021)111This effect can also take place
in a wide range of masses and separations, see e.g., Randall & Xianyu (2019);
Emami & Loeb (2020); Deme et al. (2020)., provided the oscillation timescale
is comparable to the mission lifetime. The loudest BHBs are expected to be
orbiting the $4\times 10^{6}\,\,M_{\odot}$ SMBH at the Galactic Center (Sgr
A∗), with an estimated few tens of detectable sources in the LISA band (Wang
et al., 2021; Xuan et al., 2023b). Due to their close proximity and heavy
masses, these binaries can accumulate signal-to-noise ratios (SNRs) between
$10^{2}$–$10^{4}$ over the course of a $4\,\,\mathrm{yr}$ LISA mission (Hoang
et al., 2019; Wang et al., 2021), necessitating efforts to model and subtract
them from LISA data (Littenberg et al., 2020; Littenberg & Cornish, 2023). In
this Letter, we study the GW signature of BHBs orbiting Sgr A∗, prioritizing
the initial task of detection. We calculate the GW signal in the time domain,
including the full instrument response of LISA,  by incorporating simulated
background noise to emulate a detection scenario as realistically as possible.
As a proof-of-concept, we explore the use of a wavelet decomposition, based on
the multi-resolution $Q$-transform (Chatterji et al., 2004; Chatterji, 2005),
to resolve bursts from eccentric BHBs in LISA data and constrain their time-
frequency properties. Additionally, we show that this time-frequency
information provides a means to infer the orbital period and eccentricity of
the BHB.

##  2 Simulated LISA data

We begin by examining the GW signal originating from a stellar-mass BHB
orbiting around a Sgr A∗-like SMBH, as observed by LISA. In order to generate
simulations of the expected signal, we calculate the time-domain GW
polarizations $h_{+,\times}$ for the BHB using a leading-order adiabatic
approximation described in Barack & Cutler (2004), which expresses the
waveform as a sum over harmonics of the orbital frequency:

$h_{+,\times}(t,\boldsymbol{\theta})=\frac{G^{2}}{c^{4}D}\frac{m_{1}m_{2}}{a}\sum_{n}h_{+,\times}^{(n)}(t,\boldsymbol{\theta})\,,$
(1)

where $m_{1}$ and $m_{2}$ are the component masses, $D$ is the distance from
the detector, and $\boldsymbol{\theta}$ denotes the standard (time-varying)
Keplerian orbital elements. The mathematical definitions of the waveform
harmonics $h_{+,\times}^{(n)}(t,\boldsymbol{\theta})$ are given in Barack &
Cutler (2004). We evolve the system in time by numerically integrating the
orbit-averaged equations of motion governing a hierarchical triple, including
gravitational quadrupole and octupole-level contributions (Naoz et al., 2011,
2013), general relativistic precession (Naoz et al., 2013), and gravitational
radiation (Peters, 1964). The explicit dynamical equations can be found in
Naoz (2016). We then evaluate Eq. 1 by sampling the time-evolved orbital
elements. We ensure stability by imposing that the BHB must not cross into the
Roche limit of the SMBH (Naoz & Silk, 2014; Hoang et al., 2018). Further, the
binaries we consider have hardened enough such that we expect them to resist
disruption by the local stellar environment (see e.g., Hoang et al. (2018);
Rose et al. (2020) for discussion on the survivability of wide binaries in
galactic nuclei).

Figure 2: Simulated LISA observations of a singular BHB source with component
masses $m_{1}=30\,M_{\odot}$ and $m_{2}=20\,M_{\odot}$, orbiting around a
$M_{\mathrm{SMBH}}=4\times 10^{6}\,M_{\odot}$ SMBH. The initial parameters of
the inner (outer) orbits are: semimajor axis $a=0.15\,\,\mathrm{au}$
($a_{\mathrm{out}}=250\,\,\mathrm{au}$), eccentricity $e=0.5$
($e_{\mathrm{out}}=0.9$), argument of periastron $\omega=0^{\circ}$
($\omega_{\mathrm{out}}=0^{\circ}$), and inclination $\iota=85^{\circ}$
between the inner and outer orbits. The source is placed at a distance of
$D=8\,\,\mathrm{kpc}$ with ecliptic latitude $\beta=-5.608^{\circ}$ and
longitude $\lambda=266.852^{\circ}$, corresponding to the Galactic Center. The
left panel shows the eccentricity evolution of the BHB over time. The inset
focuses on a $1\,\,\mathrm{yr}$ interval centered on the time of maximum
eccentricity ($e_{\mathrm{max}}\approx 0.9944$). The upper right panel shows
the time-domain strain measured by LISA (in terms of the TDI-$A$ channel)
during this same $1\,\,\mathrm{yr}$ period, where the signal of interest is
shown in blue, and the total strain is in grey. The noise background is
obtained from LDC-2b (Baghi, 2022). The lower right panel is a spectrogram of
the whitened strain, with the peak GW frequency (Eq. 2) represented by the
dashed red curve.

In LISA, the output data streams are the time-delay interferometry (TDI)
channels, e.g., the noise-orthogonal $\\{A,E,T\\}$ variables. These channels
are synthesized by combining time-shifted one-way phase measurements recorded
by each of LISA’s six laser links (Armstrong et al., 1999; Tinto & Dhurandhar,
2021), and are designed to suppress the otherwise dominant laser noise. We use
the instrument model implemented in LISA GW Response (Bayle et al., 2022) to
evaluate the interferometric phase shifts for each laser link given an
incident GW signal. We use PyTDI (Staab et al., 2022) to convert these
measurements into the various TDI channels. Finally, we inject the TDI signal
into mock LISA data from the Spritz edition of the LISA Data Challenges
(LDC-2b, Baghi, 2022). For simplicity, we use a “cleaned” version of the
$1\,\,\mathrm{yr}$ dataset, which is free of artifacts such as noise glitches
and data gaps but retains various instrumental and environmental background
noises.

Fig. 2 shows $1\,\,\mathrm{yr}$ of simulated LISA observations for a
representative $30\,M_{\odot}+20\,M_{\odot}$ BHB perturbed by a Sgr A∗-like
SMBH and undergoing eccentricity oscillations. The signal is characterized by
a train of repeating, short-duration pulses lasting a few minutes each, which
map onto the frequency domain as broadband bursts of excess power. The burst
cadence is equal to the orbital period of the BHB, which is about
$3\,\,\mathrm{days}$ for the example shown in Fig. 2. The majority of the GW
power is radiated near periastron, causing the spectrum of each to burst to
peak at a frequency of approximately (O’Leary et al., 2009)

$f_{\mathrm{p}}(a,e)\approx\frac{(1+e)^{1/2}}{(1-e)^{3/2}}f_{\mathrm{orb}}(a)\,,$
(2)

where $(a,e)$ are the BHB semimajor axis and eccentricity at the time the
burst was emitted, and
$f_{\mathrm{orb}}(a)=(2\pi)^{-1}\sqrt{G(m_{1}+m_{2})/a^{3}}$ is the orbital
frequency. The burst amplitude and frequency track the oscillating
eccentricity, as shown in the lower right panel of Fig. 2. The bursts also
experience complicated amplitude modulations attributed to both the source’s
evolution and the motion of LISA, the latter of which depends on the source
sky position and polarization. On account of the light travel time across
LISA’s heliocentric orbit (i.e. Roemer delay), the burst arrival times
measured by LISA deviate from the true orbital period by $\pm
500\,\,\mathrm{s}$, depending on the orbital phase of LISA.

Throughout this work, we assume the outer orbit of the triple is always
coplanar with the sky. This choice of orientation simplifies the data analysis
problem, since we can ignore phase shifts caused by the light travel time
across the outer orbit (Laeuger et al., 2024). Similarly, we do not consider
any corrections associated with the motion of the BHB (Robson et al., 2018;
Xuan et al., 2023a) or the gravitational potential of the SMBH (Kuntz & Leyde,
2023; Sberna et al., 2022), opting to focus on a simple scenario in this
proof-of-principle study and leave further generalization to future work.

##  3 Unmodeled burst detection

Motivated by their transient and burst-like nature, we adapt the $Q$-transform
(Brown, 1991) to detect highly eccentric BHBs in LISA data. Analysis pipelines
based on the $Q$-transform search for arbitrary GW transients by filtering the
data against windowed sinusoids (wavelets) in place of waveform templates
(Chatterji et al., 2004; Chatterji, 2005), and have been applied extensively
to process data from ground-based GW detectors and characterize transient
noise (Davis et al., 2021; Acernese et al., 2023). We leverage the GWpy
(Macleod et al., 2021) implementation of the $Q$-transform, which is itself
derived from earlier burst search pipelines (Chatterji et al., 2004;
Chatterji, 2005; Rollins, 2011; Robinet et al., 2020).

The $Q$-transform projects time series data onto overlapping time-frequency
planes covered by an array of rectangular tiles, with each tile representing a
wavelet. Each plane is parameterized by a quality factor, $Q\sim f/\Delta f$,
which fixes the time-frequency aspect ratio of the individual tiles for that
plane. The significance of any tile is quantified by its energy, given by
(Chatterji et al., 2004; Chatterji, 2005)

$|X(\tau,f,Q)|^{2}=\bigg{|}\int_{-\infty}^{\infty}x(t)w(t-\tau,f,Q)e^{-2\pi
ift}\,\mathrm{d}t\bigg{|}^{2}\,,$ (3)

where $(\tau,f,Q)$ are the time, frequency, and $Q$ of the tile, $x(t)$ is the
whitened strain data, and $w(t,f,Q)$ is a window function (specifically, a
bisquare window). We normalize the energy such that, in white noise,
$|X(\tau,f,Q)|^{2}$ is exponentially distributed with unit mean and variance.
The tile density is tuned to guarantee a fractional energy loss (mismatch)
between adjacent tiles no larger than a desired amount, with smaller
mismatches resulting in denser time-frequency tilings. The analysis outputs
all tiles with SNR exceeding a pre-determined threshold, called triggers,
where we define the trigger SNR as

$\hat{\rho}=\sqrt{|X(\tau,f,Q)|^{2}-1}\,.$ (4)

Because a single burst typically produces multiple triggers, we group
significant triggers that are overlapping or adjacent in time into clusters,
and report the maximum-SNR trigger for each cluster. In our implementation,
the $Q$-transform is carried out in two stages: the first stage performs a
“quick pass” of the data using a low-resolution $Q$-tiling at $15\%$ mismatch;
after clustering triggers with SNR greater than $8$, we search
$5\,\,\mathrm{hr}$ of data around each trigger with a high-resolution tiling
at $1\%$ mismatch to further refine the trigger parameters.

Prior to calculating the $Q$-transform, we whiten the data by its power
spectral density (PSD), which is estimated empirically via a Welch median
(Welch, 1967) with $2^{19}\,\,\mathrm{s}$ Hann-windowed and $50\%$ overlapping
segments. To minimize over-whitening, we calculate the PSD on a “gated”
version of the data (Zweizig & Riles, 2021; Usman et al., 2016). This gating
is done by first dividing the data into short segments of duration
$T_{\mathrm{gate}}=300\,\,\mathrm{s}$; if any of the strain values exceed some
threshold, the data are multiplied by an inverse Tukey window, which is zero
during the segment containing the excursion, and smoothly transitions from
zero to one over a duration $T_{\mathrm{gate}}/4$ on either side. Through
trial and error, we find that a gating threshold of $4.5$ times the root-mean-
square strain (averaged over the full time series) is sufficient to remove the
worst-offending bursts and prevent over-whitening.

Figure 3: Time-frequency analysis of LISA data using the $Q$-transform. The
data consists of simulated noise (LDC-2b, Baghi, 2022) with added signals from
three BHBs with masses $m_{1}=30\,M_{\odot}$ and $m_{2}=20\,M_{\odot}$
orbiting a $M_{\mathrm{SMBH}}=4\times 10^{6}\,M_{\odot}$ SMBH at the Galactic
Center. A fourth source located $50\,\,\mathrm{kpc}$ away is also shown. The
initial BHB parameters are given in the legend. The initial outer orbit is
specified by $a_{\mathrm{out}}=250\,\,\mathrm{au}$, $e_{\mathrm{out}}=0.9$,
$\iota=85^{\circ}$ for each source. The left panel is a visualization of the
raw $Q$-transform output ($15\%$ mismatch, SNR threshold $8$) around a burst,
where the time axis is given in minutes from the trigger time
($t_{\mathrm{trig}}\approx 54.02\,\,\mathrm{d}$). The SNR of each tile is
given by its color. The square in the centre is the maximum-SNR trigger for
this cluster. The right panel shows the time and SNR of all triggers ($1\%$
mismatch) with SNR above $8$, each one corresponding to a detected burst.

Fig. 3 shows all triggers with SNR above $8$ after analyzing simulated LISA
data with the $Q$-transform. For simplicity, we consider just a single TDI
channel222Our analysis could be improved by using additional TDI channels to
perform glitch rejection, as true GWs will propagate through LISA differently
from instrumental glitches (Robson & Cornish, 2019)., TDI-$A$. We search over
time-frequency planes with low $Q$-values333The lower bound of $Q=\sqrt{11}$
is an anti-aliasing condition., between $\sqrt{11}$ and $32$, specifically
targeting bursts with short durations (maximizing timing accuracy) and high
bandwidths. The searched frequency range is from $10^{-4}\,\,\mathrm{Hz}$ to
$10^{-1}\,\,\mathrm{Hz}$, covering the nominal LISA observing band. The data
contains three sources representing BHBs orbiting a SMBH at the Galactic
Center, amidst background LISA noise from LDC-2b (Baghi, 2022). We also
simulate a fourth source identical to the example from Fig. 2, but placed at a
distance of $50\,\,\mathrm{kpc}$ (roughly the distance to the Large Magellanic
Cloud) instead of $8\,\,\mathrm{kpc}$, showing how these more distant bursts
are still individually resolvable at sufficiently large eccentricities. Three
of the injected sources (red and green markers) achieve maximum eccentricity
at the midpoint of the simulated observation period, while another source
(blue markers) achieves it shortly after the end of the observation period.
Longer mission lifetimes not only increase the chances that LISA will see one
of these putative sources during a high-eccentricity state, but will also
improve LISA’s ability to study the source evolution over time and determine
whether the dynamics are consistent with EKL-induced oscillations.

An important limitation of the $Q$-transform is that it is only sensitive to
well-localized excess power, and thus cannot combine the power from multiple
bursts far apart in time. Even if the summed SNR from all bursts is above our
threshold, the source will not be detected if the SNR of individual bursts is
consistently below the detection threshold. Our approach could thus be
improved upon by stacking power from multiple bursts assuming some burst
timing model (Romero-Shaw et al., 2023; Arredondo & Loutrel, 2021), allowing
one to filter against a train of wavelets instead of isolated wavelets.

##  4 Parameter recovery

Even though the $Q$-transform only provides generic time-frequency information
about the bursts, it is possible to make rough inferences about some of the
source’s orbital parameters based on two key factors: the timing of the
bursts, and the characteristic “peak” frequency where most of the excess power
is concentrated. We outline here a simple method of estimating the orbital
period and eccentricity evolution of BHBs using only the trigger time-
frequency parameters.

###  4.1 Orbital period

Bursts from highly eccentric binaries occur once per orbit at periastron, thus
we can measure the orbital period by observing the time interval between
consecutive burst triggers. However, the $Q$-transform is also likely to
detect non-repeating events associated with, e.g., massive BHB mergers or
detector glitches, obfuscating the eccentric BHB signal. In addition, LISA may
detect bursts from more than one eccentric binary, introducing multiple
periodicities into the trigger catalog. Numerous period detection schemes have
been studied in the astronomical literature (Koen, 2016). These methods
typically involve the phase-folding of arrival time data over many trial
periods and evaluating a test statistic, e.g. the Rayleigh or $H$-test (de
Jager et al., 1989), at each trial. Here, we apply a period-finding statistic
from Nishiguchi & Kobayashi (2000) to identify repeating bursts from
individual binaries,

$\mathcal{D}(T)=\bigg{|}\sum_{k=2}^{N_{\mathrm{t}}}\sum_{l=1}^{k-1}\mathbbm{1}(|T-(t_{k}-t_{l})|<\epsilon
T)e^{2\pi it_{k}/T}\bigg{|}\,,$ (5)

where $T$ is a trial orbital period, $0<\epsilon<1$ is a small jitter
fraction, $t_{i}$ are the trigger times, $N_{\mathrm{t}}$ is the total number
of triggers, and $\mathbbm{1}(\cdot)$ is the indicator function, which has a
value of $1$ when its argument is true and $0$ otherwise. Eq. 5 essentially
iterates over all pairs of triggers $(t_{k},t_{l})$ and counts the number of
pairs separated by an interval between $T(1\pm\epsilon)$.

Figure 4: Recovered orbital period and eccentricity for three simulated BHBs
with masses $m_{1}=30\,M_{\odot}$ and $m_{2}=20\,M_{\odot}$ orbiting a
$M_{\mathrm{SMBH}}=4\times 10^{6}\,M_{\odot}$ SMBH at the Galactic Center. The
initial orbital parameters are repeated from Fig. 3. The left panel shows
periodograms of simulated LISA triggers over $0.5$ (light red) and
$1\,\,\mathrm{yr}$ (black) observing periods, calculated via Eq. 5 with period
step size $\Delta T=500\,\,\mathrm{s}$ and jitter $\epsilon=0.005$. The right
panel shows the reconstructed eccentricity evolution of each BHB inferred from
the $Q$-transform triggers. The solid lines are the eccentricity calculated
from the central frequencies of each trigger, and the shaded region represents
the uncertainty due to the bandwidth of each trigger. The true eccentricities
are shown by the dashed curves.

The $\mathcal{D}(T)$ statistic has a number of advantages that make it well-
suited for our problem: the jitter fraction $\epsilon$ can be adjusted to
admit a certain level of deviation from the expected period; the weighting by
$e^{2\pi it/T}$ suppresses higher harmonics of the fundamental period,
allowing for easier identification of distinct sources; lastly, it is robust
to one-off transients like glitches or mergers, assuming these occur
randomly444Results from LISA Pathfinder showed that the waiting time between
acceleration glitches was exponentially distributed (LISA Pathfinder
Collaboration, 2022). However, it is not inconceivable that a full-size LISA
could experience periodic glitching if the glitch source has some periodic
aspect to it. The efficacy of our proposed method will thus depend on how well
such glitches are mitigated.. The main downside to Eq. 5 is that it checks all
trigger pairs, and so the computational cost scales as $N_{\mathrm{t}}^{2}$.
Also, if there are many missing bursts, the value of $\mathcal{D}(T)$ at the
true period will be reduced.

Evaluating Eq. 5 on the BHB triggers from Fig. 3 over many trial periods
results gives a periodogram, with sharp peaks corresponding to the orbital
periods of each binary, as shown in the left panel of Fig. 4. To test the
robustness of this period search method in the presence of extraneous
triggers, we artificially add $500$ triggers distributed in time with uniform
probability between the start and end of the observing period. Even though the
number of non-periodic triggers is comparable to the number of periodic
triggers, each BHB period is clearly detected above the noise level. For
comparison, we repeat our analysis with triggers from the first
$0.5\,\,\mathrm{yr}$ of data, showing how the peak prominence rises linearly
with the number of detected bursts. Note that we fix $\epsilon=0.005$, which
limits the period deviation to $0.5\%$ of the folded trial period. Since we do
not simulate the orbit of the inner binary around the SMBH, the only source of
period deviation is due to the travel time across the orbit of LISA.
Increasing $\epsilon$ allows Eq. 5 to capture larger variations in the period,
but at the expense of a higher noise level. Larger simulations are needed to
flesh out the noise distribution of $\mathcal{D}(T)$ and determine appropriate
thresholds for claiming evidence for periodicity.

###  4.2 Eccentricity evolution

We reconstruct the BHB eccentricity by leveraging the relation between its
orbital frequency and peak GW frequency, as given in Eq. 2. For the orbital
frequency $f_{\mathrm{orb}}$ we substitute the value measured via Eq. 5, and
for the peak GW frequency $f_{\mathrm{p}}$ we substitute the central frequency
of the maximum-SNR trigger for each burst. A simple inversion of Eq. 2 yields
an estimate of the eccentricity at the time each burst was emitted, which lets
us to examine the EKL-induced eccentricity evolution over time.

We show reconstructed eccentricity tracks of simulated BHBs orbiting a SMBH at
the Galactic Center in the right panel of Fig. 4. The eccentricity estimator
roughly tracks the true eccentricity near the eccentricity peak, but is prone
to overestimating at lower eccentricities, particularly at peak GW frequencies
near $f_{\mathrm{p}}\sim 10^{-3}\,\,\mathrm{Hz}$. Further inspection shows
that the peak frequency measured by LISA deviates from the true peak frequency
at lower frequencies. This is because the bursts are effectively filtered by
the TDI response functions of LISA (e.g., Cornish & Rubbo, 2003; Flauger et
al., 2021), which have a strong frequency dependence that shifts the detector-
frame peak frequency to higher frequencies. The whitening filter follows a
similar shape to the response function and compensates somewhat for this
effect, but is not sufficient to shift the measured peak frequency back to the
source-frame value for all frequencies. The effect of the LISA/TDI response
functions on transient broadband signals should be studied in future work.

Figure 5: Eccentricity errors for simulated BHBs with masses
$m_{1}=30\,M_{\odot}$ and $m_{2}=20\,M_{\odot}$ orbiting a
$M_{\mathrm{SMBH}}=4\times 10^{6}\,M_{\odot}$ at the Galactic Center. As
before, the initial outer orbital parameters are
$a_{\mathrm{out}}=250\,\,\mathrm{au}$, $e_{\mathrm{out}}=0.9$, and
$\iota=85^{\circ}$ in all simulations. Each source is plotted at its maximum
eccentricity, and the relative error (marker colours) in $\log_{10}(1-e)$ is
averaged over a month of data centered on this time. For simulations with
fewer than two detected bursts (black cross), we cannot measure their
eccentricity. The dashed grey contours are lines of constant peak GW frequency
assuming a $30\,M_{\odot}+20\,M_{\odot}$ BHB. The maximum BHB semimajor axis
allowed for this particular outer orbit (the Roche limit) is shown by the
vertical dash-dot line.

To better understand the accuracy of our eccentricity estimator over the
parameter space, we generate $135$ galactic BHB simulations on a grid in
initial semimajor axis and eccentricity, and recover the eccentricity of each
source using our trigger-based method. Fig. 5 shows the relative eccentricity
errors in terms of

$\delta
e=1-\bigg{|}\frac{\log_{10}(1-e_{\mathrm{meas}}(t))}{\log_{10}(1-e_{\mathrm{true}}(t))}\bigg{|},$
(6)

where $e_{\mathrm{meas}}(t)$ and $e_{\mathrm{true}}(t)$ are respectively the
measured and true eccentricity at trigger time $t$. Eq. 6 is averaged over a
month of data centered on the time of maximum eccentricity for each
simulation. Overall, the relative error is less than $20\%$ across all
simulations, and the error is lowest for peak frequencies near
$f_{\mathrm{p}}\sim 10^{-2}\,\,\mathrm{Hz}$. At lower frequencies, the
eccentricity becomes overestimated for the reasons discussed above. Moreover,
the eccentricity is underestimated as the peak frequency approaches the
Nyquist limit of $0.1\,\,\mathrm{Hz}$, where our frequency resolution quickly
deteriorates on account of the logarithmic frequency scaling of the
$Q$-transform. This low sample rate hinders our ability to characterize bursts
from the most eccentric BHBs in our suite of simulations.

##  5 Summary

A detection of binary eccentricity oscillations would offer conclusive
evidence that the binary existed in a triple, providing key insight into the
formation and evolution of stellar-mass compact binaries in dense cluster
environments. In this work, we have studied the GW signature of stellar-mass
BHBs in the Galactic Center that are driven to high eccentricities by a
perturbing SMBH, and outlined techniques for detecting and identifying these
sources in LISA data. BHBs in our own galaxy offer the best chances to detect
EKL-driven dynamics, though we expect such systems to exist in other galaxies
hosting SMBHs as well. We showed that an unmodeled time-frequency analysis can
detect GW bursts emitted by highly eccentric BHBs in our galaxy, and
characterize their time-frequency properties well enough to provide a rough
estimate of the BHB orbital period and eccentricity evolution. This
information can be passed to downstream analyses which use wavelets or
waveform models to fit the signal, reducing the parameter space to be explored
by these more computationally expensive techniques. Although we have focused
on eccentric BHBs perturbed by a SMBH, our analysis techniques could be
straightforwardly applied to generically eccentric compact binaries which do
not undergo EKL oscillations.

We have not considered the effects of the BHB’s orbital motion around the
SMBH, which will introduce long-period Doppler modulations of the GW signal
(Laeuger et al., 2024). Such modulations will change the arrival time of the
bursts relative to the orbital period, allowing one to better constrain the
outer orbit of the triple, but will also complicate the periodicity search
described in Sec. 4.1. We worked with $1\,\,\mathrm{yr}$ of simulated LISA
data for this study, but emphasize that longer observing times provide better
opportunities to detect and study the long-term behaviour of EKL triples.

As this is a proof-of-concept exploration, there is much room to build upon
and optimize various aspects of our methodology. Our burst search pipeline is
based on the $Q$-transform, which uses bisquare-windowed sinusoids to maximize
detection efficiency. However, this wavelet basis is overcomplete and cannot
be used for direct signal reconstruction and subtraction, which will likely be
a necessary step for the LISA global fit (Littenberg et al., 2020; Littenberg
& Cornish, 2023). Therefore, it will be worthwhile to explore alternative
bases more suitable for burst subtraction, such as orthogonal Meyer wavelets,
as done in the cWB pipeline (Klimenko & Mitselmakher, 2004; Klimenko et al.,
2008; Drago et al., 2020), or Morlet-Gabor wavelets with a manually enforced
orthogonality condition, as done in BayesWave (Cornish & Littenberg, 2015;
Cornish et al., 2021). Another limitation of our method is that we search for
bursts in isolation and do not stack power coherently from multiple bursts.
Our sensitivity could potentially be enhanced by filtering on templates which
chain together many bursts, where the timing between bursts and their relative
amplitudes would be parameters that must be searched over.

This research was enabled in part by support provided by the Digital Research
Alliance of Canada (alliance.can.ca). A.M.K. and J.M. acknowledge funding
support from the Natural Sciences and Engineering Research Council of Canada.
A.M.K. acknowledges funding support from a Killam Doctoral Scholarship. J.M.
acknowledges support from the Canada Research Chairs program. S.N.
acknowledges the partial support from NASA ATP 80NSSC20K0505 and from NSF-AST
2206428 grant and thanks Howard and Astrid Preston for their generous support.
I.M.R.-S. acknowledges support received from the Herchel Smith Postdoctoral
Fellowship Fund.

## References

  * Abbott et al. (2023) Abbott, R., Abbott, T. D., Acernese, F., et al. 2023, Physical Review X, 13, 011048, doi: 10.1103/PhysRevX.13.011048
  * Acernese et al. (2023) Acernese, F., et al. 2023, Class. Quant. Grav., 40, 185006, doi: 10.1088/1361-6382/acd92d
  * Amaro-Seoane et al. (2017) Amaro-Seoane, P., Audley, H., Babak, S., et al. 2017, arXiv e-prints, arXiv:1702.00786, doi: 10.48550/arXiv.1702.00786
  * Amaro-Seoane et al. (2023) Amaro-Seoane, P., Andrews, J., Arca Sedda, M., et al. 2023, Living Reviews in Relativity, 26, 2, doi: 10.1007/s41114-022-00041-y
  * Antonini & Perets (2012) Antonini, F., & Perets, H. B. 2012, Astrophys. J., 757, 27, doi: 10.1088/0004-637X/757/1/27
  * Antonini & Rasio (2016) Antonini, F., & Rasio, F. A. 2016, ApJ, 831, 187, doi: 10.3847/0004-637X/831/2/187
  * Armstrong et al. (1999) Armstrong, J. W., Estabrook, F. B., & Tinto, M. 1999, The Astrophysical Journal, 527, 814, doi: 10.1086/308110
  * Arredondo & Loutrel (2021) Arredondo, J. N., & Loutrel, N. 2021, Class. Quant. Grav., 38, 165001, doi: 10.1088/1361-6382/ac1083
  * Baghi (2022) Baghi, Q. 2022, arXiv e-prints, arXiv:2204.12142, doi: 10.48550/arXiv.2204.12142
  * Barack & Cutler (2004) Barack, L., & Cutler, C. 2004, Phys. Rev. D, 69, 082005, doi: 10.1103/PhysRevD.69.082005
  * Bartos et al. (2017) Bartos, I., Kocsis, B., Haiman, Z., & Márka, S. 2017, Astrophys. J., 835, 165, doi: 10.3847/1538-4357/835/2/165
  * Bayle et al. (2022) Bayle, J.-B., Baghi, Q., Renzini, A., & Le Jeune, M. 2022, LISA GW Response, 1.1, Zenodo, doi: 10.5281/zenodo.6423436
  * Belczynski & Banerjee (2020) Belczynski, K., & Banerjee, S. 2020, A&A, 640, L20, doi: 10.1051/0004-6361/202038427
  * Brown (1991) Brown, J. C. 1991, The Journal of the Acoustical Society of America, 89, 425, doi: 10.1121/1.400476
  * Chatterji (2005) Chatterji, S. 2005, Phd thesis, Massachusetts Institute of Technology
  * Chatterji et al. (2004) Chatterji, S., Blackburn, L., Martin, G., & Katsavounidis, E. 2004, Classical and Quantum Gravity, 21, S1809, doi: 10.1088/0264-9381/21/20/024
  * Colpi et al. (2024) Colpi, M., et al. 2024.  https://arxiv.org/abs/2402.07571
  * Cornish & Littenberg (2015) Cornish, N. J., & Littenberg, T. B. 2015, Classical and Quantum Gravity, 32, 135012, doi: 10.1088/0264-9381/32/13/135012
  * Cornish et al. (2021) Cornish, N. J., Littenberg, T. B., Bécsy, B., et al. 2021, Phys. Rev. D, 103, 044006, doi: 10.1103/PhysRevD.103.044006
  * Cornish & Rubbo (2003) Cornish, N. J., & Rubbo, L. J. 2003, Phys. Rev. D, 67, 022001, doi: 10.1103/PhysRevD.67.029905
  * Davis et al. (2021) Davis, D., et al. 2021, Class. Quant. Grav., 38, 135014, doi: 10.1088/1361-6382/abfd85
  * de Jager et al. (1989) de Jager, O. C., Raubenheimer, B. C., & Swanepoel, J. W. H. 1989, A&A, 221, 180
  * Deme et al. (2020) Deme, B., Hoang, B.-M., Naoz, S., & Kocsis, B. 2020, ApJ, 901, 125, doi: 10.3847/1538-4357/abafa3
  * Drago et al. (2020) Drago, M., et al. 2020.  https://arxiv.org/abs/2006.12604
  * Emami & Loeb (2020) Emami, R., & Loeb, A. 2020, MNRAS, 495, 536, doi: 10.1093/mnras/staa1200
  * Flauger et al. (2021) Flauger, R., Karnesis, N., Nardini, G., et al. 2021, JCAP, 01, 059, doi: 10.1088/1475-7516/2021/01/059
  * Freitag et al. (2006) Freitag, M., Amaro-Seoane, P., & Kalogera, V. 2006, ApJ, 649, 91, doi: 10.1086/506193
  * Gerosa & Fishbach (2021) Gerosa, D., & Fishbach, M. 2021, Nature Astronomy, 5, 749, doi: 10.1038/s41550-021-01398-w
  * Ghez et al. (2005) Ghez, A. M., Salim, S., Hornstein, S. D., et al. 2005, ApJ, 620, 744, doi: 10.1086/427175
  * Gillessen et al. (2009) Gillessen, S., Eisenhauer, F., Trippe, S., et al. 2009, ApJ, 692, 1075, doi: 10.1088/0004-637X/692/2/1075
  * Gould (2011) Gould, A. 2011, ApJ, 729, L23, doi: 10.1088/2041-8205/729/2/L23
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hoang et al. (2019) Hoang, B.-M., Naoz, S., Kocsis, B., Farr, W. M., & McIver, J. 2019, ApJ, 875, L31, doi: 10.3847/2041-8213/ab14f7
  * Hoang et al. (2018) Hoang, B.-M., Naoz, S., Kocsis, B., Rasio, F. A., & Dosopoulou, F. 2018, ApJ, 856, 140, doi: 10.3847/1538-4357/aaafce
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Klimenko & Mitselmakher (2004) Klimenko, S., & Mitselmakher, G. 2004, Class. Quant. Grav., 21, S1819, doi: 10.1088/0264-9381/21/20/025
  * Klimenko et al. (2008) Klimenko, S., Yakushin, I., Mercer, A., & Mitselmakher, G. 2008, Class. Quant. Grav., 25, 114029, doi: 10.1088/0264-9381/25/11/114029
  * Kocsis & Levin (2012) Kocsis, B., & Levin, J. 2012, Phys. Rev. D, 85, 123005, doi: 10.1103/PhysRevD.85.123005
  * Koen (2016) Koen, C. 2016, Monthly Notices of the Royal Astronomical Society, 459, 3012, doi: 10.1093/mnras/stw848
  * Kormendy & Ho (2013) Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511, doi: 10.1146/annurev-astro-082708-101811
  * Kozai (1962) Kozai, Y. 1962, AJ, 67, 591, doi: 10.1086/108790
  * Kuntz & Leyde (2023) Kuntz, A., & Leyde, K. 2023, Phys. Rev. D, 108, 024002, doi: 10.1103/PhysRevD.108.024002
  * Laeuger et al. (2024) Laeuger, A., Seymour, B., Chen, Y., & Yu, H. 2024, Phys. Rev. D, 109, 064086, doi: 10.1103/PhysRevD.109.064086
  * Lidov (1962) Lidov, M. L. 1962, Planet. Space Sci., 9, 719, doi: 10.1016/0032-0633(62)90129-0
  * Lim & Rodriguez (2020) Lim, H., & Rodriguez, C. L. 2020, Phys. Rev. D, 102, 064033, doi: 10.1103/PhysRevD.102.064033
  * LISA Pathfinder Collaboration (2022) LISA Pathfinder Collaboration. 2022, arXiv e-prints, arXiv:2205.11938, doi: 10.48550/arXiv.2205.11938
  * Littenberg & Cornish (2023) Littenberg, T. B., & Cornish, N. J. 2023, Phys. Rev. D, 107, 063004, doi: 10.1103/PhysRevD.107.063004
  * Littenberg et al. (2020) Littenberg, T. B., Cornish, N. J., Lackeos, K., & Robson, T. 2020, Phys. Rev. D, 101, 123021, doi: 10.1103/PhysRevD.101.123021
  * Macleod et al. (2021) Macleod, D. M., Areeda, J. S., Coughlin, S. B., Massinger, T. J., & Urban, A. L. 2021, SoftwareX, 13, 100657, doi: 10.1016/j.softx.2021.100657
  * Martinez et al. (2020) Martinez, M. A. S., Fragione, G., Kremer, K., et al. 2020, ApJ, 903, 67, doi: 10.3847/1538-4357/abba25
  * McMillan & Portegies Zwart (2003) McMillan, S., & Portegies Zwart, S. 2003, in Astronomical Society of the Pacific Conference Series, Vol. 296, New Horizons in Globular Cluster Astronomy, ed. G. Piotto, G. Meylan, S. G. Djorgovski, & M. Riello, 85
  * Miralda-Escudé & Gould (2000) Miralda-Escudé, J., & Gould, A. 2000, ApJ, 545, 847, doi: 10.1086/317837
  * Morris (1993) Morris, M. 1993, ApJ, 408, 496, doi: 10.1086/172607
  * Naoz (2016) Naoz, S. 2016, ARA&A, 54, 441, doi: 10.1146/annurev-astro-081915-023315
  * Naoz et al. (2011) Naoz, S., Farr, W. M., Lithwick, Y., Rasio, F. A., & Teyssandier, J. 2011, Nature, 473, 187, doi: 10.1038/nature10076
  * Naoz et al. (2013) —. 2013, MNRAS, 431, 2155, doi: 10.1093/mnras/stt302
  * Naoz et al. (2018) Naoz, S., Ghez, A. M., Hees, A., et al. 2018, ApJ, 853, L24, doi: 10.3847/2041-8213/aaa6bf
  * Naoz et al. (2013) Naoz, S., Kocsis, B., Loeb, A., & Yunes, N. 2013, Astrophys. J., 773, 187, doi: 10.1088/0004-637X/773/2/187
  * Naoz & Silk (2014) Naoz, S., & Silk, J. 2014, ApJ, 795, 102, doi: 10.1088/0004-637X/795/2/102
  * Nishiguchi & Kobayashi (2000) Nishiguchi, K., & Kobayashi, M. 2000, IEEE Transactions on Aerospace and Electronic Systems, 36, 407, doi: 10.1109/7.845217
  * O’Leary et al. (2009) O’Leary, R. M., Kocsis, B., & Loeb, A. 2009, MNRAS, 395, 2127, doi: 10.1111/j.1365-2966.2009.14653.x
  * O’Leary et al. (2016) O’Leary, R. M., Meiron, Y., & Kocsis, B. 2016, ApJ, 824, L12, doi: 10.3847/2041-8205/824/1/L12
  * Peters (1964) Peters, P. C. 1964, Phys. Rev., 136, B1224, doi: 10.1103/PhysRev.136.B1224
  * Peters & Mathews (1963) Peters, P. C., & Mathews, J. 1963, Phys. Rev., 131, 435, doi: 10.1103/PhysRev.131.435
  * Portegies Zwart & McMillan (2000) Portegies Zwart, S. F., & McMillan, S. L. W. 2000, ApJ, 528, L17, doi: 10.1086/312422
  * Randall & Xianyu (2018) Randall, L., & Xianyu, Z.-Z. 2018, ApJ, 864, 134, doi: 10.3847/1538-4357/aad7fe
  * Randall & Xianyu (2019) —. 2019, arXiv e-prints, arXiv:1907.02283, doi: 10.48550/arXiv.1907.02283
  * Robinet et al. (2020) Robinet, F., Arnaud, N., Leroy, N., et al. 2020, SoftwareX, 12, 100620, doi: 10.1016/j.softx.2020.100620
  * Robson & Cornish (2019) Robson, T., & Cornish, N. J. 2019, Phys. Rev. D, 99, 024019, doi: 10.1103/PhysRevD.99.024019
  * Robson et al. (2018) Robson, T., Cornish, N. J., Tamanini, N., & Toonen, S. 2018, Phys. Rev. D, 98, 064012, doi: 10.1103/PhysRevD.98.064012
  * Rodriguez et al. (2018a) Rodriguez, C. L., Amaro-Seoane, P., Chatterjee, S., et al. 2018a, Phys. Rev. D, 98, 123005, doi: 10.1103/PhysRevD.98.123005
  * Rodriguez et al. (2018b) Rodriguez, C. L., Amaro-Seoane, P., Chatterjee, S., & Rasio, F. A. 2018b, Phys. Rev. Lett., 120, 151101, doi: 10.1103/PhysRevLett.120.151101
  * Rodriguez et al. (2016) Rodriguez, C. L., Chatterjee, S., & Rasio, F. A. 2016, Phys. Rev. D, 93, 084029, doi: 10.1103/PhysRevD.93.084029
  * Rollins (2011) Rollins, J. G. 2011, Phd thesis, Columbia University
  * Romero-Shaw et al. (2023) Romero-Shaw, I., Loutrel, N., & Zevin, M. 2023, Phys. Rev. D, 107, 122001, doi: 10.1103/PhysRevD.107.122001
  * Rose et al. (2020) Rose, S. C., Naoz, S., Gautam, A. K., et al. 2020, ApJ, 904, 113, doi: 10.3847/1538-4357/abc557
  * Rose et al. (2022) Rose, S. C., Naoz, S., Sari, R., & Linial, I. 2022, ApJ, 929, L22, doi: 10.3847/2041-8213/ac6426
  * Samsing (2018) Samsing, J. 2018, Phys. Rev. D, 97, 103014, doi: 10.1103/PhysRevD.97.103014
  * Sberna et al. (2022) Sberna, L., et al. 2022, Phys. Rev. D, 106, 064056, doi: 10.1103/PhysRevD.106.064056
  * Sigurdsson & Hernquist (1993) Sigurdsson, S., & Hernquist, L. 1993, Nature, 364, 423, doi: 10.1038/364423a0
  * Staab et al. (2022) Staab, M., Bayle, J.-B., & Hartwig, O. 2022, PyTDI, 1.2, Zenodo, doi: 10.5281/zenodo.6351737
  * Stephan et al. (2019) Stephan, A. P., Naoz, S., Ghez, A. M., et al. 2019, ApJ, 878, 58, doi: 10.3847/1538-4357/ab1e4d
  * Tagawa et al. (2020) Tagawa, H., Haiman, Z., & Kocsis, B. 2020, Astrophys. J., 898, 25, doi: 10.3847/1538-4357/ab9b8c
  * Tinto & Dhurandhar (2021) Tinto, M., & Dhurandhar, S. V. 2021, Living Reviews in Relativity, 24, 1, doi: 10.1007/s41114-020-00029-6
  * Usman et al. (2016) Usman, S. A., et al. 2016, Class. Quant. Grav., 33, 215004, doi: 10.1088/0264-9381/33/21/215004
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Wang et al. (2021) Wang, H., Stephan, A. P., Naoz, S., Hoang, B.-M., & Breivik, K. 2021, ApJ, 917, 76, doi: 10.3847/1538-4357/ac088d
  * Welch (1967) Welch, P. 1967, IEEE Transactions on Audio and Electroacoustics, 15, 70, doi: 10.1109/TAU.1967.1161901
  * Xuan et al. (2023a) Xuan, Z., Naoz, S., & Chen, X. 2023a, Phys. Rev. D, 107, 043009, doi: 10.1103/PhysRevD.107.043009
  * Xuan et al. (2023b) Xuan, Z., Naoz, S., Kocsis, B., & Michaely, E. 2023b.  https://arxiv.org/abs/2310.00042
  * Zweizig & Riles (2021) Zweizig, J., & Riles, K. 2021, Information on self-gating of $h(t)$ used in O3 continuous-wave searches


# MultiEarth 2022 – Multimodal Learning for Earth and Environment
Workshop and Challenge

Miriam Cha1, Kuan Wei Huang2, Morgan Schmidt2, Gregory Angelides1, Mark
Hamilton2,
Sam Goldberg2, Armando Cabrera3, Phillip Isola2, Taylor Perron2, Bill
Freeman2, Yen-Chen Lin2,
Brandon Swenson3, Jean Piou1

1MIT Lincoln Laboratory, [miriam.cha, gregangelides<EMAIL_ADDRESS>
2MIT, [kwhuang, morgansc, markth, sgoldberg, phillipi, perron, billf,
<EMAIL_ADDRESS>
3DAF MIT AI Accelerator, [armando.cabrera<EMAIL_ADDRESS>
###### Abstract

The Multimodal Learning for Earth and Environment Challenge (MultiEarth 2022)
will be the first competition aimed at the monitoring and analysis of
deforestation in the Amazon rainforest at any time and in any weather
conditions. The goal of the Challenge is to provide a common benchmark for
multimodal information processing and to bring together the earth and
environmental science communities as well as multimodal representation
learning communities to compare the relative merits of the various multimodal
learning methods to deforestation estimation under well-defined and strictly
comparable conditions. MultiEarth 2022 will have three sub-challenges: 1)
matrix completion, 2) deforestation estimation, and 3) image-to-image
translation. This paper presents the challenge guidelines, datasets, and
evaluation metrics for the three sub-challenges. Our challenge website is
available at https://sites.google.com/view/rainforest-challenge.

##  1 Introduction

Despite international efforts to reduce deforestation, the world loses, so
far, an area of forest that is equivalent to the size of 40 football fields
every minute [11]. Deforestation in the Amazon rainforest accounts for the
largest share, contributing to reduced biodiversity, habitat loss, and climate
change. Since much of the region is difficult to access, satellite remote
sensing offers a powerful tool to track changes in the Amazon. However,
obtaining a continuous time series of images is hindered by seasonal weather,
clouds, smoke, and other inherent limitations of optical sensors. Synthetic
aperture radar (SAR), which is insensitive to lighting and weather conditions,
appears to be a well suited tool for the task, but SAR images are more
difficult for humans to interpret than optical images. A key component of this
challenge is to monitor the Amazon rainforest in all weather and lighting
conditions using our multimodal remote sensing dataset, which includes a time
series of multispectral and SAR images. The 2022 Multimodal Learning for Earth
and Environment Challenge (MultiEarth 2022) will be the first competition
aimed at monitoring the Amazon rainforest and predicting deforestation using
multimodal representation learning methods.

While considerable research has been devoted to tracking changes in forests
[4, 8, 6, 9], the analyses typically rely on passive, optical sensors such as
the U.S. Geological Survey’s Landsat, NASA’s Moderate Resolution Imaging
Spectroradiometer (MODIS), and the European Space Agency’s Sentinel-2 [1].
These passive, optical sensors need an unobstructed and illuminated view of
the scene in order to capture meaningful images. This limits their
effectiveness in monitoring the Amazon rainforest with year-round cloud
coverage. Beyond the Amazon, approximately 67 percent of Earth’s surface is
typically covered by clouds [7]. Therefore, using SAR for cloud-free
observation will have a broad applicability.

MultiEarth 2022 will conduct the following sub-challenges to support the
interpretation and analysis of the rainforest at any time and any weather
conditions:

  * •

Matrix Completion Sub-Challenge: given remote sensing images taken at
different locations and dates, and in different modalities, participants are
required to predict appearance at a novel [lon, lat, date, modality] query.
Performance is measured on the following visual metrics: Peak Signal-to-Noise
Ratio (PSNR), Structural Similarity Index Measure (SSIM) [10], Learned
Perceptual Image Patch Similarity (LPIPS) [12], and Fréchet Inception Distance
(FID) [5]. Detailed challenge description is provided in Section 3.1.

  * •

Deforestation Estimation Sub-Challenge: beyond visual appearance, participants
are required to classify whether a region is deforested or not at a novel
[lon, lat, date, modality] query. Modality will be ‘deforestation’ for this
sub-challenge. Participants will be given the multimodal remote sensing
dataset along with the deforestation label maps. Performance is measured on
the following metrics: pixel accuracy, F1 score, and Intersection over Union
(IoU). Detailed description of this sub-challenge can be found in Section 3.2.

  * •

Image-to-Image Translation Sub-Challenge: participants are required to predict
a set of possible cloud-free corresponding optical images given an input SAR
image. For this sub-challenge, we provide an aligned dataset (e.g.
$[\mathbf{x},[\mathbf{y}_{1},\mathbf{y}_{2},...,\mathbf{y}_{N}]]$) where an
input SAR image $\mathbf{x}$ is paired to a set of ground truth optical images
$[\mathbf{y}_{1},\mathbf{y}_{2},...,\mathbf{y}_{N}]$. Performance is evaluated
based on $\sum_{j}\min_{i}\Arrowvert
f(\mathbf{x})_{i}-\mathbf{y}_{j}\Arrowvert$ where $f(\cdot)$ is a prediction
of $\mathbf{y}$ given $\mathbf{x}$, and $f$ may make multiple predictions
indexed by $i$. Detailed guideline of the sub-challenge is provided in Section
3.3.

To be eligible to participate in the challenge, every entry has to be
accompanied by a paper presenting the results and the methods that created
them, which will undergo peer-review. The organizers reserve the right to re-
evaluate the findings, but will not participate in the Challenge.

Figure 1: Study area in the Amazon, bounded by (4.39∘ S, 55.2∘ W), (4.39∘ S, 54.48∘ W), (3.33∘ S, 54.48∘ W) and (3.33∘ S, 55.2∘ W). Sensor | Time | Bands | Resolution (m) | # Images | Link
---|---|---|---|---|---
Sentinel-1 | 2014-2021 | VV, VH | 10 | 859,627 |  | link1
---
link2
Sentinel-2 | 2018-2021 |  | B1, B2, B3, B4, B5,
---
B6, B7, B8, B8A, B9,
B11, B12, QA60
10 | 5,395,559 |  | link1
---
link2
Landsat 5 | 1984-2012 |  | SR_B1, SR_B2, SR_B3,
---
SR_B4, SR_B5, ST_B6,
SR_B7, QA_PIXEL
30 | 3,550,368 |  | link1
---
link2
Landsat 8 | 2013-2021 |  | SR_B1, SR_B2, SR_B3,
---
SR_B4, SR_B5, SR_B6,
SR_B7, ST_B10, QA_PIXEL
30 | 2,172,564 |  | link1
---
link2
Table 1: Overview of our multimodal remote sensing dataset that includes
Sentinel-1, Sentinel-2, Landsat 5, and Landsat 8.

##  2 Datasets

###  2.1 Multimodal Remote Sensing Dataset

Participants will receive a multimodal remote sensing dataset that consists of
Sentinel-1, Sentinel-2, Landsat 5, and Landsat 8 as reported in Table 1.
Sentinel-1 uses a synthetic aperture radar (SAR) instrument, which collects in
two polarization bands: VV (vertical transmit/vertical receive) and VH
(vertical transmit/horizontal receive). Sentinel-2, Landsat 5, and Landsat 8
use optical instruments, which measure in spectral bands in the visible and
infrared spectra. We also include in the dataset the associated layers with
cloud quality for Sentinel-2, Landsat 5, and Landsat 8 (i.e. QA60 and
QA_PIXEL). Detailed band designations for each sensor can be found in Google
Earth Engine Data Catalog111Sentinel-1: https://developers.google.com/earth-
engine/datasets/catalog/COPERNICUS_S1_GRD
Sentinel-2: https://developers.google.com/earth-
engine/datasets/catalog/COPERNICUS_S2_SR
Landsat 5: https://developers.google.com/earth-
engine/datasets/catalog/LANDSAT_LT05_C02_T1_L2
Landsat 8: https://developers.google.com/earth-
engine/datasets/catalog/LANDSAT_LC08_C02_T1_L2 .

The region of interest is the rectangle bounded by the points (4.39∘ S, 55.2∘
W), (4.39∘ S, 54.48∘ W), (3.33∘ S, 54.48∘ W) and (3.33∘ S, 55.2∘ W) as shown
in Figure 1. This area is further divided into squares centered on 54 $\times$
37 = 1998 coordinate pairs by iterating in 0.02∘ increments in the latitude
and longitude directions. Our dataset consists of 11,978,118 total training
images. (Note that the image-to-image translation sub-challenge will use a
subset of the total training images.) Test images are selected from the middle
7 columns (from longitudes -54.78 to -54.90), which are withheld as test
region. The training images come from all bands and dates from the 4 sensor
collections specified in Table 1 covering the area outside of the test region.
Each image file name is in the format {Collection}_{Band}_{Longitude}_
{Latitude}_{Year}_{Month}_{Day}.tiff. Sample images of Sentinel-1, Sentinel-2,
Landsat 5, and Landsat 8 from our dataset are shown in Figure 2.

Satellite images from this dataset are downloaded using Google’s Earth Engine
platform. Each coordinate pair is projected from UTM (Universal Transverse
Mercator) to EPSG (European Petroleum Survey Group). Then a 256$\times$256
image centered at the coordinate is extracted for Sentinel-1 and Sentinel-2,
and a 85$\times$85 image centered at the coordinate is extracted for Landsat 5
and Landsat 8. The difference in pixels accounts for the 10-meter resolution
Sentinel-1 and Sentinel-2 vs the 30-meter resolution Landsat 5 and Landsat 8.
If the Landsat images are upsampled to 256$\times$256, they will be
geospatially aligned with the Sentinel images.

The dataset can be downloaded from the links provided in Table 1. Data is
freely available for development, research, or educational purposes. For more
details, please refer to Google Earth Engine License
Agreement222https://earthengine.google.com/terms/.

(a) Sentinel-1

(b) Sentinel-2

(c) Landsat 5

(d) Landsat 8

Figure 2: Image examples from our multimodal dataset (LAT/LON: -4.11/-55.14)
(a) Sentinel-1 VV on 12/30/2021, (b) Sentinel-2 RGB (B4, B3, B2) on 8/4/2021,
(c) Landsat 5 RGB (SR_B3, SR_B2, SR_B1) on 6/12/1998, and (d) Landsat 8 RGB
(SR_B4, SR_B3, SR_B2) on 12/4/2021.

###  2.2 Deforestation Labels

Labeled data of deforestation for training and testing are manually labeled
using cloud-free monthly mosaic satellite images from Planet [2]. The Planet
imagery has 3.7 m spatial resolution and consists of 3 bands (RGB) with images
from a single month joined together to achieve an image with the least amount
of cloud cover. Most are cloud-free but some have a few clouds. Eleven time
slices (months) are labeled over six years, from 2016 to 2021. They include:
08/2016; 07, 08/2017; 06, 08/2018; 07, 08/2019; 06, 08/2020; and 05, 08/2021.
The August 2016 image contains a greater number of clouds and artifacts for
the mosaic process. August 2021 is labeled by the authors and the remaining
ten are labeled by the team at Scale AI [3]. The following guidelines are used
by Scale AI to label deforested areas within polygons: 1) areas deforested for
human use in any stage of regrowth; 2) forested areas 1 hectare or larger with
unbroken forest canopy are not labeled; 3) roads are included as deforested if
clearly visible; 4) rivers going through deforested areas are left unlabeled
if their area is greater than 1 hectare; 5) clouds or image artifacts are not
labeled; 6) snapping was used during labeling to connect polygons. For the
August 2021 labeled data, our team does not use snapping but instead leave
small spaces between polygons. Therefore, in the labeled dataset, unlabeled
areas are mostly forest, but may also contain some clouds or rivers.
Shapefiles with labeled polygons are converted into rasters with 0’s for
forested/other and 1’s for deforested areas. Rasters are divided into
georeferenced image chips (256$\times$256 pixels) that correspond with the
multimodal remote sensing dataset described in Section 2.1. We show sample
images of deforestation label maps at various time slices in Figure 3.

(a) 08/2016

(b) 07/2017

(c) 06/2018

(d) 08/2020

Figure 3: Example images of deforestation label maps corresponding to the
location shown in Figure 2 (LAT/LON: -4.11/-55.14) at various time slices:
08/2016, 07/2017, 06/2018, and 08/2020.

##  3 Challenge Tasks

###  3.1 Matrix Completion Sub-Challenge

Problem Definition. Most studies of forest monitoring rely on remote sensing
data collected by optical sensors [4, 8, 6, 9]. However, it is often
difficult, particularly in moist forest areas like the Amazon rainforest, to
obtain cloud-free images. This sub-challenge focuses on filling in spatial,
temporal, and modality gaps in remote sensing data, especially gaps created by
unfavorable lighting, weather conditions, or other atmospheric factors. The
Matrix Completion Sub-Challenge represents the sensory world as a huge tensor
of measurements, in multiple modalities, sampled at different places and
times. This tensor is typically very sparse, with most measurements missing,
since satellites only photograph a tiny fraction of all possible wavelengths
at all possible positions on the globe at all possible points in time.
Participants are required to solve the tensor completion problem on this huge
yet sparse multidimensional measurement tensor, and may develop and use
methods from representation learning and generative modeling, or may consider
other approaches.

Data. Training set from the multimodal remote sensing dataset described in
Section 2.1 is provided.

Metrics. To rank all submissions to this sub-challenge, imaging results
entered by the participants will be evaluated using the following four
metrics: Peak Signal-to-Noise Ratio (PNSR), Structural Similarity Index
Measure (SSIM) [10], Learned Perceptual Image Patch Similarity (LPIPS) [12],
and Fréchet Inception Distance (FID) [5].

Submission Format. For the test input, 2000 test queries will be provided as a
list of lists _i.e_.
[[${\text{lon}_{0},\text{lat}_{0},\text{date}_{0},\text{modality}_{0}}],\text{…},[\text{lon}_{1999},\text{lat}_{1999},\text{date}_{1999}$,
$\text{modality}_{1999}]]$. Each test query is in the format [lon, lat, date,
modality]. For example, $[-55.15,-4.11,2021\\_12\\_04,$
$\text{Landsat8\\_SR\\_B2}]$ will represent
Landsat8_SR_B2_-55.15_-4.11_2021_12_04.tiff. For the test output, participants
will submit in total 2,000 256$\times$256 images, one 256$\times$256 image for
each input test query. Imagery related to the input test queries will be made
available in our website and can be used to help generate the requested
output.

###  3.2 Deforestation Estimation Sub-Challenge

Problem Definition. Beyond predicting visual appearance, this sub-challenge is
aimed at estimating deforestation from the multimodal remote sensing dataset.
As described in the Matrix Completion Sub-Challenge, the multimodal dataset
can be represented as a multidimensional data matrix with missing entries. The
goal of this sub-challenge is to perform a binary classification to predict
whether a region is deforested or not. As solutions for the Matrix Completion
Sub-Challenge can be naturally extended, we strongly encourage participants in
the Matrix Completion Sub-Challenge to submit to the Deforestation Estimation
Sub-Challenge.

Data. Participants will use the multimodal remote sensing dataset described in
Section 2.1 to predict binary deforestation label maps in Section 2.2.

Metrics. Performance is measured on the following standard metrics: pixel
accuracy, F1 score, and Intersection over Union (IoU).

Submission Format. For the test input, 1000 test queries will be provided as a
list of lists _i.e_.
[[${\text{lon}_{0},\text{lat}_{0},\text{date}_{0},\text{modality}_{0}}],\text{…},[\text{lon}_{999},\text{lat}_{999},\text{date}_{999}$,
$\text{modality}_{999}]]$. Each test query is in the format [lon, lat, date,
modality]. For example, $[-55.15,-4.11,2021\\_08\\_01,$
$\text{deforestation}]$ will represent deforestation_-55.15_-4.11_2021_08.png.
To have a consistent naming convention for the date (_i.e_. year_month_day),
we add a nominal day label of “_01” to all deforestation estimation test
queries. For the test output, participants will submit in total 1,000
256$\times$256 binary masks, one 256$\times$256 binary mask for each input
test query. Imagery related to the input test queries will be made available
in our website and can be used to help generate the requested output.

###  3.3 Image-to-Image Translation Sub-Challenge

Problem Definition. Obtaining a continuous time series of view of the Amazon
rainforest is hindered by weather, clouds, smoke, and other inherent
limitations of passive sensors (e.g. optical sensors) that rely on sunlight.
Such limitations produce a major information gap in the Amazon rainforest that
gets rain throughout the entire year. Synthetic aperture radar (SAR) is an
active sensor that can collect images with relative invariance to weather and
lighting conditions. However, visual interpretation of SAR images is not
intuitive due to the large dynamic range, low spatial correlation, and radar-
specific geometry distortion. To enhance SAR interpretability, this sub-
challenge is aimed at modeling a distribution of possible electro-optical (EO)
image outputs conditioned on a SAR input image. Here, EO image is a 3-channel
RGB image from Sentinel-2. In Sentinel-2, RGB bands are represented as B4, B3,
and B2, respectively. SAR image is a 2-channel Sentinel-1 image consisting of
VV and VH bands. In this Image-to-Image Translation Sub-Challenge,
participants need to model a distribution of potential results in a
conditional generative modeling setting.

Data. For this sub-challenge, we provide JSON files specifying which
Sentinel-2 EO images (B4, B3, and B2) correspond to which SAR images. We
provide two JSON files: one for Sentinel-1 VV
band333https://rainforestchallenge.blob.core.windows.net/dataset/sentinel_vv_image_alignment_train.json
and another for Sentinel-1VH band
444https://rainforestchallenge.blob.core.windows.net/dataset/sentinel_vh_image_alignment_train.json.
The mappings in both files are identical. The aligned dataset will have the
following format:
$[\mathbf{x},[\mathbf{y}_{1},\mathbf{y}_{2},...,\mathbf{y}_{N}]]$ where a SAR
image $\mathbf{x}$ is paired to a set of ground truth EO images
$[\mathbf{y}_{1},\mathbf{y}_{2},...,\mathbf{y}_{N}]$. For each SAR image, all
EO images of the same geographic region and which were collected within 7 days
of the SAR image timestamp will be identified. On average 3 EO images are
paired with each SAR image ($N\approx 3$).

Metrics. Performance is evaluated based on following:

$\sum_{j}\min_{i}\Arrowvert f(\mathbf{x})_{i}-\mathbf{y}_{j}\Arrowvert$ (1)

where $f(\mathbf{x})_{i}$ is a set of possible EO images translated from a
generative model $f(\cdot)$ conditioned on an input SAR image $\mathbf{x}$,
and $\mathbf{y}_{j}$ is an EO image from the corresponding ground truth set.
This metric also evaluates diversity of generated output images staying
faithful to the diversity of the ground-truth EO data.

Submission Format. For testing, we will provide 5000 Sentinel-1 SAR images,
where each SAR image is 256$\times$256$\times$2 and the 2 channels correspond
to the Sentinel-1 VV and VH bands, respectively. For the test output,
participants will submit in total 15,000 256$\times$256$\times$3 EO images,
three translated 256$\times$256$\times$3 EO images for each input Sentinel-1
SAR image. The EO image channels correspond to the Sentinel-2 RGB channels
(B4, B3, B2). This is a multimodal image-to-image translation problem where
participants will generate three possible EO images given an input SAR image.

##  4 Conclusion

We introduce MultiEarth 2022 – the first open Multimodal Learning for Earth
and Environment Challenge. It comprises three sub-challenges: 1) matrix
completion, 2) deforestation estimation, and 3) image-to-image translation.
This manuscript describes MultiEarth 2022’s challenge conditions, data,
evaluation metrics, and submission guideline. Only a few labeled, multimodal
datasets including passive and active sensors have been publicly available,
limiting the number of participants who can research and analyze Earth’s
surface at all times and in all weather conditions. We collect and disseminate
a multimodal dataset that includes a continuous time series of Sentinel-1,
Sentinel-2, Landsat 5 and Landsat 8, with deforestation labels.

##  5 Acknowledgments

The authors would like to thank Dr. M.K. Newey for helping the team labeling
the August 2021 time slice data, and Scale AI, for its efforts to label, in a
short period of time, the remaining ten time slices considered in this
MultiEarth2022.

Research was sponsored by the United States Air Force Research Laboratory and
the United States Air Force Artificial Intelligence Accelerator and was
accomplished under Cooperative Agreement Number FA8750-19-2-1000. The views
and conclusions contained in this document are those of the authors and should
not be interpreted as representing the official policies, either expressed or
implied, of the United States Air Force or the U.S. Government. The U.S.
Government is authorized to reproduce and distribute reprints for Government
purposes notwithstanding any copyright notation herein.

## References

  * [1] European Space Agency, 2015. Sentinels: Space for Copernicus.  http://esamultimedia.esa.int.
  * [2] Planet Team.  https://api.planet.com.
  * [3] Scale AI.  https://scale.com.
  * [4] Maria Antonia Brovelli, Yaru Sun, and Vasil Yordanov.  Monitoring forest change in the amazon using multi-temporal remote sensing data and machine learning classification on google earth engine.  ISPRS International Journal of Geo-Information, 2020.
  * [5] Martin Heusel, Hubert Ramsauer, Thomas Unterthiner, Bernhard Nessler, and Sepp Hochreiter.  GANs Trained by a Two Time-Scale Update Rule Converge to a Local Nash Equilibrium.  NeurIPS, 2017.
  * [6] Crismeire Isbaex and Ana Margarida Coelho.  The potential of sentinel-2 satellite images for land-cover/land-use and forest biomass estimation: A review.  In Forest Biomass. IntechOpen, Rijeka, 2021.
  * [7] Michael D. King, Steven Platnick, W. Paul Menzel, Steven A. Ackerman, and Paul A. Hubanks.  Spatial and temporal distribution of clouds observed by modis onboard the terra and aqua satellites.  IEEE Transactions on Geoscience and Remote Sensing, 2013.
  * [8] Thaís Almeida Lima, René Beuchle, Andreas Langner, Rosana Cristina Grecchi, Verena C. Griess, and Frédéric Achard.  Comparing sentinel-2 msi and landsat 8 oli imagery for monitoring selective logging in the brazilian amazon.  Remote Sensing, 2019.
  * [9] Fiona Ngadze, Kudzai Shaun Mpakairi, Blessing Kavhu, Henry Ndaimani, and Monalisa Shingirayi Maremba.  Exploring the utility of sentinel-2 msi and landsat 8 oli in burned area mapping for a heterogenous savannah landscape.  PLOS ONE, 2020.
  * [10] Zhou Wang, Alan C. Bovik, Hamid R. Sheikh, and Eero P. Simoncelli.  Image quality assessment: from error visibility to structural similarity.  IEEE Transactions on Image Processing, 2004.
  * [11] World Resources Institute.  2017 Was the Second-Worst Year on Record for Tropical Tree Cover Loss.  https://www.wri.org.
  * [12] Richard Zhang, Phillip Isola, Alexei A. Efros, Eli Shechtman, and Oliver Wang.  The unreasonable effectiveness of deep features as a perceptual metric, 2018.