# Extreme-mass-ratio inspirals into rotating boson stars:
nonintegrability, chaos, and transient resonances

Kyriakos Destounis, Federico Angeloni, Massimo Vaglio and Paolo Pani
Dipartimento di Fisica, Sapienza Università di Roma, Piazzale Aldo Moro 5,
00185, Roma, Italy  INFN, Sezione di Roma, Piazzale Aldo Moro 2, 00185, Roma,
Italy

###### Abstract

General relativity predicts that black holes are described by the Kerr metric,
which has integrable geodesics. This property is crucial to produce accurate
waveforms from extreme-mass-ratio inspirals. Astrophysical environments,
modifications of gravity and new fundamental fields may lead to nonintegrable
geodesics, inducing chaotic effects. We study geodesics around self-
interacting rotating boson stars and find robust evidence of nonintegrability
and chaos. We identify islands of stability around resonant orbits, where the
orbital radial and polar oscillation frequency ratios, known as rotation
numbers, remain constant throughout the island. These islands are generically
present both in the exterior and the interior of compact boson stars. A
monotonicity change of rotation curves takes place as orbits travel from the
exterior to the interior of the star. Therefore, configurations with neutron-
star-like compactness can support degenerate resonant islands. This anomaly is
reported here for the first time and it is not present in black holes. Such
configurations can also support extremely prolonged resonant islands that span
from the exterior to the interior of the star and are shielded by thick
chaotic layers. We adiabatically evolve inspirals using approximated post-
Newtonian fluxes and find time-dependent plateaus in the rotation curves which
are associated with island-crossing orbits. Crossings of external islands give
rise to typical gravitational-wave glitches found in non-Kerr objects.
Furthermore, when an inspiral is traversing an internal island that is
surrounded by a thick chaotic layer, a new type of simultaneous multifrequency
glitch occurs that may be detectable with space interferometers such as LISA,
and can serve as evidence of an extreme-mass-ratio inspiral around a
supermassive boson star.

##  I Introduction

Ground-based gravitational-wave (GW) interferometers are about to start a new
observation run and will continue detecting GW signals from the coalescence of
compact binaries Abbott  _et al._ (2021a, b, c), and possibly from other GW
sources in the years to come. Till date, almost a hundred of compact binary
mergers have been reported. Despite the fact that the majority of events are
well-understood as either black-hole (BH) or neutron-star (NS) or mixed BH-NS
mergers, some puzzling “mass-gap” events, such as GW190814 Abbott  _et al._
(2020a) and GW190521 Abbott  _et al._ (2020b, c), challenge standard formation
scenarios and have motivated exotic alternatives (see, e.g., Calderón Bustillo
_et al._ (2021)). With the GW event catalog ever increasing, the possibility
of detecting exotic compact objects (ECOs) other than BHs and NSs is worth
exploring Cardoso and Pani (2019); Maggio  _et al._ (2021).

We are now convinced that an important amount of non-luminous exotic matter,
known as dark matter Bertone  _et al._ (2005); Clowe  _et al._ (2006);
Bergstrom (2009), is paramount in the formation Primack (1997) and
amalgamation of galaxies Corbelli and Salucci (2000), as well as in
determining the earliest, current, and future state of the universe Kolb and
Turner (1990); Primack (1997); Carroll (2001); Del Popolo (2007); Arbey and
Mahmoudi (2021); Perivolaropoulos and Skara (2022). This raises the
possibility of new fundamental particles comprising the missing cosmological
mass. With scalar fields predominantly used to model early universe physics,
the case arises that such fields could form equilibrium condensates, held
together by their own gravity, through a mechanism known as gravitational
cooling Seidel and Suen (1994); Sanchis-Gual  _et al._ (2019). Such
prototypical class of ECOs has been dubbed boson stars (BSs), and their
conceptualization dates back to the late 1960s Bonazzola and Pacini (1966);
Feinblum and McKinley (1968); Kaup (1968); Ruffini and Bonazzola (1969) and
1980s Colpi  _et al._ (1986); van der Bij and Gleiser (1987); Frieman  _et
al._ (1988, 1989); Ferrell and Gleiser (1989); Grasso (1990); Tkachev (1991);
Liddle and Madsen (1992).

At the fundamental level, BSs are the simplest localized configurations of a
complex scalar field, governed by classical equations, thus even if they have
not yet been observed in nature, they still can serve as models for compact
objects ranging from particles to stars and less dense galactic halos. In all
these cases, BSs are endowed with a balance between the dispersive nature of
scalar matter and the gravitational pull holding them together (see Jetzer
(1992); Schunck and Mielke (2003); Liebling and Palenzuela (2023) for reviews
on various types of BSs). BSs do not have an event horizon nor a singularity,
are asymptotically flat, and may exhibit stable light rings, isolated
ergoregions, and super-extremal spins, which can lead to new GW phenomenology
Liebling and Palenzuela (2023); Kesden  _et al._ (2005); Macedo  _et al._
(2013a, b, 2016); Cardoso  _et al._ (2017); Guo  _et al._ (2019). As a
consequence, they are considered among the best models of ECOs and a proxy to
test the nature of compact objects in the extreme-gravity regime Cardoso and
Pani (2019); Berti and Cardoso (2006); Pacilio  _et al._ (2020); Vaglio  _et
al._ (2023). Another important motivation concerns the fact that BSs (and
their real-scalar counterparts, known as oscillons Seidel and Suen (1991))
could form in the early universe Coleman (1985); Kusenko and Shaposhnikov
(1998) from large overdensities, being thus (meta)stable relics of inflation
Lozanov and Amin (2019); Aurrekoetxea  _et al._ (2023) and compelling dark-
matter candidates Lozanov  _et al._ (2023); Franciolini and Pani (2023).

Static BS configurations are linearly stable to perturbations Gleiser (1988);
Gleiser and Watkins (1989); Kusmartsev  _et al._ (1991); Ryan (1997a); Guzman
(2004); Macedo  _et al._ (2013a, 2016) and can consistently form dynamically
from diffuse initial states Liddle and Madsen (1992); Levkov  _et al._ (2018);
Veltmaat  _et al._ (2018); Amin and Mocz (2019); Arvanitaki  _et al._ (2020),
as well as from collisions and binary mergers Bernal and Siddhartha Guzman
(2006); Palenzuela  _et al._ (2007, 2008); Choptuik and Pretorius (2010);
Bezares  _et al._ (2017); Palenzuela  _et al._ (2017); Bezares and Palenzuela
(2018); Calderón Bustillo  _et al._ (2021); Helfer  _et al._ (2022); Bezares
_et al._ (2022); Croft  _et al._ (2023); Evstafyeva  _et al._ (2022);
Siemonsen and East (2023). On the other hand, spinning BSs were found to be
unstable toward a bar-mode instability Sanchis-Gual  _et al._ (2019); Di
Giovanni  _et al._ (2020), unless scalar self-interaction terms in the model
are sufficiently strong Siemonsen and East (2021). This instability is not
present for rotating BSs made by a massive vector field Di Giovanni  _et al._
(2020).

Till date, a variety of bosonic potentials have been considered Colpi  _et
al._ (1986); van der Bij and Gleiser (1987); Schunck and Torres (2000);
Friedberg  _et al._ (1987); Guerra  _et al._ (2019); Collodel  _et al._
(2017); Bošković and Barausse (2022); Collodel and Doneva (2022) each one
providing a different relation between the BS maximum mass and the underlying
field-theory parameters. As a rule of thumb, strong self-interactions can make
the maximum mass parametrically larger than in the free-scalar case,
motivating alternative models for supermassive objects in galactic cores,
which can mimic the shadows Vincent  _et al._ (2016); Cunha  _et al._ (2015,
2016); Cunha and Herdeiro (2018); Rosa and Rubiera-Garcia (2022) and particle
orbits around ordinary BHs Torres  _et al._ (2000); Macedo  _et al._ (2013a).

Hence, it is of utmost importance to devise further tests in order to
distinguish BSs through orbital dynamics and GW observations. In this work we
focus on extreme-mass-ratio inspirals (EMRIs), consisting of a primary
supermassive compact object (which we assume to be a self-interacting,
compact, spinning BS) and a stellar-mass compact secondary. So far, geodesic
and inspiral studies have been carried out at the equatorial plane of non-
rotating Kesden  _et al._ (2005); Diemer  _et al._ (2013); Macedo  _et al._
(2013a, b); Brihaye  _et al._ (2014) and rotating BSs Grandclement  _et al._
(2014); Grould  _et al._ (2017); Zhang  _et al._ (2022) (see also the recent
Delgado  _et al._ (2023) for circular, equatorial EMRIs around hairy BHs
Herdeiro and Radu (2014) interpolating between a Kerr BH and a BS without
self-interactions). The emitted GWs during circular equatorial geodesics and
EMRIs were also analyzed in some BS models Kesden  _et al._ (2005); Macedo
_et al._ (2013a); Collodel  _et al._ (2022); Zhang  _et al._ (2021); Delgado
_et al._ (2023). However, a geodesic and inspiral analysis is lacking for
generic, non-equatorial and non-circular, trajectories around spinning BSs.
Generic orbits will not only help to further constrain the existence of
supermassive BSs through their considerably more intricate waveform signal,
but are also tools to examine the integrability of the underlying geodesic
equations that govern particle motion. If geodesics around BSs are integrable,
then the evolution is purely deterministic and no chaotic phenomena are
present, as in the Kerr case. In this case, geodesics and EMRIs between
rotating BSs and Kerr BHs will differ only by their different multipolar
structure Ryan (1995, 1997b, 1997c); Pacilio  _et al._ (2020); Vaglio  _et
al._ (2022, 2023). Alternatively, if geodesics around rotating BSs break
integrability then direct and indirect chaotic imprints emerge that are
clearly distinguishable both at the orbital Apostolatos  _et al._ (2009);
Lukes-Gerakopoulos  _et al._ (2010); Contopoulos  _et al._ (2011); Lukes-
Gerakopoulos (2012); Zelenka and Lukes-Gerakopoulos (2017); Zelenka  _et al._
(2020); Lukes-Gerakopoulos and Witzany (2021); Destounis  _et al._ (2020);
Deich  _et al._ (2022); Chen  _et al._ (2022) and GW level Destounis  _et al._
(2021); Destounis and Kokkotas (2021); Destounis  _et al._ (2023a).

It is the main goal of this work to study generic orbits and EMRIs around
supermassive self-interacting BSs in order to assess if their exotic
multipolar structure Ryan (1995, 1997b, 1997c); Pacilio  _et al._ (2020);
Vaglio  _et al._ (2022, 2023) breaks integrability and gives rise to imprints
of chaos. If geodesics around such objects are nonintegrable, then GW
observations from future spaceborne detectors like the Laser Interferometer
Space Antenna (LISA) Glampedakis (2005); Amaro-Seoane  _et al._ (2017);
Baibhav  _et al._ (2021); Amaro-Seoane  _et al._ (2022); Arun  _et al._
(2022); Karnesis  _et al._ (2022), TianQin Luo  _et al._ (2016) and Taiji Ruan
_et al._ (2020a, b) may lead to distinguishable effects that can break the
degeneracy between supermassive BHs and BSs in galactic centers111Note that
also EMRIs around stellar-mass compact objects would be potentially detectable
by third-generation detectors such as the Einstein Telescope if the secondary
is a subsolar compact object Barsanti  _et al._ (2022).. In what follows, we
adopt geometrized units so that $G=c=1$.

##  II Geodesics and Chaos

Generic stationary and axisymmetric spacetimes can be written as

$ds^{2}=g_{tt}dt^{2}+2g_{t\varphi}dtd\varphi+g_{rr}dr^{2}+g_{\theta\theta}d\theta^{2}+g_{\varphi\varphi}d\varphi^{2},$
(1)

where the metric tensor components are, in general, functions of $r$ and
$\theta$, and the coordinate system $\left(t,r,\theta,\varphi\right)$ can be,
e.g., of Boyer-Lindquist type Boyer and Lindquist (1967) or quasi isotropic
Misner  _et al._ (1973). The motion of a secondary point-particle orbiting the
spacetime geometry of the primary compact object (as defined by Eq. (1)) is
described by the geodesic equations,

$\ddot{x}^{\kappa}+\Gamma^{\kappa}_{\lambda\nu}\dot{x}^{\lambda}\dot{x}^{\nu}=0,$
(2)

where $\Gamma^{\kappa}_{\lambda\nu}$ are the Christoffel symbols of spacetime,
$x^{\kappa}$ is the four-position vector of the orbit, and the overdots denote
differentiation with respect to proper time $\tau$. The geodesic equation (2)
breaks down into four equations of motion, one for each coordinate component
of the particle’s trajectory, $x^{\kappa}(\tau)$. This system of second order,
coupled differential equations can be considerably simplified using spacetime
symmetries. The spacetime (1) assumed in this work possesses two Killing
vector fields resulting from stationarity and axisymmetry (the metric tensor
is $t$\- and $\varphi$-independent), yielding two conserved quantities
throughout geodesic motion, i.e. the specific energy and $z$-component of the
angular momentum of the particle,

$E/m=g_{tt}\dot{t}+g_{t\varphi}\dot{\varphi},\,\,\,\,\,\,\,\,\,L_{z}/m=g_{t\varphi}\dot{t}+g_{\varphi\varphi}\dot{\varphi},$
(3)

where $m$ is the mass of the orbiting particle. Rearranging Eq. (3) we obtain
two first-order, decoupled differential equations for the $t$ and $\varphi$
momenta as

$\dot{t}=\frac{Eg_{\varphi\varphi}+L_{z}g_{t\varphi}}{m\left(g^{2}_{t\varphi}-g_{tt}g_{\varphi\varphi}\right)},\,\,\,\,\,\,\,\,\,\dot{\varphi}=\frac{Eg_{t\varphi}+L_{z}g_{tt}}{m\left(g_{tt}g_{\varphi\varphi}-g^{2}_{t\varphi}\right)}\,,$
(4)

which can be solved once $r(\tau)$ and $\theta(\tau)$ are known. The two
remaining equations of motion for $r$ and $\theta$ are, in general, coupled
and of second differential order. Test particles in geodesic motion provide a
third constant of motion, namely the conservation of their rest mass (or
equivalently their four-velocity
$g_{\lambda\nu}\dot{x}^{\lambda}\dot{x}^{\nu}=-1$) which leads to a constraint
equation of the form

$\dot{r}^{2}+\frac{g_{\theta\theta}}{g_{rr}}\dot{\theta}^{2}+V_{\text{eff}}=0,$
(5)

with $V_{\text{eff}}$ being a Newtonian-like potential,

$V_{\text{eff}}\equiv\frac{1}{g_{rr}}\left(1+\frac{g_{\varphi\varphi}E^{2}+g_{tt}L_{z}^{2}+2g_{t\varphi}EL_{z}}{m^{2}\left(g_{tt}g_{\varphi\varphi}-g_{t\varphi}^{2}\right)}\right),$
(6)

that characterizes bound geodesic motion. When $V_{\text{eff}}=0$ the
resulting curve is called curve of zero velocity (CZV) since
$\dot{r}=\dot{\theta}=0$ there.

If a hypothetical rank-two (or higher-rank) Killing tensor field exists, then
the motion of $r(\tau)$ and $\theta(\tau)$ decouples into first-order
differential equations. A special case where a rank-two Killing tensor exists
is the Kerr solution where the role of the separation constant is played by
the famous Carter constant Carter (1968). However, since the compact object we
will assume throughout our analysis has a different multipolar structure than
that of a Kerr BH (and furthermore it is known only numerically), it is
unlikely that a separation, Carter-like, constant (or any other higher-rank
Killing tensor) exists. Thus, to evolve orbits we will use the coupled second-
order differential equation system for $r$ and $\theta$, together with Eqs.
(4) and (5), without any further symmetry assumptions besides stationarity and
axisymmetry.

The existence of Carter’s separation constant is not only a useful tool to
evolve trajectories in Kerr spacetime faster, but rather implies an important
aspect of geodesics, namely their integrability. Kerr spacetimes have four
degrees of freedom. Stationarity and axisymmetry leads to the reduction of
degrees of freedom to two. Taking into account the existence of the Carter
constant and the conservation of the rest mass of the test particle reduces
the degrees of freedom of orbital motion to zero; geodesics around Kerr BHs
are integrable, deterministic and do not present chaotic features Contopoulos
(2003). Unfortunately, Carter’s symmetry is extremely fragile and in many
occasions is broken by simply deforming the multipolar structure of spacetime
by considering accretion disks and BH environments Zipoy (1966); Voorhees
(1970); Polcar  _et al._ (2022), modifications of gravity Destounis  _et al._
(2020); Deich  _et al._ (2022), neutron stars Hartle and Thorne (1968); Manko
_et al._ (2000a); Pappas and Apostolatos (2013); Pappas (2017), exotic compact
objects Manko and Novikov (1992); Manko  _et al._ (2000b), or in a
parameterized way by introducing agnostic deformations to Kerr; a class of
metrics called bumpy or non-Kerr BHs Collins and Hughes (2004); Glampedakis
and Babak (2006); Vigeland and Hughes (2010); Johannsen and Psaltis (2011);
Vigeland  _et al._ (2011); Johannsen (2013a, b); Emparan  _et al._ (2014);
Cardoso  _et al._ (2014); Rezzolla and Zhidenko (2014); Konoplya  _et al._
(2016); Moore  _et al._ (2017). In these cases, the integrability property may
be broken, leading to chaotic effects Contopoulos (2003). Full-blown, ergodic
chaos is not expected to occur in astrophysical scenarios, such as EMRIs, but
nonintegrability-like effects are anticipated even when geodesics are
integrable Flanagan and Hinderer (2012); Berry  _et al._ (2016); Speri and
Gair (2021); Gupta  _et al._ (2022) due to dissipation. This is mainly due to
the manifestation occurring around transient orbital resonances which are
expected to affect EMRI evolution and parameter estimation, though
nonintegrable EMRIs are even more likely to amplify these effects and can
introduce clear chaotic phenomena Apostolatos  _et al._ (2009); Lukes-
Gerakopoulos  _et al._ (2010); Contopoulos  _et al._ (2011); Lukes-
Gerakopoulos (2012, 2012); Zelenka and Lukes-Gerakopoulos (2017); Zelenka  _et
al._ (2020); Destounis  _et al._ (2021); Destounis and Kokkotas (2021); Lukes-
Gerakopoulos and Witzany (2021); Destounis  _et al._ (2023a). A generic orbit
of an integrable system can be described by its revolution frequency
$\omega_{\varphi}$ and two libration-like frequencies; the frequency of
oscillation from the periapsis to the apoapsis and back, $\omega_{r}$, and the
oscillation frequency through the equatorial plane, $\omega_{\theta}$. Generic
orbits possess irrational ratios of the above frequencies and these orbits
fill densely the available phase space of a three-dimensional torus. On the
other hand, resonant (periodic) orbits, have commensurate (rational) ratios of
orbital frequencies which means that the particular orbits are returning to
their initial position after some revolutions depending on their periodicity,
thus they are not phase-space filling. When orbital resonances are encountered
during inspirals of equal mass binaries, e.g. in the kHz band of
LIGO/Virgo/KAGRA detectors, they do not affect the evolution since the
inspiral is extremely rapid. However, when the binary is asymmetric, i.e. an
EMRI, then the adiabatic nature of the secondary’s motion can experience
orbital resonances for a significant number of cycles and lead to substantial
effects, such as cumulative dephasing and putative erroneous parameter
estimation Flanagan and Hinderer (2012); Berry  _et al._ (2016); Speri and
Gair (2021).

When a nonintegrable perturbation is introduced to the system, two theorems,
namely the Kolmogorov-Arnold-Moser (KAM) Möser (1962); Arnol’d (1963) and
Poincaré-Birkhoff Birkhoff (1913) theorems, dictate how the phase space
structure is altered around resonant points. KAM theorem ensures that when the
orbits are sufficiently away from resonances, the system behaves as if it were
integrable. The trajectories in the phase space lie on a torus defined by the
integrals of motion and successive intersections of the orbits on a
perpendicular 2-dimensional surface (the Poincaré surface of section) form
curves that organize around a common, fixed, central point. These are called
KAM curves and the whole structure is known as Poincaré map, whose central
point corresponds to a planar circular orbit. Close to periodic orbits, the
KAM curves disintegrate into two sets of periodic points in the Poincaré map,
in accord to Poincaré-Birkhoff theorem; the stable ones which are surrounded
by islands of stability (resonant islands) and the unstable ones where chaotic
orbits emanate and surround the islands of stability with thin layers. The
whole structure around resonances of nonintegrable systems is called a
Birkhoff chain (see Fig. 2 in Lukes-Gerakopoulos  _et al._ (2010) for an
illustration). The crucial aspect of resonant islands, and their significance
in EMRIs, is the fact that the rational ratio of the orbital frequencies
$\omega_{r}/\omega_{\theta}$ is shared throughout the island for all geodesics
that occupy it. Put in other words, integrable EMRIs experience resonances
that occupy a zero-volume point in phase space while nonintegrable EMRIs
exhibit prolonged resonances where the secondary is locked in perfect
resonance for a rather significant amount of revolutions ($\sim 200-300$
cycles depending on the non-Kerr object Destounis  _et al._ (2021); Destounis
and Kokkotas (2021); Lukes-Gerakopoulos  _et al._ (2010); Lukes-Gerakopoulos
and Witzany (2021); Destounis  _et al._ (2023a)) without taking into account
pre- and post-resonant effects Ruangsri and Hughes (2014) and conservative
effects from gravitational self-force Barack (2009) which are essential for
integrable EMRIs to show signs of transient resonant phenomena.

The existence of ergodic chaos and islands of stability around periodic points
are direct and indirect signatures of nonintegrability, respectively,
therefore sketching a detailed Poincaré map can unravel if the system under
study is deterministic. When the perturbation introduces a slight
nonintegrable deformation to the system there are other techniques to assess
if and at which initial conditions of geodesics the aforementioned phenomena
manifest themselves. The rotation number is one of the most helpful tool to
find regions of interest in order to search for islands of stability. We
calculate it by tracking the angle $\vartheta$ between successive
intersections on KAM curves, relative to the fixed central point of the
Poincaré map. The rotation number is defined as the summation of all angles
$\vartheta$ measured between consecutive intersections, i.e. Contopoulos
(2003)

$\nu_{\vartheta}=\frac{1}{2\pi N}\sum_{i=1}^{N}\vartheta_{i},$ (7)

where $N$ is the number of angles measured. When $N\rightarrow\infty$, Eq. (7)
asymptotes to the orbital frequency ratio
$\nu_{\vartheta}=\omega_{r}/\omega_{\theta}$. Calculating consecutive rotation
numbers for different initial conditions of orbits, by smoothly varying one of
the parameters of the system while keeping the rest fixed, leads to a rotation
curve.

Integrable systems show monotonous rotation curves, while nonintegrable
systems display discontinuities in the monotonicity through the formation of
transient plateaus with a non-zero width when geodesics transverse resonant
islands. Inflection points can also appear when trajectories pass through
unstable periodic points. Nevertheless, by changing the initial conditions
properly, the orbit can be driven through the island and give rise to a
plateau. So far, all studies have dealt with compact objects that are either
non-Kerr or bumpy in nature and either possess an event horizon or have
serious causal structure pathologies. In the following sections, we will
examine the characteristic features of rotating self-interacting BSs with
geodesics and approximate EMRI evolutions in order to first assess whether
geodesic motion in such spacetimes is nonintegrable and, then, to establish
various phenomenological imprints of chaos in the associated GW signal.

##  III Self-interacting rotating BSs

The equilibrium configurations of rotating BSs can be constructed starting
from the globally $U(1)$-invariant action

${\cal S}=\int d^{4}x\sqrt{-g}\left[\frac{R}{16\pi}-{\cal L}_{\phi}\right]\ ,$
(8)

which describes the dynamics of a complex, massive scalar field $\phi$,
minimally coupled to gravity. The Lagrangian $\mathcal{L}_{\phi}$ considered
in this work is characterized by a mass parameter $\mu$, plus quartic
repulsive corrections, controlled by the coupling $\lambda$, which add
linearly to the kinetic term as

${\cal
L}_{\phi}=-\frac{1}{2}g^{\alpha\beta}\phi^{*}_{,\alpha}\phi_{,\beta}-\frac{1}{2}\mu^{2}|\phi|^{2}-\frac{1}{4}\lambda|\phi|^{4}\
.$ (9)

The requirement of stationarity and axisymmetry, leads to the following ansatz
for the scalar field

$\phi=\phi_{0}(r,\theta)e^{i(s\varphi-\Omega t)}\ ,$ (10)

where $\Omega>0$ is the field’s angular frequency, which determines the phase
evolution in time, while $s$ is an integer called the azimuthal (or
rotational) winding number and can be proven to correspond to the ratio
between the conserved angular momentum and particle number Schunck and Mielke
(1996); Yoshida and Eriguchi (1997); Schunck and Mielke (2003). The ansatz
(10) ensures that the stress-energy tensor
$T_{\alpha\beta}[g^{\alpha\beta},\phi,\partial_{\alpha}\phi]$ does not depend
on $t$ and $\varphi$ and sources a spacetime metric with time and azimuthal
symmetry.

In this work, following Ryan (1995); Vaglio  _et al._ (2022), we consider
spinning BSs with large self-interactions, i.e. characterized by
$\lambda/\mu^{2}\gg 1$. In this limit the $(r,\theta)$ scalar profile is
approximately constant in the star’s interior and one can neglect the radial
and polar derivatives of the field ($\partial_{r}\phi\sim 0$,
$\partial_{\theta}\phi\sim 0$), while assuming $\phi\sim 0$ in the exterior.
This allows expressing the stress-energy tensor of the BS as that of a perfect
fluid and to define the radius as in ordinary stars222In the general case, the
stress-energy tensor of a BS contains anisotropic terms and the scalar field,
although exponentially decreasing, extends up to infinity, so that the radius
is conventionally defined as the value of the radial coordinate enclosing a
sufficiently large fraction (typically $99\%$) of the total mass Kaup (1968);
Ruffini and Bonazzola (1969); Amaro-Seoane  _et al._ (2010). In our case,
instead, the scalar field has compact support and the radius is uniquely
defined.. Furthermore, with an appropriate redefinition of the variables and
parameters, the coupling constants in the Lagrangian can be factored out, so
that each numerical solution corresponds to a one-parameter family of
configurations, controlled by the effective mass parameter
$M_{B}=\sqrt{\lambda}/\mu^{2}$ Vaglio  _et al._ (2022). The dimensionful
physical quantities characterizing each BS, such as its mass, radius and
energy-density, can be obtained from the dimensionless ones, characterizing
each numerical solution, by multiplying them by the required power of $M_{B}$
to match their dimension in mass and scale correspondingly as it changes.

Adopting quasi-isotropic coordinates, the line element (1) can be described
through four independent functions $(\rho,\gamma,\omega,\alpha)$ of
$(r,\theta)$ as

$\displaystyle g_{tt}(r,\theta)$
$\displaystyle=-e^{\gamma(r,\theta)+\rho(r,\theta)}$
$\displaystyle+\omega(r,\theta)^{2}e^{\gamma(r,\theta)-\rho(r,\theta)}r^{2}\sin^{2}\theta\,,$
(11a) $\displaystyle g_{rr}(r,\theta)$
$\displaystyle=\frac{g_{\theta\theta}(r,\theta)}{r^{2}}=e^{2\alpha(r,\theta)}\,,$
(11b) $\displaystyle g_{t\varphi}(r,\theta)$
$\displaystyle=-\omega(r,\theta)e^{\gamma(r,\theta)-\rho(r,\theta)}r^{2}\sin^{2}\theta\,,$
(11c) $\displaystyle g_{\varphi\varphi}(r,\theta)$
$\displaystyle=e^{\gamma(r,\theta)-\rho(r,\theta)}r^{2}\sin^{2}\theta\,.$
(11d)

We constructed the solutions numerically following the method described in
Ryan (1997a); Vaglio  _et al._ (2022) which makes use of an integral
representation of the Einstein equations to set up an iterative integration
scheme. The algorithm starts with a solution corresponding to a non-rotating
BS defined on a two-dimensional grid of $r\in(10^{-6},10)M_{B}$ (which
typically corresponds to $r\in(10^{-5},100)M$ in terms of the BS mass $M$),
$\theta\in(0,\pi/2)$, and converges to a spinning configuration, close to the
initial one, within roughly $150$ iterations, with a relative error of
$\mathcal{O}(10^{-4})\%$ for all metric functions and the scalar profile. The
details of the numerical implementation are provided in Vaglio  _et al._
(2022).

The maximum mass that can be reached by these stars, scales as Colpi  _et al._
(1986); Vaglio  _et al._ (2022)

$M_{\rm
max}\sim\gamma(\chi)M_{B}=\gamma(\chi)\frac{\sqrt{\lambda\hbar}}{m_{s}^{2}}M_{p}^{3}\
,$ (12)

where $m_{s}=\mu\hbar$ is the mass of the boson, $M_{p}$ the Planck mass and
$\gamma(\chi)$ is a $\mathcal{O}(0.1)$ factor which depends on the
dimensionless spin $\chi=J/M^{2}$ (where $J$ is the angular momentum of the
solution). This means, for instance, that, for
$\lambda\sim\mathcal{O}(\hbar^{-1})$ and $m_{s}$ in the range $0.1$–$100$ MeV,
the model allows for compact stellar configurations with $M_{\rm max}$ in the
range $10$–$10^{7}\leavevmode\nobreak\ M_{\odot}$.

In what follows, we consider rapidly-rotating supermassive configurations with
$\chi\sim 0.8$, same mass and decreasing compactness, whose properties are
listed in Table 1. All the configurations have topological genus 1 Vaglio  _et
al._ (2022), at variance with their vector counterparts, spinning Proca stars
Brito  _et al._ (2016), which have instead a spherical topology even when
spinning.

|  $M/M_{\odot}$ |  $\,\,\chi$ |  $\,\,\,\,\mathcal{C}$ |  $\,\,\,\Omega$ |  $M_{B}/M_{\odot}$
---|---|---|---|---|---
Case 1 |  $10^{6}$ |  $0.8$ |  $0.25$ |  $0.73$ |  $8.5\times 10^{6}$
Case 2 |  $10^{6}$ |  $0.8$ |  $0.19$ |  $0.80$ |  $9.0\times 10^{6}$
Case 3 |  $10^{6}$ |  $0.8$ |  $0.16$ |  $0.82$ |  $1.0\times 10^{7}$
Table 1: Self-interacting, rotating BSs considered in this work. The three
configurations have the same dimensionless spin and mass but different
compactnesses and frequencies (due to the different effective coupling
$M_{B}$). The compactness is defined by ${\cal C}\equiv M/R$, where $R$ is the
star’s perimeteral radius Herdeiro and Radu (2015). The order of magnitude
chosen for $M_{B}$ corresponds, considering
$\lambda/\mu^{2}\sim\mathcal{O}(100)\gg 1$, to a mass of the boson
$m_{s}=\hbar\mu\sim 10^{-14}\,{\rm eV}$.

The configuration of Case 1 is the most compact one and corresponds to the
maximum-mass solution for $\chi=0.8$. A comparison with the Kerr solution,
sharing the same mass and spin, is shown in Fig. 1, in terms of the functions
$U_{tt}\equiv(g_{tt}-1)/2$ and $U_{rr}\equiv(g_{rr}+1)/2$ evaluated on the
equatorial plane. The corresponding Kerr event horizon and BS radius are also
shown, as well as the Newtonian gravitational potential $M/r$ for an object of
the same mass.

Figure 1: Comparison between the Newtonian potential $M/r$ (black dashed) and
the equatorial radial profiles of $U_{tt}\equiv(g_{tt}-1)/2$ and
$U_{rr}\equiv(g_{rr}+1)/2$ for the Kerr metric (red dashed), and the rotating
BS (blue solid) corresponding to Case 1. The black solid and gray dashed
vertical lines correspond to the location of the Kerr horizon and the BS
radius in the quasi-isotropic radial variable $r/M$.

All the configurations in Table 1 are in a linearly perturbative stable branch
of the mass-frequency diagram. It is known that BSs with no or weak self-
interactions are subject to a dynamical non-axisymmetric instability that
develops on short timescales Sanchis-Gual  _et al._ (2019). Such instability
has been shown to be quenched if the scalar self-interactions are sufficiently
strong Siemonsen and East (2021). This is precisely the limit in which our
solutions are obtained, with the explicit value of $\lambda/\mu^{2}$ depending
on the individual choices for $\lambda$ and $\mu$, once $M_{B}$ is fixed.
Rotating BSs in this regime have also recently been formed dynamically in
numerical simulations as a result of a binary coalescence, starting from non-
rotating components in a quasi-circular orbit Siemonsen and East (2023).

Another source of instability is potentially linked to the presence of light
rings. It has been recently shown Cunha  _et al._ (2023) that exotic
ultracompact objects, i.e. compact objects featuring a light ring Cardoso and
Pani (2019), are unstable under nonlinear perturbations due to the presence of
a stable photon orbit in their interior, which traps massless perturbation
modes Cardoso  _et al._ (2014). Light rings are found as stationary points of
the effective potentials

$V_{\pm}\propto\frac{-g_{t\varphi}\pm\sqrt{g_{t\varphi}^{2}-g_{tt}g_{\varphi\varphi}}}{g_{\varphi\varphi}},$
(13)

where the plus (minus) sign corresponds to orbits that are co-rotating
(counter-rotating) with the star. Among our three configurations, only Case 1
exhibits a pair of light rings (one stable and one unstable) in the effective
potential for counter-rotating orbits $V_{-}$, as depicted in Fig. 2. Although
this likely makes this configuration prone to instability– either through
migration to a stable non-ultracompact configuration or collapse into a BH –
the timescale of such instability is unknown due to the absence of simulations
for these particular BS models. Furthermore, the frequency of the star is
close to the critical frequency for which no light rings are featured,
corresponding to an infinite timescale. For these reasons, we chose to keep
this solution as an example of high-compactness spinning BS. In any case,
stable configurations near Case 1 are expected to have similar geodesic
properties.

Figure 2: Equatorial effective potentials $V_{-}(r)$ for the BS configurations
corresponding to our three case studies. The BS in Case 1 has two light rings,
corresponding to the stationary points $V^{\prime}_{-}(r)=0$, one of which is
unstable and outside the star, while the other is stable and inside it.
Stationary points of $V_{-}(r)$ correspond to counter-rotating circular photon
orbits, while $V_{+}$ (not shown here) is associated to co-rotating photon
orbits which are not present for these configurations. Figure 3: Vertical
section showing the energy density distribution $\epsilon$ of our
representative BSs (from left to right with decreasing compactness: Case 1, 2,
3), normalized to the largest value $\epsilon_{\text{max}}$ (attained in Case
1).

For completeness, in Fig. 3, we present the vertical section of the energy
density of the rotating BS cases considered, where the relative compactness
and star radii are evident. All energy densities have been normalized with
respect to the maximum value reached by Case 1.

##  IV Geodesic analysis of generic orbits

Before embarking into the discussion of the results obtained through geodesic
evolutions of orbits around rotating BSs, we report that the metric tensor
components of our three configurations have been constructed with varying
number of grid points on the $r$ and $\theta$ coordinates to test convergence.
The islands of stability and corresponding plateau widths at resonance in
rotation curves reported below converge to $\lesssim 1\%$ as we increased the
number of grid points from $300\times 200$ to $1000\times 700$. Interestingly,
the island and plateau widths converge from below, meaning that more grid
points lead to slightly larger islands till convergence occurs, which
guarantees their existence regardless of the grid resolution, and therefore
ensures nonintegrability. For all simulations presented herein, we find that
the constraint equation (5) is satisfied to within one part in $10^{8}$ for
the first 5 to 10 thousand revolutions (and intersections through the
equatorial plane, depending on the initial conditions and BS configuration).

Finally, in all cases we fix the mass ratio of the EMRI to $m/M=10^{-6}$ and
choose the initial energy and $z$-component of the angular momentum as
$E/m=0.95$ and $L_{z}/m=3M$, respectively.

Figure 4: Left: Poincaré map of a secondary with $m=1M_{\odot}$ orbiting
around a compact rotating BS with $M=10^{6}M_{\odot}$, $\mathcal{C}=0.25$
while the rest of the configuration quantities are stated in Table 1 (Case 1).
The secondary’s conserved energy and angular momentum are $E/m=0.95$,
$L_{z}/m=3M$, respectively. The fixed initial conditions chosen here are
$\dot{r}(0)=0,\,\theta(0)=\pi/2$, and $\dot{\theta}(0)$ is defined by the
constraint equation to guarantee bound motion, while $r(0)$ is varied. Black
curves that surround the central fixed point of the map designate
intersections of generic orbits through the equatorial plane, with different
initial $r(0)$, while colored curves designate intersections that belong to
different resonant islands of stability. Right: Zoom into the leftmost region
where the $1/2$ island of stability resides. Similar encapsulated structure is
found for the rest of the islands. Figure 5: Left: Rotation curve
corresponding to the same system as in Fig. 4. The vertical dashed line
represents the radius of the BS while the inset zooms in the surrounding
region of the $1/2$ plateau. Right: Same as in the left panel but with
$\dot{r}(0)=0.1$. The initial velocity gives access to the $2/3$ island of
stability.

###  IV.1 Case 1

Case 1 is a representative example for a compact spinning BS. The Poincaré map
shown in Fig. 4 clearly shows that rotating BSs with large compactness are
nonintegrable, since around resonances such as $\nu_{\vartheta}=1/2,\,2/3$
islands of stability form that encapsulate stable periodic points of
geodesics.

From the maps, we calculate the rotation curves for two different choices of
initial radial velocity, one with $\dot{r}(0)=0$ and one with
$\dot{r}(0)=0.1$, shown in Fig. 5. The rotation curves for this configuration
looks quite similar to those found in Lukes-Gerakopoulos  _et al._ (2010);
Destounis  _et al._ (2020, 2023a), i.e. an inflection point around the $2/3$
resonance and a plateau at $1/2$ resonance when $\dot{r}(0)=0$. The absence of
an event horizon allows for bound geodesics even inside the BS, with its
radius designated with a vertical dashed line, where surprisingly the $1/2$
plateau resides. Therefore, the first non-BH feature is the existence of
resonant islands even inside the star. Another interesting feature is the
change in monotonicity close to plunge333We refer to plunge for the orbits
that exit the CZV (the separatrix). In the BH case, these orbits end up
plunging into the horizon, while in BSs there might be cases in which plunging
orbits escape to infinity. for Case 1, which has also been observed for Kerr
BHs with soft scalar hair and BSs in the study of epicyclic frequencies
Franchini  _et al._ (2017). The inflection point at rotation number $2/3$ for
$\dot{r}(0)=0$, in the left panel of Fig. 5, designates the passage of an
orbit between two edges of resonant islands, where there is an unstable
periodic point. Changing the initial radial velocity to $\dot{r}(0)=0.1$, as
we show for the rotation curve in Fig. 5 (right panel), the phase space
trajectory crosses the island and a plateau at $2/3$ appears.

Figure 6: Left: Same as in Fig. 4 but for Case 2 (intermediate compactness
$\mathcal{C}=0.19$). Right: Zoom into the leftmost region where the $4/7$
island of stability resides. Similar encapsulated structure is found for the
rest of the islands. Figure 7: Left: Rotation curve corresponding to the same
system as in Fig. 5. The horizontal dashed line represents the radius of the
BS while the inset zooms in the surrounding region of the external $4/7$
plateau. Right: Same as in the left panel but with $\dot{r}(0)=0.058$. The
initial velocity gives access to the $2/3$ island of stability.

###  IV.2 Case 2

As we slightly decrease the compactness of the BS, further non-BH effects
arise. For an intermediate compactness rotating BS (Case 2) the Poincaré map,
shown in Fig. 6, displays a much more intricate structure around subdominant
resonances. Firstly, the radius of the star increases, therefore orbits are
allowed to exist in its interior, which is precluded for BHs, and should lead
to more prominent effects. Intriguingly, due to the aforementioned phenomenon,
here we find for the first time degenerate resonant islands – namely islands
of stability occurring for two different ranges of the radial coordinate. In
particular, Fig. 6 shows the typical exterior $2/3$ island as well as two sets
of Birkhoff chains for the $4/7$ resonance, one in the exterior and one in the
interior of the BS.

The rotation curve in Fig. 7 for $\dot{r}(0)=0$ (left panel) confirms all the
above. Decreasing the compactness leads to a more prominent change in
monotonicity when the secondary enters the star. This in turn allows for
degenerate plateaus both in the interior and the exterior of the star’s
geodesics, with the interior plateau being wider than the exterior one. Again,
choosing a different initial velocity, we can turn the inflection point of
$2/3$ resonance into a plateau, since the velocity pushes the orbit to
traverse the island.

Figure 8: Top left: Same as in Figs. 4 and 5 but for Case 3 (smaller
compactness, $\mathcal{C}=0.16$). Top right: Zoom into the region where the
leftmost $5/8$ island of stability resides. Bottom left: Zoom into the region
where the leftmost $4/7$ island of stability resides. Bottom right: Zoom into
a region of the Poincaré map on the top left where a chaotic layer is present,
shown with scattered black points that emanate from the unstable $2/3$
periodic point with $\dot{r}(0)=0$, and parts of two islands of stability,
namely the $2/3$ and $13/20$. Figure 9: Top left: Rotation curve corresponding
to the same system as in Fig. 8. The vertical dashed line represents the
radius of the BS while the inset zooms in the surrounding region of the
external $4/7$ plateau. Top right: Zoomed part of the rotation curve on the
top left figure on the subdominant resonances $\nu_{\vartheta}=5/8$ and $4/7$
in descending order. Bottom left: Zoomed part of the rotation curve on top
left figure at the chaotic layer of the unstable periodic point for $2/3$
resonance. A further zoom is shown in the inset of two extremely subdominant
islands with $\nu_{\vartheta}=15/23$ and $13/20$ in descending order. Bottom
right: Rotation curve for the same parameters and initial conditions as in the
top left figure but with $\dot{r}(0)=0.155$. The initial velocity gives access
to the $2/3$ island of stability.

###  IV.3 Case 3

Case 3 is probably the most interesting one. Its relatively smaller
compactness allows for a plethora of new phenomena in both Poincaré maps (Fig.
8) and rotation curves (Fig. 9). The radius of the BS increases even more and
the CZV enters deep inside the star where bound geodesics are still possible.
This leads to a multitude of dominant and subdominant resonances, as shown in
the top left panel of Fig. 8, as well as further degenerate Birkhoff chains.
In the top right and bottom left panels, a zoom of the interior $5/8$ and
$4/7$ islands is presented. Probably the most interesting feature is displayed
in the bottom right panel. Two islands of stability appear, namely the $2/3$
(in red) and the interior $13/20$ part of the island (in pink). Furthermore,
we find a visible chaotic layer that surrounds the $2/3$ island, designated
with scattered black points. Their source is the unstable periodic point at
$\dot{r}(0)=0$ and $r(0)\sim 3.31M$. To the best of our knowledge, this is the
first case in which a chaotic layer appears for a motivated model of compact
object without any pathologies. This layer should give rise to significant
effects on rotation curves and eventual EMRIs crossing through the particular
$2/3$ island when the fluxes will be taken into account Destounis  _et al._
(2021); Destounis and Kokkotas (2021); Destounis  _et al._ (2023a).

Figure 9 shows the significantly modified rotation curves with respect to a
typical non-Kerr BH (see, e.g., Apostolatos  _et al._ (2009); Lukes-
Gerakopoulos  _et al._ (2010); Contopoulos  _et al._ (2011); Lukes-
Gerakopoulos (2012); Zelenka and Lukes-Gerakopoulos (2017); Destounis  _et
al._ (2020); Lukes-Gerakopoulos and Witzany (2021); Deich  _et al._ (2022);
Chen  _et al._ (2022); Destounis  _et al._ (2023a)). Degenerate $4/7$ plateaus
are visible in the interior and exterior of the star, as well as a robust
internal plateau for the $5/8$ resonance. Although finding its external
counterpart requires extreme fine-tuning, it is guaranteed that such region
exists in the exterior. Closing in to the plunge (i.e., the inner boundary of
the CZV), deep inside the star, a very large inflection point appears, which
after zooming in (see Fig. 9, bottom left) reveals a thick chaotic layer with
ill-defined rotation numbers. We also find a couple of extremely subdominant
islands for the $13/20$ (also shown in Fig. 8) and $15/23$ resonances. Even
more interestingly, by assuming an appropriate initial velocity for the
geodesics, we managed to find, for the first time, an island (which leads to a
plateau) that begins from the exterior of the star and ends in the interior.
Its width ($\sim 2.6M$ in quasi-isotropic or $\sim 3.6M$ in Boyer-Lindquist
coordinates) is so wide-spread that supersedes any other plateau ever found in
non-Kerr spacetimes, where the widest one found is of order $\sim 0.05M$ (in
Boyer-Lindquist coordinates) and barely compares with the plateaus for Case 1
($2/3$ resonance) and Case 2 ($4/7$ resonance) island widths.

For completeness, in Fig. 10 we present all rotation curves obtained with
$\dot{r}(0)=0$ and varying $r/M$, together with the rotation curve of a Kerr
BH with $\chi=0.8$. The differences between BH orbits and non-compact rotating
BSs are evident and the BSs present novel features due to the absence of an
event horizon. For fixed $E$ and $L_{z}$, we find that even the most compact
configuration considered differs dramatically from that of Kerr, while the
rest completely disengage with the typical behavior of a Kerr rotation curve.
Perhaps a rotation curve in the lines of those presented in Lukes-Gerakopoulos
_et al._ (2010); Destounis  _et al._ (2020) may seem more similar to that of
Case 1 BS but still the existence of the radius, the change in monotonicity
when the orbits enters the star, the degenerate plateaus, and the existence of
an observable chaotic layer should make rotating BSs having a contrasting and
discernible behavior of geodesics when compared to those of non-Kerr
spacetimes.

Figure 10: Combined rotation curves of a secondary with $m=1M_{\odot}$
orbiting around a rotating BS with $M=10^{6}M_{\odot}$, varying compactness
$\mathcal{C}$ and $\chi=0.8$. The secondary’s conserved energy and angular
momentum are $E/m=0.95$, $L_{z}/m=3M$, respectively. The fixed initial
conditions chosen here are $\dot{r}(0)=0,\,\theta(0)=\pi/2$ and
$\dot{\theta}(0)$ is defined by the constraint equation to guarantee bound
motion, while $r(0)$ is varied. The horizontal dashed lines represents the
radii of each BS configuration. For comparison, we include the rotation curve
of a Kerr BH with the same spacetime and geodesic parameters.

##  V Adiabatic inspiral and waveforms

Inspiral in the EMRI limit can be conveniently described within BH
perturbation theory Barack (2009); Pound and Wardell (2021). Owing to the
hierarchy of scales, the dynamics can be studied as a small perturbation of
the geodesic motion of the secondary test mass around the primary object. To
leading order in the mass ratio, one could evolve geodesic quantities
adiabatically by taking into account dissipation due to radiative degrees of
freedom. Higher order corrections require including conservative and, in
general, self-force effects during the inspiral Barack (2009); van de Meent
(2016, 2018); Pound (2017); Pound and Wardell (2021).

While this program has been extremely successful for standard EMRIs around a
Kerr BH within General Relativity, going beyond the standard paradigm is much
more challenging, even at the leading order. Indeed, dissipative corrections
are computed using the Teukolsky formalism Teukolsky (1973); Teukolsky and
Press (1974); Detweiler (1978), which allows separating the perturbations of a
Kerr BH in General Relativity and computing fluxes through numerically
integration of inhomogeneous ordinary differential equations with a point-
particle source term Sasaki and Nakamura (1982); Mino  _et al._ (1997). This
technology has been widely tested in the frequency Hughes (2000, 2001); Drasco
and Hughes (2006) and time domain Poisson (2004); Martel (2004); Lopez-Aleman
_et al._ (2003); Khanna (2004). However, it heavily relies on the separability
of the perturbation equations, which does not occur if the background is not
described by the Kerr metric as in our case (which is, in addition, only known
numerically). Likewise, the gravitational self-force was computed up to
second-order in perturbations Pound (2012); Pound  _et al._ (2020); Pound and
Wardell (2021) for orbits around BHs van de Meent (2018); Spiers , but the
case of ECOs is an uncharted territory.

In the absence of a consistent framework to study EMRI dynamics around a non-
Kerr spinning object, here we use the only framework available at the moment,
namely approximate semi-relativistic inspirals Glampedakis  _et al._ (2002);
Gair and Glampedakis (2006) and waveforms with methods known as “kludge”
schemes Ruffini and Sasaki (1981); Sasaki (1983); Babak  _et al._ (2007).
Kludge waveforms are constructed through the so-called numerical kludge
scheme, namely by combining flat spacetime weak-field (PN) GW emission
together with a fully-relativistic treatment for the secondary’s motion.

###  V.1 Numerical kludge scheme

To approximate EMRIs around rotating BSs we integrate the second-order
dynamical system for $r,\,\theta$, augmented with weak-field PN fluxes for the
energy and angular momentum loss due to GW emission Glampedakis and Kennefick
(2002); Glampedakis  _et al._ (2002); Gair and Glampedakis (2006). This
treatment, though approximate and valid only for small orbital velocities,
takes into account the dominant contribution of the secondary’s radiative
backreaction to the primary’s geometry, at second PN order, and results in an
adiabatic evolution of the EMRI. During the evolution, the orbit is treated,
at small timescales, as a geodesic, while for longer timescales the trajectory
is slowly driven through consecutively damped geodesics. This scheme has been
shown to perform well when compared to Teukolsky-based waveforms of EMRIs
Babak  _et al._ (2007).

BSs have non-trivial multipolar structure which differs from the one of a Kerr
BH Ryan (1995, 1997b, 1997c); Pacilio  _et al._ (2020); Vaglio  _et al._
(2022, 2023). At second PN order, the kludge scheme Glampedakis  _et al._
(2002) involves the mass quadrupole moment $M_{2}$. Thus, to construct a more
faithful (though still approximate) inspiral around a rotating BS, we augment
the fluxes with its modified mass quadrupole moment (see Barack and Cutler
(2007); Gair  _et al._ (2008); Apostolatos  _et al._ (2009); Lukes-
Gerakopoulos  _et al._ (2010); Destounis  _et al._ (2021); Destounis and
Kokkotas (2021)). This kludge scheme, together with the modified mass
quadrupole moment, has recently been examined and found to provide results
qualitatively equivalent to evolutions with instantaneous self-force in non-
Kerr electromagnetic analogs, which indicates that this method can, in
principle, describe resonance and island-crossings in nonintegrable EMRIs with
sufficient accuracy Mukherjee  _et al._ (2023).

We employ linear variations of $E$ and $L_{z}$ in an iterative way, such that
Canizares  _et al._ (2012); Destounis  _et al._ (2021); Destounis and Kokkotas
(2021)

$\displaystyle E_{1}$
$\displaystyle=\frac{E(0)}{m}+\bigg{<}\frac{dE}{dt}\bigg{>}\bigg{|}_{t=0}N_{r}\,T_{r},$
(14) $\displaystyle L_{z,1}$
$\displaystyle=\frac{L_{z}(0)}{m}+\bigg{<}\frac{dL_{z}}{dt}\bigg{>}\bigg{|}_{t=0}N_{r}\,T_{r},$
(15)

where $E(0),\,L_{z}(0)$ are the energy and $z$-component of the angular
momentum at $t=0$, respectively. In turn,
$\langle{dE}/{dt}\rangle|_{t=0},\,\langle{dL_{z}}/{dt}\rangle|_{t=0}$ are the
averaged PN fluxes calculated at the beginning of the inspiral, through the
complicated equations outlined in Gair and Glampedakis (2006); Destounis and
Kokkotas (2021). $N_{r}$ is the number of radial periods elapsed between each
update of (14) and (15), which improves the scheme and includes cumulative
nonlinear variations, while $T_{r}$ is the radial period of the EMRI.
Equations (14) and (15) are iterated along the whole EMRI evolution with
appropriate choices of $N_{r}$ and $T_{r}$ to obtain the dissipative orbit. If
the system was integrable, then through the evolution of the Carter constant,
we would be able to also evolve the rest of the components of the angular
momentum, namely $L_{x}$ and $L_{y}$, which are imprinted in Carter’s
constant. Our case though is nonintegrable and lacks of a separation Carter-
like constant; therefore, we have no explicit way of evolving $L_{x}$ and
$L_{y}$ but rather we keep them constant.

Even though, due to dissipation, the constraint equation (5) is not anymore a
constant of motion (that can be monitored to assess the accuracy of the
numerical scheme), we have tested numerous random time instants of the
inspirals obtained with the aforementioned scheme from which we extract the
instantaneous position and velocity vectors of the orbit, $E$ and $L_{z}$, and
calculate (5) for these values. We then compare the resulting constraint with
the one obtained by evolving a geodesic with the aforementioned parameters of
the same time instants as initial conditions. For all cases, we find that the
new dissipated constraint is satisfied to within 1 part in $10^{7}$ for the
first $\sim 10^{4}$ revolutions. Through geodesic evolutions of successive
time instants of an inspiral, we can built dissipative Poincare maps with
adiabatically decreasing $E$ and $L_{z}$ and therefore calculate time-
dependent rotation curves Destounis and Kokkotas (2021); Lukes-Gerakopoulos
and Witzany (2021).

Figure 11: Left: Sustained resonance for Case 3, with mass ratio
$m/M=10^{-6}$, and simplified linear fluxes without any updates as given in
Eqs. (14) and (15). The secondary’s initial parameters and conditions are
$E(0)/m=0.95$, $L_{z}(0)/m=3M$, $r(0)=6.5M$, $\dot{r}(0)=0.155$,
$\theta(0)=\pi/2$, while $\dot{\theta}(0)$ is defined from the constraint Eq.
(5). The total evolution time is $t=5\times 10^{6}M\sim 10$ months. The
linearly dissipative inspiral spends roughly $1$ month off resonance and $9$
months in perfect resonance. Here, we have chosen a stroboscopic depiction for
the trapping to be easier identified, i.e. from the time series of the full
Poincaré map, only every third consequent point is kept. Right: Same as in
left panel but starting from a point of the inspiral where the $1$ month off
resonance have lapsed and the orbit enters the (putatively) eternal resonance.

###  V.2 GW modeling

To present the phenomenological imprints of inspiraling secondaries onto
supermassive BSs, we employ the quadrupole approximation. In the traceless and
tranverse gauge the metric perturbations read (e.g., Misner  _et al._ (1973))

$h_{ij}=\frac{2}{d}\frac{d^{2}Q_{ij}}{dt^{2}},$ (16)

where $Q_{ij}$ is the symmetric and trace-free (STF) mass quadrupole tensor,
which can be written as

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

with $t$ the coordinate time of the secondary source, and

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

Here, we employ the approximation where quasi-isotropic coordinates at
infinity are identified with spherical ones as $Z(t)=(x(t),y(t),z(t))$, where

$\displaystyle x(t)$ $\displaystyle=r(t)\sin\theta(t)\cos\phi(t),$ (19)
$\displaystyle y(t)$ $\displaystyle=r(t)\sin\theta(t)\sin\phi(t),$ (20)
$\displaystyle z(t)$ $\displaystyle=r(t)\cos\theta(t),$ (21)

and then transform them from spherical to Euclidean coordinates that describe
the secondary’s trajectory.

An incoming wave from an EMRI onto an interferometer can be projected in two
polarizations, $+$ and $\times$, by introducing two unit vectors
$\bm{\hat{p}}=\bm{\hat{n}}\times\bm{\hat{z}}/|\bm{\hat{n}}\times\bm{\hat{z}}|$
and $\bm{\hat{q}}=\bm{\hat{p}}\times\bm{\hat{n}}$ (here $\times$ is the cross
product of two vectors and should not be confused with the cross polarization
symbol of the incoming GW), which are defined in terms of a third unit vector
$\bm{\hat{n}}$ that points from the EMRI source to the detector. Finally, the
unit vector $\bm{\hat{z}}$ designates the spin direction of the BS. The
triplet $\bm{\hat{p}},\,\bm{\hat{q}},\,\bm{\hat{n}}$ forms an orthonormal
basis from which the polarization tensor components are defined as

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

and allow us to write the metric perturbation in the quadrupole approximation
as

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

where $h_{+,\times}$ are the plus and cross polarizations of the incoming GW.
The GW components are then expressed in terms of the position, $Z^{i}(t)$,
velocity, $v^{i}(t)=dZ^{i}/dt$, and acceleration,
$a^{i}(t)=d^{2}Z^{i}/dt^{2}$, vectors. One finally obtains

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

LISA’s response to an incident GW event is correlated with the antenna pattern
functions, $F^{+,\times}_{I,II}(t)$, describing the motion of the detector on
its respective spacecraft channels $I$ and $II$ (see Refs. Cutler (1998);
Barack and Cutler (2004); Destounis  _et al._ (2020) for details). The total
waveform detected by a LISA-like interferometer reads

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

We assume a detector that lies at a luminosity distance $d$ with fixed
orientation $\bm{n}=(0,0,1)$ with respect to the EMRI source, and that the
primary’s polar and azimuthal angles are fixed at the equatorial plane. The
data streams that will be considered in what follows will contain the GW
together with stationary and Gaussian noise. For simplicity, we will further
assume that the two data stream channels are uncorrelated, thus we will abide
to a single-channel approximation of our detector.

##  VI EMRIs around a spinning BS

In this section, we discuss the effects of non-integrability and chaos on EMRI
evolution and GW emission. Even though approximate, the results shown below
demonstrate the basic features of the chaotic phenomena taking place in
rotating BSs. Since the most interesting cases of the BSs constructed are
Cases 1 and 3, we will focus on them from now on. We fix $E(0)/m=0.95$,
$L_{z}(0)/m=3M$ as initial parameters of the secondary in Eqs. (14) and (15),
respectively, and perform a different numbers of updates depending on how
quick the inspiral evolves, namely how close we are to the inner boundary of
the CZV (separatrix).

Note that in some cases the inspiral would occur _inside_ the BS. This is a
striking difference with respect to the BH case in which the signal disappears
after the particle has crossed the horizon. We can, therefore, continue the
evolution as long as the orbits do not cross the separatrix, where the
inspiral plunges.

It is worth noticing that the potential inside the BS is approximately
constant (see Fig. 1) and not particularly strong, which justifies the use of
PN fluxes even inside the star. For EMRIs around Kerr, the overlap between PN
and Teukolsky-based waveforms starts to completely deteriorate typically when
the periapsis radius $r_{p}\sim 5M$ (in Boyer-Lindquist coordinates) or
smaller Babak  _et al._ (2007), but in the BH case the potential is stronger,
as shown in Fig. 1, due to the larger compactness at the horizon. Thus, our PN
approximation should be accurate enough also inside the star all the way to
the edge of the CZV, especially for the least compact configurations.
Furthermore, we stress that we are only considering radiation-reaction
effects, neglecting environmental effects such as dynamical friction and
accretion within the BS Macedo  _et al._ (2013a), as well as direct (non-
gravitational) coupling between the secondary and the scalar field, both of
which might significantly contribute to the inspiral.

###  VI.1 Sustained resonances

A first interesting result arises when we ignore the updates on the flux, Eqs.
(14) and (15), an approximation that has been considered in various analyses
Apostolatos  _et al._ (2009); Lukes-Gerakopoulos  _et al._ (2010); Lukes-
Gerakopoulos (2012). This simplification leads to the following linearly-
variated fluxes

$\displaystyle E(t)$
$\displaystyle=\frac{E(0)}{m}+\bigg{<}\frac{dE}{dt}\bigg{>}\bigg{|}_{t=0}t,$
(26) $\displaystyle L_{z}(t)$
$\displaystyle=\frac{L_{z}(0)}{m}+\bigg{<}\frac{dL_{z}}{dt}\bigg{>}\bigg{|}_{t=0}t,$
(27)

Such assumption inevitably leads to sustained resonances, that cannot exist in
Kerr van de Meent (2014), where the secondary is trapped in a resonant island
for an extremely long (potentially infinite) time. In Fig. 11 we present one
such case, where a $10$ month inspiral spends $1$ month evolving normally and
$9$ months in perfect resonance. It never eventually escapes the island, at
least for the timescale over we have evolved the EMRI. In order to identify
the trapping more easily, Fig. 11 presents a stroboscopic depiction, i.e. from
the time series of the full dissipative Poincaré map, only every third
consequent point is kept due to the multiplicity of the $2/3$ resonance.

The initial position for a sustained resonance to occur is not fine-tuned but
rather a quite large region of initial conditions $r(0)\in[6.46,6.58]M$ exists
(which corresponds to $[7.47,7.59]M$ in Boyer-Lindquist coordinates), when
fixing the rest of the parameters as $\dot{r}(0)=0.155$, $\theta(0)=\pi/2$,
and define $\dot{\theta}(0)$ through the constraint equation (5). Similar
results have recently been found in Lukes-Gerakopoulos and Witzany (2021) for
linear variations of the fluxes like Eq. (26). Even if the range of initial
conditions giving such sustained resonances is not negligible, in practice
such sustained resonances are an artifact due the assumption of low-order
energy and angular-momentum fluxes throughout the evolution. Namely, the
linear approximation in energy and angular momentum is valid as long as the
trajectory of the EMRI does not get too far from the initial geodesic, on
which the orbital elements and fluxes have been initialized.

Figure 12: External $2/3$ resonant island crossing for Case 1. To produce the
inspiral we have updated the fluxes $300$ times. The secondary’s initial
parameters and conditions are $E(0)/m=0.95$, $L_{z}(0)/m=3M$, $r(0)=5.026M$,
$\dot{r}(0)=0.1$, $\theta(0)=\pi/2$, while $\dot{\theta}(0)$ is defined from
the constraint Eq. (5). The EMRI spends $\sim 200$ cycles in resonance.

By taking into account updates (from $50$ to $150$ for some representative
cases) on the fluxes for the aforementioned radial range, which effectively
includes nonlinear terms on the fluxes, all sustained resonances vanish. Of
course, we cannot exclude the possibility that the region of sustained
resonances either shrinks significantly or hides into a different range of
initial conditions, although we consider this option as unlikely. Note also
that putative sustained resonances should persist (and not disappearing) as
the number of flux updates is increased. A similar case was recently studied
for an EMRI analog in Mukherjee  _et al._ (2023), where lower order flux
approximations gave rise to sustained resonances whereas, when higher order
terms were involved, the sustained resonances disappeared.

###  VI.2 Time-dependent rotation curves

Figures 12 and 13 show typical, time-dependent, island crossings for Cases 1
and 3. In Fig. 12, the non-zero initial velocity leads to the EMRI crossing
the $2/3$ island in the exterior of Case 1 configuration. The plateau is
finite and lasts for around 200 revolutions before it exits the island. The
fact that the rotation number decreases with time shows that the orbital
motion is outside of the star and tends towards circularization. On the other
hand, in Fig. 13, we show the dissipative rotation curve of the $4/7$ island
crossing which resides inside the BS of Case 3. The fluxes have increased at
this point significantly, therefore the EMRI evolves faster. Nevertheless, the
secondary spends $\sim 100$ cycles in perfect resonance which is extraordinary
for a subdominant resonance. Equivalently, the fact that these orbits are
inside the star and both the eccentricity and rotation numbers increase,
justifies why the dissipative rotation curve increases with time444It has been
shown that, before plunge, Kerr EMRIs enter a short phase of increasing
eccentricity Gair and Glampedakis (2006). In our case, due to the absence of
an event horizon, the eccentricity EMRI can enter the star and the
eccentricity continues increasing and leads to orbits as those found in
Grandclement  _et al._ (2014). Since the eccentricity is directly proportional
to the periapsis-apoapsis oscillation frequency $\omega_{r}$, the rotation
number increases as well..

Figure 13: Internal $4/7$ resonant island crossing for Case 3. To produce the
inspiral we have updated the fluxes $250$ times. The secondary’s initial
parameters and conditions are $E(0)/m=0.95$, $L_{z}(0)/m=3M$, $r(0)=4.06025M$,
$\dot{r}(0)=0$, $\theta(0)=\pi/2$ while the $\dot{\theta}(0)$ is defined from
the constraint Eq. (5). The EMRI spends $\sim 100$ cycles in resonance. Figure
14: Left: Interior $2/3$ resonant island crossing for Case 3. To produce the
inspiral we have updated the fluxes $200$ times. The secondary’s initial
parameters and conditions are $E(0)/m=0.95$, $L_{z}(0)/m=3M$, $r(0)=6.5M$,
$\dot{r}(0)=0.155$, $\theta(0)=\pi/2$, while $\dot{\theta}(0)$ is defined from
the constraint Eq. (5). The EMRI spends $\sim 10$ cycles on and off the
resonant island and $\sim 100$ revolutions in the vicinity the $2/3$ Birkhoff
chain until the chaotic layer and island are crossed. Right: The produced KAM
curves from the red dots shown in the inset of the left panel. All of them
(namely, $10$ points) belong to the $2/3$ resonant island.

The final, and most interesting scenario examined for Case 3 is the very large
plateau that arises when $\dot{r}(0)=0.155$, that begins from the exterior and
ends in the interior of the BS. Even though one would expect a dissipative
rotation curve with a plateau that lasts for thousands of cycles, the
existence of a thick chaotic layer around the island induces a significant
alteration to the rotation curve, which we find here for the first time. In
Fig. 14 we present the behavior of the wide $2/3$ island. On the left, we show
that instead of a plateau we encounter mostly the chaotic layer (regardless of
the initial condition $r(0)$), where the rotation numbers are ill-defined.
Nonetheless, there are points (marked in red in the zoom-in inset of the left
panel in Fig. 14) that actually enter the island. To make sure that these
points do enter the resonance, we plot on the right the corresponding KAM
curves of these points. As shown, they form islands with tips that do not
touch, therefore the EMRI is not occupying the chaotic layer, unlike the rest
of the black points in the inset of Fig. 14 that correspond to orbits residing
in the chaotic layer. What seems to occur here is a strong chaotic shielding
which does not allow the EMRI to spend continuous time intervals inside the
island. Even so, we should expect interesting GW imprints for all cases
discussed above and especially those from Case 3 which take place in the
interior of the rotating BS. We present those effects in the next section.

We stress that our results are based on PN fluxes and neglect effects other
than radiation-reaction within the BS, e.g. dynamical friction and accretion,
which can nevertheless be added into the evolution.

Figure 15: Frequency evolution of an EMRI around rotating BS from Case 1,
through the $2/3$ resonant island with parameters and initial conditions as in
Fig. 12. As a reference, the approximate GWs (fundamental frequency on the
left panel and first harmonic on the right panel) are detected at luminosity
distance $d=100\,\text{Mpc}$. Figure 16: Same as Fig. 15 but for the rotating
BS configuration of Case 3, through the $4/7$ interior resonant island with
parameters and initial conditions as in Fig. 13. Figure 17: Same as Fig. 15
but for the rotating BS configuration of Case 3, through the $2/3$ interior
Birkhoff chain with parameters and initial conditions as in Fig. 14.

###  VI.3 GW frequency evolution

After obtaining the inspirals for Case 1 and 3, it is straightforward to use
Eq. (25) in order to find the approximate GW emitted by such EMRIs. The time-
domain waveforms are Fourier transformed to the frequency domain where the
frequency evolution of the inspiral can be constructed. Here, we follow
Destounis  _et al._ (2021); Destounis and Kokkotas (2021); Destounis  _et al._
(2023a) and construct the spectrograms by performing consecutive short-time
Fourier transforms, with appropriate window sizes and offsets in order to
maximize the quality of the resulting figures.

In Fig. 15, we present the most typical glitch waveform for Case 1, where the
BS is rather compact and mimics a non-Kerr BH. Its fundamental and first-
harmonic frequency evolution display the standard modulation that occurs when
the EMRI crosses the $2/3$ island (that lies on the exterior of the rotating
BS and lasts for 200 revolutions). This should lead to a significant dephasing
with respect to an integrable EMRI evolution that crosses the same resonance.
The resemblance with glitches found in Destounis  _et al._ (2021); Destounis
and Kokkotas (2021) for non-Kerr spacetimes is remarkable. Therefore, since we
use the hybrid kludge scheme and the quadrupole formula, in this case it would
be almost impossible to distinguish between non-Kerr BHs or rotating BS based
only on these glitches. Fortunately, Case 3, the less compact configuration we
have constructed, changes this picture completely. Figure 16 shows the
frequency evolution of an EMRI through the internal $4/7$ island. The
frequency evolution becomes nonlinear due to the acceleration of the fluxes
and the EMRI frequencies after the glitch increase towards a chirp.

Finally, we examine the spectrogram of an EMRI crossing the internal $2/3$
resonant island that corresponds to the inspiral shown in Fig. 14. The GW
frequency modulation shown in Fig. 17 is the first one of its kind found in
any non-integrable spacetime so far. We do not only observe similar nonlinear
frequency evolution as a function of time for the strongest Fourier peaks
(fundamental and first harmonic) as in Fig. 16, but also linear-in-time
subdominant frequency peaks that evolve and glitch in a similar manner. We
have performed the same inspiral for two different initial conditions $r(0)$,
one that lies outside and one that lies inside the range of sustained
resonances (if linear fluxes are assumed) and the resulting spectrogram is
qualitatively the same as that shown in Fig. 17. This feature is significantly
different from all the other glitches found in previous non-integrable
spacetimes, and is associated to the fact that the plateau in the dissipative
rotation curve is replaced with an observable chaotic layer with glimpses of
island occupancy.

When an orbit is fully chaotic and resides in a chaotic sea Contopoulos
(2003), then we expect rotation curves and spectrograms to be rendered useless
since no well-defined rotation numbers and discrete GW frequencies exist
Kiuchi and Maeda (2004); Gair  _et al._ (2008). On the other hand, when an
orbit is close to a slightly chaotic region but otherwise not a fully chaotic
part of phase space, such as the chaotic layer that shields the interior $2/3$
island of Case 3, the waveform and its frequency content resemble a lot those
of a regular orbit, with discrete Fourier peaks and indistinguishable effects
of chaos in the time-domain waveform. The only key difference is that the
discrete Fourier peaks are not comprised by single harmonics (for geodesics)
but rather harmonics that are broken down to subpeaks. The spectrum,
nevertheless, remains discrete. This has been shown in Kiuchi and Maeda (2004)
for a secondary particle with spin orbiting around a Schwarzschild BH (see
Figs. 4, 10, 11 and 13 therein), although in that case this behavior
disappears as the secondary spin is treated perturbatively, as also requested
for consistency within the perturbative expansion in the mass ratio (see,
e.g., Piovano  _et al._ (2020)). It is therefore interesting that we recover
the same feature, that is a slightly chaotic region that embeds the island,
but otherwise normal generic orbits around this region, discrete Fourier
spectra, and harmonics with multiple peaks, but in rotating BS spacetimes
(where motion is not even integrable in a perturbative sense). This is the
reason why Fig. 17 has such a distinctive and peculiar nature, which can be
explained by the effects discussed in Kiuchi and Maeda (2004). If we consider
the effects on the spectrogram and the dissipative evolution through the
Birkhoff chain, the EMRIs spends around 100 cycles to cross it, and 10 cycles
in perfect resonance. Yet, a Birkhoff chain is composed of both stable and
unstable periodic points, and the chaotic layer helps to amplify the frequency
glitch in a significant way, with frequency jumps as large as $\sim 0.07$ mHz
for the most dominant fundamental Fourier peak evolution, which is almost one
order of magnitude larger than any other frequency glitches presented here,
and also higher than those found in Destounis  _et al._ (2021); Destounis and
Kokkotas (2021) for various deformation parameters that were either
exaggerated or chosen arbitrarily. Our case presents a well motivated model of
ECOs that has a precise formation scenario, thus the results presented here
can be used in the future as a testbed to understand if weak or strong (or no)
chaos exists in GW observations.

##  VII Conclusions

The detection of EMRIs by LISA Amaro-Seoane  _et al._ (2017) and other
spaceborne detectors Ruan  _et al._ (2020a, b); Luo  _et al._ (2016) will
provide for the first time an extremely accurate mapping of the primary’s
spacetime geometry, due to the large mass hierarchy that creates the
conditions for long-lasting inspirals with continuously observable GW
emission. The associated signal will present very rich phenomenology Amaro-
Seoane  _et al._ (2022) like the appearance of resonances Apostolatos  _et
al._ (2009); Contopoulos  _et al._ (2011); Lukes-Gerakopoulos (2012); Lukes-
Gerakopoulos and Witzany (2021) and, remarkably, carries information about
spacetime symmetries Destounis  _et al._ (2020); Lukes-Gerakopoulos  _et al._
(2010).

The observation of these systems, will also help to decide if BHs are the only
compact objects in the universe or if other exotic configurations exist in
nature. The different multipolar structure, relative to Kerr BHs, can lead to
fascinating phenomena around resonances such as frequency modulation and GW
glitches Destounis  _et al._ (2021), in particular for those spacetimes that
lack a Carter constant (or any other higher-rank Killing tensor) that would
guarantee the integrability of geodesics. So far, only bumpy and non-Kerr
solution of GR – that are plagued with pathologies Destounis and Kokkotas
(2021); Destounis  _et al._ (2023a) and do not come from a well motivated,
first-principle theory – have been analyzed in order to inspect their
different phenomenology.

Here, we considered for the first time one of the most well-behaved and simple
compact object that may as well exist in our universe, namely rotating BSs.
These exotic objects can serve as prototypical BH mimickers, including
supermassive ones, do not have an event horizon or singularities in their
interior and have a clear formation mechanism. They also constitute compelling
dark matter candidates, and could have formed in the early Universe, e.g. as
remnants of inflation.

While the majority of analyses for BSs have been performed for equatorial
orbits, here we took a different turn on examining supermassive self-
interacting rotating BSs in the context of EMRIs for generic orbits (without
constraints on the orbital plane or eccentricity). In an attempt to study
geodesics that are peculiar enough to be distinguishable from those around
Kerr BHs, we evolved geodesics around and inside supermassive rotating
configurations and studied their GW signatures. We have built numerically
three of such configurations with representative values of the spin and
compactness.

Our geodesic analysis reveals the existence of resonant islands (with finite
width) and Birkhoff chains in the phase space of orbits, which signals that
rotating BSs are non-integrable. We also found that compact rotating BSs
behave similarly to non-Kerr BHs up to the point where the geodesic enters the
star smoothly and the godesic structure becomes qualitatively different.
Decreasing the compactness makes things much more different between rotating
BSs and non-Kerr BHs. First, the geodesic phase space hosts many orbits that
reside in the star’s interior, where the rotation curve exhibits a change in
monotonicity and begins to increase due to the increment of eccentricity
Grandclement  _et al._ (2014). The latter effect leads to a newly-found
phenomenon, that is the existence of degenerate islands of stability, in the
interior and the exterior of the star, with the same rotation number. We have
also found regions in the phase space of geodesics where there exist not only
a plethora of resonant islands, but also a thin chaotic layer surrounding the
most dominant $2/3$ island of stability. This is, to our knowledge, the first
time that a full Birkhoff chain appears in a general relativistic setup, which
has important implications in GWs emitted by EMRIs that cross this region.
With appropriate initial conditions, we find that the $2/3$ island encounters
a plateau in the rotation curve, where the geodesic enters the island in the
exterior and exits in the interior of the star; a fascinating outcome of the
fact that the compactness of the BS of Case 3 (albeit as large as a typical
neutron star) is relatively smaller than in the other cases.

To achieve some initial estimates of the elementary structure of approximate
waveforms from such rotating BSs, and especially for generic inspirals that
cross transient resonances, we use the quadrupole approximation to model GWs
and evolve the inspiral, with augmented PN fluxes to account for the modified
multipolar structure of our configurations Vaglio  _et al._ (2022, 2023).
Despite these approximations, the phenomenology of the inspiral and GW signal
is expected to be robust.

At first sight with a linear approximation of the fluxes, we stumble upon
inspirals that enter the $2/3$ island, in a generous range of initial
conditions $r(0)$, and never exit the resonance, as also found in Lukes-
Gerakopoulos and Witzany (2021) but in a much smaller radial domain. However,
when we consistently update fluxes during the inspiral (though in a
discretized way after some amount of revolutions), higher order terms seem to
destroy such sustained resonances. Although unlikely, we cannot unequivocally
exclude the existence of other regions of phase space where such resonances
occur even when updates in the fluxes are introduced. We sketch the
dissipative rotation curves by using time instants of each EMRI as initial
condition to a geodesic evolution and find that external resonant islands are
indistinguishable with those occurring in non-Kerr objects Destounis  _et al._
(2021); Destounis and Kokkotas (2021); Destounis  _et al._ (2023a). The
plateaus that characterize islands situated in the interior of the star are
approached from below in the rotation curves, which then reprise a
monotonically increasing behavior. This corresponds to a reversed time and
radial dependence with respect to the rotation curves in the neighbourhood of
external islands, related to the non-zero energy-density distribution inside
the star. Such monotonicity reversal also pairs with the different evolution
of the eccentricity. Indeed, inside the star, eventually a point is reached
where the circularization of orbits have lapsed and the eccentricity increases
with decreasing radius, in contrast to what happens in the exterior. The most
prominent island of all found to date, the one for Case 3 that has a record
width of $2.6M$ (in quasi-isotropic coordinates), is even more complicated due
to chaotic shielding. The EMRI does not enter and stay in the island but
rather oscillates between on and off island states. It manages to do so at
least 10 times during which the evolution drives the orbit around the chaotic
layer till the secondary eventually exits the particular part of the Birkhoff
chain.

Overall, the time spent in these cases spans from $200$ cycles in external
resonances, where the evolution is adiabatically slow, to around $100$ cycles
in the interior of the star, when the Birkhoff chain is taken into account as
a part of the whole chaotic KAM curve (and not just the island with the
periodic stable point at its center).

When the EMRI crosses a resonant island, all the above effects are imprinted
onto the EMRI waveform either in a typical way, analogous to other non-Kerr
systems, or as novel signatures that designate the existence of a supermassive
BS primary. GW signals from the exterior islands are qualitatively
indistinguishable from those in non-Kerr EMRIs; a result which is expected
from the close similarity in the rotation curves discussed above. The case of
internal motion is instead quite different. Even subdominant resonant islands
introduce significant nonlinear frequency modulation because the orbits is
close to exit the CZV (see e.g. Fig. 16). The widest resonant island of Case 3
has a very special behavior imprinted in the GW when the EMRI crosses it. The
proximity of the orbits to the otherwise mild (but still observable) chaotic
layer leads to an effect in the GW Fourier peaks first found in Kiuchi and
Maeda (2004). Moderate chaotic layers do not produce continuous GW spectra, as
fully chaotic orbits do, but rather discrete sets of harmonics that are broken
down to sub-harmonics. This phenomenon is evident in Fig. 17 where two
simultaneous glitches occur under to the evolution of each sub-frequency due
to mild chaos. The underlying reason is, again, the fact that the EMRI is
close to crossing a thin chaotic region. This turns each single harmonic of
the frequency content of the approximated GW signal into multiple sub-peaks,
that arise in the spectrogram in the particular fashion shown in Fig. 17. The
most prominent (the brightest) one, which leads to a frequency glitch of
${\cal O}(0.1)\,{\rm mHz}$, renders all exterior glitches subdominant and
should affect significantly the orbital evolution and eventual parameter
estimation. Yet, we stress that our results are qualitative and should be
extended with more accurate inspiral and waveform modeling.

In this respect it would be very interesting to perform proper perturbation
theory around a spinning BS, although nonseparability of the equations makes
the flux computation quite challenging compared to the Teukolsky case for a
Kerr BH. Nevertheless our results provide solid evidence that chaos does exist
in EMRIs around rotating BSs, and the fact that this is a quite compelling
model for ECOs should drive us to further understand these objects, as well as
their effects in geodesics and EMRI dynamics.

Along these lines it would be relevant to include environmental effects such
as dynamical friction and accretion Macedo  _et al._ (2013a); Cardoso  _et
al._ (2022a, b); Destounis  _et al._ (2023b) and assess their impact on the
chaotic motion and EMRI signal. Finally, although we focused on a specific
model of BS with large quartic interactions, we expect to find the same
qualitative features in other spinning BS models, spinning Proca stars, and in
fact in all spinning ECOs where geodesic motion is most likely nonintegrable.
In particular, we expect the same (and perhaps even more prominent) signatures
of chaotic motion for EMRIs in generic orbits around BH microstates predicted
in the fuzzball scenario (see Bena  _et al._ (2022) for a recent review), due
to their rich multipolar structure and breaking of the axial and equatorial
symmetry Bena and Mayerson (2020); Bianchi  _et al._ (2020); Bena and Mayerson
(2021); Bianchi  _et al._ (2021); Fransen and Mayerson (2022).

###### Acknowledgements.

The authors are indebted to Georgios Lukes-Gerakopoulos and Kostas D. Kokkotas
for insightful discussions. We acknowledge financial support provided under
the European Union’s H2020 ERC, Starting Grant agreement no. DarkGRA–757480
and support under the MIUR PRIN (Grant 2020KR4KN2 “String Theory as a bridge
between Gauge Theories and Quantum Gravity”) and FARE (GW-NEXT, CUP:
B84I20000100001, 2020KR4KN2) programmes. K.D. is supported by the DAAD program
for the “promotion of the exchange and scientific cooperation between Greece
and Germany IKYDAAD 2022” (57628320).

## References

  * Abbott  _et al._ (2021a) R. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. X 11, 021053 (2021a), arXiv:2010.14527 [gr-qc] .
  * Abbott  _et al._ (2021b) R. Abbott _et al._ (LIGO Scientific, KAGRA, VIRGO), Astrophys. J. Lett. 915, L5 (2021b), arXiv:2106.15163 [astro-ph.HE] .
  * Abbott  _et al._ (2021c) R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), (2021c), arXiv:2111.03606 [gr-qc] .
  * Abbott  _et al._ (2020a) R. Abbott _et al._ (LIGO Scientific, Virgo), Astrophys. J. Lett. 896, L44 (2020a), arXiv:2006.12611 [astro-ph.HE] .
  * Abbott  _et al._ (2020b) R. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 125, 101102 (2020b), arXiv:2009.01075 [gr-qc] .
  * Abbott  _et al._ (2020c) R. Abbott _et al._ (LIGO Scientific, Virgo), Astrophys. J. Lett. 900, L13 (2020c), arXiv:2009.01190 [astro-ph.HE] .
  * Calderón Bustillo  _et al._ (2021) J. Calderón Bustillo, N. Sanchis-Gual, A. Torres-Forné, J. A. Font, A. Vajpeyi, R. Smith, C. Herdeiro, E. Radu, and S. H. W. Leong, Phys. Rev. Lett. 126, 081101 (2021), arXiv:2009.05376 [gr-qc] .
  * Cardoso and Pani (2019) V. Cardoso and P. Pani, Living Rev. Rel. 22, 4 (2019), arXiv:1904.05363 [gr-qc] .
  * Maggio  _et al._ (2021) E. Maggio, P. Pani, and G. Raposo,  (2021), 10.1007/978-981-15-4702-7_29-1, arXiv:2105.06410 [gr-qc] .
  * Bertone  _et al._ (2005) G. Bertone, D. Hooper, and J. Silk, Phys. Rept. 405, 279 (2005), arXiv:hep-ph/0404175 .
  * Clowe  _et al._ (2006) D. Clowe, M. Bradac, A. H. Gonzalez, M. Markevitch, S. W. Randall, C. Jones, and D. Zaritsky, Astrophys. J. Lett. 648, L109 (2006), arXiv:astro-ph/0608407 .
  * Bergstrom (2009) L. Bergstrom, New J. Phys. 11, 105006 (2009), arXiv:0903.4849 [hep-ph] .
  * Primack (1997) J. R. Primack, in  _Midrasha Mathematicae in Jerusalem: Winter School in Dynamical Systems_ (1997) arXiv:astro-ph/9707285 .
  * Corbelli and Salucci (2000) E. Corbelli and P. Salucci, Mon. Not. Roy. Astron. Soc. 311, 441 (2000), arXiv:astro-ph/9909252 .
  * Kolb and Turner (1990) E. W. Kolb and M. S. Turner,  _The Early Universe_, Vol. 69 (1990).
  * Carroll (2001) S. M. Carroll, Living Rev. Rel. 4, 1 (2001), arXiv:astro-ph/0004075 .
  * Del Popolo (2007) A. Del Popolo, Astron. Rep. 51, 169 (2007), arXiv:0801.1091 [astro-ph] .
  * Arbey and Mahmoudi (2021) A. Arbey and F. Mahmoudi, Prog. Part. Nucl. Phys. 119, 103865 (2021), arXiv:2104.11488 [hep-ph] .
  * Perivolaropoulos and Skara (2022) L. Perivolaropoulos and F. Skara, New Astron. Rev. 95, 101659 (2022), arXiv:2105.05208 [astro-ph.CO] .
  * Seidel and Suen (1994) E. Seidel and W.-M. Suen, Phys. Rev. Lett. 72, 2516 (1994), arXiv:gr-qc/9309015 .
  * Sanchis-Gual  _et al._ (2019) N. Sanchis-Gual, F. Di Giovanni, M. Zilhão, C. Herdeiro, P. Cerdá-Durán, J. A. Font, and E. Radu, Phys. Rev. Lett. 123, 221101 (2019), arXiv:1907.12565 [gr-qc] .
  * Bonazzola and Pacini (1966) S. Bonazzola and F. Pacini, Phys. Rev. 148, 1269 (1966).
  * Feinblum and McKinley (1968) D. A. Feinblum and W. A. McKinley, Phys. Rev. 168, 1445 (1968).
  * Kaup (1968) D. J. Kaup, Phys. Rev. 172, 1331 (1968).
  * Ruffini and Bonazzola (1969) R. Ruffini and S. Bonazzola, Phys. Rev. 187, 1767 (1969).
  * Colpi  _et al._ (1986) M. Colpi, S. L. Shapiro, and I. Wasserman, Phys. Rev. Lett. 57, 2485 (1986).
  * van der Bij and Gleiser (1987) J. J. van der Bij and M. Gleiser, Phys. Lett. B 194, 482 (1987).
  * Frieman  _et al._ (1988) J. A. Frieman, G. B. Gelmini, M. Gleiser, and E. W. Kolb, Phys. Rev. Lett. 60, 2101 (1988).
  * Frieman  _et al._ (1989) J. A. Frieman, A. V. Olinto, M. Gleiser, and C. Alcock, Phys. Rev. D 40, 3241 (1989).
  * Ferrell and Gleiser (1989) R. Ferrell and M. Gleiser, Phys. Rev. D 40, 2524 (1989).
  * Grasso (1990) D. Grasso, Phys. Rev. D 41, 2998 (1990).
  * Tkachev (1991) I. I. Tkachev, Phys. Lett. B 261, 289 (1991).
  * Liddle and Madsen (1992) A. R. Liddle and M. S. Madsen, Int. J. Mod. Phys. D 1, 101 (1992).
  * Jetzer (1992) P. Jetzer, Phys. Rept. 220, 163 (1992).
  * Schunck and Mielke (2003) F. E. Schunck and E. W. Mielke, Class. Quant. Grav. 20, R301 (2003), arXiv:0801.0307 [astro-ph] .
  * Liebling and Palenzuela (2023) S. L. Liebling and C. Palenzuela, Living Rev. Rel. 26 (2023), 10.1007/s41114-023-00043-4, arXiv:1202.5809 [gr-qc] .
  * Kesden  _et al._ (2005) M. Kesden, J. Gair, and M. Kamionkowski, Phys. Rev. D 71, 044015 (2005), arXiv:astro-ph/0411478 .
  * Macedo  _et al._ (2013a) C. F. B. Macedo, P. Pani, V. Cardoso, and L. C. B. Crispino, Phys. Rev. D 88, 064046 (2013a), arXiv:1307.4812 [gr-qc] .
  * Macedo  _et al._ (2013b) C. F. B. Macedo, P. Pani, V. Cardoso, and L. C. B. Crispino, Astrophys. J. 774, 48 (2013b), arXiv:1302.2646 [gr-qc] .
  * Macedo  _et al._ (2016) C. F. B. Macedo, V. Cardoso, L. C. B. Crispino, and P. Pani, Phys. Rev. D 93, 064053 (2016), arXiv:1603.02095 [gr-qc] .
  * Cardoso  _et al._ (2017) V. Cardoso, E. Franzin, A. Maselli, P. Pani, and G. Raposo, Phys. Rev. D 95, 084014 (2017), [Addendum: Phys.Rev.D 95, 089901 (2017)], arXiv:1701.01116 [gr-qc] .
  * Guo  _et al._ (2019) H.-K. Guo, K. Sinha, and C. Sun, JCAP 09, 032 (2019), arXiv:1904.07871 [hep-ph] .
  * Berti and Cardoso (2006) E. Berti and V. Cardoso, Int. J. Mod. Phys. D 15, 2209 (2006), arXiv:gr-qc/0605101 .
  * Pacilio  _et al._ (2020) C. Pacilio, M. Vaglio, A. Maselli, and P. Pani, Phys. Rev. D 102, 083002 (2020), arXiv:2007.05264 [gr-qc] .
  * Vaglio  _et al._ (2023) M. Vaglio, C. Pacilio, A. Maselli, and P. Pani, (2023), arXiv:2302.13954 [gr-qc] .
  * Seidel and Suen (1991) E. Seidel and W. M. Suen, Phys. Rev. Lett. 66, 1659 (1991).
  * Coleman (1985) S. R. Coleman, Nucl. Phys. B 262, 263 (1985), [Addendum: Nucl.Phys.B 269, 744 (1986)].
  * Kusenko and Shaposhnikov (1998) A. Kusenko and M. E. Shaposhnikov, Phys. Lett. B 418, 46 (1998), arXiv:hep-ph/9709492 .
  * Lozanov and Amin (2019) K. D. Lozanov and M. A. Amin, Phys. Rev. D 99, 123504 (2019), arXiv:1902.06736 [astro-ph.CO] .
  * Aurrekoetxea  _et al._ (2023) J. C. Aurrekoetxea, K. Clough, and F. Muia, (2023), arXiv:2304.01673 [gr-qc] .
  * Lozanov  _et al._ (2023) K. D. Lozanov, M. Sasaki, and V. Takhistov, (2023), arXiv:2304.06709 [astro-ph.CO] .
  * Franciolini and Pani (2023) G. Franciolini and P. Pani, (2023), arXiv:2304.13576 [astro-ph.CO] .
  * Gleiser (1988) M. Gleiser, Phys. Rev. D 38, 2376 (1988), [Erratum: Phys.Rev.D 39, 1257 (1989)].
  * Gleiser and Watkins (1989) M. Gleiser and R. Watkins, Nuclear Physics B 319, 733 (1989).
  * Kusmartsev  _et al._ (1991) F. V. Kusmartsev, E. W. Mielke, and F. E. Schunck, Phys. Rev. D 43, 3895 (1991), arXiv:0810.0696 [astro-ph] .
  * Ryan (1997a) F. D. Ryan, Phys. Rev. D 55, 6081 (1997a).
  * Guzman (2004) F. S. Guzman, Phys. Rev. D 70, 044033 (2004), arXiv:gr-qc/0407054 .
  * Levkov  _et al._ (2018) D. G. Levkov, A. G. Panin, and I. I. Tkachev, Phys. Rev. Lett. 121, 151301 (2018), arXiv:1804.05857 [astro-ph.CO] .
  * Veltmaat  _et al._ (2018) J. Veltmaat, J. C. Niemeyer, and B. Schwabe, Phys. Rev. D 98, 043509 (2018), arXiv:1804.09647 [astro-ph.CO] .
  * Amin and Mocz (2019) M. A. Amin and P. Mocz, Phys. Rev. D 100, 063507 (2019), arXiv:1902.07261 [astro-ph.CO] .
  * Arvanitaki  _et al._ (2020) A. Arvanitaki, S. Dimopoulos, M. Galanis, L. Lehner, J. O. Thompson, and K. Van Tilburg, Phys. Rev. D 101, 083014 (2020), arXiv:1909.11665 [astro-ph.CO] .
  * Bernal and Siddhartha Guzman (2006) A. Bernal and F. Siddhartha Guzman, Phys. Rev. D 74, 103002 (2006), arXiv:astro-ph/0610682 .
  * Palenzuela  _et al._ (2007) C. Palenzuela, I. Olabarrieta, L. Lehner, and S. L. Liebling, Phys. Rev. D 75, 064005 (2007), arXiv:gr-qc/0612067 .
  * Palenzuela  _et al._ (2008) C. Palenzuela, L. Lehner, and S. L. Liebling, Phys. Rev. D 77, 044036 (2008), arXiv:0706.2435 [gr-qc] .
  * Choptuik and Pretorius (2010) M. W. Choptuik and F. Pretorius, Phys. Rev. Lett. 104, 111101 (2010), arXiv:0908.1780 [gr-qc] .
  * Bezares  _et al._ (2017) M. Bezares, C. Palenzuela, and C. Bona, Phys. Rev. D 95, 124005 (2017), arXiv:1705.01071 [gr-qc] .
  * Palenzuela  _et al._ (2017) C. Palenzuela, P. Pani, M. Bezares, V. Cardoso, L. Lehner, and S. Liebling, Phys. Rev. D 96, 104058 (2017), arXiv:1710.09432 [gr-qc] .
  * Bezares and Palenzuela (2018) M. Bezares and C. Palenzuela, Class. Quant. Grav. 35, 234002 (2018), arXiv:1808.10732 [gr-qc] .
  * Helfer  _et al._ (2022) T. Helfer, U. Sperhake, R. Croft, M. Radia, B.-X. Ge, and E. A. Lim, Class. Quant. Grav. 39, 074001 (2022), arXiv:2108.11995 [gr-qc] .
  * Bezares  _et al._ (2022) M. Bezares, M. Bošković, S. Liebling, C. Palenzuela, P. Pani, and E. Barausse, Phys. Rev. D 105, 064067 (2022), arXiv:2201.06113 [gr-qc] .
  * Croft  _et al._ (2023) R. Croft, T. Helfer, B.-X. Ge, M. Radia, T. Evstafyeva, E. A. Lim, U. Sperhake, and K. Clough, Class. Quant. Grav. 40, 065001 (2023), arXiv:2207.05690 [gr-qc] .
  * Evstafyeva  _et al._ (2022) T. Evstafyeva, U. Sperhake, T. Helfer, R. Croft, M. Radia, B.-X. Ge, and E. A. Lim, (2022), arXiv:2212.08023 [gr-qc] .
  * Siemonsen and East (2023) N. Siemonsen and W. E. East, (2023), arXiv:2302.06627 [gr-qc] .
  * Di Giovanni  _et al._ (2020) F. Di Giovanni, N. Sanchis-Gual, P. Cerdá-Durán, M. Zilhão, C. Herdeiro, J. A. Font, and E. Radu, Phys. Rev. D 102, 124009 (2020), arXiv:2010.05845 [gr-qc] .
  * Siemonsen and East (2021) N. Siemonsen and W. E. East, Phys. Rev. D 103, 044022 (2021), arXiv:2011.08247 [gr-qc] .
  * Schunck and Torres (2000) F. E. Schunck and D. F. Torres, Int. J. Mod. Phys. D 9, 601 (2000), arXiv:gr-qc/9911038 .
  * Friedberg  _et al._ (1987) R. Friedberg, T. D. Lee, and Y. Pang, Phys. Rev. D 35, 3658 (1987).
  * Guerra  _et al._ (2019) D. Guerra, C. F. B. Macedo, and P. Pani, JCAP 09, 061 (2019), [Erratum: JCAP 06, E01 (2020)], arXiv:1909.05515 [gr-qc] .
  * Collodel  _et al._ (2017) L. G. Collodel, B. Kleihaus, and J. Kunz, Phys. Rev. D 96, 084066 (2017), arXiv:1708.02057 [gr-qc] .
  * Bošković and Barausse (2022) M. Bošković and E. Barausse, JCAP 02, 032 (2022), arXiv:2111.03870 [gr-qc] .
  * Collodel and Doneva (2022) L. G. Collodel and D. D. Doneva, Phys. Rev. D 106, 084057 (2022), arXiv:2203.08203 [gr-qc] .
  * Vincent  _et al._ (2016) F. H. Vincent, Z. Meliani, P. Grandclement, E. Gourgoulhon, and O. Straub, Class. Quant. Grav. 33, 105015 (2016), arXiv:1510.04170 [gr-qc] .
  * Cunha  _et al._ (2015) P. V. P. Cunha, C. A. R. Herdeiro, E. Radu, and H. F. Runarsson, Phys. Rev. Lett. 115, 211102 (2015), arXiv:1509.00021 [gr-qc] .
  * Cunha  _et al._ (2016) P. V. P. Cunha, J. Grover, C. Herdeiro, E. Radu, H. Runarsson, and A. Wittig, Phys. Rev. D 94, 104023 (2016), arXiv:1609.01340 [gr-qc] .
  * Cunha and Herdeiro (2018) P. V. P. Cunha and C. A. R. Herdeiro, Gen. Rel. Grav. 50, 42 (2018), arXiv:1801.00860 [gr-qc] .
  * Rosa and Rubiera-Garcia (2022) J. a. L. Rosa and D. Rubiera-Garcia, Phys. Rev. D 106, 084004 (2022), arXiv:2204.12949 [gr-qc] .
  * Torres  _et al._ (2000) D. F. Torres, S. Capozziello, and G. Lambiase, Phys. Rev. D 62, 104012 (2000), arXiv:astro-ph/0004064 .
  * Diemer  _et al._ (2013) V. Diemer, K. Eilers, B. Hartmann, I. Schaffer, and C. Toma, Phys. Rev. D 88, 044025 (2013), arXiv:1304.5646 [gr-qc] .
  * Brihaye  _et al._ (2014) Y. Brihaye, V. Diemer, and B. Hartmann, Phys. Rev. D 89, 084048 (2014), arXiv:1402.1055 [gr-qc] .
  * Grandclement  _et al._ (2014) P. Grandclement, C. Somé, and E. Gourgoulhon, Phys. Rev. D 90, 024068 (2014), arXiv:1405.4837 [gr-qc] .
  * Grould  _et al._ (2017) M. Grould, Z. Meliani, F. H. Vincent, P. Grandclément, and E. Gourgoulhon, Class. Quant. Grav. 34, 215007 (2017), arXiv:1709.05938 [astro-ph.HE] .
  * Zhang  _et al._ (2022) Y.-P. Zhang, Y.-B. Zeng, Y.-Q. Wang, S.-W. Wei, and Y.-X. Liu, Phys. Rev. D 105, 044021 (2022), arXiv:2107.04848 [gr-qc] .
  * Delgado  _et al._ (2023) J. F. M. Delgado, C. A. R. Herdeiro, and E. Radu, (2023), arXiv:2305.02333 [gr-qc] .
  * Herdeiro and Radu (2014) C. A. R. Herdeiro and E. Radu, Phys. Rev. Lett. 112, 221101 (2014), arXiv:1403.2757 [gr-qc] .
  * Collodel  _et al._ (2022) L. G. Collodel, D. D. Doneva, and S. S. Yazadjiev, Phys. Rev. D 105, 044036 (2022), arXiv:2108.11658 [gr-qc] .
  * Zhang  _et al._ (2021) Y.-P. Zhang, Y.-B. Zeng, Y.-Q. Wang, S.-W. Wei, P. A. Seoane, and Y.-X. Liu, (2021), arXiv:2108.13170 [gr-qc] .
  * Ryan (1995) F. D. Ryan, Phys. Rev. D 52, 5707 (1995).
  * Ryan (1997b) F. D. Ryan, Phys. Rev. D 56, 1845 (1997b).
  * Ryan (1997c) F. D. Ryan, Phys. Rev. D 56, 7732 (1997c).
  * Vaglio  _et al._ (2022) M. Vaglio, C. Pacilio, A. Maselli, and P. Pani, Phys. Rev. D 105, 124020 (2022), arXiv:2203.07442 [gr-qc] .
  * Apostolatos  _et al._ (2009) T. A. Apostolatos, G. Lukes-Gerakopoulos, and G. Contopoulos, Phys. Rev. Lett. 103, 111101 (2009), arXiv:0906.0093 [gr-qc] .
  * Lukes-Gerakopoulos  _et al._ (2010) G. Lukes-Gerakopoulos, T. A. Apostolatos, and G. Contopoulos, Phys. Rev. D 81, 124005 (2010), arXiv:1003.3120 [gr-qc] .
  * Contopoulos  _et al._ (2011) G. Contopoulos, G. Lukes-Gerakopoulos, and T. A. Apostolatos, Int. J. Bifurc. Chaos 21, 2261 (2011), arXiv:1108.5057 [gr-qc] .
  * Lukes-Gerakopoulos (2012) G. Lukes-Gerakopoulos, Phys. Rev. D 86, 044013 (2012), arXiv:1206.0660 [gr-qc] .
  * Zelenka and Lukes-Gerakopoulos (2017) O. Zelenka and G. Lukes-Gerakopoulos, in  _Workshop on Black Holes and Neutron Stars_ (2017) arXiv:1711.02442 [gr-qc] .
  * Zelenka  _et al._ (2020) O. Zelenka, G. Lukes-Gerakopoulos, V. Witzany, and O. Kopáček, Phys. Rev. D 101, 024037 (2020), arXiv:1911.00414 [gr-qc] .
  * Lukes-Gerakopoulos and Witzany (2021) G. Lukes-Gerakopoulos and V. Witzany,  (2021), 10.1007/978-981-15-4702-7-42-1, arXiv:2103.06724 [gr-qc] .
  * Destounis  _et al._ (2020) K. Destounis, A. G. Suvorov, and K. D. Kokkotas, Phys. Rev. D 102, 064041 (2020), arXiv:2009.00028 [gr-qc] .
  * Deich  _et al._ (2022) A. Deich, A. Cárdenas-Avendaño, and N. Yunes, Phys. Rev. D 106, 024040 (2022), arXiv:2203.00524 [gr-qc] .
  * Chen  _et al._ (2022) C.-Y. Chen, F.-L. Lin, and A. Patel, Phys. Rev. D 106, 084064 (2022), arXiv:2206.10966 [gr-qc] .
  * Destounis  _et al._ (2021) K. Destounis, A. G. Suvorov, and K. D. Kokkotas, Phys. Rev. Lett. 126, 141102 (2021), arXiv:2103.05643 [gr-qc] .
  * Destounis and Kokkotas (2021) K. Destounis and K. D. Kokkotas, Phys. Rev. D 104, 064023 (2021), arXiv:2108.02782 [gr-qc] .
  * Destounis  _et al._ (2023a) K. Destounis, G. Huez, and K. D. Kokkotas, (2023a), arXiv:2301.11483 [gr-qc] .
  * Glampedakis (2005) K. Glampedakis, Class. Quant. Grav. 22, S605 (2005), arXiv:gr-qc/0509024 .
  * Amaro-Seoane  _et al._ (2017) P. Amaro-Seoane _et al._ (LISA), (2017), arXiv:1702.00786 [astro-ph.IM] .
  * Baibhav  _et al._ (2021) V. Baibhav _et al._ , Exper. Astron. 51, 1385 (2021), arXiv:1908.11390 [astro-ph.HE] .
  * Amaro-Seoane  _et al._ (2022) P. Amaro-Seoane _et al._ , (2022), arXiv:2203.06016 [gr-qc] .
  * Arun  _et al._ (2022) K. G. Arun _et al._ (LISA), Living Rev. Rel. 25, 4 (2022), arXiv:2205.01597 [gr-qc] .
  * Karnesis  _et al._ (2022) N. Karnesis _et al._ , (2022), arXiv:2209.04358 [gr-qc] .
  * Luo  _et al._ (2016) J. Luo _et al._ (TianQin), Class. Quant. Grav. 33, 035010 (2016), arXiv:1512.02076 [astro-ph.IM] .
  * Ruan  _et al._ (2020a) W.-H. Ruan, Z.-K. Guo, R.-G. Cai, and Y.-Z. Zhang, Int. J. Mod. Phys. A 35, 2050075 (2020a), arXiv:1807.09495 [gr-qc] .
  * Ruan  _et al._ (2020b) W.-H. Ruan, C. Liu, Z.-K. Guo, Y.-L. Wu, and R.-G. Cai, Nature Astron. 4, 108 (2020b), arXiv:2002.03603 [gr-qc] .
  * Barsanti  _et al._ (2022) S. Barsanti, V. De Luca, A. Maselli, and P. Pani, Phys. Rev. Lett. 128, 111104 (2022), arXiv:2109.02170 [gr-qc] .
  * Boyer and Lindquist (1967) R. H. Boyer and R. W. Lindquist, J. Math. Phys. 8, 265 (1967).
  * Misner  _et al._ (1973) C. W. Misner, K. S. Thorne, and J. A. Wheeler,  _Gravitation_ (W. H. Freeman, San Francisco, 1973).
  * Carter (1968) B. Carter, Phys. Rev. 174, 1559 (1968).
  * Contopoulos (2003) G. Contopoulos,  _Order and Chaos in Dynamical Astronomy_ (Springer-Verlag, New York, 2003).
  * Zipoy (1966) D. M. Zipoy, Journal of Mathematical Physics 7, 1137 (1966).
  * Voorhees (1970) B. H. Voorhees, Phys. Rev. D 2, 2119 (1970).
  * Polcar  _et al._ (2022) L. Polcar, G. Lukes-Gerakopoulos, and V. Witzany, Phys. Rev. D 106, 044069 (2022), arXiv:2205.08516 [gr-qc] .
  * Hartle and Thorne (1968) J. B. Hartle and K. S. Thorne, Astrophys. J. 153, 807 (1968).
  * Manko  _et al._ (2000a) V. S. Manko, E. W. Mielke, and J. D. Sanabria-Gomez, Phys. Rev. D 61, 081501 (2000a), arXiv:gr-qc/0001081 .
  * Pappas and Apostolatos (2013) G. Pappas and T. A. Apostolatos, Mon. Not. Roy. Astron. Soc. 429, 3007 (2013), arXiv:1209.6148 [gr-qc] .
  * Pappas (2017) G. Pappas, Mon. Not. Roy. Astron. Soc. 466, 4381 (2017), arXiv:1610.05370 [gr-qc] .
  * Manko and Novikov (1992) V. S. Manko and I. D. Novikov, Classical and Quantum Gravity 9, 2477 (1992).
  * Manko  _et al._ (2000b) V. S. Manko, J. D. Sanabria-Gomez, and O. V. Manko, Phys. Rev. D 62, 044048 (2000b).
  * Collins and Hughes (2004) N. A. Collins and S. A. Hughes, Phys. Rev. D 69, 124022 (2004), arXiv:gr-qc/0402063 .
  * Glampedakis and Babak (2006) K. Glampedakis and S. Babak, Class. Quant. Grav. 23, 4167 (2006), arXiv:gr-qc/0510057 .
  * Vigeland and Hughes (2010) S. J. Vigeland and S. A. Hughes, Phys. Rev. D 81, 024030 (2010), arXiv:0911.1756 [gr-qc] .
  * Johannsen and Psaltis (2011) T. Johannsen and D. Psaltis, Phys. Rev. D 83, 124015 (2011), arXiv:1105.3191 [gr-qc] .
  * Vigeland  _et al._ (2011) S. Vigeland, N. Yunes, and L. Stein, Phys. Rev. D 83, 104027 (2011), arXiv:1102.3706 [gr-qc] .
  * Johannsen (2013a) T. Johannsen, J. Phys. Conf. Ser. 410, 012136 (2013a), arXiv:1210.0483 [gr-qc] .
  * Johannsen (2013b) T. Johannsen, Phys. Rev. D 88, 044002 (2013b), arXiv:1501.02809 [gr-qc] .
  * Emparan  _et al._ (2014) R. Emparan, P. Figueras, and M. Martínez, Journal of High Energy Physics 2014, 72 (2014), arXiv:1410.4764 [hep-th] .
  * Cardoso  _et al._ (2014) V. Cardoso, P. Pani, and J. Rico, Phys. Rev. D 89, 064007 (2014), arXiv:1401.0528 [gr-qc] .
  * Rezzolla and Zhidenko (2014) L. Rezzolla and A. Zhidenko, Phys. Rev. D 90, 084009 (2014), arXiv:1407.3086 [gr-qc] .
  * Konoplya  _et al._ (2016) R. Konoplya, L. Rezzolla, and A. Zhidenko, Phys. Rev. D 93, 064015 (2016), arXiv:1602.02378 [gr-qc] .
  * Moore  _et al._ (2017) C. J. Moore, A. J. K. Chua, and J. R. Gair, Class. Quant. Grav. 34, 195009 (2017), arXiv:1707.00712 [gr-qc] .
  * Flanagan and Hinderer (2012) E. E. Flanagan and T. Hinderer, Phys. Rev. Lett. 109, 071102 (2012), arXiv:1009.4923 [gr-qc] .
  * Berry  _et al._ (2016) C. P. L. Berry, R. H. Cole, P. Cañizares, and J. R. Gair, Phys. Rev. D 94, 124042 (2016), arXiv:1608.08951 [gr-qc] .
  * Speri and Gair (2021) L. Speri and J. R. Gair, Phys. Rev. D 103, 124032 (2021), arXiv:2103.06306 [gr-qc] .
  * Gupta  _et al._ (2022) P. Gupta, L. Speri, B. Bonga, A. J. K. Chua, and T. Tanaka, Phys. Rev. D 106, 104001 (2022), arXiv:2205.04808 [gr-qc] .
  * Möser (1962) J. Möser, Nachr. Akad. Wiss. Göttingen, II , 1 (1962).
  * Arnol’d (1963) V. I. Arnol’d, Russian Mathematical Surveys 18, 9 (1963).
  * Birkhoff (1913) G. D. Birkhoff, Transactions of the American Mathematical Society 14, 14 (1913).
  * Ruangsri and Hughes (2014) U. Ruangsri and S. A. Hughes, Phys. Rev. D 89, 084036 (2014), arXiv:1307.6483 [gr-qc] .
  * Barack (2009) L. Barack, Class. Quant. Grav. 26, 213001 (2009), arXiv:0908.1664 [gr-qc] .
  * Schunck and Mielke (1996) F. E. Schunck and E. W. Mielke, “Rotating boson stars,” in  _Relativity and Scientific Computing: Computer Algebra, Numerics, Visualization_, edited by F. W. Hehl, R. A. Puntigam, and H. Ruder (Springer Berlin Heidelberg, Berlin, Heidelberg, 1996) pp. 138–151.
  * Yoshida and Eriguchi (1997) S. Yoshida and Y. Eriguchi, Phys. Rev. D 56, 762 (1997).
  * Amaro-Seoane  _et al._ (2010) P. Amaro-Seoane, J. Barranco, A. Bernal, and L. Rezzolla, JCAP 11, 002 (2010), arXiv:1009.0019 [astro-ph.CO] .
  * Brito  _et al._ (2016) R. Brito, V. Cardoso, C. A. R. Herdeiro, and E. Radu, Phys. Lett. B 752, 291 (2016), arXiv:1508.05395 [gr-qc] .
  * Herdeiro and Radu (2015) C. Herdeiro and E. Radu, Class. Quant. Grav. 32, 144001 (2015), arXiv:1501.04319 [gr-qc] .
  * Cunha  _et al._ (2023) P. V. P. Cunha, C. Herdeiro, E. Radu, and N. Sanchis-Gual, Phys. Rev. Lett. 130, 061401 (2023), arXiv:2207.13713 [gr-qc] .
  * Franchini  _et al._ (2017) N. Franchini, P. Pani, A. Maselli, L. Gualtieri, C. A. R. Herdeiro, E. Radu, and V. Ferrari, Phys. Rev. D 95, 124025 (2017), arXiv:1612.00038 [astro-ph.HE] .
  * Pound and Wardell (2021) A. Pound and B. Wardell,  (2021), 10.1007/978-981-15-4702-7_38-1, arXiv:2101.04592 [gr-qc] .
  * van de Meent (2016) M. van de Meent, Phys. Rev. D 94, 044034 (2016), arXiv:1606.06297 [gr-qc] .
  * van de Meent (2018) M. van de Meent, Phys. Rev. D 97, 104033 (2018), arXiv:1711.09607 [gr-qc] .
  * Pound (2017) A. Pound, Phys. Rev. D 95, 104056 (2017), arXiv:1703.02836 [gr-qc] .
  * Teukolsky (1973) S. A. Teukolsky, Astrophys. J. 185, 635 (1973).
  * Teukolsky and Press (1974) S. A. Teukolsky and W. H. Press, Astrophys. J. 193, 443 (1974).
  * Detweiler (1978) S. L. Detweiler, Astrophys. J. 225, 687 (1978).
  * Sasaki and Nakamura (1982) M. Sasaki and T. Nakamura, Prog. Theor. Phys. 67, 1788 (1982).
  * Mino  _et al._ (1997) Y. Mino, M. Sasaki, M. Shibata, H. Tagoshi, and T. Tanaka, Prog. Theor. Phys. Suppl. 128, 1 (1997), arXiv:gr-qc/9712057 .
  * Hughes (2000) S. A. Hughes, Phys. Rev. D 61, 084004 (2000), [Erratum: Phys.Rev.D 63, 049902 (2001), Erratum: Phys.Rev.D 65, 069902 (2002), Erratum: Phys.Rev.D 67, 089901 (2003), Erratum: Phys.Rev.D 78, 109902 (2008), Erratum: Phys.Rev.D 90, 109904 (2014)], arXiv:gr-qc/9910091 .
  * Hughes (2001) S. A. Hughes, Phys. Rev. D 64, 064004 (2001), [Erratum: Phys.Rev.D 88, 109902 (2013)], arXiv:gr-qc/0104041 .
  * Drasco and Hughes (2006) S. Drasco and S. A. Hughes, Phys. Rev. D 73, 024027 (2006), [Erratum: Phys.Rev.D 88, 109905 (2013), Erratum: Phys.Rev.D 90, 109905 (2014)], arXiv:gr-qc/0509101 .
  * Poisson (2004) E. Poisson, Phys. Rev. D 70, 084044 (2004), arXiv:gr-qc/0407050 .
  * Martel (2004) K. Martel, Phys. Rev. D 69, 044025 (2004), arXiv:gr-qc/0311017 .
  * Lopez-Aleman  _et al._ (2003) R. Lopez-Aleman, G. Khanna, and J. Pullin, Class. Quant. Grav. 20, 3259 (2003), arXiv:gr-qc/0303054 .
  * Khanna (2004) G. Khanna, Phys. Rev. D 69, 024016 (2004), arXiv:gr-qc/0309107 .
  * Pound (2012) A. Pound, Phys. Rev. Lett. 109, 051101 (2012), arXiv:1201.5089 [gr-qc] .
  * Pound  _et al._ (2020) A. Pound, B. Wardell, N. Warburton, and J. Miller, Phys. Rev. Lett. 124, 021101 (2020), arXiv:1908.07419 [gr-qc] .
  * (183) A. Spiers,  10.5258/SOTON/T0051.
  * Glampedakis  _et al._ (2002) K. Glampedakis, S. A. Hughes, and D. Kennefick, Phys. Rev. D 66, 064005 (2002), arXiv:gr-qc/0205033 .
  * Gair and Glampedakis (2006) J. R. Gair and K. Glampedakis, Phys. Rev. D 73, 064037 (2006), arXiv:gr-qc/0510129 .
  * Ruffini and Sasaki (1981) R. Ruffini and M. Sasaki, Progress of Theoretical Physics 66, 1627 (1981).
  * Sasaki (1983) M. Sasaki, Progress of Theoretical Physics 69, 815 (1983).
  * Babak  _et al._ (2007) S. Babak, H. Fang, J. R. Gair, K. Glampedakis, and S. A. Hughes, Phys. Rev. D 75, 024005 (2007), [Erratum: Phys.Rev.D 77, 04990 (2008)], arXiv:gr-qc/0607007 .
  * Glampedakis and Kennefick (2002) K. Glampedakis and D. Kennefick, Phys. Rev. D 66, 044002 (2002), arXiv:gr-qc/0203086 .
  * Barack and Cutler (2007) L. Barack and C. Cutler, Phys. Rev. D 75, 042003 (2007), arXiv:gr-qc/0612029 .
  * Gair  _et al._ (2008) J. R. Gair, C. Li, and I. Mandel, Phys. Rev. D 77, 024035 (2008), arXiv:0708.0628 [gr-qc] .
  * Mukherjee  _et al._ (2023) S. Mukherjee, O. Kopacek, and G. Lukes-Gerakopoulos, Phys. Rev. D 107, 064005 (2023), arXiv:2206.10302 [gr-qc] .
  * Canizares  _et al._ (2012) P. Canizares, J. R. Gair, and C. F. Sopuerta, Phys. Rev. D 86, 044010 (2012), arXiv:1205.1253 [gr-qc] .
  * Cutler (1998) C. Cutler, Phys. Rev. D 57, 7089 (1998), arXiv:gr-qc/9703068 .
  * Barack and Cutler (2004) L. Barack and C. Cutler, Phys. Rev. D 69, 082005 (2004), arXiv:gr-qc/0310125 .
  * van de Meent (2014) M. van de Meent, Phys. Rev. D 89, 084033 (2014), arXiv:1311.4457 [gr-qc] .
  * Kiuchi and Maeda (2004) K. Kiuchi and K.-i. Maeda, Phys. Rev. D 70, 064036 (2004), arXiv:gr-qc/0404124 .
  * Piovano  _et al._ (2020) G. A. Piovano, A. Maselli, and P. Pani, Phys. Rev. D 102, 024041 (2020), arXiv:2004.02654 [gr-qc] .
  * Cardoso  _et al._ (2022a) V. Cardoso, K. Destounis, F. Duque, R. P. Macedo, and A. Maselli, Phys. Rev. D 105, L061501 (2022a), arXiv:2109.00005 [gr-qc] .
  * Cardoso  _et al._ (2022b) V. Cardoso, K. Destounis, F. Duque, R. Panosso Macedo, and A. Maselli, Phys. Rev. Lett. 129, 241103 (2022b), arXiv:2210.01133 [gr-qc] .
  * Destounis  _et al._ (2023b) K. Destounis, A. Kulathingal, K. D. Kokkotas, and G. O. Papadopoulos, Phys. Rev. D 107, 084027 (2023b), arXiv:2210.09357 [gr-qc] .
  * Bena  _et al._ (2022) I. Bena, E. J. Martinec, S. D. Mathur, and N. P. Warner, (2022), arXiv:2203.04981 [hep-th] .
  * Bena and Mayerson (2020) I. Bena and D. R. Mayerson, Phys. Rev. Lett. 125, 221602 (2020), arXiv:2006.10750 [hep-th] .
  * Bianchi  _et al._ (2020) M. Bianchi, D. Consoli, A. Grillo, J. F. Morales, P. Pani, and G. Raposo, Phys. Rev. Lett. 125, 221601 (2020), arXiv:2007.01743 [hep-th] .
  * Bena and Mayerson (2021) I. Bena and D. R. Mayerson, JHEP 03, 114 (2021), arXiv:2007.09152 [hep-th] .
  * Bianchi  _et al._ (2021) M. Bianchi, D. Consoli, A. Grillo, J. F. Morales, P. Pani, and G. Raposo, JHEP 01, 003 (2021), arXiv:2008.01445 [hep-th] .
  * Fransen and Mayerson (2022) K. Fransen and D. R. Mayerson, Phys. Rev. D 106, 064035 (2022), arXiv:2201.03569 [gr-qc] .