# Timescales of the chaos onset in the
general relativistic Poynting-Robertson effect

Vittorio De Falco1<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, 25121, Brescia, Italy.

###### Abstract

It has been proved that the general relativistic Poynting-Robertson effect in
the equatorial plane of Kerr metric shows a chaotic behavior for a suitable
range of parameters. As a further step, we calculate the timescale for the
onset of chaos through the Lyapunov exponents, estimating how this trend
impacts on the observational dynamics. We conclude our analyses with a
discussion on the possibility to observe this phenomenon in neutron star and
black hole astrophysical sources.

##  I Introduction

The general relativistic Poynting-Robertson (PR) effect is a phenomenon
occurring in high-energy astrophysics. The motion of accreting matter in the
vicinity of compact objects, such as black holes (BHs) and neutron stars
(NSs), is strongly affected by the gravitational field. In addition, there
could be the presence of other perturbing effects responsible to alterate the
geodesic motion of the surrounding matter (e.g., magnetic fields,
hydrodynamical/magneto-hydrodynamical processes).

The X-ray electromagnetic radiation produced by an emitting source located
close to the compact object (e.g., very hot corona around supermassive BHs in
AGNs, stellar mass BHs, and weakly magnetised NSs in X-ray binary systems
(Fabian  _et al._ , 2015; Reis and Miller, 2013; Done  _et al._ , 2007),
boundary layers around NSs or the NS surface itself (Inogamov and Sunyaev,
1999; Popham and Sunyaev, 2001), thermonuclear flashes occurring in the
outermost layers of accreting NSs (Lewin  _et al._ , 1993)) can intercept the
inflowing matter modifying its motion. Indeed, the radiation force acting on
relatively small-sized matter elements, treated as point-like particles, is
composed by the radiation pressure (being opposite to the gravitational
attraction) and the PR radiation drag force, originating when the test
particle absorbs the incoming radiation and then re-emit it isotropically in
its rest frame Poynting (1903); Robertson (1937); Bini  _et al._ (2009, 2011).
The PR effect removes very efficiently energy and angular momentum from the
affected test particle, configuring thus as a _dissipative and non-linear
dynamical system in GR_.

From a modelling perspective, there are different treatments of the general
relativistic PR effect going from Schawarzschild and Kerr to 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._ (2020a); De Falco and Wielgus
(2021). The common feature of all models is the presence of the _critical
hypersurface_ , 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).

From a theoretical point of view, the general relativistic PR effect has been
treated under a Lagrangian formalism, determining for the first time in the GR
literature the analytical form of the Rayleigh potential (De Falco  _et al._ ,
2018; De Falco and Battista, 2019; De Falco and Battista, 2020; De Falco,
2020).

There are also several high-energy astrophysical applications, like:
investigating the way in which type I X-ray bursts on the NS surface may
induce an increased mass inflow rate in the inner edge of the accretion disk
Walker (1992); Lančová  _et al._ (2017); Fragile  _et al._ (2020); studies on
the matter velocity field close to a slowly rotating NS as a result of the PR
effect for different star luminosities Abramowicz  _et al._ (1990); Lamb and
Miller (1995); Miller and Lamb (1996); Miller  _et al._ (1998); modelling the
photospheric expansion occurring during Eddington-luminosity X-ray bursts
Wielgus  _et al._ (2015); Wielgus (2016) and associated oscillations Wielgus
_et al._ (2012); Bollimpalli  _et al._ (2019); observational evidences for
changes in the inner disk properties, possibly induced by the PR effect from
high signal to noise observations of type I X-ray bursts Ballantyne and
Strohmayer (2004); Keek  _et al._ (2014a, b, 2018); diagnosing the presence of
wormholes through the detection of metric-changes occurring in strong field
regimes around BHs through the PR critical hypersurfaces De Falco  _et al._
(2020b, 2021); the PR effect can drive the dynamical evolution of unequal
supermassive BHs coalescence in galactic nuclei Chen  _et al._ (2020).

Recently, it has been shown that the general relativistic PR effect in the
equatorial plane of the Kerr metric shows chaotic dynamics for a suitable
parameter range De Falco and Borrelli (2021). In this paper, we would like to
explore the timescale of the chaos onset and check how these configurations
influence the PR dynamics within the Lyapunov exponents theory.

The paper is organized as follows: in Sec. II we briefly describe the general
relativistic PR effect model and its equations of motion; in Sec. III we
calculate the Lyapunov exponents of the general relativistic PR effect; in
Sec. IV we discuss our results proposing some observational classes of
astrophysical sources to detect the chaos in the PR dynamics and finally draw
the conclusions.

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

We consider a test particle orbiting in the equatorial plane, $\theta=\pi/2$,
around a rotating compact object, whose outside spacetime is described by the
Kerr metric. In geometrical units ($c=G=1$), its line element, ${\rm
d}s^{2}=g_{\alpha\beta}{\rm d}x^{\alpha}{\rm d}x^{\beta}$, in standard Boyer-
Lindquist coordinates is parametrized by the mass, $M$, and spin, $a$, 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=r^{2}-2r+a^{2}$, and $\rho=r^{2}+a^{2}+2a^{2}/r$.

We assume that the motion of the test particle is influenced by the X-ray
electromagnetic radiation field coming from an emitting source located in the
vicinity of the rotating compact object. Radiation is treated here as a field
superimposed on the Kerr spacetime, located in the equatorial plane, and
constituted by a coherent flux of photons, propagating outwards from the
center of the radiation source. At each time instant the test particle is hit
by a photon, moving along null geodesics in the equatorial plane on the Kerr
spacetime. The radiation field includes also the effects of non-zero angular
momentum of the photon field, $b$, that would result from the rotation of the
central radiation source and/or frame dragging due to the rotating compact
object. Therefore, the photons can be emitted either radially, $b=0$, or in
any other direction, $b\neq 0$, in the equatorial plane. It is useful to
introduce the parameter $\beta$, defined as the azimuthal angle of the photon
four-momentum measured clockwise from the $\bm{\varphi}$-axis in the local
zero angular momentum observer (ZAMO) frame, which is related to the photon
angular momentum $b$ through the formula (Bini  _et al._ , 2011)

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

where $\beta\in[0,2\pi]$ and the functions $N,N^{\varphi}$ read as

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

It is important to note that for $\sin\beta>0$ ($\sin\beta<0$) we are
considering outgoing (ingoing) photons Bini  _et al._ (2011).

We assume that the interaction between the radiation field and the test
particle occurs through coherent and isotropic scattering in the test particle
rest frame with Thomson cross section $\sigma_{\rm T}=6.7\times 10^{-25}$ cm2.
The relative luminosity of the radiation field is encoded in the parameter
$A$, which can be written as $A/M=L_{\infty}/L_{\rm Edd}$ (Bini  _et al._ ,
2009), where $L_{\infty}$ stays for the luminosity measured by an observer at
infinity and $L_{\rm Edd}\simeq 1.26\times 10^{38}\,(M/M_{\odot})$ erg s-1
represents the Eddington limiting luminosity.

The dynamical system describing the equatorial motion of the test particle
influenced by the Kerr gravitational field, the radiation pressure and the PR
effect is governed by a set of four coupled and fully general relativistic
ordinary differential equations of the first order. The first two equations
describe the test particle motion in the ZAMO frame in terms of the local
spatial velocity $\nu$ and its azimuthal angle $\alpha$ measured clockwise
with respect to the $\bm{\varphi}$-axis in the ZAMO frame. The last two
equations transform these ZAMO quantities into the radial $r$ and angular
$\varphi$ velocities. The set of equations of motion in the Kerr metric reads
as Bini  _et al._ (2009, 2011)

$\displaystyle\frac{\rm d\nu}{\rm dt}$ $\displaystyle=$
$\displaystyle-\frac{N\sin\alpha}{\gamma^{2}}[a(n)^{\hat{r}}+2\nu\cos\alpha\,\theta(n)^{\hat{r}}{}_{\hat{\varphi}}]$
$\displaystyle+A\frac{(1+bN^{\varphi})}{\gamma
Nr\sqrt{\rho}}\frac{[\cos(\alpha-\beta)-\nu][1-\nu\cos(\alpha-\beta)]}{|\sin\beta|},$
$\displaystyle\frac{\rm d\alpha}{\rm dt}$ $\displaystyle=$
$\displaystyle-\frac{N\cos\alpha}{\nu}[a(n)^{\hat{r}}+2\nu\cos\alpha\,\theta(n)^{\hat{r}}{}_{\hat{\varphi}}+\nu^{2}k_{\rm(lie)}(n)^{\hat{r}}]$
(5) $\displaystyle+A\frac{(1+bN^{\varphi})}{\gamma
Nr\sqrt{\rho}}\frac{[1-\nu\cos(\alpha-\beta)]\sin(\beta-\alpha)}{\nu|\sin\beta|},$
$\displaystyle\frac{\rm dr}{\rm dt}$ $\displaystyle=$
$\displaystyle\frac{N\nu\sin\alpha}{\sqrt{g_{rr}}},$ (6)
$\displaystyle\frac{\rm d\varphi}{\rm dt}$ $\displaystyle=$
$\displaystyle\frac{N\nu\cos\alpha}{\sqrt{g_{\varphi\varphi}}}-N^{\varphi}.$
(7)

where $\gamma=1/\sqrt{1-\nu^{2}}$ is the Lorentz factor, and
$a(n)^{\hat{r}},\theta(n)^{\hat{r}}{}_{\hat{\varphi}},k_{\rm(lie)}(n)^{\hat{r}}$
are the acceleration vector, the expansion tensor, and the relative Lie
curvature tensor, respectively, whose explicitly expressions are (Bini  _et
al._ , 2011)

$\displaystyle a(n)^{\hat{r}}$ $\displaystyle=$
$\displaystyle\frac{M[(r^{2}+a^{2})^{2}-4a^{2}Mr]}{r^{3}\rho\sqrt{\Delta}},$
(8) $\displaystyle\theta(n)^{\hat{r}}{}_{\hat{\varphi}}$ $\displaystyle=$
$\displaystyle\frac{-aM(3r^{2}+a^{2})}{r^{3}\rho},$ (9) $\displaystyle
k_{\rm(lie)}(n)^{\hat{r}}$ $\displaystyle=$
$\displaystyle\frac{-\sqrt{\Delta}(r^{3}-a^{2}M)}{r^{3}\rho}.$ (10)

###  II.1 Critical hypersurface

The general relativistic PR effect model admits as a particular solution, one
where gravitational attraction, radiation pressure, and radiation drag force
balance. This configuration defines a region, dubbed as _critical
hypersurface_ , describing a circular orbit in the equatorial plane, where the
test particle moves on it stably with constant velocity Bini  _et al._ (2011);
De Falco and Bakala (2020). Therefore, imposing $\nu=\nu_{\rm
crit},\alpha=0,r=r_{\rm crit}$ in the equations of motion, we obtain ${\rm
d}\nu/{\rm dt}={\rm d}\alpha/{\rm dt}={\rm d}r/{\rm dt}=0$, which entails Bini
_et al._ (2011)

$\displaystyle\nu_{\rm crit}=\cos\beta,$ (11) $\displaystyle
a(n)^{\hat{r}}+2\theta(n)^{\hat{r}}{}_{\hat{\varphi}}\nu_{\rm
crit}+k_{\rm(Lie)}(n)^{\rm\hat{r}}\nu_{\rm crit}^{2}$ (12)
$\displaystyle=\frac{A(1+bN^{\varphi})^{2}\rm{sgn}(\sin\beta)}{N^{2}r\sqrt{\rho}\gamma_{\rm
crit}^{3}}.$

The first equation gives the constant velocity of the test particle on the
critical hypersurface; whereas the second solved implicitly in terms of $r$
permits to determine the radius of the critical hypersurface in terms of
$M,a,A,b$.

##  III Lyapunov exponents of the general relativistic PR effect

The theory of the Lyapunov exponents is a generalization of the linear
stability theory for dynamical systems, that allows to quantify the rate of
separation of infinitesimally close trajectories. Roughly speaking, it relies
on studying the asymptotic properties of the tangent map to the dynamical
system, and this is usually achieved analysing the linearized flow. We give a
short account on the method, referring to Benettin  _et al._ (1980a, b) for
more details.

We first quickly introduce the general framework, before describing the
concrete situation we are interested in. Let $\Phi^{t}$ be a differentiable
flow on an $n$-dimensional connected Riemannian manifold $(\mathcal{M},h)$,
where $h$ is the metric on the manifold $\mathcal{M}$111In Sec. III.1.1, we
will see that the spacetime $\mathcal{M}=\mathbb{R}^{3}$ and the metric $h$ on
$\mathcal{M}$ is the usual euclidean metric., and
$d\Phi^{t}_{x}:T_{x}\mathcal{M}\to T_{x}\mathcal{M}$ the associated tangent
map. Here $T_{x}\mathcal{M}$ denotes the tangent space to $\mathcal{M}$ at
$x\in\mathcal{M}$. The Lyapunov exponents are associated with the evolution of
tangent vectors to $\mathcal{M}$. Namely, given a non-zero $\bm{v}\in
T_{x}\mathcal{M}$ they are defined as

$\chi(x,\bm{v}):=\lim_{t\to+\infty}\frac{\ln\|d\Phi^{t}_{x}\bm{v}\|_{h}}{t}\,.$
(13)

As $\bm{v}$ varies in $T_{x}\mathcal{M}$, $\chi(x,\bm{v})$ takes $s\leq n$
distinct values $\nu_{1}(x)>\ldots>\nu_{s}(x)$, with $s=s(x)$. Then one can
compute _all_ the Lyapunov exponents, letting vary vectors $\bm{v}\in
T_{x}\mathcal{M}$ in Eq. (13). However, as explained in Benettin  _et al._
(1980a), in practical computations a random choice of vectors always leads to
the _largest_ Lyapunov exponent, because the others have essentially zero
probability to be obtained.

Let $\chi_{i}(x)$, with $\chi_{1}(x)\geq\ldots\chi_{n}(x)$, be the Lyapunov
exponents at $x$ counted with their multiplicity. The sum of the first $p\geq
1$ Lyapunov exponents can be computed choosing $p$ vectors
$\bm{v_{1}},\ldots,\bm{v_{p}}\in T_{x}\mathcal{M}$ and then evaluating the
volume of the parallelepiped generated by the transformed vectors
$d\Phi^{t}\bm{v_{1}},\ldots,d\Phi^{t}\bm{v_{p}}$ through the formula:

$\lim_{t\to\infty}\frac{1}{t}\operatorname{Vol^{p}_{h}}([d\Phi^{t}_{x}\bm{v_{1}},\ldots,d\Phi^{t}_{x}\bm{v_{p}}])=\sum^{p}_{i=1}\chi_{i}(x)\,.$
(14)

As an example of the applicability of the above formula, let us suppose, for
instance, we want to compute the second Lyapunov exponent, $\chi_{2}(x)$. The
first one $\chi_{1}(x)$ is easily determined by randomly choosing a vector
$\bm{v_{1}}\in T_{x}\mathcal{M}$ and then employing Eq. (13). Therefore, now
one needs to choose just another independent vector $\bm{v_{2}}\in
T_{x}\mathcal{M}$, compute the volume of the paralleliped
$[d\Phi^{t}\bm{v_{1}},d\Phi^{t}\bm{v_{2}}]$, and finally take the limit for
$t\to\infty$:

$\lim_{t\to\infty}\frac{1}{t}\operatorname{Vol^{2}_{h}}([d\Phi^{t}_{x}v_{1},d\Phi^{t}_{x}v_{2}])=\chi_{1}(x)+\chi_{2}(x)\,.$
(15)

Subtracting from this sum the value of $\chi_{1}(x)$, one obtains
$\chi_{2}(x)$. It is now evident how to iteratively apply this method for
computing all the Lyapunov exponents $\chi_{i}(x)$. One simply needs to add
each time a new indepedendent vector $\bm{v_{i}}\in T_{x}\mathcal{M}$ and to
use Eq. (14).

However, as clarified in Benettin  _et al._ (1980b), this procedure requires a
careful implementation from a numerical point of view, as two main issues
arise. First, randomly choosing $\bm{v}\in T_{x}M$,
$\|dT^{t}_{x}\bm{v}\|_{h}$, it exponentially increases as $t\to\infty$.
Moreover, when at least two vectors are involved in (14), the angle between
two of them can rapidly become too small to be numerically handled. In
Benettin  _et al._ (1980b) the authors solve all these issues and then provide
tests and examples illustrating their computational strategy.

The concrete setting we are considering is those dynamical systems governed by
a set of ordinary differential equations of the following type

$\frac{{\rm d}X_{i}(t)}{{\rm d}t}=f_{i}(X_{1}(t),\dots,X_{n}(t)),\quad
i=1,\dots,n,$ (16)

where $f_{i}$ are supposed to be smooth functions for all $i=1,\dots,n$, and
$t$ plays the role of the time. Let us assume that a time-dependent stationary
solution $(X_{1}^{*}(t),\dots,X_{n}^{*}(t))$ exists for the dynamical system
(16), i.e., $f_{i}(X_{1}^{*}(t),\dots,X_{n}^{*}(t))=0$ for all $i=1,\dots,n$.
Let us consider now a small perturbation around the stationary solution, i.e.,
$X_{i}(t)=X_{i}^{*}(t)+\delta X_{i}(t)$ for all $i=1,\dots,n$, and linearize
thus the dynamical system (16),

$\frac{{\rm d}(\delta X_{i}(t))}{{\rm d}t}=\mathbb{A}_{ij}(t)\delta
X_{j}(t),\quad i,j=1,\dots,n,$ (17)

where $\mathbb{A}_{ij}(t)$ is the _linear stability matrix_ obtained as

$\mathbb{A}_{ij}(t)=\left.\frac{\partial f_{i}}{\partial
X_{j}}\right|_{X_{i}^{*}(t)},\quad i,j=1,\dots,n.$ (18)

The solution of the linearized equation (17) is

$\delta X_{i}(t)=L_{ij}(t)\delta X_{j}(0),\quad i,j=1,\dots,n,$ (19)

where the _evolution matrix_ $L_{ij}(t)$ satisfies

$\dot{L}_{ij}(t)=\mathbb{A}_{im}(t)L_{mj}(t),\quad i,j=1,\dots,n,$ (20)

and $L_{ij}(0)=\delta_{ij}$. We observe that, a priori, $\mathbb{A}_{im}$
depends on the time $t$. In the special case where the stationary solution is
independent from the time $t$, the Lyapunov exponents are simply the
eigenvalues of the time-independent matrix $\mathbb{A}_{ij}$. Then the
_principal Lyapunov exponent_ is the maximal eigenvalue of such a matrix.

####  III.0.1 Proprieties of the Lyapunov exponents

For _conservative_ dynamical systems the sum of all Lyapunov exponents is
zero, because a volume element in the phase space is conserved by the flow.
Instead if the dynamical system is dissipative the sum of Lyapunov exponents
is negative, as the volume element shrinks along a trajectory. Therefore, the
Lyapunov exponents provide information concerning local expansion and
contraction of phase space, formalizing thus the concept of stretching rate
along different directions Abarbanel  _et al._ (1993).

Positive Lyapunov exponents are an useful index of the _sensitive dependence
on the initial conditions_ 222A dynamical system shows sensitive dependence on
the initial conditions if tiny perturbations on the initial conditions leads
to significantly different future behaviours. Wiggins (1988); Tabor (1989);
Strogatz (2000); Ott (2002); Guckenheimer and Holmes (2002). Therefore, they
are usually taken as a possible indication of chaos provided that some other
conditions are satisfied333The widely accepted definition of _chaos_ is due to
Robert L. Devaney and it fulfils three proprieties Devaney (2018): (1)
sensitive dependence on initial conditions, (2) topologically mixing (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 (every point in the dynamical real space is approached arbitrarily
close by periodic orbits).. Lyapunov exponents are used to characterize
unstable orbits, along which chaotic dynamics can develop. As an example, the
instability of some circular orbits around a Schwarzschild BH can be
quantified by a positive principle Lyapunov exponent, although the geodesics
around a Schwarzschild BH are not chaotic Cornish (2001). If a dynamical
system shows a positive principal Lyapunov exponent, $\bar{\lambda}>0$, it is
possible to define the _Lyapunov time_ $T_{\bar{\lambda}}=1/\bar{\lambda}$.
This is a characteristic timescale on which a dynamical system is unstable or
chaotic and beyond which our predictions break down Strogatz (2000); Cornish
(2001); Cornish and Levin (2003).

Despite being helpful, Lyapunov exponents must be exploited with caution as
they present some drawbacks, especially in GR theory, which can be summarised
as:

  * •

_they vary from orbit to orbit_ , not encoding the collective behaviour of all
orbits and not usually catching generic information Strogatz (2000). In order
to get true and appropriate values, they should be averaged over many
different points on the same trajectory. Sometimes it can occur that in such a
mean process they could return zero values when the considered orbits move in
and out of unstable regions Barrow (1981, 1982); Hobill  _et al._ (1994);

  * •

_they depend on the chosen time coordinate_. Since time is relative in GR, the
same happens also for the Lyapunov exponents. If this remark is not taken
properly into account can bring to erroneous results. However, whenever a
preferred time direction exists, the uncertainty of time can be eliminated.
For example, in the Schwarzschild and Kerr spacetimes where a timelike Killing
vector exists, the coordinate time of the observer at infinity reveals to be
the most appropriate choice Cornish (2001); Cornish and Levin (2003).

###  III.1 Application of the Lyapunov theory to the general relativistic PR
effect

The general relativistic PR effect is _rotationally invariant_ (independent
from the azimuthal angle $\varphi$) and _autonomous_ (does not explicitly
depend on the time $t$). Therefore Eqs. (II), (5), (6) represent the dynamical
system to investigate. Defined $\bm{X}=(\nu,\alpha,r)$, the dynamical system
can be written as ${\rm d}\bm{X}/{\rm d}t=\bm{f}(\bm{X})$. We linearly perturb
it around the critical hypersurface $\bm{X_{0}}=(\nu_{0},0,r_{0})$ (being a
stationary solution, i.e., $\bm{f}(\bm{X_{0}})=\bm{0}$), where $\nu_{0}$ and
$r_{0}$ can be determined by exploiting Eqs. (11) and (12). We consider the
following perturbations

$\nu=\nu_{0}+\varepsilon\nu_{1},\quad\alpha=\varepsilon\alpha_{1},\quad
r=r_{0}+\varepsilon r_{1},\quad\varepsilon\ll 1,$ (21)

also written as $\bm{X}=\bm{X_{0}}+\varepsilon\bm{X_{1}}$, with
$\bm{X_{1}}=(\nu_{1},\alpha_{1},r_{1})$.

The components of the linearized $3\times 3$ matrix $\mathbb{A}=({\rm
d}\bm{f}/{\rm d}\bm{X})_{\bm{X}=\bm{X_{0}}}$ can be found in Table 1, where we
obtain the linearized dynamical system ${\rm d}\bm{X_{1}}/{\rm
d}t=\mathbb{A}\cdot\bm{X_{1}}$.

Table 1:  Explicit expressions of the coefficients of the linearized matrix
$\mathbb{A}$, where $a_{ij}=\partial f_{i}/\partial X_{j}$. The quantities
with subscript $0$ mean that they are evaluated in $\nu_{0},r_{0}$. They can
be found also in the related Mathematica notebook.
---
Coefficient | Explicit expression
$\chi_{1}$ | $2aM\rho^{3}r_{0}\left[2a^{2}r_{0}\left(-8M^{2}+Mr_{0}+6r_{0}^{2}\right)+a^{4}(4M+3r_{0})+r_{0}^{4}(9r_{0}-14M)\right]$
$\chi_{2}$ | $-3\rho_{0}^{2}\left[r_{0}^{7}\left(20a^{2}M^{2}+1\right)-Mr_{0}^{6}\left(32a^{2}M^{2}+5\right)+2r_{0}^{5}\left(12a^{4}M^{2}+a^{2}+3M^{2}\right)\right.$
|
$\left.+a^{2}r_{0}^{3}\left(-16a^{2}M^{4}+4a^{2}M^{2}+a^{2}+12M^{2}\right)-a^{2}Mr_{0}^{2}\left(5a^{2}+12M^{2}\right)+10a^{4}M^{2}r_{0}-8a^{2}Mr_{0}^{4}-2a^{6}M\right]$
$\chi_{3}$ | $2aM\rho_{0}r_{0}^{2}\left[-4\text{Ma}^{2}r_{0}^{3}\left(2a^{2}M^{2}+9\right)+r_{0}^{6}\left(44a^{2}M^{2}+3\right)-18Mr_{0}^{5}\left(4a^{2}M^{2}+1\right)+3r_{0}^{4}\left(16a^{4}M^{2}+3a^{2}+8M^{2}\right)\right.$
|
$\left.+a^{2}r_{0}^{2}\left(4a^{2}M^{2}+9a^{2}+36M^{2}\right)-2a^{2}Mr_{0}\left(8a^{4}M^{2}+9a^{2}+12M^{2}\right)+3a^{4}\left(a^{2}+4M^{2}\right)\right]$
$\chi_{4}$ | $r_{0}\left[6a^{4}r_{0}^{2}\left(-2a^{4}M^{2}+12a^{2}M^{4}+a^{2}+2M^{2}\right)+4a^{2}r_{0}^{4}\left(-8a^{4}M^{6}-2a^{4}M^{2}+3a^{2}+12M^{2}\right)\right.$
|
$\left.+2a^{2}Mr_{0}^{3}\left(16a^{6}M^{3}+6a^{4}M^{2}-24a^{2}M^{3}-15a^{2}-4M^{2}\right)+3r_{0}^{8}\left(1-16a^{4}M^{4}\right)\right.$
|
$\left.+2Mr_{0}^{7}\left(40a^{4}M^{4}+6a^{2}M^{2}-9\right)+2r_{0}^{6}\left(-24a^{6}M^{4}-6a^{4}M^{2}-12a^{2}M^{4}+5a^{2}+18M^{2}\right)\right.$
|
$\left.+2Mr_{0}^{5}\left(8a^{6}M^{4}+24a^{4}M^{2}-21a^{2}-12M^{2}\right)-6a^{6}Mr_{0}\left(4a^{2}M^{2}+1\right)+a^{8}\right]$
$\xi_{1}$ | $-6aMr_{0}^{2}\rho_{0}^{2}\left(a^{2}+3r_{0}^{2}\right)$
$\xi_{2}$ | $\rho_{0}\left[2r_{0}^{5}\left(18a^{2}M^{2}-1\right)+3a^{2}r_{0}^{3}\left(4a^{2}M^{2}-1\right)+a^{2}r_{0}\left(6M^{2}-a^{2}\right)-a^{2}Mr_{0}^{2}-2a^{4}M+3Mr_{0}^{4}\right]$
$\xi_{3}$ | $2aMr_{0}\left[3r_{0}^{5}\left(1-4a^{2}M^{2}\right)+a^{2}r_{0}^{3}\left(3-4a^{2}M^{2}\right)+2a^{2}M^{2}r_{0}-a^{2}Mr_{0}^{2}-2a^{4}M-5Mr_{0}^{4}\right]$
$\psi_{1}$ | $\nu_{0}^{2}\left[a^{2}\left(2M^{2}-Mr_{0}+r_{0}^{2}\right)+r_{0}^{3}(r_{0}-2M)\right]+M\left(2a^{2}(r_{0}-2M)+r_{0}^{3}\right)$
$\psi_{2}$ | $2a^{3}bMr_{0}^{3}\left(-2M^{2}+Mr_{0}-3r_{0}^{2}\right)+4a^{5}bM^{2}r_{0}^{2}+a^{2}r_{0}^{4}\left(-12M^{2}+2Mr_{0}+5r_{0}^{2}\right)$
|
$+a^{4}r_{0}(2M+r_{0})\left(-6M^{2}+Mr_{0}+4r_{0}^{2}\right)+a^{6}(2M+r_{0})^{2}+2abMr_{0}^{6}(5M-3r_{0})+r_{0}^{7}(2r_{0}-3M)$
$\psi_{3}$ | $-2aM\nu_{0}\left[a^{2}r_{0}^{3}\left(12M^{2}-23Mr_{0}+17r_{0}^{2}\right)+3a^{6}(M+r_{0})-a^{4}r_{0}(M-r_{0})(4M+11r_{0})+3r_{0}^{6}(3r_{0}-7M)\right]$
$\psi_{4}$ | $-M\left[2a^{4}r_{0}\left(2M^{2}+5Mr_{0}-4r_{0}^{2}\right)+a^{2}r_{0}^{4}(25M-7r_{0})-3a^{6}(M+r_{0})-2r_{0}^{7}\right]$
$\psi_{5}$ | $-a^{2}r_{0}^{4}\left(20M^{2}-9Mr_{0}+r_{0}^{2}\right)+a^{4}Mr_{0}\left(-2M^{2}-5Mr_{0}+10r_{0}^{2}\right)+3a^{6}M(M+r_{0})+r_{0}^{7}(4M-r_{0})$
$a_{11}$ | $\frac{A\gamma_{0}}{r_{0}}\frac{\left(1-2\nu_{0}^{2}\right)\left[(\rho_{0}-2abMr_{0})^{2}+b^{2}\Delta_{0}\right]-b\sqrt{\Delta_{0}}\nu_{0}\left(3\nu_{0}^{2}-1\right)(2abMr_{0}-\rho_{0})}{\sqrt{\Delta_{0}}\rho_{0}(2abMr_{0}-\rho_{0})\sqrt{1-\frac{b^{2}\Delta_{0}}{(\rho_{0}-2abMr_{0})^{2}}}}$
$a_{12}$ | $\frac{\frac{A\gamma_{0}\sqrt{\rho_{0}}r_{0}^{2}\left[\left(\nu_{0}^{2}+1\right)(2abMr_{0}+\rho_{0})-2b\sqrt{\Delta_{0}}\nu_{0}\right]}{\sqrt{\Delta_{0}}}-M\left[-4a^{2}Mr_{0}-2a\sqrt{\Delta_{0}}\nu_{0}\left(a^{2}+3r_{0}^{2}\right)+(a^{2}+r_{0}^{2})^{2}\right]}{\gamma_{0}^{2}\rho_{0}^{3/2}r_{0}^{3}}$
$a_{13}$ | $-A\frac{\left\\{b\Delta_{0}^{3/2}\left(\nu_{0}^{2}+1\right)\left[b^{3}\xi_{3}+b^{2}\xi_{2}+b\xi_{1}+\rho_{0}^{3}r_{0}\left(a^{2}+3r_{0}^{2}\right)\right]+\nu_{0}\left(b^{4}\chi_{4}+b^{3}\chi_{3}+b^{2}\chi_{2}+b\chi_{1}-r_{o}\rho^{6}\right)\right\\}}{\gamma_{0}\Delta_{0}^{3/2}\rho_{0}^{2}r_{0}^{3}\left(1-\frac{b^{2}\Delta_{0}}{(\rho_{0}-2abMr_{0})^{2}}\right)^{3/2}}$
$a_{21}$ | $\frac{\sqrt{\Delta_{0}}\left(a^{4}M+\gamma_{0}^{2}r_{0}\psi_{1}\right)-A\gamma_{0}\sqrt{\rho_{0}}\left(-2abMr_{0}-b\sqrt{\Delta_{0}}\nu_{0}^{3}+\rho\right)}{\sqrt{\Delta_{0}}\rho^{3/2}r_{0}^{2}\nu_{0}^{2}}$
$a_{22}$ | $A\frac{\nu_{0}\left[(\rho-2abMr_{0})^{2}-2b^{2}\Delta_{0}\right]+b\sqrt{\Delta_{0}}(\rho_{0}-2abMr_{0})}{\gamma_{0}\sqrt{\Delta_{0}}\rho_{0}r_{0}\nu_{0}(2abMr_{0}-\rho_{0})\sqrt{\frac{b^{2}\Delta_{0}}{(\rho-2abMr_{0})^{2}}+1}}$
$a_{23}$ | $-\frac{A\sqrt{\rho_{0}}r_{0}\left[\psi_{2}-b\Delta_{0}^{3/2}r_{0}^{2}\nu_{0}\left(a^{2}+3r_{0}^{2}\right)\right]-\gamma_{0}\Delta_{0}M\left[\sqrt{\Delta_{0}}\left(\nu_{0}^{2}\psi_{5}+\psi_{4}\right)+\psi_{3}\right]}{\gamma_{0}\Delta_{0}^{3/2}\rho_{0}^{5/2}r_{0}^{5}\nu_{0}}$
$a_{32}$ | $\frac{\Delta_{0}\nu_{0}}{\sqrt{\rho_{0}}r_{0}}$
$a_{31}=a_{33}$ | 0

As already explained, since the dynamical system is autonomous the Lyapunov
exponents coincide with the eigenvalues of the matrix $\mathbb{A}$, denoted as
$\left\\{\lambda_{1},\lambda_{2},\lambda_{3}\right\\}$. The _characteristic
eigenvalue equation_ in terms of $\lambda$ is

$c_{0}+c_{1}\ \lambda+c_{2}\ \lambda^{2}-\lambda^{3}=0,$ (22)

where

$\displaystyle c_{0}={\rm det}\mathbb{A},\qquad c_{2}={\rm Tr}\mathbb{A},$
(23) $\displaystyle c_{1}=\sum_{\begin{subarray}{c}i,j=1\\\ i\neq
j\end{subarray}}^{3}\left(\frac{a_{ij}a_{ji}-a_{ii}a_{jj}}{2}\right).$

Since the eigenvalues are the zeroes of a polynomial of third order, we can
analytically determine them Press  _et al._ (2002).

####  III.1.1 Lyapunov and PR timescales

We can calculate the Lyapunov timescale for the chaos onset as
$T_{\bar{\lambda}}=1/\bar{\lambda}$, where $\bar{\lambda}$ is the principal
Lyapunov exponent. This time must be compared with the PR timescale $T_{\rm
PR}$, defined as the time from the start of the numerical simulation until the
test particle reaches for the first time the critical hypersurface for moving
then on that stably444There are some examples, where the test particle can
cross the critical hypersurface and then not moving anymore on that (see
figures in Sec. 3.4 of Ref. De Falco (2019), for more details).. This time is
calculated numerically for each simulation. In order to understand whether the
chaotic behavior is influential on the global PR dynamics, we compute the
ratio $T_{\bar{\lambda}}/T_{\rm PR}$, and if it is smaller (greater) than one,
then the chaos is important (unimportant) Cornish and Levin (2003).

As examples we consider selected data from our previous numerical simulations,
see Fig. 1 (see also Fig. 7 in Ref. De Falco and Borrelli (2021), for
comparison). The specific data for these simulations are reported in Table 2,
see simulations 1 – 15. We added also other five examples not reported in
previous studies in order to show some variability of parameters, see
simulations 16 – 20 in Table 2.

Figure 1: Parameter space $(r_{u},a,b)$, where $r_{u}$ is the homoclinic
orbit’s parameter, and $b$ is fixed to $b=3$. This space is delimited by the
innermost bound circular orbit, $r_{\rm IBCO}$, and the innermost stable
circular orbit, $r_{\rm ISCO}$. All the (blue and red) dots represent the
parameter subspace in correspondence of which the chaotic dynamic occurs. The
red dots are the selected examples reported in simulations 1 – 15 in Table 2.
Table 2:  Numerical simulations performed using as test particle’s initial
conditions $(r_{p},\alpha_{p},\nu_{p})$ and the PR effect model parameters
$b,A,a$ together with the critical hypersurface radius $r_{0}$. As a result we
obtain the Lyapunov timescale $T_{\lambda}$, the PR timescale $T_{\rm PR}$,
and the the ratio between these two times in order to see whether the PR
effect is important.
---
SIM. | $\quad\bm{\nu_{p}}\quad$ | $\quad\bm{\alpha_{p}}\quad$ | $\quad\bm{r_{p}}\quad$ | $\quad\bm{b}\quad$ | $\quad\bm{A}\quad$ | $\quad\bm{a}\quad$ | $\quad\bm{r_{0}}\quad$ | $\quad\bm{T_{\bar{\lambda}}}\quad$ | $\quad\bm{T_{\rm PR}}\quad$ | $\quad\bm{T_{\bar{\lambda}}/T_{\rm PR}}\quad$
$\bm{\\#}$ |  |  | $\bm{(M)}$ | $\bm{(M)}$ |  |  | $\bm{(M)}$ | $\bm{(M)}$ | $\bm{(M)}$ |
1 | 0.82 | 0.52 | 2.84 | 3 | 0.0001 | 0 | 2.00 | $7.97\times 10^{6}$ | 178703.57 | 44.59
2 | 0.72 | 0.30 | 3.51 | 3 | 0.0001 | 0 | 2.00 | $7.97\times 10^{6}$ | 178703.57 | 43.60
3 | 0.89 | 0.76 | 2.39 | 3 | 0.001 | 0.1 | 1.9950 | 587441.09 | 47363.77 | 12.40
4 | 0.92 | 0.91 | 2.25 | 3 | 0.001 | 0.1 | 1.9950 | 587441.09 | 47397.65 | 12.39
5 | 0.79 | 0.45 | 2.80 | 3 | 0.01 | 0.2 | 1.9800 | 14352.60 | 5164.82 | 2.78
6 | 0.75 | 0.35 | 3.03 | 3 | 0.01 | 0.2 | 1.9800 | 14352.60 | 5187.24 | 2.77
7 | 0.75 | 0.34 | 2.92 | 3 | 0.01 | 0.3 | 1.9541 | 2470.75 | 5347.96 | 0.46
8 | 0.73 | 0.29 | 3.06 | 3 | 0.01 | 0.3 | 1.9541 | 2470.75 | 5350.97 | 0.46
9 | 0.77 | 0.36 | 2.70 | 3 | 0.01 | 0.4 | 1.9167 | 15.79 | 5558.25 | $2.84\times 10^{-3}$
10 | 0.73 | 0.29 | 2.88 | 3 | 0.01 | 0.4 | 1.9167 | 15.79 | 5583.53 | $2.83\times 10^{-3}$
11 | 0.29 | 0.73 | 14.31 | 3 | 0.1 | 0.5 | 4.55 | 24.14 | 3397.87 | $7.10\times 10^{3}$
12 | 0.40 | 0.33 | 7.35 | 3 | 0.1 | 0.5 | 4.55 | 24.14 | 2691.15 | 0.01
13 | 0.62 | 0.04 | 3.37 | 3 | 0.1 | 0.6 | 5.62 | 213.32 | 3733.26 | 0.06
14 | 0.62 | 0.03 | 3.42 | 3 | 0.1 | 0.6 | 5.62 | 213.32 | 3727.37 | 0.06
15 | 0.14 | 0.24 | 21.31 | 3 | 0.1 | 0.7 | 6.16 | 264.20 | 264.20 | 0.06
16 | 0.70 | 0.22 | 3.29 | 3.1 | 0.1 | 0.3 | 1.9698 | 4.83 | 1029.39 | $4.69\times 10^{-3}$
17 | 0.67 | 0.11 | 3.17 | 3.1 | 0.1 | 0.5 | 5.88 | 619.76 | 4193.14 | 0.15
18 | 0.63 | 0.03 | 3.61 | 3.2 | 0.1 | 0.5 | 6.82 | 557.15 | 5227.02 | 0.11
19 | 0.61 | 0.05 | 4.03 | 3.3 | 0.1 | 0.4 | 7.22 | 415.77 | 5608.68 | 0.07
20 | 0.60 | 0.05 | 4.36 | 3.4 | 0.1 | 0.3 | 7.66 | 499.52 | 6182.17 | 0.08

##  IV Discussions and Conclusions

The strategy to detect chaos in the PR dynamics relies on computing the
_Melnikov integral_ $\mathscr{M}$, whose integrating function is the Poisson
brackets of the unperturbed Kerr Hamiltonian function and the PR dissipative
perturbations. This is a function of mass $M$, spin $a$, homoclinic orbit
parameter $r_{u}$, photon impact parameter $b$, and test particle’s initial
condition $r_{0}$, i.e., $\mathscr{M}=\mathscr{M}(M,a,r_{u},b;r_{0},)$ De
Falco and Borrelli (2021). If there is a combination of these parameters which
nullifies $\mathscr{M}$, this means that in correspondence of these values
chaotic dynamics occurs. The general relativistic PR effect shows chaos for a
suitable range of parameters, see Fig. 1 for an example555For discovering
other ranges of parameters for which chaos occurs, we developed a code written
in Mathematica, which permits to facilitate this research (see Ref. [47] in De
Falco and Borrelli (2021), for more details)..

As subsequent analysis, in this work we have investigated whether the chaotic
dynamics is observationally relevant on the PR dynamics. To this end, we have
exploited the theory of Lyapunov exponents, which allows to estimate the
timescale of the chaos onset. The procedure is essentially based on perturbing
the general relativistic PR equations of motion (II) – (6) for low
luminosities in terms of the parameter $\varepsilon=A/M\ll 1$ around the
critical hypersurface values (configuration of equilibrium for the PR
dynamical system), see Eq. (21). Therefore, after performing these
calculations, we obtain the linearized matrix $\mathbb{A}$, which is
numerically determined once $M,a,A,b,r_{0}$ have been assigned. Since the PR
model is an autonomous dynamical system, the eigenvalues of $\mathbb{A}$ are
exactly the _Lyapunov exponents_
$\left\\{\lambda_{1},\lambda_{2},\lambda_{3}\right\\}$. Finally, considering
the maximum of the real parts of the eigenvalues, we determine the _principal
Lyapunov exponent_ $\bar{\lambda}$, whose inverse value corresponds to the
_Lyapunov timescale_ $T_{\bar{\lambda}}=1/\bar{\lambda}$ for estimating the
chaos onset. Another fundamental information is encoded in the PR timescale
$T_{\rm PR}$, defined as the time from the start of the numerical simulation
until the test particle reaches for the first time the critical hypersurface
for then moving on it stably. In order to understand how the chaotic behaviour
impacts on the PR dynamics, we consider the ratio between the Lyapunov
$T_{\bar{\lambda}}$ and PR $T_{\rm PR}$ timescales: if
$T_{\bar{\lambda}}/T_{\rm PR}\leq 1$ it means that the chaos is
observationally significant, whereas if $T_{\bar{\lambda}}/T_{\rm PR}>1$ is
unimportant.

We performed 20 numerical simulations, whose detailed values are reported in
Table 2. We note that the chaos is significant for simulations 7 – 20, while
there is any influence for simulations 1 – 6. However, we checked by
performing other numerical simulations that varying the luminosity parameter
$A$, namely making it smaller and smaller, it is possible to have also for
simulations 1 – 6 a ratio lower than one. Once we calculate the Lyapunov
timescale, we have an indication on how to tune the parameters for having a PR
timescale such that it is smaller or greater than the time for the chaos
onset. We note also that fixed the value of the photon impact parameter $b$,
and chosen, for hypothesis, a small value of the luminosity parameter $A$, it
follows that an important role is played by the spin $a$. Indeed, the test
particle’s initial conditions are also fundamental, but they are strongly
related to $\left\\{M,a,b\right\\}$ via the Melnikov integral $\mathscr{M}$ De
Falco and Borrelli (2021).

From an observational point of view, it is significant to understand how to
identify the astrophysical systems where chaos in the PR dynamics can be
detected. The requirement of low luminosities permits to provide a first
stringent criterion. The mass $M$ and spin $a$ of a compact object can be
normally estimated by means of several strategies (see e.g., Falanga  _et al._
(2015); Middleton (2016), for more details). In addition, the surrounding
accreting matter can be found distributed almost everywhere with generally
different velocities, including some initial configurations for having chaotic
dynamics Frank  _et al._ (2002). Therefore, only the photon impact parameter
$b$ must be estimated. It cannot be measured directly from the observational
data, but it can be linked to the emitting surface radius $R_{\star}$
(supposed to be a spherical region) and angular velocity $\Omega_{\star}$
(assuming that the emitting surface rigidly rotates) through the formula (see
Ref. Bakala  _et al._ (2019), for more details)

$b\equiv\left[-\frac{g_{t\varphi}+g_{\varphi\varphi}\Omega_{\star}}{g_{tt}+g_{t\varphi}\Omega_{\star}}\right]_{r=R_{\star}}=\frac{a^{2}R_{\star}\Omega_{\star}\rho(R_{\star})}{R_{\star}+2M(a\Omega_{\star}-1)}.$
(24)

In addition, we have that $\Omega_{\star}\in[\Omega_{-},\Omega_{+}]$, being

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

Since $R_{\star}$ can be estimated from the observations, we can relate
$\Omega_{\star}$ in terms of $b,a,M,R_{\star}$, via Eq. (24), as

$\Omega_{\star}=\frac{2aM+b(R_{\star}-2M)}{\rho(R_{\star})-2Mab}.$ (26)

The last condition imposes thus a further constraint to single out the
astrophysical systems exhibiting chaos.

It is important now to distinguish the physics of BH and NS systems. For a
_standard NS_ of mass $M=1.4M_{\odot}$ and radius $R_{\star}=6M$, if we
consider the NS surface as emitting region, then the spin $a$ can be expressed
as a function of the NS angular velocity $\Omega_{\star}$ through
$a=\mathcal{C}\Omega_{\star}/z$, where $\mathcal{C}$ depends on the NS
structure and equation of state, which in our case amounts to be
$\mathcal{C}=1.1\times 10^{-4}$ s/rad Bakala  _et al._ (2012), and
$z=(1.4GM_{\odot}/c^{3})/(2\pi)$ is the conversion gravitational factor. In
this way, we have $a=a(b)$, which for $b\sim 3$ we obtain $a\sim 10^{-6}$,
corresponding thus to extremely slowly rotating NSs, namely
$\Omega_{\star}\sim 0.06\ {\rm rad/s}$ or spin period $T_{\star}\equiv
2\pi/\Omega_{\star}\sim 100$ s. In Fig. 2, we plot how the NS angular velocity
$\Omega_{\star}$ and spin period $T_{\star}$ changes in terms of the photon
impact parameter $b$. In addition, since NSs are very small compact objects,
they have a very low luminosity. Assuming that the NS surface temperature is
$T=10^{6}$ K, its luminosity can be calculated through the Stefan-Boltzmann
law in terms of the Sun’s luminosity (where Sun’s temperature $T_{\odot}=5800$
K and radius $R_{\odot}=7\times 10^{5}$ Km), obtaining thus $L/L_{\odot}=0.2$.
The Eddington luminosity is $L_{\rm Edd}=4.48\times 10^{4}L_{\odot}$,
therefore $A\leq 4.46\times 10^{-}6$. In rows 1, 2, 3 in Table 3, we report
the data of some NSs, in which our study can be applied.

Figure 2: Plot of the NS surface’s angular velocity $\Omega_{\star}$ and spin
period $T_{\star}$ in terms of photon impact parameter $b$.

Instead, if we consider either a boundary layer around a NS or a hot corona
around a BH, $a$ and $\Omega_{\star}$ are now independent. In the case of a
_boundary layer around a NS_ we know that it is located very close to the NS
surface, $R_{\star}\sim(5-7)M$, and it is rigidly rotating with Keplerian
angular velocity $\Omega_{\star}=\Omega_{K}(R_{\star})\equiv
M/(a+R_{\star})^{3/2}$ Popham and Sunyaev (2001). The additional requirement
of low luminosities sees these configurations hosted in atoll sources,
characterised by $A\lesssim 0.06$ Ludlam  _et al._ (2019). In Fig. 3 we plot
the photon impact parameter $b$ in terms of the NS spin $a$, the only free
parameter in this case. Finally in rows 3, 4, 5 in Table 3, there are the data
of some astrophysical examples.

Figure 3: Plot of the NS boundary layer’s photon impact parameter $b$ in terms
of the NS spin $a$. The continuous ($R_{\star}=5M$) and dashed
($R_{\star}=7M$) lines delimit the light grey shaded area, which includes all
admissible physical cases.

In the BH case, the emitting surface is represented by the hot corona, which
is located in the range of $R_{\star}\sim(3-10)M$. We assume that the matter
is rotating with Keplerian velocity $\Omega_{\star}=\Omega_{K}(R_{\star})$.
There are however different models of the hot corona’s angular velocity
proposed in the literature (see Ref. Niedźwiecki (2005), for more details). We
impose also that these sources must be characterised by very small
luminosities $A\lesssim 0.01$. In Fig. 4 we plot the photon impact parameter
$b$ in terms of the BH spin values using different models of the hot corona’s
angular velocity as reported in Ref. Niedźwiecki (2005), in order to show how
the $b$ range values change in terms of different approaches. In rows 7, 8 ,9
in Table 3 we show the data of some astrophysical sources.

Figure 4: Plot of the BH hot corona’s photon impact parameter $b$ in terms of
the BH spin $a$ for different locations (light grey shaded area), delimited by
$R_{\star}=3M$ (continuous line) and $R_{\star}=10M$ (dashed line) in terms of
different models of angular velocity $\Omega_{\star}$: Keplerian
$\Omega_{K}(R_{\star})$ coincident with the slab-like in the equatorial plane
(black lines), and ZAMO $\Omega_{\rm
ZAMO}=[-g_{t\varphi}/g_{\varphi\varphi}]_{r=R_{\star}}$ (blue line). Looking
at Eq. (24), it is understandable why to $\Omega_{\rm ZAMO}$ corresponds
always $b=0$ for all possible values of $a$ and $R_{\star}$.

The above discussions together with the related formulas, plots, and examples
provide some basic strategies, which could be more extensively investigated
and improved from an observational point of view through the analysis of
several data related to astrophysical BH and NS sources. In addition, our
approach is not only restricted to the PR effect, but it can be further
extended to study the timescales of other relevant phenomena occurring in
high-energy accretion physics (see e.g., Gaspari  _et al._ (2013); King and
Nixon (2015); Suková  _et al._ (2016); Blinova  _et al._ (2016)).

Table 3:  Information about six examples (ex.) of some astrophysical sources
divided in three types analysed in Sec. IV.
---
Type | EX. | SOURCE | $\quad\bm{M}\quad$ | $\quad\bm{R_{\star}}\quad$ | $\quad\bm{T_{\star}}\quad$ | $\quad\bm{A}\quad$ | Ref.
| $\bm{\\#}$ | id. | $\bm{{\rm(M_{\odot})}}$ | $\bm{{\rm(M)}}$ | $\bm{{\rm(s)}}$ |  |
NS surface | 1 | A0535+26 | 1.50 | 4.52 | 103.00 | $4.46\times 10^{-6}$ | Ikhsanov (2001)
2 | GX 1+4 | 1.35 | 5.03 | 121.00 | $4.46\times 10^{-6}$ | González-Galán  _et al._ (2012)
3 | Vela X-1 | 1.88 | 3.61 | 283.00 | $4.46\times 10^{-6}$ | Quaintrell  _et al._ (2003)
Boundary layer around a NS | 4 | GX 3+1 | 1.40 | 6.67 | $0.74\times 10^{-3}$ | 0.058 | Ludlam  _et al._ (2019)
5 | 4U 1702–429 | 1.40 | 5.35 | $0.53\times 10^{-3}$ | 0.058 | Ludlam  _et al._ (2019)
6 | GX 301–2 | 1.40 | 5.36 | $0.54\times 10^{-3}$ | 0.058 | Ludlam  _et al._ (2019)
Hot corona around a BH666The sources we have chosen have all extreme spin values, namely $0.95<a<1$. | 7 | NGC 5506 | $10^{8}$ | 10.00 | $1.01\times 10^{-3}$ | 0.002 | Fabian  _et al._ (2015)
8 | MCG-6-30-15 | $0.18\times 10^{-3}$ | 2.90 | $0.07$ | 0.03 | Fabian  _et al._ (2015)
9 | Cyg A | $2.51\times 10^{9}$ | 10.00 | $1.01\times 10^{-3}$ | 0.003 | Fabian  _et al._ (2015)

## 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.

## References

  * Fabian  _et al._ (2015) A. C. Fabian, A. Lohfink, E. Kara, M. L. Parker, R. Vasudevan, and C. S. Reynolds, Monthly Notices of the Royal Astronomical Society 451, 4375 (2015), arXiv:1505.07603 [astro-ph.HE] .
  * Reis and Miller (2013) R. C. Reis and J. M. Miller, ApJL 769, L7 (2013), arXiv:1304.4947 [astro-ph.HE] .
  * Done  _et al._ (2007) C. Done, M. Gierliński, and A. Kubota, A&A Rev. 15, 1 (2007), arXiv:0708.0148 .
  * Inogamov and Sunyaev (1999) N. A. Inogamov and R. A. Sunyaev, Astronomy Letters 25, 269 (1999), astro-ph/9904333 .
  * Popham and Sunyaev (2001) R. Popham and R. Sunyaev, The Astrophysical Journal 547, 355 (2001), astro-ph/0004017 .
  * Lewin  _et al._ (1993) W. H. G. Lewin, J. van Paradijs, and R. E. Taam, Space Science Reviews 62, 223 (1993).
  * 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._ (2020a) V. De Falco, P. Bakala, and M. Falanga, Phys. Rev. D 101, 124031 (2020a), arXiv:2006.01452 [gr-qc] .
  * De Falco and Wielgus (2021) V. De Falco and M. Wielgus, arXiv e-prints , arXiv:2103.17165 (2021), arXiv:2103.17165 [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] .
  * 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] .
  * Walker (1992) M. A. Walker, The Astrophysical Journal 385, 642 (1992).
  * Lančová  _et al._ (2017) D. Lančová, P. Bakala, K. Goluchová, M. Falanga, V. De Falco, and L. Stella, in  _Proceedings of RAGtime 17-19: Workshops on black holes and neutron stars, 17-19/23-26 Oct., 1-5 Nov. 2015/2016/2017, Opava, Czech Republic, Z. Stuchlík, G. Török and V. Karas editors, Silesian University in Opava, 2017, ISBN 978-80-7510-257-7, ISSN 2336-5676, p. 127-136_ (2017) pp. 127–136.
  * Fragile  _et al._ (2020) P. C. Fragile, D. R. Ballantyne, and A. Blankenship, Nature Astronomy , 6 (2020).
  * Abramowicz  _et al._ (1990) M. A. Abramowicz, G. F. R. Ellis, and A. Lanza, The Astrophysical Journal 361, 470 (1990).
  * Lamb and Miller (1995) F. K. Lamb and M. C. Miller, The Astrophysical Journal 439, 828 (1995).
  * Miller and Lamb (1996) M. C. Miller and F. K. Lamb, The Astrophysical Journal 470, 1033 (1996).
  * Miller  _et al._ (1998) M. C. Miller, F. K. Lamb, and D. Psaltis, The Astrophysical Journal 508, 791 (1998), astro-ph/9609157 .
  * Wielgus  _et al._ (2015) M. Wielgus, W. Kluźniak, A. Sadowski, R. Narayan, and M. Abramowicz, MNRAS 454, 3766 (2015), arXiv:1505.06099 [astro-ph.HE] .
  * Wielgus (2016) M. Wielgus, in  _Star Clusters and Black Holes in Galaxies across Cosmic Time_, IAU Symposium, Vol. 312, edited by Y. Meiron, S. Li, F.-K. Liu, and R. Spurzem (2016) pp. 131–134.
  * Wielgus  _et al._ (2012) M. Wielgus, A. Stahl, M. Abramowicz, and W. Kluźniak, Astronomy & Astrophysics 545, A123 (2012), arXiv:1208.2939 [astro-ph.HE] .
  * Bollimpalli  _et al._ (2019) D. A. Bollimpalli, M. Wielgus, D. Abarca, and W. Kluźniak, MNRAS 487, 5129 (2019), arXiv:1812.01299 [astro-ph.HE] .
  * Ballantyne and Strohmayer (2004) D. R. Ballantyne and T. E. Strohmayer, The Astrophysical Journal Letters 602, L105 (2004), astro-ph/0401154 .
  * Keek  _et al._ (2014a) L. Keek, D. R. Ballantyne, E. Kuulkers, and T. E. Strohmayer, ApJj 789, 121 (2014a), arXiv:1403.4940 [astro-ph.HE] .
  * Keek  _et al._ (2014b) L. Keek, D. R. Ballantyne, E. Kuulkers, and T. E. Strohmayer, ApJL 797, L23 (2014b), arXiv:1410.5434 [astro-ph.HE] .
  * Keek  _et al._ (2018) L. Keek, Z. Arzoumanian, P. Bult, E. M. Cackett, D. Chakrabarty, J. Chenevez, A. C. Fabian, K. C. Gendreau, S. Guillot, T. Güver, J. Homan, G. K. Jaisawal, F. K. Lamb, R. M. Ludlam, S. Mahmoodifar, C. B. Markwardt, J. M. Miller, G. Prigozhin, Y. Soong, T. E. Strohmayer, and M. T. Wolff, ApJL 855, L4 (2018), arXiv:1808.06441 [astro-ph.HE] .
  * De Falco  _et al._ (2020b) V. De Falco, E. Battista, S. Capozziello, and M. De Laurentis, Phys. Rev. D 101, 104037 (2020b), arXiv:2004.14849 [gr-qc] .
  * De Falco  _et al._ (2021) V. De Falco, E. Battista, S. Capozziello, and M. De Laurentis, Phys. Rev. D 103, 044007 (2021).
  * Chen  _et al._ (2020) X. Chen, D. N. C. Lin, and X. Zhang, ApjL 893, L15 (2020), arXiv:2003.11770 [astro-ph.HE] .
  * De Falco and Borrelli (2021) V. De Falco and W. Borrelli, Phys. Rev. D 103, 064014 (2021).
  * Benettin  _et al._ (1980a) G. Benettin, L. Galgani, A. Giorgilli, and J.-M. Strelcyn, Meccanica 15, 044007, pag. 9 – 20 (1980a).
  * Benettin  _et al._ (1980b) G. Benettin, L. Galgani, A. Giorgilli, and J.-M. Strelcyn, Meccanica 15, 044007, pag. 21 – 30 (1980b).
  * Abarbanel  _et al._ (1993) H. D. Abarbanel, M. I. Rabinovich, and M. M. Sushchik,  _Introduction to nonlinear dynamics for physicists_ , Vol. 53 (World Scientific, 1993).
  * Wiggins (1988) S. Wiggins,  _Global Bifurcations and Chaos: Analytical Methods_, Applied Mathematical Sciences Series No. v. 73 (Springer-Verlag, 1988).
  * Tabor (1989) M. Tabor,  _Chaos and integrability in nonlinear dynamics: an introduction_, Wiley-Interscience publication (Wiley, 1989).
  * Strogatz (2000) S. H. Strogatz,  _Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry and Engineering_ (Westview Press, 2000).
  * Ott (2002) E. Ott,  _Chaos in Dynamical Systems_ (Cambridge University Press, 2002).
  * 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).
  * Devaney (2018) R. Devaney,  _An introduction to chaotic dynamical systems_ (CRC Press, 2018).
  * Cornish (2001) N. J. Cornish, Phys. Rev. D 64, 084011 (2001), arXiv:gr-qc/0106062 [gr-qc] .
  * Cornish and Levin (2003) N. J. Cornish and J. Levin, Classical and Quantum Gravity 20, 1649 (2003), arXiv:gr-qc/0304056 [gr-qc] .
  * Barrow (1981) J. D. Barrow, Phys. Rev. Lett. 46, 963 (1981).
  * Barrow (1982) J. D. Barrow, Physics Reports 85, 1 (1982).
  * Hobill  _et al._ (1994) D. Hobill, A. Burd, and A. Coley,  _Deterministic Chaos in General Relativity_, Nato Science Series B: (Springer US, 1994).
  * Press  _et al._ (2002) W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery,  _Numerical recipes in C++ : the art of scientific computing_ (2002).
  * De Falco (2019) V. De Falco, arXiv e-prints , arXiv:1904.01013 (2019), arXiv:1904.01013 [gr-qc] .
  * Falanga  _et al._ (2015) M. Falanga, T. Belloni, P. Casella, M. Gilfanov, P. Jonker, and A. King, eds.,  _The Physics of Accretion onto Black Holes_ (Springer-Verlag New York, 2015).
  * Middleton (2016) M. Middleton, “Black Hole Spin: Theory and Observation,” in  _Astrophysics of Black Holes: From Fundamental Aspects to Latest Developments_, Astrophysics and Space Science Library, Vol. 440, edited by C. Bambi (Springer Berlin Heidelberg, Berlin, Heidelberg, 2016) p. 99.
  * Frank  _et al._ (2002) J. Frank, A. King, and D. J. Raine,  _Accretion Power in Astrophysics, by Juhan Frank and Andrew King and Derek Raine, pp. 398. ISBN 0521620538\. Cambridge, UK: Cambridge University Press, February 2002._ (2002) p. 398.
  * Bakala  _et al._ (2012) P. Bakala, M. Urbanec, E. Šrámková, Z. Stuchlík, and G. Török, Classical and Quantum Gravity 29, 065012 (2012).
  * Ludlam  _et al._ (2019) R. M. Ludlam, J. M. Miller, D. Barret, E. M. Cackett, B. M. Coughenour, T. Dauser, N. Degenaar, J. A. García, F. A. Harrison, and F. Paerels, ApJ 873, 99 (2019), arXiv:1902.00520 [astro-ph.HE] .
  * Niedźwiecki (2005) A. Niedźwiecki, MNRAS 356, 913 (2005), arXiv:astro-ph/0409710 [astro-ph] .
  * Gaspari  _et al._ (2013) M. Gaspari, M. Ruszkowski, and S. P. Oh, MNRAS 432, 3401 (2013), arXiv:1301.3130 [astro-ph.CO] .
  * King and Nixon (2015) A. King and C. Nixon, MNRAS 453, L46 (2015), arXiv:1507.05960 [astro-ph.HE] .
  * Suková  _et al._ (2016) P. Suková, M. Grzedzielski, and A. Janiuk, A&A 586, A143 (2016), arXiv:1506.02526 [astro-ph.HE] .
  * Blinova  _et al._ (2016) A. A. Blinova, M. M. Romanova, and R. V. E. Lovelace, MNRAS 459, 2354 (2016), arXiv:1501.01948 [astro-ph.SR] .
  * Ikhsanov (2001) N. R. Ikhsanov, A&A 367, 549 (2001).
  * González-Galán  _et al._ (2012) A. González-Galán, E. Kuulkers, P. Kretschmar, S. Larsson, K. Postnov, A. Kochetkova, and M. H. Finger, A&A 537, A66 (2012), arXiv:1111.6791 [astro-ph.HE] .
  * Quaintrell  _et al._ (2003) H. Quaintrell, A. J. Norton, T. D. C. Ash, P. Roche, B. Willems, T. R. Bedding, I. K. Baldry, and R. P. Fender, A&A 401, 313 (2003), arXiv:astro-ph/0301243 [astro-ph] .