# Gravitational-wave glitches: resonant islands and frequency jumps in non-
integrable extreme-mass-ratio inspirals

Kyriakos Destounis<EMAIL_ADDRESS>Theoretical
Astrophysics, IAAT, University of Tübingen, 72076 Tübingen, Germany  Kostas D.
Kokkotas  Theoretical Astrophysics, IAAT, University of Tübingen, 72076
Tübingen, Germany

###### Abstract

The detection of gravitational waves from extreme-mass-ratio inspirals with
upcoming space-borne detectors will allow for unprecedented tests of general
relativity in the strong-field regime. Aside from assessing whether black
holes are unequivocally described by the Kerr metric, such detections may
place constraints on the degree of spacetime symmetry. In particular,
depending on exactly how a hypothetical departure from the Kerr metric
manifests, the Carter symmetry, which implies the integrability of the
geodesic equations, may be broken. Here, we examine the gravitational
waveforms associated with non-integrable extreme-mass-ratio inspirals
involving a small-mass companion and a supermassive compact object of general
relativity, namely the Manko-Novikov spacetime. We show that the waveforms
displays sudden frequency jumps, when the companion transverses resonant
islands. These findings demonstrate that such abrupt manifestations in the
gravitational-wave frequencies are generic, have a genuine astrophysical
origin and function as a distinctive signature of chaotic phenomena in
extreme-mass-ratio binaries.

##  I Introduction

The direct detection of gravitational waves (GWs) by ground-based detectors,
such as the Laser Interferometer Gravitational-wave Observatory (LIGO) and
Virgo [1, 2], has unlocked a new realm into the exploration of coalescing
compact objects in the strong-field regime. The successors to ground-based
detectors, the Laser Interferometer Space Antenna (LISA) and Taiji [3, 4],
will redefine our understanding of GW astronomy by probing the low-frequency
regime of binary mergers [5]. One of their primary targets are extreme-mass-
ratio inspirals (EMRIs).

These binaries involve a supermassive compact object (the primary) and a
small-mass companion (the secondary). The fact that one of the binary
components has, essentially, negligible mass compared to the other, implies
that the secondary object transverses the gravitational field generated by the
primary following geodesic motion, at zeroth order. Including higher-order
self-force corrections to the equations of motion of the secondary, obtained
from black-hole (BH) perturbation theory [6, 7], leads to an adiabatic drift
between geodesics due to radiation reaction. Taking into consideration that
EMRIs evolve slowly and complete $\sim 10^{5}$ cycles in the strong-field
regime means that they will provide unmatched information about the governing
theory of gravity [8, 9, 10, 11, 12, 13].

While a plethora of BH solutions currently exist in general relativity (GR),
astrophysical BHs are expected to be described by the Kerr metric [14, 15].
The governing equations of an object of negligible, but conserved, rest mass
orbiting around a Kerr BH are integrable, due the stationarity, axisymmetry
and the existence of a hidden symmetry; the Carter constant [16]. Since
integrability implies the absence of chaotic phenomena [17, 18], an EMRI
involving a Kerr primary (or any other compact object with integrable
geodesics) will produce a waveform with slowly drifting fundamental
frequencies and harmonics. If, however, the primary’s geometry is
significantly distorted due to the presence of accretion disks [19, 20]
(though their influence is expected to be negligible [21]), or GR does not
provide an exact description of geometries surrounding compact objects, then
the Kerr metric may not adequately describe supermassive BHs. Including
slight, or otherwise, modifications to the Kerr description provides a
reasonable framework to deal with the aforementioned subtleties with a
drawback; integrability may be broken (though see [22, 23, 24]).

When a non-integrable perturbation is introduced into the Hamiltonian, the
resulting system may exhibit regions where chaotic phenomena arise, in accord
with the Kolmogorov-Arnold-Moser (KAM) and Poincaré-Birkhoff theorems [25, 26,
27, 28]. Such ergodic regions may either surround resonances [29] or form
close to the plunge [30]. In EMRIs, though, chaotic motion itself is not
expected to be a prominent effect, therefore one must search for alternative
generic features that clearly distinguish the existence of chaos. It turns out
that a solid indicator of non-integrability is the passage through resonant
islands [31, 32, 33, 34, 35, 36].

A resonance occurs when two or more characteristic frequencies of a system
match in integer-number ratios. When this is the case, the resonant condition
is fulfilled [18] and the orbit becomes periodic. The resonant orbit no longer
samples the available phase space, in contrast with generic orbits which cover
it densely. Although resonances generically occur in integrable systems, such
as Kerr, non-integrable Kerr-like compact objects (henceforth referred to as
‘non-Kerr’ objects) possess regimes of non-zero phase-space volume called
resonant islands (‘Birkhoff islands’), within which geodesics share the same
rational ratio of orbital frequencies [27, 28].

Resonant orbits are ubiquitous in EMRIs and pose quite a challenge in data
analysis, by virtue of their topological structure which differs significantly
from that of generic orbits. When the secondary object crosses a resonance,
the flux of GWs can be diminished or enhanced, leading to a shift in its
trajectory and the phase of the resulting waveform [37, 38, 39]. Inspiraling
objects may spend, roughly, a few hours [40, 41] or a few dozens of orbital
cycles ‘near’ resonance [42] (the resonant condition is satisfied within an
arbitrarily short range) for supermassive BHs with masses $\sim
10^{6}M_{\odot}$, though the observation of a very long-lived (sustained)
resonance is highly unlike, if not practically impossible [43].

Resonances surrounded by Birkhoff islands of non-Kerr objects, on the other
hand, lead to the entrapment of small-mass companions in ‘perfect’ resonance
(the resonant condition is precisely satisfied) for a time comparable to the
width of the island, and should introduce even greater difficulties in
parameter estimation. At the orbital level, occupancy in a resonant island
translates to a dynamical plateau of the ratio of radial and polar frequencies
[18] which can last for up to a week or a few hundred orbital cycles [32, 34].
Surprisingly, even sustained resonances are likely to occur in non-Kerr EMRIs,
where the trapping time of the companion inside a resonant island is, roughly,
comparable to the inspiral timescale [44].

Recently, it was shown that the appearance of a dynamical plateau in the
orbital frequencies of a non-Kerr EMRI induces a substantial ‘glitch’ to the
frequency evolution of a GW detected by LISA [45]. These abrupt frequency
jumps can be understood as direct imprints of chaotic phenomena in EMRIs,
which in turn may either designate the violation of the ‘no-hair’ theorem [32,
45], the existence of an $N>2$-body inspiral [46, 47] or the presence of spin
in the secondary inspiraling object [48, 49, 50].

Here, we consider an EMRI consisting of a small-mass secondary object
inspiraling into a non-Kerr compact object described by a particular subclass
of the Manko-Novikov spacetime [51] which is an exact solution to the vacuum
Einstein equations and includes a generic deviation to the quadrupole mass
moment. Although the orbital dynamics of such bumpy spacetime has been
extensively analyzed [30, 32], we show that the passage through a resonant
island corresponds to a sudden jump in the orbital elements of the inspiral.
This abrupt phenomenon manifests itself into the gravitational waveform as a
frequency glitch, in an identical manner to those found in [45]. Hence, we
demonstrate that GW glitches are generic features of non-integrable EMRIs, are
intrinsically linked to the fundamental spacetime symmetries of the system and
signify the presence of chaotic phenomena during the inspiral.

##  II The Manko-Novikov spacetime

The goal of this study is to examine the universality of GW glitches in non-
integrable EMRIs, such as those found recently in [45]. To achieve this, we
utilize a metric, namely the Manko-Novikov (MN) spacetime [51], which is an
exact solution to the vacuum Einstein field equations and generalizes the
Kerr-Newman metric into a multiparametric family. Although this solution
depends, in general, on an infinite amount of deformation parameters related
to the multipole moments of the gravitational field, we will adopt a
particular subclass [30, 32] which admits one extra parameter, besides the
mass $M$ and spin $S$, which controls the integrability of the geodesic
equations. This new deformation parameter, $q$, is a dimensionless number that
measures the deviation of the mass quadrupole moment of the gravitational
field from that of the corresponding Kerr metric. The first three, non-
vanishing, mass and current multipole moments are, then, characterized by
$M,\,S,\,q$ as [30, 32]

$\displaystyle M_{0}$ $\displaystyle=M,\,\,\,S_{1}=\chi
M^{2},\,\,\,M_{2}=-M^{3}\left(\chi^{2}+q\right),$ (1)

where $\chi\equiv S/M^{2}$ is the dimensionless spin parameter, while higher
order moments are characterized by higher order polynomials of $q$. A general
stationary, axisymmetric, vacuum spacetime can be written in the Weyl-
Papapetrou coordinates $(t,\rho,z,\phi)$ as [30, 32]

$ds^{2}=-f(dt-\omega
d\phi)^{2}+\frac{1}{f}\left[e^{2\gamma}(d\rho^{2}+dz^{2})+\rho^{2}d\phi^{2}\right].$
(2)

For the MN spacetime, $f,~{}\omega,~{}\gamma$ are functions of the prolate
spheroidal coordinates $(x,y)$ with the following form:

$\displaystyle f$ $\displaystyle=e^{2\psi}\frac{A}{B},$ (3)
$\displaystyle\omega$
$\displaystyle=2ke^{-2\psi}\frac{C}{A}-4k\frac{\alpha}{1-\alpha^{2}},$ (4)
$\displaystyle e^{2\gamma}$
$\displaystyle=e^{2\gamma^{\prime}}\frac{A}{(x^{2}-1)(1-\alpha^{2})^{2}},$ (5)
$\displaystyle A$ $\displaystyle=(x^{2}-1)(1+a~{}b)^{2}-(1-y^{2})(b-a)^{2},$
(6) $\displaystyle B$
$\displaystyle=[(x+1)+(x-1)a~{}b]^{2}+[(1+y)a+(1-y)b]^{2},$ (7) $\displaystyle
C$ $\displaystyle=(x^{2}-1)(1+a~{}b)[(b-a)-y(a+b)]$
$\displaystyle+(1-y^{2})(b-a)[(1+a~{}b)+x(1-a~{}b)],$ (8) $\displaystyle\psi$
$\displaystyle=\beta\frac{P_{2}}{R^{3}},$ (9) $\displaystyle\gamma^{\prime}$
$\displaystyle=\ln{\sqrt{\frac{x^{2}-1}{x^{2}-y^{2}}}}+\frac{3\beta^{2}}{2R^{6}}(P_{3}^{2}-P_{2}^{2})$
$\displaystyle+\beta\left(-2+\displaystyle{\sum_{\ell=0}^{2}}\frac{x-y+(-1)^{2-\ell}(x+y)}{R^{\ell+1}}P_{\ell}\right),$
(10) $\displaystyle a$
$\displaystyle=-\alpha\exp{\left[-2\beta\left(-1+\displaystyle{\sum_{\ell=0}^{2}}\frac{(x-y)P_{\ell}}{R^{\ell+1}}\right)\right]},$
(11) $\displaystyle b$
$\displaystyle=\alpha\exp{\left[2\beta\left(1+\displaystyle{\sum_{\ell=0}^{2}}\frac{(-1)^{3-\ell}(x+y)P_{\ell}}{R^{\ell+1}}\right)\right]},$
(12)

where $R=\sqrt{x^{2}+y^{2}-1}$ and $P_{\ell}=P_{\ell}(xy/R)$ are the Legendre
polynomials of order $l$ given by

$P_{\ell}(w)=\frac{1}{2^{\ell}\ell!}\left(\frac{d}{dz}\right)^{\ell}(w^{2}-1)^{\ell}.$
(13)

The line element (2) is related to the parameters $\alpha,\,k,\,\beta$ which
are defined as

$\alpha=\frac{-1+\sqrt{1-\chi^{2}}}{\chi},\,\,\,k=M\frac{1-\alpha^{2}}{1+\alpha^{2}},\,\,\,\beta=q\frac{M^{3}}{k^{3}}.$
(14)

Consequently, choosing $M,\,\chi$ and $q$ completely defines the MN metric.
The relation between prolate spheroidal coordinates $(x,y)$ and quasi-
cylindrical coordinates $(\rho,z)$ is given by [32, 52]

$\rho=k\sqrt{(x^{2}-1)(1-y^{2})},\quad z=kxy,$ (15)

while the inverse transformation reads

$\displaystyle x$
$\displaystyle=\frac{1}{2k}\left(\sqrt{\rho^{2}+\left(z+k\right)^{2}}+\sqrt{\rho^{2}+\left(z-k\right)^{2}}\right),$
(16) $\displaystyle y$
$\displaystyle=\frac{1}{2k}\left(\sqrt{\rho^{2}+\left(z+k\right)^{2}}-\sqrt{\rho^{2}+\left(z-k\right)^{2}}\right).$
(17)

In turn, the relation between $(x,y)$, $(\rho,z)$ and the standard Boyer-
Lindquist coordinates $(r,\theta)$ is given by

$r=kx+M,\quad\theta=\cos^{-1}y$ (18)

and

$\displaystyle r/M$
$\displaystyle=1+\sqrt{-\sigma_{2}+\sqrt{\sigma_{2}^{~{}2}+\sigma_{1}(z/M)^{2}}},$
(19) $\displaystyle\theta$
$\displaystyle=\cos^{-1}\left({\frac{z/M}{(r/M)-1}}\right),$ (20)

where

$\sigma_{1}=\chi^{2}-1,\quad\sigma_{2}=\frac{\sigma_{1}-(\rho/M)^{2}-(z/M)^{2}}{2}.$

The MN spacetime, having a multipolar structure different from that of Kerr,
is not a BH solution, in accordance with the no-hair theorem. In fact, it
lacks a compact event horizon and contains closed time-like curves, where
$g_{\phi\phi}<0$, exterior to its horizon which is broken in the equatorial
plane by a circular line singularity. For non-zero spin parameters, the metric
contains an ergoregion, determined by $g_{tt}=0$, which generally admits a
multiple lobed structure [30]. For $q=0$, the Kerr metric is recovered, while
a choice of $q>0\,(q<0)$ represents an oblate (prolate) perturbation of the
Kerr metric.

Regardless of its pathological nature, the MN metric is a prime candidate to
describe the gravitational field of a stationary axisymmetric compact object
with a well-defined multipolar structure which deviates from that of Kerr. In
particular, it does not admit a Carter constant, as opposed to Kerr, and thus
provides a proper testing ground for the no-hair theorem, and its consequences
on the waveforms of non-integrable inspirals. In what follows, we will
restrict our numerical evolutions into regions where the near-horizon
pathological features of the MN spacetime are not present.

##  III Short timescale EMRI modeling: geodesic motion

Due to GW emission, the orbital dynamics that governs an EMRI is not purely
geodesic. Nevertheless, at zeroth-order and for short timescales, we can
neglect the structure of the secondary and adopt a point-particle
approximation, with its motion on the gravitational field created by the
primary being described by the geodesic equation

$\frac{d^{2}x^{\kappa}}{d\tau^{2}}+\Gamma^{\kappa}_{\lambda\nu}\frac{dx^{\lambda}}{d\tau}\frac{dx^{\nu}}{d\tau}=0,$
(21)

where $\Gamma^{\alpha}_{\lambda\nu}$ are the Christoffel symbols, $x^{\kappa}$
the $4$-position vector and $\tau$ the proper time. Equivalently, one could
operate in the Lagrange formulation to describe geodesic motion. The equations
of motion for a secondary object with rest mass $\mu$ can be derived, via the
Euler-Lagrange equations, by considering a Lagrangian

$\mathcal{L}=\frac{1}{2}\mu\,g_{\lambda\nu}\dot{x}^{\lambda}\dot{x}^{\nu},$
(22)

where $g_{\lambda\nu}$ are the metric tensor components of (2) and an overhead
dot denotes differentiation with respect to proper time.

The MN spacetime, being stationary and axisymmetric, immediately leads to the
conservation of the secondary’s temporal and azimuthal momenta
$p_{t},\,p_{\phi}$, with $p_{\lambda}=\mu\,g_{\lambda\nu}\dot{x}^{\nu}$,
through the equations of motion $\dot{p}_{t}=0=\dot{p}_{\phi}$. There are,
therefore, two constant of motion; the energy $E$ and the $z$-component of the
angular momentum $L_{z}$, which are given by

$\displaystyle E$
$\displaystyle=-\mu\left(g_{tt}\dot{t}+g_{t\phi}\dot{\phi}\right)=\mu\,f\left(\dot{t}-\omega\dot{\phi}\right),$
(23) $\displaystyle L_{z}$
$\displaystyle=\mu\left(g_{t\phi}\dot{t}+g_{\phi\phi}\dot{\phi}\right)=\mu\left(f\omega\left(\dot{t}-\omega\dot{\phi}\right)+\frac{\rho^{2}}{f}\dot{\phi}\right).$
(24)

Together with the preservation of the secondary’s rest mass $\mu$, which
results from the conservation of $4$-velocity along a geodesic
$g_{\lambda\nu}\dot{x}^{\lambda}\dot{x}^{\nu}=-1$, one obtains three constants
of motion. Unfortunately, these integrals do not suffice to decouple the
secondary’s motion in the $(\rho,z)$-plane, in contrast to geodesics in Kerr
($q=0$) which possess a fourth integral of motion, the Carter constant, and
the equations of motion in the $(\rho,z)$-plane are separable. Consequently,
to obtain the motion in MN spacetime, one requires use of the coupled second-
order differential equations (21):

$\displaystyle
2\,g_{\rho\rho}\,\ddot{\rho}+2\,\dot{\rho}\,\dot{z}\,\partial_{z}g_{\rho\rho}+\dot{\rho}^{2}\,\partial_{\rho}g_{\rho\rho}-\dot{t}^{2}\,\partial_{\rho}g_{tt}-2\,\dot{t}\,\dot{\phi}\,\partial_{\rho}g_{t\phi}$
$\displaystyle-\dot{z}^{2}\,\partial_{\rho}g_{zz}-\dot{\phi}^{2}\,\partial_{\rho}g_{\phi\phi}=0,$
(25) $\displaystyle
2\,g_{zz}\,\ddot{z}+2\,\dot{\rho}\,\dot{z}\,\partial_{\rho}g_{zz}-\dot{\rho}^{2}\,\partial_{z}g_{\rho\rho}-\dot{t}^{2}\,\partial_{z}g_{tt}-2\,\dot{t}\,\dot{\phi}\,\partial_{z}g_{t\phi}$
$\displaystyle+\dot{z}^{2}\,\partial_{z}g_{zz}-\dot{\phi}^{2}\,\partial_{z}g_{\phi\phi}=0,$
(26)

where

$\displaystyle g_{tt}=-f,\,g_{t\phi}=\omega
f,\,g_{\rho\rho}=g_{zz}=\frac{e^{2\gamma}}{f},\,g_{\phi\phi}=\frac{\rho^{2}}{f}-\omega^{2}f,$
(27)

and the $t$, $\phi$-momenta are given by

$\displaystyle\dot{t}$
$\displaystyle=\frac{Eg_{\phi\phi}+L_{z}g_{t\phi}}{\mu(g_{t\phi}^{2}-g_{tt}g_{\phi\phi})}=\frac{\omega
f}{\mu\,\rho^{2}}\left[{L_{z}}+{\omega
E}\left(\frac{\rho^{2}}{\omega^{2}f^{2}}-1\right)\right],$ (28)
$\displaystyle\dot{\phi}$
$\displaystyle=-\frac{Eg_{t\phi}+L_{z}g_{tt}}{\mu(g_{t\phi}^{2}-g_{tt}g_{\phi\phi})}=\frac{f}{\mu\,\rho^{2}}(L_{z}-\omega
E).$ (29)

Eqs. (25) and (26) are sufficient to fully describe the geodesic motion of the
secondary in the MN spacetime, while the evolution of the $t,\,\phi$
coordinates can be obtained by integrating (28) and (29). The combination of
(28), (29) and the conservation of $4$-velocity leads to the constraint
equation [32]

$\dot{\rho}^{2}+\dot{z}^{2}+V_{\textrm{eff}}(\rho,z)=0,$ (30)

with

$\displaystyle
V_{\textrm{eff}}(\rho,z)\equiv\frac{e^{-2\gamma}}{\mu^{2}}\left[\mu^{2}f-E^{2}+\left(\frac{f}{\rho}(L_{z}-\omega
E)\right)^{2}\right]$ (31)

a $2$-dimensional effective potential. The motion in the $(\rho,z)$-plane may,
therefore, be thought of as bound motion in the effective potential (31). The
curve defined by $V_{\text{eff}}=0$ is called the curve of zero velocity
(CZV), since whenever a geodesic orbit reaches it the velocity components in
the $(\rho,z)$-plane vanish. Equation (30), not only provides an understanding
of the general properties of bound orbits in the MN spacetime, but also can be
used as a cross-check to verify the convergence of geodesic evolution. For an
exhaustive analysis of the properties of the effective potential as well as
its dependence on various parameters, we refer the reader to [30, 32, 34, 33].

To numerically integrate the second-order equations (25) and (26), we
generally need to prescribe appropriate initial conditions of the form
$(\rho(0),\dot{\rho}(0),z(0),\dot{z}(0))$, besides the parameters describing
the particulars of the system, which establish some initial position and
velocity to the secondary, at $\tau=0$, within the allowed region of bound
motion on the $(\rho,z)$-plane defined by the CZV. We reduce the initial
conditions to $(\rho(0),\dot{\rho}(0),z(0))$, with the remaining condition
$\dot{z}(0)$ being inferred from Eq. (30). Despite the fact that the equations
of motion are written with respect to proper time for simplicity, GW detectors
operate with respect to the inertial time $t$ at very large distances from the
EMRI source. This is merely an inconvenience since one can recast (25), (26)
and (30) to evolve with respect to $t$ by using (28) and the chain rule, e.g.
$\dot{\rho}=d\rho/d\tau=\dot{t}\,d\rho/dt$.

##  IV Long timescale EMRI modeling: radiation reaction and the kludge scheme

GR predicts that astrophysically-captured compact objects emit radiation in
accordance with the Einstein equations. In the presence of gravitational
radiation, the integrals of motion cannot be considered constant anymore, but
rather evolve as energy and angular momentum is gradually carried away by GWs.
Owing to their small mass ratio, EMRIs radiate adiabatically during very long
timescales. As discussed above, on short timescales the orbit is approximately
undergoing geodesic motion. By taking into account the structure of the
secondary, which induces a gravitational self-force, one can calculate
gravitational perturbations at the primary’s horizon and infinity to obtain
the gravitational waveform and the corresponding fluxes of energy and angular
momentum. The gravitational self-force, then, leads to some evolution
equations for $E$ and $L_{z}$ and forms the basis for the adiabatic
approximation introduced in [53, 54], which is valid only when the change in
the momenta is sufficiently small over a single orbit. Over short timescales,
the system follows an orbit that neglects radiative backreaction, while over
long timescales, the inspiral behaves as a ‘flow’ of the system through a
sequence of successively damped geodesics.

Although these fluxes can be completely understood through full-blown self-
force calculations, a description of gravitational self-force, up to first-
order terms, has been achieved only recently for Kerr [7, 55]. To extend these
results to bumpy BHs remains a rather difficult and demanding task. Therefore,
to simulate the dissipative effect of $E$ and $L_{z}$ we employ the hybrid
kludge scheme, introduced in [56, 57], which mixes a notion of short-timescale
motion with an approximate description of the long-timescale radiation
effects. This scheme is based on combining exact relativistic expressions for
the evolution of the orbital elements with approximate post-Newtonian (PN)
formulae for energy and angular momentum fluxes. Hence, radiation reaction is
implicitly included by the particulars of the PN coefficients, and all we need
to do is solve the “instantaneous” geodesic equations (21) augmented with the
appropriate PN formulae.

The kludge approximation has been thoroughly tested, is widely utilized for
the study of EMRIs and has been proven to be remarkably faithful, thus
provides a trustworthy mechanism to generate trajectories of generic EMRIs
particularly quickly [57].

Here, since we are dealing with a spacetime with an anomalous quadrupole
moment (see Eq. (1)), we augment the fluxes with an additional PN term which
represents the effect of $q$ on the evolution of $E$ and $L_{z}$ [12, 30, 32],
through $M_{2}$. Up to second PN order, the kludge equations for the evolution
of an inspiral read [57]

$\displaystyle\bigg{<}\frac{dE}{dt}\bigg{>}$
$\displaystyle=-\frac{32}{5}\frac{\mu^{2}}{M^{2}}\left(\frac{M}{p}\right)^{5}(1-e^{2})^{3/2}\Bigg{[}g_{1}(e)-\chi\left(\frac{M}{p}\right)^{3/2}g_{2}(e)\cos\iota-\left(\frac{M}{p}\right)g_{3}(e)+\pi\left(\frac{M}{p}\right)^{3/2}g_{4}(e)$
$\displaystyle-\left(\frac{M}{p}\right)^{2}g_{5}(e)-(M_{2}/M^{3})\left(\frac{M}{p}\right)^{2}g_{6}(e)+\frac{527}{96}(M_{2}/M^{3})\left(\frac{M}{p}\right)^{2}\sin^{2}\iota\Bigg{]},$
(32) $\displaystyle\bigg{<}\frac{dL_{z}}{dt}\bigg{>}$
$\displaystyle=-\frac{32}{5}\frac{\mu^{2}}{M}\left(\frac{M}{p}\right)^{7/2}(1-e^{2})^{3/2}\Bigg{[}g_{9}(e)\cos\iota+\chi\left(\frac{M}{p}\right)^{3/2}\\{g_{10}^{a}(e)-\cos^{2}\iota
g_{10}^{b}(e)\\}-\left(\frac{M}{p}\right)g_{11}(e)\cos\iota$
$\displaystyle+\pi\left(\frac{M}{p}\right)^{3/2}g_{12}(e)\cos\iota-\left(\frac{M}{p}\right)^{2}g_{13}(e)\cos\iota-(M_{2}/M^{3})\left(\frac{M}{p}\right)^{2}\cos\iota\left(g_{14}(e)-\frac{45}{8}\,\sin^{2}\iota\right)\Bigg{]},$
(33) $\displaystyle g_{1}(e)$
$\displaystyle=1+\frac{73}{24}e^{2}+\frac{37}{96}e^{4},\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad
g_{2}(e)=\frac{73}{12}+\frac{823}{24}e^{2}+\frac{949}{32}e^{4}+\frac{491}{192}e^{6},\
$ $\displaystyle g_{3}(e)$
$\displaystyle=\frac{1247}{336}+\frac{9181}{672}e^{2},\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\,\,\,\,g_{4}(e)=4+\frac{1375}{48}e^{2},\
$ $\displaystyle g_{5}(e)$
$\displaystyle=\frac{44711}{9072}+\frac{172157}{2592}e^{2},\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\,g_{6}(e)=\frac{33}{16}+\frac{359}{32}e^{2},\
$ $\displaystyle g_{9}(e)$
$\displaystyle=1+\frac{7}{8}e^{2},\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\,\,\,g_{10}(e)=\frac{61}{12}+\frac{119}{8}e^{2}+\frac{183}{32}e^{4},\
$ (34) $\displaystyle g_{11}(e)$
$\displaystyle=\frac{1247}{336}+\frac{425}{336}e^{2},\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\,\,\,\,\,g_{12}(e)=4+\frac{97}{8}e^{2},\
$ $\displaystyle g_{13}(e)$
$\displaystyle=\frac{44711}{9072}+\frac{302893}{6048}e^{2},\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\,\,\,\,\,g_{14}(e)=\frac{33}{16}+\frac{95}{16}e^{2},\
$ $\displaystyle g_{10}^{a}(e)$
$\displaystyle=\frac{61}{24}+\frac{63}{8}e^{2}+\frac{95}{64}e^{4},\
\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\,\,\,\,g_{10}^{b}(e)=\frac{61}{8}+\frac{91}{4}e^{2}+\frac{461}{64}e^{4},\
$

where the orbital elements $e,p,\iota$ (where $p$ is _not_ to be confused with
the momentum), corresponding to the eccentricity, semi-latus rectum and
inclination, respectively, are given by

$e=\frac{r_{+}-r_{-}}{r_{+}+r_{-}},\quad\quad
p=\frac{2r_{+}r_{-}}{r_{+}+r_{-}},\quad\quad\iota=\frac{\theta_{+}-\theta_{-}}{2},$
(35)

where $r_{+}(r_{-})$ and $\theta_{+}(\theta_{-})$ denote the periapsis
(apoapsis) and maximum (minimum) angular values, along the orbit,
respectively. At this point we need to admit the following: Eqs. (IV) and (IV)
are only valid for Kerr BHs. The augmentation of those fluxes with additional
terms to account for the anomalous quadrupoles of MN spacetime [12, 30, 31,
32], which is also a solution to GR, serves as a model to understand the
qualitative behavior of orbits but does not capture the exact dynamics of the
EMRI. Furthermore, we note that we do not attempt to evolve the other
components of the angular momentum or some generalization of the Carter
constant in MN spacetime, but only introduce dissipation to the energy and
azimuthal angular momentum as in [30, 31, 32, 45, 44]. Even so, since we will
only evolve the EMRIs for small intervals of time compared to the radiation
reaction timescale, and since we are interested only in particular
phenomenological imprints, we do not expect that the inclusion of the
aforementioned effects will dramatically obscure the following results.

To construct the inspiral we numerically integrate Eqs. (25) and (26),
utilizing (28) and (29), augmented with the PN fluxes (IV) and (IV) and assume
linear variations for the energy and angular momentum as in [58]

$\displaystyle
E_{1}=\frac{E_{0}}{\mu}+\bigg{<}\frac{dE}{dt}\bigg{>}\bigg{|}_{0}N_{r}\,T_{r},$
(36) $\displaystyle
L_{z,1}=\frac{L_{z,0}}{\mu}+\bigg{<}\frac{dL_{z}}{dt}\bigg{>}\bigg{|}_{0}N_{r}\,T_{r},$
(37)

where $E_{0},\,L_{z,0}$ are the initial energy and angular momentum,
respectively, $\langle{dE}/{dt}\rangle|_{0},\,\langle{dLz}/{dt}\rangle|_{0}$
are the fluxes calculated at the beginning of the inspiral and $T_{r}$ is the
time that the orbit takes to travel from the periapsis to apoapsis and back.
Eqs. (35), and thus the fluxes (36) and (37), are updated every $N_{r}$ cycles
for the whole EMRI evolution. To further guarantee numerical accuracy when
radiation reaction is included, we calculate the $4$-velocity in each update
of the fluxes and check its conservation along a geodesic evolution with
initial conditions the energy, angular momentum, position and velocity at
every update time step. For all simulations presented herein, we find that
they are satisfied to within a part in $\sim 10^{9}$ for the first $10^{4}$
crossings through the equatorial plane $z=0$.

##  V Imprints of non-integrability

Due to stationarity and axisymmetry, the Kerr metric admits two integrals of
motion which respectively imply the conservation of energy and angular
momentum of relativistic particles. Together with the preservation of the
particle’s rest mass and the Carter constant, geodesic motion in Kerr is
integrable and the orbits are bound in $2$-dimensional tori of the
$4$-dimensional phase space. In general, orbits oscillate between the two
available degrees of freedom with characteristic frequencies $\omega_{1}$ and
$\omega_{2}$. The ratio of these frequencies,
$\nu_{\theta}=\omega_{1}/\omega_{2}$, is called the rotation number, and
carries pertinent information about the features of the orbit. If the ratio of
the frequencies is an irrational number, the corresponding orbit winds around
its torus, covering it densely. Such orbits are called quasiperiodic. In
special cases, the ratio of frequencies becomes a rational number $n/m$, where
$n,\,m\in\mathbb{N}$, and the orbit returns to its starting point after $m$
windings. These orbits are called periodic and the corresponding torus is
called resonant. The phase-space trajectory of a periodic orbit is not space-
filling, in contrast to quasiperiodic orbits.

When monitoring the trajectory’s successive intersections on a particular
$2$-dimensional surface of section within the phase space (Poincaré map), one
can define closed curves for each torus, named invariant curves [59]. A
quasiperiodic orbit intersects the surface of section repeatedly and covers
densely the corresponding invariant curve which surrounds a fixed central
point of the map. On the other hand, the intersections of periodic orbits form
resonant invariant curves consisting of an infinite number of $m$-multiplets
of periodic points. When the system in study is non-integrable, the resonant
invariant curves nest around periodic stable points, thus forming islands of
stability (Birkhoff islands), in accordance with the Kolmogorov-Arnold-Moser
(KAM) [25, 26] and Poincaré-Birkhoff theorems[27, 28].

By tracking the angles $\vartheta$ formed between subsequent intersections on
a surface of section, relative to some fixed point, we can define a sequence
of rotation numbers,

$\nu_{\theta,N}=\frac{1}{2\pi N}\sum^{N}_{i=1}\vartheta_{i},$ (38)

where in the limit $N\rightarrow\infty$, the sequence (38) converges to
$\nu_{\theta}$. Generally, integrable systems exhibit a monotonous change in
consecutive rotation numbers (rotation curve), while non-integrable systems
exhibit plateaus [18, 32, 36] when the orbit is crossing a Birkhoff island.

Figure 1: Dissipative rotation curves for the EMRI cases (i) (left) and (i)
(right). Both plateaus appear at $\nu_{\theta}=2/3$, designating the crossing
of the $2/3$-resonant island.

The rotation curve has been utilized to distinguish whether various EMRI
systems have integrable equations of motion (see [31, 32, 60, 34, 33, 35, 50,
36]). At the geodesic level, to form a rotation curve, one needs to calculate
successive rotation numbers, with respect to a fixed point
$\boldsymbol{u}_{0}$, by smoothly increasing one of the phase space
parameters, e.g. the radius, while keeping the rest fixed. As explained above,
the calculation of rotation numbers involves the geodesic evolution of an
orbit whilst keeping track of all crossings through a particular surface of
section. The longer we evolve the geodesic, the more intersection points we
obtain, with the number of intersections being proportional to the precision
of the rotation number (e.g. $\sim 10^{4}$ crossings gives $4$ digits of
precision to the resulting rotation number). A typical graph of a rotation
curve in MN spacetime can be found in [32] (see Figure 3 therein).

Rotation curves are particularly helpful to determine the existence of non-
integrability in EMRIs even when radiation reaction is included. After
evolving an inspiral to obtain the trajectory of the secondary around a
primary MN compact object, one can analyze the $\rho$ and $z$ time-series in
the frequency domain with Fourier analysis, and infer the rotation numbers in
successive time segments [32]. This leads to a rotation curve which evolves
with time, and exhibits a plateau during the passage through a resonant
island. An alternative way of producing the rotation curve of an inspiral,
which we utilize here, was recently proposed in [44]. The main idea is to take
various points from the dissipative evolution of $\rho$, $z$, $E$ and $L_{z}$
as initial conditions and evolve those particular points as geodesics for
enough time in order to obtain the respective rotation number of each point.

In what follows, we consider two inspirals, with the following parameters and
initial conditions, so that they pass through the $2/3$-resonance

  1. (i).

$\mu/M=10^{-5},\,[\rho(0),\dot{\rho}(0),z(0)]=(4.23M,0,0)$,

  2. (ii).

$\mu/M=10^{-6},\,[\rho(0),\dot{\rho}(0),z(0)]=(4.129M,0,0)$,

where we fix the parameters $E_{0}=0.95\mu$, $L_{z,0}=3M\mu$, $\chi=0.9$ and
consider a light secondary object of $\mu=1M_{\odot}$ inspiraling into a
primary MN supermassive compact object with deformation parameter $q=0.95$,
where its mass $M$ is defined by the respective mass ratio
($M=10^{5}M_{\odot}$ for (i) and $M=10^{6}M_{\odot}$ for (ii)). As detailed in
Sec. III, the remaining initial condition $\dot{z}(0)$ is fixed by the
constraint equation (30). The total evolution time of the EMRIs are
$t_{\text{evol}}=2\times 10^{6}M\simeq 3.5$ months for case $(\text{i})$ and
$t_{\text{evol}}=10^{6}M\simeq 5.5$ days for case $(\text{ii})$. Both cases
evolve with an approximate orbital period of $T_{r}\sim 2\times 10^{2}M$,
which translates to $\sim 5\times 10^{3}$ and $\sim 10^{4}$ total cycles, for
EMRI $(\text{i})$ and $(\text{ii})$, respectively.

In Fig. 1 we present the dissipative rotation curves for the cases
$(\text{i})$ and $(\text{ii})$ discussed above. A clear plateau appears in
both EMRIs, at roughly half of their respective $t_{\text{evol}}$. The
plateaus are fixed at $\nu_{\theta}=2/3$, which designates the time it takes
for the secondary to cross the $2/3$-resonant island. Occupancy in the island
lasts for $\sim 200$ cycles, in both cases, and should lead to a rather
significant effect on the characteristics of the orbit.

The evolution of the orbital elements during the inspiral is shown in Fig. 2.
Indeed, a very abrupt change occurs, during passage through the resonant
island, with the semi-latus rectum, eccentricity and inclination undergoing a
sudden shift of $\sim 0.6\%$, $1\%$ and $2\%$, respectively. These ‘kicks’ are
quite sudden and occur at entry into the resonant island. Note that away from
resonant islands, we update the fluxes every $10$ cycles, i.e. $N_{r}=10$.
Close to and during passage through an island, we increase the update rate to
$1$ per cycle to properly capture such dynamical phenomenon. This leads to
$\sim 250$ updates around the vicinity of a resonance. We further note that
for the calculation of the orbital elements (and the subsequent fluxes per
update) we evolved geodesics with appropriate initial conditions for
$2.5\times 10^{3}$ and $10^{4}$ cycles and obtained convergence of up to $7$
digits. Since the sudden jumps in the orbital elements occur on around the
third digit, we are reassured that the behavior portrayed in Fig. 2 is
qualitatively valid.

Figure 2: Eccentricity (top row), semi-latus rectum (middle row) and
inclination (bottom row) evolution for the EMRI cases (i) (left column) and
(i) (right column). The abrupt changes occur during the $2/3$-resonant island
crossing for both cases.

We expect that the corresponding kicks resulting from a hypothetical
instantaneous self-force evolution, which includes also conservative terms
affecting the immediate motion, would be smoother, in accordance with the
findings in [37, 39]. Even so, the fact that we can observe any effect in the
orbital elements during passage through a resonant island with the adiabatic
approximation is quite fascinating, although such approach is inadequate in
capturing the full extent of this phenomenon.

In the following sections, we will demonstrate that the abrupt shifts in the
orbital elements, caused during passage through a resonant island, is directly
imprinted in the gravitational waveform of a non-integrable EMRI.

##  VI EMRI Waveform modeling and detector response

Although crossing a resonant island serves as a clear signature of non-
integrability in EMRIs, it is quite challenging for one to extract the
fundamental orbital frequencies from a waveform in order to sketch a rotation
curve. GWs, indeed, carry information regarding the inspiral’s fundamental
orbital frequencies, though in most cases these quantities are hidden within
the waveform as linear combinations of all orbital constituents. In this
section, we investigate the GW emission of the aforementioned demonstrative
EMRIs. We will show that the frequency evolution of an incoming GW, detected
by LISA, contains a clear imprint of the passage through a resonant island, in
the form of an abrupt, but consistent, frequency ‘glitch’, in a similar manner
as those observed in [45].

To model the gravitational waveform of island-crossing EMRIs we adopt the
‘quick and dirty’ numerical kludge scheme, which combines exact particle
trajectories with approximate GW emission [61]. A more sophisticated approach
could involve some (appropriately generalized, see e.g. [62, 63]) Teukolsky
equations to deduce directly the GW characteristics through the Weyl scalars.
However, such an analysis is considerably more complicated and does not help
to elucidate the main new features presented in this work: the universality of
‘glitches’ in non-integrable EMRIs. The numerical kludge is perfectly-suited
for phenomenology and has demonstrated remarkable agreement ($\sim 95\%$) with
Teukolsky-based waveforms [61].

The radiative component of a metric perturbation at large luminosity distance
$d$ from a source $\boldsymbol{T}$, within the Einstein-quadrupole
approximation, can be written, in the transverse-traceless (TT) gauge, as [64,
61, 58]

$h^{\text{TT}}_{ij}=\frac{2}{d}\frac{d^{2}I_{ij}}{dt^{2}}$ (39)

where $I_{ij}$ is the symmetric and trace-free (STF) mass quadrupole
associated to the perturbation,

$I^{ij}=\left[\int d^{3}x\,x^{i}x^{j}\,T^{tt}(t,x^{i})\right]^{\text{STF}},$
(40)

with $t$ the inertial time measured by the detector. The $tt-$component of the
stress-energy tensor for a point-particle with trajectory $\boldsymbol{Z}(t)$
reads [65]

$T^{tt}(t,x^{i})=\mu\delta^{(3)}\left[x^{i}-Z^{i}(t)\right].$ (41)

where $\delta^{(3)}$ denotes the $3$-dimensional Dirac delta distribution.
Utilizing Eqs. (19) and (20), one can transform the quasi-cylindrical
coordinates used in Eq. (2), $(\rho,z)$, to the standard Boyer-Lindquist
coordinates $(r,\theta)$. These, asymptotically reduce to spherical
coordinates, so that suitable pseudo-flat coordinates (not to be confused with
the prolate spheroidal coordinates $(x,y)$)

$\displaystyle
x=r\sin\theta\cos\phi,\,\,\,\,\,\,y=r\sin\theta\sin\phi,\,\,\,\,\,\,z=r\cos\theta,$
(42)

can be identified, and related to the position of a space-borne detector at
infinity. In reality, the detector is located at some finite $d$, so this
prescription is not strictly speaking valid, though has been shown to model
EMRI waveforms quite well when compared to those produced using more
sophisticated approaches [61, 66].

Figure 3: Power spectral density of case (i) (blue) and (ii) (red) EMRIs at a
luminosity distance $d=50\,\text{Mpc}$, over-plotted on the LISA power
spectral density sensitivity curve (dashed black curve). The most prominent
Fourier peaks are presented in the zoomed inlays, while the remaining peaks
correspond to higher harmonics.

An incoming GW can be projected onto its mutually orthogonal $+$ and $\times$
polarization states by introducing two vectors
$\boldsymbol{p}=\boldsymbol{n}\times{\boldsymbol{z}}/|\boldsymbol{n}\times{\boldsymbol{z}}|$
and $\boldsymbol{q}=\boldsymbol{p}\times\boldsymbol{n}$, which are defined in
terms of a unit vector $\boldsymbol{n}$ which points from the source to the
direction of the detector and the unit vector $\boldsymbol{z}$ designating the
spin direction, $\boldsymbol{S}=aM\boldsymbol{z}$. In terms of the
polarization tensors

$\epsilon_{+}^{ij}=p^{i}p^{j}-q^{i}q^{j},\,\,\,\,\,\,\epsilon_{\times}^{ij}=p^{i}q^{j}+p^{j}q^{i},$
(43)

constructed by the spatial orthonormal basis
$(\boldsymbol{n},\boldsymbol{p},\boldsymbol{q})$, the corresponding GW metric
perturbation is

$h^{ij}(t)=\epsilon_{+}^{ij}h_{+}(t)+\epsilon_{\times}^{ij}h_{\times}(t),$
(44)

with

$h_{+}(t)=\frac{1}{2}\epsilon_{+}^{ij}h_{ij}(t),\,\,\,\,\,\,\,h_{\times}(t)=\frac{1}{2}\epsilon_{\times}^{ij}h_{ij}(t).$
(45)

To express the GW components in terms of the position, $Z^{i}(t)$, velocity,
$v^{i}(t)=dZ^{i}/dt$, and acceleration vectors $a^{i}(t)=d^{2}Z^{i}/dt^{2}$,
we use Eqs. (39) and (41) to obtain [58]

$h_{+,\times}(t)=\frac{2\mu}{d}\epsilon^{+,\times}_{ij}\left[a^{i}(t)Z^{j}(t)+v^{i}(t)v^{j}(t)\right].$
(46)

In general, LISA’s response to an incident GW is determined by a vector
$h_{\alpha}$ which depends on the antenna pattern response functions
$F^{+,\times}_{I,II}$ and the $+$ and $\times$ polarization states through
[67]

$h_{\alpha}(t)=\frac{\sqrt{3}}{2}\left[F^{+}_{\alpha}(t)h_{+}(t)+F^{\times}_{\alpha}(t)h_{\times}(t)\right],$
(47)

where $\alpha=(I,II)$ is an index representing the different antenna pattern
functions [68, 69, 67, 24]

$\displaystyle F^{+}_{I}$
$\displaystyle=\frac{1}{2}(1+\cos^{2}\Theta)\cos(2\Phi)\cos(2\Psi)$
$\displaystyle-\cos\Theta\sin(2\Phi)\sin(2\Psi),$ (48) $\displaystyle
F^{\times}_{I}$
$\displaystyle=\frac{1}{2}(1+\cos^{2}\Theta)\cos(2\Phi)\cos(2\Psi)$
$\displaystyle+\cos\Theta\sin(2\Phi)\sin(2\Psi),$ (49) $\displaystyle
F^{+}_{II}$
$\displaystyle=\frac{1}{2}(1+\cos^{2}\Theta)\sin(2\Phi)\cos(2\Psi)$
$\displaystyle+\cos\Theta\cos(2\Phi)\sin(2\Psi),$ (50) $\displaystyle
F^{\times}_{II}$
$\displaystyle=\frac{1}{2}(1+\cos^{2}\Theta)\sin(2\Phi)\sin(2\Psi)$
$\displaystyle-\cos\Theta\cos(2\Phi)\cos(2\Psi),$ (51)

where $(\Theta,\Phi)$ is the sky location and $\Psi$ the polarization of the
source in a precessing detector-based coordinate system. Assuming a fixed
orientation $\boldsymbol{n}=(0,0,1)$ and that the supermassive compact
object’s spin polar and spin azimuthal angles remain fixed at the equatorial
plane for simplicity, the angles $(\Theta,\Phi,\Psi)$ introduced in (48)-(51)
read (see Ref. [58] for formulae in the general case) $\Theta(t)=\pi/3$,
$\Phi(t)=2\pi t/T+\pi/2$, and $\Psi=-2\pi t/T$, where $T=1$ year is the
orbital period of the Earth around the Sun. Although the LISA data stream
consists of two linearly independent channels, we will omit one channel for
simplicity and neglect any noise in the data stream to demonstrate the main
characteristics of detection, therefore, we assume a data stream
$s(t)=h_{\alpha}(t)$.

##  VII Detectability and Frequency evolution

Figure 4: Frequency evolution diagrams for the EMRI cases (i) and (ii)
depicted in the inlays of Fig. 3.

Here, we consider the EMRI sources $(\text{i})$ and $(\text{ii})$ from Sec. V
and assume a luminosity distance $d=50\,\text{Mpc}$ from the LISA detector. To
test the detectability of GWs from these sources we utilize the square root of
the power spectral density (PSD). When discussing a detector, we can define
the square of the PSD $S_{h}(f)$ for the source as [70]

$\sqrt{S_{h}(f)}=2\sqrt{f}|\tilde{h}(f)|,$ (52)

where $f$ the frequency range and $|\tilde{h}(f)|$ the magnitude of the
Fourier transform of the signal. The root PSD is the most frequently plotted
quantity in the literature and has units of $\text{Hz}^{-1/2}$.

Fig. 3 displays the root PSD of EMRIs $(\text{i})$ and $(\text{ii})$ on top of
a typical LISA sensitivity curve. A plethora of peaks are present in the
spectrum which correspond to the fundamental frequencies and higher harmonics
of the GW signals. As expected, such sources produce clearly observable
signals, in the vicinity of LISA’s maximum sensitivity ($\sim 10^{-3}-10^{-2}$
Hz). Although a Kerr (or, in general, any integrable) EMRI produces peaks
which smoothly vary over the available frequency range (see Figure 1 in [45]),
non-integrable EMRIs have richer Fourier structures and form amplitude valleys
whilst crossing the $2/3$-resonant island, where the PSD peaks drop by up to
$2$ orders of magnitude (see inlays in Fig. 3). Consequently, the formation of
valleys in the Fourier peaks of a GW signal designates the passage through an
island and the non-integrability of the EMRI equations of motion, in a precise
manner.

Despite the fact that the Fourier peak structure of non-integrable EMRIs
implies the existence of transient frequency modulation during island
crossings, the frequency dynamics are better captured through spectrograms. To
construct a spectrogram we use the short-time Fourier transform (STFT); a
typical method used to determine the frequency content of local sections of a
signal as it changes over time. By dividing the waveform into shorter segments
of equal length, we compute the Fourier transform separately on each segment,
thus revealing its frequency content. We then plot the change of the
frequencies from one temporal segment to another, eventually obtaining the
frequency evolution. To achieve a smoother transition between segments, we
employ a window of fixed size with an offset, which enables it to slide over
the available signal. At each slide the Fourier transform is evaluated. This
method allows overlapping of time segments and leads to smoother frequency
transitions. One of the pitfalls of the STFT is that it has a fixed
resolution. The width of the window is directly related to how one wants to
represent the signal. A wide window (wideband transform) gives better
frequency resolution though poor time resolution. On the other hand, a narrow
window (narrowband transform) gives good time resolution but poor frequency
resolution. Therefore, the choice of window size is related to an uncertainty
principle as the standard deviation in frequency and time is limited. In what
follows, we employ a window with appropriate offset so that we represent the
frequency evolution of the given signals as optimal as possible, though paying
more attention to frequency resolution.

In Fig. 4 we present the evolution of the fundamental frequencies of the EMRI
cases $(\text{i})$ and $(\text{ii})$ with respect to time. Both cases exhibit
an abrupt, but consistent, glitch when the valley of the peaks is observed.
Smaller mass ratios (such as the one in case $(\text{ii})$) tend to display
more prominent frequency jumps, as a result of the slower passage of the orbit
through successive geodesics. In actuality, even though the glitch in case
$(\text{ii})$ seems more distinguishable than that of case $(\text{i})$, both
inspirals spend a similar amount of cycles in the island and span over a
comparable frequency range. These phenomenon is shared throughout the whole GW
spectrum, though higher harmonics tend to display less prominent effects. The
results presented here, are coherent with those in [45] and demonstrate that
glitches are universal and signify the presence of Birkhoff islands and
chaotic phenomena in non-integrable EMRIs, even when the primary is a solution
to GR.

##  VIII Discussion

In this work, we have studied EMRI orbits and waveforms of a non-Kerr
supermassive compact object with anomalous multipolar structure, with respect
to that of Kerr, namely the MN spacetime. By assuming that such deformations
from a Kerr description may be present in the EMRI composition, we are left
with a non-integrable system of equations of motion. In turn, non-
integrability implies the existence of chaos [18].

Nevertheless, full-blown chaos is not expected to be a prominent effect in
EMRIs since the occurrence of ergodic orbits is rather rare. The phase-space
structure of non-integrable EMRIs is rather similar to that of an integrable
one in the majority of the available phase space [26]. Even so, when the
inspiral is close to a resonance, the structure of phase space is
fundamentally altered. Islands of stability surround resonances [27, 28],
within which the ratio of the radial and polar orbital frequencies, namely the
rotation number, is rational and shared throughout the island. Since
resonances are abundant in astrophysically-captured bodies [40], the
identification of ‘smoking-gun’ phenomena associated with non-integrability
will shed light into the nature of EMRIs and the role of fundamental spacetime
symmetries.

We find that at the orbital level, a clear imprint evinces when the test body
transverses a resonant island; besides the fact that the rotation number
displays a dynamical plateau [31, 32, 33, 34, 35, 36], the orbital elements
undergo rapid changes when the secondary enters the resonant island. We expect
that these abrupt manifestations would be smoother when one considers an
instantaneous self-force evolution [37, 38, 39]. Recently, these plateaus were
associated with an abrupt jump in the GW frequency evolution, thus providing a
definitive signature of non-integrability at the waveform level [45].

We have extended such waveform analysis to an exact vacuum solution of the
Einstein field equations, such as the MN spacetime [51]. The existence of GW
glitches during resonant-island crossings appears to be a generic attribute of
non-integrable EMRIs, and not just an artifact of theory-agnostic spacetimes.
These abrupt jumps in the GW frequencies are sourced by the rapid changes on
the orbital elements, even at the adiabatic-approximation level. Despite the
fact that Kerr EMRIs also undergo resonant kicks [42], their effect can only
be revealed with instantaneous self-force calculations [39], and evince
themselves in the form of accumulated dephasing [37]. One can extrapolate how
strong will be the impact of resonant islands in parameter estimation, with
respect to a typical Kerr resonance, since the orbit experiences prolonged
resonant motion for hundreds of cycles, without taking into account near-
resonance motion.

Ultimately, the adiabatic approximation of EMRIs around resonances should be
replaced by full self-force computations, though such task is highly
challenging, especially for bumpy EMRIs with complicated multipolar structure.
A first step towards this direction could be the augmentation of the numerical
kludge scheme utilized here, with an effective resonance model, such as the
one recently proposed in [77]. Another direction would be to utilize Newtonian
analogs, such as the one in [78], to understand the phenomenology of island
crossings under instantaneous self-force dissipation, though a Newtonian
analog prescription of a non-Kerr spacetime is still lacking.

We have argued that GW glitches are universal in non-Kerr EMRIs and have a
direct astrophysical origin. In our particular case, they function as a
signature of non-Kerrness, though we expect that similar effects should occur
in $N>2$-body EMRIs [46, 47] or when the multipole structure of the small-mass
companion is taken into account [48, 49, 50]. Interestingly, instrumental
glitches have been spotted in LISA pathfinder data [79], though such events in
the acceleration noise are rather abrupt and have a distinct structure with
respect to the glitches shown here. Therefore, our results serve as a
phonomenological imprint of non-integrability and portray a “smoking-gun” of
chaotic phenomena in EMRIs.

## Acknowledgements

The authors are indebted to Arthur G. Suvorov for helpful discussions.

## References

  * [1] B. P. Abbott et al. [LIGO Scientific and Virgo], “GWTC-1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs,” Phys. Rev. X 9, no.3, 031040 (2019)
  * [2] R. Abbott et al. [LIGO Scientific and Virgo], “GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run,” Phys. Rev. X 11, 021053 (2021)
  * [3] P. Amaro-Seoane et al. [LISA Collaboration], “Laser Interferometer Space Antenna,” [arXiv:1702.00786].
  * [4] W. R. Hu and Y. L. Wu, “The Taiji Program in Space for gravitational wave physics and the nature of gravity,” (2017): 685-686. Natl. Sci. Rev. 4, 685 (2017).
  * [5] E. Berti, E. Barausse, I. Cholis, J. García-Bellido, K. Holley-Bockelmann, S. A. Hughes, B. Kelly, E. D. Kovetz, T. B. Littenberg and J. Livas, et al. “Tests of General Relativity and Fundamental Physics with Space-based Gravitational Wave Detectors,” [arXiv:1903.02781 [astro-ph.HE]].
  * [6] E. Poisson, A. Pound and I. Vega, “The Motion of point particles in curved spacetime,” Living Rev. Rel. 14, 7 (2011)
  * [7] L. Barack and A. Pound, “Self-force and radiation reaction in general relativity,” Rept. Prog. Phys. 82, no.1, 016904 (2019)
  * [8] F. D. Ryan, “Gravitational waves from the inspiral of a compact object into a massive, axisymmetric body with arbitrary multipole moments,” Phys. Rev. D 52, 5707 (1995).
  * [9] E. Berti, A. Buonanno, and C. M. Will, “Testing general relativity and probing the merger history of massive black holes with LISA,” Classical and Quantum Gravity 22.18 (2005): S943.
  * [10] K. Glampedakis and S. Babak, “Mapping spacetimes with LISA: inspiral of a test body in a ‘quasi-Kerr’field,” Classical and Quantum Gravity 23.12 (2006): 4167.
  * [11] B. F. Schutz, “Fundamental physics with LISA,” Classical and Quantum Gravity 26.9 (2009): 094020.
  * [12] L. Barack and C. Cutler, “Using LISA EMRI sources to test off-Kerr deviations in the geometry of massive black holes,” Phys. Rev. D 75, 042003 (2007).
  * [13] E. Barausse, E. Berti, T. Hertog, S. A. Hughes, P. Jetzer, P. Pani, T. P. Sotiriou, N. Tamanini, H. Witek and K. Yagi, et al. “Prospects for Fundamental Physics with LISA,” Gen. Rel. Grav. 52, no.8, 81 (2020)
  * [14] R. P. Kerr, “Gravitational field of a spinning mass as an example of algebraically special metrics,” Phys. Rev. Lett. 11, 237-238 (1963)
  * [15] D. C. Robinson, “Uniqueness of the Kerr black hole,” Phys. Rev. Lett. 34.14 (1975): 905.
  * [16] B. Carter, “Global structure of the Kerr family of gravitational fields,” Phys. Rev. 174, 1559 (1968).
  * [17] D. C. Wilkins, “Bound geodesics in the Kerr metric,” Phys. Rev. D 5, 814 (1972).
  * [18] G. Contopoulos, _Order and Chaos in Dynamical Astronomy_ (Springer, New York, 2002).
  * [19] E. Barausse, L. Rezzolla, D. Petroff, and M. Ansorg, “Gravitational waves from extreme mass ratio inspirals in nonpure Kerr spacetimes,” Phys. Rev. D 75, 064026 (2007).
  * [20] P. Suková and O. Semerák, “Free motion around black holes with discs or rings: between integrability and chaos–III,” Mon. Not. Roy. Astron. Soc. 436, 978 (2013).
  * [21] E. Barausse, V. Cardoso and P. Pani, “Can environmental effects spoil precision gravitational-wave astrophysics?,” Phys. Rev. D 89, no.10, 104059 (2014)
  * [22] T. Johannsen, “Regular black hole metric with three constants of motion,” Phys. Rev. D 88, 044002 (2013).
  * [23] G. O. Papadopoulos and K. D. Kokkotas, “Preserving Kerr symmetries in deformed spacetimes,” Class. Quant. Grav. 35, 185014 (2018).
  * [24] C. J. Moore, A. J. K. Chua and J. R. Gair, “Gravitational waves from extreme mass ratio inspirals around bumpy black holes,” Class. Quant. Grav. 34, no.19, 195009 (2017)
  * [25] J. Moser, “On invariant curves of area-preserving mappings of an annulus,” Nachr. Akad. Wiss. Göttingen Math. Phys. Kl. II, 1 (1962).
  * [26] V. I. Arnold, “Proof of a theorem of A. N. Kolmogorov on the invariance of quasi-periodic motions under small perturbations of the Hamiltonian,” Russian Math. Survey 18, 9 (1963).
  * [27] H. Poincaré, “Sur un théorème de géometrie,” Rend. Circ. Mat. Palermo 33, 375 (1912).
  * [28] G. D. Birkhoff, “Proof of Poincaré’s geometric theorem,” Trans. Am. Math. Soc. 14, 14 (1913).
  * [29] V. I. Arnold, _Mathematical Methods of Classical Mechanics_ , (Springer, New York, 1989), 2nd ed.
  * [30] J. R. Gair, C. Li and I. Mandel, “Observable Properties of Orbits in Exact Bumpy Spacetimes,” Phys. Rev. D 77, 024035 (2008).
  * [31] T. Apostolatos, G. Lukes-Gerakopoulos, and G. Contopoulos, “How to observe a non-kerr spacetime using gravitational waves,” Physical review letters 103.11 (2009): 111101.
  * [32] G. Lukes-Gerakopoulos, T. A. Apostolatos, and G. Contopoulos, “Observable signature of a background deviating from the Kerr metric,” Phys. Rev. D 81, 124005 (2010).
  * [33] G. Contopoulos, G. Lukes-Gerakopoulos, and T. A. Apostolatos, “Orbits in a non-Kerr dynamical system,” International Journal of Bifurcation and Chaos 21, 2261 (2011).
  * [34] G. Lukes-Gerakopoulos, G. Contopoulos and T. A. Apostolatos, “Non-Linear Effects in Non-Kerr spacetimes,” Springer Proc. Phys. 157, 129-136 (2014).
  * [35] O. Zelenka and G. Lukes-Gerakopoulos, “Chaotic motion in the Johannsen-Psaltis spacetime,” [arXiv:1711.02442 [gr-qc]].
  * [36] K. Destounis, A. G. Suvorov, and K. D. Kokkotas, “Testing spacetime symmetry through gravitational waves from extreme-mass-ratio inspirals,” Phys. Rev. D 102, 064041 (2020).
  * [37] E. E. Flanagan and T. Hinderer, “Transient resonances in the inspirals of point particles into black holes,” Phys. Rev. Lett. 109, 071102 (2012)
  * [38] É. E. Flanagan, S. A. Hughes and U. Ruangsri, “Resonantly enhanced and diminished strong-field gravitational-wave fluxes,” Phys. Rev. D 89, no.8, 084028 (2014)
  * [39] C. P. L. Berry, R. H. Cole, P. Cañizares and J. R. Gair, “Importance of transient resonances in extreme-mass-ratio inspirals,” Phys. Rev. D 94, no.12, 124042 (2016)
  * [40] J. Brink, M. Geyer and T. Hinderer, “Orbital resonances around Black holes,” Phys. Rev. Lett. 114, no.8, 081102 (2015)
  * [41] J. Brink, M. Geyer and T. Hinderer, “Astrophysics of resonant orbits in the Kerr metric,” Phys. Rev. D 91, no.8, 083001 (2015)
  * [42] U. Ruangsri and S. A. Hughes, “Census of transient orbital resonances encountered during binary inspiral,” Phys. Rev. D 89, no.8, 084036 (2014)
  * [43] M. van de Meent, “Conditions for Sustained Orbital Resonances in Extreme Mass Ratio Inspirals,” Phys. Rev. D 89, no.8, 084033 (2014)
  * [44] G. Lukes-Gerakopoulos and V. Witzany, “Non-linear effects in EMRI dynamics and their imprints on gravitational waves,” [arXiv:2103.06724 [gr-qc]].
  * [45] K. Destounis, A. G. Suvorov and K. D. Kokkotas, “Gravitational-wave glitches in chaotic extreme-mass-ratio inspirals,” Phys. Rev. Lett. 126, 141102 (2021)
  * [46] E. Barausse, L. Rezzolla, D. Petroff and M. Ansorg, “Gravitational waves from Extreme Mass Ratio Inspirals in non-pure Kerr spacetimes,” Phys. Rev. D 75, 064026 (2007)
  * [47] P. Amaro-Seoane, P. Brem, J. Cuadra and P. J. Armitage, “The butterfly effect in the extreme-mass ratio inspiral problem,” Astrophys. J. Lett. 744, L20 (2012)
  * [48] K. Kiuchi and K. Maeda,“Gravitational waves from a chaotic dynamical system,” Phys. Rev. D 70, 064036 (2004).
  * [49] O. Zelenka, G. Lukes-Gerakopoulos and V. Witzany, “Recurrence analysis of spinning particles in the Schwarzschild background,” [arXiv:1903.00360 [gr-qc]].
  * [50] O. Zelenka, G. Lukes-Gerakopoulos, V. Witzany and O. Kopáček, “Growth of resonances and chaos for a spinning test particle in the Schwarzschild background,” Phys. Rev. D 101, no.2, 024037 (2020)
  * [51] Manko V.S., Novikov I.D., “Generalizations of the Kerr and Kerr-Newman metrics possessing an arbitrary set of mass-multipole moments,” Class. Quant. Grav., 9, 2477, (1992)
  * [52] C. Bambi and G. Lukes-Gerakopoulos, “Testing the existence of regions of stable orbits at small radii around black hole candidates,” Phys. Rev. D 87, no.8, 083006 (2013)
  * [53] Y. Mino, Phys. Rev. D 67, 084027 (2003).
  * [54] D. V. Gal’tsov, “Radiation reaction in the Kerr gravitational field,” Journal of Physics A: Mathematical and General 15, 3737 (1982).
  * [55] M. van de Meent, “Gravitational self-force on generic bound geodesics in Kerr spacetime,” Phys. Rev. D 97, no.10, 104033 (2018)
  * [56] L. Barack and C. Cutler, “LISA capture sources: Approximate waveforms, signal-to-noise ratios, and parameter estimation accuracy,” Phys. Rev. D 69, 082005 (2004).
  * [57] J. R. Gair and K. Glampedakis, “Improved approximate inspirals of test bodies into Kerr black holes,” Phys. Rev. D 73, 064037 (2006).
  * [58] P. Canizares, J. R. Gair, and C. F. Sopuerta, “Testing Chern-Simons modified gravity with gravitational-wave detections of extreme-mass-ratio binaries,” Phys. Rev. D 86, 044010 (2012).
  * [59] Lichtenberg A. J. , Lieberman M. A. , “Regular and Chaotic Dynamics,” Springer-Verlag, (1992)
  * [60] G. Lukes-Gerakopoulos, “The non-integrability of the Zipoy-Voorhees metric,” Phys. Rev. D 86, 044013 (2012)
  * [61] S. Babak et al., ““Kludge” gravitational waveforms for a test-body orbiting a Kerr black hole,” Phys. Rev. D 75, 024005 (2007).
  * [62] A. G. Suvorov, “Gravitational perturbations of a Kerr black hole in f (R) gravity,” Phys. Rev. D 99, 124026 (2019).
  * [63] A. G. Suvorov, “Signatures of extended theories of gravity in black hole oscillations,” Acta Phys. Pol. B Proc. Suppl. 13, 283 (2020).
  * [64] K. S. Thorne, “Multipole expansions of gravitational radiation,” Reviews of Modern Physics 52, 299 (1980).
  * [65] P. C. Peters, “Relativistic gravitational bremsstrahlung,” Phys. Rev. D 1, 1559 (1970).
  * [66] C. F. Sopuerta and N. Yunes, “New Kludge Scheme for the Construction of Approximate Waveforms for Extreme-Mass-Ratio Inspirals,” Phys. Rev. D 84, 124060 (2011)
  * [67] L. Barack and C. Cutler, “LISA capture sources: Approximate waveforms, signal-to-noise ratios, and parameter estimation accuracy,” Phys. Rev. D 69, 082005 (2004).
  * [68] T. A. Apostolatos, C. Cutler, G. J. Sussman and K. S. Thorne, “Spin induced orbital precession and its modulation of the gravitational wave forms from merging binaries,” Phys. Rev. D 49, 6274-6297 (1994)
  * [69] C. Cutler, “Angular resolution of the LISA gravitational wave detector,” Phys. Rev. D 57, 7089-7102 (1998)
  * [70] C. J. Moore, R. H. Cole and C. P. L. Berry, “Gravitational-wave sensitivity curves,” Class. Quant. Grav. 32, no.1, 015014 (2015)
  * [71] D. Lynden-Bell, “Galactic nuclei as collapsed old quasars,” Nature 223, 690 (1969)
  * [72] J. Kormendy and D. Richstone, “Inward bound: The Search for supermassive black holes in galactic nuclei,” Ann. Rev. Astron. Astrophys. 33, 581 (1995)
  * [73] m. Miyoshi, J. Moran, J. Hernstein, L. Greenhill, N. Nakai, P. Diamond and M. Inoue, “Evidence for a black hole from high rotation velocities in a subparsec region of NGC4258,” Nature 373, 127-129 (1995)
  * [74] A. M. Ghez, B. L. Klein, M. Morris and E. E. Becklin, “High proper motion stars in the vicinity of Sgr A*: Evidence for a supermassive black hole at the center of our galaxy,” Astrophys. J. 509, 678-686 (1998)
  * [75] A. Eckart and R. Genzel, “Stellar proper motions in the central 0.1 PC of the galaxy,” Mon. Not. Roy. Astron. Soc. 284, 576-598 (1997)
  * [76] R. Abuter et al. [GRAVITY], “Detection of the gravitational redshift in the orbit of the star S2 near the Galactic centre massive black hole,” Astron. Astrophys. 615, L15 (2018)
  * [77] L. Speri and J. R. Gair, Phys. Rev. D 103, no.12, 124032 (2021)
  * [78] A. Eleni and T. A. Apostolatos, “Newtonian analogue of a Kerr black hole,” Phys. Rev. D 101, no.4, 044056 (2020)
  * [79] M. C. Edwards, P. Maturana-Russel, R. Meyer, J. Gair, N. Korsakova and N. Christensen, “Identifying and Addressing Nonstationary LISA Noise,” Phys. Rev. D 102, no.8, 084062 (2020)