New approaches to the general relativistic Poynting-Robertson effect

Objectives: A systematic study on the general relativistic Poynting-Robertson effect has been developed so far by introducing different complementary approaches, which can be mainly divided in two kinds: (1) improving the theoretical assessments and model in its simple aspects, and (2) extracting mathematical and physical information from such system with the aim to extend methods or results to other similar physical systems of analogue structure. Methods/Analysis: We use these theoretical approaches: relativity of observer splitting formalism; Lagrangian formalism and Rayleigh potential with a new integration method; Lyapunov theory os stability. Findings: We determined the three-dimensional formulation of the general relativistic Poynting-Robertson effect model. We determine the analytical form of the Rayleigh potential and discuss its implications. We prove that the critical hypersurfaces (regions where there is a balance between gravitational and radiation forces) are stable configurations. Novelty /Improvement: Our new contributions are: to have introduced the three-dimensional description; to have determined the general relativistic Rayleigh potential for the first time in the General Relativity literature; to have provided an alternative, general and more elegant proof of the stability of the critical hypersurfaces.


INTRODUCTION
The last four years have been witness of revolutionary discoveries in astrophysics, which have seen as protagonists these two significant events: (1) the first detection of gravitational waves from the binary black hole (BH) GW151226 [1,2] and then from the binary neutron star (NS) GW170817 [3] thanks to the LIGO and VIRGO collaborations; (2) the first imaging of the matter motion around the supermassive BH in the center of M87 Galaxy [4][5][6][7][8][9] thanks to the strong efforts spent on building the Event Horizon Telescope (EHT) and the synergetic cooperation between EHT and Black Hole Cam project. The achievement of such scientific milestones have increasingly motivated all research groups to improve the actual theoretical models, to validate Einstein theory or possible extension of it, when benchmarked with these new amount of powerful observational data.
Generally, the motion of the matter around massive compact objects, as stellar NSs or BHs, or supermassive BHs, is approximated to be mainly geodetic. However, in view of the actual powerful observational capacities and facilities, it is important to take into account other small perturbing effects. In particular, when we consider the motion of relatively small-sized test particles (e.g., dust grains, gas clouds, meteors, accretion disk matter elements) around electromagnetic radiating sources (like type-I X-ray bursts on NS polar caps, boundary layer around a NS, or a hot corona around a BH) an important effect to be taken into account is the Poynting-Robertson (PR) effect [10,11]. The forces acting on the * vittorio.defalco@physics.cz test particle are: the gravitational field, directed toward the compact object and opposite to the radiation pressure, pointing outward, and also the PR effect. This phenomenon is triggered each time the radiation field invests the test particle, raising up its temperature, which for the Stefan-Boltzmann law starts re-emitting radiation. In this model, the test particle is considered as an ideal black body in thermal equilibrium, meaning that all the absorbed energy is isotropically re-emitted.
This process of absorption and remission of radiation generates a recoil force opposite to the test body orbital motion. This can be interpreted as an aberration effect in the test particle's frame or also as an anisotropic reemission in the star reference frame. It is important to note that radiation pressure and PR effect can be split in the classical frame, while in GR frame they constitute one single function, which must satisfy the relativistic covariance principle, in order not to run into paradoxes. However, the PR effect can be seen as the action of an electromagnetic field on a moving body. Such mechanism removes thus very efficiently angular momentum and energy from the test particle, forcing it to spiral inward or outward depending on the radiation field intensity.
Such effect has been initially introduced in classical physics by Poynting in 1903 [10], and then extended in special relativity by Robertson in 1937 [11], with several applications to the Solar system [12]. From 2009, it has been extended in GR by Bini and collaborators within the two-dimensional (2D) equatorial plane of the Kerr spacetime [13,14]. Recently, it has been extended also in the three-dimensional (3D) space [15][16][17].
Our aim is to investigate such phenomenon under different perspectives, in order to extract more peculiar information, and for developing new valuable mathematical tools, which can be broadly applied to other dissipative systems in GR. in this paper, we would like to revise our new approaches and the consequent implications. The article is structured as follows: in Sec. 2, we concentrate on the PR effect model, showing how to pass from the 2D description to the 3D case, discussing the further implications and new advantages; in Sec. 3, we treat the PR effect as a dissipative system under a Lagrangian formalism, determining analytically the Rayleigh dissipation function through a new procedure; in Sec. 4, we study the PR effect as a dynamical system, proving that the critical hypersurfaces (regions where gravitational and radiation forcese balance) are stable configurations within Lyapunov theory. Finally, in Sec. 5 the conclusions are drawn.

Geometry and strategy
We consider a rotating compact object, whose outside spacetime is described by the Kerr metric. We use the signature (−, +, +, +) for the metric, and geometrical units for gravitational constant G, and speed of light c (c = G = 1). The metric line element, ds 2 = g αβ dx α dx β , expressed in Boyer-Lindquist coordinates, parameterized by mass M and spin a, reads as [18] 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 strategy for determining the equations of motion of a test particle under the general relativistic PR effect is obtained by following the initial approach followed by Poynting and Robertson. We write first the equations in the test particle rest frame through the zero angular momentum observers (ZAMOs), considered at the position occupied by the test particle at each instant of time. Then, we transform them in the static observer frame located at infinity. The successful technique exploited to achieve such objective is reached through the relativity of observer splitting formalism. This represents a powerful method in GR to distinguish the gravitational effects from the fictitious forces arising from the relative motion of two non-inertial observers [19][20][21][22]. Such formalism allows us to derive the test particle equations of motion in the reference frame of the static observer located at infinity as a set of coupled first order ordinary and highly-linear differential equations [13][14][15][16].
The orthonormal frame adapted to the ZAMOs is [13][14][15][16] where N = (−g tt ) −1/2 and N ϕ = g tϕ /g ϕϕ . The nonzero ZAMO kinematical quantities in the decomposition of the ZAMO congruence 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 [15], for their explicit expressions). We denote scalar and tensors measured in the ZAMO frame respectively followed by (n) and by a superposed hat.

Radiation field
The radiation field is modeled as a coherent flux of photons traveling along null geodesics on the Kerr metric. We consider that at each instant of time a single photon, from an emitting surface around the central compact object, reaches the test particle in its position. The related stress-energy tensor T µν is given by [15,16] where I is a parameter linked to the radiation field intensity and k is the photon four-momentum field, where the last two equations express the condition of null geodesic. In Kerr spacetime, we have that the energy E = −k t , the angular momentum with respect to the polar axis L z = k ϕ , the Carter constant Q, and the module of the photons k µ k µ = 0 are conserved along its trajectory. Splitting k with respect to the ZAMO frame (see Fig.  1), we obtain [13][14][15][16] k = E(n)[n +ν(k, n)], (4) ν(k, n) = sin ζ sin β er + cos ζ eθ + sin ζ cos β eφ, (5) where E(n) is the photon energy measured in the ZAMO frame,ν(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. The radiation field in the 3D model is governed by the two impact parameters (b, q), associated respectively with the two emission angles (β, ζ). The radiation field photons are emitted from a spherical rigid surface having a radius R centered at the origin of the Boyer-Lindquist coordinates, and rotating rigidly with angular velocity Ω . In the 2D model, the motion occurs only in the equatorial plane θ = π/2, where we have only one impact parameter b related to the only emission angle β in the ZAMO frame. We say that for b = 0 the radiation field is radial, otherwise it is a general radiation field (see Fig. 1). The photon impact parameters are [14,16] b Both relations are evaluated at the emitting surface radius R . For the azimuthal photon impact parameter b, we have that b = b(R , Ω , θ, a), where θ is the polar angle occupied by the test particle during its motion. Instead, for the latitudinal photon impact parameter q, it fixes the Carter constant Q for a given photon trajectory, and depends only on the test particle's polar angle θ and the value assumed by b.
By imposing that the photon must hit the test particle in the equatorial plane of the ZAMO frame, even at infinity, we substantially reduce the complexity of the problem, because everything is expressed in terms of one single parameter, which is astrophysically realistic. In both 2D and 3D models, from the conservation of the stress-energy tensor conditions (Bianchi identities), namely ∇ µ T µν = 0, we are able to determine the parameter I, which has the following expression [16] where I 0 is I evaluated at the emitting surface.
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. We can split the photon four momentum (4) in terms of the velocity U as [16] where E(U ) is the photon energy measured by the test particle. The radiation force can be written as [16] where m is the test particle mass and the termσ[IE(U )] 2 reads as [16] σ and with A =σ[I 0 E] 2 being the luminosity parameter, which can be equivalently written as 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, whose expressions are [16] Gathering all information together, it is possible to derive the resulting equations of motion for a test particle moving in a 3D space, which are [13][14][15][16] where τ is the affine parameter (proper time) along the test particle trajectory. Naturally, these equations reduce to the 2D case, when ψ = θ = π/2.

Critical hypersurfaces and test particle's trajectories
The general relativistic PR effect, defined by Eqs. (20)-(25), exhibits, both in the 2D and 3D models, a critical hypersurface around the compact object. This is a region, where there exists a balance among gravitational and radiation forces, see Fig. 2 for some examples. Such structures are obtained by the requirement that the test particle moves on them (α = 0, π) with constant velocity (ν = const) with respect to the ZAMO frame, and the polar axis is orthogonal to the critical hypersurface (ψ = ±π/2), which in turn implies that dν/dτ = dα/dτ = 0 [13][14][15][16] where the first condition means that the test particle moves on the critical hypersurface with constant velocity equal to the azimuthal photon velocity, see Eq. (5); whereas the second condition determine through an implicit equation the shape of the critical hypersurface in terms of the critical radius r (crit) as a function of the polar angle, once metric background (i.e., a) and radiation field proprieties (i.e., A, R , Ω ) are assigned. When the test particle reaches the critical hypersurface, we can have to possible behaviors: • latitudinal drift, due to the interplay of gravitational and radiation actions in the polar direction, which brings definitively the test particle on the equatorial plane [15,16]. This corresponds to the condition dψ/dτ = 0, because the ψ angle change during the test particle motion on the critical hypersurface; • suspended orbits, the test particle does not drift down to the equatorial plane, but moves on a pure circular orbit at θ =θ = π/2. This condition is mathematically achieved by imposing that dψ/dτ = 0, which for b = 0 reads as [16] a(n)θ + k (Lie) (n)θ ν 2 + 2ν sin ψ θ(n)θφ which is an implicit equation in terms of ψ. Instead for b = 0 we obtain ψ = ±π/2 [15]. It is important to note that in the Schwarzschild case (a = 0) for b = 0 the test particle stops on a point on the critical hypersurface, without moving on a purely circular orbit [15].
In Fig. 3 we display some selected test particle trajectories in 2D and 3D cases in order to show how the PR effect alters the matter motion around a compact object [13][14][15][16]. It is important to note that in both cases, such effect is strongly sensitive from the initial conditions, therefore it requires that the integration error should be very low, otherwise its propagation could lead to not realistic trajectories. In both models, the test particle has two possible endings, strongly depending on the initial conditions and the parameters defining the geometrical structure and the radiation field: orbiting on the critical hypersurface, or departing at infinity.
where X = (t, r, θ, ϕ) and U = (U t , U r , U θ , U ϕ ), and the unknown functions to be determined are: the Lagrangian function L(X, U ) (including the kinetic energy and all the conservative and generalised forces) and the Rayleigh dissipative potential V (X, U ) (encompassing the dissipative forces). This is a well posed mathematical problem and it is known in the literature as the inverse problem in the calculus of variations [23][24][25].
Our aim is to derive the analytical form of both functions, without recurring to numerical simulations or codes. We note that the Lagrangian function is already known in the classical literature of GR [18], connected with the pure gravitational structure of the background spacetime (i.e., a = 0), and given by The challenging issue is of course to determine the analytical form of the Rayleigh potential, connected to the dissipative effects in GR, where generally the dissipative force is highly-nonlinear, because it strongly couples with the curved geometrical background.

The method
The general relativistic PR effect represents the first system in GR, where we have been able to analytically determine the Rayleigh potential [30? , 31]. To obtain such result, we have exploited two important ideas: the integrating factor to make a differential (semi-basic) oneform closed in its simply connected domain (i.e., exact), and an integration strategy based on the energy dissipated by the system. In this section, we explain into details the followed procedure.
Before to start, we set up some preliminary considerations. The motion of the test particle occurs in M, a simply connected domain (the region outside of the compact object including the event horizon). We denote with T M the tangent bundle of M, whereas T * M stands for the cotangent bundle over M. Let ω : T M → T * M be a smooth differential semi-basic one-form [26][27][28], then the radiation force components (13) can be seen as the components of ω, namely The vertical exterior derivative d V is an operator, . Left lower panel: Test particle trajectories around a NS of spin a = 0.41, radius R = 6M , angular velocity Ω = 0.031, and luminosity parameter A = 0.8, starting at the position (r0, θ0) = (15M, 10 • ) with the initial velocity ν0 = 0.01 oriented in the azimuthal corotating direction direction (orange) and oriented radially towards the emitting surface (red). Right lower panel:Test particle trajectories around a NS of spin a = 0.07, radius R = 6M , angular velocity Ω = 0.005, and luminosity parameter A = 0.85, starting at the position (r0, θ0) = (15M, 10 • ) with the initial velocity ν0 = 0.01 oriented in the azimuthal corotating direction direction (orange) and oriented radially towards the emitting surface (red). whose local expression on a smooth differential semi-basic one-form β = β α dX α is given by [27][28][29] whose components are the cross derivatives of β α with respect to U β and ∧ is the wedge product to lift the forms. The differential semi-basic one-form ω is closed under the vertical exterior derivative d V if d V ω = 0. However, it can be checked that the cross derivatives of (33) are not equal to zero if applied to Eq. (13). Due to the non-linear dependence of the radiation force on the test particle velocity field, the semi-basic one-form turns out to be not exact [22]. However, the PR phenomenon exhibits the peculiar propriety according to which ω(X, U ) becomes exact through the introduction of the integrating factor µ = (E/E) 2 [22,30,31], where For the Poincaré lemma (generalised to the vertical differentiation) the closure condition and the simply connected domain T M guarantee that µω is exact [32]. Therefore, it exists a smooth 0-form V (X, U ) such that −d V V = µω, which in local coordinates reads as Substituting all the occurrences of E in F (rad) (X, U ) α , see Eq. (13), we obtain [30,31] We consider the velocity derivative operator with respect to the energy E through the chain rule, having Equation (34) reads explicitly in such case as [30,31] µF α (rad) = k α ∂V ∂E .
Considering the scalar product of both members of Eq. (37) by U α , we obtain a differential equation for V , which yields at the following integral [30,31] where f (X, U ) is constant with respect to E, i.e., ∂f (X, U )/∂E = 0. Such term can be calculated by using the iterative process of integration of exact differential one-forms, after having rewritten E as −k α U α in order to have all coherently expressed in terms of (X, U ). Such step is fundamental to understand our strategy, because it permits not only to reduce the calculations from four in terms of the velocity field U to just one in terms of the energy E, but permits to analytically obtain the Rayleigh potential in terms of the energy E whenever the calculations to obtain f (X, U ) become complicate.
Integrating Eq. (38), the final result is [30,31] V =σI 2 ln where the constant term 1/2 − ln(E) has been determined considering the classical limit, where it is possible to match the general relativistic solution with the classical description [10,11,30,31].

Discussion of the results
Our new developed strategy is based on two stages.
(1) Use of an integrating factor to make a differential semi-basic one-form exact (closed in its simply connected domain of integration). This approach permits to have more dissipative systems admitting their differential semi-basic one-forms closed, since in GR we deal with dissipative forces highly nonlinear with respect to the velocity field. The closure condition translates in solving a partial differential equation for the integrating factor µ, which in the general relativistic PR case can be easily solved by separation of variables [31]; (2) Writing the dissipative force in terms of the dissipated energy E (see Eq. (35)) and passing the derivative operatore from the velocity field to the dissipated energy E by applying the chain rule (see Eq. (36)), we can finally have through some simple algebraic calculations an analytical form of the Rayleigh potential in terms of the dissipated energy E (see Eq. (38)). This integral is defined up to a constant function f (X, U ) with respect to the dissipated energy E, i.e., ∂f (X, U )/∂E = 0. In this term is contained our ignorance on how the dissipative action occurs in the physical system, but at least we know how to act in energetic terms.
It is important to note that E is considered as the dissipated energy, because we can see our radiation field made by photon-bullets, which are shot against the test particle-target. Faster the test particle moves, less are the photons absorbed, or in other words, more are the photons dissipated. Indeed, if the test particle is at rest, E = E, which is the maximum attained energy; while instead if the test particle moves close to a speed of light, E = 0 is the minimum reachable energy [31].
The analytical form of the Rayleigh potential related to the general relativistic PR effect is very important because it represents the first example in the GR literature. The Rayleigh potential (38) is a valuable tool to investigate the proprieties of the general relativistic PR effect and more in general the radiation processes in highenergy astrophysics. In addition, such function entails two main important consequences: • the Rayleigh potential contains a logarithm of the energy, which, in view of the interpretation of the dissipated energy E, can be physically interpreted as the absorbed energy from the test particle, while the other function U α U α represents the re-emission process, which is in agreement with the underlying hypothesis of the PR effect model, i.e., the test particle behaves as an ideal black body in thermal equilibrium, which re-emits radiation at constant rate, isotropically, and independent from the velocity field U α [31].
• the Rayleigh potential permits to directly connect the theory with the observations. This last sentence can be better understood by looking at Fig.  4, where in panel a) the test particle trajectory is displayed (i.e., what can be observed) and in panels b) − f ) the Rayleigh potential in terms of the coordinates r, ϕ, t,ṙ, rφ, respectively (i.e., what derives from the theory). In other words, detecting the test particle motion, it is possible to infer the analytic structure of the Rayleigh potential; viceversa by using different functional form of the Rayleigh potential to investigate other kinds of radiation processes in high-energy astrophysics, it can be possibile to numerically simulate the test particle trajectory (see Ref. [31], for details).

DYNAMICAL SYSTEM APPROACH: STABILITY OF THE CRITICAL HYPERSURFACES
The general relativistic PR effect can be also analysed under the dynamical system point of view for proving the stability of the critical hypersurfaces. To this end, we consider only those initial conditions, where the test particle ends its motion on the critical hypersurfaces without escaping at infinity. It is not possible to formally characterise these configurations, because the general relativistic PR effect models show a sensitive dependence from the initial data [? ]. Once the stability issue has been proved, it immediately follows that the critical equatorial ring is a stable attractor (region where the test particle is attracted for ending its motion), and the whole critical hypersurface is a basin of attraction [33].

Previous approach: linearization theory
In the previous approach pursued by Bini and collaborators, they have been able to prove such statement only in the Schwarzschild case within the linear stability theory (see Appendix in Ref. [14]). This method relies on the linearization of a dynamical systemẋ = f (x) towards the critical point Then searching for the eigenvalues of the matrix A, and checking whether they are all negatives, we can finally obtain that the dynamical system is stable toward the critical point x 0 . At glance, this approach seems theoretically simple, but practically it implies computationallyexpensive calculations (especially in the Kerr case).

New approach: Lyapunov theory
We have introduced a new, simpler in term s of calculations, and more physical approach based on the Lyapunov theory [34]. The dynamical system (20)- (25), x = f (x), is defined in the domain D (spacetime outside the compact object including the event horizon), and we call with H the critical hypersurface. Let Λ = Λ(x) be a smooth and real valued function, continuously differentiable in all points of D, then Λ is a Lyapunov function forẋ = f (x) if it satisfies the following conditions: Once a Lyapunov function Λ is determined for all points belonging to the critical hypersurface H, a theorem due to Lyapunov assures that H is stable [34]. The conditions (I) and (II) are very easy to be proved, while the last condition (III) requires to perform some calculations, that anyway can be easily carried out by hand. Such approaches has the great advantage to study the behavior of a dynamical system without knowing its analytical solution or performing any approximation. Nevertheless, some limits should be also taken into account, like: there is any fixed rule or recipe to determine it, but sometimes only physical intuitions can help; it is not unique, but it could be possible to determine more than one or sometimes can even not exist.

Three Lyapunov functions for the general relativistic PR effect
For the general relativistic PR effect, three different Lyapunov functions have been determined. An idea of the proof that they are Lyapunov functions is based on expanding all the kinematic terms of the equations of motion with respect to the radius, estimating thus their magnitude, and then considering that the test particle's orbit is confined in a bounded box (since by hypothesis we consider only those configurations reaching the critical hypersurface), see Ref. [33], for further details.
• The relative classical mechanical energy of the test particle with respect to the critical hypersurface measured in the ZAMO frame is where ν crit (θ) = [cos β] r=rcrit(θ) , which includes as a particular case the velocity ν eq = [cos β] r=rcrit(π/2) in the equatorial ring. Its derivative iṡ where sgn(x) is the signum function.
• The relative classical angular momentum of the test particle measured in the ZAMO frame is L = m(rν sin ψ cos α − r crit ν crit ).
Its derivative is given bẏ +r dν dτ cos α sin ψ + ν(ṙ cos α sin ψ −r sin α sin ψα + r sin α cos ψψ) . (46) • The relative general relativistic Rayleigh dissipation function is (see Sec. 3 for details) where E crit is the energy E evaluated on the critical hypersurface, given by Its derivative iṡ In Fig. 5 we show the trajectory of a test particle orbiting in the equatorial plane around a BH and influenced by a radiation field together with the general relativistic PR effect. The test particle ends its motion on the critical hypersurface, which is the configuration we focussed on. In the other panels of Fig. 5, the trend of the three proposed functions have been displayed (i.e., K, L, F), together with their derivatives (i.e.,K,L,Ḟ), to graphically prove that they verify the three proprieties to be Lyapunov functions. We note that the the first two Lyapunov functions (mechanical energy and angular momentum) are classical definition, while the third example (Rayleigh dissipation function) is a pure general relativistic case. The former two functions are not in contradiction with the latter, rather they are very useful to substantially reduce the calculations with respect to their general relativistic versions. In addition, it must be said that another great advantage of Lyapunov theory is that the Lyapunov functions should not need to have a physical meaning, they can be also mathematical function, which should respect the three conditions stated above, so that they can prove the stability of the critical hypersurfaces.

CONCLUSIONS
In this work we have presented three different and complementary approaches to study the general relativistic PR effect. They can be summarised as: • working on the model by improving it in its ingenuous aspects (see Sec. 2). Indeed, we have improved the PR effect model from the 2D equatorial plane to the 3D space in Kerr geometry. The emitted photons are in general parametrized by two impact parameters (b, q), since we are in the 3D case. However, imposing that the radiation field must lie in the equatorial plane of the ZAMO frame (even at infinity), we not only reduce the radiation field to only one parameter, but we also develop an astrophysical coherent model. The resulting equations of motion represent a system of six coupled ordinary and highly nonlinear differential equations of first order, see Eqs. (20)-(25). Such dynamical system admits the existence of a critical hypersurface, regions where the gravitational attraction is balanced by the radiation force. The test particle can end its motion on it, moving over there stably (see Sec. 2 2.4). The main difference of the 3D model with the 2D case is in having the phenomena of latitudinal drift and suspended orbits.
• the general relativistic PR effect can be seen as a dissipative system (see Sec. 3). The equations of motion can be treated within the theory of inverse problem of calculus of variations, where the unknown functions are the Lagrangian and the Rayleigh potential. The former is easily found (31), while for the latter we have developed a strategy to determine it analytically (39). Such approach is based first on making a differential semi-basic one-form exact through the introduction of an integrating factor, and then to substantially reduce the calculations for obtaining its analytic expression at least in terms of the dissipated energy (38). The full complete analytical form of the Rayleigh potential permitted to discover a new functional class related to absorption processes in high-energy astrophysics (see Sec. 3 3.3), and to develop a strategy to closely relate observations and theory (see Fig. 4).
• another perspective to analyse the general relativistic PR effect is under the dynamical system point of view (see Sec. 4). We have proved the stability of the critical hypersurfaces by introducing a new method. The previous approach was based on the linearization theory, where the calculations revealed to be very demanding (see Sec. 4 4.1). Therefore, we have thought to introduce a more straightforward method within the Lyapunov theory (see Sec. 4 4.2). The determination of three different Lyapunov functions (i.e., classical mechanical energy (43) and angular momentum (45) of the test particle and general relativistic Rayleigh po-tential (48)) permitted to prove in an equivalent way and under different physical aspects the stability of the critical hypersurfaces.
As future projects, we plan to further investigate the results obtained in the three different approaches, namely: (1) improve the actual theoretical assessments employed to model the radiation field in some elementary aspects, like: the momemntum-transfer cross section will be not anymore constant, but it will depend on the angle and frequency of the incoming radiation field, the radiation field is not emitted anymore by a point-like source, but from a finite extended source; (2) generalise the strategy for finding the analytical form of the Rayleigh potential to other dissipative systems in GR, like for example the gravitational wave theory; (3) apply the Lyapunov functions to all the extensions of the general relativistic PR effect model with the due modifications.
Anyway we would like to find other alternative approaches to the same problem in order to extract several other interesting information, and in the same time to develop new formalisms, which can be applied to other dissipative systems in GR. We would like also to apply the general relativistic PR effect model to describe some astrophysical problems, like: accretion disk model, type-I X-ray burst, photospheric radius expansion.

Funding and conflict of interest
The author thanks the Silesian University in Opava and Gruppo Nazionale di Fisica Matematica of Istituto Nazionale di Alta Matematica for support. The results contained in the present paper have been partially presented at the conference WASCOM 2019, Maiori (Italy).