Analytical study of light ray trajectories in Kerr spacetime in the presence of an inhomogeneous anisotropic plasma

We calculate the exact solutions to the equations of motion that govern the light ray trajectories as they travel in a Kerr black hole's exterior that is considered to be filled with an inhomogeneous and anisotropic plasmic medium. This is approached by characterizing the plasma through conceiving a radial and an angular structure function, which are let to be constant. The description of the motion is carried out by using the Hamilton-Jacobi method, that allows defining two effective potentials, characterizing the evolution of the polar coordinates. The elliptic integrals of motion are then solved analytically, and the evolution of coordinates is expressed in terms of the Mino time. This way, the three-dimensional demonstrations of the light ray trajectories are given respectively.


II. LIGHT PROPAGATION IN AN INHOMOGENEOUS PLASMA IN KERR SPACETIME
To study the light propagation in a medium, we need to consider the phase space dynamics, by introducing an appropriate Hamiltonian. For light rays propagating in a non-magnetized pressure-less plasma, this Hamiltonian is characterized as in the coordinate system x ≡ x µ = (x 0 , x 1 , x 2 , x 3 ), where p ≡ p µ = (p 0 , p 1 , p 2 , p 3 ) is the canonical momentum covector.
Here, ω p (x) is the plasma electron frequency given by (for a review see Ref. [52]) in terms of the electron charge e, electron mass m e , and the number density N p . Defining the plasmic refractive index [53] with u as the observer's four-velocity, we can recast the Hamiltonian as Here, the quantity p ν u ν = −ω gives the effective energy of the photons of frequency ω, as measured by the observer. This way, the refractive index n ≡ n( x, ω) constructs the optical metric o µν of the plasmic medium, such that H = 1 2 o µν p µ p ν . In fact, the Hamilton-Jacobi approach, requires the contribution of the Jacobi action S, given by with λ as the curve parametrization. Now, if the observer is located on the x 0 curves, it has the four-velocity u µ = δ µ 0 / √ −g 00 , which, employing Eqs. (5), yields the Hamilton-Jacobi equation as [53,54] in which, we have used the identity p ν u ν = p 0 / √ g 00 = −ω. Therefore, the general equations governing the light ray trajectories

A. Light propagating in Kerr spacetime
The Kerr black hole spacetime is described by the line element ds 2 = − 1 − 2M r ρ 2 dt 2 + ρ 2 ∆ dr 2 + ρ 2 dθ 2 + sin 2 θ r 2 + a 2 + 2M ra 2 sin 2 θ ρ 2 dφ 2 − 4M ra sin 2 θ ρ 2 dtdφ, which, here, is given in terms of the Boyer-Lindquist coordinates x µ = (t, r, θ, φ). In the above line element with M and a = J/M being, respectively, the mass and the spin parameter of the black hole, where J is the black hole's angular momentum. The Kerr black hole spacetime admits for a Cauchy and an event horizon (notated respectively by r − and r + ), whose surfaces are determined by solving ∆ = 0, and are given by The other significant hypersurface in Kerr spacetime, where g tt = 0, corresponds to a region, inside which, no static observes can exist and stationary observers are in the state of corotation with the black hole. This surface of static limit is located at and together with the event horizon, forms the black hole's ergosphere in the region r + < r < r SL . Throughout this paper, we restrict our study to the domain of outer communications, i.e., the domain outside the event horizon (r > r + ), and we consider the case of a 2 ≤ M 2 , that corresponds to a black hole rather than a naked singularity. The famous method of separation of the Jacobi action, is based on the definition [5,55] in which E, L and m are, respectively, the energy, angular momentum, and the mass associated with the test particles (here, m = 0). The first two parameters are known as the constants of motion which are specified by the Hamilton equations. Along with the method of separation of the Hamilton-Jacobi equation introduced in Ref. [43], we recast the Hamiltonian (1) in the Kerr spacetime as Assuming ω p ≡ ω p (r, θ), it is then straightforward to see that ∂H/∂t = 0 = ∂H/∂φ = 0. Therefore, taking into account Eqs. (5a) and (12), we can specify the constants of motion as Physically, the angular momentum component p φ corresponds, to the axial symmetry of the spacetime. The nature of ω 0 , on the other hand, becomes clear if we specify the light rays' frequency. In fact, the special case of ω p (r, θ) corresponds to the three constants of motion, H = 0, ω 0 and L. Accordingly, we can apply the Carter's method of separation of the Hamilton-Jacobi equation, through which, the light ray trajectories are given in terms of integrable equations [5,55]. Using Eqs. (13) and (14), the Hamilton-Jacobi equation, H(x, p) = 0, yields The separability property of the Hamilton-Jacobi equation, demands that Eq. (15) can be divided into separated r-dependent and θ-dependent segments. In Ref. [43], this has been made possible by defining for some functions f r (r) and f θ (θ). Now, the identity together with Eqs. (15) and (16), separates the Hamilton-Jacobi equation as in which, Q is the so-called Carter's constant. Recasting the above equations, yields Insertion of Eqs. (19) and (20) into the Hamiltonian (13), and then using Eq. (7a), provides the first order differential equations of motion as where Finally, defining the dimension-less Mino time 1 , γ, as ρ 2 dγ = M dλ [56], the equations of motion are now rewritten as In the next section, we continue our discussion by analyzing the above equations of motion, in order to find analytically exact solutions to the light ray trajectories travelling in an inhomogeneous anisotropic plasma around the black hole.

III. GENERAL ANALYSIS OF THE EQUATIONS OF MOTION
The separation condition in Eq. (16) characterizes the plasma, regarding its geometric distribution in the spacetime surrounding the black hole. Therefore, the characteristic functions f r (r) and f θ (θ), play important roles in defining the plasma's configuration. In this section, we confine the light rays to an inhomogeneous anisotropic plasma, by adopting the particular choices This way, Eq. (16) can be recast as 1 Note that, in the geometric units we use, the parameter time has the dimensions of length. Same holds for the curve parameter λ. where R is the mean radius of the gravitating object, and Ω 0 is a positive constant. Within the text, we use, frequently, the conventions in order to simplify the analysis. In what follows, the temporal evolution of the spacetime coordinates are analyzed separately for this plasma, and the relevant exact solutions are presented.
A. The evolution of the radial distance (the r-motion) In the study of particle trajectories in curved spacetimes, it is of crucial importance to know how the particles approach and recede from the source of gravity. Based on the nature of the interactions, this study is traditionally done by calculating the radial, effective gravitational potential, that acts on the particles [57]. Hence, to scrutinize the r-motion for the light ray trajectories in the context under consideration, we focus on the radial equation of motion (26), and rewrite the expression in Eq. (25a) as where P(r) = r 4 + a 2 r 2 + 2M a 2 r, and the radial gravitational potentials are given by taking into account the condition (30). The negative branch is not of our interest, since it has no classical interpretations 2 . We therefore, choose the positive branch of Eq. (36) as the effective potential, i.e.
In Fig. 1, this effective potential has been demonstrated for different values of a and L. According to the diagrams, no stable orbits are expected since the effective potentials do not possess any minimums. Note that, the photons can also pursue a motion along the opposite direction of the black hole's spin, which is the significance of the retrograde motion. As seen in the left panel of Fig. 1, the effective potential corresponding to the retrograde motion (plotted for a < 0), exhibits a lower maximum energy for the the same angular momentum. Therefore, photons on the retrograde motion encounter a remarkably smoother gravitational potential.
A typical effective potential, plotted for a = 0.8M , L = 1M , and Q + fr = 10M 2 (which are the values that will be taken into account for all the forthcoming diagrams in this paper). The categorization of orbits is done by comparing the photon energy ω0 with that for photons on the UCO, i.e. ωU (for the above values, ωU ≈ 0.755232). Photons with the energy ω0 > ωU , will experience an inevitable fall onto the black hole's event horizon. For the special case of ω0 < ωU , the photons encounter two turning points rD and rF .
The possible trajectories are then categorized based on the turning points. However, before proceeding with the determination of the turning points, let us rewrite Eq. (26) as in terms of the characteristic polynomial where and η r = f r /ω 2 0 .

The turning points
Let us consider a typical effective potential as demonstrated in Fig. 2. The possible types of motion are categorized regarding the photon frequency (energy) ω 0 , compared with its value ω U , for photons on the unstable circular orbits (UCO). When ω 0 > ω U , the characteristic polynomial P(r) has four complex roots, which for ω 0 = ω U , reduces to two complex roots and a (degenerate) positive real root. This latter corresponds to the UCO at the orbital radius r U . For the case of ω 0 < ω U , the polynomial has two complex and two positive reals roots, r D and r F , that correspond to the turning points of the photon orbits. In fact, Eqs. (26) and (29) also inform about the characteristics of the turning points. Defining the coordinate velocity v c (r) = dr/dt, then the turning points, r t , are where v c (r t ) = 0. Based on the definition given in Eq. (37), this condition is equivalent to P(r t ) = 0, which together with Eq. (38), results in the two radii and Ξ = 2 in which, The turning points r D and r F correspond to different fates for the trajectories. Respecting Fig. 2, photons that approach the black hole at r D , will deflect to infinity by pursuing a hyperbolic type of motion (deflection of the first kind (DFK)). The equation of motion (37), if solved at the vicinity of the deflection point r D , yields the DFK as (appendix A) in which where ℘(x) ≡ ℘(x;g 2 ,g 3 ) is the elliptic Weierstraßian ℘ function (see appendix B), here, associated with the invariants given that α 1 = C 2 /C 3 , α 2 = C 1 /C 3 and α 3 = C 0 /C 3 , and the respected coefficients defined as On the other hand, r F is the point of no return, from which, the photons will be captured bu the black hole (deflection of the second kind (DSK)). Applying the same method we pursued to calculate the DFK, at the approaching point r F , the DSK is found to obey the equation where with the Weierstraß coefficients having the same form of expressions as those in Eqs. (47), andα 1 =C 2 /C 3 ,α 2 =C 1 /C 3 and α 3 =C 0 /C 3 , that relate respectively to the same coefficients in Eqs. (48), by doing the exchange r D → r F . In Fig. 3, the DFK and DSK have been plotted for the light rays approaching the black hole, for definite dynamical parameters. As it is expected, the DFK initiates from the turning point r D and escape to infinity, and can therefore, reach a distant observer. This concept, if treated for the trajectories with angular components, has the significance of bending of light in curved spacetimes. It is worth notifying the difference between the two types of the DFKs, as indicated in panel (a) of the figure. As it is expected, the more ω 0 increases toward its critical value ω U , the more the trajectories are inclined towards the black hole. In order to find a switching value for ω 0 , at which the DFKs change their deflecting character, let us recall that r D ≡ r D (ω 0 ) and r F ≡ r F (ω 0 ), whose behaviors have been plotted in Fig. 4. As it is observed in the figure, there is a region of the steepest declination of r D , located at the vicinity of ω U . This region starts in accordance with a switching value ω e 0 , from which, the trajectories start to change their deflecting character. On the other hand, the DSK starts from r F and ends in falling onto the singularity. Hence, light rays that are engaged in this process, will never go beyond the distance r F from the black hole and cannot reach an observer at infinity. Note that, in the case of ω 0 = ω U , the light rays will approach at the point r F < r U < r D , that satisfy the equation V (r U ) = 0. As mentioned above, this corresponds to the UCO. In Fig. 5, this kind of orbit has been plotted, regarding its first and the second kinds (UCOFK and UCOSK), by exploiting the equations of motion obtained above, applied for the radial distance r U . The UCOFK, in particular, is responsible for the formation of the black hole shadow.

The capture zone
In addition to the photons with ω 0 < ω U that approach the black hole from the radial distance r F , those incident photons with ω 0 > ω U , also become completely unstable in the region dominant by the effective potential, and are captured by the black hole  (see Fig. 2). The form of the equation of motion for such photons, is the same as those in Eqs. (45) and (49), but the point of approach can be any point r I > r + . In Fig. 6, an example of this kind of orbit has been plotted.
B. The evolution of the polar angle (the θ-motion) The θ-motion is governed by Eqs. (27) and (25b), for which, the condition (31) implies that can be recast as where is the angular gravitational potential felt by the light rays in the plasma. Essentially, this potential has the general form as defined, for example, for the case of null geodesics in Ref. [59]. This general form can be recovered by letting f θ = 0 in Eq. (53). We defineη for more convenience. Now, performing the change of variable z = cos θ, Eq. (27) can be recast as where and the condition Θ z > 0 is required. Clearly, the characteristics of the motion depend directly on the nature of Θ z . Accordingly, we assess the equation of motion (55), separately, for the casesη > 0,η = 0, andη < 0. These cases completely categorize the types of the θ-motion.

The case ofη > 0
In this case, the expression in Eq. (56) remains unchanged, and the condition Θ z > 0, confines the z parameter in the domain −z min ≤ z ≤ z max , where in which, χ 0 = ξ 2 +η − a 2 > 0. The mentioned domain, corresponds to the polar range θ min ≤ θ ≤ θ max , given that θ min = arccos (z max ) and θ max = arccos (−z max ). This range defines a cone, to which, the test particles' motion is confined.
Having this in mind, we can solve the equation of motion (55) by direct integration, resulting in where with and are the respected Weierstraß invariants. Using the analytical solution (58), the temporal evolution of θ(γ) for the case ofη > 0 has been shown in Fig. 7, together with the behavior of W (θ). As it can be observed from the behavior of W (θ), the values of energy that rely in the region ω 0 ≤ ω U are allowed. The light rays, therefore, can opt all kinds of possible orbits that were discussed in the previous section. Accordingly, and applying the analytical solutions of the radial coordinate, the cross-sectional behaviors of the above orbits have been plotted in Fig. 8, for the case ofη > 0, inside the cone of confinement in the z-x plane (i.e. for φ = 0).

The case ofη = 0
The parameter z, in this case, is confined to the domainz min ≤ z ≤z max withz min = 0 andz max = 1 − (ξ/a) 2 . This domain corresponds to the coneθ min ≤ θ ≤ π/2, whereθ min = arccos (z max ). In this case, the effective angular potential (53) takes the form It is straightforward to see that W (θ) = 0 givesθ 0 = π/2 =θ max , according to which, the minimum of the angular effective potential, W min = W (π/2), is obtained. One can therefore infer that the case ofη = 0 also allows for a stable polar equatorial which implies a 2 > ξ 2 , and we have defined κ 1 = ω 0 a 2 − ξ 2 /M . Accordingly, the temporal evolution of θ(γ) and the corresponding angular effective potential have been plotted in Fig. 9. As it can be observed, in contrast with the case ofη > 0, the allowed energies for the case ofη = 0 are higher than their critical value (i.e. W min > ω 2 U ), and therefore, only the capturing trajectories are possible, which in Fig. 10, has been demonstrated inside the cone of confinement.

The case ofη < 0
Under this condition, the expression in Eq. (56) takes the form and Θ z > 0 requires |η| + a 2 − ξ 2 > 0, which is satisfied inside the domainz min ≤ z ≤z max , wherē (66b) The corresponding particle-cone is therefore confined toθ min ≤ θ ≤θ max , whereθ min = arccos (z max ) andθ max = arccos (z min ), and the respected effective potential is Once again, to determine the possible stable polar orbits we solve W (θ) = 0, which yields This value satisfiesθ min ≤θ 0 ≤θ max , and corresponds to the minimum of the angular effective potential for the case ofη < 0.
To find the analytical solution for θ(γ), we pursue the same method as in the case ofη > 0, that provides where and the corresponding Weierstraß invariants arē As in the previous cases, we have plotted the respected angular effective potential and the temporal evolution of the θ-coordinate, in Fig. 11, for the case ofη < 0, which similar to the case ofη = 0, indicates that W min > ω 2 U . Hence, only the capturing trajectories are allowed (see Fig. 12).

C. The evolution of the azimuth angle (the φ-motion)
A feasible approach to calculate the φ-motion, is to divide the integral equation (28) into θ-dependent and r-dependent parts, in a way that we get with in which, the minimum value of the θ coordinate has been assumed to coincide with the initial point r i , which can be set as either of the turning points. According to the fact that the general form of the equation of motion has been considered, this assumption does not affect the final analytical results for the φ-motion. Direct integration of the integral (73a) results in the following cases: • Forη > 0 we get in which, κ 0 has been given in Eq. (59a), and (with j = 1, 2) where, ℘ (υ) ≡ d dυ ℘(υ; g 2 , g 3 ), which here, is given in terms of the corresponding Weierstraß invariants g 2 and g 3 , as in Eqs. (61). Furthermore, in which, ß(x) ≡ ℘ −1 (x; g 2 , g 3 ) is the inverse Weierstraßian ℘ function, z max and ψ 0 are given in Eqs. (57a) and (59b), The function in Eq. (75), also depends on the Weierstraßian Zeta and Sigma functions (see appendix B). Note that, for the sake of simplicity in the demonstration of the trajectories, we will set φ(θ min ) = 0 as the initial condition. The complete γ-dependent expression for Φ θ (γ) is then obtained by substituting θ → θ(γ) in the above relations.
The r-dependent integral (73b) provides in which given that where the coefficient α 1 and the Weierstraß invariants, are the same as those given in Eqs. (47), considering r D → r i . Furthermoreα Now, having in hand the analytical expressions for all the spatial coordinates, we are able to simulate the possible orbits, based on the simultaneous evolution of these coordinates. Respecting the allowed values of ω 0 , in Fig. 13, the orbits that we have previously illustrated in Figs. 8, 10, and 12, are demonstrated in the three-dimensional form, by applying the following Cartesian correspondents of the Boyer-Lindquist coordinates [60] x(γ) = r 2 (γ) + a 2 sin θ(γ) cos φ(γ), (89a) In presenting the figures, we have also considered the case of the vacuum Kerr spacetime, for which, the characteristic functions f r (r) and f θ (θ) in Eq. (16), vanish identically. In fact, in the three-dimensional treatment of the trajectories, a fixed positive Carter's constant Q, does not allow for the construction of the casesη = 0 andη < 0 in the vacuum Kerr spacetime. Hence, the possible comparison between the light propagation in the plasmic and vacuum Kerr spacetimes can only be done in the context ofη > 0 which corresponds to the first six diagrams of Fig. 13. As indicated in these diagrams, the sensible differences between these media, are seen in the DFKs. While the light ray trajectories in plasma (black curves) occupy a wider spatial range around the black hole, those in the vacuum (blue curves) are more close to the event horizon. Passing the plasma, will therefore, change the amount of light deflection and this can be inspected through the process of gravitational lensing (see below). The other types of trajectories do not show any sensible differences. Hence, in what follows we continue with the equatorial lens equation for the black hole.

Gravitational lensing
Let us, for now, confine ourselves to the equatorial plan (with θ = π/2), on which we have Q = f θ . Therefore, by means of Eqs. (26) and (28), the differential equation that governs the gravitational lensing is according to which, the lens equation is written asθ whereθ is the deflection angle and r D is that in Eq. (40). This provideŝ in which, F ± are given in Eq. (86), and Υ ± are the same as those in Eq. (87d) with r i → r D . Applying the data given in Fig. 13, one obtainsθ = 47.137 • and ϑ = 110.869 • , respectively, for ω 0 = 0.755 and ω 0 = 0.70. This is while for a vacuum background (i.e. f r = f θ = 0), these values change toθ = 27.276 • andθ = 52.997 • , for the same initial frequencies.
• Forη < 0: The r-dependent integral (95b) provides the solution t r (γ) = τ 0 5 j=1 τ j T j r(γ) − T j (r i ) , effective potential does not offer any planetary orbits, so that the light rays, beside being able to form photon rings on the UCO, can only escape from or be captured by the black hole. These kinds of trajectories, as demonstrated within the text, are bounded to cones of definite vertex angles. Based on the specific frequency (energy) of the light rays, ω 0 , the intensity of the deflection in the escaping trajectories can vary, and the rays may travel on more fast-changing hyperbolic curves as ω 0 approaches its critical value, ω U . Note that, since we express the plasma frequency ω p in terms of the black hole's physical characteristics in Eq. (32), the impacts of plasma are then included, indirectly, through the variations of the effective potential and the consequent types of orbit. In the presented analytical solutions, this was done by introducing the effective impact parameterη and its dependent constants. Although the temporal evolution of the coordinates have appeared to be of rather complicated mathematical expressions, they however, can establish the basement of further studies, specified to certain models of plasma that are in connection with black holes' parameters. Such studies can offer astrophysical applications and are left to our future works.
Furthermore, in Eq. (B2c), it is defined that There is also another integral relation that we have exploited in this paper: