# Detection of chaos in the general relativistic
Poynting-Robertson effect: Kerr equatorial plane

Vittorio De Falco1<EMAIL_ADDRESS><EMAIL_ADDRESS>William Borrelli2<EMAIL_ADDRESS>1Department of Mathematics and
Applications “R. Caccioppoli”, University of Naples Federico II, Via Cintia,
80126 Naples, Italy,
2Dipartimento di Matematica e Fisica, Università Cattolica del Sacro Cuore,
Via dei Musei 41, Brescia, Italy.

###### Abstract

The general relativistic Poynting-Robertson effect is a dissipative and non-
linear dynamical system obtained by perturbing through radiation processes the
geodesic motion of test particles orbiting around a spinning compact object,
described by the Kerr metric. Using the Melnikov method we find that, in a
suitable range of parameters, chaotic behavior is present in the motion of a
test particle driven by the Poynting-Robertson effect in the Kerr equatorial
plane.

##  I Introduction

Chaos is a widespread feature in many physical non-linear dynamical systems. A
chaotic system hides behind the visible randomness of the complex dynamics,
some underlying rich mathematical structures, such as: constant feedback
loops, self-similarities, fractals, and self-organization Wiggins (1988); Ott
(2002). Although an universally accepted formal definition of chaos does not
exist, the one due to Robert L. Devaney is widely accepted, and it is based on
the following three proprieties Devaney (2018):

  * •

$(1)$ _sensitive dependence on initial conditions_ , i.e., tiny perturbations
on the initial conditions leads to significantly different future behaviors,

  * •

$(2)$ _topologically mixing_ , i.e., any given region or open set of the phase
space eventually overlaps with any other given region in the phase space;

  * •

$(3)$ _presence of a dense set of periodic orbits_ , i.e., every point in the
dynamical real space is approached arbitrarily close by periodic orbits.

General Relativity (GR), being a non-linear theory, can potentially exhibit
chaotic phenomena Hobill  _et al._ (1994). The studies on chaos in GR can be
mainly divided in two branches: (1) problems of geodesic/non-geodesic motion
of a particle in a given gravitational field; (2) evolution of cosmological
models. Regarding works on the first class, it is important to mention: the
motion in spaces with negative curvature Arnol’d and Avez (1989), motion
around two fixed black holes (BHs) Contopoulos (1990, 1991); Hobill  _et al._
(1994), relativistic restricted three-body problem Wanex (2002), Schwarzschild
BH affected by high-frequency periodic perturbations Bombelli and Calzetta
(1992), spinning particle motion around a Kerr and Schwarzschild BH Suzuki and
Maeda (1997); Lukes-Gerakopoulos (2018), gravitational waves from spinning
compact binaries Cornish (2001); Cornish and Levin (2002, 2003). Moreover,
studies on chaos in cosmology include: the model of Belinski-Khalatnikov-
Lifshitz dealing with the dynamic evolution of the universe near the initial
gravitational singularity, described by an anisotropic, chaotic solution of
the Einstein field equation of gravitation. Belinskij  _et al._ (1970a, b),
Bianchi IX (“mixmaster universe”) Barrow and Sirousse-Zia (1989); Burd  _et
al._ (1991); Contopoulos  _et al._ (1999), Friedmann-Robertson-Walker (FRW)
plus a massive scalar field Calzetta and El Hasi (1993), and the non-linear
interaction among dark matter, dark energy, normal matter, and radiation on
the FRW spacetime Aydiner (2016).

In high-energy astrophysics, dealing with electromagnetic radiation processes
around compact objects, like neutron stars (NSs) or BHs, relatively small-
sized test particles can drastically depart from their geodesic motion. The
gravitational pull is contrasted by the radiation pressure, and in the process
of absorption and reemission of radiation from the test particle an additional
radiation torque appears, acting as a drag force opposite to the test particle
orbital motion Poynting (1903); Robertson (1937). _This is the so-called
Poynting-Robertson (PR) effect, which configures thus as a dissipative non-
linear dynamical system efficiently removing energy and angular momentum from
the affected test particle_. There are several models of the general
relativistic PR effect in Kerr and also other spacetimes from the two
dimensional (2D) Bini  _et al._ (2009, 2011, 2015) to the three dimensional
(3D) formulations De Falco  _et al._ (2019); Bakala  _et al._ (2019); Wielgus
(2019); De Falco  _et al._ (2020); De Falco (2019). They all exhibit the
existence of a critical hypersurface, a region where gravitational and
radiation forces balance and the test particle moves on it stably Bini  _et
al._ (2011); De Falco and Bakala (2020); De Falco (2020).

Here, we focus our attention on the _general relativistic PR effect in the
equatorial plane around a Kerr compact object_. To search for chaotic
behavior, we employ the _Melnikov method_ Guckenheimer and Holmes (2002);
Wiggins (1988), which is an independent diagnostic procedure, complementary to
other numerical and analytical methods Tabor (1989); Ott (2002); Guckenheimer
and Holmes (2002). Its strength relies on the fact, that _it requires only the
knowledge of few elements without having any insight into the solution of the
perturbed dynamics, i.e.: (1) invariant subsets in the phase space of the
unperturbed dynamics (homoclinic orbits); (2) explicit expression of the
perturbations_.

The paper is organized as follows: in Sec. II we briefly recall the general
relativistic PR effect model, underlining also how to derive its dissipative
perturbations; in Sec. III the homoclinic orbits in the equatorial plane of
the Kerr spacetime are described; in Sec. IV we apply the Melnikov method to
the general relativistic PR effect; in Sec. V we discuss our results and draw
the conclusions.

##  II General relativistic Poynting-Robertson effect in the equatorial plane
of the Kerr metric

The general relativistic PR effect in the Kerr metric describes the motion of
a test particle influenced by the gravitational field, the radiation pressure
and the radiation drag force. The radiation field is modeled by photons
stemming out from a spherical and rigidly rotating radiation source, which
permits to calculate their impact parameter. They move along null geodesics of
the Kerr metric and hit the test particle at each time instant, modifying thus
its timelike geodesic trajectory (see Sec. II.1). We underline the ranges of
the model parameters, which will be useful in Sec. IV. This model can be
recast in Hamiltonian form, convenient to extract the dissipative PR
perturbations (see Sec. II.2). Finally, we discuss some a-priori indications
of possible chaotic behaviour in the general relativistic PR effect (see Sec.
II.3).

###  II.1 The model

We consider a central compact object, whose outside spacetime is described by
the Kerr metric with signature $(-,+,+,+)$. In geometrical units ($c=G=1$),
the line element of the Kerr spacetime,
$ds^{2}=g_{\alpha\beta}dx^{\alpha}dx^{\beta}$, in Boyer-Lindquist coordinates,
parameterized by mass $M$ (set equal to unity, $M=1$) and spin $a$, settled in
the equatorial plane $\theta=\pi/2$, reads as

$\displaystyle\mathrm{d}s^{2}$
$\displaystyle=\left(\frac{2}{r}-1\right)\mathrm{d}t^{2}-\frac{4a}{r}\mathrm{d}t\,\mathrm{d}\varphi+\frac{r^{2}}{\Delta}\,\mathrm{d}r^{2}+\rho\mathrm{d}\varphi^{2},$
(1)

where $\Delta\equiv r^{2}-2r+a^{2}$, and $\rho\equiv r^{2}+a^{2}+2a^{2}/r$. We
introduce the zero angular momentum observers (ZAMOs), whose adapted
orthonormal frame is given by111The hat over the indices indicates that the
corresponding vector or tensor quantity is calculated in the ZAMO frame. Bini
_et al._ (2009, 2011)

$\displaystyle\bm{e_{\hat{t}}}\equiv\bm{n}=\frac{\bm{\partial_{t}}-N^{\varphi}\bm{\partial_{\varphi}}}{N},\
\bm{e_{\hat{r}}}=\frac{\bm{\partial_{r}}}{\sqrt{g_{rr}}},\
\bm{e_{\hat{\varphi}}}=\frac{\bm{\partial_{\varphi}}}{\sqrt{g_{\varphi\varphi}}},$
(2)

where $\left\\{\bm{\partial_{t}},\ \bm{\partial_{r}},\
\bm{\partial_{\varphi}}\right\\}$ is the orthonormal frame adapted to the
static observer at infinity, $N=(-g^{tt})^{-1/2}$ is the time lapse function
and $N^{\varphi}=g_{t\varphi}/g_{\varphi\varphi}$ the spatial shift vector
field, whose explicit expressions are Bini  _et al._ (2011)

$N=\sqrt{\frac{\Delta}{\rho}},\qquad N^{\varphi}=-\frac{2a}{r\Delta}.$ (3)

The radiation field is constituted by a coherent flux of photons traveling
along null geodesics in the Kerr geometry. The related stress-energy tensor is
Bini  _et al._ (2009, 2011); De Falco  _et al._ (2019); Bakala  _et al._
(2019)

$T^{\mu\nu}=\Phi^{2}k^{\mu}k^{\nu}\,,\qquad k^{\mu}k_{\mu}=0,\qquad
k^{\mu}\nabla_{\mu}k^{\nu}=0,$ (4)

where $\bm{k}$ is the photon four-momentum field, and $\Phi$ is a parameter
linked to the radiation field intensity, whose explicit expression is given by
Bini  _et al._ (2011) 222The radial radiation function $\mathcal{R}_{\rm
rad}(r)$ can be equivalently written as $\mathcal{R}_{\rm
rad}(r)=(r^{2}+a^{2}-ab)^{2}-\Delta(a-b)^{2}$ De Falco  _et al._ (2019);
Bakala  _et al._ (2019).

$\Phi^{2}=\frac{\Phi_{0}^{2}}{\sqrt{\mathcal{R}_{\rm
rad}(r)}},\qquad\mathcal{R}_{\rm rad}(r)=rN|b\tan\beta|,$ (5)

where $\Phi_{0}$ is $\Phi$ evaluated at the emitting surface. Splitting
$\bm{k}$ with respect to the ZAMOs, we obtain Bini  _et al._ (2009, 2011)

$\displaystyle\bm{k}=E(n)[\bm{n}+\bm{\hat{\nu}}],\qquad\bm{\hat{\nu}}=\sin\beta\
\bm{e_{\hat{r}}}+\cos\beta\ \bm{e_{\hat{\varphi}}},$ (6)

where $\bm{\hat{\nu}},\ \beta,\ E(n)=E_{p}(1+bN^{\varphi})/N$ with
$E_{p}=-k_{t}$ is the conserved photon energy along its trajectory are
respectively the photon spatial unit relative velocity, the angle in the
azimuthal direction, and the photon energy, where all quantities are measured
in the ZAMO frame Bini  _et al._ (2009, 2011), The radiation field is governed
by the impact parameter $b$, associated with the emission angle $\beta$.

The photons of the radiation field are emitted from a spherical surface having
radius ${\mathrm{R_{\star}}}$ centered at the origin of the Boyer-Lindquist
coordinates, and rigidly rotating with angular velocity
$\Omega_{\mathrm{\star}}\geq 0$. Defined the event horizon $R_{\rm
H}=1+\sqrt{1-a^{2}}$ and the static limit $R_{\rm SL}=2$ radii in the
equatorial plane, we have that $R_{\star}\in(R_{\rm H}(a),\bar{R}_{\star}]$,
where $\bar{R}_{\star}<\infty$. Once $R_{\star}$ has been chosen, we want that
$\Omega_{\star}\in[\Omega_{\rm min},\Omega_{\rm
max}]=[\Omega_{-},\Omega_{+}]\cap[0,\Omega_{+}]$, where Bakala  _et al._
(2019)

$\Omega_{\pm}=\frac{-g_{t\varphi}\pm\sqrt{g_{t\varphi}^{2}-g_{\varphi\varphi}g_{tt}}}{g_{\varphi\varphi}}.$
(7)

The photon impact parameter is given by Bakala  _et al._ (2019)

$\displaystyle
b=-\left[\frac{\mathrm{g_{t\varphi}}+\mathrm{g_{\varphi\varphi}}\Omega_{\star}}{\mathrm{g_{tt}}+\mathrm{g_{t\varphi}}\Omega_{\star}}\right]_{r={\mathrm{R_{\star}}}},$
(8)

which in these premises ranges in $[b_{\rm min},b_{\rm
max}]\subseteq\mathbb{R}$ (see Fig. 2 in Ref. Bakala  _et al._ (2019), for
more details). The related photon angle in the ZAMO frame is Bakala  _et al._
(2019)

$\cos\beta=\frac{bN}{\sqrt{g_{\varphi\varphi}}(1+bN^{\varphi})},$ (9)

where $\beta\in[0,2\pi]$. For $\sin\beta>0$ ($\sin\beta<0$) we are considering
outgoing (ingoing) photons, see Ref. Bini  _et al._ (2011).

A test particle moves with a timelike four-velocity $\bm{U}$ and a spatial
three-velocity with respect to the ZAMOs, $\bm{\nu}(U,n)$, which both read as
Bini  _et al._ (2009, 2011); De Falco  _et al._ (2019); Bakala  _et al._
(2019)

$\displaystyle\bm{U}=\gamma[\bm{n}+\bm{\nu}],\qquad\bm{\nu}=\nu(\sin\alpha\bm{e_{\hat{r}}}+\cos\alpha\bm{e_{\hat{\varphi}}}),$
(10)

where $\gamma=1/\sqrt{1-||\bm{\nu}||^{2}}$ is the Lorentz factor,
$\nu=||\bm{\nu}||$ is the magnitude of the test particle spatial velocity
$\bm{\nu}(U,n)$, $\alpha$ is the azimuthal angle of the vector $\bm{\nu}$
measured clockwise from the positive $\bm{\hat{\varphi}}$ direction in the
$\bm{\hat{r}}-\bm{\hat{\varphi}}$ tangent plane in the ZAMO frame. The energy
absorbed by the test particle affected by the incoming photon is
$E(U)=-k_{\mu}U^{\mu}$, which can be related to the photon energy $E(n)$ in
the ZAMO frame through Bini  _et al._ (2009, 2011); De Falco  _et al._ (2019);
Bakala  _et al._ (2019)

$E(U)=\gamma E(n)[1-\nu\sin\psi\cos(\alpha-\beta)].$ (11)

We assume that the radiation test particle interaction occurs through Thomson
scattering, characterized by a constant momentum-transfer cross section
$\sigma$, independent from direction and frequency of the radiation field. The
radiation force is given by Bini  _et al._ (2009, 2011); De Falco  _et al._
(2019); Bakala  _et al._ (2019)

${\mathcal{F}}_{\rm(rad)}(U)^{\hat{\alpha}}=\sigma\,[\Phi
E(U)]^{2}\,\hat{\mathcal{V}}(k,U)^{\hat{\alpha}}\,.$ (12)

where the term $\tilde{\sigma}[\Phi E(U)]^{2}$ reads as Bini  _et al._ (2009,
2011); De Falco  _et al._ (2019); Bakala  _et al._ (2019)

$\displaystyle\frac{\tilde{\sigma}[\Phi E(U)]^{2}}{A}$ $\displaystyle=$
$\displaystyle\frac{\gamma^{2}(1+bN^{\varphi})^{2}\mathbb{A}^{2}}{N^{2}\sqrt{\mathcal{R}_{\rm
rad}(r)}},$ (13) $\displaystyle\mathbb{A}$ $\displaystyle=$
$\displaystyle\left[\gamma-\frac{p_{r}}{\sqrt{g_{rr}}}\sin\beta-\frac{p_{\varphi}}{\sqrt{g_{\varphi\varphi}}}\cos\beta\right].$
(14)

The term $A=\tilde{\sigma}[\Phi_{0}E_{p}]^{2}$ is the luminosity parameter,
which can be equivalently written as $A=L/L_{\rm Edd}\in[0,1]$, with $L$ the
emitted luminosity at infinity and $L_{\rm Edd}$ the Eddington luminosity. We
have that $\tilde{\sigma}=\sigma/m$, where $m$ is the test particle mass,
which for easing the notations we set equal to unity, $m=1$. The terms
$\hat{\mathcal{V}}(k,U)^{\hat{\alpha}}$ are the radiation field components,
which are Bini  _et al._ (2009, 2011); De Falco  _et al._ (2019); Bakala  _et
al._ (2019)

$\displaystyle\hat{\mathcal{V}}_{\hat{r}}$ $\displaystyle=$
$\displaystyle\frac{1}{\mathbb{A}}\left[\sin\beta-\frac{p_{r}}{\sqrt{g_{rr}}}\mathbb{A}\right],$
(15) $\displaystyle\hat{\mathcal{V}}_{\hat{\varphi}}$ $\displaystyle=$
$\displaystyle\frac{1}{\mathbb{A}}\left[\cos\beta-\frac{p_{\varphi}}{\sqrt{g_{\varphi\varphi}}}\mathbb{A}\right],$
(16) $\displaystyle\hat{\mathcal{V}}_{\hat{t}}$ $\displaystyle=$
$\displaystyle\frac{1}{\mathbb{A}}\left[1-\gamma\mathbb{A}\right],$ (17)

###  II.2 Hamiltonian formulation

The general relativistic PR effect in the Lagrangian formalism has been
already treated in De Falco  _et al._ (2018); De Falco and Battista (2019); De
Falco and Battista (2020); De Falco (2020), and we now pass to its Hamiltonian
formulation. In the geodesic case, we consider the mass shell constraint
$g^{\alpha\beta}p_{\alpha}p_{\beta}=-1$, where the momentum $p_{\alpha}$ is
canonically conjugate to $x^{\alpha}$ through the _Legendre transform_
$p_{\alpha}=g_{\alpha\beta}\dot{x}^{\beta}$. Here the dot stands for the
derivative with respect to the affine parameter $\tau$. Therefore, the
Hamiltonian is
$\mathcal{H}(\bm{p},\bm{x})=g^{\alpha\beta}p_{\alpha}p_{\beta}/2$ and, the
Hamilton equations are

$\dot{x}^{\mu}=\frac{\partial\mathcal{H}}{\partial
p_{\mu}},\qquad\dot{p}_{\mu}=-\frac{\partial\mathcal{H}}{\partial x^{\mu}}.$
(18)

Such formulation can be also extended to a dissipative system, where the
perturbations $\bm{f}(\bm{p},\bm{x})=(f_{1}^{\mu},f_{2,\mu})$ are not of
Hamiltonian type, therefore Eqs. (18) become

$\displaystyle\dot{x}^{\mu}=\frac{\partial\mathcal{H}}{\partial
p_{\mu}}+\epsilon
f_{1}^{\mu},\qquad\dot{p}_{\mu}=-\frac{\partial\mathcal{H}}{\partial
x^{\mu}}+\epsilon f_{2,\mu},$ (19)

where $\epsilon\ll 1$ is a small parameter.

####  II.2.1 General relativistic PR perturbations

The test particle velocity components are Bini  _et al._ (2009, 2011)

$\displaystyle
U^{\hat{r}}\equiv\frac{dr}{d\tau}=\frac{\gamma\nu\sin\alpha}{\sqrt{g_{rr}}},$
(20) $\displaystyle
U^{\hat{\varphi}}\equiv\frac{d\varphi}{d\tau}=\frac{\gamma\nu\cos\alpha}{\sqrt{g_{\varphi\varphi}}}-\frac{\gamma
N^{\varphi}}{N},$ (21) $\displaystyle
U^{\hat{t}}\equiv\frac{dt}{d\tau}=\frac{\gamma}{N},$ (22)

where $\tau$ is the affine parameter (proper time) along the test particle
trajectory, see Eqs. (10). In the PR effect case, the conjugate momenta
$p_{\mu}$ to the $x^{\mu}=(t,r,\varphi)$ are

$\displaystyle\frac{p_{r}}{\sqrt{g_{rr}}}=\gamma\nu\sin\alpha,\quad\frac{p_{\varphi}}{\sqrt{g_{\varphi\varphi}}}=\gamma\nu\cos\alpha,\quad
p_{t}=\frac{\gamma}{N}.$ (23)

In such formalism, we have that $\nu$ and $\gamma$ read as

$\nu=\sqrt{\frac{\left(\frac{p_{r}^{2}}{g_{rr}}+\frac{p_{\varphi}^{2}}{g_{\varphi\varphi}}\right)}{1+\left(\frac{p_{r}^{2}}{g_{rr}}+\frac{p_{\varphi}^{2}}{g_{\varphi\varphi}}\right)}},\quad\gamma=\sqrt{1+\frac{p_{r}^{2}}{g_{rr}}+\frac{p_{\varphi}^{2}}{g_{\varphi\varphi}}}.$
(24)

Using the radiation force components (15) – (17), we obtain
$\tilde{F}_{\mu}=\tilde{\sigma}[\Phi
E(U)]^{2}\hat{\mathcal{V}}_{\hat{\mu}}/A$, where $\epsilon=A\equiv L/L_{\rm
Edd}\ll 1$, namely low luminosities. The PR dissipative perturbations are
$(f^{\mu}_{1},f_{2,\mu})=(0,\tilde{F}_{\mu})$. We note that $f^{\mu}_{1}=0$,
because the radiation field, including radiation pressure and PR drag force,
affects only the accelerations and not the velocity components, see Refs. Bini
_et al._ (2011) for more details.

###  II.3 A-priori indications of chaotic behavior

The main motivations for the present study are explained in this section.
During the investigation of the general relativistic PR effect, a series of
a-priori indications of possible chaotic dynamics can be found:

  * •

it is a dissipative and non-linear dynamical system in GR De Falco and
Battista (2020), which makes the Kerr geodesic motion not _integrable_
Strogatz (2000); Mori  _et al._ (2013);

  * •

it has been analytically and numerically confirmed that such effect generally
behaves as a _forced harmonic oscillator_ Bini  _et al._ (2009, 2011) endowed
with a non-linear driven force (close to a Duffing oscillator), responsible to
potentially create _resonance effects_ Tabor (1989);

  * •

a test particle under the general relativistic PR effect can end its motion
either on the critical hypersurface or escaping at infinity. It has been
already formally proved that the critical hypersurface behaves as a _stable
attractor_ De Falco and Bakala (2020), and the same holds also for the spatial
infinity (never returning back);

  * •

as proved in Bini  _et al._ (2009), such effect admits positive _Lyapunov
exponents_ , which measure the mean rate of exponential separation of
neighboring trajectories Tabor (1989). This is an useful index indicating that
a dynamical system shows _sensitive dependence on the initial conditions_.
This propriety has been further confirmed by numerical simulations.

##  III Homoclinic orbits

The notion of homoclinic orbits for a dynamical system is based on the
research of _recurrent invariant sets_ $\Lambda$ Wiggins (1988); Guckenheimer
and Holmes (2002), such as fixed points, periodic orbits, or $n$-dimensional
invariant tori. The set of all trajectories approaching an invariant set
$\Lambda$ asymptotically in the infinite future (past) is a submanifold of the
phase space termed _stable (unstable) manifold of $\Lambda$_, usually
indicated by $W^{s}(\Lambda)$ ($W^{u}(\Lambda)$). An invariant set $\Lambda$
possessing both stable and unstable manifolds is called _hyperbolic_ 333The
given definition has a clear dynamical meaning. However, the rigorous
definition of an hyperbolic point $p$ for a $C^{1}$ vector field
$\bm{F}:\mathbb{R}^{n}\to\mathbb{R}^{n}$ is the following: $p$ is a critical
point for $\bm{F}$, i.e., $\bm{F}(p)=0$, and the Jacobian matrix of $\bm{F}$
at $p$, $\bm{J}=(\nabla\bm{F})(p)$, has no eigenvalues with zero real parts
Wiggins (1988); Guckenheimer and Holmes (2002). The stable (unstable) manifold
of $p$ consists of points $q$ such that $\phi_{t}(q)\to p$ as $t\to+\infty$
($t\to-\infty$), where $\phi_{t}$ is the _flow_ associated with $\bm{F}$..

A trajectory is defined to be homoclinic to a hyperbolic invariant set
$\Lambda$ if it approaches $\Lambda$ in the infinite future as in the infinite
past, i.e., $W^{u}(\Lambda)\cap W^{s}(\Lambda)\cap\Lambda\neq\emptyset$
Wiggins (1988); Guckenheimer and Holmes (2002). _Therefore, for determining
the class of the homoclinic orbits of a dynamical system we need to identify
the intersections of their stable and unstable manifolds on their hyperbolic
invariant sets_.

###  III.1 Homoclinic orbits in the equatorial plane of Kerr spacetime

We consider the following dynamical system represented by the motion of a
timelike test particle governed only by gravity and no other perturbing
effects in the equatorial plane of the Kerr metric Misner  _et al._ (1973);
Levin and Perez-Giz (2009)

$\displaystyle\dot{t}=\frac{r\rho
E-2aL_{z}}{r\Delta},\quad\dot{r}=\pm\frac{\sqrt{R(r)}}{r^{2}},\quad\dot{\varphi}=\dot{\varphi}(a,r),$
(25)

where $\dot{\varphi}(a,r)=[2aE+L_{z}(r-2)]/(r\Delta)$, 444It is possible to
factorize $R(r)=-(1-E^{2})r^{4}+2r^{3}-[a^{2}(1-E^{2})+L_{z}^{2}]r^{2}+2(aE-
L_{z})^{2}r$ as reported in Eq. (26).

$\displaystyle R(r)=-(1-E^{2})r(r-r_{u})^{2}(r-r_{a}),$ (26)

$E=-p_{t}$, and $L_{z}=p_{\varphi}$ are respectively the energy and angular
momentum with respect to the $\bm{z}$-axis (orthogonal to the equatorial
plane) conserved along the test particle trajectory, and $r_{u}$ and $r_{a}$
are respectively the periastron and apastron radii of the homoclinic orbit.
Throughout the paper the signs $\pm$ refers to prograde and retrograde orbits,
respectively.

The invariant sets are the circular orbits (defined by the conditions $R(r)=0$
and $dR(r)/dr=0$), while the hyperbolic invariant sets coincide with the
unstable circular orbits (defined by circular orbit condition, and
$d^{2}R(r)/dr^{2}<0$, which corresponds to the maximum of $dR/dr=0$). Among
these trajectories, the homoclinic orbits are the unstable circular orbits
energetically bounded ($E<1$) Levin and Perez-Giz (2009), that we describe
through the periastron and apastron radii $(r_{p},r_{a})$. _Homoclinic orbits
are in a one-to-one correspondence with bound energy values $E<1$, and
therefore constitute a one-parameter family specified by the (periastron)
radius $r_{u}=r_{p}$_.

Figure 1: Homoclinic orbit (left panel) and its phase portrait (right panel),
for $a=0.1$, $r_{u}=4.40$. The other parameters are $r_{a}=13.61$, $E=0.95$,
and $L_{z}=3.52$. The red dashed line represents the circular orbit of radius
$r_{u}$ centered in $(0,0)$ towards which the homoclinic orbit moves. The
hyperbolic fixed red point in the phase space has coordinates
$(r,p_{r})=(r_{u},0)$.

The one-parameter family of homoclinic orbits in the equatorial plane of the
Kerr spacetime $\mathcal{O}^{\rm hc}(r_{u})$, see Fig. 1 as an example555To
plot the homoclinic orbit in the equatorial plane of Kerr metric, we use Eq.
(26c) in Ref. Levin and Perez-Giz (2009) for describing the azimuthal
coordinate $\varphi$, while the radial coordinate $r$ ranges in
$[r_{u},r_{a}]$., is characterized by Levin and Perez-Giz (2009)

$\displaystyle E$ $\displaystyle=$
$\displaystyle\frac{r_{u}^{3/2}-2r_{u}^{1/2}\pm
a}{r_{u}^{3/4}\sqrt{r_{u}^{3/2}-3r_{u}^{1/2}\pm 2a}}<1,$ (27) $\displaystyle
L_{z}$ $\displaystyle=$ $\displaystyle\frac{r_{u}^{2}\mp
2ar_{u}^{1/2}+a^{2}}{r_{u}^{3/4}\sqrt{r_{u}^{3/2}-3r_{u}^{1/2}\pm 2a}},$ (28)
$\displaystyle r_{a}$ $\displaystyle=$ $\displaystyle\frac{2(aE-
L_{z})^{2}}{r_{u}^{2}(1-E^{2})}\equiv\frac{2r_{u}(a\mp\sqrt{r_{u}})^{2}}{r_{u}^{2}-4r_{u}\pm
4a\sqrt{r_{u}}-a^{2}},$ (29)

where $r_{u}$ ranges between the innermost bound circular orbit (IBCO), and
the innermost stable circular orbit (ISCO), i.e., $r_{u}\in[r_{\rm
IBCO},r_{\rm ISCO}]$ (see Fig. 2), with

$\displaystyle r_{\rm IBCO}$ $\displaystyle=$ $\displaystyle 2\mp
a+2\sqrt{1\mp a},$ (30) $\displaystyle r_{\rm ISCO}$ $\displaystyle=$
$\displaystyle 3+Z_{2}\mp\sqrt{(3-Z_{1})(3+Z_{1}+2Z_{2})},$ (31)
$\displaystyle Z_{1}$ $\displaystyle=$ $\displaystyle
1+{}^{3}\sqrt{1-a^{2}}\left[{}^{3}\sqrt{1+a}+{}^{3}\sqrt{1-a}\right],$ (32)
$\displaystyle Z_{2}$ $\displaystyle=$ $\displaystyle\sqrt{3a^{2}+Z_{1}^{2}}.$
(33)

Figure 2: The shaded area, delimited by IBCO and ISCO radii, is where $r_{u}$
can range in terms of the spin $a$. The vertical dashed red line represent the
Schwarzschild limit $(a=0)$, the horizontal dashed green line is the static
limit radius, and the continuous blue line is the event horizon radius.

##  IV Melnikov integral

The Melnikov method is a powerful mathematical tool to identify the occurrence
of chaos in 2D and even higher-dimensional dynamical systems affected by
Hamiltonian periodic or non-Hamiltonian perturbations Wiggins (1988);
Guckenheimer and Holmes (2002); Holmes and Marsden (1982a, b).

Let $\bm{\Phi}:\mathbb{R}^{2n}\to\mathbb{R}^{2n}$ be an Hamiltonian integrable
dynamical system, which for the Liouville theorem is an area-preserving map in
the phase space, possessing a hyperbolic fixed point $P$ and a homoclinic
orbit $\mathcal{O}$. Such a Hamiltonian system is affected by dissipative
perturbations like Eqs. (19). In these hypothesis, the Melnikov method goes in
search of _homoclinic tangles_ Guckenheimer and Holmes (2002); Wiggins (1988);
Bombelli and Calzetta (1992), see Fig. 3.

Figure 3: Example of homoclinic tangle.

In such structures, once $W^{s}(P)$ and $W^{u}(P)$ intersect each other, they
will continue to intersect infinitely in a discrete number of points, like
$\left\\{Q,Q^{\prime},Q^{\prime\prime},Q^{\prime\prime\prime},\cdots,P\right\\}$.
They cannot touch the same point twice, otherwise they will be trapped in a
cycle and will not reach the point $P$, and $P$ is not touched in a finite
number of steps, since $P$ is a fixed point (has no image or pre-image of a
point other than itself). Since $\bm{\Phi}$ is an area-preserving map, the
areas formed by the intersection of $W^{s}(P)$ and $W^{u}(P)$ (i.e.,
$\mathcal{A},\mathcal{A}^{\prime},\mathcal{A}^{\prime\prime},\cdots$) are all
equivalent. The _Smale-Birkhoff theorem_ claims that the dynamics produced by
$W^{s}(P)$ and $W^{u}(P)$ in approaching the point $P$ gives rise to the
_Smale horseshoe’s map_ , which is a _chaotic map_ Guckenheimer and Holmes
(2002); Wiggins (1988).

In order to find the homoclinic tangles, we have to find a time $t_{0}$ such
that $W^{s}(P_{t_{0}})$ and $W^{u}(P_{t_{0}})$ intersect transversally. To
this end, we fix an arbitrary initial time $t_{0}$ or Poincaré section (where
we follow the dynamics) corresponding to the hyperbolic fixed point
$P_{t_{0}}$ (see Fig. 4).

Figure 4: Sketch of Poincaré sections for unperturbed (left cartoon) and
perturbed (right cartoon) dynamical systems.

We consider a point $Q\in\mathcal{O}$, and then we define the distance from
$W^{s}(P_{t_{0}})$ to $W^{u}(P_{t_{0}})$ along a transversal direction to
$\mathcal{O}$ in $Q$, which intersects $W^{s}(P_{t_{0}})$ and
$W^{u}(P_{t_{0}})$ respectively in $Q^{s}$ and $Q^{u}$, see Fig. 5,

Figure 5: Cartoon to explain the Melnikov integral.

$d(t_{0})\approx\epsilon\frac{\mathcal{M}(t_{0})}{||\nabla\mathcal{H}(P_{t_{0}})||}+O(\epsilon^{2}),$
(34)

where $||\nabla\mathcal{H}(P_{t_{0}})||\neq 0$. Here $\mathcal{M}(t_{0})$ is
the _Melnikov integral_ defined as Guckenheimer and Holmes (2002); Wiggins
(1988); Holmes and Marsden (1982a, b); Bombelli and Calzetta (1992)

$\displaystyle\mathcal{M}(t_{0})$
$\displaystyle=\int_{-\infty}^{+\infty}\left\\{\mathcal{H},\bm{f}\right\\}\
dt$ (35)
$\displaystyle=\int_{-\infty}^{+\infty}\sum_{\mu=1}^{n}\left(\frac{\partial\mathcal{H}}{\partial
p_{\mu}}f_{2,\mu}+\frac{\partial\mathcal{H}}{\partial
x^{\mu}}f_{1}^{\mu}\right)\ dt,$

where $\left\\{\cdot,\cdot\right\\}$ are the Poisson brackets, $\bm{f}$ are
the perturbations of Eqs. (19), and the integral is taken along the
unperturbed homoclinic orbit $\mathcal{O}$. Depending on the values assumed by
the Melnikov integral, we have:

  * •

if $\mathcal{M}(t_{0})$ admits odd order zeros, there is chaos;

  * •

if $\mathcal{M}(t_{0})$ is bounded away from zero, there is no occurrence of
chaos in the perturbed dynamics;

  * •

if $\mathcal{M}(t_{0})$ is identically zero or admits even order zeros, the
method cannot predict anything.

###  IV.1 Application to the general relativistic PR effect perturbing the
equatorial Kerr dynamics

We apply the Melnikov method to the general relativistic PR effect, where the
perturbations are $\bm{f}=(0,\tilde{F}_{\mu})$, see Sec. II.2.1. The Melnikov
integral (35) reads as

$\displaystyle\mathcal{M}\equiv\mathcal{M}(a,r_{u},R_{\star},\Omega_{\star};t_{0})$
(36)
$\displaystyle=\int_{-\infty}^{+\infty}\left[\left(g^{tt}p_{t}+g^{t\varphi}p_{\varphi}\right)_{t=t_{0}}\tilde{F}_{t}(t-t_{0})\right.$
$\displaystyle\left.+\left(g^{\varphi\varphi}p_{\varphi}+g^{t\varphi}p_{t}\right)_{t=t_{0}}\tilde{F}_{\varphi}(t-t_{0})\right.$
$\displaystyle\left.+\left(g^{rr}p_{r}\right)_{t=t_{0}}\tilde{F}_{r}(t-t_{0})\right]\
dt.$

This integral is evaluated along the homoclinic orbit
$\mathcal{O}\in\mathcal{O}^{\rm hc}(r_{u})$ at the time $t_{0}$. It is
important to note that if there exists an intersection for some $t_{0}$, then
there will be one for every $t_{0}$ Bombelli and Calzetta (1992). Considering
$p_{r}=\dot{r}g_{rr}$, $p_{\varphi}=L_{z}$, and $p_{t}=-E$ (see Eqs. (25),
(28), and (27), respectively), the explicit expressions of
$\hat{\mathcal{V}}^{\hat{\mu}}$ (see Eqs. (13) – (17)), and passing from the
coordinate time $t$ to the coordinate radius $r$ integration, cf. Eqs. (25),
we have

$\mathcal{M}=2\int_{r_{u}}^{r_{a}}\Psi_{1}\Psi_{2}dr,$ (37)

where

$\displaystyle\frac{\dot{t}}{\dot{r}}$ $\displaystyle=$
$\displaystyle\frac{r(r\rho E-2aL_{z})}{\Delta\sqrt{R(r)}},$ (38)
$\displaystyle\Psi_{1}$ $\displaystyle=$
$\displaystyle\frac{\gamma^{2}(1+bN^{\varphi})^{2}}{N^{2}\sqrt{R_{\rm
rad}}}\frac{\mathbb{A}}{r_{0}\Delta_{0}}\frac{\dot{t}}{\dot{r}},$ (39)
$\displaystyle\Psi_{2}$ $\displaystyle=$
$\displaystyle(\rho_{0}r_{0}E-2aL_{z})(1-\gamma\mathbb{A})$
$\displaystyle+\left[(r_{0}-2)L_{z}+2aE\right]\left(\cos\beta-\mathcal{B}\mathbb{A}\right)$
$\displaystyle+\frac{\sqrt{R_{0}}\Delta_{0}}{r_{0}}\left(\sin\beta-\mathcal{A}\mathbb{A}\right).$

In order to simplify the notations we have defined

$\mathcal{A}=\frac{p_{r}}{\sqrt{g_{rr}}}\equiv\frac{1}{r}\sqrt{\frac{R}{\Delta}}\geq
0,\quad\mathcal{B}=\frac{L_{z}}{\sqrt{g_{\varphi\varphi}}}\equiv\frac{L_{z}}{\sqrt{\rho}}>0.$
(41)

This implies that Eqs. (24) – (14) reads respectively as

$\gamma=\sqrt{1+\mathcal{A}^{2}+\mathcal{B}^{2}},\quad\mathbb{A}=\gamma-\mathcal{A}\sin\beta-\mathcal{B}\cos\beta.$
(42)

The quantities with a subscript zero means that they are evaluated at the
initial time $t_{0}$. Since the general relativistic PR effect dynamics does
not depend explicitly on the time, we can set without loss of generality and
for simplicity $t_{0}=0$. We will show that $\mathcal{M}$ has not a defined
sign for all parameters ranging in their intervals. To achieve this goal, we
will study the signs of each components of the integrating function by
performing either analytical calculations or numerical simulations, where the
functions are difficult to handle analytically.

Resuming what has been discussed in the previous sections, we have that the
set of parameters
$\left\\{a,r_{u},R_{\star},\Omega_{\star},r_{0},r,\beta\right\\}$ range over
the following intervals

$\displaystyle a\in[0,1),\qquad r_{u}\in[r_{\rm IBCO}(a),r_{\rm ISCO}(a)],$
(43) $\displaystyle R_{\star}\in(r_{\rm
H}(a),\bar{R}_{\star}],\qquad\Omega_{\star}\in[\Omega_{\rm min},\Omega_{\rm
max}],$ $\displaystyle r_{0}\in(r_{\rm H}(a),\bar{R}],\quad
r\in[r_{u},r_{a}(a,r_{u})],\quad\beta\in[0,2\pi],$

where $\bar{R}_{\star},\bar{R}$ are finite values. In addition, we know

$\frac{\sqrt{3}}{3}\leq E<1,\quad\frac{2\sqrt{3}}{3}\leq L_{z}\leq
4,\quad\gamma\geq 1.$ (44)

We prove that $\Psi_{1}>0$, because

$\frac{\gamma^{2}(1+bN^{\varphi})^{2}}{N^{2}\sqrt{R_{\rm
rad}}}>0,\quad\frac{\mathbb{A}}{r_{0}\Delta_{0}}>0,\quad\frac{\dot{t}}{\dot{r}}>0.$
(45)

The first term is composed by positive quantities. The second term is
positive, because $\mathbb{A}>0$ (see Appendix A.1 and Eqs. (42)). The term
$\dot{t}/\dot{r}$ has a positive denominator, as well as the numerator (see
the proof in Appendix A.2).

We prove that $\Psi_{2}$ has not a defined sign. Indeed, we have that
$\rho_{0}r_{0}E-2aL_{z}$ is equal to the numerator of $\dot{t}/\dot{r}$, which
is non-negative, $\sqrt{R_{0}}\Delta_{0}/r_{0}$ is non-negative and

$(r_{0}-2)L_{z}+2aE>(1-2)\frac{2\sqrt{3}}{3}+\frac{2\sqrt{3}}{3}=0.$ (46)

Figure 6: Ranges of the $f,g,h$ functions (blue surfaces). The orange surface
corresponds to zero value of the functions.

The signs of $\Psi_{2}$ depends therefore only by (defining
$f=\cos\beta-\mathcal{B}\mathbb{A}$, $g=\sin\beta-\mathcal{A}\mathbb{A}$,
$h=1-\gamma\mathbb{A}$)

$\displaystyle f(\mathcal{A},\mathcal{B},\beta)$
$\displaystyle=\cos\beta(1+\mathcal{B}^{2})+\mathcal{B}(\mathcal{A}\sin\beta-\gamma),$
(47) $\displaystyle g(\mathcal{A},\mathcal{B},\beta)$
$\displaystyle=\sin\beta(1+\mathcal{A}^{2})+\mathcal{A}(\mathcal{B}\cos\beta-\gamma),$
(48) $\displaystyle h(\mathcal{A},\mathcal{B},\beta)$
$\displaystyle=-(\mathcal{A}^{2}+\mathcal{B}^{2})+\mathcal{A}\sqrt{1+\mathcal{A}^{2}+\mathcal{B}^{2}}.$
(49)

After having found the ranges of $\mathcal{A},\mathcal{B}$ (see Appendix A.3),
we see that for $\beta\in[0,2\pi]$ the functions $f,g,h$ does not have a
definite sign as can be seen in Fig. 6

This result implies that the Melnikov integral may admit zero values. Due to
the behaviors of the $f,g,h$ functions it is very difficult to analytically
describe the set of parameters for which the Melnikov integral vanishes.
Therefore, we resort to numerical simulations to investigate this issue. We
develop a code in Mathematica 12.1.1.0, which permits to numerically check
whether there are values of $r_{0}$ such that vanish the Melnikov integral. We
calculate also the derivative of the Melnikov integral with respect to
$r_{0}$, proving that it is non-zero at the value of $r_{0}$ for which the
Melnikov integral vanishes, thus assuring that the zero is simple and chaos
occurs (see Theorem 4.5.2 in Ref. Guckenheimer and Holmes (2002), for more
details) 666We note that our dynamical system is autonomous with respect to
the time $t$. Since the model is set in the equatorial plane all the functions
depend only by $r(t)$. Therefore, we should find the value of $r_{0}$ such
that $\mathcal{M}(r_{0})=0$. In addition the condition
$\frac{\partial\mathcal{M}}{\partial t_{0}}\neq 0,$ can be substituted by
$\frac{\partial\mathcal{M}}{\partial r_{0}}\left[\frac{\partial r(t)}{\partial
t}\right]_{t=t_{0}}.$ .

In our numerical simulations we found chaotic behavior for certain parameter
values. In Fig. 7, we display the region of the parameter space where chaos
occurs.

Figure 7: Parameter space $(a,r_{u})$ delimited by the curves $r_{\rm IBCO}$
and $r_{\rm ISCO}$ for the photon impact parameter value $b=3$. The blue dots
are the values found from our numerical simulations in Mathematica,
corresponding to chaotic dynamics.

We performed several numerical simulations, and we found that for photon
impact parameter $b=3$ there is the occurrence of chaos for almost each spin
value, while for $b=0,1,2$ chaos is not present. _Therefore, we conclude that
for radial radiation field, $b=0$, chaos does not manifest_. In addition, we
checked that for values closer to $b=3$, chaos still reveals its presence.
Once $(a^{*},r_{u}^{*})$, and the radius $r_{0}$ for having chaotic dynamics
have been found, we calculate the initial conditions on the test particle
trajectory by calculating first
$\mathcal{A}^{*}=\mathcal{A}(a^{*}r_{u}^{*},r_{0})$ and
$\mathcal{B}^{*}=\mathcal{B}(a^{*}r_{u}^{*},r_{0})$, see Eqs. (42), and then
we can calculate the test particle initial velocity conditions

$\nu_{0}=\sqrt{\frac{\mathcal{A}^{*}{}^{2}+\mathcal{B}^{*}{}^{2}}{1+\mathcal{A}^{*}{}^{2}+\mathcal{B}^{*}{}^{2}}},\qquad\alpha_{0}=\arctan\left(\frac{\mathcal{A}^{*}{}}{\mathcal{B}^{*}{}}\right),$
(50)

obtained by employing Eqs. (23), (24), and (42). As an example we plot in Fig.
8 a chaotic orbit together with a normal dynamics to highlight the different
behaviours.

Figure 8: Trajectories of two distinct test particles for $a=0.5$, $A=L/L_{\rm
Edd}=0.1$, $b=3.1M$ starting both at $r_{0}=3.17M$ with angular velocity
$\alpha_{0}=0.11$, but with different initial velocities, $\nu_{0}=0.67$
(black and chaotic orbit) and $\nu_{0}=0.60$ (green orbit). The red circle is
the critical hypersurface for the black orbit located at $r_{\rm crit}=5.88M$,
while the critical hypersurface for the green one is located at $r_{\rm
crit}=1.88M$, very close to the event horizon $R_{\rm H}(0.5)=1.87M$.

##  V Conclusions

We have analysed the general relativistic PR effect in the equatorial plane of
Kerr spacetime from a dynamical system point of view (see Sec. II). We have
employed the Melnikov method to investigate whether it admits chaotic
behaviors. The idea behind such investigation relies on the existence of some
_a-priori_ indications of chaos, which are: $(i)$ non-integrability, $(ii)$
presence of a stable attractor (critical hypersurface), $(iii)$ strong analogy
with a forced harmonic oscillator (Duffing equations), $(iv)$ sensitive
dependence on the initial conditions (see Sec. II.3).

The Melnikov method is based on the knowledge of unperturbed Hamiltonian Kerr
metric, general relativistic PR dissipative perturbations (see Sec. II.2.1),
and homoclinic orbits in the equatorial plane of Kerr spacetime, parametrized
by the periastron $r_{u}$ (see Sec. III.1). The aim is for the existence of
homoclinic tangles in the phase space, whose dynamics reproduce that of the
Smale horseshoe’s map, which is a chaotic map. This reduces to determining
whether the Melnikov integral, see Eq. (36), admits zeros in terms of its
parameters $(a,r_{u},b(R_{\star},\Omega_{\star}))$ and initial condition
$r_{0}$, see Sec. IV. We proved that the Melnikov integral admits simple zeros
and therefore chaos is present in the dynamics of the general relativistic PR
effect in the Kerr equatorial plane for low luminosities.

This result is relevant, because we discovered that the general relativistic
PR effect can admit chaotic orbits for a suitable range of parameters provided
by the Melnikov method. Although several numerical simulations of the PR
trajectories have been performed in the literature (see Refs. Bini  _et al._
(2009, 2011); De Falco (2019), for further details), to our knowledge the
existence of chaotic trajectories has never been reported in previous works on
the PR effect. These configurations are useful for astrophysical purposes,
because they can be exploited as a valuable tool for lighting up the compact
object around which a test particle orbits, being thus a further source of
information Levin (1999). Instead, for observational goals it could be better
to avoid chaotic orbits for not compromising the detections of phenomena
related to the general relativistic PR effect Cornish (2001).

## Acknowledgements

V.D.F. thanks Gruppo Nazionale di Fisica Matematica of Istituto Nazionale di
Alta Matematica for support. W.B. acknowledges support from Gruppo Nazionale
per l’Analisi Matematica, la Probabilità e le loro Applicazioni of Istituto
Nazionale di Alta Matematica. V.D.F. and W.B. are grateful to Dr. Paolo
Giulietti for fruitful discussions on the Melnikov method and on Chaos Theory.

##  Appendix A signs and ranges of parameters

This Appendix is devoted to prove the sign or range of some parameters,
appearing in the terms $\Psi_{1},\Psi_{2}$.

###  A.1 Sign and range of $\mathbb{A}$

The expression of $\mathbb{A}$, see Eq. (14), can be equivalently written in
terms of Eqs. (42) as

$\mathbb{A}=\sqrt{1+\mathcal{A}^{2}+\mathcal{B}^{2}}-\mathcal{A}\sin\beta-\mathcal{B}\cos\beta,$
(51)

where $\gamma=\sqrt{1+\mathcal{A}^{2}+\mathcal{B}^{2}}$, cf. Eqs. (24). We
claim that $\mathbb{A}>0$. If $\sin\beta,\ \cos\beta\leq 0$ this derives
immediately from (51). Assuming that $\sin\beta>0$ and $\cos\beta\leq 0$ we
get

$\mathbb{A}\geq\sqrt{1+\mathcal{A}^{2}+\mathcal{B}^{2}}-\mathcal{A}\sin\beta>\sqrt{1+\mathcal{A}^{2}}-\mathcal{A}>0\,$

and the same argument holds exchanging the role of $\cos\beta$ and
$\sin\beta$. Finally, if $\cos\beta,\ \sin\beta>0$, we obtain

$1+\mathcal{A}^{2}+\mathcal{B}^{2}>(\mathcal{A}\sin\beta+\mathcal{B}\cos\beta)^{2},$
(52)

from which we obtain

$1+(\mathcal{A}\cos\beta-\mathcal{B}\sin\beta)^{2}>0.$ (53)

###  A.2 Sign of $\dot{t}/\dot{r}$

The term $\dot{t}/\dot{r}$ is the Jacobian of coordinate transformation, cf.
Eq. (38), therefore it must be non-zero. Since the denominator is always
positive, we focus only on the sign of the numerator, which is estimated
through

$r\rho E-2aL_{z}\geq
r\rho\frac{\sqrt{3}}{3}-\frac{4\sqrt{3}}{3}>\frac{4\sqrt{3}}{3}-\frac{4\sqrt{3}}{3}=0.$
(54)

###  A.3 Ranges of $\mathcal{A}$ and $\mathcal{B}$

The analytical expression of $\mathcal{A}$ is, see Eqs. (26) – (42),

$\mathcal{A}=\sqrt{\frac{(1-E^{2})(r-r_{u})^{2}(r_{a}-r)}{r\Delta}}.$ (55)

where the numerator has a maximum at $r_{m}=(r_{u}+2r_{a})/3$, while the
denominator is a monotone increasing function for $r\in[r_{u},r_{a}]$. We
numerically checked that $\mathcal{A}$ attains its maximum for $r_{u}=R_{\rm
IBCO}(a)$. It becomes a constant function independent from $a$, such that
$\mathcal{A}(a,R_{\rm IBCO}(a))\approx 0.7$, showing thus that
$\mathcal{A}<1$.

The range of $\mathcal{B}$ can be obtained through

$0\leq\mathcal{B}\leq\frac{L_{z}}{r_{u}}<2.$ (56)

## References

  * Wiggins (1988) S. Wiggins,  _Global Bifurcations and Chaos: Analytical Methods_, Applied Mathematical Sciences Series No. v. 73 (Springer-Verlag, 1988).
  * Ott (2002) E. Ott,  _Chaos in Dynamical Systems_ (Cambridge University Press, 2002).
  * Devaney (2018) R. Devaney,  _An introduction to chaotic dynamical systems_ (CRC Press, 2018).
  * Hobill  _et al._ (1994) D. Hobill, A. Burd, and A. Coley,  _Deterministic Chaos in General Relativity_, Nato Science Series B: (Springer US, 1994).
  * Arnol’d and Avez (1989) V. Arnol’d and A. Avez,  _Ergodic Problems of Classical Mechanics_, Advanced book classics (Addison-Wesley, 1989).
  * Contopoulos (1990) G. Contopoulos, Proceedings of the Royal Society of London Series A 431, 183 (1990).
  * Contopoulos (1991) G. Contopoulos, Proceedings of the Royal Society of London Series A 435, 551 (1991).
  * Wanex (2002) L. F. Wanex,  _Chaotic amplification in the relativistic restricted three-body problem_ , Ph.D. thesis, UNIVERSITY OF NEVADA, RENO (2002).
  * Bombelli and Calzetta (1992) L. Bombelli and E. Calzetta, Classical and Quantum Gravity 9, 2573 (1992).
  * Suzuki and Maeda (1997) S. Suzuki and K.-I. Maeda, PRD 55, 4848 (1997), arXiv:gr-qc/9604020 [gr-qc] .
  * Lukes-Gerakopoulos (2018) G. Lukes-Gerakopoulos, in  _Fifteenth Marcel Grossmann Meeting - MG15_ (2018) pp. 1960–1965, arXiv:1606.09430 [gr-qc] .
  * Cornish (2001) N. J. Cornish, Phys. Rev. D 64, 084011 (2001), arXiv:gr-qc/0106062 [gr-qc] .
  * Cornish and Levin (2002) N. J. Cornish and J. Levin, Phys. Rev. Lett. 89, 179001 (2002), arXiv:gr-qc/0207020 [gr-qc] .
  * Cornish and Levin (2003) N. J. Cornish and J. Levin, Phys. Rev. D 68, 024004 (2003).
  * Belinskij  _et al._ (1970a) V. A. Belinskij, E. M. Lifshits, and I. M. Khalatnikov, Uspekhi Fizicheskikh Nauk 102, 463 (1970a).
  * Belinskij  _et al._ (1970b) V. A. Belinskij, I. M. Khalatnikov, and E. M. Lifshits, Advances in Physics 19, 525 (1970b).
  * Barrow and Sirousse-Zia (1989) J. D. Barrow and H. Sirousse-Zia, PRD 39, 2187 (1989).
  * Burd  _et al._ (1991) A. B. Burd, N. Buric, and R. K. Tavakol, Classical and Quantum Gravity 8, 123 (1991).
  * Contopoulos  _et al._ (1999) G. Contopoulos, N. Voglis, and C. Efthymiopoulos, Celestial Mechanics and Dynamical Astronomy 73, 1 (1999).
  * Calzetta and El Hasi (1993) E. Calzetta and C. El Hasi, Classical and Quantum Gravity 10, 1825 (1993), arXiv:gr-qc/9211027 [gr-qc] .
  * Aydiner (2016) E. Aydiner, arXiv e-prints , arXiv:1610.07338 (2016), arXiv:1610.07338 [gr-qc] .
  * Poynting (1903) J. H. Poynting, Monthly Notices of the Royal Astronomical Society 64, 1 (1903).
  * Robertson (1937) H. P. Robertson, Monthly Notices of the Royal Astronomical Society 97, 423 (1937).
  * Bini  _et al._ (2009) D. Bini, R. T. Jantzen, and L. Stella, Classical and Quantum Gravity 26, 055009 (2009), arXiv:0808.1083 [gr-qc] .
  * Bini  _et al._ (2011) D. Bini, A. Geralico, R. T. Jantzen, O. Semerák, and L. Stella, Classical and Quantum Gravity 28, 035008 (2011), arXiv:1408.4945 [gr-qc] .
  * Bini  _et al._ (2015) D. Bini, A. Geralico, and A. Passamonti, MNRAS 446, 65 (2015), arXiv:1410.3099 [astro-ph.HE] .
  * De Falco  _et al._ (2019) V. De Falco, P. Bakala, E. Battista, D. Lančová, M. Falanga, and L. Stella, Phys. Rev. D 99, 023014 (2019).
  * Bakala  _et al._ (2019) P. Bakala, V. De Falco, E. Battista, K. Goluchová, D. Lančová, M. Falanga, and L. Stella, PRD 100, 104053 (2019).
  * Wielgus (2019) M. Wielgus, MNRAS 488, 4937 (2019), arXiv:1907.11268 [astro-ph.HE] .
  * De Falco  _et al._ (2020) V. De Falco, P. Bakala, and M. Falanga, Phys. Rev. D 101, 124031 (2020), arXiv:2006.01452 [gr-qc] .
  * De Falco (2019) V. De Falco, arXiv e-prints , arXiv:1904.01013 (2019), arXiv:1904.01013 [gr-qc] .
  * De Falco and Bakala (2020) V. De Falco and P. Bakala, PRD 101, 024025 (2020), arXiv:1911.03649 [hep-th] .
  * De Falco (2020) V. De Falco, arXiv e-prints , arXiv:2006.01462 (2020), arXiv:2006.01462 [gr-qc] .
  * Guckenheimer and Holmes (2002) J. Guckenheimer and P. Holmes,  _Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields_, Applied Mathematical Sciences (Springer New York, 2002).
  * Tabor (1989) M. Tabor,  _Chaos and integrability in nonlinear dynamics: an introduction_, Wiley-Interscience publication (Wiley, 1989).
  * De Falco  _et al._ (2018) V. De Falco, E. Battista, and M. Falanga, Physical Review D 97, 084048 (2018), arXiv:1804.00519 [gr-qc] .
  * De Falco and Battista (2019) V. De Falco and E. Battista, EPL 127, 30006 (2019), arXiv:1907.13354 [gr-qc] .
  * De Falco and Battista (2020) V. De Falco and E. Battista, Phys. Rev. D 101, 064040 (2020), arXiv:2003.04416 [gr-qc] .
  * Strogatz (2000) S. H. Strogatz,  _Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry and Engineering_ (Westview Press, 2000).
  * Mori  _et al._ (2013) H. Mori, G. Paquette, and Y. Kuramoto,  _Dissipative Structures and Chaos_ (Springer Berlin Heidelberg, 2013).
  * Misner  _et al._ (1973) C. W. Misner, K. S. Thorne, and J. A. Wheeler,  _San Francisco: W.H. Freeman and Co., 1973_ (1973).
  * Levin and Perez-Giz (2009) J. Levin and G. Perez-Giz, PRD 79, 124013 (2009), arXiv:0811.3814 [gr-qc] .
  * Holmes and Marsden (1982a) P. J. Holmes and J. E. Marsden, Journal of Mathematical Physics 23, 669 (1982a).
  * Holmes and Marsden (1982b) P. J. Holmes and J. E. Marsden, Communications in Mathematical Physics 82, 523 (1982b).
  * Levin (1999) J. Levin, Phys. Rev. D 60, 064015 (1999), arXiv:astro-ph/9811213 [astro-ph] .