Photon motion in Kerr-de Sitter spacetimes

We study general non-equatorial motion of photons in the Kerr-de Sitter black hole and naked singularity spacetimes. The motion is governed by the impact parameter $X$ related to axial symmetry of the spacetime, and impact parameter $q$ related to its hidden symmetry. Appropriate 'effective potentials' governing the latitudinal and radial motion are introduced and their behaviour is examined by 'Chinese boxes' technique giving regions allowed for the motion in terms of the impact parameters. Limits on the impact parameter $X$ and $q$ are established in dependence on the spacetime parameters $M, \Lambda, a$. The motion can be of orbital type (crossing the equatorial plane, $q>0$) and vortical type (tied above or bellow the equatorial plane, $q<0$). It is shown that for negative values of $q$ the reality conditions imposed on the latitudinal motion yield stronger constraints on the parameter $X$ than that following from the reality condition of motion in the radial direction, which, e. g., exclude possibility of existence of vortical motion of constant radius. Several other consequences are determined and a classification of Kerr-de Sitter spacetimes reflecting the behaviour of the effective potentials is given.


I. INTRODUCTION
In the framework of inflationary paradigm [45], recent cosmological observations indicate that a very small relict vacuum energy (equivalently, repulsive cosmological constant Λ > 0), or, generally, a dark energy demonstrating repulsive gravitational effect, has to be introduced to explain dynamics of the recent Universe [3,7,16,40,41,49,92]. These conclusions are supported strongly by the observations of distant Ia-type supernova explosions indicating that starting at the cosmological redshift z ≈ 1 expansion of the Universe is accelerated [55]. The total energy density of the Universe is very close to the critical energy density ρ crit corresponding to almost flat universe predicted by the inflationary scenario [62], and the dark energy represents about 70% of the energy content of the observable universe [15,63]. These conclusions are confirmed by recent measurements of cosmic microwave background anisotropies by the space satellite observatory PLANCK [1,52].
The dark energy equation of state is very close to those corresponding to the vacuum energy [15]. Therefore, it is relevant to study the astrophysical consequences of the effect of the observed cosmological constant implied by the cosmological tests to be Λ ≈ 1.3 × 10 −56 cm −2 , and the related vacuum energy ρ vac ∼ 10 −29 g/cm 3 , close to the critical density of the universe. The repulsive cosmological constant changes significantly the asymptotic structure of black-hole, naked singularity, or any compact-body backgrounds as such backgrounds become asymptotically de Sitter spacetimes, and an event horizon (cosmological horizon) always exists, behind which * Electronic address: daniel.charbulak@fpf.slu.cz † Electronic address: zdenek.stuchlik@fpf.slu.cz the geometry is dynamic.
The role of the cosmological constant can be significant for both the geometrically thin Keplerian accretion discs [46,61,69,72,86] and the toroidal accretion discs [6,42,53,54,60,87,89] orbiting supermassive black holes (Kerr superspinars) in the central parts of giant galaxies. Both high-frequency quasiperiodic oscillations and jets originating at the accretion discs can be reflected by current carrying string loops in SdS and KdS spacetimes [28,33,77,78,93]. In the spherically symmetric spacetimes, the Keplerian and toroidal disc structures can be precisely described the Pseudo-Newtonian potential of Paczynski type [79,88] that appears to be useful also in studies of motion of interacting arXiv:1702.07850v3 [gr-qc] 19 Feb 2018 galaxies [56,81,82] demonstrating relation of the gravitationally bound galactic systems to the so called static radius of the SdS or KdS spacetimes [4,5,23,24,66,67]. This idea has been confirmed by the recent study of general relativistic static polytropic spheres in spacetimes with the repulsive cosmological constant [75,85].
The present paper is devoted to detailed study of properties of the photon motion in the KdS black hole and naked singularity spacetimes. We concentrate attention to the behavior of the effective potentials determining the regions allowed for the photon motion. Such a study is necessary for full understanding of the optical phenomena occuring in the black hole or naked singularity spacetimes with the repulsive cosmological constant. We generalize the previous work concentrated on the properties of the photon motion in the equatorial plane [73], discussing properties of the effective potential of the latitudinal motion in terms of the motion constant related to the equatorial plane, and then continuing by study of the effective potential of the radial motion. We concentrate our study on the spherical photon orbits representing a natural generalization of the photon circular geodesics that enables a natural classification of the KdS spacetimes according to the properties of the null geodesics representing the photon motion.

II. KERR-DE SITTER SPACETIME AND CARTER'S EQUATIONS OF GEODESIC MOTION
A. Kerr-de Sitter geometry The line element describing the Kerr-de Sitter geometry is in the standard Boyer-Lindquist coordinates, using geometric system of units (c = G = 1), given by ds 2 = − ∆ r I 2 ρ 2 dt − a sin 2 θ dφ 2 + ∆ θ sin 2 θ I 2 ρ 2 a dt − r 2 + a 2 dφ where Here, as usual, we denoted by M the mass of the central gravitating body, by a its specific angular momentum (a = J/M ) and by Λ the cosmological constant. In order to simplify the discussion of the following equations, it is convenient to introduce a new cosmological parameter y = 1 3 ΛM 2 , and use dimensionless quantities, redefining them such that s/M → s, t/M → t, r/M → r, a/M → a, which is equivalent to putting M = 1. The above expressions then read ∆ r = 1 − yr 2 r 2 + a 2 − 2r, (6) ∆ θ = 1 + a 2 y cos 2 θ, (7) I = 1 + a 2 y, with equation (5) being left unchanged. The physical singularity is located, as in the Kerr geometry, at the ring r = 0, θ = π/2. The black hole horizons are determined by the condition ∆ r = 0 (9) and their loci can be determined by the relation y = y h (r; a 2 ) ≡ r 2 − 2r + a 2 r 2 (r 2 + a 2 ) .
The zeros of y h (r; a 2 ), determining the loci of black hole horizons in pure Kerr spacetimes, are given by the relation the loci of its extrema are given by the functions where the function a 2 ex(h)− (r) < 0 in its whole definition range, hence is irrelevant. The functions y h (r; a 2 ), a 2 z(h) (r) and a 2 ex(h)± (r) will be needed in the section devoted to the discussion of the radial motion. Three event horizons, two black hole r − , r + , and the cosmological horizon r c , (r − < r + < r c ) exist for y min(h) (a 2 ) < y < y max(h) (a 2 ), where the limits y min/max(h) (a 2 ) correspond to local minimum or local maximum of the function y h (r; a 2 ), respectively, for given rotational parameter a. For 0 < y < y min(h) (a 2 ) or y > y max(h) (a 2 ) Kerr-de Sitter naked singularity spacetimes exist. The limit case y = y min(h) (a 2 ) corresponds to an extreme black hole spacetime, when the two black hole horizons coalesce. If y = y max(h) (a 2 ), the outer black hole and cosmological horizon merge. There exists critical value of the rotational parameter a 2 crit = 1.21202, for which the two local extrema of the function y h (r; a 2 ) coalesce in an inflection point at r crit = 1.61603 with the critical value y crit = 0.0592. Thus, for a 2 > a 2 crit only Kerr-de Sitter naked singularity can exist for any y > 0.
Properties of the event horizons for the more general case of the Kerr-Newman-de Sitter spacetimes can be found in [73].

B. Carter's equations of geodesic motion
The motion of test particles and photons following its geodesics in the Kerr-de Sitter spacetime is described by the well known Carter equations [17] and Here E and Φ are the constants of motion connected respectively with the time and axial symmetry of the Kerrde Sitter geometry, and K is the fourth Carter constant of motion connected with the hidden symmetry of the Kerr-de Sitter geometry. Another constant of motion is the rest mass m (energy) of the test particle; for photons m = 0. Recall that E and Φ cannot be interpreted as energy and axial component of the angular momentum of the test particle at infinity, because, due to the presence of the cosmological Λ term, the geometry is not asymptotically flat, but de Sitter [72].
Detailed discussion of the equatorial motion of photons in the Kerr-Newman-de Sitter spacetimes has been published in [73]. Circular motion of test particles in the Kerr-de Sitter spacetimes has been presented in [86].
Here we restrict our attention on the general motion of photons in the Kerr-de Sitter spacetimes.
In fact, the motion of photons is independent of the constant of motion E and depends only on the ratio Φ/E (E = 0), usually referred to as the impact parameter , and on the parameter K/E 2 . For our general discussion it is convenient to use Q = K−I 2 (Φ−aE) 2 that vanishes for the equatorial motion. For our purposes it is, however, following the paper [73], convenient to introduce a new constants of motion X ≡ − a. Further the constant q ≡ Q/I 2 E 2 will be applied. Then the relations (17) and (18) simplify to the form R(r; X, q, y, a) ≡ I 2 E 2 r 2 − aX 2 − ∆ r X 2 + q .
(20) Following the work [73] we study the general photon motion in terms of the parameter X. However, since we consider the non-equatorial motion here, it is also necessary to find out the restrictions to be imposed on parameter X that follows from the reality conditions of the latitudinal motion. The latitudinal motion in the Kerr-de Sitter spacetimes has been investigated yet [66]; however, the discussion has been related to the motion constant K.
Here we give the discussion of the effective potential of the latitudinal motion related to the motion constant Q, as it is convenient for the purposes of our study.

III. LATITUDINAL MOTION
Because it is more convenient to work with algebraic functions instead of trigonometric ones, we introduce a new variable This implies replacing the equation (13) by where M (m; X, q, y, a) ≡ Note that dm/dλ = 0 does not necessarily imply dθ/dλ = 0, since it can mean just transit through the equatorial plane or polar axis. Therefore, in some cases, in order to avoid any doubts, we rather discuss the behaviour of the function (19). The reality condition M (m; a, y, X, q) ≥ 0 can be expressed by the relations in regions where ∆ m − a 2 y > 0, i. e., equivalently, is the solution of the equation and by the relations in regions where ∆ m − a 2 y < 0, i. e., m < m d , which requires y > 1/a 2 . The functions X θ ± (m; q, y, a), regarded as 'effective potentials' governing the latitudinal motion, are defined by The functions X θ ± (m; q, y, a) thus determine the regions allowed for the latitudinal motion, conditions X = X θ ± (m; q, y, a) give the turning points. In order to understand the behaviour of the functions X θ ± (m; q, y, a), it is necessary to find the reality regions, and loci of its local extrema and divergencies. Following [73], we shall perform this analysis using the well known procedure called 'Chinese boxes technique' and adopting labelling of the appropriate characteristic functions in similar way. The parameters are of various significance -q is a constant of motion, whereas a, y govern the geometry. The natural choice is therefore to give the properties of the potentials X θ ± (m; q, y, a) by family of functions q(m; y, a), and properties of these functions by another families of functions of variable m with parameters lowered by one, with spacetime parameters excluded at last.
In the following analysis the relevant range of variable m is, of course, 0 ≤ m ≤ 1, but somewhere, in order to better understand the behaviour of the characteristic functions, we formally permit m ∈ R. First we shall determine the reality region of X ± (m). It is given by where Of course, this function also determines the common points of the potentials X θ − (m; q, y, a) and X θ + (m; q, y, a), which values are then Of particular importance, if defined, is the value X (±) (1; y, a) = −a (see bellow).
The divergency points of the functions q θ r (m; y, a 2 ), X (±) (m; y, a) and X θ − (m; q, y, a) are determined by Both the functions X θ ± (m; q, y, a) can diverge, if well defined, for m = 0, another divergencies are given by the function y θ d (m; a 2 ) for the potential X θ − (m; q, y, a), but there are no other divergencies for the potential X θ + (m; q, y, a), as can be seen, if we rewrite the definition (29) in an alternative form X θ ± (m; q, y, a) = (34) The function y θ d (m; a 2 ) → ∞ for m → 1 from the left. There are no local extrema of this function and for 0 ≤ m < 1 it is increasing. For m = 0 we get y θ d (0; a 2 ) = 1/a 2 . The point m d given by the definition (25) determines the loci where the functions q θ r (m; y, a 2 ), X (±) (m; y, a) and X θ − (m; q, y, a) diverge; it occurs at relevant interval (0; 1) for y > 1/a 2 and m d → 1 for a 2 y → ∞. In such case, q θ r (m; y, a 2 ) → +∞ (−∞) for m → m d from the left (right).
From the equality one can see that the function q θ r (m; y, a 2 ) has no local extrema and is decreasing for y < 1/a 2 , or piecewise increasing with the discontinuity point m d for y > 1/a 2 , i. e., q θ r (m; y, a 2 ) → +∞ (−∞) for m → m d from the left (right). It always holds q θ r (m = 0; y, a 2 ) = 0 and q θ r (m = 1; y, a 2 ) = −a 2 .
Based on the conditions (29), (30) and the above characteristic functions, we can complete setting the definition range of the potentials X θ ± (m; q, a, y), which we leave to the end of this section. Now we shall determine the loci of local extrema of the effective potentials X θ ± (m; q, y, a). They can be derived from the condition ∂X θ ± /∂m = 0, which implies the equation (a 2 m 2 +q)(∆ m −a 2 y)[a 2 m 2 I 2 +q(1−a 2 y +2a 2 my) 2 ] = 0.
(36) It can be verified that the function X θ + (m; q, y, a) has local extrema given by the relation A discussion of this function is trivial, so we only note that it is independent of the cosmological parameter y and renders the loci of extrema only for −a 2 ≤ q ≤ 0, while, as we shall see below, they can exist even for q < −a 2 . The character of these extrema reveals inserting this expression into the second derivative, which yields clearly they must be maxima. From the equation (36) we find that another extrema of the potentials X θ ± (m; q, y, a) are determined by the condition The divergencies of the functions q θ ex(±) (m; y, a 2 ) are determined by the relation .
From the properties of the function y θ d(ex±) (m; a 2 ) we deduce that the function q θ ex(±) (m; y, a 2 ) can diverge only if y > 1/a 2 , at In the following we shall decide about the monotony and possible existence of local extrema of the function q ex(±) . From it is clear that there are no local extrema of q θ ex(±) (m; y, a 2 ) in the interval m ∈ (0; 1). The derivative changes its sign at the divergent point m d(ex) , which reflects the behaviour given by (41). For y < 1/a 2 , there is ∂q θ ex(±) /∂m < 0 for m ∈ 0; 1 , that is, q θ ex(±) (m; y, a 2 ) is decreasing. In the limit case y = 1/a 2 , we get q θ ex(±) (m; y = 1/a 2 , a 2 ) = −a 2 = q θ r (m; y = 1/a 2 , a 2 ).
In the next step we shall characterize the extrema given by q θ ex(±) (m; y, a 2 ). First we find that If we now require ∂X θ + /∂m(m; q = q θ ex(±) , y, a) = 0 somewhere at 0 < m < 1, we obtain a condition m > m d(ex) for y > 1/a 2 , which ensures Considering the previous results, we can conclude that the functions q θ ex(±) (m; y, a 2 ) determine local maxima of the potential X θ + (m; q, y, a) for q < −a 2 that occur on this curve in the case y > 1/a 2 in the interval m ∈ (m d(ex) ; 1). Proceeding the same way with the function X θ − (m; q, y, a), we first find that the equation ∂X θ − /∂m(m; q = q θ ex(±) , y, a) = 0 has always solution for some m ∈ (0; 1) in the case y < 1/a 2 , but for y ≥ 1/a 2 this solution must fulfil m < m d(ex) . Substituting q = q θ ex(±) into the second derivative of X θ − (m; q, y, a) yields the same expression as that in (45), but now with the above conditions we have indicating local minima. Therefore, the function q θ ex(±) gives local minima of X θ − (m; q, y, a) for y < 1/a 2 at the whole interval m ∈ (0; 1), and for y > 1/a 2 at m ∈ (0; m d(ex) ). In the special case y = 1/a 2 , the function q θ ex(±) reduces to the form q θ ex(±) (m; y = 1/a 2 , a 2 ) = −a 2 ; we can easily convince ourself that the function X θ − has no local extremum in such case, and the extrema of X θ + are given by the function q θ ex+ (m; y, a). The conditions (24), (27) ensuring the allowance of the latitudinal motion must be complemented by case when the functions X θ ± are not defined. Their definition range is given by relations (29), (30), but it can be shown that the violation of the latter one imply M (m; X, q, y, a) > 0. In such case, the latitudinal motion is allowed for any impact parameter X (see the details in the discussion bellow). All characteristic functions are depicted in Fig.  1. and the graphs of the potentials in Fig. 2 for selected representative values of parameters.
Now we are able to discuss the behaviour of the potentials X θ ± (m; q, y, a) for various representative values of its parameters. The intersections of a line X = const. with the curves X θ ± (m; q, y, a) represents the turning points in variable m. From the knowledge of these functions we can thus get qualitative insight into the character of the latitudinal motion. This entitles us to following classification of Kerr-de Sitter spacetimes and brief description of the latitudinal motion. The basic division apparently consists of cases y < 1/a 2 , y = 1/a 2 , y > 1/a 2 : the definition range of the potentials is an empty set; the latitudinal motion is not possible; the potentials X θ ± (m; q, y, a) are defined only for m = 1, where photons with such values of parameters are the special case of the so called PNC photons 'radially' moving along the spin axis [10]; • −a 2 < q < 0 (Fig. 2a, 2b) both the potentials are defined for m ∈ m l ; 1 , where the lower limit is the solution of the equation the limits of the interval are the common points of the potentials, where is given by the equation (39), the loci m max(+) of maximum X θ max(+) is determined by relation (37); if X takes one of these extremal values, then the trajectory of such photon lies entirely on cones θ = arccos √ m ex , there are two solutions m 1 < m 2 of each of the two equations X = X θ ± (m; q, y, a), implying that photon executes so called vortical motion, which is restricted between two pairs of cones, symmetrically placed relative to equatorial plane: FIG. 2: Graphs of potentials X θ + (m; q, y, a) (full curve) and X θ − (m; q, y, a) (dashed curve) depicted for given typical values a 2 , y corresponding successively to cases y < 1/a 2 (top row), y = 1/a 2 (middle row), y > 1/a 2 (bottom row), and with q representing significant cases q < 0, q = 0, q > 0. Shading demarcates the region where the latitudinal motion is forbidden.
in the special case X = −a one of the turning points is m 2 = 1, which represents transit through the spin axis; such photon therefore oscillate above one of the poles in cone which is delimited by the angle θ = arccos √ m 1 ; from the preceding discussion it follows that we can expect that the case X = −a represents a change in azimuthal direction with respect to some privileged family of observers; • q = 0 ( Fig. 2c) the expression in the definition (29) can be reduced to which validity can be enlarged, without any repercussion on the correctness of the analysis, even for m = 0; the definition range of the potentials is thus 0; 1 ; from the equality W (θ = π/2; X, q, y, a) = q it follows that at least in the equatorial plane the (radial) motion always exists for q = 0, where it can be both stable or unstable (see bellow); for q > 0 the equatorial plane is crossed, for q < 0 it can not be reached; there are no extrema of the potentials -X θ + (m; q, y, a) is decreasing, X θ − (m; q, y, a) is increasing; the permissible values of X for which dθ/dλ > 0 are still confined to an interval with limits X θ max(+) = X θ + (m = 0; q = 0, y, a) = 0, (50) where X θ min(−) < X θ max(+) ; if X ≤ X θ min(−) or X ≥ X θ max(+) then the requirement W (θ) ≥ 0 is fulfilled only if θ = π/2, and in such case dθ/dλ = 0, thus the motion is stably confined to the equatorial plane; for X θ min(−) < X < X θ max(+) photon initially released in the direction off the equatorial plane is once reflected at θ = arccos √ m 0 or θ = π−arccos √ m 0 respectively, where m 0 denotes the only solution of X = X θ ± (m; q, y, a); another point where dθ/dλ = 0 is now in the equatorial plane, however the equality d 2 θ/dλ 2 = 0 implies halting in the latitudinal direction; the function W (θ) has at θ = π/2 local minimum, which indicates, as follows from perturbation analysis, instability in the equatorial plane; if specially X = −a then m 0 = 1, thus photon initially directed off the equatorial plane crosses the spin axis and finally is captured in the equatorial plane; • q > 0 (Fig. 2d) the potentials are defined for m ∈ (0, 1 ; they are monotonous in the same manner as in the case q = 0, but X θ + (m; q, a, y) → +∞ and X θ − (m; q, y, a) → −∞ as m → 0; from the behaviour of the potentials it follows that for X = −a photon is forced to oscillate in θ-direction through the equatorial plane between two cones governed by arccos with m 0 of the same meaning as above; case X = −a represents the motion above both poles; the foregoing conclusion is a reason to have a suspicion that cases X < −a and X > −a differ in the azimuthal direction relative to some family of stationary observers, it corresponds to > 0 and < 0; 2. Case y = 1/a 2 the potentials simplify into the form • q < −a 2 the potentials are not defined, thus the latitudinal motion is not allowed; the curves X = X θ ± (m; q = −a 2 , y = 1/a 2 , a) coalesce, since for X ≤ −a there is one solution of the equation X = X θ (±) (m; a), which gives m = m (±) ≡ −a/X; this corresponds to PNC photons moving along for X → −∞ the cones approach the equatorial plane; if specially X = −a the cones degenerate to spin axis, therefore, such PNC photons move along the spin axis; for X > −a there is no motion allowed; • −a 2 < q < 0 (Fig. 2e, 2f) the potentials are both defined for m ∈ (0; 1 ; there is one local maximum X θ max(+) given by (37) of the function X θ + (m; q, y, a) and no ex- , the vortical motion exists; for X = −a the 'inner' cones coalesce with the spin axis, thus the vortical motion involves crossing the poles; for X = X θ max(+) both the 'inner' and 'outer' cones coalesce, giving thus rise to PNC photons; if X > X θ max(+) , no motion is allowed; • q = 0 ( Fig. 2g) the same discussion holds as in the case y < 1/a 2 , except that the motion exists for X arbitrarily small; • q > 0 ( Fig. 2h) the same conclusions holds as in the case y < 1/a 2 , ; 3. Case y > 1/a 2 • q < −a 2 (Fig. 2i) the definition range of both potentials is an interval (0; m u (see the purple curve in Fig. 1d), where the upper limit m u < 1 is given as m l in the previous case by (46); there is X θ + (m; q, y, a) → −∞ and X θ − (m; q, y, a) → +∞ as m → 0, moreover, X θ − (m; q, y, a) now diverge at m = m d , which is the solution of (33), and X θ − (m; q, y, a) → +∞ (−∞) as m → m d from the left (right); there are thus two regions of permissible values X in the (m, X)-plane for which the motion can exist; the lower one bounded by the graph of X θ + and the lower branch of X θ − , which at m = m u join into continuous curve, and the upper region given by the upper branch of X θ − ; the motion is therefore allowed for X ≤ X θ max(+) < −a or X ≥ X θ min(−) > 0, where the loci of local extrema X θ max(+) , X θ min(−) are given by (39) (see the blue curve in Fig. 1d); , the inner cones delimiting the vortical motion are the narrowest; for X → −∞ or X → +∞ the outer cones given by angles approach the equatorial plane since m 1 → 0; for the inner cones there is • q = −a 2 (Fig. 2j) there is no local extremum of the function X θ + (m; q, y, a), which is now increasing; it holds m u = 1, X θ − (m u ) = X θ + (m u ) = X θ +(max) = −a, hence for X = −a both the inner and outer cones coalesce with the spin axis, which again corresponds to 'axial' PNC photon; another PNC photons exist for X = X θ min(−) > 0; there are no other qualitative differences from the case q < −a 2 ; • −a 2 < q < 0 (Fig. 2k) the definition range is an interval (0; 1 and the divergencies of the potentials are the same as above; the function X θ + (m; q, y, a) has now local maximum X θ max(+) , −a < X θ max(+) < 0, X θ max(+) → 0 for q → 0, determined by equation (37); case X = −a now corresponds to vortical motion above the poles -the inner cones have coalesced with the spin axis, the outer ones stay open; the vortical motion exists as in the previous cases and above that for −a < X < X θ max(+) ; • q = 0 (Fig. 2l) the definition (48) holds, the functions X θ +(−) (m; q, y, a) are defined at 0, 1 ( 0, 1 \ {m d }); the values for m = 0 are given by (50), but now X θ max(+) < X θ min(−) ; the potential X θ + (m; y, a) is decreasing in its whole definition range, X θ − (m; y, a) is piecewise increasing because of the divergent point m d ; if X ≤ X θ max(+) = 0 or X ≥ X θ min(−) then the same conclusions can be made as in the case y < 1/a 2 for X θ min(−) ≤ X ≤ X θ max(+) ; for X θ max(+) < X < X θ min(−) it holds W (θ = π/2; X, q = 0, y, a) = 0 again, otherwise W (θ; X, q = 0, y, a) < 0, therefore photons can radially move in the equatorial plane; where m l is given by (46) with the difference that now X θ (±) (m l ) > 0; the graphs of both functions now form a single open curve, which intersects a line X = const. at a single point; in the interval m ∈ 0; m l the latitudinal motion is allowed for arbitrarily large or small value of the motion constant X; for arbitrary X = −a there exists oscillatory motion through the equatorial plane as described in the case y < 1/a 2 , q > 0; if X = X (±) (m l ) the boundary cones are closest to equatorial plane, they are given by angles the case X = −a corresponds to orbits above both poles crossing also the equatorial plane; there is no vortical motion or PNC photons; We finish this section with setting the allowed region in the (X, q)-plane delimiting such combinations of the motion constants, for which the latitudinal motion is possible, in dependence on the spacetime parameters a, y.
From the requirement that the function M (m; a, y, X, q) defined in the relation (22) has to be non-negative somewhere in the interval m ∈ 0; 1 , one can derive that the allowed region of the (X − q)plane is determined by the condition where q min (X, y, a) is defined using functions and as follows (see Fig. 3 and Fig. 4): • Case y < 1/a 2 (Fig. 3a) • Case y = 1/a 2 (Fig. 3b) • Case y > 1/a 2 (Fig. 3c) The case y < 1/a 2 qualitatively corresponds to both black hole and naked singularity spacetimes, the other two cases y = 1/a 2 and y > 1/a 2 describe the naked singularity spacetimes (see Fig. 7 in the next section). The q = const.-slices of the function q min (X, y, a) give for q < 0 extremal values X min(−) , X max(+) of the potentials X θ ± (m; q, y, a) discussed in the text.

IV. RADIAL MOTION
From the equation (14) it is clear that the radial motion can exist if R(r) ≥ 0, where the equality gives the turning points of the radial motion. This condition can be rewritten in terms of an 'effective potential' X ± in the form if a 2 − ∆ r > 0 (and X − < X + ), The allowed region in the motion constant plane (X − q) depicted for some representative values of spacetime parameters a, y, successively corresponding to cases y < 1/a 2 , y = 1/a 2 , y > 1/a 2 . An intersections of a line q = const. < 0 with the border curves q1(X), q2(X; y, a) determines the extremal values X min(−) , X max(+) of the potentials X θ ± (m; q, y, a) introduced in the discussion above. or where We start the analysis by determining the reality region of the effective potential X ± . From the expression (61) it follows that this function is well defined for where we have introduced the reality function . (63) There are thus two different types of the boundary points of the definition range of X ± (r; q, a, y). The points of the first type lie stably for given spacetime parameters on the borders of the static regions determined by the relation 9, i. e. at the event horizons (r = r h ). At these horizons, if they exist, for arbitrary parameter q, the functions X ± have common values (c. f. [73]). The points of the second type, which are also common points of X ± , depend on the value of parameter q and are given by the equality q = q r (r; y, a 2 ). If we denote them r = r q then it holds The divergencies of the function q r (r; y, a 2 ), which are incident with divergent points of X + (r; q, y, a), are located at radii where which one can express by the relation .
The function X − (r; q, y, a) can not diverge at radii r d given by (66), since using an alternative expression X ± (r; q, y, a) ≡ r 4 − q∆ r it can be shown that it has finite value Another point where the functions X ± (r; q, y, a) diverge is r = 0, with X ± (r; q, y, a) → ±∞ as r → 0 for q > 0, but for q = 0 it holds X ± (r; q = 0, y, a) → 0. The character of the function y d (r; a 2 ) has been discussed in [70,73], therefore we briefly repeat that the only zero of y d (r; a 2 ) is at r = 2, the extrema, which for a 2 > 0 must be maxima, yields the relation .
The only zero of q r (r; y, a 2 ) is at r = 0. For r → ∞ it holds q r (r; y, a) → −1/y. Its extrema are determined by The divergency of y ex(r) (r; a 2 ) is at r = 0 and y ex(r) (r; a 2 ) → −∞ for r → 0. For r → ∞ it approaches the line 1/a 2 from bellow. The zero is at r = 3 and its extrema do not exist, the function is purely increasing. Now we shall specify the local extrema of the effective potential, which determine the radii of spherical photon orbits. They are given by the condition ∂X ± /∂r = 0, which implies or This can be rewritten in terms of parameter q as and q = q ex (r; y, a 2 ) (75) Note that the function q ex1 (r; a 2 ) is independent of the cosmological parameter. Both the functions q ex1 (r; a 2 ), q ex (r; y, a 2 ) have common points determined by and y = 1 i. e., they are located at event horizons and so called static radius r s = 1/ 3 √ y, where the gravitational attraction is just compensated by cosmological repulsion [66,72]. The function q ex1 (r; a 2 ) is negative valued and hence, as we shall see bellow, the extrema of the potentials X ± determined by this function lie in regions forbidden by conditions for the reality of latitudinal motion. The divergencies of q ex (r; y, a 2 ) are determined by the relation its asymptotic behaviour is given by q ex (r; y, a 2 ) → −(I/2ay) 2 as r → ∞.
The function y d(ex) (r; a 2 ) diverges for r = 0 and y d(ex) (r; a 2 ) → −∞ as r → 0. For r → ∞ it holds y d(ex) (r; a 2 ) → 0. The zero of this function is at r = 1 and its local extrema are determined by the relation where the label 'max' indicates that at relevant range r ≥ 3/2, these extrema must be maxima. The zero point of the function q ex (r; y, a 2 ) is at r = 0, another zeros determine the loci of the circular equatorial photon orbits. They are given by the relation Since the function y z(ex)− (r; a 2 ) < 0 for r > 0, it is irrelevant in our discussion. The function y z(ex)+ (r; a 2 ) is real valued for all r > 0 and it diverges at r = 0, with y z(ex)+ (r; a 2 ) → ∞ as r → 0. For r → ∞, we find y z(ex)+ (r; a 2 ) → −1/a 2 . Its zeros represent the equatorial circular photon orbits in the Kerr spacetimes, being determined by the relation [65] The extrema of the function y z(ex)+ (r; a 2 ) are determined by the equation hence the loci of extrema of the functions y z(ex)+ (r; a 2 ) and y h (r; a 2 ) coalesce. The function a 2 ex(z(ex)+)− (r) should be excluded from further analysis since for r > 0 there is a 2 ex(z(ex)+)− (r) < 0. It remains to determine loci of the local extrema of the function q ex (r; y, a 2 ). Proceeding the usual way we find that their occurrence is governed by the relations and y = y ex(ex)± (r; a 2 ) ≡ Using the relation (83), one can show that the extrema of both functions q ex (r; y, a 2 ), q r (r; y, a 2 ) coalesce. At this point let us add that another common points of functions q ex (r; y, a 2 ), q r (r; y, a 2 ), as well as of the functions q ex1 (r; y, a 2 ), q r (r; y, a 2 ), are also given by i. e. they are located at the event horizons. The reality conditions of the functions y ex(ex)± (r; a 2 ) read if 0 < r ≤r, and a 2 ≤ a 2 r(ex(ex±)) (r) or a 2 r(ex(ex±))+ (r) ≤ a 2 , (87) and The marginal radiusr has valuê and it holds a 2 r(ex(ex)±)+ (r) = a 2 r(ex(ex)±) (r) = a 2 crit = 1.21202, (91) where a 2 crit corresponds to local maximum of the function a 2 ex(h)+ (r) (see e. g. [73] for details). The functions y ex(ex)± (r; a 2 ) have the divergency point at r = 0 and y ex(ex)± (r; a 2 ) → ±∞ for r → 0. For r → ∞ we find that y ex(ex)+ (r; a 2 ) → ∞ and y ex(ex)− (r; a 2 ) → 0 from above. The zero point of y ex(ex) (r; a 2 ) is at r = 3 and the function is increasing for all r > 0. Zeros of the functions y ex(ex)± (r; a 2 ) are given by The condition for stationary points ∂y ex(ex)± (r; a 2 )/∂r = 0 leads to a 4 + a 2 r(2r − 1) + r 3 (r − 3) = 0, which can be solved with respect to a 2 with the same result as given by (12). However, substitution into the second derivative concurrently with the requirement y ex(ex)± (r; a 2 ) > 0 implies ∂ 2 y ex(ex)± (r; a 2 = a 2 ex(h)+ (r)) = 0, therefore the function determines the loci of the inflex points of the functions y ex(ex)± (r; a 2 ). If we compare the asymptotic behaviour of all characteristic functions y(r; a 2 ), we find that following inequality is satisfied: In Fig.5 we present all the characteristic functions related to spin parameter governing the effective potential on the lowest level: • a 2 max(d(ex)) (r) • a 2 z(z(ex)+) (r) • a 2 r(ex(ex)±)+ (r) • a 2 r(ex(ex)±) (r) • a 2 z(ex(ex)±) (r).
These functions determine the behaviour of the characteristic functions related to the cosmological parameter, characterizing the functions q(r; y, a 2 ) and then effective potentials on the higher level: • y h (r; a 2 ) • y d (r; a 2 ) • y ex(r) (r; a 2 ) = y ex(ex) (r; a 2 ) • y d(ex) (r; a 2 ) • y z(ex)+ (r; a 2 ) • y ex(ex)± (r; a 2 ). From the significance of the individual characteristic functions a 2 (r) depicted in Fig. 5, one can infer that there are just two values of a 2 being of particular importance and leading to qualitatively different behaviour of the functions y(r; a 2 ) : a 2 = 1 -the common local maximum of the functions a 2 z(h) (r) and a 2 z(z(ex)) (r) at r = 1, which coincides with the inflection point of the function a 2 z(ex(ex)±) (r; a 2 ) and with the intersection with the curve a 2 ex(h)+ (r) a 2 = a 2 crit = 1.21202 -the local maximum a 2 ex(h)+ (r) which is the intersection of the curves a 2 r(ex(ex)±)+ (r), a 2 r(ex(ex)±) (r) and a 2 max(d(ex)) (r).
The graphs of characteristic functions y(r; a 2 ) depicted for some values of spin parameter a representing cases 0 < a 2 < 1, 1 < a 2 < a 2 crit and a 2 crit < a 2 are presented in Fig. 6.
In general, behaviour of the characteristic functions q r (r; y, a 2 ) and q ex (r; y, a 2 ) will be qualitatively different, if for the parameter a being fixed we take the yvalues from different intervals, which are limited by intersections and/or extrema of the characteristic functions y(r; a 2 ) that are demonstrated in Fig. 6. We therefore need to determine curves y(a 2 ) that separate the (a 2 -y)-plane into regions that correspond to that different behaviour of the characteristic functions q r (r; y, a 2 ) and q ex (r; y, a 2 ). The number of these functions is substantially lowered by the fact that all the local extrema are multiple intersections with other curves and coincide with other extrema. Moreover, as explained bellow, the behaviour of the characteristic functions q r (r; y, a 2 ) and q ex (r; y, a 2 ) in their negative values we can omit as irrelevant for the character of the photon motion. The functions we need are the following: • y max(d(ex)) (a 2 ) = y d(ex)-ex(ex)− (a 2 ) • y ex(ex)-(ex(ex)+) (a 2 ) Here the dashes between two labels denote affiliation to intersection of appropriate functions (it can be proved that there are no other intersections of these functions than that shown in Fig.6). These functions are projections of extremal values or intersections of characteristic functions y(r; a 2 ) into (a 2 -y)-plane and they are demonstrated in Fig. 7.
The functions y ex(h) (a 2 ) divide the parameter plane (a 2 -y) into regions describing Kerr-de Sitter black hole and naked singularity spacetimes, the curve y max(d) (a 2 ) divides spacetimes with so called divergent and restricted repulsive barrier of photon motion. A detailed discussion of these functions have been performed e. g. in [73,86] and will not be repeated here. The significance of the remaining functions can be understood from the depiction of the characteristic functions q r (r; y, a 2 ), q ex (r; y, a 2 ) in Fig 8. They are given parametrically by appropriate functions a 2 (r), y(r; a 2 (r)) with r being the parameter: • the functions y max(h) (a 2 ) and y min(h) (a 2 ) are both determined by a 2 ex(h)+ (r) and y h (r; a 2 = a 2 ex(h)+ (r)); • y max(d) (a 2 ) we obtain from a 2 max(d) (r) with y d (r; a 2 = a 2 max(d) (r)); • the curve y max(d(ex)) (a 2 ) is given by functions a 2 max(d(ex)) (r) and y d(ex) (r; a 2 = a 2 max(d(ex)) (r));

FIG. 7:
The functions y max(h) (a 2 ) (bold full curve), y min(h) (a 2 ) (bold dashed curve), y d-z(ex)+ (a 2 ) (bold dash-dotted curve), y max(d) (a 2 ) (full curve), y max(d(ex)) (a 2 ) (dashed curve), y ex(ex)-ex(ex)+ (a 2 ) (dash-dotted) and y = 1/a 2 (bold dotted curve). In order to clearly display the asymptotic behaviour of these functions we give their plots with the a 2 -axis in logarithmic scale. These curves represent such qualitative changes in behaviour of the characteristic functions qr(r; y, a 2 ), qex(r; y, a 2 ), which are resulting in different character of the photon motion, and divide the (a 2 − y)-plane into regions distinguished by different Roman numerals. The bold grey line is function y (ex(ex)+)-(ex(ex)−)2 (a 2 ) that separates spacetimes I-III, IVa, VIa endowed with both prograde and retrograde spherical photon orbits, as seen by family of the locally non-rotating observers (see bellow), that are separated by 'polar' spherical photon orbit, from the spacetimes IVb, V, VIb,VII, VIII possessing just retrograde spherical orbits. The function y = 1/a 2 represents additional division of parameter plane reflecting different character of the latitudinal motion as discussed in previous section.
• y d-z(ex)+ (a 2 ) is determined by and y d (r; a 2 = a 2 d-z(ex)+ (r)), where the function a 2 d-z(ex)+ (r) is a solution of y d (r; a 2 ) = y z(ex)+ (r; a 2 ) with respect to parameter a 2 ; all such functions are obtained by analogous manner; • y ex(ex)-ex(ex)+ (a 2 ) are constructed from and y ex(ex) (r; a 2 = a 2 ex(ex)-ex(ex)+ (r)); There exist another functions y(a 2 ), corresponding to intersections of the characteristic functions y(r; a 2 ), which are not displayed in Fig.7. The reason is that all the functions y(a 2 ) lie under the curve y = 1/a 2 , and thus we have to take into account the restriction q ≥ −a 2 (see Section 1). Therefore, the changes of the characteristic functions q r (r; y, a 2 ), q ex (r; y, a 2 ) in values under this limit can be omitted as irrelevant. Moreover, we can easily show that in the case q < 0, the restrictions (53) imposed on the latitudinal motion yields stronger constraints on the value X than that given by the relations (59), (60), (61) conditioning the reality of the radial motion. Indeed, for any triad (q < 0, y, a 2 ) there is no intersection of the curves X = X ± (r; q, y, a) with the lines X θ min(−) , X θ max(+) , where X θ min(−) , X θ max(+) are extrema of the functions X θ (m; q, y, a) introduced in Sec.1 (see Fig. 9 e-g,ω). To verify this, it is convenient to regard the curves X = X ± (r; q, y, a) as q = const-slices of the surface q = q max (r; X, a, y), where is an alternative expression of the reality condition R(r; X, q, y, a) ≥ 0,, and search instead for intersections of surfaces q = q max (r; X, y, a) and q = q min (X, y, a), defined by relations (54) -(58). We therefore solve two equations -q max (r; X, y, a 2 ) = q 1 (X), with the result and q max (r; X, a, y) = q 2 (X; a, y), which gives The solution (96) yields X > 0, which, however, does not apply to the case q < 0 for y < 1/a 2 . Moreover, this solution represents touching points of the surface q r (r; X, y, a 2 ) with the parabolic surface q 1 (X) at X = + √ −q, and hence can be omitted even in the case y ≥ 1/a 2 , since theses values lie in the region forbidden by the relations (57)- (58). The solutions (97) are evidently irrelevant, since in stationary regions ∆ r > 0 they are imaginary. The above analysis shows that in the case q < 0 the 'potentials' X ± (r; q, y, a) have values in regions forbidden by reality conditions of the latitudinal motion and hence play no role at all. The limits for impact parameter X of photons with q < 0 are thus given by relation (53); photons satisfying the relation (53) have thus no turning points of the radial motion. In the rest of this treatise we can thus focus on the behaviour of the characteristic functions for q ≥ 0. In Fig. 8 we present all possible variants of behaviour of the characteristic functions q(r; y, a 2 ). These variants involve cases III: y max(d) (a 2 ) ≤ y ≤ y max(h) (a 2 ) for a 2 ≤ 1.08316, or y min(h) (a 2 ) ≤ y ≤ y max(h) (a 2 ) for 1.08316 ≤ a 2 ≤ 1.21202 = a 2 crit ; IVa: y ≤ y min(h) (a 2 ) for 1 ≤ a 2 ≤ 1.08316, or y ≤ y max(d) (a 2 ) for 1.08316 ≤ a 2 ≤ 1.28282, or y ≤ y (ex(ex)+)-(ex(ex)−)2 (a 2 ) for 1.28282 ≤ a 2 ≤ 6 √ 3 − 9 = 1.3923; IVb: y (ex(ex)+)-(ex(ex)−)2 (a 2 ) ≤ y ≤ y max(d) (a 2 ) for 1.28282 ≤ a 2 ≤ 1.3923, or y ≤ y max(d) (a 2 ) for 1.3923 ≤ a 2 ≤ 9, or y ex(ex)-ex(ex)+ (a 2 ) ≤ y ≤ y max(d) (a 2 ) for a 2 ≥ 9; V: y ≤ y ex(ex)-ex(ex)+ (a 2 ) for a 2 ≥ 9; VIa: y max(d) (a 2 ) ≤ y ≤ y min(h) (a 2 ) for 1.08316 ≤ a 2 ≤ 1.21202, or y max(d) (a 2 ) ≤ y ≤ y (ex(ex)+)-(ex(ex)−)2 (a 2 ) for 1.21202 ≤ a 2 ≤ 1.28282; VIb: y (ex(ex)+)-(ex(ex)−)2 (a 2 ) ≤ y ≤ y max(d(ex)) (a 2 ) for 1.21202 ≤ a 2 ≤ 1.28282, or y max(d) (a 2 ) ≤ y ≤ y max(d(ex)) (a 2 ) for a 2 ≥ 1.28282; VII: y max(h) (a 2 ) ≤ y ≤ 1/a 2 for a 2 ≤ 1.21202, or y max(d(ex)) (a 2 ) ≤ y ≤ 1/a 2 for a 2 ≥ 1.21202; Now it remains to assign to each region of the (a 2 -y)plane functions q(y, a 2 ), which by themselves represent marginal values of the parameter q corresponding to some qualitative shift in the behaviour of the potentials X ± (r; q, y, a).
In the regions I, II, which describe black hole spacetimes with divergent repulsive barrier, we have to compare the two local maxima q max(ex+) (y, a 2 ) located under the inner horizon and q max(ex) (y, a 2 ) = q min(r) (y, a 2 ) between the outer and cosmological horizons respectively (see Figs. 8a,b). The function q max(ex+) (y, a 2 ) is given parametrically by functions y ex(ex)+ (r; y, a 2 ) and q ex (r; y = y ex(ex)+ (r; y, a 2 ), a 2 ) with r being the parameter, similarly q max(ex) (y, a 2 ) is given by y ex(ex) (r; y, a 2 ) and q ex (r; y = y ex(ex) (r; y, a 2 ), a 2 ). The extrema function q max(ex) (y, a 2 ) diverges at the curve y max(d) (a 2 ), which forms boundary between regions II-III and IV-VI, i. e. q max(ex) (y = y max(d) (a 2 ), a 2 ) → +∞ (c. f. Figs. 8b, 8c and 8d, 8f). In the region III, corresponding to black hole spacetimes with the restricted repulsive barrier, only local maximum q max(ex+) (y, a 2 ) located under the inner horizon remains. As can be seen from the behaviour of the characteristic functions in Fig. 6b, for y → y min(h) (a 2 ) from above, the 'inner' local maximum of q (ex) (r; y, a 2 ), determined by y ex(ex)+ (r; y, a 2 ), approaches from the left its divergency point given by y = y d(ex) (r; a 2 ), where q ex → −∞ (Fig. 8b), so that q ex (r; y = y min(h) (a 2 ), a 2 ) becomes continuous. For y ≤ y min(h) (a 2 ), the divergency of the function q ex (r; y, a 2 ) appears again with q ex → +∞ and a local minimum has formed on the right (cf. Figs. 8b, 8d). Hence the curve y min(h) (a 2 ) forms a boundary on which the local maxima q max(ex+) (y, a 2 ) convert into local minima. We denote them q min(ex±) (y, a 2 ), since, as follows from relations (86)-(91), for 1.125 ≤ a 2 ≤ a 2 crit and y ≤ y (ex(ex) or a 2 crit ≤ a 2 ≤ 1.3923 and they are given by y ex(ex)− (r; y, a 2 ). The function y (ex(ex)+)-(ex(ex)−)1 (a 2 ) is given parametrically by a 2 r(ex(ex)±) (r) and, e.g., by y ex(ex)− (r; y, a 2 = a 2 r(ex(ex)±) (r)), and the function y (ex(ex)+)-(ex(ex)−)2 (a 2 ) by a 2 r(ex(ex)±)+ (r) and y ex(ex)− (r; y, a 2 = a 2 r(ex(ex)±)+ (r)). The analytical expressions in (98), (99) can be then derived by eliminating the radius r. Both these functions have their relevant parts entirely in regions IV,VI.
The function y (ex(ex)+)-(ex(ex)−)2 (a 2 ) corresponds to the local minima of the potential X − (r; q, y, a) reaching the value X = −a, i. e., = 0. The photons corresponding to these minima, and having appropriate motion constant q, persist on 'spherical' orbits with r = const, which are crossing the spacetime rotation axis alternately above both poles. In the next we shall call them 'polar' spherical orbits -in the following section we shall see that such polar spherical orbits form a border surface between prograde and retrograde spherical photon orbits, as related to the locally non-rotating observers. The function y (ex(ex)+)-(ex(ex)−)2 (a 2 ) has therefore an important meaning, since it represents a boundary between regions of qualitatively different KdS spacetimes in the (a 2 − y)-plane. From this point of view, the function y (ex(ex)+)-(ex(ex)−)2 (a 2 ) creates another qualitative shift in the parameter plane (a 2 −y) with regard to character of the photon motion, however, no qualitative shift in mathematical properties of characteristic functions q r (r; y, a 2 ) and q ex (r; y, a 2 ) in their relevant values q ≥ 0. The parts of the (a 2 −y)-plane corresponding to different behaviour of the characteristic functions are in Fig. 7 distinguished by Roman numerals, the curve y (ex(ex)+)-(ex(ex)−)2 (a 2 ) then induces an additional division a/b. Further we have to relate the minima q min(ex±) (y, a 2 ) with the maxima q max(ex) (y, a 2 ) = q min(r) (y, a 2 ). In the region V, the minima of q ex (r; y, a 2 ) coalesce with the minima of q r (r; y, a 2 ) (Fig.  8e). We therefore have to compare the minima function q min(ex) (y, a 2 )=q min(r) (y, a 2 ) determined by y ex(ex) (r; y, a 2 ) and, e.g., q ex (r; y = y ex(ex) (r; y, a 2 ), a 2 ), with the maxima function q max(ex+) (y, a 2 ) parametrized by y ex(ex)+ (r; y, a 2 ) and q ex (r; y = y ex(ex)+ (r; y, a 2 ), a 2 ). The boundary curve y ex(ex)-ex(ex)+ (a 2 ) then represents such combinations of parameters a 2 , y for which the local extrema of q ex (r; y, a 2 ) have coalesced into an inflection point.
For parameters from the region VI, corresponding to naked singularity spacetimes with restricted repulsive barrier (as well as from the remaining regions VII, VIII), the function q ex (r; y, a 2 ) has one local minimum (Fig. 8f), and we therefore construct a function q min(ex±) (y, a 2 ) determined by functions y ex(ex)± (r; a 2 ) and q ex (r; y = y ex(ex)± (r; a 2 ), a 2 ), where the minus sign has to be chosen for 1.17007 ≤ a 2 ≤ a 2 crit and y max(d) (a 2 ) ≤ y ≤ y (ex(ex)+)-(ex(ex)−)1 (a 2 ), or a 2 crit ≤ a 2 ≤ 1.2828 and y max(d) (a 2 ) ≤ y ≤ y (ex(ex)+)-(ex(ex)−)2 (a 2 ) (see Fig.7). For y → y max(d(ex)) (a 2 ) and a 2 ≥ a 2 crit there is q min(ex+) (y, a 2 ) → +∞,, and for y > y max(d(ex)) (a 2 ), i.e. in the region VII, it converts into the local maximum (c. f. Figs 8f, 8g). The transition into the region VII from the region III can be inferred from comparison of Fig. 8c with Fig. 8g. Therefore, in the region VII we have to follow up the values of the function q max(ex+) (y, a 2 ) determined by the functions y ex(ex)+ (r; a 2 ) and q ex (r; y = y ex(ex)+ (r; a 2 ), a 2 ). The functions q(y, a 2 ) are demonstrated in Fig. 9.
With the knowledge of the behaviour of the extremal values q min/max(ex) (y, a 2 ) at each region of the (a 2 -y)plane, we can finally construct all qualitatively different types of the behavior of the effective potentials X ± (r; q, y, a). They are presented in Figure 10 for appropriately chosen representative combinations (q, y, a 2 ).

V. SPHERICAL PHOTON ORBITS AND CLASSIFICATION OF THE KERR-DE SITTER SPACETIMES DUE TO PROPERTIES OF THE PHOTON MOTION
In the following section we demonstrate by using the behaviour of the effective potentials X ± (r; q, a, y) that the null geodesics create qualitatively different structures in the various cases of Kerr-de Sitter spacetimes with the spacetime parameters chosen from different parts of the (a 2 -y)-plane labelled by numerals I-VIII. Hence the regions of the spacetime parameter space of these labels can be considered as representatives of the classification of the Kerr-de Sitter spacetimes due to the photon motion (null geodesics). Similarly as in [73], there are three (four) criteria used -the main criterion for the classification is the existence (number) of the event horizons. The other differentiating factors follow from the nature of photon motion. First, there is some kind of repulsive barrier defending a light to reach the ring singularity, which is always created in its vicinity for photons with q > 0. However, a similar barrier can emerge between the outer black hole horizon and the cosmological horizon in black hole and naked singularity spacetimes, repelling photons towards one of these horizons. In the naked singularity spacetimes, occurrence of an additional barrier, which reflects photons towards the ring singularity, leads to occurence of the phenomenon of bound photon orbits. Such bound photon orbits are not present in the case of the black hole spacetimes. The presence and character of this barrier we take as another criterion in the following classification. The other aspect that authorizes us to make such distinction between the KdS spacetimes will be the existence and character of the spherical photon orbits. In the KdS naked singularity spacetimes the bound orbits are concentrated around the stable spherical photon orbits.  r; q, y, a). The dynamic region is highlighted by shading, the grey vertical bar demarcates loci of the static radius. Note that the last scheme corresponding to case y > 1/a 2 , which is introduced due to changes in the latitudinal motion, brings no new features in character of the radial motion, since it qualitatively coincide with case VII and is shown for completeness. -slices corresponding to values between outstanding limits a 2 = 1, a 2 rrb(max) = 1.08316, a 2 crit = 1.21202 and a 2 (ex(ex)-ex(ex)+)min = 9. The positive branch of the full curve is the graph of the function q max(ex) (y, a 2 ), its asymptote intersects the (a 2 -y)-plane in the curve y max(d) (a 2 ). The dashed curve in Fig. (b) describes the local maxima q max(ex+) (y, a 2 ) under the inner horizon. In Figs. (c), (d) it represents local minima q min(ex+) (y, a 2 ) (q min(ex−) (y, a 2 ) eventually) for y < y min(h) (a 2 ), and maxima q max(ex+) (y, a 2 ) for y > y min(h) (a 2 ), where the critical value q max(ex+) (y = y min(h) (a 2 ), a 2 ) can be identified as the inflection point of this curve. The minima diverge as y → y max(d(ex)) (a 2 ) (Fig. e), the descending part then describes the local maxima q max(ex+) (y, a 2 ) in region VII. Fig. (f) demonstrates intersection of both curves at y = y ex(ex)−ex(ex+) (a 2 ). The full curve for y < y ex(ex)−ex(ex+) (a 2 ) then corresponds to the local minima q min(ex) (y, a 2 ) = q min(r) (y, a 2 ), while the dashed one now matches the maxima q max(ex+) (y, a 2 ). , and X−(r; q, a, y) (dashed curve), with spacetime parameters chosen successively from regions corresponding to the classes I-VIII in the (a 2 − y)-plane, and some representative values of parameter q divided by the extrema q(y, a 2 ). Vertical dashed lines demarcate loci of the event horizons, the grey line demarcate the static radius, shading highlights the forbidden region. Figs. (e)-(g) depict the case q < 0 with the limiting values X θ min(−) , X θ max(+) as horizontal dashed lines. The graphs X±(r; q, a, y) lie entirely in forbidden region. The behaviour of these functions in other cases differing by values a, y with q < 0 is qualitatively the same. Figs. (γ)-(δ) describes the case y > 1/a 2 , q < 0. The potentials have values in forbidden region again; for q > 0 the behaviour of the potentials qualitatively corresponds to class VII. In the classes IVb, VIb, the local minima of the potentials X−(r; q, a, y) (see Figs. (q),(y)) have values X min(−) < −a, contrary to the case IVa,VIa, which are not displayed, since there are no other qualitative differences.

A. Spherical photon orbits
The spherical photon orbits are determined by the conditions R(r) = 0 and dR/dr = 0 that have to be solved simultaneously. The physically acceptable solution is governed by the relations for the photon motion constants X and q that are expressed as functions of the radius r and the spacetime parameters a, y, and take the form where the function q ex (r; y, a 2 ) is defined by the relation (75). These solutions governing the spherical photon orbits are allowed in the interval of radii limited by the equatorial photon circular orbits. Stability of the spherical photon orbits relative to radial perturbations is determined by the sign of the expression evaluated at appropriate radii. It can be shown that local maxima of the potential X + and local minima of X − correspond to unstable orbits in the black hole spacetimes. However, local minima of X + and local maxima of X − represent stable orbits, which occur in the naked singularity spacetimes.
It is useful to relate the parameters (wave vector components) of photons orbiting along the spherical null geodesics to the locally non-rotating frames (LNRFs) that are the most convenient frames for description of physical processes in the Kerr(dS) spacetimes [9]. The LNRF tetrad of differential one-forms is given by the relations the corresponding tetrad of dual vectors reads The wave-vector components related to the LNRFs are then determined by the relations hence and is the angular velocity of the LNRFs related to distant static observers. In order to determine the orientation of the spherical orbits, we have chosen as the azimuthal direction indicator the sign of the ratio k (φ) /k (t) . If we define the directional angle Ψ in such a way that Ψ = 0 for motion in the direction of the latitudinal tetrad vector e (θ) , while Ψ = π/2 for motion in the direction of the azimuthal tetrad vector e (φ) , then k (φ) /k (t) = sin Ψ and we find the relation If the sign of sinΨ is positive, we call the spherical orbit prograde, if it is negative, we call the spherical orbit retrograde. Special case of limiting spherical orbits corresponds to the equatorial circular orbits that are again corotating (prograde), respectively counter-rotating (retrograde). It can be shown that the sign of the directional angle remains fixed at any latitude of any particular spherical orbit, i.e., the locally non-rotating observers see the photon motion in fixed azimuthal direction. [94] Since the functions A, Ω LN RF are positive [73], it is clear that all photons with X < −a ( < 0) are retrograde. However, photons with can be retrograde as well. Considering in such a case the relation for the tetrad LNRF component (112), we can see that in order to keep for k (t) the standard physical meaning, i.e., k (t) > 0, we have to put E < 0. [95] In order to find conditions under which such a situation occurs, it is convenient to express from the alternate relation the impact parameter in the form and reverse the problem by searching for conditions, under which there is > 0. Such a relation is evidently fulfilled if sin ψ > 0, i.e., the positive impact parameters pertain to prograde photons. However, there is another possibility, to have sin ψ < 0 together with from which it follows However, the last inequality can be written in the form which implies Hence, such a situation can occur only in the ergosphere. Of course, the impact parameter of such photons must fulfil the condition (119). The function has for ∆ r = 0 common points with the potentials X ± given by the relation (64). There are no other intersections with the potentials, hence, the reality condition of the radial motion together with (119) imply X > X + > 0. Therefore, the motion of photons with negative energy E, which appear to be retrograde in the LNRFs, is governed by effective potentials X + (r; q, a, y) with positive values.
Using the properties of the effective potentials X ± (r; q, a, y), we can identify the radii r = r 0 of the spherical photon orbits as loci of the local extrema of the effective potentials and determine their stability and orientation as described above. At each allowed radius r 0 , located between the radii of the equatorial photon circular orbits, we can assign corresponding limits θ min , θ max on the latitudinal motion by solving the equation which due to he results of Section 3 has one real positive root m 0 , since q sph (r 0 ) ≥ 0 (q sph (r 0 ) = 0 for r 0 = r ph± , i.e., equatorial circular co-rotating or counter-rotating photon orbit). The marginal latitudes (turning points of the latitudinal motion) then read for details see the discussion of the latitudinal motion in Section 3. We can thus easily determine for a spherical orbit at an allowed radius r 0 the impact parameters of the orbit, the extension of the latitudinal motion, and the orientation of the azimuthal motion.

B. Classification
In the following classification we introduce ten classes of the Kerr-de Sitter spacetimes and demonstrate properties of the photon motion using the spherical photon orbits that serve as crucial characteristic for the classification. We give the loci of the spherical photon orbits and their extension in latitude, stability against radial perturbations, and orientation of their azimuthal motion. The classification is represented by family of characteristic figures corresponding to the separated classes of the KdS spacetimes. For easy interpretation of the family of the figures representing the classification, we introduce an auxiliary Fig.11 commented with detailed explanatory notes. In order to fully and clearly characterize the KdS spacetimes and their horizon and ergosphere structure, and to demonstrate the spheroidal character of the applied coordinate system, we use now the so called Kerr-Schild coordinates x, y, z that are connected to the Boyer-Lindquist coordinates r, θ by the relations x 2 + y 2 = (r 2 + a 2 ) sin 2 θ, z 2 = r 2 cos 2 θ.
In the figures we, of course, use the meridional sections of y = 0. The characteristics of the classes of the KdS spacetime according to the photon orbits are presented as follows.
Class I : Black hole spacetimes with the divergent repulsive barrier of the radial photon motion, having The ergosphere is denoted by shading; note however that the shading is canceled by colours corresponding to various types of the spherical photon orbits located at the ergosphere (here and in the following figures representing the classification). The outer ergosphere and the cosmological horizon is outside of the drawing in this figure. Thin dashed red ellipse denotes the sphere with radius r = r ph+ = 1.98, which intersects the equatorial plane in the co-rotating photon circular orbit, displayed as red bold points. Similarly, thin dashed blue ellipse represents the r = r ph2 = 3.43-sphere, and the bold blue points mark the equatorial counter-rotating circular orbit. Bold red arcs represent a prograde spherical photon orbit at r = 2.26, selected just to touch the ergosurface in the equatorial plane. Thick red dashed arcs denote rest of the selected sphere. Accordingly, full parts of the bold blue ellipse mark some spherical retrograde photon orbit, selected at r = 3.1. The endpoints of the bold arcs then correspond to marginal latitudes, among which the photons oscillate. All such endpoints thus create a surface, here depicted by full thin red or blue curves, limiting regions filled by the prograde or retrograde, respectively, spherical photon orbits, coloured in light red or blue, respectively. In this case, the regions of opposite orientations are separated by purple ellipse representing the polar spherical orbit at which a photon with X = −a ( = 0) crosses the spaxcetime rotation axis alternately above both poles. Light/full hues we reserved to regions of unstable/stable orbits, the latter ones not present in this case. In some KdS spacetimes, there exist retrograde spherical orbits of photons with X > −a. They are confined exclusively inside the ergosphere and occupied by photons with negative energy E. These regions will be depicted in green.
one equatorial counter-rotating circular unstable orbit with negative energy located under the inner black hole horizon (0 < r < r − ), which is limiting the range of the spherical photon orbits with negative energy. There exist stable orbits, corresponding to local minima X min(+) of the effective potential X + at 0 < r < r max(ex)1 , and unstable orbits, corresponding to local maxima X max(+) of X + at r max(ex)1 < r < r z(ex)1 for 0 < q < q max(ex+) (y, a 2 ) (Fig. 12a). We denote as r min/max(ex) the local extrema, and as r z(ex) the zero point, of the function q ex (r; y, a 2 ) hereafter. Such a structure is present under the inner horizon of any KdS black hole spacetime. Outside the ergosphere, one unstable co-rotating equatorial circular orbit, located at r = r ph+ = r z(ex)2 , and polar spherical orbit with r = r pol , r ph+ < r pol , limit the range of unstable prograde spherical orbits given by the local minima X min(−) of the effective potential X − , for which X min(−) > −a. The radius of the polar spherical orbit is found by solving X − (r pol ; q ex (r pol )) = −a. The counter-rotating equatorial circular orbit at r = r ph− = r z(ex)3 gives the limit of region of unstable retrograde spherical orbits, given by the local minima X min(−) < −a, and maxima X max(+) < −a for 0 < q < q max(ex) (y, a 2 ), such that r pol < r ph− .
Class II : Black hole spacetimes with the same features as in the class I, but now the ergosphere enters the region of the spherical photon orbits (Fig. 12b). No spherical orbit is fully immersed in the ergosphere and photons at all the spherical orbits have positive energy. The presence of the ergosphere in region of the spherical photon orbits influences character of the light escape cones [18].
Class III : Black hole spacetimes with the restricted repulsive barrier of the radial photon motion. The ergosphere spreads over all radii. The prograde spherical orbits are given by the local minima X min(−) > −a at r ph+ < r < r pol , while the retrograde spherical orbits with E > 0 are given by the minima X min(−) < −a at r pol < r < r d(ex) . The spherical orbits given by the local maxima X max(+) (see Figs 10k-n) at r d(ex) < r < r ph− , where r d(ex) denotes the divergence point of q ex (r; y, a 2 ) (see Fig. 8c), are fully immersed in the ergosphere. Such areas are drawn in Figs 12 in green and the spheres with r = r d(ex) as full/dashed green ellipses. Photons in such regions have E < 0.
Class IVa : Naked singularity spacetimes with divergent repulsive barrier of the photon motion. At radii 0 < r < r d(ex) (Fig.8d can be used for illustration), there are local minima of the potential X + (for illustration use Figs. 10p-r) corresponding to the stable retrograde spherical orbits with negative energy (E < 0) (Fig. 12d). The stable retrograde orbits with positive energy corresponding to the local maxima of X − are at r d(ex) < r < r pol1 . These maxima exceed the value X max(−) = −a at r pol1 < r < r min(ex) , where they yield stable prograde spherical orbits. At radii r min(ex) < r < r pol2 , there are the local minima of X − with values X min(−) > −a, -these radii are thus occupied by the unstable prograde orbits. The local minima of X − , and the local maxima of X + at r > r pol2 , correspond to the unstable retrograde orbits. There are thus two polar spherical orbits enclosing region of prograde orbits -the inner at the radius r = r pol1 being stable, the outer at the radius r = r pol2 being unstable.
Class IVb : Naked singularity spacetimes with the same features as in the class IVa, but the two polar orbits have coalesced, therefore, there are no prograde spherical orbits (Fig. 12e).
Class V : Naked singularity spacetimes having the structure of the spherical orbits corresponding to the previous case (Fig. 12f), but with is a small region of bound orbits for photons with motion constants q min(ex) (y, a 2 ) < q < q max(ex+) (y, a 2 ) and X between the appropriate local extrema of X + (c. f. Figs. 10q, u), which is not contained in the other cases.
Class VIa : Naked singularity spacetimes with the restricted repulsive barrier of the radial photon motion. For 0 < r < r d(ex)1 (the function q ex (r; y, a 2 ) has two divergence points r d(ex)1 , r d(ex)2 , -see Fig. 8f for preview) the minima of X + correspond to stable retrograde orbits with E < 0; for r d(ex)1 < r < r pol1 , there are the local maxima of X − with values X max(−) < −a giving retrograde orbits with E > 0. The local minima of X − at r pol1 < r < r min(ex) give the stable prograde orbits. At radius r = r pol1 the stable polar orbit is located. For r min(ex)<r<r pol2 , the function X − has minima with values X min(−) > −a giving the unstable prograde orbits. For r pol2 < r < r d(ex)2 , they correspond to the unstable retrograde orbits. At radius r = r pol2 , the unstable polar orbit exists. The local maxima of the function X + at r d(ex)2<r<r ph− correspond to the retrograde unstable spherical orbits with E < 0.
Class VIb : The structure of the spherical orbits corresponds to the class VIa with an exception that the local extrema of the potential X − have values X < −a,, implying that there are no polar spherical orbits, neither the prograde spherical orbits (Fig.  12h).
Class VII : Naked singularity spacetimes with the restricted repulsive barrier of the radial photon motion having stable retrograde spherical orbits at 0 < r max(ex) corresponding to local minima of X + (Fig. 10α), and unstable retrograde spherical orbits at r max(ex) < r < r ph− corresponding to local maxima of X + . All these spherical orbits, including the counter-rotating equatorial circular orbit at r = r ph− , correspond to photons with E < 0.
Class VIII : Special class of the naked singularity spacetimes demonstrating the same features of the radial motion of photons with q ≥ 0 as the class VII, but differing from all previous cases by the existence of null geodesics for arbitrary q < 0. The allowed values of the impact parameter X are for q < 0 confined to the intervals X < X θ max(+) < 0 or X > X θ min(−) > 0 (see Section 3). The potentials governing the radial photon motion are fully immersed in the forbidden region (Figs. 10(γ), (δ)), thus in the radial direction the photons with such parameters move freely in the whole range between the ring singularity and the cosmological horizon.

VI. CONCLUSIONS
We can summarize our results by the following concluding remarks.
1. In any kind of the black hole spacetimes, there are no radially bound null geodesics in the stationary region, i.e., the trajectory of a photon has at most one turning point in radial direction between the outer and cosmological horizon, or the photons can move freely between the outer black hole and the cosmological horizons. However, such bounded photon orbits exist in each naked singularity spacetime for photons with parameters q > 0 and X chosen appropriately.
2. No photons with q > 0 can reach the ring singularity at r = 0 in any of the Kerr-de Sitter spacetimes.
3. In the Kerr-de Sitter spacetimes of classes I-VII, i.e., with the spacetime parameters satisfying the condition y < 1/a 2 , there is a lower limit q = −a 2 of the parameter q < 0, for which the photon motion is allowed. The range of the allowed values of the impact parameter X is then an interval given by the relations (53) - (56). Photons with such tuned parameters have no turning point in radial direction, since the effective potential lies entirely in the forbidden region (Figs 10(e)-(g)). Further, by the results of Section 3, only such photons execute the vortical motion, or their trajectory lies completely on the cones of θ = constant. We can therefore reject possibility of existence of vortical photon motion of constant radius, or off-equatorial circular photon orbits.
4. In the Kerr-de Sitter spacetimes of class VIII (y > 1/a 2 ), the photon motion is allowed for any q < 0.
The permissible values of the parameter X are then two disjunct unlimited intervals determined by the relation (58.) In the extreme case y = 1/a 2 , it must be q ≥ −a 2 again and for negative q, the parameter X can take less than certain negative value given by (57). The consequences for photon motion are then the same as in previous note (Figs 10(γ)-(δ)).
5. In the Kerr-de Sitter spacetimes with the divergent repulsive barrier of the radial photon motion, there exists a critical value q max(ex) (y, a 2 ), for which this barrier becomes impermeable between the outer black hole horizon and cosmological horizon, or, in naked singularity spacetimes, between the ring singularity and cosmological horizon, for photons with any impact parameter X. In spacetimes with the restricted repulsive barrier of the radial photon motion, the height of this barrier slowly grows with increasing parameter q, but stays finite for any q > 0 (Figs 10(n), 10(y)).
6. In the Kerr-de Sitter spacetimes of classes I-III, IVa, VIa, there exist spherical photon orbits, which can be both prograde or retrograde as seen by the family of locally non-rotating observers. Additionally, each of the two types can be stable or unstable with respect to radial perturbations. The regions of spherical orbits of different orientations are separated by the so called polar spherical orbit, at which photons cross the spacetime rotation axis alternately above both poles. In the naked singularity spacetimes of class IVa, VIa, there are two polar spherical orbits, the inner one being stable, the outer one being unstable.
7. In the Kerr-de Sitter spacetimes of classes IVb, V ,VIb, VII, VIII there are no prograde or polar spherical orbits. 8. In each class of the Kerr-de Sitter spacetimes, there exist region, where the effective potential X + have positive values. Photons with impact parameter X exceeding this values appear to move in retrograde direction as seen in LNRFs. This region must be located inside the ergosphere, and photons with such impact parameters must have negative energy, E < 0. In the black hole spacetimes with the divergent barrier of the radial photon motion, the ergosphere has two parts above the black hole outer event horizon -the inner one, which is limited to the outer vicinity of the outer event horizon, and the outer one, limited to the inner vicinity of the cosmological horizon. In the black hole spacetimes with restricted repulsive barrier of the radial photon motion, the two regions of the ergosphere merge in the equatorial plane, and they spreads at any radii except for certain region in the vicinity of the rotation axis.
9. In the LNRFs, trajectories of photons moving along any spherical orbit have no turning point of the azimuthal motion.
We have thus demonstrated a variety of very extraordinary phenomena related to the photon motion in the KdS spacetimes, of both black hole and naked singularity types. Especially relevant effects are found in the case of spherical photon orbits that can be directly related to the observational phenomena. It is quite interesting that we could expect another interesting phenomena related to the charged Kerr-Newman or Kerr-Newman-de Sitter naked singularity spacetimes (with both the standard electric charge, or the tidal charge of the braneworld models), especially in the case of the so called mining Kerr-Newman spacetimes [12], containing a special type of equatorial stable photon orbits.