Photon motion in Kerr-de Sitter spacetimes
Daniel Charbulak Zdenek Stuchlik








Schwarzschild
Schwarzschild–de Sitter
de Sitter
Schwarzschild–anti-de Sitter
Kerr–de Sitter
Kerr–Newman–de Sitter

    dd E Kma-xi-mumpa-ra-me-tersga-la-xiesgeo-met-ry sign daniel.charbulak@fpf.slu.czzdenek.stuchlik@fpf.slu.czInstitute of Physics and Research Centre of Theoretical Physics and Astrophysics, 
        	Faculty of Philosophy and Science, 
        	Silesian university in Opava, 
        	Bezručovo nám. 13, CZ-746 01 Opava, Czech Republic
        We study general motion of photons in the Kerr-de Sitter black hole and naked singularity spacetimes. The motion is governed by the impact parameters X, related to axial symmetry of the spacetime, and 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. Restrictions on the impact parameters X and q are established in dependence on the spacetime parameters M, Λ, 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 the radial motion, excluding existence of vortical motion of constant radius. Properties of the spherical photon orbits of the orbital type are determined and used along with properties of the effective potentials as criteria of classification of the KdS spacetimes according to the properties of the photon motion.    
        Photon motion in   spacetimes
            Zdeněk Stuchlík
            December 30, 2023
        =============================
        
        
        
        § INTRODUCTION
        
        
        In the framework of inflationary paradigm <cit.>, 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 <cit.>. 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 <cit.>. 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 <cit.>, and the dark energy represents about 70% of the energy content of the observable universe <cit.>. These conclusions are confirmed by recent measurements of cosmic microwave background anisotropies by the space satellite observatory PLANCK <cit.>.
        
        The dark energy equation of state is very close to those corresponding to the vacuum energy <cit.>. 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 , naked singularity, or any compact-body backgrounds as such backgrounds become asymptotically  spacetimes, and an event horizon (cosmological horizon) always exists, behind which the geometry is dynamic. 
        
        Substantial influence of the repulsive cosmological constant has been demonstrated for astrophysical situations related to active galactic nuclei and their central supermassive black holes <cit.>. The black hole spacetimes with the Λ term are described in the spherically symmetric case by the vacuum (SdS)
        
         geometry <cit.>, while the internal, uniform density SdS spacetimes are given in <cit.>. The axially symmetric, rotating black holes are determined by the Kerr-de Sitter (KdS) geometry <cit.>. In the spacetimes with the repulsive cosmological term, motion of photons was extensively investigated in many papers <cit.>. The motion of massive test particles was studied in <cit.>. The KdS geometry can be relevant also for the so called Kerr superspinars representing an alternative to black holes <cit.>, breaking the black hole bound on the dimensionless spin and exhibiting a variety of unusual physical phenomena <cit.>. It is worth to note that the SdS and KdS spacetimes are equivalent to some solutions of the f(R) gravity representing black holes and naked singularities <cit.>. 
        
        The role of the cosmological constant can be significant for both the geometrically thin Keplerian accretion discs <cit.> and the toroidal accretion discs <cit.> 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 <cit.>. In the spherically symmetric spacetimes, the Keplerian and toroidal disc structures can be precisely described the Pseudo-Newtonian potential of Paczynski type <cit.> that appears to be useful also in studies of motion of interacting galaxies <cit.> demonstrating relation of the gravitationally bound galactic systems to the so called static radius of the SdS or KdS spacetimes  <cit.>. This idea has been confirmed by the recent study of general relativistic static polytropic spheres in spacetimes with the repulsive cosmological constant <cit.>. 
        
        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 <cit.>, 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. 
        
        
        
        
        
        
        §  SPACETIME AND CARTER'S EQUATIONS OF GEODESIC MOTION
        
        
         §.§  geometry
        
        The line element describing the  geometry is in the standard Boyer-Lindquist coordinates, using geometric system of units (c = G = 1), given by
         s^2 =
             -   Δ_r/I^2 ρ^2( t - a sin^2θϕ)^2
           +   Δ_θsin^2θ/I^2ρ^2[a  t - (r^2 +   a^2)ϕ]^2 
           +   ρ^2/Δ_r r^2 + ρ^2/Δ_θθ^2, 
        where
        Δ_r   =   (1 - 1/3Λ r^2)(r^2 + a^2) - 2Mr,
        Δ_θ   =   1 + 1/3Λ a^2 cos^2θ,
        
        I   =   1 + 1/3Λ a^2,
        ρ^2   =   r^2 + a^2cos^2θ. 
         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 - y r^2)(r^2 + a^2) - 2r,
        Δ_θ   =   1 + a^2 y cos^2θ,
        
        I   =   1 + a^2 y,
        with equation (<ref>) 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 
        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
        
        a^2=a^2_z(h)(r)≡ 2r-r^2, 
        the loci of its extrema are given by the functions
        
        a^2=a^2_ex(h)±(r)≡r(1-2r±√(1+8r))/2,
        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) 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.212 02, for which the two local extrema of the function y_h(r; a^2) coalesce in an inflection point at r_crit=1.616 03 with the critical value y_crit=0.0592. Thus, for a^2>a^2_crit only  naked singularity can exist for any y>0.
        
        Properties of the event horizons for the more general case of the  spacetimes can be found in <cit.>.
        
        
        
         §.§ Carter's equations of geodesic motion
               
        The motion of test particles and photons following its geodesics in the  spacetime is described by the well known Carter equations <cit.>ρ^2 θ/λ   =   ±√(W(θ; E, Φ, 𝒦, y, a)), 
        ρ^2  r/λ   =   ±√(R(r; E, Φ, 𝒦, y, a)),
        ρ^2 φ/λ   =   aI^2[E(r^2+a^2)-aΦ]/Δ_r
           -   I^2[aEsin^2θ-Φ]/Δ_θsin^2θ,
        ρ^2  t/λ   =   I^2(r^2+a^2)[E(r^2+a^2)-aΦ]/Δ_r
           -   aI^2[aEsin^2θ-Φ]/Δ_θ,
        where
        
        W(θ; E, Φ, 𝒦, y, a) = 𝒦Δ_θ-I^2(aEsin^2θ-Φ)^2/sin^2θ
        and
         R(r; E, Φ, 𝒦, y, a) = [IE(r^2 + a^2) - IaΦ]^2 - Δ_r 𝒦.
        Here E and Φ are the constants of motion connected respectively with the time and axial symmetry of the  geometry, and 𝒦 is the fourth Carter constant of motion connected with the hidden symmetry of the  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 <cit.>.
        Detailed discussion of the equatorial motion of photons in the  spacetimes has been published in <cit.>. Circular motion of test particles in the  spacetimes has been presented in <cit.>. Here we restrict our attention on the general motion of photons in the  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 𝒦/E^2. For our general discussion it is convenient to use Q = 𝒦 - I^2(Φ - a E)^2 that vanishes for the equatorial motion. For our purposes it is, however, following the paper <cit.>, convenient to introduce a new constants of motion X ≡ℓ - a. Further the constant q ≡ Q/I^2E^2 will be applied. Then the relations (<ref>) and (<ref>) simplify to the form
         
         W(θ; X, q, y, a) ≡ I^2E^2[(X^2+q)Δ_θ-(acos^2θ+X)^2/sin^2θ],
         R(r; X, q, y, a)≡ I^2E^2[(r^2 - aX)^2 -Δ_r(X^2 + q)].
         Following the work <cit.> 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  spacetimes has been investigated yet <cit.>; however, the discussion has been related to the motion constant 𝒦. 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.
         
        
        § LATITUDINAL MOTION
        
         Because it is more convenient to work with algebraic functions instead of trigonometric ones, we introduce a new variable 
         
            m = cos ^2 θ,
        
            m=2(θ-π/2)√(m(1-m))dθ.
        
         This implies replacing the equation (<ref>) by
        ρ^2  m/λ = ± 2 √(M(m; X, q, y, a)),
        where
        
        M(m; X, q, y, a)≡ I^2E^2m[(1-m)(X^2+q)Δ_m-(am+X)^2] 
        with notation Δ_m=1+a^2ym. Note that m/λ =0 does not necessarily imply θ/λ=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 (<ref>).
        
        The reality condition M(m; a, y, X, q)≥0 can be expressed by the relations
        
         X^θ_-(m; q, y, a)≤ X≤ X^θ_+(m; q, y, a) 
          in regions where Δ_m-a^2y > 0, i. e., equivalently, 
            m > m_d,
         where 
          
          m_d = 1-1/a^2y 
           is the solution of the equation
          Δ_m-a^2y = 0,
          
          and by the relations
         
          X≤ X^θ_+(m; q, y, a),   X^θ_-(m; q, y, a)≤ X, 
         in regions where Δ_m-a^2y < 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
                X^θ_±(m; q, y, a) ≡
              -am ±√(m(1-m)Δ_m [a^2m+q(Δ_m-a^2y)])/m(Δ_m-a^2y).
         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 <cit.>, 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
         
         q   ≥    q^θ_r(m; y, a^2)    Δ_m-a^2y > 0,
        
         q   ≤    q^θ_r(m; y, a^2)    Δ_m-a^2y < 0, 
         where
         
         q^θ_r(m; y, a^2)≡a^2m/a^2y-Δ_m.
         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
         
         X^θ_c=X_(±)(m; q=q^θ_r, y, a) = X_(±)(m; y, a)≡a/a^2y-Δ_m.
         
         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
         
         y=y^θ_d(m; a^2)≡1/a^2(1-m).
         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 (<ref>) in an alternative form
                X^θ_±(m; q, y, a)=
              a^2m^2-q(1-m)Δ_m/-am∓√(m(1-m)Δ_m[a^2m+q(Δ_m-a^2y)]).
         
         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 (<ref>) 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^2y→∞. In such case, q^θ_r(m; y, a^2)→ +∞ (-∞) for m→ m_d from the left (right). 
         From the equality
         ∂ q^θ_r/∂ m = a^2(a^2y-1)/(Δ_m-a^2y)^2
         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.
        
         In the special case y=1/a^2 we get
         
            q^θ_r(m; y=1/a^2, a^2)=const.= -a^2     m≠ 0
        
         with
         
            lim_m→ 0 q^θ_r(m; y=1/a^2, a^2)= -a^2.
        
         Based on the conditions (<ref>), (<ref>) 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^2m^2+q)(Δ_m-a^2y)[a^2m^2I^2+q(1-a^2y+2a^2my)^2]=0.
         
         It can be verified that the function X^θ_+(m; q, y, a) has local extrema given by the relation
          q=q^θ_ex(+)(m; a^2)≡ -a^2m^2. 
          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
         ∂^2X^θ_+/∂ m^2(m; q=q_ex(+), y, a)=-a/m(1-m)Δ_m,
         
         clearly they must be maxima. 
         From the equation (<ref>) we find that another extrema of the potentials X^θ_±(m; q, y, a)  are determined by the condition
          q=q^θ_ex(±)(m; y, a^2)≡-a^2m^2I^2/[Δ_m-a^2y(1-m)]^2.
         The divergencies of the functions q^θ_ex(±)(m; y, a^2) are determined by the relation
         
         y=y^θ_d(ex±)(m; a)≡1/a^2(1-2m).
         The function y^θ_d(ex±)(m; a^2) is positively valued at 0≤ m<0.5, where y^θ_d(ex±)(m; a^2)→ +∞ for m→ 0.5 from the left. For m=0, there is y^θ_d(ex±)(0; a^2)=y^θ_d(r)(0; a^2)=1/a^2.
        
         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 
            m=m_d(ex)≡ 0.5 (1-1/a^2y)=0.5 m_d
        
         located such that 0≤ m_d(ex) <0.5. Obviously
         
         q^θ_ex(±)(m; y, a^2) → -∞     m→ m_d(ex). 
         In the following we shall decide about the monotony and  possible existence of local extrema of the function q_ex(±). From 
         ∂ q^θ_ex(±)/∂ m = 2m(a^2y-1)a^2I^2/[Δ_m-a^2y(1-m)]^3
         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 (<ref>).
        
         For y<1/a^2, there is
         
            ∂ q^θ_ex(±)/∂ m < 0     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).
        
         Comparing both the functions q^θ_ex+(m; a^2) and q^θ_ex(±)(m; y, a^2), we find that
         
            q^θ_ex(±)(m; y, a^2)≤ q^θ_ex+(m; a^2)≤ 0
        
         and have common points at m = 0, 1 with
         
            q^θ_ex+(0; a^2)= q^θ_ex(±)(0; y, a^2)=0
        
         and
         
            q^θ_ex+(1; a^2)= q^θ_ex(±)(1; y, 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 
               ∂ X^θ_±/∂ m(m; q=q^θ_ex(±), y, a)=
              a^3y{[1-a^2y(1-2m)]±(1-a^2y)}/[1-a^2y(m-1)]^3. 
         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)     y>1/a^2,
        
         which ensures
         ∂^2X^θ_+/∂ m^2(m; q=q^θ_ex(±), y, a)   =   
        aI^2/m(1-m)Δ_m(1-a^2y)[Δ_m-a^2y(1-m)]   <   0. 
         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 (<ref>), but now with the above conditions we have
         
            ∂^2X^θ_-/∂ m^2(m; q=q^θ_ex(±), y, a)>0,
        
         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 (<ref>), (<ref>) 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 (<ref>), (<ref>), 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 = 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  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:
          *  Case y<1/a^2
          * q<-a^2
          – the definition range of the potentials is an empty set; the latitudinal motion is not possible;
            
          * q=-a^2
          – the potentials X^θ_±(m; q, y, a) are defined only for m=1, where
              
            X^θ_+(1; q, y, a)=X^θ_-(1; q, y, a)=-a;
        
               photons with such values of parameters are the special
                case of the so called PNC photons 'radially' moving along the spin axis <cit.>;
            
          * -a^2<q<0 (Fig. 2a, 2b)
            
          –both the potentials are defined for m∈⟨ m_l;1⟩, where the lower limit
            
            m_l = q(a^2y-1)/a^2(qy+1)>0 
            is the solution of the equation
            
            q=q^θ_r(m; y, a)  ;
         
                        
            the limits of the interval are the common points of the potentials, where
            
            X^θ_±(m=m_l; q, y, a)=X^θ_(±)(m_l)=a(1+qy)/a^2y-1<0; 
          –the latitudinal motion is allowed for values of the parameter X between some local minimum X^θ_min(-)=X^θ_-(m_min(-); q, y, a) and maximum X^θ_max(+)=X^θ_+(m_max(+); q, y, a), for which
                
            X^θ_min(-)<-a<X^θ_max(+)<0;
        
                the loci m_min(-) of minimum X^θ_min(-) is given by the equation (<ref>), the loci m_max(+) of maximum X^θ_max(+) is determined by relation (<ref>);
              
          – if X takes one of these extremal values, then the trajectory of such photon lies entirely on cones θ=arccos√(m_ex),θ=π -arccos√(m_ex), where m_ex∈{m_min(-),m_max(+)}; such photons are called PNC photons <cit.>;
              
          – for X^θ_min(-)<X<X^θ_max(+) 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: 
            0 < arccos√( m_2)≤θ≤arccos√(m_1) < π/2
         and 
            π/2 < π  - arccos√(m_1)≤θ≤π-arccos√(m_2) < π;
        
          – 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 (<ref>) can be reduced to
            
            X^θ_±(m; y, a)=-a(1∓√((1-m)Δ_m))/Δ_m-a^2y, 
            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 θ/λ>0 are still confined to an interval with limits
                  
                  X^θ_min(-)   =   X^θ_-(m=0; q=0, y, a) 
           =   2a/(a^2y-1)
        
                  X^θ_max(+)   =   X^θ_+(m=0; q=0, y, a)=0, 
                  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 θ/λ=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 θ/λ=0 is now in the equatorial plane, however the equality ^2θ/λ ^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√(m_0)≤θ≤π-arccos√(m_0), 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;
            
          *  Case y=1/a^2
          –the potentials simplify into the form
            
            X^θ_±(m; q, a)=-a±√((1-m^2)(q+a^2))/m;
          * q<-a^2
          –the potentials are not defined, thus the latitudinal motion is not allowed;
            
          * q=-a^2
          –the curves X=X^θ_±(m; q=-a^2, y=1/a^2, a) coalesce, since
            
            X^θ_+(m;q=-a^2,y=1/a^2,a)   =   
        
            X^θ_-(m;q=-a^2,y=1/a^2,a)   =    X^θ_(±)(m; a)≡-a/m;
          –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 cones θ=arccos√(m_(±)),θ=π-arccos√(m_(±));
          –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 (<ref>) of the function X^θ_+(m; q, y, a) and no extremum of X^θ_-(m; q, y, a); it holds X^θ_-(m; q, y, a)<X^θ_+(m; q, y, a)<0 and X^θ_-(m; q, y, a),X^θ_+(m; q, y, a)→ -∞ as m → 0 from the right;
            
          –if X<-a or -a<X<X^θ_max(+), 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,;
          *  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 (<ref>);
                            
          –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 (<ref>), 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 (<ref>) (see the blue curve in Fig. 1d);
                            
          –if X<X^θ_max(+) or X>X^θ_min(-) photon executes vortical motion, cases X=X^θ_max(+), X=X^θ_min(-) correspond to PNC photons;
                            
          –for X=X^θ_-(m_u)=X^θ_+(m_u)=a(1+qy)/(a^2y-1), the inner cones delimiting the vortical motion are the narrowest;
                            
          –for X→ -∞ or X→ +∞ the outer cones given by angles 
            θ=arccos√(m_1),  θ=π-arccos√(m_1)
         approach the equatorial plane since m_1→ 0; for the inner cones 
            θ=arccos√(m_2),  θ=π-arccos√(m_2),
         there is 
            m_2→ m_d=1-1/a^2y;
        
          * 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 (<ref>);
                              
          –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 (<ref>) 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 (<ref>), 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;
                        
          * q>0 (Fig.2m, 2n)
             
          –the function X^θ_+(m; q, y, a) is defined at ⟨ m_l,1⟩, the function X^θ_-(m; q, y, a) at
                          ⟨ m_l,1⟩∖{m_d}, where m_l is given by (<ref>) 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
                          
            θ=arccos√(m_l),  θ=π-arccos√(m_l);
        
          –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 (<ref>) 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
         
        q≥ q_min(X, y, a),
        where q_min(X, y, a) is defined using functions
        
        q_1(X)≡ -X^2 
        and 
        
        q_2(X; y, a)≡ -I^-2[(1-a^2y)X+2a]^2
        as follows (see Fig. 3 and Fig. 4):
        
          *  Case y<1/a^2 (Fig. 3a)
          
            q_min(X, y, a) ≡{[  0,     X < 2 a/a^2 y - 1     X > 0;;                                     ; q_2(X; y, a),    2 a/a^2y-1≤ X < -a;;                                     ;            q_1(X),     -a ≤  X ≤  0;;                                      ].
        
          *  Case y=1/a^2 (Fig. 3b)
          
            q_min(X, y, a) ≡{[       -a^2,     X ≤ -a;;                        ; q_1(X),     -a ≤ X ≤ 0;;                        ;            0,     X≥ 0; ].
        
          *  Case y>1/a^2 (Fig. 3c)
          
            q_min(X, y, a) ≡{  [     q_2(X; y, a),     X < -a;                             ;           2 a/a^2 y - 1 < X;;                             ;      q_1(X),     -a ≤ X ≤ 0;;                             ; 0,     0 < X ≤2 a/a^2 y - 1; ].
        
        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. <ref> 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.
        
        
           
         
        
        
        § RADIAL MOTION
        
        
         From the equation (<ref>) 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
         
         X≤ X_-  or   X≥ X_+, 
         if
         
            a^2 - Δ_r > 0   (and   X_- < X_+),
        
         or
         
         X_+≤ X ≤ X_-  if   a^2 - Δ_r < 0, 
         where
         
         X_±(r;q, y, a)≡ar^2±√(Δ_r[r^4+q(a^2-Δ_r)])/a^2-Δ_r.
        
         We start the analysis by determining the reality region of the effective potential X_±. From the expression (<ref>) it follows that this function is well defined for
          q {[ ≤ q_r(r; y, a^2)     a^2 < Δ_r     Δ_r ≤ 0;          ≥  q_r(r; y, a^2)     0≤Δ_r < a^2 ].,
        where we have introduced the reality function
        
        q_r(r; y, a^2)≡r^4/Δ_r-a^2= r^3/r-2-yr(r^2+a^2).
        
        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 <ref>,  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  X_+(r_h) = X_-(r_h) = r_h^2/a (c. f.  <cit.>). 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  X_+(r_q) = X_-(r_q) = -a q/r_q^2.
        
        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
        Δ_r=a^2,  
        which one can express by the relation
        
        y=y_d(r; a^2)≡r-2/r(r^2+a^2). 
        The function X_-(r; q, y, a) can not diverge at radii r_d given by (<ref>), since using an alternative expression
        
        X_±(r; q, y, a)≡r^4-qΔ_r/ar^2∓√(Δ_r[r^4+q(a^2-Δ_r)])
        it can be shown that it has finite value 
        
        X_-(r_d; q, y, a)=r_d^4-qa^2/2ar_d^2.
        
         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 <cit.>, 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
        
        a^2=a^2_max(d)(r)≡ r^2(r-3)..
        
         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
        
        y = y_ex(r)(r;  a^2) ≡r-3/a^2r.
        
        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
        
        r^4+a^2q=0,
        
        or
              qa^2[2yr^3+(ya^2-1)r+1]^2+
              r^3[y^2a^4r^3+2ya^2r^2(r+3)+ r(r-3)^2-4a^2]=0.
        
        
        This can be rewritten in terms of parameter q as
        
        q=q_ex1(r; a^2)≡ -r^4/a^2, 
        and
        
        q   =   q_ex(r; y, a^2)
           ≡    -r^3/a^2 y^2a^4r^3+2ya^2r^2(r+3)+r(r-3)^2-4a^2/[2yr^3+(ya^2-1)r+1]^2.  
        
         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
        
        y=r^2-2r+a^2/r^2(r^2+a^2)=y_h(r; a^2),
        
        and
        
        y=1/r^3,
        
        i. e., they are located at event horizons and so called static radius r_s=1/√(y), where the gravitational attraction is just compensated by cosmological repulsion <cit.>. 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
        
        y=y_d(ex)(r; a^2)≡r-1/r(2r^2+a^2),
        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
        
        a^2=a^2_max(d(ex))(r)≡ 2r^2(2r-3),
        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
        
        y=y_z(ex)±(r; a^2)≡-r(r+3)± 2√(r(3r^2+a^2))/a^2r^2.
        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 <cit.>
        a^2=a^2_z(z(ex)+)(r)≡r(r-3)^2/4.
        The extrema of the function y_z(ex)+(r; a^2) are determined by the equation
        
        a^2=a^2_ex(z(ex)+)±(r)≡r(1-2r±√(1+8r))/2=a^2_ex(h)±(r),
        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
        
        y=y_ex(ex)(r; a^2)≡r-3/a^2r=y_ex(r)(r;  a^2) 
        and
        
            y=y_ex(ex)±(r; a^2)≡
        3r^2√(r)-a^2√(r)(3+2r)±√((4a^2-3r)(a^4+6a^2r^2-3r^4))/2a^4√(r^3).
        
        Using the relation (<ref>), 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
        
        y=y_h(r; a^2),
        
        i. e. they are located at the event horizons.
        
        The reality conditions of the functions y_ex(ex)±(r; a^2) read
        
        a^2≤ a^2_r(ex(ex±))+(r)     a^2_r(ex(ex±))(r)≤ a^2, 
            0<r≤r̂,
        
        and
        
        a^2≤ a^2_r(ex(ex±))(r)     a^2_r(ex(ex±))+(r)≤ a^2,
        
            r̂≤ r,
        
        where
         a^2_r(ex(ex)±)+(r)≡ (+2√(3)-3)r^2
        and
         a^2_r(ex(ex)±)(r)≡3/4r. 
        The marginal radius r̂ has value
        r̂=2√(3)+3/4=1.61603=r_crit
        and it holds
         a^2_r(ex(ex)±)+(r̂)=a^2_r(ex(ex)±)(r̂)=a^2_crit=1.21202, 
        where a^2_crit corresponds to local maximum of the function a^2_ex(h)+(r) (see e. g. <cit.> 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
        
        a^2=a^2_z(ex(ex)±)(r)≡ r(r^2-3r+3).
        The condition for stationary points 
            ∂ y_ex(ex)±(r; a^2)/∂ r=0
         leads to
        
            a^4+a^2r(2r-1)+r^3(r-3)=0,
        
        which can be solved with respect to a^2 with the same result as given by (<ref>). 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
        
        a^2_inf(ex(ex)±)+(r)≡ a^2_ex(h)+(r)
        
        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:
        1/a^2>y_ex(ex)(r; a^2)>y_ex(ex)-(r; a^2)>y_h(r; a^2)>y_d(r; a^2)>y_d(ex)(r; a^2)>0>y_z(ex)(r; a^2)>-1/a^2 as r→∞.
        
        In  Fig.5 we present all the characteristic functions related to spin parameter governing the effective potential on the lowest level:
        
          * a^2_z(h)(r)
          * a^2_ex(h)+(r)=a^2_ex(z(ex))+(r)=a^2_inf(ex(ex)±)+(r)
          * a^2_max(d)(r)
          * 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 y- values 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^2y)-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(h)(a^2)=y_max(z(ex)+)(a^2)=y_inf(ex(ex)-)(a^2)=
        	y_d(ex)h(z(ex)+)ex(ex)-(a^2)
          * y_min(h)(a^2)=y_min(z(ex)+)(a^2)=y_inf(ex(ex)+)(a^2)=
        	y_d(ex)h(z(ex)+)ex(ex)+(a^2)
          * y_max(d)(a^2)=y_dd(ex)(z(ex)+)ex(ex)(a^2)
          * y_d(z(ex)+)(a^2)
          * 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^2y)-plane and they are demonstrated in Fig. 7.
        
        
        
        The functions y_ex(h)(a^2) divide the parameter plane (a^2y) into regions describing  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 <cit.> 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));
          * y_dz(ex)+(a^2) is determined by 
        	
            a^2_dz(ex)+(r)≡r/8(1-4r+√(40r+1))
        
        	and y_d(r; a^2=a^2_dz(ex)+(r)), where the function a^2_dz(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 
        	
            a^2_ex(ex) ex(ex)±(r) ≡
        
            r/2(4r^2-12r+3±√(16r^4-96r^3+156r^2-36r+9))
        
        	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 (<ref>) imposed on the latitudinal motion yields stronger constraints on the value X than that given by the relations (<ref>), (<ref>), (<ref>) 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
        
        q ≤ q_max(r; X, a, y)≡(r^2-a X)^2/Δ_r-X^2 
        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 (<ref>) - (<ref>). 
        
        We therefore solve two equations - q_max(r; X, y, a^2)=q_1(X), with the result
        
        X=r^2/a, 
        and q_max(r;X, a, y)=q_2(X; a, y), which gives
        
        X_1,2=I^2r^2+2(a^2y-1)Δ_r∓ 2I√(-2rΔ_r)/a(I^2-4yΔ_r). 
        The solution (<ref>)  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 (<ref>)-(<ref>). The solutions (<ref>) 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 (<ref>); photons satisfying the relation (<ref>) 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
        
          I:y ≤ y_dz(ex)+(a^2) for a^2≤ 0.5
          II:y_dz(ex)+(a^2)≤ y≤ y_max(d)(a^2) for a^2≤ 0.5,
        
        	or
        y≤ y_max(d)(a^2) for 0.5 ≤ a^2≤ 1,
        
        	or
        y_min(h)(a^2)≤ y≤ y_max(d)(a^2) for 1≤ a^2≤ 1.08316;
          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;
          VIII:y≥ 1/a^2;
        
        
        Now it remains to assign to each region of the (a^2y)-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 (<ref>)-(<ref>), for 1.125≤ a^2 ≤ a^2_crit and 
         y≤ y_(ex(ex)+)(ex(ex)-)1(a^2)≡8a^2-9/8a^4 , 
        or a^2_crit≤ a^2 ≤ 1.3923 and
         y≤ y_(ex(ex)+)(ex(ex)-)2(a^2)≡√(3(2√(3)-3)/a^6)-1/a^2 , 
        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 (<ref>), (<ref>) 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^2y)-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).
           
         
        
        
        § SPHERICAL PHOTON ORBITS AND CLASSIFICATION OF THE  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  spacetimes with the spacetime parameters chosen from different parts of the (a^2y)-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  spacetimes due to the photon motion (null geodesics). Similarly as in <cit.>, 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. 
        
        
         §.§ 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 
         
           X =X_sph(r)≡r[(1-a^2y)r^2-3r+2a^2]/a[yr(2r^2+a^2)-r+1], 
        
           q = q_sph(r)≡ q_ex(r; y, a^2),
        
        where the function q_ex(r; y, a^2) is defined by the relation (<ref>).
        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 
        d^2R/d^2r = 12r^2[1+y(X^2+q)]+2(X^2+q)(a^2y-1)-4aX, 
        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 <cit.>.
        The LNRF tetrad of differential one-forms is given by the relations
        
        ω^(t)   =   √(Δ_rΔ_θρ^2/I^2 A) t, 
        ω^(r)   =   √(ρ^2/Δ_r) r,
        ω^(θ)   =   √(ρ^2/Δ_θ)θ,
        ω^(ϕ)   =   √(A sin ^2 θ/I^2 ρ^2)(ϕ-Ω_LNRF t), 
        
        the corresponding tetrad of dual vectors reads 
        
             e_(t)   =   √(I^2A/Δ_r Δ_θρ^2)( ∂/∂ t+Ω_LNRF∂/∂ϕ) ,
        
             e_(r)   =   √(Δ_r/ρ^2)∂/∂ r,
        
             e_(θ)   =   √(Δ_θ/ρ^2)∂/∂θ,
        
             e_(ϕ)   =   √(I^2 ρ^2/A sin^2θ)∂/∂ϕ. 
        The wave-vector components related to the LNRFs are then determined by the relations 
         k^(a)= ω^(a)_μk^μ,    k_(b)=k_ν e_(b)^ν, 
        hence
        
        k^(t)   =   I E √(A/Δ_r Δ_θρ^2) (1-Ω_LNRF(X+a)),
         
        k^(r)   =   ±I E/√(Δ_rρ^2)[(r^2-aX)^2-Δ_r(X^2+q)],
         
        k^(ϕ)   =   I E √(ρ^2/A sin^2θ)(X+a), 
        
        k^(θ)   =   ±I E/√(Δ_θρ^2)×
              √((X^2+q)Δ_θ-(acos^2θ+X)^2/sin^2θ),
        
        where
        
        A = (r^2+a^2)^2-a^2Δ_rsin^2θ,
        
        and
        Ω_LNRF = a[(r^2+a^2)Δ_θ-Δ_r]/A
        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 
        
        sinΨ = ρ^2√(Δ_rΔ_θ)/A sinθX+a/1-Ω_LNRF(X+a).
         
        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 co-rotating (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. [However, we have to note that, similarly to the case of Kerr black holes <cit.>, the sign of the variation of the azimuthal coordinate can be changed at some latitude, if related to distant observers.] Since the functions A, Ω_LNRF are positive <cit.>, it is clear that all photons with X<-a (ℓ<0) are retrograde. However, photons with 
        
        X>1/Ω_LNRF-a    (ℓ>1/Ω_LNRF>0) 
        can be retrograde as well. Considering in such a case the relation for the tetrad LNRF component (<ref>), 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.[Keeping E>0 means k^(t)<0, i.e., a photon in negative-root state with time evolution directed to past – for details see <cit.>.] In order to find conditions under which such a situation occurs, it is convenient to express from the alternate relation  
        sinΨ = ρ^2√(Δ_rΔ_θ)/A sinθℓ/1-Ω_LNRFℓ 
        the impact parameter ℓ in the form
        ℓ=Asinψsinθ/√(Δ_rΔ_θ)ρ^2+AΩ_LNRFsinψsinθ 
        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 
        
        -1≤sinψ < -ρ^2 √(Δ_rΔ_θ)/A Ωsinθ,
        
        from which it follows
        
        AΩsinθ - ρ^2 √(Δ_rΔ_θ)≥0.
        	
        However, the last inequality can be written in the form
        
         A I^2ρ^2 g_tt≥0,
        
        which implies 
        
        g_tt≥0.
        
        Hence, such a situation can occur only in the ergosphere. Of course, the impact parameter of such photons must fulfil the condition (<ref>). The function
        1/Ω-a=r^2(r^2+a^2)Δ_θ+a^2Δ_rcosθ^2/a[(r^2+a^2)Δ_θ-Δ_r]
        has for Δ_r=0 common points with the potentials X_± given by the relation (<ref>). There are no other intersections with the potentials, hence, the reality condition of the radial motion together with (<ref>) 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
        
        M(m; X_sph(r_0), q_sph(r_0), y, a)=0, 
        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
        θ_min=arccos√(m_0),   θ_max=π-arccos√(m_0) ,
        
        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. 
        
        
        
         §.§ Classification
        
        In the following classification we introduce ten classes of the  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.<ref> 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^2cos^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 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. <ref>a). 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. <ref>b). 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 <cit.>.
        	 
        		
          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 <ref> 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. <ref>d). 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. <ref>e). 
                
                
          Class V:
                Naked singularity spacetimes having the structure of the spherical orbits corresponding to the previous case (Fig. <ref>f), 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. <ref>h). 
        
                
          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. 
        
        
        § CONCLUSIONS
        
        
        We can summarize our results by the following concluding remarks.
        
        	
          *  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.
        		 
        		
          *  No photons with q>0 can reach the ring singularity at r=0 in any of the  spacetimes.
        		
        		
          *  In the  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 (<ref>) - (<ref>). 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.
        		
        		
          *  In the  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 (<ref>.) 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 (<ref>). The consequences for photon motion are then the same as in previous note (Figs 10(γ)-(δ)).
        		
        		
          *  In the  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)).
                
                
          *  In the  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.
                
                
          *  In the   spacetimes of classes IVb, V ,VIb, VII, VIII there are no prograde or polar spherical orbits. 
                
                
          *  In each class of the  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.
        
                
          *  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 <cit.>, containing a special type of equatorial stable photon orbits. 
        
        
        
        
        
        
        § ACKNOWLEDGMENTS
        
        
        Z.S. acknowledges the Albert Einstein Centre for Gravitation and Astrophysics supported by the Czech
        Science Foundation Grant No. 14-37086G. D.Ch. acknowledges the Silesian University in Opava Grant No. SGS/14/2016.
        
        
        abbrv10Ada-etal:2013:ASTRA:DifLiYoClG
        	C. Adami, F. Durret, L. Guennou, and C. Da Rocha.
        	Diffuse light in the young cluster of galaxies CL J1449+0856 at
        		z=2.07.
        	 Astronomy and Astrophysics, 551:A20 (7 pages), Mar. 2013.
        	
        	Ali:2007:PHYSR4:EMPropKadS
        	A. N. Aliev.
        	Electromagnetic properties of Kerr–anti-de Sitter black holes.
        	 Phys. Rev. D, 75(8):084041, Apr. 2007.
        	
        	ArP-Muk-Ste:2000:PHYRL:
        	C. Armendariz-Picon, V. Mukhanov, and P. J. Steinhardt.
        	 Dynamical solution to the problem of a small cosmological constant
        	and late-time cosmic acceleration.
        	 Phys. Rev. Lett., 85(21):4438, 2000.
        	
        	Arra:2014:PHYSR4:
        	I. Arraut.
        	Komar mass function in the de Rham-Gabadadze-Tolley nonlinear theory
        		of massive gravity.
        	 Phys. Rev. D, 90:124082, Dec 2014.
        	
        	Arra:2017:Universe:
        	I. Arraut.
        	 The astrophysical scales set by the cosmological constant, black-hole
        	thermodynamics and non-linear massive gravity.
        	 Universe, 3(2), 2017.
        	
        	Asc:2008:CHIAA:MassSpinBHQPO
        	B. Aschenbach.
        	Measurement of Mass and Spin of Black Holes with QPOs.
        	 Chinese Journal of Astronomy and Astrophysics Supplement,
        	8:291–296, Oct. 2008.
        	
        	Bah-etal:1999:SCIEN:
        	N. Bahcall, J. P. Ostriker, S. Perlmutter, and P. J. Steinhardt.
        	 The cosmic triangle: Revealing the state of the universe.
        	 Science, 284:1481–1488, 1999.
        	
        	Bak-etal:2007:CEURJP:
        	P. Bakala, P. Čermák, S. Hledík, Z. Stuchlík, and
        	K. Truparová.
        	Extreme gravitational lensing in vicinity of
        		Schwarz­schild–de Sitter black holes.
        	 Central European J. Phys., 5(4):599–610, Dec. 2007.
        	
        	Bardeen:1973
        	J. M. Bardeen.
        	Timelike and null geodesics in the Kerr metric. In C. Dewitt and B. S. Dewitt, editors,  Black Holes (Les
        		Astres Occlus), pages 215–239, 1973.
        	
        	BS76:
        	J. Bičák and Z. Stuchlík.
        	 On the latitudinal and radial motion in the field of a rotating black
        	hole.
        	 Bull. Astronom. Inst. Czechoslovakia, 27(3):129–133, 1976.
        	
        	Bic-Stu-Bal:1989:BAC:
        	J. Bičák, Z. Stuchlík, and V. Balek.
        	The motion of charged particles in the field of rotating charged
        		black holes and naked singularities.
        	 Bulletin of the Astronomical Institutes of Czechoslovakia,
        	40:65–92, Mar. 1989.
        	
        	Bla-Stu:2016:PHYSR4:
        	M. Blaschke and Z. Stuchlík.
        	 Efficiency of the keplerian accretion in braneworld kerr-newman
        	spacetimes and mining instability of some naked singularity spacetimes.
        	 Phys. Rev. D, 94:086006, Oct 2016.
        	
        	Boh:2004:GENRG2:
        	C. G. Böhmer.
        	Eleven Spherically Symmetric Constant Density Solutions with
        		Cosmological Constant.
        	 Gen. Relativity Gravitation, 36:1039–1054, May 2004.
        	
        	Boy-etal:2003:PHYSR4:HoloProtChron
        	E. K. Boyda, S. Ganguli, P. Hořava, and U. Varadarajan.
        	Holographic Protection of Chronology in Universes of the Gödel
        		Type.
        	 Phys. Rev. D, 67:106003, 2003.
        	
        	Cal-Kam:2009:NATURE:CosDarkMat
        	R. Caldwell and M. Kamionkowski.
        	Cosmology: Dark matter and dark energy.
        	 Nature, 458(7238):587–589, Apr. 2009.
        	
        	Cal-Dav-Ste:1998:PHYRL:
        	R. R. Caldwell, R. Dave, and P. J. Steinhardt.
        	 Cosmological imprint of an energy component with general equation of
        	state.
        	 Phys. Rev. Lett., 80(8):1582, 1998.
        	
        	Car:1973:BlaHol:
        	B. Carter.
        	Black hole equilibrium states. In C. Dewitt and B. S. Dewitt, editors,  Black Holes (Les
        		Astres Occlus), pages 57–214, 1973.
        	
        	l_e_cones:
        	D. Charbulák and Z. Stuchlík.
        	Light escape cones in local reference frames of Kerr- de Sitter
        		black hole spacetimes and related black hole shadows,(to be published in The
        		European Physical Journal C).
        	
        	Che:2008:CHINPB:DkEnGeoMorSchw
        	J.-H. Chen and Y.-J. Wang.
        	Influence of dark energy on time-like geodesic motion in
        		Schwarzschild spacetime.
        	 Chinese Physics B, 17(4):1184, 2008.
        	
        	Cru-Oli-Vil:2005:CLAQG:GeoSdSBH
        	N. Cruz, M. Olivares, and J. R. Villanueva.
        	The geodesic structure of the Schwarzschild anti-de Sitter black
        		hole.
        	 Classical Quantum Gravity, 22(6):1167–1190, Mar. 2005.
        	
        	deFel:1974:ASTRA:
        	F. de Felice.
        	 Repulsive phenomena and energy emission in the field of a naked
        	singularity.
        	 Astronomy and Astrophysics, 34:15–19, 1974.
        	
        	deFel:1978:NATURE:InstabNS
        	F. de Felice.
        	Classical instability of a naked singularity.
        	 Nature, 273:429–431, June 1978.
        	
        	Far:2016:PDU:
        	V. Faraoni.
        	Turnaround radius in modified gravity.
        	 Physics of the Dark Universe, 11:11–15, Mar. 2016.
        	
        	Far-Lap-Pra:2015:JCAP:
        	V. Faraoni, M. Lapierre-Léonard, and A. Prain.
        	Turnaround radius in an accelerated universe with quasi-local mass.
        	 Journal of Cosmology and Astroparticle Physics,
        	2015(10):013–013, Oct. 2015.
        	
        	Gib-Haw:1977:PHYSR4:
        	G. W. Gibbons and S. W. Hawking.
        	 Cosmological event horizons, thermodynamics, and particle creation.
        	 Phys. Rev. D, 15:2738–2751, May 1977.
        	
        	Gim-Hor:2004:hep-th0405019:GodHolo
        	E. G. Gimon and P. Hořava.
        	Over-Rotating Black Holes, Gödel Holography and the Hypertube,
        	2004.
        	
        	Gim-Hor:2009:PHYLB:AstVioSignStr
        	E. G. Gimon and P. Hořava.
        	Astrophysical Violations of the Kerr Bound as a Possible Signature
        		of String Theory.
        	 Phys. Lett. B, 672:299, 2009.
        	
        	Gu-Cheng:2007:GENRG2:CircLoopKdS
        	Z. Gu and H. Cheng.
        	The circular loop equation of a cosmic string in Kerr–de Sitter
        		spacetimes.
        	 Gen. Relativity Gravitation, 39(1):1–7, Jan. 2007.
        	
        	Hac-etal:2010:PHYSR4:KerrBHCoStr:
        	E. Hackmann, B. Hartmann, C. Lämmerzahl, and P. Sirimachan.
        	Test particle motion in the space-time of a Kerr black hole pierced
        		by a cosmic string.
        	 Phys. Rev. D, 82(4):044024, Aug. 2010.
        	
        	Hio-Mae:2009:PHYSR4:KerrSpinMeas
        	K. Hioki and K.-i. Maeda.
        	Measurement of the Kerr spin parameter by observation of a compact
        		object's shadow.
        	 Phys. Rev. D, 80(2):024042 (9 pages), July 2009.
        	
        	Ior:2009:NEWASTR:CCDGPGrav
        	L. Iorio.
        	Constraining the cosmological constant and the DGP gravity with the
        		double pulsar PSR J0737-3039.
        	 New Astronomy, 14(2):196–199, Feb. 2009.
        	
        	Kag-Kun-Lam:2006:PHYLB:SolarSdS
        	V. Kagramanova, J. Kunz, and C. Lammerzahl.
        	Solar system effects in Schwarzschild–de Sitter space-time.
        	 Phys. Lett. B, 634(5–6):465–470, Mar. 2006.
        	
        	Kol-Stu:2010:PHYSR4:CurCarStrLoops
        	M. Kološ and Z. Stuchlík.
        	Current-carrying string loops in black-hole spacetimes with a
        		repulsive cosmological constant.
        	 Phys. Rev. D, 82(12):125012 (21 pages), Dec. 2010.
        	
        	Kot:1918:ANNPH2:PhyBasEinsGr
        	F. Kottler.
        	Über die physikalischen Grundlagen der Einsteinschen
        		Gravitationstheorie.
        	 Annalen der Physik, 361(14):401–462, 1918.
        	
        	Kra:2005:DARK:CCPerPrec
        	G. V. Kraniotis.
        	Precise theory of orbits in general relativity, the cosmological
        		constant and the perihelion precession of Mercury.
        	 pages 469–479.
        	
        	Kra:2004:CLAQG:
        	G. V. Kraniotis.
        	Precise relativistic orbits in Kerr and Kerr–(anti-)de Sitter
        		spacetimes.
        	 Classical Quantum Gravity, 21:4743–4769, 2004.
        	
        	Kra:2007:CLAQG:Periapsis
        	G. V. Kraniotis.
        	Periapsis and gravitomagnetic precessions of stellar orbits in Kerr
        		and Kerr–de Sitter black hole spacetimes.
        	 Classical Quantum Gravity, 24:1775–1808, 2007.
        	
        	Kra:2011:CLAQG:
        	G. V. Kraniotis.
        	Precise analytic treatment of Kerr and Kerr-(anti) de Sitter black
        		holes as gravitational lenses.
        	 Class. Quant. Grav., 28:085021, 2011.
        	
        	Kra:2014:GRG:
        	G. V. Kraniotis.
        	Gravitational lensing and frame dragging of light in the Kerr-Newman
        		and the Kerr-Newman-(anti) de Sitter black hole spacetimes.
        	 Gen. Rel. Grav., 46(11):1818, 2014.
        	
        	Kra:1998:ASTRJ2:
        	L. M. Krauss.
        	 The end of the age problem, and the case for a cosmological constant
        	revisited.
        	 Astrophys. J., 501(2):461–466, 1998.
        	
        	Kra-Tur:1995:GENRG2:
        	L. M. Krauss and M. S. Turner.
        	The cosmological constant is back.
        	 Gen. Relativity Gravitation, 27(11):1137–1144, Nov. 1995.
        	
        	Kuc-Sla-Stu:2011:JCAP:ToroPerFlRNadS:
        	H. Kučáková, P. Slaný, and Stuchlík.
        	 Toroidal configurations of perfect fluid in the
        	reissner-nordström-(anti-)de sitter spacetimes.
        	 Journal of Cosmology and Astroparticle Physics, 2011(01):033,
        	2011.
        	
        	Lak:2002:PHYSR4:BendLiCC
        	K. Lake.
        	Bending of light and the cosmological constant.
        	 Phys. Rev. D, 65(8, B):087301, Apr. 2002.
        	
        	Lak-Zan:2016:PHYSR4:
        	K. Lake and T. Zannias.
        	Global structure of Kerr -de Sitter spacetimes.
        	 Phys. Rev. D, 92:084003, Oct 2015.
        	
        	Lin:1990:InfCos:
        	A. D. Linde.
        	 Particle Physics and Inflationary Cosmology.
        	 Gordon and Breach, New York, 1990.
        	
        	Mul-Asch:2007:CLAQG:NonMonoVel
        	A. Müller and B. Aschenbach.
        	Non-monotonic orbital velocity profiles around rapidly rotating Kerr
        		(anti-)de Sitter black holes.
        	 Classical and Quantum Gravity, 24:2637–2644, May 2007.
        	
        	Mul:2008:GENRG2:FallSchBH
        	T. Müller.
        	Falling into a Schwarzschild black hole.
        	 Gen. Relativity Gravitation, pages 56–+, Feb. 2008.
        	
        	Oli-etal:2011:MODPLA:ChaParRNadS:
        	M. Olivares, J. Saavedra, C. Leiva, and J. R. Villanueva.
        	Motion of charged particles on the Reissner–Nordström
        		(anti)–de Sitter black hole spacetime.
        	 Modern Phys. Lett. A, 26(39):2923–2950, Dec. 2011.
        	
        	Ost-Ste:1995:NATURE:
        	J. P. Ostriker and P. J. Steinhardt.
        	The observational case for a low-density universe with a nonzero
        		cosmological constant.
        	 Nature, 377(6550):600–602, Oct. 1995.
        	
        	Per-Rom-PeB:2013:ASTRA:AccDiBHModGra:
        	D. Pérez, G. E. Romero, and S. E. Bergliaffa.
        	Accretion discs around black holes in modified strong gravity.
        	 Astronomy and Astrophysics, 551:A4(15 pages), Mar. 2013.
        	
        	Cha-Har:2012:PHYSR4:BEConGRStar
        	C. Pierre-Henri and T. Harko.
        	Bose–Einstein Condensate general relativistic star.
        	 Phys. Rev. D, 86(6):064011, Sept. 2012.
        	
        	Ade-etal:2014:ASTRA:Planck Collaboration, P. A. R. Ade, N. Aghanim, C. Armitage-Caplan,
        	M. Arnaud, M. Ashdown, F. Atrio-Barandela, J. Aumont,
        	C. Baccigalupi, A. J. Banday, and et al.
        	Planck 2013 results. XII. Diffuse component separation.
        	 Astronomy and Astrophysics, 571:A12, Nov. 2014.
        	
        	Pug-Stu:2015:ApJS:
        	D. Pugliese and Z. Stuchlík.
        	Ringed Accretion Disks: Equilibrium Configurations.
        	 Astrophys. J. Suppl., 221(2):25, dec 2015.
        	
        	Rez-Zan-Fon:2003:ASTRA:
        	L. Rezzolla, O. Zanotti, and J. A. Font.
        	Dynamics of thick discs around Schwarzschild–de Sitter black
        		holes.
        	 Astronomy and Astrophysics, 412(3):603–613, Dec. 2003.
        	
        	Rie-etal:2004:ASTRJ2:
        	A. G. Riess et al.
        	Type Ia Supernova Discoveries at z>1 From the Hubble Space
        		Telescope: Evidence for Past Deceleration and Constraints on Dark Energy
        		Evolution.
        	 Astrophys. J., 123:145, 2004.
        	
        	Sche-Stu-Pet:2013:JCAP:
        	J. Schee, Z. Stuchlík, and M. Petrásek.
        	Influence of the cosmic repulsion on the MOND model of the
        		Magellanic Cloud motion in the field of Milky Way.
        	 Journal of Cosmology and Astroparticle Physics, 12:026, Dec.
        	2013.
        	
        	Sch-Zai:2008:0801.3776:CCTimeDelay
        	T. Schücker and N. Zaimen.
        	Cosmological constant and time delay.
        	 Astronomy and Astrophysics, 484(1):103–106, June 2008.
        	
        	Ser:2008:PHYSR4:CCLens
        	M. Sereno.
        	On the influence of the cosmological constant on gravitational
        		lensing in small systems.
        	 Phys. Rev. D, 77(4):043004, 2008.
        	
        	Zhou-Chen:2011:AstrSpSc:
        	Z. Sheng, C. Ju-Hua, and W. Yong-Jiu.
        	Time-like geodesic structure of a spherically symmetric black hole
        		in the brane-world.
        	 Chinese Physics B, 20(10):100401, 2011.
        	
        	Sla-Stu:2005:CLAQG:
        	P. Slaný and Z. Stuchlík.
        	Relativistic thick discs in the Kerr–de Sitter backgrounds.
        	 Classical Quantum Gravity, 22(17):3623–3651, 2005.
        	
        	Sla-Stu:2008:CLAQG:CmtNoMonKadS
        	P. Slaný and Z. Stuchlík.
        	 Comment on 'non-monotonic orbital velocity profiles around rapidly
        	rotating kerr–(anti-)de sitter black holes'.
        	 Classical and Quantum Gravity, 25(3):038001, 2008.
        	
        	Spe-etal:2007:ASTJS:3yrWMAP
        	D. N. Spergel, R. Bean, O. Dore, M. R. Nolta, C. L. Bennett, J. Dunkley,
        	G. Hinshaw, N. Jarosik, E. Komatsu, L. Page, H. V. Peiris, L. Verde,
        	M. Halpern, R. S. Hill, A. Kogut, M. Limon, S. S. Meyer, N. Odegard, G. S.
        	Tucker, J. L. Weiland, E. Wollack, and E. L. Wright.
        	Three-year Wilkinson Microwave Anisotropy Probe (WMAP) observations:
        		Implications for cosmology.
        	 Astrophys. J. Suppl., 170(2):377–408, June 2007.
        	
        	Spe-etal:2007:ApJSuppl:
        	D. N. Spergel, R. Bean, O. Dore, M. R. Nolta, C. L. Bennett, J. Dunkley,
        	G. Hinshaw, N. Jarosik, E. Komatsu, L. Page, H. V. Peiris, L. Verde,
        	M. Halpern, R. S. Hill, A. Kogut, M. Limon, S. S. Meyer, N. Odegard, G. S.
        	Tucker, J. L. Weiland, E. Wollack, and E. L. Wright.
        	Three‐Year Wilkinson Microwave Anisotropy Probe ( WMAP )
        		Observations: Implications for Cosmology.
        	 Astrophys. J. Suppl., 170(2):377–408, jun 2007.
        	
        	Stu:1980:BULAI:
        	Z. Stuchlík.
        	Equatorial circular orbits and the motion of the shell of dust in
        		the field of a rotating naked singularity.
        	 Bull. Astronom. Inst. Czechoslovakia, 31:129–144, 1980.
        	
        	Stu:1981:BAIC:Rad.mot.ph.Kerr
        	Z. Stuchlík.
        	The radial motion of photons in Kerr metric.
        	 Bulletin of the Astronomical Institutes of Czechoslovakia,
        	32:40–52, 1981.
        	
        	Stu:1983:BULAI:
        	Z. Stuchlík.
        	 The motion of test particles in black-hole backgrounds with non-zero
        	cosmological constant.
        	 Bull. Astronom. Inst. Czechoslovakia, 34(3):129–149, 1983.
        	
        	Stu:1984:BULAI:
        	Z. Stuchlík.
        	 An Einstein–Strauss–de Sitter model of the universe.
        	 Bull. Astronom. Inst. Czechoslovakia, 35(4):205–215, 1984.
        	
        	Stu:2000:ACTPS2:
        	Z. Stuchlík.
        	 Spherically symmetric static configurations of uniform density in
        	spacetimes with a non-zero cosmological constant.
        	 Acta Phys. Slovaca, 50(2):219–228, Mar. 2000.
        	
        	Stu:2005:MODPLA:
        	Z. Stuchlík.
        	Influence of the Relict Cosmological Constant on Accretion Discs.
        	 Modern Phys. Lett. A, 20(8):561–575, Mar. 2005.
        	
        	Stu-Ba_Ost:1998:KNdSrest.rep.bar.epm
        	Z. Stuchlík, G. Bao, E. Østgaard, and S. Hledík.
        	Kerr-Newman-de Sitter black holes with a restricted repulsive
        		barrier of equatorial photon motion.
        	 Phys. Rev. D, 58(8):084003, Oct. 1998.
        	
        	Stu-Cal:1991:GENRG2:
        	Z. Stuchlík and M. Calvani.
        	Null geodesics in black-hole metrics with nonzero cosmological
        		constant.
        	 Gen. Relativity Gravitation, 23(5):507–519, May 1991.
        	
        	Stu-Hle:1999:PHYSR4:
        	Z. Stuchlík and S. Hledík.
        	Some properties of the Schwarz­schild–de Sitter and
        		Schwarz­schild–anti-de Sitter spacetimes.
        	 Phys. Rev. D, 60(4):044006 (15 pages), Aug. 1999.
        	
        	Stu-Hle:2000:CLAQG:
        	Z. Stuchlík and S. Hledík.
        	Equatorial photon motion in the Kerr–Newman spacetimes with a
        		non-zero cosmological constant.
        	 Classical Quantum Gravity, 17(21):4541–4576, Nov. 2000.
        	
        	Stu-Hle:2002:ACTPS2:
        	Z. Stuchlík and S. Hledík.
        	Properties of the Reissner–Nordström spacetimes with a nonzero
        		cosmological constant.
        	 Acta Phys. Slovaca, 52(5):363–407, Oct. 2002.
        	
        	Stu-Hle-Nov:2016:PHYSR4:
        	Z. Stuchlík, S. Hledík, and J. Novotný.
        	 General relativistic polytropes with a repulsive cosmological
        	constant.
        	 Phys. Rev. D, 94:103513, Nov 2016.
        	
        	Stu-Hle-Tru:2011:CLAQG:
        	Z. Stuchlík, S. Hledík, and K. Truparová.
        	Evolution of Kerr superspinars due to accretion counterrotating thin
        		discs.
        	 Classical Quantum Gravity, 28(15):155017, Aug. 2011.
        	
        	Stu-Kol:2012:PHYSR4:AccStringLoops
        	Z. Stuchlík and M. Kološ.
        	Acceleration of string loops in the Schwarzschild–de Sitter
        		geometry.
        	 Phys. Rev. D, 85(6):065022 [13 pages], 2012.
        	
        	Stu-Kol:2012:JCAP:StringLoops:
        	Z. Stuchlík and M. Kološ.
        	String loops in the field of braneworld spherically symmetric black
        		holes and naked singularities.
        	 Journal of Cosmology and Astroparticle Physics, 2012:008, 2012.
        	
        	Stu-Kov:2008:INTJMD:PsNewtSdS:
        	Z. Stuchlík and J. Kovář.
        	Pseudo-Newtonian gravitational potential for
        		Schwarzschild–de Sitter spacetimes.
        	 INTJMD, 17(11):2089–2105, 2008.
        	
        	Stu-Sch:2010:CLAQG:AppKepDiOrKerrSSp
        	Z. Stuchlík and J. Schee.
        	Appearance of Keplerian discs orbiting Kerr superspinars.
        	 Classical Quantum Gravity, 27(21):215017 (39 pages), Nov. 2010.
        	
        	Stu-Sch:2011:JCAP:CCMagOnCloud:
        	Z. Stuchlík and J. Schee.
        	Influence of the cosmological constant on the motion of Magellanic
        		Clouds in the gravitational field of Milky Way.
        	 Journal of Cosmology and Astroparticle Physics, 9:018–018,
        	Sept. 2011.
        	
        	Stu-Sch:2012:INTJMD:GRvsPsNewtMagClou:
        	Z. Stuchlík and J. Schee.
        	Comparison of general relativistic and pseudo-Newtonian description
        		of Magellanic-clouds motion in the field of Milky Way.
        	 Internat. J. Modern Phys. D, 21(4):1250031, Apr. 2012.
        	
        	0264-9381-29-6-065002
        	Z. Stuchlík and J. Schee.
        	Observational phenomena related to primordial Kerr superspinars.
        	 Classical and Quantum Gravity, 29(6):065002, 2012.
        	
        	Stu-Sch:2013:CLAQG:UHEKerrGeo
        	Z. Stuchlík and J. Schee.
        	Ultra-high-energy collisions in the superspinning Kerr geometry.
        	 Classical Quantum Gravity, 30(7):075012, Apr. 2013.
        	
        	Stu-etal:2017:JCAP:
        	Z. Stuchlík, J. Schee, B. Toshmatov, J. Hladík, and
        	J. Novotný.
        	 Gravitational instability of polytropic spheres containing region of
        	trapped null geodesics: a possible explanation of central supermassive black
        	holes in galactic halos.
        	 Journal of Cosmology and Astroparticle Physics, 2017(06):056,
        	2017.
        	
        	Stu-Sla:2004:PHYSR4:
        	Z. Stuchlík and P. Slaný.
        	Equatorial circular orbits in the Kerr–de Sitter spacetimes.
        	 Phys. Rev. D, 69:064001, 2004.
        	
        	Stu-Sla-Hle:2000:ASTRA:
        	Z. Stuchlík, P. Slaný, and S. Hledík.
        	Equilibrium configurations of perfect fluid orbiting
        		Schwarz­schild–de Sitter black holes.
        	 Astronomy and Astrophysics, 363(2):425–439, Nov. 2000.
        	
        	Stu-Sla-Kov:2009:CLAQG:PseNewSdS
        	Z. Stuchlík, P. Slaný, and J. Kovář.
        	Pseudo-Newtonian and general relativistic barotropic tori in
        		Schwarzschild–de Sitter spacetimes.
        	 Classical Quantum Gravity, 26(21):215013 (34 pp), Nov. 2009.
        	
        	Stu-etal:2005:PHYSR4:AschenUnexpTopo:
        	Z. Stuchlík, P. Slaný, G. Török, and M. A.
        	Abramowicz.
        	Aschenbach effect: Unexpected topology changes in the motion of
        		particles and fluids orbiting rapidly rotating Kerr black holes.
        	 Phys. Rev. D, 71(2):024037, Jan. 2005.
        	
        	Teo:2003:GenRelGrav:
        	E. Teo.
        	Spherical Photon Orbits Around a Kerr Black Hole.
        	 General Relativity and Gravitation, 35:1909–1926, Nov. 2003.
        	
        	Vil-etal:2013:ASTSS1:PhMoChgAdS:
        	J. R. Villanueva, J. Saavedra, M. Olivares, and N. Cruz.
        	Photons motion in charged Anti-de Sitter black holes.
        	 Astrophys. and Space Sci., 344(2):437–446, Dec. 2012.
        	
        	Wan-etal:2000:ASTRJ2:
        	L. Wang, R. R. Caldwell, J. P. Ostriker, and P. J. Steinhardt.
        	 Cosmic concordance and quintessence.
        	 Astrophys. J., 530(1):17–35, 2000.
        	
        	Wan-Che:2012:PHYLB:CirLoopPerTens
        	L. Wang and H. Cheng.
        	The evolution of circular loops of a cosmic string with periodic
        		tension.
        	 Phys. Lett. B, 713(1):59–62, 2012.
