Stable attractors in the three-dimensional general relativistic Poynting-Robertson effect

We prove the stability of the critical hypersurfaces associated with the three-dimensional general relativistic Poynting-Robertson effect. The equatorial ring configures to be as a stable attractor and the whole critical hypersurface as a basin of attraction for this dynamical system. We introduce a new, simpler (in terms of calculations), and more physical approach within the Lyapunov theory. We propose three different Lyapunov functions, each one carrying important information and very useful for understanding such phenomenon under different aspects.


INTRODUCTION
In the study of radiation processes in high-energy astrophysics around compact objects, like neutron stars (NSs) and black holes (BHs), it is of utmost importance to accurately describe the motion of the surrounding matter to then test the model in strong field regimes. In particular, when we deal with relatively small-sized test particles, like dust grains, plasma or gas elements invested by an electromagnetic radiation field, their motion can be considerably altered by the general relativistic Poynting-Robertson (PR) effect [1,2]. The radiation field exerts not only a force outward the compact object contrasting the gravitational pull, but also a radiation drag force opposite to the test particle orbital velocity. The PR effect configures thus as a pure relativistic dissipative effect, which efficiently removes energy and angular momentum from the affected test particle.
The general relativistic treatments, from the two dimensional (2D) [3,4] until the three dimensional (3D) formulations [5][6][7], show all the existence of a critical hypersurface, region where gravitational and radiation forces balance. From selected test particle orbits, it graphically results that the critical hypersurfaces are stable, namely once the test particle reaches such region and moves on that, it will remain there forever. The test particle can either spiral down towards the equatorial ring (latitudinal drift) or move in suspended orbits (see Refs. [5,6], for details). This implication must be formally proved, because it contains fundamental information not only on the PR effect, but also on the system under study.
We propose a new method to prove the stability of the critical hypersurfaces within the Lyapunov theory. Such approach carries important information on the physical * vittorio.defalco@physics.cz † pavel.bakala@physics.cz system under study and it substantially reduces the calculations with respect to a previous method employed in the literature [4]. The idea of a Lyapunov function has been proposed in 1956 [8], and since then it has been exploited in several and disparate contexts in physics and in mathematics [9][10][11][12]. There are also a wealth of applications in astrophysics and cosmology, such as: in accretion disc theory to control the large-and small-scale in/stability of such continuous-like structures (see, e.g., Refs. [13][14][15][16]); in celestial mechanics to study the motion of bodies under the influence of a gravitational (or other kinds of) forces (see, e.g., Refs. [17][18][19]); in cosmology to understand the stability of the models, to analyse dark energy's origin and implications, and to investigate modified gravity scenarios (see e.g., Refs [20,21]). The paper is structured as follows: in Sec. 2 we introduce the formal aspects of the 3D general relativistic PR effect model, including a detailed description of the critical hypersurfaces' derivation and proprieties. In Sec. 3 we review what has been done in the literature so far, underlying the limiting aspects, and then we present the power and advantages of our new approach. Finally, in Sec. 4 we draw our conclusions.

GENERAL RELATIVISTIC 3D PR EFFECT 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 αβ dx α dx β , in Boyer-Lindquist coordinates, parameterized by mass M and spin a, reads as where Σ ≡ r 2 + a 2 cos 2 θ, ∆ ≡ r 2 − 2M r + a 2 , and ρ ≡ r 2 + a 2 + 2M a 2 r sin 2 θ/Σ. The determinant of the Kerr metric is g = −Σ 2 sin 2 θ. We introduce the zero angular momentum observers (ZAMOs), whose adapted orthonormal frame is given by where where {∂ t , ∂ r , ∂ θ , ∂ ϕ } is the orthonormal frame adapted to the static observer at infinity, N = (−g tt ) −1/2 is the time lapse function and N ϕ = g tϕ /g ϕϕ the spatial shift vector field. All the vector and tensor indices (e.g., v α , T αβ ) associated to the ZAMO frame will be labeled by a hat (e.g., vα, Tαβ), instead all the scalar quantities measured in the ZAMO frame (e.g., f ) will be followed by (n) (e.g., f (n)). In the kinematical decomposition of the ZAMO congruence, we have that the nonzero ZAMO kinematical quantities are acceleration a(n) = ∇ n n, expansion tensor along theφ-direction θφ(n), and the relative Lie curvature vector k (Lie) (n) (see Table 1 in [5], for their explicit expression). 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 [5,6] where I is a parameter linked to the radiation field intensity and k is the photon four-momentum field. Splitting k with respect to the ZAMO frame, we obtain ν(k, n) = sin ζ sin β er + cos ζ eθ + sin ζ cos β eφ, (5) whereν(k, n) is the photon spatial unit relative velocity with respect to the ZAMOs, β and ζ are the two angles measured in the ZAMO frame in the azimuthal and polar direction, respectively. In addition, E(n) is the photon energy measured in the ZAMO frame, which is [5,6] where E p = −k t is the conserved photon energy along its trajectory. The radiation field is governed by the two photon impact parameters (b, q), associated respectively with the two emission angles (β, ζ). The photons of the radiation field are emitted from a spherical surface having a radius R centered at the origin of the Boyer-Lindquist coordinates, and rigidly rotating with angular velocity Ω . The photon impact parameters have the following expressions [6] b = − g tϕ + g ϕϕ Ω g tt + g tϕ Ω r=R , The related photon angles in the ZAMO frame are [6] cos β = bN √ g ϕϕ (1 + bN ϕ ) , ζ = π/2.
The parameter I is given by [6] where I 0 is I evaluated at the emitting surface. A test particle moves with a timelike four-velocity U and a spatial three-velocity with respect to the ZAMOs, ν(U, n), which both read as [5] ν = ν(sin ψ sin αer + cos ψeθ + sin ψ cos αeφ), (12) where γ(U, n) ≡ γ = 1/ 1 − ||ν(U, n)|| 2 is the Lorentz factor, ν = ||ν(U, n)||, γ(U, n) = γ. We have that ν represents the magnitude of the test particle spatial velocity ν(U, n), α is the azimuthal angle of the vector ν(U, n) measured clockwise from the positiveφ direction in thê r −φ tangent plane in the ZAMO frame, and ψ is the polar angle of the vector ν(U, n) measured from the axis orthogonal to ther−φ tangent plane in the ZAMO frame.
We assume that the radiation test particle interaction occurs through Thomson scattering, characterized by a constant momentum-transfer cross section σ, independent from direction and frequency of the radiation field. The radiation force is given by [5,6] where the termσ[IE(U )] 2 reads as [5,6] (14) E(U ) is the photon energy absorbed by the test particle, which can be related to E(n) through [5,6] The term A =σ[I 0 E p ] 2 is the luminosity parameter, which can be equivalently written A/M = L/L EDD ∈ [0, 1] with L the emitted luminosity at infinity and L EDD the Eddington luminosity. The termsV(k, U ) α are the radiation field components, which are [5,6] Collecting all the information together, it is possible to derive the resulting equations of motion for a test particle moving in a 3D space, which are [5,6] +2ν cos α sin ψ θ(n)rφ + cos ψ a(n)θ +2ν cos α sin 2 ψ θ(n)θφ − sin α cos ψ a(n)r where τ is the affine (proper time) parameter along U.

Critical hypersurfaces
The dynamical system defined by Eqs. (20)-(25) exhibits an axially symmetric hypersurface outside around the compact object, where there exists a balance among gravitational and radiation forces. We impose that on such region the test particle must move in purely circular orbits with constant velocity (ν = const) with respect to the ZAMO frame (α = 0, π), and the polar axis orthogonal to the critical hypersurface (ψ = ±π/2). These requirements entail dν/dτ = dα/dτ = 0, from which we determine the following conditions [5,6] where Eq. (27) means that the test particle moves on the critical hypersurface with constant velocity equal to the azimuthal photon velocity; whereas Eq. (28) determine the critical radius r crit as a function of the polar angle θ through an implicit equation, once the radiation field proprieties are assigned (i.e., the radius R and the angular velocity Ω of the emitting surface together with the luminosity parameter A, see Ref. [6]). It is important to note that Eq. (28) might admit three different solutions with precise locations: one inside the emitting surface (considered to be non physical), one close to the emitting surface (the solution we actually study and plot), and another one very far from the emitting surface (that we do not take into account) [4,6]. In general we have dψ/dτ = 0, because the ψ angle change during the test particle motion on the critical hypersurface, having what we termed latitudinal drift. This effect, occurring for the interplay of gravitational and radiation actions in the polar direction, brings definitively the test particle on the equatorial plane [5,6]. For ψ = θ = π/2, we have dψ/dτ = 0, corresponding to the equatorial ring. However, we can have dψ/dτ = 0, also for θ =θ = π/2, having what we termed as suspended orbits. The condition for this last configuration in the case of b = 0 can be expressed as [6] which permits to be solved in terms of ψ. Instead in the case of b = 0 we obtain ψ = ±π/2 [5,6]. The critical points are either the suspended orbits or the equatorial ring, where the test particle ends its motion. In the Schwarzschild case, Eq. (29) is an identity, because the test particle either stops on a point (for b = 0) or move on a purely circular orbit in the equatorial ring (for b = 0) of the critical hypersurface (see Refs. [3][4][5][6], for details).
The test particle comoves with the local corototaing observer frames in bound quasicircular orbits and in the equatorial plane in circular orbits [4,5]. In the Kerr case, the critical hypersurface assumes a quasi-ellipsoid shape, depending on the radiation emitting source parameters (R , Ω , A). Therefore, the critical radius is a function of r crit = r crit (θ, A, R , Ω ). In the Schwarzschild case instead, the critical hypersurface becomes a sphere, where such radius can be easily computed at the equator through Eq. (2.33) of Ref. [4]. However, it is important to note that the non-vanishing angular velocity of the emitting surface Ω = 0 breaks the spherical symmetry of the Schwarzschild metric [6], so the 2D model is not anymore valid, and only the 3D case must be employed.
In addition, there are also the conditions dr/dτ = dθ/dτ = 0 (based on the critical hypersurface definition). This means that at fixed radius r, the polar angle θ is constant (dθ/dτ = 0), while at fixed polar angle θ the radius r is constant (dr/dτ = 0). In other words, such conditions require that each parallel section to the equatorial plane intersecting the critical hypersurface gives a circular ring. For the azimuthally rotational simmetry, the critical hypersurfaces are surfaces of revolution symmetric with respect to the polar axis (see Fig. 1).

STABILITY OF THE CRITICAL HYPERSURFACES
For performing our proof, we set before our notations. The dynamical system (20)-(24) is represented by 1 whereẋ = dx/dτ is defined in D, representing the spatial region outside the compact object together with the admissible velocity field field of the test particle, namely The function Φ : D×R −→ D is called the flow associated to the dynamical system (30). Called x i the initial conditions, we have that the solution at timeτ , indicated by x(τ ), can be also written as We consider only those initial configurations, where the test particle ends its motion on the critical hypersurface without escaping at infinity. Unfortunately, it is not simple to mathematically characterize such class of solutions, because the dynamical system under study is too complex. In addition, there is a strong dependence not only on the input parameters determining the radiation field, but also on the test particle initial data. Indeed, this dynamical system shows a sensitive dependence on the initial conditions, whose propriety is to exhibit extremely different behaviors with only tiny changes in the initial conditions 2 (see Fig. 2, as example).
Once the stability has been proven, it immediately follows that the critical equatorial ring is a stable attractor, i.e., a region toward which the test particle tends to be attracted and to end its motion. More formally, a set A is an attractor for the dynamical system (30) if [11,22] 1 We do not include Eq. (25) because it is not vanishing at the critical hypersurface. In addition, being the dynamical system symmetric with respect to ϕ-rotations, it is possible to remove such equation without any lost of generality. Equation (26) is also not added, because it expresses only the conversion between proper time τ and coordinate time t. 2 In Ref. [5], it is explicitly stated that the 3D case is more sensitive from the initial conditions with respect to the 2D case. Indeed, to solve such issue there are more controls on the integration process reaching thus an average relative accuracy of ∼ 10 −14 .

1.
A is forward invariant under Φ τ , namely if x 0 ∈ A then also Φ τ (x 0 ) ∈ A for all τ > 0; 2. there exists a neighborhood of A, called the basin of attraction for A, denoted by B(A), which consists of all points x that enter A in the limit τ → ∞. More formally, B(A) is the set of all points x in the phase space with the following property: for any open neighborhood U of A, there is a positive constant T such that f (t, x) ∈ U for all real τ > T ; 3. there is no proper (non-empty) subset of A having the first two properties. In other words, A = ∩ τ ≥0 Φ τ (U), where U is a basin of attraction.
The first and second propriety are basically linked to the proof of the stability of the critical hypersurface, while the third shows that the critical hypersurface H = A, because once the test particle moves on H, it will not leave such region, configuring thus as the smallest basin of attraction for the dynamical system (30) 3 . Another important propriety of the critical hypersurfaces is they are compact sets (see Eq. (B5), for more details).
In the next sections, we show how to formally prove the stability of the critical hypersurfaces by recalling what has been done in the literature (see Sec. 3 3.1) and then by introducing our new contributions (see Sec. 3 3.2).

Linear stability theory
Bini and collaborators [3,4] have previously presented the proof of the stability by performing calculations within the linear stability theory [11]. This method consists in linearizing the non-linear dynamical system towards the critical points of the critical hypersurface, i.e., where A = (∂ x f )(x 0 ) is a linear operator. Then, after having diagonalized the matrix A, one looks at its eigenvalues and check whether they are negative or have real part negative for inferring the stability of the critical points. Such procedure holds whenever the critical points are not hyperbolic, meaning that the matrix A has no eigenvalues with real part equal to zero. Indeed, the Hartman-Grobman theorem (or also known in the literature as linearization theorem) states that [23,24] there exists a neighborhood U of x 0 and a homeomorphism h ≡ A · (x − x 0 ) : U → R n with h(x 0 ) = 0 such that in the neighborhood U the dynamical system (30) is topologically homeomorphic to the dynamical system (33) through the map h.
For the difficulty of calculations Bini and collaborators have only shown the stability of the 2D critical hypersurfaces in the Schwarzschild case (see Appendix in Ref. [4]). This method, albeit simple in its theoretical explanation, practically requires to develop several calculations (especially in the Kerr case), therefore we are looking for a new, simpler, and more physical approach.

Lyapunov function
We propose a new method framed within the Lyapunov theory [11], which is easier both in terms of calculations and gives more physical insight into the problem under investigation. Let Λ = Λ(ν, ψ, α, r, θ) be a real valued function of the test particle position and velocity fields, continuously differentiable in all points outside of the compact object. Then Λ is a Lyapunov function for the set H, if it fulfills the following conditions: Once the Lyapunov function Λ has been found for all points belonging to the critical hypersurface H, a theorem due to Lyapunov assures that H is stable [11,25]. In addition, if the third condition (34) is replaced by then H is asymptotically stable [11,25]. The great advantage of such a method relies on the fact that it can be applied without solving the differential equations (30). It is important to note that the Lyapunov function is not unique at all, indeed there could be the cases where it is possible to find just one, more than one, or even anything. Unfortunately, there is no a mathematical recipe for determining a Lyapunov function, because it is usually a matter of ingenuity, trials, or luck in each case. However, sometimes there are natural functions to try, like for example the associated first integrals (see Ref. [11], for examples). In our case we were able to determine three different Lyapunov functions, having all important physical meanings.

Energy
We propose as first Lyapunov function the test particle relative energy, where both kinetic and potential energies related to the radiation pressure and the gravitational force are all measured in the ZAMO frame, i.e., where m is the the test particle mass and ν crit (θ) = [cos β] r=rcrit(θ) , which includes as a particular case the velocity in the equatorial ring ν eq = [cos β] r=rcrit(π/2) . By definition it is defined positive outside the critical hypersurface, because the second term is the product of two negative terms (since A/M ∈ [0, 1] and r ≥ r crit ), satisfying thus condition (I), and it is identically zero on the critical hypersurface, as requested by condition (II). Regarding the condition (III) we need to calculate the τ -derivative of K, which is given bẏ where sgn(x) is the signum function. In Appendix C we prove thatK is definite non-positive. In Fig. 3 we show an example of test particle trajectory and the related mechanical energy, together with its τ -derivative (see upper right panel). As we can graphically see, K fulfills all the conditions to be a Lyapunov function.

Angular momentum
Now, we consider as second Lyapunov function the relative angular momentum of the test particle measured in the ZAMO frame, i.e., L = m(rν sin ψ cos α − r crit ν crit ).
By definition it is defined positive outside the critical hypersurface, satisfying thus condition (I), and it is identically zero on the critical hypersurface, as requested by condition (II). Regarding the condition (III) we need to calculate the τ -derivative of L, which is given bẏ ν(ṙ cos α sin ψ − r sin α sin ψα + r sin α cos ψψ) .
In Appendix D we prove thatL is definite non-positive. In Fig. 3 we show the angular momentum and its τderivative behaviors (see lower left panel). Therefore, graphically we see that also the function L fully respects the conditions to be a Lyapunov function.

Rayleigh potential
In the study of the general relativistic PR effect, it has been proved it admits a Lagrangian formulation, albeit it is a highly-non-linear dissipative system in GR [26]. This can be realised through the help of an integrating factor µ = E 2 p /E 2 [27,28], where E p is the photon energy and E ≡ E(U ) = −k α U α , see Eq. (15).
Using the energy formalism [27] it was possible to determine the explicit formula for the Rayleigh potential F related to the radiation force µF (rad) (U ) α (see Eq. (13)). Therefore, the third Lyapunov function is the relative Rayleigh potential 4 , i.e., where E crit is the energy E evaluated on the critical hypersurface, and it is simply given by where the subscript "crit" means to evaluate a quantity at the critical hypersurface. The energy absorbed by the test particle E is always minor than the photon energy E p , and its maximum is attained at E p only when the test particle is in rest [27]. Therefore, it is everywhere positive, satisfying thus the condition (I), and the condition (II) on the critical hypersurface. The τ -derivative of F iṡ In Appendix E we prove thatḞ is definite non-positive, verifying thus condition (III). In Fig. 3 we show how the profile of the Rayleigh potential respects the conditions to be a Lyapunov function (see lower right panel).

CONCLUSIONS
In this paper, we have shown how to formally prove the stability of the critical hypersurfaces for the general relativistic (both 2D and 3D) PR effect models. We have proposed a new approach based on the Lyapunov functions, which is more elegant, easier in terms of calculations, and also contains important physical information on the system under study. Previously in the literature, it has been exploited the linearized theory around the critical points, but it revealed to be not so powerful, because it requires strong computational efforts in linearizing the dynamical system and then in finding the eigenvalues and studying the sign, especially for the 3D case. For such reasons, only the Schwarzschild (2D) case has been completely proved, since in the Kerr spacetime everything becomes extremely more difficult to handle [4].
Our method is able to fully prove the stability issue, without recurring to any symbolic program or dedicated software to carry out our analysis. All the estimations and calculations reported in this paper can be relatively easily handled. The Lyapunov theory applied to our problem reveals to be very ingenious, clear in the calculation process, and more powerful than a numericalprogramming approach. In addition, we understand the contribution played by each single term present in the dynamical equations (see Appendices C, D, E).
In addition, we were able to find and propose three different Lyapunov functions with a well precise meaning, proving thus the stability in different ways. Since the PR effect removes energy and angular momentum from the test particle, those were the inspiring ideas, which led us to build up the first two Lyapunov functions (see Secs. 3 3.2 3.2.1 and 3 3.2 3.2.2). They represent the classical version and not the general relativistic expression. Even if we used not the proper definition, they permit to easily carry out the calculations and achieve the stability results. There is no contradictions with the definition of Lyapunov function and its application. We note that even a mathematical function, with no physical meaning connected with the system under study, but verifying the Lyapunov conditions would be a good candidate for proving the stability of the critical hypersurfaces.
The third Lyapunov function is less intuitive than the previous cases, because it stems out from the Lagrangian approach to the general relativistic PR effect [26][27][28]. Indeed, employing the energy formalism, it is possible to analytically derive the general relativistic Rayleigh potential, which contains the radiation field absorption processes affecting the test particle motion. Such a function involves the logarithm of the test particle absorbed energy. We thought about this function by looking at the plots and physical meaning reported in Refs. [28]. Through these valuable results, we have understood that the critical hypersurfaces are basin of attraction and the equatorial ring is a stable attractor.
This approach shows also another great potentiality, because it can be exploited to prove the stability of the critical hypersurfaces in further possible extensions of the general relativistic PR effect model, naturally with the due modifications, and still keeping its good performances. In a future work, we would like to deeply investigate the proprieties of this effect under the dynamical system point of view. Such kind of research, which has never been studied or proposed so far in the literature, will shed new light on the PR effect. The difficulties of such approach rely on mathematically formalising some notions used intuitively so far. These efforts will permit to develop new methods and techniques apt to infer not only new interesting results on the PR effect, but more in general on dissipative systems in GR. FIG. 3. We show a test particle orbit and the related three Lyapunov functions, which graphically prove the stability of the critical hypersurface. Upper left panel: test particle moving around a rotating compact object with mass M = 1, spin a = 0.3, luminosity parameter A = 0.2, and photon impact parameter b = 0. The test particle starts its motion at the position (r0, ϕ0) = (30M, 0) with velocity (ν0, α0) = ( M/r0, 0). The critical hypersurface is a circle with radius rcrit = 2.07M . The energy (see Eqs. in a box during its evolution, which means to impose the following reasonable limits on the test particle position and velocity parameters 5 : r ∈ [R ,R], ϕ ∈ [0, 2π], θ ∈ 0, π 2 , ν ∈ [0, 1], α ∈ [0, −π] , ψ ∈ 0, π 2 , (A1) 5 We prove the stability for test particle set outside the critical hypersurface. The proof can be easily extended also to a test particle set inside the critical hypersurface (but outside the emitting surface), see also the argument of footnote 3.
The terms γ sin α sin ψ and sin α sin ψ γ − 1 γ are both non-positive. We distinguish two cases: • if ν > x, it is obvious that A is non-positive, because it is the sum of three negative terms; • if ν ≤ x, it is not evident, but we need to perform some calculations, where we have It is important to note that the sign function sgn(ν 2 − cos 2 β) does not change the final sign ofK, because the former multiplies only the kinetic term, while the negative dominant contribution from the gravitational and radiation potential persists with its own sign. Therefore, we have finally proved thatK ≤ 0, both on the critical hypersurface and at the equatorial ring.