kds short=KdS, long=Kerr-de Sitter, qnm short=QNMs, long=Quasinormal Modes,
rkds short=RKdS, long=Kerr-de Sitter Revisited, sds short=SdS,
long=Schwarzschild-de Sitter, egb short=EGB, long=Einstein-Gauss-Bonnet, egbds
short=EGBdS, long=Einstein-Gauss-Bonnet de Sitter, gr short=GR, long=General
Relativity, vlbi short=VLBI, long=Very Long Baseline Interferometry, eht
short=EHT, long=Event Horizon Telescope, sgrA short=Sgr$A^{*}$,
long=Sagittarius $A^{*}$, zamo short=ZAMO, long=Zero Angular Momentum Orbit,

# Black hole lensing in Kerr-de Sitter spacetimes

Eunice Omwoyo<EMAIL_ADDRESS>PPGCosmo, CCE, Universidade
Federal do Espírito Santo (UFES)
Av. Fernando Ferrari, 540, CEP 29.075-910, Vitória, ES, Brazil.  Humberto
Belich<EMAIL_ADDRESS>Núcleo Cosmo-ufes & Departamento de Física,
Universidade Federal do Espírito Santo (UFES)
Av. Fernando Ferrari, 540, CEP 29.075-910, Vitória, ES, Brazil.  Júlio C.
Fabris<EMAIL_ADDRESS>Núcleo Cosmo-ufes & Departamento de
Física, Universidade Federal do Espírito Santo (UFES)
Av. Fernando Ferrari, 540, CEP 29.075-910, Vitória, ES, Brazil.  National
Research Nuclear University MEPhI, Kashirskoe sh. 31, Moscow 115409, Russia
Hermano Velten  §<EMAIL_ADDRESS>Departamento de Física,
Universidade Federal de Ouro Preto (UFOP), Campus Universitário Morro do
Cruzeiro, 35.400-000, Ouro Preto, Brazil

###### Abstract

We have derived analytical solutions using Jacobi elliptic functions for bound
and nearly bound photon orbits in kds and rkds spacetimes. Leveraging our
obtained solutions, we have conducted an analytic ray-tracing in both
spacetimes. We have obtained direct images, lensing rings and photon rings for
equatorial disks considering inclined locally static observers. Images
corresponding to $n=(2,3)$ exhibit a significantly closer resemblance to the
critical curve as compared to the $n=1$ image. This highlights the remarkable
potential of these higher-order images as robust testing grounds for gr.
Furthermore in both spacetimes, we have obtained analytical solutions for the
critical parameters governing the structure of the photon ring and analyzed
these parameters in details.

##  I Introduction

It is well-established that in the vicinity of rotating black holes there
exists a photon region, also known as a photon sphere for static black holes,
which consists of orbits where photons move at a constant radial distance,
referred to as bound photon orbits. When these orbits are unstable to radial
perturbations, photons on nearby orbits that are not precisely in the photon
region/sphere may undergo numerous half orbits close to the bound orbit before
either plunging into the black hole or escaping. These orbits are known as
nearly bound photon orbits. In the image plane of an observer, the photons on
nearly bound orbits manifest as an annular feature that is located in close
proximity to the critical curve, and this feature is known as the photon ring
[[1], [2],[3]]. The term critical curve describes the locus on the image plane
that arises from the photons moving along bound orbits and its shape follows
from General Relativity (GR). Recently, the eht collaboration has provided
interferometric observations of the black holes M87 [4] and sgrA [5] at the
center of the galaxies Messier 87 and Milky Way, respectively, thus providing
valuable insights into the emission structure on horizon scales. Calculations
within the framework of gr, [3], show that the photon ring is not a single
continual entity. Rather, it is made up of an infinite sequence of self-
similar subrings converging to the critical curve. Each subring is uniquely
identified by the number of half orbits, denoted as $n=0,1,2,..\infty$, which
nearly bound photons execute before escaping or plunging into the black hole.
The $n=0$ image, also known as the primary or direct image, is formed by
photons on weakly lensed orbits and its characteristics strongly depend on the
astrophysical source profile around the black hole. On the other hand, the
$n=1$ image is referred to as the lensed image [1]. The photon ring consists
of subrings with $n\geq 2$ resulting from photons on orbits that have
undergone extreme lensing due to the strong gravity of the black hole, causing
them to execute multiple half orbits around the black hole. Probing the photon
ring is thus a gravity probe because the presence of it is due to strong
lensing and is insensitive of the astrophysical source profile’s details.

As a photon on a nearly bound orbit escapes or falls into a black hole, the
distance between it and a bound orbit grows exponentially. The Lyapunov
exponent of the respective bound orbit measures the rate of this exponential
deviation. Furthermore, when these photons are projected onto the celestial
sphere of an observer, they generate subrings whose thickness is controlled by
the Lyapunov exponent [1]. With each additional half orbit that photons on
nearly bound orbits execute, there is a corresponding change in the azimuthal
angle, which causes subsequent subrings to rotate. Moreover, the detection of
subsequent subrings by an observer occurs at different times due to the extra
half orbits executed by the photons, thereby causing a delay in their
detection. Thus the photon ring is controlled by three critical parameters,
the change in azimuthal angle $\delta$, Lyapunov exponent $\gamma$, and time
delay $\zeta$ [1]. The time delay parameter controls the arrival time of
subsequent subrings by determining how much time photons take to complete half
orbits. Meanwhile, the change in azimuthal angle parameter controls the
rotation of the subrings by specifying the amount of change in the azimuthal
angle per half orbit. Finally, the Lyapunov exponent parameter measures the
exponential demagnification of subsequent subrings. These critical parameters
are defined on the critical radii, which is the radial coordinate of the bound
photon orbits and are independent of the emitting matter around the black
hole. Thus, to decipher the photon ring alongside these critical parameters, a
comprehensive understanding of bound photon orbits is vital.

In this work, we begin with performing an analysis of bound and nearly bound
photon orbits in Kerr-de Sitter (KdS) [6] and Kerr-de Sitter Revisited (RKdS)
[7] spacetimes. The rkds is a newly proposed solution that has been obtained
through gravitational decoupling for axially symmetric systems ([8],[9],[10]).
Unlike the kds solution, the rkds solution is neither a $\Lambda$ Vacuum
solution - it does not belong to the Plebanski-Demianski class of metrics -
nor is it a constant curvature solution. The salient feature of the rkds
solution is the emergence of a warped curvature induced by rotational effects
which is absent in kds solution. This deformation is likely to be closely
related to the warped vacuum energy in the direct proximity of extreme
gravitational sources [11]. It is possible that the deformation of the vacuum
energy in the vicinity of the rkds black hole could be of relevant
astrophysical significance and thus extensive study of this solution is
necessary.

In Ref.[12], analytic solutions to geodesics in kds spacetime have been
obtained using Weierstrass’ elliptic $\wp$ function and Weierstrass’ $\zeta$
and $\sigma$ functions. The Weierstrass’ elliptic $\wp$ function and $\sigma$
function are multivalued while the Weierstrass’ $\zeta$ function is single-
valued. A multivalued function, also known as a set-valued function, is a
mapping that assigns a set of output values to each input value in its domain
while a single-valued function can be defined mathematically as a mapping that
assigns a unique output value to each input value in its domain. Nearly bound
orbits involve photons winding around the black hole multiple times before
escaping or falling into the black hole, thus, there can be multiple values of
the angular coordinate at certain radial positions. To capture this winding
behaviour explicitly, it is important to obtain solutions that are solely in
terms of multivalued functions. As a result, the first goal of this work is to
obtain solutions to null geodesic equations of kds and rkds spacetime in terms
of Jacobi elliptic functions, which are all multivalued in nature.
Furthermore, solutions obtained in terms of Jacobi elliptic functions do not
require gluing and instead extend through various regions, providing explicit
complete trajectories. The Jacobi elliptic functions approach has been
considered to fully characterize null geodesics in the exterior of a Kerr
black hole [13], so we refer the reader to this reference in the case of
vanishing cosmological constant. Using our solutions, we perform an analytic
ray-tracing to obtain direct images, lensing rings and photon rings in kds and
rkds spacetime. Besides, we also obtain analytic forms of the critical
parameters in both spacetimes.

This work is structured as follows: In section II, a comprehensive analysis of
kds and rkds spacetimes is presented, accompanied by the derivation of the
solutions to their respective null geodesic equations utilizing Jacobi
elliptic functions. This section establishes the groundwork upon which the
subsequent sections are built. Moreover, the classification and investigation
of the general structure of photon orbits, integral to this work, are examined
in this section. Section III capitalizes on the derived solutions to perform
an analytical ray-tracing in kds and rkds spacetimes. By considering inclined
locally static observers, we investigate images of equatorial disks around kds
and rkds black holes, encompassing direct images, lensing rings and photon
rings. Building upon the properties of the aforementioned images, section IV
delves into the critical parameters that govern the structure of the photon
ring. Analytical expressions for these parameters are derived and their
detailed investigation is provided. Lastly but not least, in section V we give
a conclusion of this work. The common derivation steps involved in this work
are elucidated in appendix A, (B), (C), (D) and (E).

##  II Solutions of the Equations of Motion

In this section, we commence by providing an overview of kds and rkds
spacetimes. Leveraging the general forms derived for the radial and angular
integrals, as outlined in Appendix C, we proceed to derive solutions for the
$r$, $\theta$, $\phi$, and $t$ components encompassing the trajectory of the
photons. The primary objective of this section is to unravel the general
structure characterizing bound and nearly bound photon orbits thus we keep the
solutions simple. In so doing, we make the assumption that a turning point has
yet to be reached [13]. The periodic nature of Jacobi elliptic functions
allows for this assumption. Hence, in spite of this assumption, the solutions
will still provide trajectories that go through the turning points, and the
physical effects on the turning points, such as gravitational lensing, will
still hold. Throughout our analysis, the terms ”source” and ”observer” are
frequently employed. Herein, the ”source” refers to the location of photon
emission, while the ”observer” denotes the location of the photon detector. It
is imperative to emphasize that our considerations are confined to observers
situated within the domain of outer communication.

###  II.1 Kerr-de Sitter spacetime

kds metric is a solution of Einstein’s field equations that describes a
rotating black hole in a universe with a positive cosmological constant. In
the Boyer-Lindquist coordinate system the metric takes the form [6],

$\text{ds}^{2}=\left(\frac{a^{2}\Delta_{\theta}\sin^{2}(\theta)}{L^{2}\Sigma}-\frac{\text{$\Delta_{r}$}}{L^{2}\Sigma}\right)\text{dt}^{2}-\frac{2a\sin^{2}(\theta)\left(\Delta_{\theta}\left(a^{2}+r^{2}\right)-\text{$\Delta_{r}$}\right)}{L^{2}\Sigma}\text{dt}\text{d$\phi$}\\\
+\frac{\sin^{2}(\theta)\left(\Delta_{\theta}\left(a^{2}+r^{2}\right)^{2}-a^{2}\text{$\Delta_{r}$}\sin^{2}(\theta)\right)}{L^{2}\Sigma}\text{d$\phi$}^{2}+\frac{\Sigma}{\Delta_{\theta}}\text{d$\theta$}^{2}+\frac{\Sigma}{\text{$\Delta_{r}$}}\text{dr}^{2}.$
(1)

The parameters $\Delta_{r}$, $\Delta_{\theta}$,$L$ and $\Sigma$ are defined
as,

$\displaystyle\Delta_{r}=\left(1-\dfrac{\Lambda
r^{2}}{3}\right)(r^{2}+a^{2})-2Mr,\Delta_{\theta}=\left(1+\dfrac{a^{2}\Lambda\cos^{2}\theta}{3}\right),L=\left(1+\dfrac{a^{2}\Lambda}{3}\right),\Sigma=r^{2}+a^{2}\cos^{2}\theta.$
(2)

Furthermore, the coordinates $t$ and $r$ range over all $\mathbb{R}$ while
$\theta\in[0,\pi]$ and $\phi\in[0,2\pi)$. The coefficients of metric eq. 1 are
independent of $t$ and $\phi$ thus $\partial_{t}$ and $\partial_{\phi}$ are
Killing vector fields. The roots of $\Delta_{r}=0$ and $\Sigma=0$ are
singularities of (1). Additionally, one can immediately see that when $a=0$,
eq. 1 reduces to the sds metric and when $\Lambda=0$, it reduces to the Kerr
metric with the Schwarzschild metric being the special case when both
$\Lambda$ and $a$ are set to zero.

Equations motion of photons are determined by four fundamental quantities: the
Lagrangian, the angular momentum in the $z$ axis $L_{z}$, the energy $E$, and
Carter’s constant $K$. Specifically, $L_{z}$ and $E$ are associated to the
axial symmetry and stationarity, while $K=K^{\mu\nu}p_{\nu}p_{\nu}$, the
Carter’s constant [14], is as a result of a hidden symmetry in kds spacetime.
The discussion in this work will be based on the modified form of the Carter’s
constant, $Q=K-L(aE-L_{z})^{2}$ [12],[15]. For photons, the quantities
$L_{z}$, $K$ and $E$ can be reparametrized as,

$\displaystyle\lambda=\frac{L_{z}}{E},\eta=\frac{Q}{E^{2}}.$ (3)

Null geodesics in kds spacetime follow Carter’s equations given by [6],

$\displaystyle\dfrac{\Sigma}{E}p^{r}=\pm_{r}\sqrt{R(r)},\dfrac{\Sigma}{E}p^{\theta}=\pm_{\theta}\sqrt{\Theta(\theta)},$
(4)
$\displaystyle\dfrac{\Sigma}{E}p^{\phi}=\dfrac{aL^{2}}{\Delta_{r}}(a(a-\lambda)+r^{2})-\dfrac{L^{2}}{\Delta_{\theta}\sin^{2}\theta}(a\sin^{2}\theta-\lambda),\dfrac{\Sigma}{E}p^{t}=\dfrac{L^{2}}{\Delta_{r}}((r^{2}+a^{2})^{2}-a\lambda(a^{2}+r^{2}))-\dfrac{aL^{2}}{\Delta_{\theta}}(a\sin^{2}\theta-\lambda).$
(5)

$p^{\mu}=\frac{dx^{\mu}}{d\sigma}$ is the photon’s four momentum while
$\pm_{r}$ and $\pm_{\theta}$ are the signs of $p^{r}$ and $p^{\theta}$
respectively. We have as well defined the radial and angular potentials as,

$\displaystyle
R(r)=L^{2}(r^{2}+a^{2}-a\lambda)^{2}-\Delta_{r}(\eta+L^{2}(\lambda-a)^{2}),$
(6)
$\displaystyle\Theta(\theta)=a^{2}\Delta_{\theta}L^{2}+a^{2}L^{2}\cos^{2}(\theta)-a^{2}L^{2}-2a\Delta_{\theta}\lambda
L^{2}+2a\lambda
L^{2}+\Delta_{\theta}\eta+\Delta_{\theta}\lambda^{2}L^{2}-\lambda^{2}L^{2}\cot^{2}(\theta)-\lambda^{2}L^{2}.$
(7)

The physically allowed ranges of $r$ and $\theta$ are determined by the
conditions $R(r)\geq 0$ and $\Theta(\theta)\geq 0$ respectively.
To obtain analytic solutions to eq. 4 and eq. 5, one can introduce the Mino
parameter $\tau$, which helps to decouple eq. 4 and is defined as [16],

$\displaystyle\frac{\Sigma}{E}p^{\mu}=\frac{\text{dx}^{\mu}}{d\tau}.$ (8)

The integral forms for eq. 4 and (5) are shown in eq. 120-(124).

####  II.1.1 Latitudinal Motion

In this part we examine the angular potential and its roots. The equation for
the photon’s trajectory in the $\theta$ direction will then be deduced.

Substituting the variable $u=\cos^{2}\theta$ in the angular potential gives,

$\displaystyle\Theta(u)=\frac{1}{1-u}\left(\eta+(\frac{1}{3}a^{2}\Lambda\left(L^{2}(a-\lambda)^{2}+\eta\right)+L^{2}(a-\lambda)(a+\lambda)-\eta)u+(-\frac{1}{3}a^{2}\left(L^{2}\left(\Lambda(a-\lambda)^{2}+3\right)+\eta\Lambda\right))u^{2}\right).$
(9)

Equation 9 is a second-order polynomial and, as a consequence, it possesses
two roots, $u_{\pm}$,

$\displaystyle
u_{\pm}=\chi\pm\sqrt{\frac{3\eta}{a^{2}\left(L^{2}\left(\Lambda(a-\lambda)^{2}+3\right)+\eta\Lambda\right)}+\chi^{2}},\quad\chi=\frac{1}{2}-\frac{3\left(\eta+\lambda^{2}L^{2}\right)}{2a^{2}\left(L^{2}\left(\Lambda(a-\lambda)^{2}+3\right)+\eta\Lambda\right)}.$
(10)

The inversion of $u=\cos^{2}\theta$ leads to the expression
$\theta=\arccos(\pm\sqrt{u})$. Consequently, the roots of 7 can be represented
in terms of $u_{\pm}$,

$\displaystyle\theta_{1}=\arccos(\sqrt{u_{+}}),$ (11)
$\displaystyle\theta_{2}=\arccos(\sqrt{u_{-}}),$ (12)
$\displaystyle\theta_{3}=\arccos(-\sqrt{u_{-}}),$ (13)
$\displaystyle\theta_{4}=\arccos(-\sqrt{u_{+}}).$ (14)

The ranges of $u$ and the nature of the roots in equations (11)-(14) are
determined by the sign of $\eta$. This can be investigated by the relation,

$\displaystyle
u_{+}u_{-}=-\frac{3\eta}{a^{2}\left(L^{2}\left(\Lambda(a-\lambda)^{2}+3\right)+\eta\Lambda\right)}.$
(15)

Examining the expression in (15) for different $\eta$, namely $\eta>0$,
$\eta<0$, and $\eta=0$, results in the following insights,

  1. 1.

For $\eta>0$, it can be inferred that $u_{-}u_{+}<0$. Given that $u_{+}>0$,
then $u_{-}<0$. Therefore, the roots in eq. 12 and (13) will be complex
conjugates, while those in eq. 11 and (14) will be real. Specifically,
$\theta_{1}$ will lie in the first quadrant and $\theta_{4}$ will be in the
second quadrant, indicating that $\theta_{1}<\pi/2<\theta_{4}$. Whereby, the
term quadrant refers to the region of the two-dimensional Cartesian coordinate
system that is delimited by the intersection of its horizontal and vertical
axes hence demarcating four distinct sections (the first quadrant, the second
quadrant and the third quadrant). The photon will oscillate between
$\theta_{1}$ and $\theta_{4}$, crossing the equatorial plane each time. This
type of orbit is referred to as an ordinary geodesic.

  2. 2.

For the case where $\eta<0$, we can observe that $u_{-}u_{+}>0$. In order for
the angular potential to be non-negative, the condition

$\displaystyle
a^{2}\left(L^{2}\left(\Lambda(a-\lambda)^{2}+3\right)+\eta\Lambda\right)-3\left(\eta+\lambda^{2}L^{2}\right)>0,$
(16)

has to be satisfied. As both $u_{+}$ and $u_{-}$ are positive, it follows that
all four roots given by eq. 11 to eq. 14 will be real. Thus, $\theta_{1}$ and
$\theta_{2}$ lie in the first quadrant, while $\theta_{3}$ and $\theta_{4}$
are in the second quadrant. This means that
$\theta_{1}<\theta_{2}<\pi/2<\theta_{3}<\theta_{4}$. As a result, the photon
will oscillate between $\theta_{1}$ and $\theta_{2}$ or between $\theta_{3}$
and $\theta_{4}$. The photon will be confined within either the northern or
southern hemisphere, and cannot cross the equatorial plane. Such orbits are
referred to as vortical geodesics.

  3. 3.

When $\eta=0$, we have $u_{+}u_{-}=0$. This implies two scenario : either
($u_{+}=0$ $u_{-}=2\chi$) or ($u_{-}=0$ $u_{+}=2\chi$). In our analysis, we
find that the possible scenario is $u_{+}=0$ and $u_{-}=2\chi$ such that
$2\chi<0$. Therefore, $\theta_{1}$ and $\theta_{4}$ both equal $\pi/2$, while
$\theta_{2}$ and $\theta_{3}$ are complex conjugates. This implies that the
orbit is confined to the equatorial plane.

We then derive the $\theta$ component of the trajectory by considering the
path of a photon emitted from a source point $s$ and detected by an observer
$o$. We will consider the case for which $\eta\geq 0$ and will show in the
next subsection that this is because for bound photon orbits, $\eta<0$ is not
physically relevant. The source is located at a latitude $\theta_{s}$, while
the observer is located at $\theta_{o}$.

By utilizing eq. 128c, we obtain,

$\displaystyle\tau=I_{\theta}=\dfrac{1}{\sqrt{-u_{-}\Xi}}\left[\nu_{\theta}F(x_{s}|k)-\nu_{\theta}F(x_{o}|k)\right].$
(17)

We invert eq. 17 using the relations,

$\displaystyle am(F(x,k),k)=x,$ (18)
$\displaystyle\sin(am(x\mid|k))=sn(x\mid|k),$ (19)

which results in,

$\displaystyle\frac{\cos\theta_{o}}{\sqrt{u_{+}}}=-\nu_{\theta}sn\left[\left(\tau\sqrt{-u_{-}\Xi}-\nu_{\theta}F\left(\varphi_{s}|k\right)\right)|k\right].$
(20)

Using equation (20), we can easily obtain $\theta_{o}$, which represents the
latitude of the observer,

$\displaystyle\theta_{o}(\tau)=\arccos\left(-\sqrt{u_{+}}\nu_{\theta}sn\left[\left(\tau\sqrt{-u_{-}\Xi}-\nu_{\theta}F\left(\varphi_{s}|k\right)\right)|k\right]\right).$
(21)

Accordingly, this equation encapsulates the motion of the photon in the
latitudinal direction in terms of Mino parameter from a point $\theta_{s}$ to
$\theta_{o}$. Regardless of the fact that it is assumed that a turning point
has not yet been reached, eq. 21 is in terms of the Jacobi elliptic sine
function which oscillates between $-1$ and $+1$ with period $4K(k)$. As a
result, the solution in eq. 21 will indeed oscillate between the turning
points, giving rise to full trajectories.

####  II.1.2 Radial Motion

In this subsection, we will investigate the conditions for bound orbits and
obtain the constants of motion characterizing them as well as the radii at
which they occur in the exterior of a kds black hole. Furthermore, for photons
on unbound orbits, we will derive the radial component of their trajectory by
taking into account the nature of the roots of the radial potential.

Bound photon orbits are characterized by the condition,

$\displaystyle R(r)=0,\quad\frac{dR}{dr}=0.$ (22)

Simultaneously solving eq. 22 for $\eta$ and $\lambda$ results in two cases of
equations,

  1. 1.

$\displaystyle\eta=-\frac{L^{2}r^{4}}{a^{2}},\quad\lambda=\frac{a^{2}+r^{2}}{a}.$
(23)

  2. 2.

$\displaystyle\eta=-\frac{L^{2}r^{3}\left(6a^{2}\left(\Lambda
r^{2}(3M+r)-6M\right)+a^{4}\Lambda^{2}r^{3}+9r(r-3M)^{2}\right)}{a^{2}\left(r\left(a^{2}\Lambda+2\Lambda
r^{2}-3\right)+3M\right)^{2}},$ (24)
$\displaystyle\lambda=\frac{r\left(a^{2}\left(6-\Lambda
r^{2}\right)+3r(r-3M)\right)}{a\left(r\left(a^{2}\Lambda+2\Lambda
r^{2}-3\right)+3M\right)}+a.$ (25)

$\eta$ is clearly negative for the first case. Substituting the parameters of
the first case into eq. 16 results in :
$a^{2}\left(L^{2}\left(\Lambda(a-\lambda)^{2}+3\right)+\eta\Lambda\right)-3\left(\eta+\lambda^{2}L^{2}\right)=-6L^{2}r^{2}<0$.
For bound photon orbits, this case is then not physical since it violates the
condition in eq. 16 . From this analysis, it follows that vortical motion is
not possible for kds bound photon orbits. Bound photon orbits in Kerr
spacetime are also found to lack the presence of vortical motion, as shown by
Teo [17]. It should be emphasized that the absence of vortical motion is
specific to bound orbits, unbound photon orbits may still exhibit such
behavior. However, we shall refrain from delving into the intricacies of
unbound vortical orbits in this work.
The second case is only relevant for $\eta>0$ and when $\eta$ vanishes the
photon orbit will be confined to the equatorial plane. Solving for $r$ in
$\eta=0$ gives the radial coordinate of the equatorial circular photon orbits
[18],

$\displaystyle
r_{ph+}=-\frac{2M(\mathcal{Y}-1)}{(\mathcal{Y}+1)^{2}}+2\sqrt{\frac{M^{2}((\mathcal{Y}-14)\mathcal{Y}+1)}{(\mathcal{Y}+1)^{4}}}\cos\left(\frac{\kappa}{3}+\frac{4\pi}{3}\right),$
(26) $\displaystyle
r_{ph-}=-\frac{2M(\mathcal{Y}-1)}{(\mathcal{Y}+1)^{2}}+2\sqrt{\frac{M^{2}((\mathcal{Y}-14)\mathcal{Y}+1)}{(\mathcal{Y}+1)^{4}}}\cos\left(\frac{\kappa}{3}\right).$
(27)

We have also defined

$\displaystyle\kappa=\arccos\left(\frac{M\left(2a^{2}(\mathcal{Y}+1)^{4}+M^{2}(\mathcal{Y}-1)(\mathcal{Y}(\mathcal{Y}+34)+1)\right)}{(\mathcal{Y}+1)^{6}\left(\frac{M^{2}((\mathcal{Y}-14)\mathcal{Y}+1)}{(\mathcal{Y}+1)^{4}}\right)^{3/2}}\right),\quad\mathcal{Y}=\frac{a^{2}\Lambda}{3}.$
(28)

Thus $r_{ph+}\leq r\leq r_{ph-}$ delineates the radii at which bound photon
orbits exist in the exterior of a kds black hole, hence the photon region.

$\eta$ increases from $r_{ph+}$ and attains a maximum at $r_{\eta
max}=-\frac{9M}{a^{2}\Lambda-3}$ where it begins to decrease towards
$r_{ph-}$.

$\lambda$ on the other hand is monotonic and it decreases from positive values
at $r_{ph+}$ towards negative values at $r_{ph-}$. This parameter vanishes at,

$\displaystyle
r_{zamo}=\frac{3M}{a^{2}\Lambda+3}+2\sqrt{\frac{9M^{2}}{\left(a^{2}\Lambda+3\right)^{2}}-\frac{a^{2}}{3}}\cos\left(\dfrac{1}{3}\arccos\left[-\frac{9\sqrt{3}M\left(a^{2}\left(a^{2}\Lambda+3\right)^{2}-9M^{2}\right)}{\left(a^{2}\Lambda+3\right)^{3}\left(\frac{27M^{2}}{\left(a^{2}\Lambda+3\right)^{2}}-a^{2}\right)^{3/2}}\right]\right).$
(29)

$\lambda$ is related to the angular momentum of the photons about the
$z$-axis. Thus, the case where $\lambda=0$ corresponds to photons that possess
no angular momentum in the $z$ direction, and can therefore traverse the
entire range of the angular coordinate, $\theta\in[0,\pi]$, allowing them to
cross the symmetry axis of the black hole.

A discussion of the roots of the radial potential has been given in appendix A
with the roots in eq. 107-(108). In this work we will consider the case of
double roots whose physical interpretation implies that the photon orbits are
bound at a constant radial coordinate $r_{3}=r_{4}$. In addition, we will
briefly examine solutions with $(r_{1}<r_{2}<r_{3}<r_{4}<r_{s})$, which
represents orbits that go through a turning point at $r_{4}$ and get detected
by an observer. In some instances all these roots can be less than the Cauchy
horizon hence the orbits do not have a turning point in the exterior of a
black hole. Therefore, the configuration of four distinct real roots,
$(r_{1}<r_{2}<r_{3}<r_{4}<r_{s})$, can represent orbits either escaping or
falling into the black hole. We will also consider orbits with
$r_{1}<r_{2}<r_{s}$ and $r_{3}=\bar{r}_{4}$, which represent photon
trajectories that always plunge into the black hole as they never encounter a
turning point in the exterior region. Our main objective for considering these
two scenarios is to investigate the structure of photon orbits that are nearly
bound since they play a vital role in the formation of the photon ring.
Nearly bound photon orbits refer to those orbits that, with a slight change in
initial conditions from bound photon orbits, result in trajectories that
remain close to the bound trajectories for a certain period before
exponentially diverging and escaping or plunging into the black hole. We will
solve the radial integrals for the two aforementioned cases of escape and
plunge orbits.

  * •

Escape orbits:
Making use of eq. 151 results in,

$\displaystyle\tau=I_{r}=\nu_{r}g_{E}(F(\arcsin\sqrt{\frac{(r_{3}-r_{1})(r_{o}-r_{4})}{(r_{4}-r_{1})(r_{o}-r_{3})}}|k_{E})-F(\varphi_{E,s}|k_{E})).$
(30)

Inverting for $r_{o}$ gives,

$\displaystyle
r_{o}^{E}[\tau]=\frac{r_{3}\left(r_{1}-r_{4}\right)\text{sn}^{2}\left(\left.\frac{\tau}{g_{E}}+\nu_{r}F(\varphi_{E,s}|k_{E})\right|k_{E}\right)+\left(r_{3}-r_{1}\right)r_{4}}{(r_{1}-r_{4})\text{sn}^{2}\left(\left.\frac{\tau}{g_{E}}+\nu_{r}F(\varphi_{E,s}|k_{E})\right|k_{E}\right)-r_{1}+r_{3}}.$
(31)

Equation 31 describes the radial component of a photon’s trajectory that is
emitted as a point $r_{s}$ and has a turning point at $r=r_{4}$ where it
escapes towards an observer at located at $r_{o}$. It is important to note
that this case can also represent a plunge orbit in the case that
$r_{1}<r_{2}<r_{3}<r_{4}<r_{-}<r_{+}<r_{s}$, where $r_{-}$ and $r_{+}$ are the
Cauchy and event horizons respectively.

  * •

Plunge orbits:
Making use of eq. 155 results in,

$\displaystyle\tau=I_{r}=\nu_{r}g_{P}\left(F\left(\arccos\left[\frac{(A-B)r_{o}+r_{2}B-r_{1}A}{(A+B)r_{o}-r_{2}B-r_{1}A}\right]|k_{P}\right)-F(\varphi_{P,s}|k_{P})\right).$
(32)

Inverting for $r_{o}$ gives,

$\displaystyle
r_{o}^{P}[\tau]=\frac{(r_{2}B+Ar_{1})\text{cn}\left(\left.\nu_{r}F(\varphi_{P,s}|k_{P})+\frac{\tau}{g_{P}}\right|k_{P}\right)+r_{2}B-Ar_{1}}{(A+B)\text{cn}\left(\left.\nu_{r}F(\varphi_{P,s}|k_{P})+\frac{\tau}{g_{P}}\right|k_{P}\right)-A+B}.$
(33)

Equation 33 represents the radial component of a photon’s trajectory which is
emitted at the radial coordinate $r_{s}$ and eventually plunges into the black
hole as it has no turning point in the exterior region of the black hole

####  II.1.3 Azimuthal Motion

The equation of motion in the $\phi$ direction, eq. 148, incorporates both
radial and angular integrals. This subsection discusses the solution to these
integrals.

Making use of the general solutions presented in eq. 130 and eq. 131c, we
obtain the angular integrals in the azimuthal direction as follows:

$\displaystyle
I_{\phi}=\dfrac{1}{\sqrt{-u_{-}\Xi}}\left[\nu_{\theta}\Pi(-u_{+}\mathcal{Y},\varphi_{s},k)-\nu_{\theta}\Pi(-u_{+}\mathcal{Y},\varphi_{\tau},k)\right],$
(34)
$\displaystyle\bar{I}_{\phi}=\dfrac{1}{\sqrt{-u_{-}\Xi}}\dfrac{1}{1+\mathcal{Y}}\left[\nu_{\theta}(\Pi(u_{+},\varphi_{s},k)+\mathcal{Y}\Pi(-u_{+}\mathcal{Y},\varphi_{s},k))-\nu_{\theta}(\Pi(u_{+},\varphi_{\tau}|k)+\mathcal{Y}\Pi(-u_{+}\mathcal{Y},\varphi_{\tau}|k))\right],$
(35)

where $\varphi_{\tau}=\varphi_{o}$ can also be expressed in terms of only the
Mino parameter and the initial conditions through ,

$\displaystyle\varphi_{\tau}=-\nu_{\theta}am\left[\left(\tau\sqrt{-u_{-}\Xi}-\nu_{\theta}F\left(\varphi_{s}|k\right)\right)|k\right].$
(36)

Equations (34) and (35) are the solutions to the angular integrals in the
azimuthal direction.

For the radial integral $I_{x}$, which we have shown in eq. 150, we proceed to
solve for both escape and plunge orbits.

  * •

Escape orbits:
By utilizing the solution provided in eq. 153, we obtain the expression for
$I_{x}$ as,

$\displaystyle
I_{x}=-\frac{2(r_{4}-r_{3})(\Pi_{o}-\Pi_{s})}{(r_{x}-r_{4})(r_{x}-r_{3})\sqrt{(r_{4}-r_{2})(r_{3}-r_{1})}}-\frac{2(F_{o}-F_{s})}{(r_{x}-r_{3}).\sqrt{(r_{4}-r_{2})(r_{3}-r_{1})}}.$
(37)

$\Pi_{i}$ and $F_{i}$ have been defined in eq. 154 where $i$ represents $s$ or
$o$. Furthermore, we re-define $\varphi_{E,o}$ as,

$\displaystyle\varphi_{E,o}=\arcsin\sqrt{\frac{(r_{3}-r_{1})(r_{o}-r_{4})}{(r_{4}-r_{1})(r_{o}-r_{3})}}=am(\frac{\tau+\nu_{r}g_{E}F(\varphi_{E,s}|k_{E})}{g_{E}}|k_{E}).$
(38)

The second equality has been obtained by inverting for
$\arcsin\sqrt{\frac{(r_{3}-r_{1})(r_{o}-r_{4})}{(r_{4}-r_{1})(r_{o}-r_{3})}}$
in eq. 30 using eq. 18.

Equations (34), (35) and (37) will then be substituted into eq. 148 to obtain
the azimuthal motion of photons on unbound geodesics that escape or plunge
into the black hole in the case that
$r_{1}<r_{2}<r_{3}<r_{4}<r_{-}<r_{+}<r_{s}$.

  * •

Plunge orbits:
Utilizing eq. 156, we obtain $I_{x}$ as,

$\displaystyle
I_{x,P}=\frac{(g_{P}(B-A))\left((\alpha-\alpha_{2})(R_{1,o}-\nu_{r}R_{1,s})+\alpha_{2}(F(\varphi_{P,o}\mid
k_{P})-\nu_{r}F(\varphi_{P,s}\mid
k_{P}))\right)}{r_{2}B+r_{1}A-Ar_{x}-Br_{x}}.$ (39)

The various parameters in eq. 39 have been defined in (159)-(161). We also
define $\varphi_{P,o}$ as,

$\displaystyle\varphi_{P,o}=\arccos\left[\frac{(A-B)r_{o}+r_{2}B-r_{1}A}{(A+B)r_{o}-r_{2}B-r_{1}A}\right]=am\left[\frac{\tau+\nu_{r}g_{P}F(\varphi_{P,s}\mid
k_{P})}{g_{P}}\right].$ (40)

Specifically, we obtained the second equality in (40) by inverting for
$\arccos\left[\frac{(A-B)r_{o}+r_{2}B-r_{1}A}{(A+B)r_{o}-r_{2}B-r_{1}A}\right]$
in eq. 32 using (18).

Similarly, equations (34), (35) and (39) will then be substituted into eq. 148
to obtain the $\phi$ component of a photon’s trajectory that plunges into the
black hole.

####  II.1.4 Motion in $t$ direction

Radial and angular integrals are additionally incorporated in the equation of
$t$ direction, eq. 149. The radial integral is still $I_{x}$, and we have
solutions for both escape and plunge orbits in eq. 37 and eq. 39,
respectively. The solution of the angular integral is obtained from eq. 134
as,

$\displaystyle
I_{t}=\dfrac{u_{+}}{\sqrt{-u_{-}\Xi}}\left[\nu_{\theta}J(-u_{+}\mathcal{Y},\varphi_{s}|k)-\nu_{\theta}J(-u_{+}\mathcal{Y},\varphi_{\tau}|k)\right],$
(41)

where $\varphi_{\tau}=\varphi_{o}$ is defined in eq. 36. These solutions will
then be substituted into eq. 149 to obtain the $t$ component of the photon’s
escape and plunge trajectory.

Using these solutions, we have illustrated examples of bound, escape and
plunge orbits of kds spacetime in fig. 1. Photons on the escape and plunge
geodesics in the first column are orbiting close to a bound orbit at
$r=r_{zamo}=2.56M$. Photons on the escape and plunge geodesics in the second
column on the other hand are orbiting close to a bound orbit at $r=3.00M$. The
photons have been emitted from a source located at a radial coordinate
$r_{s}=5M$ and latitude $\theta_{s}=\pi/2$. Furthermore, the observer is
situated within the domain of outer communication.

|
---|---

(a) A bound photon orbit at $r=r_{zamo,KdS}=2.56M$ with $\lambda=0$ and
$\eta=23.3576792$.

(b) A bound photon orbit at $r=3.00M$ with $\lambda=-1.8$ and $\eta=27$.

|
---|---

(c) The escape orbit has $r_{4}=2.749M$, $\lambda=0$ and $\eta=23.591256$
whereas the plunge orbit has $\lambda=0$ and $\eta=23.124102$.

(d) The escape orbit is such that $r_{4}=3.2206M$,$\lambda=-1.818$ and
$\eta=27.27$ while the plunge orbit has $\lambda=-1.782$ and $\eta=26.73$.

|
---|---

(e) The escape orbit is such that $r_{4}=2.560052765M$, $\lambda=0$,
$\eta=23.3576794$ while the plunge orbit has $\lambda=0$ and
$\eta=23.3576789$.

(f) The escape orbit is such that $r_{4}=3.0000647M$,$\lambda=-1.800000018$
and $\eta=27.00000027$. The plunge orbit is such that $\lambda=-1.799999982$
and $\eta=26.99999973$.

Figure 1:  Illustration of bound, escape (solid red) and plunge (dotted red)
orbits in kds spacetime with $a=0.9$ and $\Lambda=1.11\times 10^{-52}m^{-2}$.
The blue disk is an equatorial source.

###  II.2 Kerr de Sitter Revisited

The rkds metric is a newly proposed solution obtained through gravitational
decoupling for axially symmetric systems ([8], [9], [10]). The metric is
expressed in the Boyer-Lindquist coordinates ($t$, $r$, $\theta$, $\phi$) [7],

$\displaystyle\text{ds}^{2}=-\left(\dfrac{\hat{\Delta}-a^{2}\sin^{2}\theta}{\rho^{2}}\right)\text{dt}^{2}+\dfrac{\rho^{2}}{\hat{\Delta}}\text{dr}^{2}+\rho^{2}\text{d$\theta$}^{2}+\dfrac{((r^{2}+a^{2})^{2}-\hat{\Delta}a^{2}\sin^{2}\theta)\sin^{2}\theta}{\rho^{2}}\text{d$\phi$}^{2}-\dfrac{2a\sin^{2}\theta}{\rho^{2}}(r^{2}+a^{2}-\hat{\Delta})\text{dt}\text{d$\phi$},$
(42)

with the following definitions,

$\displaystyle\hat{\Delta}=r^{2}-2Mr+a^{2}-\dfrac{\Lambda
r^{4}}{3},\rho^{2}=r^{2}+a^{2}\cos^{2}\theta,$ (43)

where $a$, $M$, and $\Lambda$ represent the black hole spin, mass, and
cosmological constant, respectively. The rkds metric exhibits both
stationarity and axial symmetry, and the coordinates are such that
$r,t\in\mathbb{R}$ while $\theta\in[0,\pi]$ and $\phi\in[0,2\pi)$.
Singularities arise in the rkds geometry when $\hat{\Delta}=0$ and
$\rho^{2}=0$. rkds metric is independent of $t$ and $\phi$ thus $\partial_{t}$
and $\partial_{\phi}$ are Killing vector fields. When $a=0$, eq. 42 reduces to
the sds metric and when $\Lambda=0$, it reduces to the Kerr metric with the
Schwarzschild metric being the special case when both $\Lambda$ and $a$ are
set to zero.

Unlike the kds metric, the rkds metric is not a $\Lambda$ Vacuum solution and
satisfies the standard Einstein’s equation, $R_{\mu\nu}-1/2Rg_{\mu\nu}=8\pi
GT_{\mu\nu}$, where $T_{\mu\nu}$ generating this metric is given by:

$\displaystyle
T^{\mu\nu}=l^{\mu}l^{\nu}p_{r}+m^{\nu}m^{\nu}p_{\phi}+n^{\mu}n^{\nu}p_{\theta}+\epsilon
u^{\mu}u^{\nu}.$ (44)

Here, the energy density and pressures satisfy

$\displaystyle\epsilon=\frac{\Lambda r^{4}}{\rho^{4}}=-p_{r},\quad
p_{\theta}=p_{\phi}=\epsilon-\frac{2\Lambda r^{2}}{\rho^{2}}.$ (45)

Moreover, the scalar curvature of the rkds solution is expressed as [[7],
[11]],

$\displaystyle R=\frac{-4\Lambda r^{2}}{r^{2}+a^{2}\cos^{2}{\theta}}.$ (46)

Thus, besides the rkds metric not being a $\Lambda$ vacuum solution, it is
also not a constant curvature solution, but rather it exhibits a warped
curvature everywhere except on the equatorial plane ($\theta=\pi/2$). Besides,
this deformation is significant near the black hole i.e. $r\sim a$ and
disappears away from the black hole where the scalar curvature approaches a
constant value, $R\sim-4\Lambda$. This deformation is likely to be closely
related to the warped vacuum energy in the direct proximity of extreme
gravitational sources [11].

Due to the symmetries of rkds spacetime, there exists the conserved
quantities,

$\displaystyle\hat{p}_{t}=-g_{tt}\hat{p}^{t}-g_{t\phi}p^{\phi}=:\hat{E},\hat{p}_{\phi}=g_{t\phi}p^{t}+g_{\phi\phi}p^{\phi}=:\hat{L}_{z}.$
(47)

Furthermore, the Lagrangian and the Carter’s constant
$\hat{K}=\hat{K}^{\mu\nu}\hat{p}_{\nu}\hat{p}_{\nu}$ [14] are also fundamental
parameters determining photon motion in rkds spacetime. In this work, our
discussion focuses on the modified Carter constant
$\hat{Q}=\hat{K}-(a\hat{E}-\hat{L}_{z})^{2}$. For photons, these parameters
can be reparametrized as

$\displaystyle\hat{\lambda}=\frac{\hat{L}_{z}}{\hat{E}},\hat{\eta}=\frac{\hat{Q}}{\hat{E}^{2}}.$
(48)

Null geodesic equations in rkds spacetime are given by [18],

$\displaystyle\dfrac{\rho^{2}}{E}\hat{p}^{r}=\pm_{r}\sqrt{\hat{R}(r)},\dfrac{\rho^{2}}{E}\hat{p}^{\theta}=\pm_{\theta}\sqrt{\hat{\Theta}(\theta)},$
(49)
$\displaystyle\dfrac{\rho^{2}}{E}\hat{p}^{\phi}=\dfrac{(ar^{2}+a^{3}-a\hat{\Delta}-a^{2}\hat{\lambda})}{\hat{\Delta}}+\dfrac{\hat{\lambda}}{\sin^{2}\theta},\dfrac{\rho^{2}}{E}\hat{p}^{t}=\dfrac{(r^{2}+a^{2})(r^{2}+a^{2}-a\hat{\lambda})}{\hat{\Delta}}+a\hat{\lambda}-a^{2}\sin^{2}\theta,$
(50)

where

$\displaystyle\hat{R}(r)=(r^{2}+a^{2}-a^{2}\hat{\lambda})^{2}-\hat{\Delta}(\hat{\eta}+(\hat{\lambda}-a)^{2}),$
(51)
$\displaystyle\hat{\Theta}(\theta)=\hat{\eta}+a^{2}\cos^{2}\theta-\hat{\lambda}^{2}\cot^{2}\theta,$
(52) $\displaystyle\hat{\Delta}=r^{2}-2Mr+a^{2}-\dfrac{\Lambda r^{4}}{3}.$
(53)

The integral forms of eq. 49 and (50) are shown in eq. 162-(165). To perform
an analysis of bound and nearly bound photon orbits in this spacetime, we use
the same approach as in kds to generalize our calculations to rkds in a
quicker way.

####  II.2.1 Latitudinal Motion

In this subsection, we will perform an analysis of the angular potential of
rkds spacetime, analytically determining its roots. In addition, we will
derive the $\theta$ component of the photon’s trajectory.

Substituting $\hat{u}=\cos^{2}\theta$ to the angular potential gives,

$\displaystyle\hat{\Theta}(\hat{u})=\frac{1}{1-u}(\hat{\eta}+(a^{2}-\hat{\eta}-\hat{\lambda}^{2})\hat{u}-a^{2}\hat{u}^{2}),$
(54)

with its two zeros expressed as,

$\displaystyle\hat{u}_{\pm}=\frac{(a^{2}-\hat{\eta}-\hat{\lambda}^{2})\pm\sqrt{\left(-a^{2}+\hat{\eta}+\hat{\lambda}^{2}\right)^{2}+4a^{2}\hat{\eta}}}{2a^{2}}.$
(55)

From eq. 55, we arrive at the roots of eq. 52,

$\displaystyle\theta_{1}=\arccos(\sqrt{\hat{u}_{+}}),$ (56)
$\displaystyle\theta_{2}=\arccos(\sqrt{\hat{u}_{-}}),$ (57)
$\displaystyle\theta_{3}=\arccos(-\sqrt{\hat{u}_{-}}),$ (58)
$\displaystyle\theta_{4}=\arccos(-\sqrt{\hat{u}_{+}}).$ (59)

The nature of the roots in eq. 56 to (59) can be assessed by investigating the
relation,

$\displaystyle\hat{u}_{+}\hat{u}_{-}=-\frac{\hat{\eta}}{a^{2}}.$ (60)

From eq. 60, it is clear that when $\eta>0$, then $\hat{u}_{+}\hat{u}_{-}<0$.
Since $\hat{u}_{+}>0$ it implies that $\hat{u}_{-}<0$. As a result, eq. 56 and
(59) will be real such that $\theta_{1}<\pi/2<\theta_{4}$. Thus, the photon
will oscillate between $\theta_{1}$ and $\theta_{4}$ crossing the equatorial
plane each time exhibiting ordinary motion. On the other hand, when $\eta<0$,
it follows that $\hat{u}_{+}\hat{u}_{-}>0$. For this case, the angular
potential will be non-negative under the condition,

$\displaystyle a^{2}-\hat{\eta}-\hat{\lambda}^{2}>0.$ (61)

Consequently, eq. 56-(59) will be real such that,
$\theta_{1}<\theta_{2}<\pi/2<\theta_{3}<\theta_{4}$. In this case, the photon
oscillates between $(\theta_{1},\theta_{2})$ or $(\theta_{3},\theta_{4})$,
hence will be confined in the northern or southern hemisphere, thus,
exhibiting vortical motion.
When $\eta=0$, $\hat{u}_{+}\hat{u}_{-}=0$, thus $\hat{u}_{+}=0$,
$\hat{u}_{-}<0$. This gives $\hat{\theta}_{1}=\hat{\theta}_{4}=\pi/2$, hence
the photon will be confined to the equatorial plane.

We restrict our consideration to cases where $\eta\geq 0$. This choice will be
further shown in the next subsection, where we demonstrate that bound photon
orbits do not display vortical motion. It is important to note that there are
instances where vortical motion can occur for unbound photon orbits. However,
for the purpose of our current discussion, we will not delve into the
specifics of this scenario.

To obtain the $\theta$ component of the trajectory, we make use of eq. 168
which results in,

$\displaystyle\tau=\dfrac{1}{\sqrt{-u_{-}a^{2}}}\left(\nu_{\theta}F(\hat{\psi}_{s}|\hat{k})-\nu_{\theta}F(\hat{\psi}_{o}|\hat{k})\right),\hat{\psi}_{i}=\arcsin\left(\frac{\cos\theta_{i}}{\sqrt{\hat{u}_{+}}}\right).$
(62)

By utilizing the inversion relations presented in equations eq. 18 and eq. 19
on (62), we solve for $\theta_{o}$,

$\displaystyle\theta_{o}[\tau]=\arccos\left(-\sqrt{\hat{u}_{+}}\nu_{\theta}sn\left[\tau\sqrt{-\hat{u}_{-}a^{2}}-\nu_{\theta}F(\hat{\psi}_{s}|\hat{k})|\hat{k}\right]\right).$
(63)

Equation 63 represents $\theta$ component of the trajectory in terms of Mino
parameter. The equation is in terms of the Jacobi elliptic sine function which
oscillates between $-1$ and $+1$ with a period of $4K(\hat{k})$. Thus,
regardless of the assumption that a turning point has not yet been reached,
the equation goes through turning points providing complete trajectories.

####  II.2.2 Radial Motion

In this part, we will begin our investigation by analysing the condition for
the existence of bound photon orbits. We will obtain the constants of motion
that characterize these orbits, as well as the precise radii at which they are
confined. Furthermore, we will derive the equation for the radial component of
photons on unbound orbits.

Bound photon orbits in rkds spacetime satisfy the condition
$\hat{R}(r)=\hat{R}^{\prime}(r)=0$. The simultaneous solution of this
condition leads to the two distinct parameter cases,

  1. 1.

$\displaystyle\hat{\eta}=-\frac{r^{4}}{a^{2}},\hat{\lambda}=\frac{a^{2}+r^{2}}{a}.$
(64)

  2. 2.

$\displaystyle\hat{\eta}=-\frac{3r^{3}\left(4a^{2}\left(\Lambda
r^{3}-3M\right)+3r(r-3M)^{2}\right)}{a^{2}\left(3M+2\Lambda
r^{3}-3r\right)^{2}},\quad\hat{\lambda}=\frac{3a^{2}M+2a^{2}\Lambda
r^{3}+3a^{2}r-9Mr^{2}+3r^{3}}{a\left(3M+2\Lambda r^{3}-3r\right)}.$ (65)

In the first case, it is evident that $\hat{\eta}$ takes on a negative value.
Upon substituting the parameters of this case into eq. 61, the resulting
expression, $-2r^{2}$, violates the requirement that eq. 61 must be strictly
positive. Consequently, for bound photon orbits, vortical motion is not
physically admissible. However, it should be noted that vortical motion is
possible for unbound orbits in the rkds spacetime.
The second case is relevant for $\hat{\eta}>0$ and in the limit
$\hat{\eta}=0$, the photon is confined to the equatorial plane. Furthermore,
it attains a maximum at $r=3M$. The zeros of $\hat{\eta}$ will then give the
radii of equatorial circular photon orbits in the exterior of a rkds black
hole [18],

$\displaystyle\hat{r}_{ph+}=\frac{6M}{4a^{2}\Lambda+3}+6\sqrt{\frac{M^{2}\left(1-4a^{2}\Lambda\right)}{\left(4a^{2}\Lambda+3\right)^{2}}}\cos\left(\frac{\tilde{\kappa}}{3}+\frac{4\pi}{3}\right),$
(66)
$\displaystyle\hat{r}_{ph-}=\frac{6M}{4a^{2}\Lambda+3}+6\sqrt{\frac{M^{2}\left(1-4a^{2}\Lambda\right)}{\left(4a^{2}\Lambda+3\right)^{2}}}\cos\left(\frac{\tilde{\kappa}}{3}\right).$
(67)

$\tilde{\kappa}$ is defined via the expression,

$\displaystyle\tilde{\kappa}=\arccos\left(\frac{9M^{2}-2a^{2}\left(2\Lambda\left(8a^{4}\Lambda+12a^{2}-27M^{2}\right)+9\right)}{9M\left(4a^{2}\Lambda-1\right)\left(4a^{2}\Lambda+3\right)\sqrt{\frac{M^{2}\left(1-4a^{2}\Lambda\right)}{\left(4a^{2}\Lambda+3\right)^{2}}}}\right).$
(68)

Hence, the bound photon orbits in the exterior of a rkds black hole exist in
the range $\hat{r}_{ph+}<r<\hat{r}_{ph-}$. The angular momentum
($\hat{\lambda}$) of the photons on these orbits, monotonically decreases from
positive values at $\hat{r}_{ph+}$ to negative values at $\hat{r}_{ph-}$ and
vanishes at a radial coordinate,

$r_{zamo,RKdS}=2\sqrt{\frac{-2a^{4}\Lambda-3a^{2}+9M^{2}}{\left(2a^{2}\Lambda+3\right)^{2}}}\cos\left(\frac{1}{3}\arccos\left[\frac{3M\sqrt{\frac{\left(2a^{2}\Lambda+3\right)^{2}}{-2a^{4}\Lambda-3a^{2}+9M^{2}}}\left(2a^{6}\Lambda^{2}+9a^{4}\Lambda+9a^{2}-9M^{2}\right)}{\left(2a^{2}\Lambda+3\right)\left(2a^{4}\Lambda+3a^{2}-9M^{2}\right)}\right]\right)\\\
+\frac{3M}{2a^{2}\Lambda+3}.$ (69)

In the case of unbound photon orbits, we will derive the radial components of
the trajectories by considering the nature of the roots given in eq. 110 and
(111). Similar to the situation in the kds case, we will focus on roots of the
form $(\hat{r}_{1}<\hat{r}_{2}<\hat{r}_{3}<\hat{r}_{4}<\hat{r}_{s})$ and
$\hat{r}_{1}<\hat{r}_{2}<\hat{r}_{s}$, where $\hat{r}_{3}=\bar{\hat{r}}_{4}$.
These two cases correspond to escape and plunge orbits, respectively. By
utilizing eq. 182 and (185), we obtain the $r$ component of the trajectories
as follows,

  * •

Escape orbits:
Making use of eq. 182 gives,

$\displaystyle\hat{r}_{o}^{E}[\tau]=\frac{\hat{r}_{3}\hat{r}_{1,4}\text{sn}^{2}\left(\left.\frac{\tau}{\hat{g}_{E}}+\nu_{r}F(\hat{\varphi}_{E,s}|\hat{k}_{E})\right|\hat{k}_{E}\right)+\hat{r}_{4}\hat{r}_{3,1}}{\hat{r}_{1,4}\text{sn}^{2}\left(\left.\frac{\tau}{\hat{g}_{E}}+\nu_{r}F(\hat{\varphi}_{E,s}|\hat{k}_{E})\right|\hat{k}_{E}\right)-\hat{r}_{1}+\hat{r}_{3}},$
(70)

where the various parameters in this equation are to be obtained from eq. 184
by setting $i$ to $s$ or $o$ and $\hat{r}_{a,b}=\hat{r}_{a}-\hat{r}_{b}$.

  * •

Plunge orbits:
For plunge orbits, we make use of eq. 185 for which we obtain,

$\displaystyle\hat{r}_{o}^{P}[\tau]=\frac{(\hat{r}_{2}\hat{B}+\hat{A}\hat{r}_{1})\text{cn}\left(\left.\nu_{r}F(\hat{\varphi}_{P,s}|\hat{k}_{P})+\tau\sqrt{\hat{A}\hat{B}}\right|\hat{k}_{P}\right)+\hat{r}_{2}\hat{B}-\hat{A}\hat{r}_{1}}{(\hat{A}+\hat{B})\text{cn}\left(\left.\nu_{r}F(\hat{\varphi}_{P,s}|\hat{k}_{P})+\tau\sqrt{\hat{A}\hat{B}}\right|\hat{k}_{P}\right)-\hat{A}+\hat{B}}.$
(71)

The various parameters can be obtained from the definitions in eq. 188-(191).

As for the bound motion, the radial potential has double roots and hence, the
$r$ component of the trajectory needs to be evaluated at a constant radial
coordinate.

####  II.2.3 Azimuthal motion

As evidenced by eq. 164, the equation governing motion in $\phi$ direction
encompasses both the radial and angular integrals. In this subsection we will
obtain the complete solutions for these integrals which will then be
substituted into eq. 164 to obtain the photons full parameterized trajectories
in $\phi$ direction.

From eq. 170 we obtain the solution of the angular integral in the azimuthal
direction as,

$\displaystyle\hat{I}_{\phi}=\dfrac{1}{\sqrt{-u_{-}a^{2}}}\left(\nu_{\theta}\Pi(\hat{u}_{+};\hat{\psi}_{s}|\hat{k})-\nu_{\theta}\Pi(\hat{u}_{+};\hat{\psi}_{\tau}|\hat{k})\right),\hat{\psi}_{\tau}=-\nu_{\theta}am\left[\tau\sqrt{-\hat{u}_{-}a^{2}}-\nu_{\theta}F(\hat{\psi}_{s}|\hat{k})|\hat{k}\right].$
(72)

Using eq. 183, we obtain the radial integral $\hat{I}_{x}$ for escape orbits
as,

$\displaystyle\hat{I}_{x,E}=-\frac{2}{\hat{r}_{x,3}.\sqrt{\hat{r}_{4,2}\hat{r}_{3,1}}}\left(F\left(\hat{\varphi}_{E,o}|\hat{k}_{E}\right)-F\left(\hat{\varphi}_{E,s}|\hat{k}_{E}\right)+\frac{\hat{r}_{4,3}(\hat{\Pi}_{o}-\hat{\Pi}_{s})}{\hat{r}_{x,4}}\right).$
(73)

The parameters $\Pi_{o}$, $\Pi_{s}$ and $\varphi_{E,s}$ can be directly
obtained from eq. 184 by replacing $i$ with $s$ or $o$. We define
$\varphi_{E,o}$ as,

$\displaystyle\hat{\varphi}_{E,o}=\arcsin\sqrt{\frac{\hat{r}_{3,1}\hat{r}_{o,4}}{\hat{r}_{4,1}\hat{r}_{o,3}}}=am(\frac{\tau+\nu_{r}\hat{g}_{E}F(\hat{\varphi}_{E,s}|\hat{k}_{E})}{\hat{g}_{E}}|\hat{k}_{E}).$
(74)

For plunge orbits we utilize eq. 186 which gives,

$\displaystyle\hat{I}_{x,P}=\frac{((\hat{B}-\hat{A}))\left((\hat{\alpha}-\hat{\alpha}_{2})(\hat{R}_{1,o}-\nu_{r}\hat{R}_{1,s})+\hat{\alpha}_{2}(F(\hat{\varphi}_{P,o}\mid\hat{k}_{P})-\nu_{r}F(\hat{\varphi}_{P,s}\mid\hat{k}_{P}))\right)}{\sqrt{\hat{A}\hat{B}}(\hat{r}_{2}\hat{B}+\hat{r}_{1}\hat{A}-\hat{A}\hat{r}_{x}-\hat{B}\hat{r}_{x})}.$
(75)

Similarly, the parameters $\hat{R}_{1,s},\hat{R}_{1,o},\hat{\varphi}_{P,s}$
are directly obtained from eq. 189-(191) by replacing $i$ with $s$ or $o$. We
define $\hat{\varphi}_{P,o}$ as,

$\displaystyle\hat{\varphi}_{P,o}=\arccos\left[\frac{(\hat{A}-\hat{B})r_{o}+\hat{r}_{2}\hat{B}-\hat{r}_{1}A}{(\hat{A}+\hat{B})r_{o}-\hat{r}_{2}\hat{B}-\hat{r}_{1}\hat{A}}\right]=am\left[\tau\sqrt{\hat{A}\hat{B}}+\nu_{r}F(\hat{\varphi}_{P,s}\mid\hat{k}_{P})\right].$
(76)

####  II.2.4 Motion in $t$ direction

The motion in $t$ direction also entails both the radial and angular
integrals. The radial integrals is still $\hat{I}_{x}$ for which we have
obtained the solutions for both escaping and plunge orbits in eq. 73 and (75)
respectively. The angular integral in $t$ direction is obtained through eq.
171 which results in,

$\displaystyle\hat{I}_{t}=\dfrac{\hat{u}_{+}}{\hat{k}\sqrt{-u_{-}a^{2}}}\left(\nu_{\theta}[F(\hat{\psi}_{s}|\hat{k})-E(\hat{\psi}_{s}|\hat{k})]-\nu_{\theta}[F(\hat{\psi}_{\tau}|\hat{k})-E(\hat{\psi}_{\tau}|\hat{k})]\right),$
(77)

where $\varphi_{\tau}$ has been defined in eq. 72.

Equation 77 and eq. 73 or (75) can then be substituted in eq. 180 in order to
get the $t$ component of a photon on an escape or plunge orbit.

Examples of bound, escape and plunge orbits in rkds spacetime are as shown in
fig. 2. The photons on escape and plunge geodesics in the first column are
orbiting close to a bound orbit at $r=r_{zamo,RKdS}=2.56M$. Photons on the
escape and plunge geodesics in the second column on the other hand are
orbiting close to a bound orbit at $r=3.00M$. The photons have been emitted
from a latitude $\theta_{s}=\pi/2$ and at a radial coordinate $r_{s}=5M$ with
the observer located in the domain of outer communication.
In (fig. 1(a), fig. 1(b)) and (fig. 2(a), fig. 2(b)) we give examples of
selected bound orbits around a kds and rkds black hole. Bound orbits are
characterized by specific constants of motion that determine their behavior
and nearly bound orbits can be constructed by making deviations from these
constants. For the rest of the diagrams in fig. 1 and fig. 2, we examine
escape and plunge orbits for different deviations on constants of motion of
the bound orbits. The photons on these orbits are all emitted from an
equatorial source (blue disk on the diagrams) at a radial coordinate
$r_{s}=5M$. A relatively large deviation from the constants of motion is made
for the illustrations in (fig. 1(c), fig. 1(d)) and (fig. 2(c), fig. 2(d)). As
can be seen, the emitted photons make fewer orbits close to the bound orbit
before escaping or plunging into the black hole. In (fig. 1(e), fig. 1(f)) and
(fig. 2(e), fig. 2(f)), we make a very small deviation on the constants of
motion. Here, the photons make many orbits close to the bound orbit before
escaping or plunging into the black hole. Decreasing the deviation even
further would result in the photons making even more orbits close to the bound
orbit before escaping or plunging into the black hole. The rate of exponential
deviation experienced by the photons as they orbit close to the bound orbits
is measured by the Lyapunov exponent. Additionally, it is evident that photons
on escape orbits have higher constants of motion than the photons on the
corresponding bound photon orbits, while plunge orbits have less. This can be
attributed to the fact that the photon region constitutes the boundary of
escape and plunge orbits. These nearly bound orbits result in subrings of
varying thickness, contingent upon the number of half orbits they complete and
the Lyapunov exponent.

|
---|---

(a) A bound photon orbit at $r=2.56M$ with $\hat{\lambda}=0$ and
$\hat{\eta}=23.3576792$.

(b) A bound photon orbit at $r=3.00M$ with $\hat{\lambda}=-1.8$ and
$\hat{\eta}=27$.

|
---|---

(c) The escape orbit has $r_{4}=2.749M$, $\hat{\lambda}=0$ and
$\hat{\eta}=23.591256$ whereas the plunge orbit has $\hat{\lambda}=0$ and
$\hat{\eta}=23.124102$.

(d) The escape orbit is such that $r_{4}=3.2206M$,$\hat{\lambda}=-1.818$ and
$\hat{\eta}=27.27$ while the plunge orbit has $\hat{\lambda}=-1.782$ and
$\hat{\eta}=26.73$.

|
---|---

(e) The escape orbit is such that $r_{4}=2.560052765M$, $\hat{\lambda}=0$,
$\hat{\eta}=23.3576794$ while the plunge orbit has $\hat{\lambda}=0$ and
$\hat{\eta}=23.3576789$.

(f) The escape orbit is such that
$r_{4}=3.0000647M$,$\hat{\lambda}=-1.800000018$ and $\hat{\eta}=27.00000027$.
The plunge orbit is such that $\hat{\lambda}=-1.799999982$ and
$\hat{\eta}=26.99999973$.

Figure 2:  Illustration of bound, escape (solid magenta) and plunge (dotted
magenta) orbits in rkds spacetime with $a=0.9$ and $\Lambda=1.11\times
10^{-52}m^{-2}$. The blue disk is an equatorial source.

##  III Direct images, lensing rings and photon rings.

In section II we have gained insights into the general structure of bound and
nearly bound photon orbits. We have seen that the number of orbits executed by
photons before escaping or plunging into the black hole varies depending on
the proximity of their trajectory to a bound orbit. Consequently, these
diverse photon paths give rise to distinct images. The central objective of
this section is to analyze these images in detail, namely direct images,
lensing rings, and photon rings. Direct images emerge from photons following
weakly lensed orbits, while lensing rings result from photons completing one
half orbit around the black hole. As for the photon ring, it arises when
photons undergo two or more half orbits. It is important to note that in
section II, our primary focus was to gain insights into the general structure
of bound and nearly bound orbits. Consequently, the assumption that a turning
point had not yet been reached was efficient and simple to achieve that
objective. However, in this section, our primary objective is to investigate
the images of a disk surrounding the black hole through tracing the photon
trajectories analytically from an observer to various source points on the
disk. Therefore, it becomes essential to precisely account for the multiple
times that photons intersects the disk in order to accurately obtain the
corresponding images. Consequently, in this section we will utilize our
solutions in section II and incorporate the turning points into the solutions.
Specifically, the solutions will be in terms of the number of angular turning
points $m$ that the photons encounter as they orbit around the black hole.
Moreover, we will see that this quantity is directly linked to the number of
distinct images that can be generated for each specific value of $m$.

We will consider locally static observers located in a frame of the form [19]:

$\displaystyle\hat{e}_{t}=\sqrt{\frac{g_{\phi\phi}}{g_{t\phi}^{2}-g_{tt}g_{\phi\phi}}}\left(\partial_{t}-\frac{g_{t\phi}}{g_{\phi\phi}}\partial_{\phi}\right),\hat{e}_{r}=\frac{1}{\sqrt{g_{rr}}}\partial_{r},\hat{e}_{\theta}=\frac{1}{\sqrt{g_{\theta\theta}}}\partial_{\theta},\hat{e}_{\phi}=\frac{1}{\sqrt{g_{\phi\phi}}}\partial_{\phi},$
(78)

where $\hat{e}_{\mu}$ are basis vectors, with $\hat{e}_{t}$ being the timelike
vector within the domain of outer communication. This vector can be assigned
as the four-velocity of the observer. The remaining three space-like vectors
correspond to the observer’s motion in the three spatial dimensions.
Additionally, the observer has zero angular momentum, since
$\hat{e}_{t}\partial_{\phi}=0$, hence the frame defined by eq. 78 is also
referred to as the zero angular momentum frame.

The four-momentum $p^{\mu}$ can be projected onto the basis defined by eq. 78.
This projection gives the quantities that can be directly measured by the
observer in their local frame [19],

$\displaystyle p^{t}=-p_{\mu}\hat{e}_{t}^{\mu},\quad
p^{\mu}=-p_{\mu}\hat{e}_{i}^{\mu},\quad i=r,\theta,\phi.$ (79)

Figure 3:  The projection of the photon’s four-momentum onto the observer’s
frame and the solid angles $(\mathfrak{a},\mathfrak{b})$

The observation angles $(\mathfrak{a},\mathfrak{b})$, see fig. 3, are defined
as ([19],[20]),

$\displaystyle p^{r}=|\vec{P}|\cos(\mathfrak{a})\cos(\mathfrak{b}),\quad
p^{\theta}=|\vec{P}|\sin(\mathfrak{a}),\quad
p^{\phi}=|\vec{P}|\cos(\mathfrak{a})\sin(\mathfrak{b}).$ (80)

Where,
$|\vec{P}|^{2}=\left(p^{\phi}\right)^{2}+\left(p^{\theta}\right)^{2}+\left(p^{r}\right)^{2}$
for which $|\vec{P}|=p^{t}$ owing to the photon being massless.

The Cartesian coordinates $(x,y)$ that represent the apparent position of an
observer located at a distance $r_{o}$ from a black hole on the plane of the
sky are mathematically defined as [21],

$\displaystyle x=-r_{o}\cos(\mathfrak{a})\sin(\mathfrak{b}),\quad
y=r_{o}\sin(\mathfrak{a}).$ (81)

Making use of eq. 80, we have that,
$\cos(\mathfrak{a})\sin(\mathfrak{b})=p^{\phi}/p^{t}$ and
$\sin(\mathfrak{a})=p^{\theta}/p^{t}$ which when substituted into eq. 81 gives
$x=-r_{o}p^{\phi}/p^{t}$ and $y=r_{o}p^{\theta}/p^{t}$. For the case of kds we
make use of eq. 4 and (5) to obtain,

$\displaystyle
x=\frac{r_{o}\left(a\Delta_{\theta}\left(a^{2}-a\lambda+r^{2}\right)+\Delta_{r}\left(\lambda\csc^{2}(\theta)-a\right)\right)}{a\Delta_{r}\left(a\sin^{2}(\theta)-\lambda\right)-\left(a^{2}+r^{2}\right)\Delta_{\theta}\left(a^{2}-a\lambda+r^{2}\right)},$
(82) $\displaystyle
y=\frac{r_{o}\sqrt{\Delta_{\theta}\left(L^{2}(a-\lambda)^{2}+\eta\right)-L^{2}(\lambda\csc(\theta)-a\sin(\theta))^{2}}}{\mathbb{T}}.$
(83)

Here, we have defined,

$\displaystyle\mathbb{T}=L^{2}\left(\frac{\left(a^{2}+r^{2}\right)\left(a^{2}-a\lambda+r^{2}\right)}{\Delta_{r}}+\frac{a\left(\lambda-a\sin^{2}(\theta)\right)}{\Delta_{\theta}}\right).$
(84)

Inserting equations (24) and (25) into (82) and (83), gives the critical curve
of kds spacetime, which forms a boundary on the image plane delineating the
region of escape and plunge orbits. Meaning that photons that escape the black
hole region have points outside the critical curve on an observer’s image
plane, whereas photons that fall into the black hole will have points inside
the critical curve. It is convenient to express $\lambda$ and $\eta$ in terms
of $x$ and $y$ to obtain such points.

By inverting (82) and (83), $\eta$ and $\lambda$, are then expressed in terms
of the Cartesian coordinates ($x,y$) as,

$\displaystyle\lambda=\frac{\sin^{2}(\theta)\left(\left(a^{2}+r^{2}\right)\Delta_{\theta}\left(x\left(a^{2}+r^{2}\right)+ar_{o}\right)-a\Delta_{r}\left(ax\sin^{2}(\theta)+r_{o}\right)\right)}{a\Delta_{\theta}\sin^{2}(\theta)\left(x\left(a^{2}+r^{2}\right)+ar_{o}\right)-\Delta_{r}\left(ax\sin^{2}(\theta)+r_{o}\right)},$
(85)
$\displaystyle\eta=\frac{L^{2}r_{o}^{2}\left((\lambda\csc(\theta)-a\sin(\theta))^{2}-(a-\lambda)^{2}\Delta_{\theta}\right)+\mathbb{T}^{2}y^{2}}{\Delta_{\theta}r_{o}^{2}}.$
(86)

Performing the same calculation for rkds spacetime, we obtain $\hat{\lambda}$
and $\hat{\eta}$ in terms of the Cartesian coordinates as,

$\displaystyle\hat{\lambda}=\frac{\sin^{2}(\theta)\left(ar_{o}\left(a^{2}-\Delta+r^{2}\right)+x\left(\left(a^{2}+r^{2}\right)^{2}-a^{2}\Delta\sin^{2}(\theta)\right)\right)}{r_{o}\left(a^{2}\sin^{2}(\theta)-\Delta\right)+ax\sin^{2}(\theta)\left(a^{2}-\Delta+r^{2}\right)},$
(87)
$\displaystyle\hat{\eta}=-a^{2}\cos^{2}(\theta)+\frac{y^{2}\left(2a^{4}+a^{2}\Delta\cos(2\theta)-a^{2}\Delta-2a\lambda\left(a^{2}-\Delta+r^{2}\right)+4a^{2}r^{2}+2r^{4}\right)^{2}}{4\Delta^{2}r_{o}^{2}}+\lambda^{2}\cot^{2}(\theta).$
(88)

Given a locally static observer located at
($r_{o},\theta_{o},\phi_{o},t_{o}$), our focus lies in finding photon orbits
connecting $(x,y)$ on the observers screen to distinct source points
($r_{s},\theta_{s},\phi_{s},t_{s}$) on a disk around the black hole. This task
involves solving the null geodesic equations, starting from an observer $o$ to
the source $s$. Since our goal is to do this analytically, we will directly
use the solutions that we obtained in section II and modify the integrals by
simply swapping the roles of the source and observer.

Furthermore, we will express the solutions in terms of the number of turning
points that the photons encounter in the $\theta$ direction. We deem it useful
to include them so as to precisely track the emission paths of photons and
enable the generation of accurate source images. Besides, as we proceed with
our analysis, we will see the close connection between these turning points
and the multiplicity of images that can be produced for a given source.

In the case that the radial potential has four real roots, the radial
component of a trajectory traced from an observer into the kds black hole
geometry is given by,

$\displaystyle
r_{s}^{E}=\frac{r_{3}\left(r_{1}-r_{4}\right)\text{sn}^{2}\left(\left.\frac{\mathbb{I}_{\theta}}{g_{E}}+\nu_{r}F(\varphi_{E,o}|k_{E})\right|k_{E}\right)+\left(r_{3}-r_{1}\right)r_{4}}{(r_{1}-r_{4})\text{sn}^{2}\left(\left.\frac{\mathbb{I}_{\theta}}{g_{E}}+\nu_{r}F(\varphi_{E,o}|k_{E})\right|k_{E}\right)-r_{1}+r_{3}},$
(89)

while for the case of two real and two complex roots, the radial component of
the trajectory will be given by ,

$\displaystyle
r_{s}^{P}=\frac{(r_{2}B+Ar_{1})\text{cn}\left(\left.\nu_{r}F(\varphi_{P,o}|k_{P})+\frac{\mathbb{I}_{\theta}}{g_{P}}\right|k_{P}\right)+r_{2}B-Ar_{1}}{(A+B)\text{cn}\left(\left.\nu_{r}F(\varphi_{P,o}|k_{P})+\frac{\mathbb{I}_{\theta}}{g_{P}}\right|k_{P}\right)-A+B}.$
(90)

The parameter $\mathbb{I}_{\theta}$ has been defined in eq. 195 and is in
terms of the number of turning points $m$ in the $\theta$ direction. Further,
all other parameters that depend on $\eta$ and $\lambda$ will be calculated
through eq. 85 and (86). The $\phi$ and $t$ components of the trajectory will
be obtained by swapping $s$ and $o$ in eq. 148 and eq. 149. Furthermore, the
angular integrals as a function of the number of angular turning points are
expressed in eq. 195-(198).

Likewise, in rkds, the radial component for the two cases will be given by,

$\displaystyle\hat{r}_{s}^{E}=\frac{\hat{r}_{3}\hat{r}_{1,4}\text{sn}^{2}\left(\left.\frac{\hat{\mathbb{I}}_{\theta}}{\hat{g}_{E}}+\nu_{r}F(\hat{\varphi}_{E,o}|\hat{k}_{E})\right|\hat{k}_{E}\right)+\hat{r}_{4}\hat{r}_{3,1}}{\hat{r}_{1,4}\text{sn}^{2}\left(\left.\frac{\hat{\mathbb{I}}_{\theta}}{\hat{g}_{E}}+\nu_{r}F(\hat{\varphi}_{E,o}|\hat{k}_{E})\right|\hat{k}_{E}\right)-\hat{r}_{1}+\hat{r}_{3}},$
(91)
$\displaystyle\hat{r}_{s}^{P}=\frac{(\hat{r}_{2}\hat{B}+\hat{A}\hat{r}_{1})\text{cn}\left(\left.\nu_{r}F(\hat{\varphi}_{P,o}|\hat{k}_{P})+\hat{\mathbb{I}}_{\theta}\sqrt{\hat{A}\hat{B}}\right|\hat{k}_{P}\right)+\hat{r}_{2}\hat{B}-\hat{A}\hat{r}_{1}}{(\hat{A}+\hat{B})\text{cn}\left(\left.\nu_{r}F(\hat{\varphi}_{P,o}|\hat{k}_{P})+\hat{\mathbb{I}}_{\theta}\sqrt{\hat{A}\hat{B}}\right|\hat{k}_{P}\right)-\hat{A}+\hat{B}}.$
(92)

$\hat{\mathbb{I}}_{\theta}$ has been defined in eq. 199 and all the parameters
depending on $\hat{\lambda}$ and $\hat{\eta}$ are to be obtained through eq.
87 and (88). In the same way, the $t$ and $\phi$ component of the trajectory
is obtained by swapping $s$ and $o$ in eq. 179 and (180). The respective
angular integrals in terms of the photon’s number of angular turning points
are expressed in eq. 199-(201).

The $r$, $\phi$, and $t$ components will play a crucial role in mapping
observer directions to spacetime points where light rays intersect a disk
around a black hole, hence, providing a means of generating and analyzing
images. These functions have been referred to as transfer functions [1], [22].
We will focus on equatorial sources and that the photons have a radial
coordinate turning point, which implies that $\theta_{s}=\pi/2$ and
$\nu_{r}=-1$. Furthermore, $\mathbb{I}_{\theta}/\hat{\mathbb{I}}_{\theta}$,
has a discontinuity along the $x$ axis thus rendering the image curves
discontinuous. To remedy this, we reparametrize $m=n+H(y)$, where $H$
represents the Heaviside function, see also [[1], [22]]. This transformation
eliminates the aforementioned discontinuity.

For even values of $n$, the photons cross the disk from the front and for odd
values of $n$, on the other hand, the photons cross the disk from the back.
The front of the disk is the side of the disk facing the observer while the
back is the side opposite to the observer. An illustration of this is as we
show in fig. 4. Therefore, as $n$ varies, the photons alternate between
showing the front and back side of the disk.

The $n=0$ image will be a direct image of the equatorial disk’s front side. In
fig. 5, we show direct images ($n=0$) of an equatorial disk around a kds and
rkds black hole. Each of the contours on these plots is a ring with a constant
Boyer-Lindquist radial coordinate $r_{s}$, and this produces a single image
for each emission from the equatorial disk. In addition, $r_{s}$ increases
monotonically outwards, with the central region in the plots representing the
region beneath the event horizon. Direct images are strongly sensitive to the
astrophysical source profile surrounding the black hole since they are
generated by photons on weakly lensed orbits.

When $n=1$, photons cross the equatorial plane twice, resulting in two images
for each source on the equatorial disk. These images will show the disk’s
front side and backside. The second image will be more lensed in comparison to
the first because it appears after the photons have completed an additional
half orbit around the black hole. Figure 6 is an illustration of $n=1$ images.
It is apparent that in this instance, two regions emerge, implying that there
are two images for each source point on the equatorial disk. The first image
has been wrapped on itself, and the second lensed image has been superposed on
this first image. The second image is oftentimes referred to as the lensing
ring [23], [24]. When the images are merged with the critical curve, the
lensed ring gets closer to the critical curve slightly. The region beneath the
critical curve corresponds to the captured photons, which are then traced back
to the event horizon. Hence, the contours inside the critical curve are an
infinite unfolding of the event horizon. Contours outside the critical curve
correspond to sources outside the photon region, whereas the critical curve
corresponds to points within the photon region.

For $n=2$, the photons cross the equatorial plane thrice hence there will be a
nested sequence of three images for each source. The third image will be
extremely lensed as compared to the first and second. Furthermore, the third
image will be superposed on the first and second image and it will move even
more closer to the critical curve. This image and other higher order images
constitute the photon ring. In fig. 7, we illustrate images for $n=2$ in a kds
and rkds black hole. It is now clear from fig. 7(a) and (7(b)) that there is
an extra black curve superposed within these images as compared to the $n=1$
images. The extremely lensed region (black curves) is the third image in both
plots. Merging these images with the critical curve shows that the third image
is very close to the critical curve and tends to take the shape of the
critical curve as compared to the first and second images.

Finally, we analyze the images for $n=3$, which will be a sequence of four
images, two from the front and two from the backside of the disk. Because the
extremely lensed images are very close to each other and hard to distinguish
in a single plot, we have chosen to separate the front and backside images as
illustrated in fig. 8 and fig. 9. Furthermore, we have extracted out the the
region of photon capture to make our illustration clear. The two images on the
front are related to the photons’ first and third equatorial crossings, whilst
the two images on the back correspond to the photons’ second and fourth
equatorial crossings. Thus, there will be four images for a particular source.
Notably, these higher-order images exhibit extreme lensing, manifesting as
thin black curves in the plots. We already saw in the $n=2$ images that the
third image was extremely close to the critical curve. For $n=3$ we also merge
the forth image with the critical curves. It is clear that the fourth image
merges with the critical curve, resulting in an indistinguishable appearance
in both the kds and rkds spacetimes (see Figures 8(c) and 9(c)). It is evident
that the images composing the photon ring exhibit a striking resemblance to
the critical curve.

Concluding our analysis with $n=3$, we state that the pattern remains
consistent as $n$ increases: each increment of $n$ generates $n+1$ images, all
of which exhibit extreme lensing effects. Furthermore, the higher-order images
continue to appear in close proximity to the critical curve of their
respective spacetimes.

Moreover, as $n$ increases, the photons undergo a change in azimuthal angle.
Consequently, subsequent images exhibit a rotation. In fig. 10, we illustrate
the azimuthal angle shift for images corresponding to $n=1$. Additionally, we
provide a color palette indicating the measured radians for each segment.
Remarkably, in both spacetimes, it is evident that the variation in azimuthal
angle across the outer region (representing the first image) differs from that
observed in the second image. Hence, the second image manifests a distinct
rotation relative to the first image. To enhance clarity, we have extracted
out the region of photon capture in our illustration.

In summary, this section has provided a comprehensive analysis of images for
with equatorial disks around kds and rkds black holes. Subsequent images are
extremely lensed and approach the critical curve exponentially
$(\propto\exp(n\gamma))$. $\gamma$ is the Lyapunov exponent and will be
analyzed in details in the next section. Furthermore, through the change in
azimuthal angle, subsequent images undergo a rotation. The amount of this
rotation due to change in azimuthal angle is calculated by
$(\Delta\phi_{m+1}-\Delta\phi_{m})$. Additionally, the amount of delay in
detection of subsequent images is also obtained through $(\Delta
t_{m+1}-\Delta t_{m})$. We have obtained explicit forms of
($\Delta\phi_{m+1}-\Delta\phi_{m}$) and ($\Delta t_{m+1}-\Delta t_{m}$) in
appendix E. In the next section, we shall see that they are related to the
critical parameters controlling the structure of the photon ring.

|
---|---

(a) $n=0$: These photons have made zero half orbits and are weakly lensed.
They terminate on the front side of the disk resulting in one image of the
front side of the disk.

(b) $n=1$: The photons on this diagram have made one half orbit around the
black hole. Their first equatorial crossing is from the front of the disk,
while the second is from the back, resulting in two images of the equatorial
disk being detected by an observer (the front and back image).

---

(c) $n=2$: The photons on these orbits have made two half orbits around the
black hole. Their first equatorial crossing is from the front of the disk, the
second is from the back of the disk and finally the photons terminate on the
front side of the disk. Thus an observer will detect two images of the front
of the disk and one image of the back side.

Figure 4:  Illustration of the equatorial crossings of the photons traced from
an observer back into the black hole geometry for even and odd values of $n$.
It is evident that for every $n$, the photon makes contact with the equatorial
disk precisely $n+1$ times. As such, each value of $n$ corresponds to $n+1$
distinct images of the equatorial disk.

|
---|---

(a) Direct ($n=0$) image of the equatorial disk around a kds black hole.

(b) Direct ($n=0$) image of the equatorial disk around a rkds black hole.

Figure 5:  The observer is at $r_{o}=50$, $\theta_{o}=80^{\circ}$ with black
hole spin $a=0.5$ and $\Lambda=1.11\times 10^{-52}m^{-2}$. Each contour
represents a constant radial coordinate $r_{s}$ increasing from the event
horizon outwards. The photons generating these images have made zero half
orbits around the black hole and are weakly lensed.

|
---|---

(a) Images of the equatorial disk around a kds black hole for $n=1$.

(b) Images of the equatorial disk around a rkds black hole for $n=1$.

|
---|---

(c) Images of the equatorial disk around a kds black hole merged with the
critical curve (the red solid curve).

(d) Images of the equatorial disk around a rkds black hole merged with the
critical curve (the cyan solid curve).

Figure 6:  The observer is at $r_{o}=50$, $\theta_{o}=80^{\circ}$ with black
hole spin $a=0.5$ and $\Lambda=1.11\times 10^{-52}m^{-2}$.

|
---|---

(a) $n=2$ images of the equatorial disk around a kds black hole.

(b) $n=2$ images of the equatorial disk around a rkds black hole.

|
---|---

(c) $n=2$ images of the equatorial disk around a kds black hole merged with
the critical curve.

(d) $n=2$ images of the equatorial disk around a rkds black hole merged with
the critical curve.

Figure 7:  The observer is at $r_{o}=50$, $\theta_{o}=80^{\circ}$ with black
hole spin $a=0.5$ and $\Lambda=1.11\times 10^{-52}m^{-2}$.

|
---|---

(a) Two images (the colored region and the thin black curve) of the front side
of the equatorial disk

(b) Two images (the colored region and the thin black curve in the inner
boundary) for the backside of the disk

---

(c) Merging the forth image with the critical curve.

Figure 8: Illustration of $n=3$ images for a kds equatorial disk. The observer
is at $r_{o}=50$, $\theta_{o}=80^{\circ}$ with black hole spin $a=0.5$ and
$\Lambda=1.11\times 10^{-52}m^{-2}$.

|
---|---

(a) Two images (the colored region and the thin black curve) of the front side
of the equatorial disk

(b) Two images (the colored region and the thin black curve in the inner
boundary) for the backside of the disk

---

(c) Merging the forth image with the critical curve.

Figure 9: Illustration of $n=3$ images for a rkds equatorial disk. The
observer is at $r_{o}=50$, $\theta_{o}=80^{\circ}$ with black hole spin
$a=0.5$ and $\Lambda=1.11\times 10^{-52}m^{-2}$.

|
---|---

(a) Change in azimuthal angle in images of a kds black hole equatorial disk.

(b) Change in azimuthal angle in images of a rkds black hole equatorial disk.

Figure 10:  The observer is at $r_{o}=50$, $\theta_{o}=80^{\circ}$ with black
hole spin $a=0.5$ and $\Lambda=1.11\times 10^{-52}m^{-2}$.

##  IV Analysis of Critical Parameters

In section III, we have analyzed direct images, lensed rings and photon rings
of the equatorial disks around a kds and rkds blackhole. By varying the number
of half orbits undertaken by the photons around the black hole, we have
observed a nested sequence of images, each undergoing significant lensing
effects. Notably, the images comprising the photon ring exhibit a remarkable
proximity and resemblance to the critical curve. Additionally, we have
observed the change in the azimuthal angle between the first and second
images. Furthermore, as the photons make additional half orbits around the
black hole, subsequent images exhibit a delay in detection. Consequently, the
photon ring encompasses a series of subsequently demagnified images that
appear rotated and are detected at different time intervals. Our focus in this
section is to examine the critical parameters governing this demagnification,
rotation, and delay in detection of successive images. These critical
parameters encompass the Lyapunov exponent, azimuthal angle parameter, and
time delay and are evaluated on the radial coordinate of bound photon orbits.

###  IV.1 Time Delay

The critical parameter time delay is defined as the amount of time photons
take to complete half orbits around the black hole. This parameter controls
the arrival time of subsequent images. Form eq. 136 and eq. 173 we obtain the
relations for time delay in kds and rkds as,

$\displaystyle\zeta=\left(\dfrac{L^{2}}{\Delta_{r}}((r^{2}+a^{2})^{2}-a\lambda(a^{2}+r^{2}))\Gamma_{\theta}+a^{2}L^{2}\Gamma_{t}-aL^{2}(a-\lambda)\Gamma_{\phi}\right),$
(93)
$\displaystyle\hat{\zeta}=\left(\frac{(r^{2}+a^{2})(r^{2}+a^{2}-a\hat{\lambda})}{\hat{\Delta}}+a\hat{\lambda}-a^{2}\right)\hat{\Gamma}_{\theta}+a^{2}\hat{\Gamma}_{t},$
(94)

Where $\Gamma_{\phi},\Gamma_{t},\hat{\Gamma}_{t}$ have been derived in eq.
137c-(140) while $\hat{\Gamma}_{\theta}$ and $\hat{\Gamma}_{t}$ have been
derived in eq. 174 and (175), respectively.

Upon comparing these equations with eq. 203 and eq. 205, a clear
correspondence emerges: $\Delta t_{m+1}-\Delta t_{m}\approx\zeta$ and
$\hat{\Delta}\phi_{m+1}-\hat{\Delta}\phi_{m}\approx\hat{\zeta}$. Thus the
delay in the detection of successive images in the photon ring is
approximately equal to the time delay evaluated on radial points within the
photon region. We then proceed to analyze eq. 93 and (94).

Each point in the photon region is delayed differently such that time delay
decreases monotonically from
$r_{ph+}/\hat{r}_{ph+}$ to $r_{zamo,KdS}/r_{zamo,RKdS}$ where it begins to
increase monotonically towards $r_{ph-}/\hat{r}_{ph-}$. This behaviour can be
seen in the plots of fig. 11. As a result, there are two local maxima for time
delay, at $r_{ph+}/\hat{r}_{ph+}$ and $r_{ph-}/\hat{r}_{ph-}$. The maxima at
$r_{ph+}/\hat{r}_{ph+}$, however, is greater than that at
$r_{ph-}$/$\hat{r}_{ph-}$, hence, the equatorial circular prograde photon
orbit experiences the greatest amount of time delay. Time delay obtains a
minima at $r_{zamo,KdS}/r_{zamo,RKdS}$, hence, the zero angular momentum orbit
experiences the least amount of delay. We attribute this small delay to the
fact that the photon lacks angular momentum at this radial coordinate, and
thus its orbit is solely due to the dragging of inertial frames around the
black hole. Furthermore, as $a\longrightarrow 1$ and
$\theta\longrightarrow\pi/2$, the equatorial circular prograde photon orbit
starts to experience very large time delays. We further express the time delay
in standard units $(\zeta_{dim}=t_{g}\zeta)$ so as to have dimensional
estimates, with $t_{g}=\frac{GM}{c^{3}}$.
Figure 11 illustrates the time delays of our results applied to sgrA and M87.
In sgrA, when $a=0.94,0.8,0.5$, the time delay range within
$\zeta_{dim}\in[4.98352,5.5589],[5.10068,5.34345]$ and $[5.26041,5.32105]$
minutes respectively. In M87, for the same values of black hole spin, the time
delays range within $\zeta_{dim}=[5.62376,6.27306],[5.75597,6.02994]$ and
$[5.93623,6.00465]$ days respectively. As a result, sgrA exhibits very rapid
changes as compared to M87, which can remain unchanged for many days. The eht
collaboration’s recent observation of M87 also indicates that the image
features of this black hole are equivalent across all four days [25].
Furthermore, unlike M87, the eht results of sgrA show that the image features
of sgrA revolve in a $4-30$ minutes range [5].

In the case of a rapidly rotating black hole, such as $a=0.99999$ with an
observer located in the equatorial plane, the time delay for sgrA ranges
within $\zeta_{dim}=[4.92024,463.102]$ minutes, whereas M87 will be
$\zeta_{dim}=[5.55235,522.598]$ days. As a result, in such a scenario, sgrA’s
equatorial circular prograde photon orbit experiences a time delay of
approximately $463.102$ minutes, while M87’s experiences a time delay of
approximately $522.598$ days. Thus, for very high values of black hole spin
and large angles of inclination, M87 can remain unchanged even up to a period
of a year. However, such large values of spin and inclinations have been ruled
out in recent observations[18].
For the case of vanishing $\Lambda$, which represents a Kerr black hole, we
find that when $a=0.94,0.8,0.5$, the time delay applied to sgrA gives
$\zeta_{dim,Kerr}\in[4.98213,5.55707],[5.10001,5.34276]$ and
$[5.26032,5.32098]$ minutes respectively. Moreover, considering M87, for the
same values of black hole spin we obtain a time delay of
$\zeta_{dim,Kerr}\in[5.62219,6.271],[5.75522,6.02916]$ and $[5.93612,6.00457]$
days respectively. Thus, we observe that the time delays when
$\Lambda=1.11\times 10^{-52}m^{-2}$ and $\Lambda=0$ in both sgrA and M87 are
approximately the same.
Increasing the value of $\Lambda$ begins to have a noticeable effect on
$\zeta$ when $\Lambda\approx 10^{-2}m^{-2}$ where the values of time delay
begins to increase. This is evident from fig. 12 where the curves representing
for $\Lambda=0,1.11\times 10^{-52}m^{-2},10^{-5}m^{-2}$ are indistinguishable.
However, as $\Lambda$ begins to have a noticeable effect we observe that an
increase in $\Lambda$ results in an increase on the amount of time photons
take to complete their respective half orbits in both spacetimes.
Comparing time delays in kds and rkds for large values of the cosmological
constant, we observe that they are no longer indistinguishable but rather the
amount of time photons take to complete half orbits is larger in kds black
hole. This can is also illustrated in fig. 13.
In the limit $a=0$, which represents a sds black hole, we obtain $\zeta_{SdS}$
$=3\sqrt{3}\pi\sqrt{\frac{M^{2}}{1-9\Lambda M^{2}}}$. For $\Lambda=0$, this
relation reduces to $3\sqrt{3}\pi M=16.3242$ which is the time delay for a
Schwarzschild black hole. For $\Lambda=1.11\times
10^{-52}m^{-2},10^{-5}m^{-2},10^{-2}m^{-2},10^{-1}m^{-2}$ we obtain the values
of time delay as $16.3242,16.3249,17.1124$ and $51.6216$ respectively. We then
observe that just as the rotating case, the cosmological constant begins to
have a significant effect when $\Lambda\approx 10^{-2}m^{-2}$. Furthermore,
when $\Lambda$ has a significant effect, we observe that time delay increases
with the increasing value of $\Lambda$.

|
---|---

(a) Time delay in sgrA.

(b) Time delay in M87.

Time delay for a kds black hole applied to sgrA and M87

|
---|---

(c) Time delay in sgrA.

(d) Time delay in M87

Time delay for a rkds black hole applied to sgrA and M87 Figure 11:  Time
delay in kds and rkds applied to sgrA and M87 black holes for different values
of spin. The observer is inclined at $30^{\circ}$ and $\Lambda=1.11\times
10^{-52}m^{-2}$.

(a) Time delay for a kds black hole

(b) Time delay for a rkds black hole.

Figure 12: Time delay in both kds and rkds for different values of $\Lambda$
with $a=0.5$ and the observers are inclined at $30^{\circ}$. The curves for
$\Lambda=0,1.11\times 10^{-52}m^{-2},10^{-5}m^{-2}$ are all indistinguishable.

(a) Time delay for $\Lambda=0.02m^{-2}$.

(b) Time delay for $\Lambda=0.01m^{-2}$.

Figure 13: Comparison of time delay in kds and rkds with $a=0.5$ and the
observers are inclined at $30^{\circ}$.

###  IV.2 Lyapunov Exponent

Small changes in the initial conditions can cause the distances between two
orbits to diverge or converge exponentially. Our reference orbits in this work
are bound photon orbits, and our goal is to investigate the rate of
exponential deviation experienced by photons that come very close to the
photon region. The Lyapunov exponents measure the rate of this exponential
divergence/convergence. The Lyapunov exponent thus embodies the rate of
exponential deviation of nearly bound photons from bound photon orbits after a
number of half orbits around the black hole. Besides, the Lyapunov exponent
controls the exponential demagnification of successive images. From eq. 145
and (176) we obtain the Lyapunov exponent in kds and rkds spacetime as,

$\displaystyle\gamma=\sqrt{\frac{1}{3}\Lambda\left(a^{2}+6r^{2}\right)\left(L^{2}(a-\lambda)^{2}+\eta\right)+L^{2}\left(a^{2}-\lambda^{2}+6r^{2}\right)-\eta}\Gamma_{\theta},$
(95)
$\displaystyle\hat{\gamma}=\sqrt{a^{2}+2\Lambda\text{r}^{2}\left((a-\hat{\lambda})^{2}+\hat{\eta}\right)-\hat{\eta}-\hat{\lambda}^{2}+6\text{r}^{2}}\hat{\Gamma}_{\theta}.$
(96)

Analysis of eq. 95 and (96) show that the Lyapunov exponents are always
positive in both spacetimes. This occurs because bound photon orbits around
kds and rkds black holes are unstable to radial perturbations. As a result,
the distance between these orbits and any nearly bound photon orbit will
diverge exponentially with time before escaping or plunging into the black
hole. Furthermore, from fig. 14, we observe that the general behaviour of
Lyapunov exponent greatly depends on the observer’s angle of inclination. That
is to say, for observers inclined at $90^{\circ}$, the Lyapunov exponent
decreases from $\pi$ at $r_{ph+}$ towards a minimum and begins to increase to
$\pi$ at $r_{ph-}$. Thus at this angle of inclination fig. 14(b), orbits close
to the retrograde and prograde equatorial circular orbits experience the same
rate of exponential deviation. The rate of exponential deviation from
equatorial circular photon orbits is no longer equal as the angle of
inclination decreases, i.e. at $30^{\circ}$ fig. 14(a). Rather, the equatorial
retrograde circular orbit exhibits more instability than the prograde orbit.
Furthermore, as black hole spin increases, the rate of exponential deviation
decreases. Nearly bound photon orbits that escape the black hole will then
generate a series of subrings whose thickness is controlled by these Lyapunov
exponents.

In the case of a Kerr black hole with $\Lambda=0$, the Lyapunov exponent
values are approximately equal to those obtained with an extremely small
cosmological constant of $\Lambda=1.11\times 10^{-52}m^{-2}$. However, as the
cosmological constant increases, a significant impact on the Lyapunov exponent
becomes apparent when $\Lambda\approx 10^{-2}m^{-2}$. This can be seen in fig.
15, where the curves for $\Lambda=0,1.11\times 10^{-52}m^{-2},10^{-5}m^{-2}$
cannot be distinguished in both spacetimes. However, when the effect of the
cosmological constant becomes noticeable, larger values result in the rate of
exponential deviation from equatorial circular prograde orbits to increase,
and leads to a decrease on the rate of deviation from equatorial circular
retrograde orbits in a kds black hole, as demonstrated in fig. 15(a) for
observers inclined at $30^{\circ}$. In contrast, for rkds fig. 15(b), an
increase in the cosmological constant leads to an increase in the rate of
exponential deviation from both the equatorial circular prograde and
retrograde orbits. We attribute this difference in the general behavior of the
Lyapunov exponent in rkds to the presence of warped curvature in the vicinity
of this black hole, which can cause changes in the behavior of the rate of
exponential deviation of orbits. As demonstrated in fig. 15(c) and fig. 15(d),
when observers are located in the equatorial plane, the Lyapunov exponent
displays a similar overall behavior in both spacetimes. The Lyapunov exponent
in both spacetimes exhibits a monotonic decrease from $\pi$ at the equatorial
circular prograde orbit to a local minimum, followed by a monotonic increase
to a value of $\pi$ at the equatorial circular retrograde orbit. It is now
evident that the general behavior of the rate of exponential deviation is
similar in both spacetimes. We attribute this similarity to the fact that the
rkds solution has constant curvature on the equatorial plane, resulting in a
similar structure of the general behavior of the rate of exponential deviation
as that of kds, which also has constant curvature.
Furthermore, a comparison between rkds and kds black holes, fig. 16, shows
that for large values of the cosmological constant, the rate of exponential
deviation for photons about rkds is higher than that experienced by photons
about kds black hole.
In the static case, we obtain that the Lyapunov exponent is equal to $\pi$,
indicating that the cosmological constant has no effect on the rate of
exponential deviation from sds photon spheres unlike in kds/rkds where
$\Lambda$ has an effect on the rate of exponential deviation. Negative
Lyapunov exponents may indicate stability, and thus exponential convergence of
orbits. This scenario can be possible in black hole spacetimes which have
stable bound photon orbits.

|
---|---

(a) kds Lyapunov exponent for observers inclined at $30^{\circ}$

(b) kds Lyapunov exponent for observers inclined at $90^{\circ}$

|
---|---

(c) rkds Lyapunov exponent for observers inclined at $30^{\circ}$

(d) rkds Lyapunov exponent for observers inclined at $90^{\circ}$

Figure 14:  Lyapunov exponent for different values of black hole spin with
$\Lambda=1.11\times 10^{-52}m^{-2}$.

(a) Lyapunov exponent for kds black hole for an observer inclined at
$30^{\circ}$.

(b) Lyapunov exponent for rkds black hole for an observer inclined at
$30^{\circ}$.

(c) Lyapunov exponent for kds black hole for an observer inclined at
$90^{\circ}$.

(d) Lyapunov exponent for rkds black hole for an observer inclined at
$90^{\circ}$.

Figure 15: Lyapunov exponent for different values of $\Lambda$ in kds and rkds
black hole spin with $a=0.5$. The curves for $\Lambda=0,1.11\times
10^{-52}m^{-2},10^{-5}m^{-2}$ cannot be distinguished from each other.

(a) Lyapunov exponent for both kds and rkds black hole with
$\Lambda=0.02m^{-2}$

(b) Lyapunov exponent for both kds and rkds black hole with
$\Lambda=0.01m^{-2}$

Figure 16:  Comparing the Lyapunov exponent in kds and rkds black hole with
$a=0.5$ and observers are inclined at $30^{\circ}$.

###  IV.3 Change in azimuthal angle

The change in azimuthal angle parameter measures the amount with which the
azimuthal angle changes per half orbit. This effect is reflected on the photon
ring structure as it causes a rotation on the subsequent subrings. Thus, the
change in azimuthal angle parameter governs the rotation of subrings. Using
eq. 135 and (172) we obtain the change in azimuthal angle parameter in both
kds and rkds as,

$\displaystyle\delta=\left(\frac{\left(aL^{2}\right)\left(a(a-\lambda)+\text{r}^{2}\right)}{\text{$\Delta_{r}$}}\Gamma_{\theta}-L^{2}a\Gamma_{\phi}+L^{2}\lambda\bar{\Gamma}_{\phi}\right),\hat{\delta}=\left(\frac{ar^{2}+a^{3}-a\hat{\Delta}-a^{2}\hat{\lambda}}{\hat{\Delta}}\right)\hat{\Gamma}_{\theta}+\hat{\lambda}\hat{\Gamma}_{\phi}.$
(97)

Upon comparing these equations with eq. 202c and (204), it becomes evident
that $\Delta\phi_{m+1}-\Delta\phi_{m}\approx\delta$ and
$\hat{\Delta}\phi_{m+1}-\hat{\Delta}\phi_{m}\approx\hat{\delta}$. Thus, the
change in azimuthal angle parameter between successive images of the photon
ring is approximately equivalent to evaluating the change in azimuthal angle
at the critical radial coordinate. We then proceed to analyze eq. 52.

This change in azimuthal angle parameter decreases from
$r_{ph+}$/$\hat{r}_{ph+}$ to $r_{ph-}$/$\hat{r}_{ph-}$, with a discontinuity
at $r_{zamo,KdS}$/$\hat{r}_{zamo}$. The discontinuity occurs when the orbits
are transitioning from prograde to retrograde and can be discarded by adding
$2\pi H(r-r_{zamo})/2\pi H(r-\hat{r}_{zamo})$ [1] to eq. 97 ,where $H$ is a
Heaviside function. This parameter specifies the amount of azimuthal lapse for
the orbits and, as a result, encodes the factor through which subsequent
images will appear rotated on the observer’s screen.
From fig. 17, the change in azimuthal angle parameter is monotonic. Besides
that, we have done away with the discontinuity that appears at $r_{zamo,KdS}$.
The equatorial prograde circular orbits experience the greatest azimuthal
shift, while the retrograde orbits have the least. Furthermore, as the black
hole spin increases, so does the change in azimuthal angle parameter. Thus,
the subrings for black holes with larger spin will experience greater rotation
as compared to those with small spin.

When $\Lambda=0$, we observe that the change in azimuthal angle parameter is
still indistinguishable with the result of $\Lambda=1.11\times
10^{-52}m^{-2}$. Moreover, increasing $\Lambda$ starts to have a considerable
effect on the change in azimuthal angle when $\Lambda\approx 10^{-2}m^{-2}$.
This is apparent from fig. 18 where we observe that the curves for
$\Lambda=0,1.11\times 10^{-52}m^{-2},10^{-5}m^{-2}$ overlap. However, when
$\Lambda$ is approximately equal or greater than $10^{-2}m^{-2}$, the change
in azimuthal angle takes on larger values. Comparing this parameter in a kds
and rkds for larger values of $\Lambda$, fig. 19, we observe that photons
around a kds black hole undergo a greater change is azimuthal angle than those
around a rkds black hole.

In the static case, we find that the change in azimuthal angle equals $\pi$
with the exception of the bound photon orbits going through the poles which
experience zero azimuthal shift.

|
---|---

(a) kds Change in azimuthal angle for observers inclined at $30^{\circ}$

(b) kds Change in azimuthal angle for observers inclined at $90^{\circ}$

|
---|---

(c) rkds Change in azimuthal angle for observers inclined at $30^{\circ}$

(d) rkds Change in azimuthal angle for observers inclined at $90^{\circ}$

Figure 17:  Change in azimuthal angle with $\Lambda=1.11\times
10^{-52}m^{-2}$.

(a) Change in azimuthal angle for kds black hole.

(b) Change in azimuthal angle for rkds black hole.

Figure 18: Change in azimuthal angle for different values of the cosmological
constant with $a=0.5$ and observers inclined at $30^{\circ}$. The curves for
$\Lambda=0,1.11\times 10^{-52}m^{-2},10^{-5}m^{-2}$ are all indistinguishable.

(a) Change in azimuthal angle for $\Lambda=0.02m^{2}$.

(b) Change in azimuthal angle for $\Lambda=0.1m^{2}$.

Figure 19: Comparing the change in azimuthal angle in kds and rkds black hole
with $a=0.5$ and observers inclined at $30^{\circ}$.

##  V Conclusion

This work is structured into three sections, each building upon the preceding
one. We commence by deriving analytical solutions for bound and nearly bound
photon orbits, strategically chosen as the starting point to comprehend the
overall structure of these orbits in both kds and rkds spacetimes. Given that
these orbits serve as the fundamental basis for the analysis in this work,
gaining a comprehensive understanding of their general structure facilitates
more intuitive insights into the subsequent sections. Leveraging our obtained
solutions, we have conducted an analytic ray-tracing in kds and rkds
spacetimes. We have embarked on a detailed analysis of equatorial disk images,
considering the varying number of half orbits taken by photons around the
black hole. An increase in the number of half orbits results in the emergence
of additional lensed images, with each successive image exhibiting greater
lensing effects and approaching the critical curve. Notably, the $n\geq 2$
images demonstrates a closer approximation to the critical curve compared to
the $n=1$ image, affirming their potential as robust testing ground for gr.
Moreover, an analysis of the azimuthal angle variations in the $n=1$ images
has showcased the rotation between the first and second images. Our
observations among other works in the literature, provide compelling evidence
for the profound capabilities of using the photon ring to test gr in the
strong gravity regimes as the critical curve itself is not experimentally
observable. Subsequently, we have proceeded to derive and analyze the critical
parameters that govern the intricate structure of the photon ring.

We observe that time delay initially decreases from the equatorial circular
prograde photon orbit to the zamo, and then increases towards the equatorial
circular retrograde orbit. Photons in equatorial circular prograde orbits
experience the longest delays, while those in zamo orbits experience the
shortest delays. The shorter delays in zamo orbits can be attributed to the
Lense-Thirring effect, which arises from the dragging of inertial frames. We
have applied our findings to the black holes sgrA and M87, considering various
spin values and a cosmological constant of $\Lambda=1.11\times
10^{-52}m^{-2}$. For example, in sgrA, for spin values of $a=0.94,0.8,0.5$,
the time delay is approximately 5 minutes. In contrast, for the same spin
values in M87, the time delay is on the scale of 5 days. This implies that M87
can appear unchanged for days, while sgrA may exhibit noticeable changes
within minutes. For large values of black hole spin and angle of inclination
(e.g., $a=0.99999$ and $\theta=\pi/2$), the time delay for sgrA ranges from
approximately $5-460$ minutes, while for M87, it ranges from about $5-522$
days. However, large values of black hole spin and angle of inclination have
been shown not to pass the constraints of recent eht observations. We observe
that these results are approximately equal to those of a rkds and Kerr black
hole, $\Lambda=0$. Increasing the cosmological constant initiates a
discernible influence on time delay, primarily when $\Lambda$ values approach
or exceed $10^{-2}m^{-2}$. For such large values, the time delay increases as
$\Lambda$ increases. In addition, in circumstances where $\Lambda$ is
considerable, photons orbiting a kds black hole will be delayed more than
those orbiting a rkds black hole.

Photons on the equatorial circular prograde orbit undergo the maximum
azimuthal shift, while those following the equatorial circular retrograde
orbit exhibit the least azimuthal shift per half-orbit. For
$\Lambda=1.11\times 10^{-52}m^{-2}$, the change in azimuthal angle is
indistinguishable in kds, rkds and Kerr ($\Lambda=0$). Additionally, the
change in the azimuthal angle increases with an increase in black hole spin.
Notably, when the value of the cosmological constant approaches or exceeds
$10^{-2}m^{-2}$, it starts to manifest an observable effect on the change in
the azimuthal angle. For such significant values, increasing the cosmological
constant results in an increase in the azimuthal shift experienced by photons
per half-orbit. Furthermore, for these large values of $\Lambda$, photons
traversing a kds black hole experience a more substantial azimuthal shift than
those orbiting a rkds black hole.

The Lyapunov exponent, on the other hand, does not reveal a global behavior
but rather depends on the observer’s inclination. This parameter exhibits
positive values in both spacetimes, indicating that the trajectories
approaching the bound photon orbits tend to exponentially deviate with time.
For $\Lambda=1.11\times 10^{-52}m^{-2}$, both the kds and rkds exhibit the
same instability rate for the circular prograde and retrograde orbits, when
observed from the equatorial plane. However, this behavior drastically changes
when observed away from the equatorial plane, where the rate of exponential
deviation from the equatorial circular prograde orbit is less than the rate of
deviation from the equatorial circular retrograde orbit. As we increase the
value of $\Lambda$, we begin to observe a notable impact on the Lyapunov
exponent, particularly when $\Lambda$ approaches and exceeds $10^{-2}m^{-2}$.
For such values, the behavior of the Lyapunov exponent in kds and rkds
manifests distinct general behaviors for observers inclined away from the
equatorial plane. For kds, large values of $\Lambda$ increase the rate of
exponential deviation from the equatorial circular prograde photon orbit,
while decreasing the rate of exponential deviation from the equatorial
circular retrograde photon orbit. On the other hand, for rkds, increasing
$\Lambda$ increases the rate of exponential deviation from both the equatorial
circular prograde and retrograde photon orbits. However, when the observer is
inclined at $\theta=\pi/2$, the Lyapunov exponent in both spacetime exhibits a
similar general behavior, where it decreases from a value of $\pi$ at the
equatorial circular photon orbit to a local minimum, where it again increases
to $\pi$ on the equatorial circular retrograde photon orbit. We attribute this
difference in the general behavior of the Lyapunov exponent to the warped
curvature of the rkds spacetime in the vicinity of the black hole. However, on
the equatorial plane, the warped curvature vanishes, and the solution becomes
a constant curvature solution, resulting in a similar general behavior of the
rate of exponential deviation to that of a kds black hole. Overall, this work
highlights the Lyapunov exponent as the sole critical parameter that exhibits
different general behaviors in kds and rkds black holes.

Thus, the accepted cosmological value $\Lambda\sim 10^{-52}m^{-2}$ plays no
relevant role in the physical situations studied here. Even though large
values of the cosmological constant may not be astrophysically relevant, we
note that in rkds spacetime, such large values are feasible, as $\Lambda$ has
the interpretation of vacuum energy and can be deformed in the vicinity of
this black hole [7],[11]. Furthermore, it is worth mentioning that in strong
gravitational field environments there is the possibility of an amplification
of the $\Lambda$ value due to the quantum effects [26]. This kind of mechanism
may justify a high $\Lambda$ value in the vicinity of a black hole compared
with its asymptotic value. From our work, we observe that detection of the
critical parameters analysed in this work might still not distinguish a Kerr
black hole from kds or rkds since such parameters are approximately the same
for the typical cosmological relevant value of the cosmological constant.

##  Appendix A  Calculating roots using Ferrari’s method

The radial potential, eq. 6, can be re-expressed as,

$\displaystyle R(r)=r^{4}+\mathcal{E}r^{2}+\mathcal{F}r+\mathcal{G},$ (98)

for which we have defined the coefficients as,

$\displaystyle\mathcal{E}=a^{2}-\frac{3\left(\eta+\lambda^{2}L^{2}\right)}{L^{2}\left(\Lambda(a-\lambda)^{2}+3\right)+\eta\Lambda},\mathcal{F}=\frac{6M\left(L^{2}(a-\lambda)^{2}+\eta\right)}{L^{2}\left(\Lambda(a-\lambda)^{2}+3\right)+\eta\Lambda},\mathcal{G}=-\frac{3a^{2}\eta}{L^{2}\left(\Lambda(a-\lambda)^{2}+3\right)+\eta\Lambda}.$
(99)

Equation 98 is already a depressed quartic and the general form of its roots
are given by Ferrari’s method [27] as,

$\displaystyle
r=\pm_{1}\sqrt{\frac{p_{1}}{2}}\pm_{2}\sqrt{-\left(\frac{\mathcal{E}}{2}+\frac{p_{1}}{2}\pm_{1}\sqrt{\frac{2}{p_{1}}}\frac{\mathcal{F}}{4}\right)},$
(100)

$p_{1}$ is a root of the resolvent cubic of this quartic polynomial. The
resolvent cubic is obtained as follows. Firstly, eq. 98 can be re-written as,

$\displaystyle\left(r^{2}+\frac{\mathcal{E}}{2}\right)^{2}=-\mathcal{F}r-\mathcal{G}+\frac{\mathcal{E}^{2}}{4},$
(101)

A parameter $p$ is then introduced into eq. 101 by adding
$2r^{2}p+\mathcal{E}p+p^{2}$ on the left and right hand side such that the
right side is perfect square. This results in,

$\displaystyle\left(r^{2}+p+\frac{\mathcal{E}}{2}\right)^{2}=-\mathcal{F}r-\mathcal{G}+\frac{\mathcal{E}^{2}}{4}+2r^{2}p+\mathcal{E}p+p^{2}.$
(102)

Recall that a quadratic equation is a perfect square if and only if it’s
discriminant equal zero. Thus, obtaining the discriminant in $r$ for the right
hand side of eq. 102 and equating to zero results in,

$\displaystyle
8p^{3}+8\mathcal{E}p^{2}+(2\mathcal{E}^{2}-8\mathcal{G})p-\mathcal{F}^{2}=0.$
(103)

Equation 103 is the resolvent cubic of eq. 98 whose roots are obtained by
first expressing it as a depressed cubic through the relation
$p=v-\mathcal{E}/3$. This results in,

$\displaystyle
v^{3}+\mathcal{P}t+\mathcal{U}=0,\mathcal{U}=\mathcal{G}-\frac{\mathcal{E}^{3}}{108}-\frac{\mathcal{F}^{2}}{8},\mathcal{T}=-\mathcal{G}-\frac{\mathcal{E}^{2}}{12}$
(104)

Cardano’s method [27] gives the roots of eq. 104, whereby we judiciously
select the real root,

$\displaystyle
v=\mathcal{Z}_{+}+\mathcal{Z}_{-},\mathcal{Z}_{\pm}=\sqrt[3]{\pm\sqrt{\frac{\mathcal{U}^{2}}{4}+\frac{\mathcal{T}^{3}}{27}}-\frac{\mathcal{U}}{2}}.$
(105)

Recall that we used the substitution $p=v-\mathcal{E}/3$. Plugging $v$ into
this relation gives,

$\displaystyle p_{1}=\mathcal{Z}_{+}+\mathcal{Z}_{-}-\mathcal{E}/3,$ (106)

Substituting eq. 106 into (100) and distributing the signs $\pm_{1}$/$\pm_{2}$
results in the four roots of the radial potential,

$\displaystyle
r_{1}=-\sqrt{\frac{p_{1}}{2}}-\sqrt{-\frac{\mathcal{E}}{2}+\frac{\mathcal{F}}{2\sqrt{2p_{1}}}-\frac{p_{1}}{2}},r_{2}=-\sqrt{\frac{p_{1}}{2}}+\sqrt{-\frac{\mathcal{E}}{2}+\frac{\mathcal{F}}{2\sqrt{2p_{1}}}-\frac{p_{1}}{2}},$
(107) $\displaystyle
r_{3}=-\sqrt{\frac{p_{1}}{2}}-\sqrt{-\frac{\mathcal{E}}{2}-\frac{\mathcal{F}}{2\sqrt{2p_{1}}}-\frac{p_{1}}{2}},r_{4}=\sqrt{\frac{p_{1}}{2}}+\sqrt{-\frac{\mathcal{E}}{2}+\frac{\mathcal{F}}{2\sqrt{2p_{1}}}-\frac{p_{1}}{2}}.$
(108)

The roots in eq. 108 can be analyzed using the discriminant of eq. 98,

$\displaystyle\Delta_{0}=16\mathcal{E}^{4}\mathcal{G}-4\mathcal{E}^{3}\mathcal{F}^{2}-128\mathcal{E}^{2}\mathcal{G}^{2}+144\mathcal{E}\mathcal{F}^{2}\mathcal{G}-27\mathcal{F}^{4}+256\mathcal{G}^{3}.$
(109)

$\Delta_{0}<0$ results in two real $(r_{1}<r_{2}<r_{-}<r_{+})$ and two complex
conjugate roots $(r_{3}=\bar{r}_{4})$. Further, $\Delta_{0}>0$ gives four real
roots $(r_{1}<r_{2}<r_{3}<r_{4})$ if $\mathcal{E}<0$ and only complex roots
$(r_{1}=\bar{r}_{2},r_{3}=\bar{r}_{4})$ if $\mathcal{E}>0$. $\Delta_{0}=0$
produces double roots $(r_{1}<r_{2}<r_{3}=r_{4})$. Generally, irrespective of
the nature of the roots, the following relations are satisfied:
$r_{1}+r_{2}+r_{3}+r_{4}=0$ and $r_{1}r_{2}r_{3}r_{4}=\mathcal{G}$.
The same approach applies to the radial potential of the rkds metric which
results in,

$\displaystyle\hat{r}_{1}=-\sqrt{\frac{\hat{p}_{1}}{2}}-\sqrt{-\frac{\hat{\mathcal{E}}}{2}+\frac{\hat{\mathcal{F}}}{2\sqrt{2\hat{p}_{1}}}-\frac{\hat{p}_{1}}{2}},\hat{r}_{2}=-\sqrt{\frac{\hat{p}_{1}}{2}}+\sqrt{-\frac{\hat{\mathcal{E}}}{2}+\frac{\hat{\mathcal{F}}}{2\sqrt{2\hat{p}_{1}}}-\frac{\hat{p}_{1}}{2}},$
(110)
$\displaystyle\hat{r}_{3}=-\sqrt{\frac{\hat{p}_{1}}{2}}-\sqrt{-\frac{\hat{\mathcal{E}}}{2}-\frac{\hat{\mathcal{F}}}{2\sqrt{2\hat{p}_{1}}}-\frac{\hat{p}_{1}}{2}},\hat{r}_{4}=\sqrt{\frac{\hat{p}_{1}}{2}}+\sqrt{-\frac{\hat{\mathcal{E}}}{2}+\frac{\hat{\mathcal{F}}}{2\sqrt{2\hat{p}_{1}}}-\frac{\hat{p}_{1}}{2}},$
(111)

where we have defined,

$\displaystyle\hat{\mathcal{E}}=\frac{3\left(a^{2}-\hat{\eta}-\hat{\lambda}^{2}\right)}{\Lambda\left((a-\hat{\lambda})^{2}+\hat{\eta}\right)+3},\hat{\mathcal{F}}=\frac{6M\left((a-\hat{\lambda})^{2}+\hat{\eta}\right)}{\Lambda\left((a-\hat{\lambda})^{2}+\hat{\eta}\right)+3},\hat{\mathcal{G}}=-\frac{3a^{2}\hat{\eta}}{\Lambda\left((a-\hat{\lambda})^{2}+\hat{\eta}\right)+3}.$
(112)

To avoid crowding equations, we state that the parameter $\hat{p}_{1}$ takes
the same form as eq. 106 and only need to add the hat symbol to $p_{1}$ and
its components.

##  Appendix B  Definition of the elliptic integrals

The functions’ $F(x|k)$, $E(x|k)$, and $\Pi(n,x|k)$, respectively, represent
incomplete elliptic integrals of the first, second, and third kinds, where
$0\leq k<1$ is the elliptic modulus and $x$ is the argument of the integral.
These functions’ equations are as follows:

The incomplete elliptic integral of the first kind, $F(x|k)$, is defined as
follows,

$\displaystyle F(x|k)=\int_{0}^{x}\frac{1}{\sqrt{1-k\sin^{2}t}}\text{dt}.$
(113)

The incomplete elliptic integral of the second kind, $E(x|k)$, is defined as:

$\displaystyle E(x|k)=\int_{0}^{x}\sqrt{1-k\sin^{2}t}\text{dt}.$ (114)

The incomplete elliptic integral of the third kind, $\Pi(n,x|k)$, is defined
as:

$\displaystyle\Pi(n,x|k)=\int_{0}^{x}\frac{1}{(1-n\sin^{2}t)\sqrt{1-k\sin^{2}t}}\text{dt}.$
(115)

where $n\in\mathbb{R}$ is the characteristic of the elliptic integral.

The incomplete elliptic integrals become complete when $x=\frac{\pi}{2}$. As a
result, we have:

$\displaystyle
F\left(\frac{\pi}{2}|k\right)=K(k),E\left(\frac{\pi}{2}|k\right)=E(k),\Pi(n,\frac{\pi}{2}|k)=\Pi(n,k).$
(116)

Moreover, a linear combination of the the elliptic integral of first and third
kind results in the so called associate incomplete elliptic integral of the
third kind [28], $J(n,x|k)$,

$\displaystyle\Pi(n;x|k)-F(x|k)=nJ(n,x|k),$ (117)

where $J(x,n|k)$ is related to Carlson’s $R_{j}$ through the relation,

$\displaystyle
J(n,x|k)=\dfrac{1}{3}\sin^{3}(x)R_{J}\left(\cos^{2}(x),1-k\sin^{2}(x),1,1-n\sin^{2}(x)\right).$
(118)

Carlson’s $R_{j}$ is a Carlson’s elliptic integral of the first kind and is
expressed as;

$\displaystyle
R_{J}(x,y,z,\rho)=\frac{3}{2}\int_{0}^{\infty}\frac{\text{dt}}{(t+x)^{-1/2}(t+\mathcal{Y})^{-1/2}(t+z)^{-1/2}(t+\rho)^{-1}}.$
(119)

##  Appendix C Using Jacobi’s Elliptic Functions to Solve Differential
Equations

In this section we give the general approach that we use to obtain solutions
of the integrals in terms of Jacobi elliptic functions.

###  C.1 Kerr de Sitter

We express the integral form of eq. 4 and (5) as,

$\displaystyle
I_{\theta}=\int_{\theta_{s}}^{\theta_{o}}\frac{\text{d$\theta$}}{\pm_{\theta}\sqrt{\Theta(\theta)}},$
(120) $\displaystyle
I_{r}=\int_{r_{s}}^{r_{o}}\frac{\text{dr}}{\pm_{r}\sqrt{R(r)}},$ (121)
$\displaystyle\phi_{o}-\phi_{s}=\int_{r_{s}}^{r_{o}}\frac{\left(aL^{2}\right)\left(a(a-\lambda)+\text{r}^{2}\right)\text{dr}}{\pm_{r}\text{$\Delta_{r}$}\sqrt{R(r)}}-L^{2}aI_{\phi}+L^{2}\lambda\bar{I}_{\phi},$
(122) $\displaystyle
t_{o}-t_{s}=\int_{r_{s}}^{r_{o}}\dfrac{L^{2}((r^{2}+a^{2})^{2}-a\lambda(a^{2}+r^{2}))\text{dr}}{\pm_{r}\Delta_{r}\sqrt{R(r)}}+a^{2}L^{2}I_{t}-aL^{2}(a-\lambda)I_{\phi},$
(123)

with,

$\displaystyle
I_{\phi}=\int_{\theta_{s}}^{\theta_{o}}\dfrac{\text{d$\theta$}}{\pm_{\theta}\sqrt{\Theta(\theta)}\Delta_{\theta}},\quad\bar{I_{\phi}}=\int_{\theta_{s}}^{\theta_{o}}\dfrac{\text{d$\theta$}}{\pm_{\theta}\sqrt{\Theta(\theta)}\Delta_{\theta}\sin^{2}\theta},\quad
I_{t}=\int_{\theta_{s}}^{\theta_{o}}\dfrac{\cos^{2}\theta\text{d$\theta$}}{\pm_{\theta}\sqrt{\Theta(\theta)}\Delta_{\theta}},$
(124)

where the subscripts $o,s$ denote the observer and the source respectively.
The integral form and Mino parameter approach are related by [13],

$\displaystyle\tau=I_{\theta}=I_{r}.$ (125)

We proceed to obtain general solutions of the angular integrals, eq. 124, as
follows.

To begin with, we express the angular potential in terms of its roots,$u_{+}$
and $u_{-}$ as,

$\displaystyle\Theta(u)=\frac{1}{1-u}[\Xi(u_{+}-u)(u-u_{-})].$ (126)

The parameter $\Xi$ has been defined as,

$\displaystyle\Xi=\left(\frac{1}{3}a^{2}\eta\Lambda+\frac{1}{3}a^{2}\lambda^{2}\Lambda
L^{2}-\frac{2}{3}a^{3}\lambda\Lambda L^{2}+\frac{1}{3}a^{4}\Lambda
L^{2}+a^{2}L^{2}\right).$ (127)

Recall that $u=\cos^{2}\theta$, hence,
$\text{\text{du}}=-2\sqrt{u}\sqrt{1-u}\text{d$\theta$}$. We use this
substitution to transform the integrals,

  1. 1.

$\displaystyle\left|\int_{\pi/2}^{\theta_{i}}\frac{\text{d$\theta$}}{\sqrt{\Theta(\theta)}}\right|=\int_{0}^{u_{i}}\frac{\sqrt{1-u}\text{du}}{2\sqrt{u}\sqrt{1-u}\sqrt{\Xi(u_{+}-u)(u-u_{-})}},$
(128a)
$\displaystyle=\frac{1}{2\sqrt{\Xi}}\int_{0}^{u_{i}}\frac{\text{du}}{\sqrt{u(u_{+}-u)(u-u_{-})}}=\frac{1}{2\sqrt{\Xi}}\int_{0}^{x_{i}}\frac{2u_{+}\cos
x\sin
x\text{dx}}{\sqrt{-u_{-}u_{+}^{2}\sin^{2}x(1-\sin^{2}x)(1-\frac{u_{+}}{u_{-}}\sin^{2}x)}},$
(128b)
$\displaystyle=\dfrac{1}{\sqrt{-u_{-}\Xi}}\int_{0}^{x_{i}}\dfrac{\text{dx}}{\sqrt{1-k\sin^{2}x}}=\dfrac{1}{\sqrt{-u_{-}\Xi}}F\left(x_{i}|k\right).$
(128c)

In eq. 128b we have utilized the relation $u=u_{+}\sin^{2}{x}$. We define
$x_{i}$ and $k$ in eq. 128c as,

$\displaystyle
x_{i}=\arcsin\sqrt{\frac{u_{i}}{u_{+}}}=\arcsin\left[\frac{\cos{\theta_{i}}}{\sqrt{u_{+}}}\right],\quad
k=\frac{u_{+}}{u_{-}}.$ (129)

  2. 2.

$\displaystyle\left|\int_{\pi/2}^{\theta_{i}}\frac{\text{d$\theta$}}{\Delta_{\theta}\sqrt{\Theta_{\theta}}}\right|=\dfrac{1}{\sqrt{-u_{-}\Xi}}\int_{0}^{x_{i}}\dfrac{\text{dx}}{(1+u_{+}\mathcal{Y}\sin^{2}x)\sqrt{1-k\sin^{2}x}}=\dfrac{1}{\sqrt{-u_{-}\Xi}}\Pi\left(-u_{+}\mathcal{Y};x_{i}|k\right).$
(130)

  3. 3.

$\displaystyle\left|\int_{\pi/2}^{\theta_{i}}\dfrac{\text{d$\theta$}}{\sqrt{\Theta(\theta)}\Delta_{\theta}\sin^{2}\theta}\right|=\dfrac{1}{\sqrt{-u_{-}\Xi}}\int_{0}^{x_{i}}\dfrac{\text{dx}}{(1-u_{+}\sin^{2}x)(1+u_{+}\mathcal{Y}\sin^{2}x)\sqrt{1-k\sin^{2}x}},$
(131a)
$\displaystyle=\dfrac{1}{\sqrt{-u_{-}\Xi}}\int_{0}^{t_{i}}\dfrac{\text{dt}}{(1-u_{+}t^{2})(1+u_{+}\mathcal{Y}t^{2})\sqrt{(1-t^{2})(1-kt^{2})}},$
(131b)
$\displaystyle=\dfrac{1}{\sqrt{-u_{-}\Xi}}\dfrac{1}{1+\mathcal{Y}}[\Pi(u_{+};\arcsin
t_{i}|k)+\mathcal{Y}\Pi(-u_{+}\mathcal{Y};\arcsin t_{i}|k)].$ (131c)

Where $t_{i}$ is given by the relation,

$\displaystyle t_{i}=\sqrt{\frac{u_{i}}{u_{+}}}.$ (132)

  4. 4.

$\displaystyle\left|\int_{\pi/2}^{\theta_{i}}\frac{\cos^{2}\theta\text{d$\theta$}}{\Delta_{\theta}\sqrt{\Theta_{\theta}}}\right|=\dfrac{1}{\sqrt{-u_{-}\Xi}}\int_{0}^{x_{i}}\dfrac{u_{+}\sin^{2}x\text{dx}}{(1+u_{+}\mathcal{Y}\sin^{2}x)\sqrt{1-k\sin^{2}x}}$
(133a)
$\displaystyle=\dfrac{u_{+}}{\sqrt{-u_{-}\Xi}}\int_{0}^{t_{i}}\dfrac{t^{2}\text{dt}}{(1+u_{+}\mathcal{Y}t^{2})\sqrt{(1-t^{2})\left(1-kt^{2}\right)}}=\dfrac{3}{a^{2}\Lambda\sqrt{-u_{-}\Xi}}(F(\arcsin
t_{i}|k)-\Pi(-u_{+}\mathcal{Y};\arcsin t_{i}|k)).$ (133b)

We have observed the equivalence of $F(x|k)$ and $\Pi(n;x|k)$ when $n=0$.
Thus, when the value of $n$ is extremely small, this combination may not be
suitable due to potential round-off errors. For example, in eq. 133b,
considering the astrophysically relevant value of the cosmological constant,
$\mathcal{Y}=\frac{a^{2}\Lambda}{3}$ becomes very small. Consequently, in our
numerical computations, the term $(F(\arcsin
t_{i}|k)-\Pi(\frac{-u_{+}a^{2}\Lambda}{3};\arcsin t_{i}|k))$ suffers round off
errors arising from this smallness. To circumvent the round off errors, we
employ an alternative representation of this linear combination using the so
called associate incomplete elliptic integral of third kind, $J(n,x|k)$ [28],

$\displaystyle\dfrac{3}{a^{2}\Lambda\sqrt{-u_{-}\Xi}}(F(\arcsin
t_{i}|k)-\Pi(-u_{+}\mathcal{Y};\arcsin
t_{i}|k))=\frac{u_{+}J(-u_{+}\mathcal{Y},\arcsin t_{i}|k)}{\sqrt{-u_{-}\Xi}}.$
(134)

The angular component of the photon trajectories can then be determined
through utilizing the above solutions by substituting the subscript $i$ with
either $s$ or $o$, representing the source or observer respectively, along
with the corresponding transformed coordinate.

To derive the solutions for the integrals used in determining critical
parameters, we adopt the following approach. Critical parameters are evaluated
on the critical radius which is the radial coordinate of bound photon orbits.
As a consequence, the integral over the radial coordinate is treated as a
constant. Moreover, these parameters are defined over half orbits. A photon is
deemed to have completed half an orbit around the black hole when it
transitions between its maximum and minimum inclinations with respect to the
equatorial plane, or vice versa. This corresponds to integrating the angular
integrals from one turning point to the next. Our analysis focuses on ordinary
motion, leading us to evaluate the integrals for half orbits over $\theta_{1}$
to $\theta_{4}$ or vice versa.

As a result, the change in azimuthal angle parameter and time delay will be
obtained by evaluating the integrals in eq. 122 and eq. 123 over half orbits
and at constant $r$,

$\displaystyle\delta=\frac{\left(aL^{2}\right)\left(a(a-\lambda)+\text{r}^{2}\right)}{\text{$\Delta_{r}$}}\left|\int_{\theta_{1}}^{\theta_{4}}\frac{\text{d$\theta$}}{\sqrt{\Theta(\theta)}}\right|-L^{2}a\left|\int_{\theta_{1}}^{\theta_{4}}\frac{\text{d$\theta$}}{\Delta_{\theta}\sqrt{\Theta_{\theta}}}\right|+L^{2}\lambda\left|\int_{\theta_{1}}^{\theta_{4}}\dfrac{\text{d$\theta$}}{\sqrt{\Theta(\theta)}\Delta_{\theta}\sin^{2}\theta}\right|,$
(135)
$\displaystyle\zeta=\dfrac{L^{2}((r^{2}+a^{2})^{2}-a\lambda(a^{2}+r^{2}))}{\Delta_{r}}\left|\int_{\theta_{1}}^{\theta_{4}}\frac{\text{d$\theta$}}{\sqrt{\Theta(\theta)}}\right|+a^{2}L^{2}\left|\int_{\theta_{1}}^{\theta_{4}}\frac{\cos^{2}\theta\text{d$\theta$}}{\Delta_{\theta}\sqrt{\Theta_{\theta}}}\right|-aL^{2}(a-\lambda)\left|\int_{\theta_{1}}^{\theta_{4}}\frac{\text{d$\theta$}}{\Delta_{\theta}\sqrt{\Theta_{\theta}}}\right|.$
(136)

Due to the constancy of the integral over the radial coordinate, we have
employed eq. 125 to replace the Mino parameter with the angular integral in
the first part of eq. 135 and (136).

Utilizing eq. 128c and the definition of $x_{i}$ in eq. 129, we have that,

$\displaystyle\left|\int_{\theta_{1}}^{\theta_{4}}\frac{\text{d$\theta$}}{\sqrt{\Theta(\theta)}}\right|=\left|\dfrac{1}{\sqrt{-u_{-}\Xi}}F\left(\arcsin\left[\frac{\cos{\theta_{4}}}{\sqrt{u_{+}}}\right]|k\right)-\dfrac{1}{\sqrt{-u_{-}\Xi}}F\left(\arcsin\left[\frac{\cos{\theta_{1}}}{\sqrt{u_{+}}}\right]|k\right)\right|,$
(137a)
$\displaystyle=\left|-\dfrac{1}{\sqrt{-u_{-}\Xi}}F\left(\arcsin\left[\frac{\sqrt{u_{+}}}{\sqrt{u_{+}}}\right]|k\right)-\dfrac{1}{\sqrt{-u_{-}\Xi}}F\left(\arcsin\left[\frac{\sqrt{u_{+}}}{\sqrt{u_{+}}}\right]|k\right)\right|,$
(137b)
$\displaystyle=\dfrac{2}{\sqrt{-u_{-}\Xi}}F\left(\frac{\pi}{2}|k\right)=\dfrac{2}{\sqrt{-u_{-}\Xi}}K(k)=\Gamma_{\theta}.$
(137c)

In eq. 137b, we have used the definitions from eq. 11 and (14) where
$\sqrt{u_{+}}=\cos{\theta_{1}}$ and $-\sqrt{u_{+}}=\cos{\theta_{4}}$.
Furthermore, we have made use of the property, $F(-x\mid k)=-F(x\mid k)$ and
$\arcsin{1}=\pi/2$. If we still integrated from $\theta_{4}$ to $\theta_{1}$,
the result will remain the same.

For the subsequent integrals, we make use of the general solutions eq.
130-(134), and directly substitute
$\arcsin\left[\frac{\cos{\theta_{4}}}{\sqrt{u_{+}}}\right]=-\frac{\pi}{2}$ and
$\arcsin\left[\frac{\cos{\theta_{1}}}{\sqrt{u_{+}}}\right]=\frac{\pi}{2}$ as
we have explained how this comes up which transforms the incomplete integrals
to be complete. This gives,

$\displaystyle\left|\int_{\theta_{1}}^{\theta_{4}}\frac{\text{d$\theta$}}{\Delta_{\theta}\sqrt{\Theta_{\theta}}}\right|=\dfrac{2}{\sqrt{-u_{-}\Xi}}\Pi\left(-u_{+}\mathcal{Y};|k\right)=\Gamma_{\phi},$
(138)
$\displaystyle\left|\int_{\theta_{1}}^{\theta_{4}}\dfrac{\text{d$\theta$}}{\sqrt{\Theta(\theta)}\Delta_{\theta}\sin^{2}\theta}\right|=\dfrac{2}{\sqrt{-u_{-}\Xi}}\dfrac{1}{1+\mathcal{Y}}[\Pi(u_{+};|k)+\mathcal{Y}\Pi(-u_{+}\mathcal{Y};|k)]=\bar{\Gamma}_{\phi},$
(139)
$\displaystyle\left|\int_{\theta_{1}}^{\theta_{4}}\frac{\cos^{2}\theta\text{d$\theta$}}{\Delta_{\theta}\sqrt{\Theta_{\theta}}}\right|=\frac{2u_{+}J(-u_{+}\mathcal{Y},|k)}{\sqrt{-u_{-}\Xi}}=\Gamma_{t}.$
(140)

Equation 137c-(140) will then be inserted into (135) and (136) to obtain the
change in azimuthal angle and time delay.

The Lyapunov exponent will be calculated in the following way. The radial
coordinate of a nearly bound photon orbit is very close to the radial
coordinate of a bound photon orbit, with just a small difference between them.
We therefore write the radial coordinate of the nearly bound orbit as
$r=r_{b}+\delta r$, where $r_{b}$ is the radial coordinate of the bound photon
orbit and $0<\delta r<1$ is a small deviation from this radial coordinate.
Performing a Taylor expansion on the radial potential $R(r)$ with respect to
this deviation results in,

$\displaystyle
R(r)=R(r_{b})+(r-r_{b})R^{\prime}(r_{b})+\frac{1}{2}(r-r_{b})^{2}R^{\prime\prime}(r_{b})+\frac{1}{6}(r-r_{b})^{3}R^{(3)}(r_{b})+O\left((r-r_{b})^{4}\right).$
(141)

Equation 141 is an approximation of the radial potential for a nearly bound
photon orbit. We then proceed by keeping terms up to second order in the
deviation $\delta r$ because higher order terms will become smaller as we move
away from the bound orbit. Further, recall that for a bound photon orbit,
$R(r)=R^{\prime}(r)=0$, hence the first two terms of eq. 141 will vanish
leaving us with only the third term. Let $\delta r_{1}$ be the initial
deviation and $\delta r_{n}$ be the deviation after $n$ half orbits. Thus
after $n$ half orbits, we have that,

$\displaystyle\int_{r_{b}+\delta r_{1}}^{r_{b}+\delta
r_{n}}\frac{\text{dr}}{\sqrt{R(r)}}\approx\int_{r_{b}+\delta
r_{1}}^{r_{b}+\delta
r_{n}}\frac{\sqrt{2}\text{dr}}{\sqrt{R^{\prime\prime}(r_{b})(r-r_{b})^{2}}}=\frac{\sqrt{2}(\ln{(\delta
r_{n})}-\ln{(\delta r_{1})})}{\sqrt{R^{\prime\prime}(r_{b})}}.$ (142)

The Mino period for half an orbit in the latitudinal direction has been
defined in eq. 137c. Therefore for $n$ half orbits the photon will have a Mino
period of $n\Gamma_{\theta}$. Equation 125 gives us the equality of the radial
and latitudinal integral and by this definition we have that after $n$ half
orbits,

$\displaystyle\frac{\sqrt{2}(\ln{(\delta r_{n})}-\ln{(\delta
r_{1})})}{\sqrt{R^{\prime\prime}(r_{b})}}=n\Gamma_{\theta}.$ (143)

Simplifying eq. 143 gives,

$\displaystyle\delta r_{n}=\delta
r_{1}\exp{\left(\sqrt{\frac{R^{\prime\prime}(r_{b})}{2}}n\Gamma_{\theta}\right)}.$
(144)

Equation 144 tells us that after $n$ half orbits, a nearly bound photon orbit
exponentially deviates from a bound orbit at a rate given by,

$\displaystyle\gamma=\sqrt{\frac{R^{\prime\prime}(r_{b})}{2}}\Gamma_{\theta}.$
(145)

Thus, eq. 145 defines the Lyapunov exponent per half orbit. A similar relation
has been obtained in [29] in the case of a Kerr black hole.
For nearly bound photon orbits we also need to obtain the general solutions of
the radial integrals. Firstly, we express the radial potential and the
parameter $\Delta_{r}$ in terms of their roots,

$\displaystyle
R(r)=(r-r_{1})(r-r_{2})(r-r_{3})(r-r_{4}),\Delta_{r}=-\frac{3}{\Lambda}(r-r_{-c})(r-r_{-})(r-r_{+})(r-r_{c}).$
(146)

Where $r_{-c},r_{-},r_{+},r_{c}$ are the kds horizons. Applying the notation
eq. 146,we re-express eq. 121\- (123) as,

$\displaystyle
I_{r}=\int_{r_{s}}^{r_{o}}\frac{\text{dr}}{\pm\sqrt{(r-r_{1})(r-r_{2})(r-r_{3})(r-r_{4})}},$
(147)
$\displaystyle\phi_{o}-\phi_{s}=-\frac{3aL^{2}}{\Lambda}\left(\frac{\Omega_{c}}{r_{c,-c}r_{c,-}r_{c,+}}-\frac{\Omega_{-c}}{r_{c,-c}r_{-c,-}r_{-c,+}}-\frac{\Omega_{-}}{r_{c,-}r_{-,-c}r_{-,+}}-\frac{\Omega_{+}}{r_{c,+}r_{+,-c}r_{+,-}}\right)-L^{2}aI_{\phi}+L^{2}\lambda\bar{I}_{\phi},$
(148) $\displaystyle
t_{o}-t_{s}=\frac{3L^{2}}{\Lambda}\left(-\frac{\Phi_{c}}{r_{c,-c}r_{c,-}r_{c,+}}+\frac{\Phi_{-c}}{r_{c,-c}r_{-c,-}r_{-c,+}}+\frac{\Phi_{-}}{r_{c,-}r_{-,-c}r_{-,+}}+\frac{\Phi_{+}}{r_{c,+}r_{+,-c}r_{+,-}}-\tau\right)+a^{2}L^{2}I_{t}-aL^{2}(a-\lambda)I_{\phi},$
(149)
$\displaystyle\Phi_{x}=\left(a^{2}+r_{x}^{2}\right)\Omega_{x},\Omega_{x}=\left(a^{2}-a\lambda+r_{x}^{2}\right)I_{x},I_{x}=\int_{r_{s}}^{r_{o}}\frac{\text{dr}}{\pm_{r}(r-r_{x})\sqrt{R(r)}}.$
(150)

The subscripts $x$ correspond to the horizons $-c,-,+,c$, where $r_{c}$,
$r_{-}$, $r_{+}$, and $r_{-c}$ represent the cosmological horizon, Cauchy
horizon, event horizon, and the dual of the cosmological horizon,
respectively. Additionally, $r_{a,b}=r_{a}-r_{b}$. We will solve the radial
integrals $I_{r}$ and $I_{x}$ for two scenarios; when the radial potential has
four real roots and that which when the radial potential has only two real
roots with the other two being complex conjugates. The solutions will be
obtained by utilizing the substitutions in Ref.[30].

  * •

$r_{1}<r_{2}<r_{3}<r_{4}<r_{s}$
Utilizing part $258.00$ of Ref. [30], an integral of the type $I_{r}$ has a
solution of the form,

$\displaystyle\int_{r_{4}}^{r_{i}}\frac{\text{dr}}{\sqrt{(r-r_{1})(r-r_{2})(r-r_{3})(r-r_{4})}}=g_{E}F(\varphi_{E,i}\mid
k_{E}),$ (151)

where $i$ denotes the source $s$ or observer $o$ and,

$\displaystyle
k_{E}=\frac{(r_{4}-r_{1})(r_{3}-r_{2})}{(r_{4}-r_{2})(r_{3}-r_{1})},\quad
g_{E}=\frac{2}{\sqrt{(r_{4}-r_{2})(r_{3}-r_{1})}},\quad\varphi_{E,i}=\arcsin\sqrt{\frac{(r_{3}-r_{1})(r_{i}-r_{4})}{(r_{4}-r_{1})(r_{i}-r_{3})}}.$
(152)

Furthermore, making use of part $258.39$ and $340.01$ of Ref. [30], the
solution to the integral of the form $I_{x}$ is expressed as,

$\displaystyle\int_{r_{4}}^{r_{i}}\frac{\text{dr}}{(r-r_{x})\sqrt{(r-r_{1})(r-r_{2})(r-r_{3})(r-r_{4})}}=-\frac{2[(r_{4}-r_{3})\Pi_{i}-(r_{x}-r_{4})F_{i}]}{(r_{x}-r_{4})(r_{x}-r_{3})\sqrt{(r_{4}-r_{2})(r_{3}-r_{1})}},$
(153)

where,

$\displaystyle\Pi_{i}=\Pi\left(\frac{(r_{4}-r_{1})(r_{x}-r_{3})}{(r_{x}-r_{4})(r_{3}-r_{1})};\varphi_{E,i}|k_{E}\right),F_{i}=F\left(\varphi_{E,i}|k_{E}\right).$
(154)

  * •

$r_{1}<r_{2}<r_{s},r_{3}=\bar{r}_{4}$
For these nature of roots, the appropriate substitution to use in solving
integrals of the form $I_{r}$ and $I_{x}$ in Ref.[30] are parts $260.00$ and
$260.04$ respectively. Using these substitutions we obtain,

$\displaystyle\int_{r_{2}}^{r_{i}}\frac{\text{dr}}{\sqrt{(r-r_{1})(r-r_{2})(r-r_{3})(r-r_{4})}}=g_{P}F(\varphi_{P,i}\mid
k_{P}),$ (155)
$\displaystyle\int_{r_{2}}^{r_{i}}\frac{\text{dr}}{(r-r_{x})\sqrt{(r-r_{1})(r-r_{2})(r-r_{3})(r-r_{4})}}=\frac{(g_{P}(B-A))\left((\alpha-\alpha_{2})R_{1,i}+\alpha_{2}F(\varphi_{P,i}\mid
k_{P})\right)}{r_{2}B+r_{1}A-Ar_{x}-Br_{x}},$ (156)

where the parameters $A,g_{P},k_{P}$ and $B$ have been defined as,

$\displaystyle
A=\sqrt{(r_{2}-b_{1}^{2})^{2}+a_{1}^{2}},B=\sqrt{(r_{1}-b_{1})^{2}+a_{1}^{2}},\quad
g_{P}=\frac{1}{\sqrt{AB}},k_{P}=\frac{(A+B)^{2}-(r_{2}-r_{1})^{2}}{4AB},$
(157) $\displaystyle a_{1}=\sqrt{-\frac{(r_{3}-r_{4})^{2}}{4}},\quad
b_{1}=\frac{r_{3}+r_{4}}{2},$ (158)

we have also defined,

$\displaystyle
R_{1,i}=\frac{1}{1-\alpha^{2}}\left(\Pi\left[\frac{\alpha^{2}}{\alpha^{2}-1};\varphi_{P,i}\mid
k_{P}\right]-\alpha
f_{1,i}\right),\quad\alpha=\frac{Br_{2}+Ar_{1}-r_{x}(A+B)}{Br_{2}-r_{1}A+r_{x}(A+B)},\quad\alpha_{2}=\frac{A+B}{B-A},$
(159) $\displaystyle
f_{1,i}=\sqrt{\frac{\alpha^{2}-1}{k_{P}+(1-k_{P})\alpha^{2}}}\ln{\left|\frac{\sqrt{k_{P}+(1-k_{P})\alpha^{2}}dn(u_{i}\mid
k_{P})+\sqrt{\alpha^{2}-1}sn(u_{i}\mid
k_{P})}{\sqrt{k_{P}+(1-k_{P})\alpha^{2}}dn(u_{i}\mid
k_{P})-\sqrt{\alpha^{2}-1}sn(u_{i}\mid k_{P})}\right|},$ (160) $\displaystyle
u_{i}=cn^{-1}\left[\cos(\varphi_{P,i})\mid
k_{P}\right],\varphi_{P,i}=\arccos\left[\frac{(A-B)r_{i}+r_{2}B-r_{1}A}{(A+B)r_{i}-r_{2}B-r_{1}A}\right].$
(161)

The complete solutions for radial integrals will be derived using the above
solutions by substituting $i$ with $s$ or $o$ which denotes the source or the
observer respectively.

###  C.2 Kerr de Sitter Revisited

Equation 49 and (50) can then be expressed in integral forms as,

$\displaystyle\hat{I}_{\hat{\theta}}=\int_{\hat{\theta}_{s}}^{\hat{\theta}_{o}}\frac{\text{d$\theta$}}{\sqrt{\hat{\Theta}(\hat{\theta})}},$
(162)
$\displaystyle\hat{I}_{r}=\int_{r_{s}}^{r_{o}}\frac{\text{dr}}{\sqrt{\hat{R}(r)}},$
(163)
$\displaystyle\hat{\phi}_{o}-\hat{\phi}_{s}=\int_{r_{s}}^{r_{o}}\left(\frac{ar^{2}+a^{3}-a\hat{\Delta}-a^{2}\hat{\lambda}}{\hat{\Delta}}\right)\frac{\text{dr}}{\sqrt{\hat{R}(r)}}+\hat{\lambda}\hat{I}_{\phi},$
(164) $\displaystyle
t_{o}-t_{s}=\int_{r_{s}}^{r_{o}}\left(\frac{(r^{2}+a^{2})(r^{2}+a^{2}-a\lambda)}{\hat{\Delta}}+a\hat{\lambda}-a^{2}\right)\frac{\text{dr}}{\sqrt{\hat{R}(r)}}+a^{2}\hat{I}_{t},$
(165)

for which we have defined,

$\displaystyle\hat{I}_{\phi}=\int_{\theta_{s}}^{\theta_{o}}\frac{\text{d$\theta$}}{\sin^{2}\theta\sqrt{\hat{\Theta}(\theta)}},\hat{I}_{t}=\int_{\theta_{s}}^{\theta_{o}}\frac{\cos^{2}\theta\text{d$\theta$}}{\sqrt{\hat{\Theta}(\theta)}}.$
(166)

The subscripts $s$ and $o$ denotes the source and observer respectively. To
obtain the solutions of eq. 166, we proceed by first expressing the angular
potential in terms of its roots $\hat{u_{+}}$ and $\hat{u_{-}}$,

$\displaystyle\hat{\Theta}_{\hat{u}}=\frac{a^{2}}{1-\hat{u}}(\hat{u}_{+}-\hat{u})(\hat{u}-\hat{u}_{-}).$
(167)

The angular integrals still unpack as eq. 193. We obtain the general solutions
as,

  1. 1.

$\displaystyle\left|\int_{\pi/2}^{\theta_{i}}\frac{\text{d$\theta$}}{\sqrt{\hat{\Theta}(\theta)}}\right|=\dfrac{1}{\sqrt{-\hat{u}_{-}a^{2}}}F\left(\hat{x}_{i}|\hat{k}\right),$
(168)

where,

$\displaystyle\hat{x}_{i}=\arcsin\left(\sqrt{\frac{\cos^{2}\theta_{i}}{\hat{u}_{+}}}\right),\hat{k}=\frac{\hat{u}_{+}}{\hat{u}_{-}}.$
(169)

  2. 2.

$\displaystyle\left|\int_{\pi/2}^{\theta_{i}}\frac{\text{d$\theta$}}{\sin^{2}\theta\sqrt{\hat{\Theta}(\theta)}}\right|=\dfrac{1}{\sqrt{-u_{-}a^{2}}}\int_{0}^{\hat{x}_{i}}\dfrac{\text{dx}}{(1-u_{+}\sin^{2}x)\sqrt{1-k\sin^{2}x}}=\dfrac{1}{\sqrt{-u_{-}a^{2}}}\Pi[\hat{u}_{+};\hat{x}_{i}|\hat{k}].$
(170)

  3. 3.