# Slowly-rotating compact objects: the nonintegrability of Hartle-Thorne
particle geodesics

Kyriakos Destounis  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<EMAIL_ADDRESS>and Kostas D.
Kokkotas  Theoretical Astrophysics, IAAT, University of Tübingen, 72076
Tübingen, Germany  Section of Astrophysics, Astronomy, and Mechanics,
Department of Physics, National and Kapodistrian University of Athens,
Panepistimiopolis Zografos GR15783, Athens, Greece

###### Abstract

X-ray astronomy provides information regarding the electromagnetic emission of
active galactic nuclei and X-ray binaries. These events provide details
regarding the astrophysical environment of black holes and stars, and help us
understand gamma-ray bursts. They produce estimates for the maximum mass of
neutron stars and eventually will contribute to the discovery of their
equation of state. Thus, it is crucial to study these configurations to
increase the yield of X-ray astronomy when combined with multimessenger
gravitational-wave astrophysics and black hole shadows. Unfortunately, an
exact solution of the field equations does not exist for neutron stars.
Nevertheless, there exist a variety of approximate compact objects that may
characterize massive or neutron stars. The most studied approximation is the
Hartle-Thorne metric that represents slowly-rotating compact objects, like
massive stars, white dwarfs and neutron stars. Recent investigations of photon
orbits and shadows of such metric revealed that it exhibits chaos close to
resonances. Here, we thoroughly investigate particle orbits around the Hartle-
Thorne spacetime. We perform an exhaustive analysis of bound motion, by
varying all parameters involved in the system. We demonstrate that chaotic
regions, known as Birkhoff islands, form around resonances, where the ratio of
the radial and polar frequency of geodesics, known as the rotation number, is
shared throughout the island. This leads to the formation of plateaus in
rotation curves during the most prominent $2/3$ resonance, which designate
nonintegrability. We measure their width and show how each parameter affects
it. The nonintegrability of Hartle-Thorne metric may affect quasiperiodic
oscillations of low-mass X-ray binaries, when chaos is taken into account, and
improve estimates of mass, angular momentum and multipole moments of
astrophysical compact objects.

##  I Introduction

X-ray observations have been providing astrophysical explanations for various
questions regarding the exterior of astrophysical compact objects such as
black holes (BHs) and neutron stars (NSs), particularly related to their
surrounding astrophysical environments, including active galactic nuclei and
low-mass X-ray binaries [1, 2, 3, 4].

The motion of matter around compact objects can be influenced by
electromagnetic fields, which originate from a NS’s boundary layer, an X-ray
burst on the NS’s polar caps, a hot corona around a BH, or an accretion disk
in general (as described by the Poynting-Robertson effect; see, for example,
[5, 6]). Observations from the Event Horizon Telescope that probe the
environmental vicinity and shadow of M87* [7] and Sgr A* [8] and
gravitational-wave (GW) observations from the LIGO/Virgo/KAGRA collaboration
[9, 10] of BH-NS and NS binaries [11, 12] have paved the way for multi-
messenger astrophysics. For binaries, this approach can combine both
electromagnetic and gravitational knowledge [13]. By incorporating the
experiments mentioned above, there are endless possibilities to explore new
fundamental physics, gain a deeper understanding of the composition of
astrophysical environments surrounding compact objects, and investigate the
effects of dark matter halos in galaxies [14, 15, 16, 17, 18].

The binary NS merger GW170817 [11], was followed by a gamma-ray burst and a
kilonova [19]. This observation was made through the use of electromagnetic
telescopes across the spectra, and it confirmed that binary NSs are the
progenitors of short gamma-ray bursts. Furthermore, it constrained the
difference between the speed of light and the propagation speed of GWs,
established new bounds on the violation of Lorentz invariance, and presented a
new test for the equivalence principle.

In general, it is difficult to find any astrophysical object or system that is
not spinning. As a result, they are not spherical and have specific multipolar
structures. The degree of deviation from sphericity depends on their spin,
environment, and other possible modifications to their composition. If the
object is a spinning BH, then the spacetime is fully described by an exact
solution of the Einstein field equations, the Kerr metric [20]. In general,
things are quite simple when dealing with Kerr BHs, which are described only
by their mass and spin. On the other hand, compact rotating NSs are much more
complicated and cannot be efficiently approximated by a Kerr exterior [21]. As
of today, there is no analytical description of such spacetimes, and the
geometry surrounding a rotating NS is much more complex because it directly
depends on the star’s internal structure and its spin. Hence, the assumption
that the exterior of a NS is similar to the Kerr metric is excessively
simplistic and will likely lead to inaccurate conclusions regarding
astrophysical processes occurring in their close proximity [22].

In order to achieve this goal, the scientific community has invested a
significant effort in finding more precise approximations of NS exteriors, as
evidenced by studies such as those referenced in [23, 24, 25, 26, 27, 28, 29,
30, 31], with the eventual aim of using numerical methods to construct
relativistic, highly accurate, and realistic NSs [32, 33]. However, despite
the advantages of numerical solutions of the field equations, such as
increased reliability and adaptability to different symmetries and equation of
state assumptions, analytic expressions remain preferable.

The first and most extensively tested approximate solution to the field
equations for the exterior of slowly-rotating, massive stars, white dwarfs or
NSs, has been constructed perturbatively by Hartle and Thorne (HT) [34, 35,
36, 37, 38, 39]. This solution incorporates the Lense-Thirring effect [40,
41], exhibits frame dragging, leads to important precession effects [42] and
thus it has been specifically used to infer several astrophysical phenomena.

Quasiperiodic oscillations, observed as flickering in the kHz X-ray luminosity
spectra of BHs and NSs at specific frequencies, assist astronomers in
comprehending the innermost regions of accretion disks [43, 44, 45].
Additionally, they provide estimations of the masses, radii, and spin periods
of white dwarfs, BHs, and NSs. These astrophysical investigations include an
accreting rotating NS, which is described by the HT metric that engulfs matter
from its donor star. This, in turn, aids in further understanding the
frequency properties of radial and vertical epicyclic motion [46], as well as
geodesic motion [42, 47]. However, potential chaotic phenomena resulting from
the altered multipolar structure of these objects, compared to Kerr, have not
been examined yet, neither the effects of such phenomena on geodesic motion,
particle frequencies, and the accretion process.

The Kerr metric, which is currently providing a good description of the shadow
of M87* as measured by the Event Horizon Telescope, has integrable geodesics.
This means that the motion of particles and photons around it can be separated
into four first-order differential equations with the help of the Carter
constant [48]. Integrability ensures a deterministic evolution of the
equations of motion and the absence of chaos in null and timelike geodesics.
However, compact objects that are not Kerr BHs have a different multipolar
structure, which often leads to the absence of a Carter-like constant [49,
50]. This may be occur in alternative theories of gravity [51, 52, 53], in
non-vacuum BHs surrounded by astrophysical environments [54], in many body
systems [55], and in asymmetric binaries where the stellar-mass spinning
secondary orbits a supermassive BH [56, 57, 58]. As a result, a plethora of
direct and indirect imprints of chaos occur. One recent example is the
realization of chaos in the shadows of highly deformed spacetimes, such as the
Hartle-Thorne metric, as discussed in [59]. These spacetimes form photon-
trapping pockets that give rise to chaos, where photons may become trapped
near an orbital resonance between their radial and polar oscillations. This
prolongs the trapping time and can produce chaotic fractal features in the
corresponding shadows.

In this paper, we explore a significant portion of the parameter space of
geodesics for massive particles around a HT compact object, up to second-order
in rotation, denoted as $\mathcal{O}(\Omega^{2})$ [35, 34]. We do not restrict
our study to circular or equatorial orbits, which have been examined in
previous works. Instead, we investigate generic orbits to search for chaotic
phenomena. Our findings confirm the existence of indirect chaos, such as the
formation of Birkhoff chains in phase space [60]. Additionally, we observe
islands of stability around periodic orbits and plateaus in the geodesic’s
radial ($\omega_{r}$) and polar ($\omega_{\theta}$) oscillation frequencies,
which define the rotation curve. We analyzed the dependence of the width of
resonant islands with respect to the $2/3$ resonance, which one of the most
prominent ones in extreme-mass-ratio inspirals [61, 62], as well as the
particle’s energy, $z$-component of angular momentum, the spacetime’s spin and
mass quadrupole deformation.

We observe that fixing all parameters and then varying the particle’s energy
(or $z$-component of its angular momentum) results in wider (or narrower)
$2/3$ resonant islands and corresponding plateaus in the rotation curves.
Additionally, we have found that increasing the spin and mass quadrupole
deformation generally leads to wider islands, which agrees with previous
research on non-Kerr metrics [63, 64, 65, 58, 66], parameterized metrics [67,
68, 69], and rotating boson stars [70]. Accounting for these chaotic effects
could provide further insight into the astrophysical environment of low-mass
X-ray binaries, improve mass, spin, and multipole moment estimates of the BH
or NS that consumes the donor star, and aid in the search for an ultimate
equation of state for the supranuclear matter. We use geometric units in the
following discussion where the gravitational constant $G$ and the speed of
light $c$ are set $G=c=1$.

##  II The Hartle-Thorne metric

The HT metric describes stationary and axisymmetric spacetimes and can be used
to approximately describe the interior and exterior of a rigidly, slowly-
rotating compact object, such as a massive star, a white dwarf or a NS [35,
34]. It is not an exact solution of field equations but rather it is
constructed perturbatively in terms of the rotation rate of the object. This
methodology is also used to construct other rotating compact configurations,
such as boson stars [71, 72, 73]. The explicit form of the HT metric is

$\displaystyle
ds^{2}=-e^{\nu}\left(1+2h\right)dt^{2}+e^{\lambda}\left(1+\frac{2\mu}{r-2m}\right)dt^{2}$
$\displaystyle+r^{2}\left(1+2k\right)\\{d\theta^{2}+\sin^{2}\theta\left[d\phi-\left(\Omega-\omega\right)dt\right]^{2}\\}+\mathcal{O}(\Omega^{3}),$
(1)

where $\Omega$ is the rotation rate of the configuration, while the metric
functions $h(r,\theta)$, $\mu(r,\theta)$, $k(r,\theta)$ and $\omega(r,\theta)$
can be expressed in terms of Legendre polynomials $P_{\ell}(\cos\theta)$ as

$\displaystyle h(r,\theta)$ $\displaystyle=h_{0}(r)+h_{2}(r)P_{2},$ (2)
$\displaystyle\mu(r,\theta)$ $\displaystyle=\nu_{0}(r)+\mu_{2}(r)P_{2},$ (3)
$\displaystyle k(r,\theta)$ $\displaystyle=k_{2}(r)P_{2},$ (4)
$\displaystyle\omega(r,\theta)$ $\displaystyle=\omega_{1}(r)P^{\prime}_{1}.$
(5)

Here $\chi=a/M=J/M^{2}$ is the dimensionless spin parameter with $a$ the spin,
$J$ the angular momentum and $M$ the mass of the non-rotating object.
Furthermore, we can define a quadrupole deviation $\delta q$ from the Kerr
quadrupole that leads to $Q=-\chi^{2}M^{3}(1-\delta q)$ [35, 34, 74, 75, 59].
The term $\left(1-\delta q\right)$ means that the larger $\delta q$ is, the
more prolate the object becomes along the direction of the spin axis, leading
to a more positive quadrupole than that of Kerr. In the special case when
$\delta q=1$, $Q=0$. Nevertheless, the object is not spherical due to its non-
zero rotation $\chi\Rightarrow g_{t\phi}\neq 0$. Stated differently, the
metric indicates that there is still a prolate morphology because of the
presence of a non-zero (positive) quadrupole deformation. The object’s spin,
together with the prolateness introduced by $\delta q=1$, counterbalance each
other in the mass quadrupole to give $Q=0$, but the object is still spinning,
which is a typical behavior of prolate rotating compact objects.

Furthermore, we define the second-order correction to the mass $\delta m$ due
to the deformation and spin, and a dimensional radial coordinate $x=r/M$ to
obtain the HT metric as follows [76, 77, 59]:

$\displaystyle m$
$\displaystyle=M,\,\,\,\,\,\,e^{\nu}=e^{-\lambda}=1-\frac{2}{x},\,\,\,\,\,\,\omega_{1}=\Omega-\frac{2\chi}{Mx^{3}},$
$\displaystyle\mu_{0}$ $\displaystyle=M\chi^{2}\left(\delta
m-x^{-3}\right),\,\,\,\,\,\,h_{0}=\frac{\chi^{2}}{x-2}\left(x^{-3}-\delta
m\right),$ $\displaystyle h_{2}$ $\displaystyle=\frac{5\chi^{2}}{16}\delta
q\left(1-\frac{2}{x}\right)\left[3x^{2}\log\left(1-2/x\right)\right.$
$\displaystyle\left.+\frac{(2-2/x)}{x(1-2/x)^{2}}(3x^{2}-6x-2)+\frac{\chi^{2}}{x^{3}}\left(1+\frac{1}{x}\right)\right],$
$\displaystyle k_{2}$
$\displaystyle=-\frac{\chi^{2}}{x^{3}}\left(1+\frac{1}{x}\right)-\frac{5\chi^{2}}{8}\delta
q$
$\displaystyle\left[3\left(1+x-\frac{2}{x}-3\left(1-\frac{x^{2}}{2}\right)\right)\log\left(1-\frac{2}{x}\right)\right],$
$\displaystyle\mu_{2}$ $\displaystyle=-\frac{5Mx\chi^{2}}{16}\delta
q\left(1-\frac{2}{x}\right)^{2}\left[3x^{2}\log\left(1-\frac{2}{x}\right)\right.$
$\displaystyle\left.+\frac{(2-2/x)}{x(1-2/x)^{2}}(3x^{2}-6x-2)-\frac{\chi^{2}}{x^{2}}\left(1-\frac{7}{x}+\frac{10}{x^{2}}\right)\right].$

For simplicity, we assume that the total mass of spacetime is $M\equiv
m+\delta m$ that includes the correction due to $\delta q$ and the small spin
$\chi=a/M$ and basically amounts to setting $\delta m=0$ in the aforementioned
equations. Lastly, the composition of matter that constructs the HT compact
object is assumed to be such that usual issues with photon emission are
evaded, so that the metric can be considered as a BH mimicker [78, 79]. The
metric described by (II) reduces to the approximate Kerr spacetime, up to
$\Omega^{2}$, when $\delta q=0$ and an appropriate transformation is imposed
[75, 77, 80]. It also possesses an ergosphere when $g_{tt}=0$ [42] while when
the spin is set to zero, it reduces to Schwarzschild, since $Q=0$,
$\omega_{1}=\Omega$ and $h=\mu=k=\omega=0$ in (II). For more details regarding
regions with potential pathological features of the HT metric (which we avoid
in this work) see [59].

##  III Geodesics and chaotic imprints

The motion of test-particles on a general, stationary and axisymmetric
spacetime of the form

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

is described by the geodesic equations

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

where $\Gamma^{\kappa}_{\lambda\nu}$ are the Christoffel symbols, $x^{\kappa}$
is the four-position of the timelike orbit, and the overdots denote
differentiation with respect to proper time $\tau$111Even though we will not
be dealing with GW emission in this work, one can tranform from proper to
coordinate time $t$ with a simple chain rule such that
$r^{\prime}=dr/dt=(dr/d\tau)(d\tau/dt)=\dot{r}/\dot{t}$.. From stationarity
and axisymmetry, the spacetime axiomatically possesses two Killing vector
fields, namely the conserved energy $E$ and $z$-component of the angular
momentum $L_{z}$, described through Eq. (7), as

$E/\mu=g_{tt}\dot{t}+g_{t\phi}\dot{\phi},\,\,\,\,\,\,\,\,\,L_{z}/\mu=g_{t\phi}\dot{t}+g_{\phi\phi}\dot{\phi},$
(8)

where $\mu$ is the conserved mass of the test particle (not to be confused
with the functions of spacetime $\mu(r,\theta)$). From Eq. (8), we obtain two
first-order differential equations for the $t$ and $\phi$ momenta as

$\dot{t}=\frac{Eg_{\phi\phi}+L_{z}g_{t\phi}}{\mu\left(g^{2}_{t\phi}-g_{tt}g_{\phi\phi}\right)},\,\,\,\,\,\,\,\,\,\dot{\phi}=\frac{Eg_{t\phi}+L_{z}g_{tt}}{\mu\left(g_{tt}g_{\phi\phi}-g^{2}_{t\phi}\right)}.$
(9)

Test particles in geodesic motion conserve, besides $E$ and $L_{z}$, their
mass (or equivalently their four-velocity
$g_{\lambda\nu}\dot{x}^{\lambda}\dot{x}^{\nu}=-1$), which leads to a
constraint equation that defines bound motion:

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

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

$V_{\text{eff}}\equiv\frac{1}{g_{rr}}\left(1+\frac{g_{\phi\phi}E^{2}+g_{tt}L_{z}^{2}+2g_{t\phi}EL_{z}}{\mu^{2}\left(g_{tt}g_{\phi\phi}-g_{t\phi}^{2}\right)}\right).$
(11)

When $V_{\text{eff}}=0$ the resulting curve is called curve of zero velocity
(CZV) since $\dot{r}=\dot{\theta}=0$ need to be satisfied. In general, the
radial and polar motion is described by second-order coupled differential
equations which are solved simultaneously. In the special case where there
exists a fourth constant of motion, such as the Carter constant in Kerr
spacetime, the system is characterized as integrable. All four equations for
the motion of $(t,r,\theta,\phi)$ are described by first-order, decoupled
equations, rendering the system deterministic. However, integrability is a
fragile property and can be easily broken.

Figure 1: Top row: Equatorial effective potential $V_{\text{eff}}$ and
corresponding CZV of the HT metric for varying $E/\mu$ and fixed
$L_{z}/\mu=3.5M$, $\chi=0.3$, $\delta q=1$. Bottom row: Equatorial effective
potential $V_{\text{eff}}$ and corresponding CZV for varying $L_{z}/\mu$ and
fixed $E/\mu=0.97$, $\chi=0.3$, $\delta q=1$.

Nonintegrability implies that there might be mild or strong chaotic phenomena
in the system’s evolution depending on its parameters [81, 60]. In these
cases, the Kolmogorov-Arnold-Moser (KAM) theorem [82, 83] states that _if the
orbits lie sufficiently away[60] from periodic points then its evolution is
slightly shifted from that of the corresponding integrable system but presents
no chaotic behavior_. For these orbits, curves that define successive
intersections on a surface of section called a Poincaré map [60], namely the
KAM curves, form closed curves and surround a periodic central point of the
Poincaré map. _On the other hand, the Poincaré-Birkhoff theorem[84] dictates
that the structure of a resonant KAM curve is disintegrated into an even
number of stable and unstable periodic points that form a Birkhoff chain when
the orbit is periodic_. Around stable periodic points, a series of resonant
KAM curves form, defining Birkhoff islands of stability (resonant islands).
These KAM curves do not surround the central point of the map but rather the
stable periodic points of the resonance. The unstable periodic orbits give
rise to chaotic orbits that shield the stability islands with thin chaotic
layers. The structure, as a whole, discussed above for resonant orbits of
nonintegrable systems is called a Birkhoff chain.

The most important aspect of resonant islands is the fact that particle orbits
that reside inside them share the same periodicity with the central periodic
orbit of the island, namely the ratio of the radial over the polar frequencies
$\omega_{r}/\omega_{\theta}$, or the rotation number. Integrable systems do
not exhibit Birkhoff chains, but still have resonant points. The difference
between a resonant island of a nonintegrable system and a resonant orbit of an
integrable system lies to the fact that the integrable resonance occupies a
zero-volume of phase space while the non-integrable resonant island induces a
finite volume of phase space where the rotation number is shared for all
geodesics inside the island and the trajectory is locked in perfect resonance.
This gives rise to an indirect imprint of chaos (due to the underlying but
mostly indistinguishable chaotic layers, though see [70]) and have been found
to affect the orbital evolution [63, 64] and GW emission [67, 68, 69, 66, 70]
from test-particles orbiting around nonintegrable spacetimes.

Figure 2: Left: Poincaré map of geodesics in HT spacetime with $E/\mu=0.97$,
$L_{z}/\mu=3.5M$, $\chi=0.3$, $\delta q=1$. The black curves correspond to
typical KAM curves while the red ones belong to the the resonant KAM curves
that surround the period stable point. Right: Zoom in to the leftmost island
from Left. Figure 3: Rotation curve of geodesics in HT metric from Fig. 1. The
inset zooms in the region where the $2/3$ plateau is clearly visible.

The two tools we will utilize in this work to conclude if chaos exists also in
particle orbits, besides photon ones [59], are the Poincaré map and rotation
curve [60]. Poincaré maps are formally defined as an accumulation of KAM
curves for successive geodesics, by fixing all parameters and varying one of
them, such as $r$ or $\theta$. The definition explicitly requires that the
intersection through the chosen surface of section to be one directional and
perpendicular to the surface. By sketching a high resolution Poincaré map, one
can find resonant islands and Birkhoff chains depending on the system and the
deformation from Kerr used. In HT spacetime, which we will operate, the
deformation parameter is $\delta q$, but nevertheless we will also explore the
dependence of the islands width with respect to all the spacetime parameters
and particle geodesics.

The rotation number is helpful to find regions where islands of stability
exist. We calculate it by tracking the angle $\vartheta$ between successive
intersections of KAM curves, relative to the fixed periodic central point of
the Poincaré map. The rotation number is then defined as the series of all
angles $\vartheta$ measured between one directional consecutive intersections,
i.e. [60]

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

where $N$ is the number of angles measured. When $N\rightarrow\infty$, Eq.
(12) asymptotes to $\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. The more angles measured, which corresponds
to more evolution time of the geodesic, the better we approximate the rotation
number. Empirical tests have shown that the rotation number’s accuracy is
inversely proportional to the number of intersections through the surface of
section, i.e. evolving an orbit for $10^{4}$ cycles, which leads to $\sim
10^{4}$ intersections, estimates the rotation number with order
$\mathcal{O}(10^{-4})$ accuracy (four decimal points). We stress that this is
an empirical law that seems to hold for gravitational systems and is not a
systematic rule for all possible scenarios.

Integrable systems have monotonous rotation curves, while nonintegrable ones
display discontinuities in the monotonicity through the formation of transient
plateaus with non-zero widths when geodesics traverse resonant islands (the
rotation number is shared throughout the island). Inflection points can also
appear when trajectories pass through unstable periodic points. Nevertheless,
by varying the initial conditions, we can drive the orbit through any island
and obtain a plateau.

##  IV Potential and Curves of Zero Velocity

Figure 4: Top left: Rotation curves of geodesics in HT metric with fixed
$L_{z}/\mu=3.5M$, $\delta q=1$, $\chi=0.3$ and varying $E/\mu$. The zoom-in
insets portray the plateaus appearing in $2/3$ islands for some of the cases.
Top right: Same us top left with fixed $E/\mu=0.97$, $\delta q=1$, $\chi=0.3$
and varying $L_{z}/\mu$. Bottom left: Same as top left with fixed
$E/\mu=0.97$, $L_{z}/\mu=3.5M$, $\chi=0.3$ and varying $\delta q$. Bottom
right: Same as top left with fixed $E/\mu=0.97$, $L_{z}/\mu=3.5M$, $\delta
q=1$ and varying $\chi$.

To gain an understanding of the region of bounded motion for test particles in
the HT metric, one can begin by comprehending the effective potential
$V_{\text{eff}}$ and the corresponding CZVs for various combinations of
parameters. This study should provide information about the regions where
bound motion is permitted, as well as the motion’s form. For instance, it
could clarify whether the CZV is entirely closed and geodesics starting inside
it can never escape, or if it is open toward the surface of the compact object
and can lead to plunging orbits. Since the HT metric is symmetric in the
equatorial plane, it is enough to plot the potential for $\theta=\pi/2$, where
it achieves its deepest valley (if it exists).

Fig. 1 demonstrates the behavior of both the equatorial potential and the CZVs
for varying $E$ and $L_{z}$. We observe that as the energy of the particle
increases, the potential forms a bigger valley and this leads to a larger
volume in the CZV where bound motion can occur. Notice that for the particular
choice of fixed $L_{z}$ in Fig. 1, all CZV are closed so there are no plunging
orbits that begin from inside the region of bound motion and abide to the
appropriate initial conditions and constraint equation. On the contrary, we
may fix $E$ and vary $L_{z}$. Then by increasing $L_{z}$ we observe that the
potential changes from one that intersects $V_{\text{eff}}=0$ only once to one
that forms a valley. Therefore, the increment of $L_{z}$ changes the nature of
the corresponding CZVs from open, that allow plunging orbits, to closed where
only bound motion occurs. It is curious that, in the case of the largest
$L_{z}$ presented in Figure 1, there is a secondary region very close to the
compact object where only plunging orbits occur. Therefore, we will avoid
these regimes in phase space when they appear.

##  V Results

To understand the appearance of indirect chaos on a Poincaré map, we present a
typical map of the HT metric with a quadrupole deformation of $\delta q=1$ in
Figure 2. Although the majority of the curves formed surround the central
periodic point in an ordered fashion, as if the system were integrable, a
thorough examination reveals the presence of a $2/3$ resonant island
(highlighted in red). Upon closer inspection of the leftmost island of the
Birkhoff chain, we observe the encapsulating structure of the island and how
the resonant KAM curves surround the stable periodic point of the resonance
rather than that of the map. This is an indirect manifestation of chaos, and
the system being studied is nonintegrable, at least up to second order in the
rotation rate, denoted as $\mathcal{O}(\Omega^{2})$. In a recent study [59],
the authors analyzed photon orbits for the HT metric at the next leading
order, i.e., $\mathcal{O}(\Omega^{3})$, and still observed chaotic effects.
Therefore, we do not expect that including even higher orders of rotation will
make the system integrable, as has been observed in Kerr spacetime [85].
Figure 3 completes our basic demonstration of nonintegrability. For the same
parameters as in Fig. 2, we sketch the rotation curve from the left edge of
the CZV up to the central point of the Poincaré map. We find that the
existence of the $2/3$ island breaks the monotonicity of the rotation curve
and at $\nu_{\vartheta}=2/3$ it presents a plateau where the rotation number
remains constant for a range of initial conditions $r(0)/M$. To determine the
width of the island, we run geodesics near the entrance and exit of the
plateau and extract the initial conditions that result in a $2/3$ island in
the Poincaré map. This ensures that we are inside the island upon entry and
just before exit.

Figure 5: Top left: $2/3$ resonant island widths $w/M$ of geodesics in HT
metric with fixed $L_{z}/\mu=3.5M$, $\delta q=1$, $\chi=0.3$ and varying
$E/\mu$. Top right: Same us top left with fixed $E/\mu=0.97$, $\delta q=1$,
$\chi=0.3$ and varying $L_{z}/\mu$. Bottom left: Same as top left with fixed
$E/\mu=0.97$, $L_{z}/\mu=3.5M$, $\chi=0.3$ and varying $\delta q$. Bottom
right: Same as top left with fixed $E/\mu=0.97$, $L_{z}/\mu=3.5M$, $\delta
q=1$ and varying $\chi$.

###  V.1 Varying $E/\mu$

As demonstrated in Fig. 1, the relationship between the potential and CZV is
proportional to the magnitude of $E$. The rotation curves displayed in Fig. 4
(top left) exhibit the anticipated behavior, whereby an increase in the
particle’s energy results in rotation curves spanning a greater range of
$r(0)/M$ due to the expansion of the CZV volume. If $E/\mu=0.95$, the $2/3$
resonance is not satisfied by any bound geodesic. This can be misleading and
result in incorrect assumptions about integrability since there are no
plateaus. However, changing $E$ will eventually disrupt the smoothness of the
rotation curve. The insets in Fig. 4 (top left) zoom in on the regions where
the plateau occurs. We observe that for the cases presented, the width of the
plateau’s width increases with $E$. In fact, Fig. 5 (top left) shows the
dependence of the width of the $2/3$ island on $E$, which clearly demonstrates
that the width increases when the particle is initialized with more energy.

###  V.2 Varying $L_{z}/\mu$

According to Fig. 1, if we increase $L_{z}$, then the width of the island is
expected to decrease as the volume of the CZV is reduced, even though it
becomes closed. This means that the geodesic has access to fewer initial
conditions in general. In this case, the dependence of the potential and CZV
is inversely proportional to the size of $L_{z}$. The rotation curves in Fig.
4 (top right) demonstrate this behavior, as increasing the momentum of the
particles leads to rotation curves that cover smaller ranges of $r(0)/M$. The
insets shown in the top right of Fig. 4 zoom onto the regions where the
plateau occurs. For the cases presented, the plateau widths decrease as
$L_{z}$ increases. Figure 5 (top right) confirms this trend, as the width of
the $2/3$ island decreases with increasing $L_{z}$. Thus, as the particle’s
$z$-component of angular momentum is reduced, the island’s width increases.
However, from Fig. 5, it is also apparent that by decreasing $L_{z}$ even
further while keeping the other parameters constant, we can reach islands with
widths of around $\mathcal{O}(10^{-1}M)$, which are quite large for slowly-
rotating compact objects. This is especially true when compared to rapidly-
rotating non-Kerr and bumpy spacetimes [67, 68, 69]. It should be noted that
rotating boson stars, depending on their compactness, can reach unprecedented
levels of island widths of order $\mathcal{O}(M)$ [70].

###  V.3 Varying $\chi$

The expectation that the increase in spin will expand the width of the island
has been thoroughly demonstrated in [64]. Our results align with this analysis
even though we are limited in how large the HT spin can be due to the
perturbative nature of the metric. Fig. 4 (bottom right) shows exactly what we
expected: the increase in $\chi$ results in larger $2/3$ islands but smaller
CZVs as the range of $r(0)/M$ is reduced. Figure 5 is consistent with this and
other analyses found in the literature [64, 86, 58].

###  V.4 Varying $\delta q$

The final and most crucial parameter that we have varied is the quadrupole
deformation $\delta q$. We kept $\delta q=1$ as a reference and increased it
accordingly. Figures 4 and 5 (bottom left) show that the further we deviate
from the Kerr mass quadrupole, the more noticeable the chaotic phenomenon
becomes. This finding aligns with non-Kerr spacetimes as reported in [64].
Island widths increase significantly as the deformation increases, reaching
widths on the order of $\mathcal{O}(10^{-1}M)$ for a sufficiently large
$\delta q$. This is a fascinating result, particularly considering that the
spacetime is slowly rotating. However, it should be noted that the HT metric
is only an approximate solution to the field equations, and the quadrupole
deformation must be very large for prominent chaotic phenomena to occur. Based
on the analysis presented above, it is clear that particle geodesics around
the HT metric are nonintegrable to at least $\mathcal{O}(\Omega^{2})$ spin,
and possibly even higher orders in $\Omega$ (as discussed in [59]).

##  VI Conclusions and Discussion

In this work, we have studied the geodesics of massive particles around the HT
metric. We have explored the available parameter space of approximate HT
solutions up to $\mathcal{O}(\Omega^{2})$ in rotation order by evolving non-
spherical, generic particle orbits within the allowed regions of bound motion.
Our analysis not only reveals how the structure of the metric and the region
of bound geodesics changes when we vary the parameters of spacetime and the
particle but also the fact that geodesics are nonintegrable. This observation
was also found for photon orbits in HT spacetime [59]. This means that the
evolution of particles is not deterministic, and chaotic effects occur.

We demonstrated how chaotic features arise using two tools: the Poincaré map
and the corresponding rotation curve of successive geodesics with varying
initial conditions. Poincaré maps reveal that most of the phase space seems to
generate integrable geodesics around the compact object, with KAM curves
enclosing the central periodic point of the map, until a resonance is reached.
In these regions, KAM curves break down into periodic stable orbits surrounded
by resonant KAM curves and encapsulated by islands of stability. The indirect
indication of chaos is that the rotation number is shared with that of the
central periodic point (in our case, $\nu_{\vartheta}=2/3$) for all geodesics
that occupy the island. This phenomenon does not appear in integrable
spacetimes like Kerr. The shared rotation number translates to a plateau in
the rotation curve. We expect other subdominant resonances to exhibit the same
features, but only when the deformation parameter is much larger.
Alternatively, the deformation can be chosen to be oblate, meaning $\delta
q<0$. However, in [59] it is argued that oblate deformations do not generate
chaotic effects as strong as those for prolate deformations. We did not
explore extreme mass quadrupole deformations to avoid deviating too much from
General Relativity and to satisfy the pathological bounds of the HT metric
[75].

Our extensive analysis reveals the following:

(i) As we increase the particle’s energy ($E$), the volume of CZVs increases,
allowing more geodesics to be bound around the HT spacetime. Consequently, the
widths of $2/3$ resonant islands increase with $E$.

(ii) Increasing the particle’s $z$-component of angular momentum has the
opposite effect to that in (i). CZVs change from an open configuration, which
allows for plunging orbits, to a closed and smaller-volume one. Thus, the
island widths decrease with $L_{z}$.

(iii) The widths of $2/3$ resonant islands are proportional to both spacetime
parameters, namely quadrupole deformation ($\delta q$) and spin ($\chi$).
Specifically, for a suitable combination of large $E$, $\delta q$, and spin,
as well as small $L_{z}$, island widths of order $\mathcal{O}(w/M)\gtrsim
10^{-1}$, as in [64], can be achieved. This is of interest since the HT metric
is slowly-rotating.

Our results clarify that it is necessary to perform a similar analysis for
orders higher than $\Omega^{2}$. This is because similar results, as presented
in this study, can be expected, as observed in photon orbits [59]. Moreover,
to construct a solution for rapidly rotating NSs and investigate their
integrability, it is essential to develop numerical methods, as outlined in
[32, 33, 21], and to evolve geodesics around the numerically generated
solution, as discussed in [70], which is a much more demanding task.
Additionally, it may be worth exploring whether resonant islands and chaotic
effects impact quasiperiodic oscillations. The ultimate goal would be to
estimate the mass, angular momentum, and multipole moments of accreting
compact objects associated with low-mass X-ray binaries to distinguish between
BHs and NSs with multimessenger astrophysics [87, 88, 89].

## Acknowledgments

The authors warmly thank Prof. G. Pappas, Prof. Th. Apostolatos, Prof. P. Pani
and Dr. P. Manoharan for critical comments and helpful discussions. This work
was supported by the DAAD program for the “Promotion of the exchange and
scientific cooperation between Greece and Germany IKYDAAD 2022” (57628320).
K.D. and K.K. are grateful for hospitality provided by the Department of
Physics of the University of Athens, Greece. K.D. acknowledges financial
support provided under the European Union’s H2020 ERC, Starting Grant
agreement no. DarkGRA–757480 and the MIUR PRIN and FARE programmes (GW-NEXT,
CUP: B84I20000100001).

## References

  * Taylor  _et al._ [1979] J. H. Taylor, L. A. Fowler, and P. M. McCulloch, Measurements of general relativistic effects in the binary pulsar PSR 1913+16, Nature 277, 437 (1979).
  * van der Klis [2000] M. van der Klis, Millisecond oscillations in x-ray binaries, Ann. Rev. Astron. Astrophys. 38, 717 (2000), arXiv:astro-ph/0001167 .
  * Kluzniak and Abramowicz [2003] W. Kluzniak and M. Abramowicz, General Relativity and Gravitation. Millisecond oscillators in accreting neutron stars and black holes, in  _12th Workshop on General Relativity and Gravitation_ (2003) arXiv:astro-ph/0304345 .
  * Lewin and van der Klis [2006] W. H. G. Lewin and M. van der Klis,  _Compact Stellar X-ray Sources_ , Vol. 39 (2006).
  * Poynting [1903] J. H. Poynting, “Radiation in the Solar System; its Effect on Temperature and its Pressure on Small Bodies.”, Monthly Notices of the Royal Astronomical Society 64, 1a (1903), https://academic.oup.com/mnras/article-pdf/64/1/1a/3438550/mnras64-0001a.pdf .
  * Robertson and Russell [1937] H. P. Robertson and H. N. Russell, Dynamical Effects of Radiation in the Solar System, Monthly Notices of the Royal Astronomical Society 97, 423 (1937), https://academic.oup.com/mnras/article-pdf/97/6/423/3169085/mnras97-0423.pdf .
  * Akiyama  _et al._ [2019] K. Akiyama _et al._ (Event Horizon Telescope), First M87 Event Horizon Telescope Results. I. The Shadow of the Supermassive Black Hole, Astrophys. J. Lett. 875, L1 (2019), arXiv:1906.11238 [astro-ph.GA] .
  * Akiyama  _et al._ [2022] K. Akiyama _et al._ (Event Horizon Telescope), First Sagittarius A* Event Horizon Telescope Results. I. The Shadow of the Supermassive Black Hole in the Center of the Milky Way, Astrophys. J. Lett. 930, L12 (2022).
  * Abbott  _et al._ [2016] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Observation of Gravitational Waves from a Binary Black Hole Merger, Phys. Rev. Lett. 116, 061102 (2016), arXiv:1602.03837 [gr-qc] .
  * Abbott  _et al._ [2021a] R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run, (2021a), arXiv:2111.03606 [gr-qc] .
  * Abbott  _et al._ [2017a] B. P. Abbott _et al._ (LIGO Scientific, Virgo), GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral, Phys. Rev. Lett. 119, 161101 (2017a), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ [2021b] R. Abbott _et al._ (LIGO Scientific, KAGRA, VIRGO), Observation of Gravitational Waves from Two Neutron Star–Black Hole Coalescences, Astrophys. J. Lett. 915, L5 (2021b), arXiv:2106.15163 [astro-ph.HE] .
  * Barack  _et al._ [2019] L. Barack _et al._ , Black holes, gravitational waves and fundamental physics: a roadmap, Class. Quant. Grav. 36, 143001 (2019), arXiv:1806.05195 [gr-qc] .
  * Barausse  _et al._ [2014] E. Barausse, V. Cardoso, and P. Pani, Can environmental effects spoil precision gravitational-wave astrophysics?, Phys. Rev. D 89, 104059 (2014), arXiv:1404.7149 [gr-qc] .
  * Cardoso  _et al._ [2022a] V. Cardoso, K. Destounis, F. Duque, R. P. Macedo, and A. Maselli, Black holes in galaxies: Environmental impact on gravitational-wave generation and propagation, 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, Gravitational Waves from Extreme-Mass-Ratio Systems in Astrophysical Environments, Phys. Rev. Lett. 129, 241103 (2022b), arXiv:2210.01133 [gr-qc] .
  * Destounis  _et al._ [2023a] K. Destounis, A. Kulathingal, K. D. Kokkotas, and G. O. Papadopoulos, Gravitational-wave imprints of compact and galactic-scale environments in extreme-mass-ratio binaries, Phys. Rev. D 107, 084027 (2023a), arXiv:2210.09357 [gr-qc] .
  * Cheung  _et al._ [2022] M. H.-Y. Cheung, K. Destounis, R. P. Macedo, E. Berti, and V. Cardoso, Destabilizing the Fundamental Mode of Black Holes: The Elephant and the Flea, Phys. Rev. Lett. 128, 111103 (2022), arXiv:2111.05415 [gr-qc] .
  * Abbott  _et al._ [2017b] B. P. Abbott _et al._ (LIGO Scientific, Virgo, Fermi-GBM, INTEGRAL), Gravitational Waves and Gamma-rays from a Binary Neutron Star Merger: GW170817 and GRB 170817A, Astrophys. J. Lett. 848, L13 (2017b), arXiv:1710.05834 [astro-ph.HE] .
  * Kerr [1963] R. P. Kerr, Gravitational field of a spinning mass as an example of algebraically special metrics, Phys. Rev. Lett. 11, 237 (1963).
  * Berti and Stergioulas [2004] E. Berti and N. Stergioulas, Approximate matching of analytic and numerical solutions for rapidly rotating neutron stars, Mon. Not. Roy. Astron. Soc. 350, 1416 (2004), arXiv:gr-qc/0310061 .
  * Pappas [2012] G. Pappas, What can quasi-periodic oscillations tell us about the structure of the corresponding compact objects?, Mon. Not. Roy. Astron. Soc. 422, 2581 (2012), arXiv:1201.6071 [astro-ph.HE] .
  * Manko and Sibgatullin [1993] V. S. Manko and N. R. Sibgatullin, Construction of exact solutions of the Einstein-Maxwell equations corresponding to a given behaviour of the Ernst potentials on the symmetry axis, Classical and Quantum Gravity 10, 1383 (1993).
  * Manko  _et al._ [1995] V. S. Manko, J. Martín, and E. Ruiz, Six-parameter solution of the Einstein–Maxwell equations possessing equatorial symmetry, Journal of Mathematical Physics 36, 3063 (1995).
  * Manko  _et al._ [2000a] V. S. Manko, J. D. Sanabria-Gomez, and O. V. Manko, Nine parameter electrovac metric involving rational functions, Phys. Rev. D 62, 044048 (2000a).
  * Manko  _et al._ [2000b] V. S. Manko, E. W. Mielke, and J. D. Sanabria-Gomez, Exact solution for the exterior field of a rotating neutron star, Phys. Rev. D 61, 081501 (2000b), arXiv:gr-qc/0001081 .
  * Stute and Camenzind [2002] M. Stute and M. Camenzind, Towards a self-consistent relativistic model of the exterior gravitational field of rapidly rotating neutron stars, Mon. Not. Roy. Astron. Soc. 336, 831 (2002), arXiv:astro-ph/0301466 .
  * Pachon  _et al._ [2006] L. A. Pachon, J. A. Rueda, and J. D. Sanabria-Gomez, Realistic exact solution for the exterior field of a rotating neutron star, Phys. Rev. D 73, 104038 (2006), arXiv:gr-qc/0606060 .
  * Teichmuller  _et al._ [2011] C. Teichmuller, M. B. Fröb, and F. Maucher, Analytical approximation of the exterior gravitational field of rotating neutron stars, Class. Quant. Grav. 28, 155015 (2011), arXiv:1102.5252 [gr-qc] .
  * Pappas and Apostolatos [2013] G. Pappas and T. A. Apostolatos, An all-purpose metric for the exterior of any kind of rotating neutron star, Mon. Not. Roy. Astron. Soc. 429, 3007 (2013), arXiv:1209.6148 [gr-qc] .
  * Pappas [2017] G. Pappas, An accurate metric for the spacetime around rotating neutron stars, Mon. Not. Roy. Astron. Soc. 466, 4381 (2017), arXiv:1610.05370 [gr-qc] .
  * Stergioulas and Friedman [1995] N. Stergioulas and J. L. Friedman, Comparing models of rapidly rotating relativistic stars constructed by two numerical methods, Astrophys. J. 444, 306 (1995), arXiv:astro-ph/9411032 .
  * Stergioulas [2003] N. Stergioulas, Rotating Stars in Relativity, Living Rev. Rel. 6, 3 (2003), arXiv:gr-qc/0302034 .
  * Hartle and Thorne [1968] J. B. Hartle and K. S. Thorne, Slowly Rotating Relativistic Stars. II. Models for Neutron Stars and Supermassive Stars, ApJ 153, 807 (1968).
  * Hartle [1967] J. B. Hartle, Slowly Rotating Relativistic Stars. I. Equations of Structure, ApJ 150, 1005 (1967).
  * Hartle and Thorne [1969] J. B. Hartle and K. S. Thorne, Slowly Rotating Relativistic Stars. III. Static Criterion for Stability, ApJ 158, 719 (1969).
  * Hartle [1970] J. B. Hartle, Slowly-Rotating Relativistic Stars.IV. Rotational Energy and Moment of Inertia for Stars in Differential Rotation, ApJ 161, 111 (1970).
  * Hartle  _et al._ [1972] J. B. Hartle, K. S. Thorne, and S. M. Chitre, Slowly Rotating Relativistic Stars.VI. Stability of the Quasiradial Modes, ApJ 176, 177 (1972).
  * Hartle [1975] J. B. Hartle, Slowly rotating relativistic stars. IIIA. The static stability criterion recovered., ApJ 195, 203 (1975).
  * Lense and Thirring [1918] J. Lense and H. Thirring, On the Influence of the Proper Rotation of Central Bodies on the Motions of Planets and Moons According to Einstein’s Theory of Gravitation, Physikalische Zeitschrift 19, 156 (1918).
  * Vieira  _et al._ [2022] H. S. Vieira, K. Destounis, and K. D. Kokkotas, Slowly-rotating curved acoustic black holes: Quasinormal modes, Hawking-Unruh radiation, and quasibound states, Phys. Rev. D 105, 045015 (2022), arXiv:2112.08711 [gr-qc] .
  * Abramowicz  _et al._ [2003] M. A. Abramowicz, G. J. E. Almergren, W. Kluzniak, and A. V. Thampan, Circular geodesics in the Hartle-Thorne metric, (2003), arXiv:gr-qc/0312070 .
  * van der Klis  _et al._ [1985] M. van der Klis, F. Jansen, J. van Paradijs, W. H. G. Lewin, E. P. J. van den Heuvel, J. E. Trumper, and M. Szatjno, Intensity-dependent quasi-periodic oscillations in the X-ray flux of GX5-1, Nature 316, 225 (1985).
  * Abramowicz  _et al._ [2002] M. A. Abramowicz, G. J. E. Almergren, W. Kluzniak, A. V. Thampan, and F. Wallinder, Holonomy invariance, orbital resonances, and kilohertz QPO(s), Class. Quant. Grav. 19, L57 (2002), arXiv:gr-qc/0202020 .
  * Ingram  _et al._ [2016] A. Ingram, M. van der Klis, M. Middleton, C. Done, D. Altamirano, L. Heil, P. Uttley, and M. Axelsson, A quasi-periodic modulation of the iron line centroid energy in the black hole binary H1743$-$322, Mon. Not. Roy. Astron. Soc. 461, 1967 (2016), arXiv:1607.02866 [astro-ph.HE] .
  * Urbancová  _et al._ [2019] G. Urbancová, M. Urbanec, G. Török, Z. Stuchlík, M. Blaschke, and J. C. Miller, Epicyclic Oscillations in the Hartle–Thorne External Geometry, Astrophys. J. 877, 66 (2019), arXiv:1905.00730 [astro-ph.HE] .
  * Sulieva  _et al._ [2022] G. Sulieva, K. Boshkayev, G. Nurbakyt, H. Quevedo, A. Taukenova, A. Tlemissov, Z. Tlemissova, and A. Urazalina, Adiabatic theory of motion of bodies in the Hartle-Thorne spacetime, (2022), arXiv:2205.04217 [gr-qc] .
  * Carter [1968] B. Carter, Global structure of the Kerr family of gravitational fields, Phys. Rev. 174, 1559 (1968).
  * Papadopoulos and Kokkotas [2018] G. O. Papadopoulos and K. D. Kokkotas, Preserving Kerr symmetries in deformed spacetimes, Class. Quant. Grav. 35, 185014 (2018), arXiv:1807.08594 [gr-qc] .
  * Papadopoulos and Kokkotas [2021] G. O. Papadopoulos and K. D. Kokkotas, On Kerr black hole deformations admitting a Carter constant and an invariant criterion for the separability of the wave equation, Gen. Rel. Grav. 53, 21 (2021), arXiv:2007.12125 [gr-qc] .
  * Cornish [2001] N. J. Cornish, Chaos and gravitational waves, Phys. Rev. D 64, 084011 (2001), arXiv:gr-qc/0106062 .
  * Cornish and Levin [2003] N. J. Cornish and J. J. Levin, Lyapunov timescales and black hole binaries, Class. Quant. Grav. 20, 1649 (2003), arXiv:gr-qc/0304056 .
  * Verhaaren and Hirschmann [2010] C. Verhaaren and E. W. Hirschmann, Chaotic orbits for spinning particles in Schwarzschild spacetime, Phys. Rev. D 81, 124034 (2010), arXiv:0912.0031 [gr-qc] .
  * Barausse  _et al._ [2007] E. Barausse, L. Rezzolla, D. Petroff, and M. Ansorg, Gravitational waves from Extreme Mass Ratio Inspirals in non-pure Kerr spacetimes, Phys. Rev. D 75, 064026 (2007), arXiv:gr-qc/0612123 .
  * Barausse and Rezzolla [2008] E. Barausse and L. Rezzolla, The Influence of the hydrodynamic drag from an accretion torus on extreme mass-ratio inspirals, Phys. Rev. D 77, 104027 (2008), arXiv:0711.4558 [gr-qc] .
  * Suzuki and Maeda [1997] S. Suzuki and K.-i. Maeda, Chaos in Schwarzschild space-time: The motion of a spinning particle, Phys. Rev. D 55, 4848 (1997), arXiv:gr-qc/9604020 .
  * Zelenka  _et al._ [2020] O. Zelenka, G. Lukes-Gerakopoulos, V. Witzany, and O. Kopáček, Growth of resonances and chaos for a spinning test particle in the Schwarzschild background, Phys. Rev. D 101, 024037 (2020), arXiv:1911.00414 [gr-qc] .
  * Lukes-Gerakopoulos  _et al._ [2014] G. Lukes-Gerakopoulos, G. Contopoulos, and T. A. Apostolatos, Non-Linear Effects in Non-Kerr spacetimes, Springer Proc. Phys. 157, 129 (2014), arXiv:1408.4697 [gr-qc] .
  * Kostaros and Pappas [2022] K. Kostaros and G. Pappas, Chaotic photon orbits and shadows of a non-Kerr object described by the Hartle–Thorne spacetime, Class. Quant. Grav. 39, 134001 (2022), arXiv:2111.09367 [gr-qc] .
  * Contopoulos [2003] G. Contopoulos,  _Order and Chaos in Dynamical Astronomy_ (Springer-Verlag, New York, 2003).
  * Glampedakis [2005] K. Glampedakis, Extreme mass ratio inspirals: LISA’s unique probe of black hole gravity, Class. Quant. Grav. 22, S605 (2005), arXiv:gr-qc/0509024 .
  * Amaro-Seoane [2018] P. Amaro-Seoane, Relativistic dynamics and extreme mass ratio inspirals, Living Rev. Rel. 21, 4 (2018), arXiv:1205.5240 [astro-ph.CO] .
  * Apostolatos  _et al._ [2009] T. A. Apostolatos, G. Lukes-Gerakopoulos, and G. Contopoulos, How to Observe a Non-Kerr Spacetime Using Gravitational Waves, 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, Observable signature of a background deviating from the Kerr metric, Phys. Rev. D 81, 124005 (2010), arXiv:1003.3120 [gr-qc] .
  * Lukes-Gerakopoulos [2012] G. Lukes-Gerakopoulos, The non-integrability of the Zipoy-Voorhees metric, Phys. Rev. D 86, 044013 (2012), arXiv:1206.0660 [gr-qc] .
  * Destounis  _et al._ [2023b] K. Destounis, G. Huez, and K. D. Kokkotas, Geodesics and gravitational waves in chaotic extreme-mass-ratio inspirals: The curious case of Zipoy-Voorhees black-hole mimickers, (2023b), arXiv:2301.11483 [gr-qc] .
  * Destounis  _et al._ [2020] K. Destounis, A. G. Suvorov, and K. D. Kokkotas, Testing spacetime symmetry through gravitational waves from extreme-mass-ratio inspirals, Phys. Rev. D 102, 064041 (2020), arXiv:2009.00028 [gr-qc] .
  * Destounis  _et al._ [2021] K. Destounis, A. G. Suvorov, and K. D. Kokkotas, Gravitational-wave glitches in chaotic extreme-mass-ratio inspirals, Phys. Rev. Lett. 126, 141102 (2021), arXiv:2103.05643 [gr-qc] .
  * Destounis and Kokkotas [2021] K. Destounis and K. D. Kokkotas, Gravitational-wave glitches: Resonant islands and frequency jumps in nonintegrable extreme-mass-ratio inspirals, Phys. Rev. D 104, 064023 (2021), arXiv:2108.02782 [gr-qc] .
  * Destounis  _et al._ [2023c] K. Destounis, F. Angeloni, M. Vaglio, and P. Pani, Extreme-mass-ratio inspirals into rotating boson stars: nonintegrability, chaos, and transient resonances, (2023c), arXiv:2305.05691 [gr-qc] .
  * Pacilio  _et al._ [2020] C. Pacilio, M. Vaglio, A. Maselli, and P. Pani, Gravitational-wave detectors as particle-physics laboratories: Constraining scalar interactions with a coherent inspiral model of boson-star binaries, Phys. Rev. D 102, 083002 (2020), arXiv:2007.05264 [gr-qc] .
  * Vaglio  _et al._ [2022] M. Vaglio, C. Pacilio, A. Maselli, and P. Pani, Multipolar structure of rotating boson stars, Phys. Rev. D 105, 124020 (2022), arXiv:2203.07442 [gr-qc] .
  * Vaglio  _et al._ [2023] M. Vaglio, C. Pacilio, A. Maselli, and P. Pani, Bayesian parameter estimation on boson-star binary signals with a coherent inspiral template and spin-dependent quadrupolar corrections, (2023), arXiv:2302.13954 [gr-qc] .
  * De Falco  _et al._ [2020] V. De Falco, P. Bakala, and M. Falanga, Three-dimensional general relativistic Poynting-Robertson effect. III. Static and nonspherical quadrupolar massive source, Phys. Rev. D 101, 124031 (2020), arXiv:2006.01452 [gr-qc] .
  * De Falco and Wielgus [2021] V. De Falco and M. Wielgus, Three-dimensional general relativistic Poynting-Robertson effect. IV. Slowly rotating and nonspherical quadrupolar massive source, Phys. Rev. D 103, 084056 (2021), arXiv:2103.17165 [gr-qc] .
  * Benhar  _et al._ [2005] O. Benhar, V. Ferrari, L. Gualtieri, and S. Marassi, Perturbative approach to the structure of rapidly rotating neutron stars, Phys. Rev. D 72, 044028 (2005), arXiv:gr-qc/0504068 .
  * Glampedakis and Pappas [2019] K. Glampedakis and G. Pappas, Modification of photon trapping orbits as a diagnostic of non-Kerr spacetimes, Phys. Rev. D 99, 124041 (2019), arXiv:1806.09333 [gr-qc] .
  * Cardoso and Pani [2019] V. Cardoso and P. Pani, Testing the nature of dark compact objects: a status report, Living Rev. Rel. 22, 4 (2019), arXiv:1904.05363 [gr-qc] .
  * Maggio  _et al._ [2020] E. Maggio, P. Pani, and G. Raposo, Testing the nature of dark compact objects with gravitational waves, in  _Handbook of Gravitational Wave Astronomy_, edited by C. Bambi, S. Katsanevas, and K. D. Kokkotas (Springer Singapore, Singapore, 2020) pp. 1–37.
  * Bini  _et al._ [2013] D. Bini, K. Boshkayev, R. Ruffini, and I. Siutsou, Equatorial Circular Geodesics in the Hartle-Thorne Spacetime, Nuovo Cim. C 036, 31 (2013), arXiv:1306.4792 [gr-qc] .
  * Contopoulos  _et al._ [2011] G. Contopoulos, G. Lukes-Gerakopoulos, and T. A. Apostolatos, Orbits in a non-Kerr Dynamical System, Int. J. Bifurc. Chaos 21, 2261 (2011), arXiv:1108.5057 [gr-qc] .
  * Arnol’d [1963] V. I. Arnol’d, Proof of a theorem of a. n. kolmogorov on the invariance of quasi-periodic motions under small perturbations of the hamiltonian, Russian Mathematical Surveys 18, 9 (1963).
  * Möser [1962] J. Möser, On invariant curves of area-preserving mappings of an annulus, Nachr. Akad. Wiss. Göttingen, II , 1 (1962).
  * Birkhoff [1913] G. D. Birkhoff, Proof of poincaré’s geometric theorem, Transactions of the American Mathematical Society 14, 14 (1913).
  * Deich  _et al._ [2022] A. Deich, A. Cárdenas-Avendaño, and N. Yunes, Chaos in quadratic gravity, Phys. Rev. D 106, 024040 (2022), arXiv:2203.00524 [gr-qc] .
  * Lukes-Gerakopoulos and Contopoulos [2013] G. Lukes-Gerakopoulos and G. Contopoulos, Mind the Resonances: Final stages of accretion into bumpy black holes, J. Phys. Conf. Ser. 453, 012005 (2013), arXiv:1304.7612 [gr-qc] .
  * Boshkayev  _et al._ [2014] K. Boshkayev, D. Bini, J. Rueda, A. Geralico, M. Muccino, and I. Siutsou, What can we extract from quasi-periodic oscillations?, Grav. Cosmol. 20, 233 (2014), arXiv:1412.8214 [astro-ph.HE] .
  * Boshkayev  _et al._ [2015] K. Boshkayev, J. Rueda, and M. Muccino, Extracting multipole moments of neutron stars from quasi-periodic oscillations in low mass X-ray binaries, Astronomy Reports 59, 441 (2015).
  * Boshkayev  _et al._ [2017] K. Boshkayev, J. A. Rueda, and M. Muccino, Main parameters of neutron stars from quasi-periodic oscillations in low mass X-ray binaries, in  _14th Marcel Grossmann Meeting on Recent Developments in Theoretical and Experimental General Relativity, Astrophysics, and Relativistic Field Theories_, Vol. 4 (2017) pp. 3433–3440, arXiv:1604.02398 [astro-ph.HE] .