Black hole lensing in Kerr-de Sitter spacetimes

Abstract

We have derived analytical solutions using Jacobi elliptic functions for bound and nearly bound photon orbits in Kerr-de Sitter and Kerr-de Sitter revisited 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 general relativity. 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.

Appendices

Appendix A: Calculating roots using Ferrari’s method

The radial potential, Eq. (6), can be re-expressed as,

$$\begin{aligned} R(r)= r^{4}+{\mathcal {E}}r^{2}+{\mathcal {F}}r+{\mathcal {G}}, \end{aligned}$$

(A1)

for which we have defined the coefficients as,

$$\begin{aligned} {\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{6 M \left( L^2 (a-\lambda )^2+\eta \right) }{L^2 \left( \Lambda (a-\lambda )^2+3\right) +\eta \Lambda },{\mathcal {G}}=-\frac{3 a^2 \eta }{L^2 \left( \Lambda (a-\lambda )^2+3\right) +\eta \Lambda }. \end{aligned}$$

(A2)

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

$$\begin{aligned} 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) }, \end{aligned}$$

(A3)

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

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

(A4)

A parameter p is then introduced into Eq. (A4) by adding $2 r^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,

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

(A5)

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. (A5) and equating to zero results in,

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

(A6)

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

$$\begin{aligned} 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} \end{aligned}$$

(A7)

Cardano’s method [27] gives the roots of Eq. (A7), whereby we judiciously select the real root,

$$\begin{aligned} v={\mathcal {Z}}_++{\mathcal {Z}}_-, {\mathcal {Z}}_{\pm }=\root 3 \of {\pm \sqrt{\frac{{\mathcal {U}}^{2}}{4}+\frac{{\mathcal {T}}^{3}}{27}}-\frac{{\mathcal {U}}}{2}}. \end{aligned}$$

(A8)

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

$$\begin{aligned} p_1={\mathcal {Z}}_++{\mathcal {Z}}_--{\mathcal {E}}/3, \end{aligned}$$

(A9)

Substituting Eq. (A9) into (A3) and distributing the signs $\pm _1$/$\pm _2$ results in the four roots of the radial potential,

$$\begin{aligned} r_{1}&= -\sqrt{\frac{p_{1}}{2}}-\sqrt{-\frac{{\mathcal {E}}}{2}+\frac{{\mathcal {F}}}{2\sqrt{2 p_1}}-\frac{p_{1}}{2}}, r_{2}=-\sqrt{\frac{p_{1}}{2}}+\sqrt{-\frac{{\mathcal {E}}}{2}+\frac{{\mathcal {F}}}{2\sqrt{2 p_1}}-\frac{p_{1}}{2}}, \end{aligned}$$

(A10)

$$\begin{aligned} r_{3}&= -\sqrt{\frac{p_{1}}{2}}-\sqrt{-\frac{{\mathcal {E}}}{2}-\frac{{\mathcal {F}}}{2\sqrt{2 p_1}}-\frac{p_{1}}{2}}, r_{4}= \sqrt{\frac{p_{1}}{2}}+\sqrt{-\frac{{\mathcal {E}}}{2}+\frac{{\mathcal {F}}}{2\sqrt{2 p_1}}-\frac{p_{1}}{2}}. \end{aligned}$$

(A11)

The roots in Eq. (A11) can be analyzed using the discriminant of Eq. (A1),

$$\begin{aligned} \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. \end{aligned}$$

(A12)

$\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,

$$\begin{aligned} {\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}}, \end{aligned}$$

(A13)

$$\begin{aligned} {\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}}, \end{aligned}$$

(A14)

where we have defined,

$$\begin{aligned} \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{6 M \left( (a-{\hat{\lambda }} )^2+{\hat{\eta }} \right) }{\Lambda \left( (a-{\hat{\lambda }} )^2+{\hat{\eta }} \right) +3},\hat{{\mathcal {G}}}=-\frac{3 a^2 {\hat{\eta }} }{\Lambda \left( (a-{\hat{\lambda }} )^2+{\hat{\eta }}\right) +3}. \end{aligned}$$

(A15)

To avoid crowding equations, we state that the parameter ${\hat{p}}_1$ takes the same form as Eq. (A9) 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\le 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,

$$\begin{aligned} F(x|k) = \int _0^x \frac{1}{\sqrt{1-k\sin ^2 t}} \text {dt}. \end{aligned}$$

(B1)

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

$$\begin{aligned} E(x|k) = \int _0^x \sqrt{1-k\sin ^2 t} \text {dt}. \end{aligned}$$

(B2)

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

$$\begin{aligned} \Pi (n,x|k) = \int _0^x \frac{1}{(1-n\sin ^2 t)\sqrt{1-k\sin ^2 t}} \text {dt}. \end{aligned}$$

(B3)

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:

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

(B4)

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),

$$\begin{aligned} \Pi (n;x|k)-F(x|k)= n J(n,x | k), \end{aligned}$$

(B5)

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

$$\begin{aligned} 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) . \end{aligned}$$

(B6)

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

$$\begin{aligned} 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}}. \end{aligned}$$

(B7)

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.

1.1 1. Kerr de Sitter

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

$$\begin{aligned} I_{\theta }&= \int _{\theta _{s}}^{\theta _{o}} \frac{\text {d}\theta }{\pm _{\theta } \sqrt{\Theta (\theta )}} , \end{aligned}$$

(C1)

$$\begin{aligned} I_{r}&= \int _{r_{s}}^{r_{o}} \frac{\text {dr}}{\pm _{r} \sqrt{R(r)}}, \end{aligned}$$

(C2)

$$\begin{aligned} \phi _{o}-\phi _{s}&= \int _{r_{s}}^{r_{o}} \frac{\left( a L^2\right) \left( a (a-\lambda )+\text {r}^2\right) \text {dr}}{\pm _{r} {\Delta _{r} } \sqrt{R(r)}} -L^{2}a I_{\phi }+ L^{2}\lambda \bar{I}_{\phi }, \end{aligned}$$

(C3)

$$\begin{aligned} 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}-a L^2(a-\lambda )I_{\phi }, \end{aligned}$$

(C4)

with,

$$\begin{aligned} 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 }}, \end{aligned}$$

(C5)

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

$$\begin{aligned} \tau =I_{\theta }=I_{r}. \end{aligned}$$

(C6)

We proceed to obtain general solutions of the angular integrals, Eq. (C5), as follows.

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

$$\begin{aligned} \Theta (u)=\frac{1}{1-u}[\Xi (u_{+}-u)(u-u_{-})]. \end{aligned}$$

(C7)

The parameter $\Xi$ has been defined as,

$$\begin{aligned} \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) . \end{aligned}$$

(C8)

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

• 1.
  $$\begin{aligned} \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_-)}}, \end{aligned}$$

  (C9a)
  $$\begin{aligned}&= \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{2 u_{+}\cos x\sin x \text {dx}}{\sqrt{-u_{-}u_{+}^{2}\sin ^{2}x(1-\sin ^{2}x)(1-\frac{u_{+}}{u_{-}}\sin ^{2}x)}}, \end{aligned}$$

  (C9b)
  $$\begin{aligned}&= \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) . \end{aligned}$$

  (C9c)

  In Eq. (C9b) we have utilized the relation $u=u_{+}\sin ^2{x}$. We define $x_{i}$ and k in Eq. (C9c) as,

  $$\begin{aligned} x_i=\arcsin \sqrt{\frac{u_i}{u_+}}=\arcsin \left[ \frac{\cos {\theta _{i}}}{\sqrt{u_{+}}}\right] , \quad k=\frac{u_{+}}{u_{-}}. \end{aligned}$$

  (C10)
• 2.
  $$\begin{aligned} \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) . \end{aligned}$$

  (C11)
• 3.
  $$\begin{aligned} \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 }}, \end{aligned}$$

  (C12a)
  $$\begin{aligned}&=\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}) }}, \end{aligned}$$

  (C12b)
  $$\begin{aligned}&=\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) ] . \end{aligned}$$

  (C12c)

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

  $$\begin{aligned} t_{i}=\sqrt{\frac{u_{i}}{u_{+}}}. \end{aligned}$$

  (C13)
• 4.
  $$\begin{aligned} \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 }} \end{aligned}$$

  (C14a)
  $$\begin{aligned}&= \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) }} \nonumber \\&= \dfrac{3}{a^2 \Lambda \sqrt{-u_{-}\Xi }} (F(\arcsin t_i|k)-\Pi ( -u_+ {\mathcal {Y}};\arcsin t_i|k)). \end{aligned}$$

  (C14b)

  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. (C14b), 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],

  $$\begin{aligned} \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 }}. \end{aligned}$$

  (C15)

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. (C3) and Eq. (C4) over half orbits and at constant r,

$$\begin{aligned} \delta&= \frac{\left( a L^2\right) \left( a (a-\lambda )+\text {r}^2\right) }{ {\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| , \end{aligned}$$

(C16)

$$\begin{aligned} \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| -a L^2(a-\lambda ) \left| \int _{\theta _1}^{\theta _4} \frac{\text {d}\theta }{\Delta _ {\theta } \sqrt{\Theta _{\theta }}}\right| . \end{aligned}$$

(C17)

Due to the constancy of the integral over the radial coordinate, we have employed Eq. (C6) to replace the Mino parameter with the angular integral in the first part of Eq. (C16) and (C17).

Utilizing Eq. (C9c) and the definition of $x_i$ in Eq. (C10), we have that,

$$\begin{aligned} \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| , \end{aligned}$$

(C18a)

$$\begin{aligned}&= \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| , \end{aligned}$$

(C18b)

$$\begin{aligned}&=\dfrac{2}{\sqrt{-u_{-}\Xi }} F\left( \frac{\pi }{2}| k \right) =\dfrac{2}{\sqrt{-u_{-}\Xi }} K(k)=\Gamma _\theta . \end{aligned}$$

(C18c)

In Eq. (C18b), we have used the definitions from Eqs. (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 Eqs. (C11–C15), 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,

$$\begin{aligned} \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 , \end{aligned}$$

(C19)

$$\begin{aligned} \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 , \end{aligned}$$

(C20)

$$\begin{aligned} \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. \end{aligned}$$

(C21)

Equations (C18c–C21) will then be inserted into (C16) and (C17) 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,

$$\begin{aligned} R(r)=R(r_b)+(r-r_b) R'(r_b)+\frac{1}{2} (r-r_b)^2 R''(r_b)+\frac{1}{6} (r-r_b)^3 R^{(3)}(r_b)+O\left( (r-r_b)^4\right) . \end{aligned}$$

(C22)

Equation (C22) 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'(r)=0$, hence the first two terms of Eq. (C22) 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,

$$\begin{aligned} \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''(r_b)(r-r_b)^2}}=\frac{\sqrt{2}(\ln {(\delta r_n)}-\ln {(\delta r_1)})}{\sqrt{R''(r_b)}}. \end{aligned}$$

(C23)

The Mino period for half an orbit in the latitudinal direction has been defined in Eq. (C18c). Therefore for n half orbits the photon will have a Mino period of $n \Gamma _\theta$. Equation (C6) gives us the equality of the radial and latitudinal integral and by this definition we have that after n half orbits,

$$\begin{aligned} \frac{\sqrt{2}(\ln {(\delta r_n)}-\ln {(\delta r_1)})}{\sqrt{R''(r_b)}}=n \Gamma _\theta . \end{aligned}$$

(C24)

Simplifying Eq. (C24) gives,

$$\begin{aligned} \delta r_n=\delta r_1 \exp {\left( \sqrt{\frac{R''(r_b)}{2}} n \Gamma _\theta \right) }. \end{aligned}$$

(C25)

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

$$\begin{aligned} \gamma =\sqrt{\frac{R''(r_b)}{2}} \Gamma _\theta . \end{aligned}$$

(C26)

Thus, Eq. (C26) 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,

$$\begin{aligned} 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}). \end{aligned}$$

(C27)

Where $r_{-c},r_{-},r_{+},r_{c}$ are the KdS horizons. Applying the notation Eq. (C27),we re-express Eqs. (C2–C4) as,

$$\begin{aligned} I_{r}&= \int _{r_{s}}^{r_{o}} \frac{\text {dr}}{\pm \sqrt{(r-r_{1})(r-r_{2})(r-r_{3})(r-r_{4})}}, \end{aligned}$$

(C28)

$$\begin{aligned} \phi _{o}-\phi _{s}&=-\frac{3 a L^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}a I_{\phi }+ L^{2}\lambda \bar{I}_{\phi } , \end{aligned}$$

(C29)

$$\begin{aligned} t_{o}-t_{s}&=\frac{3 L^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}-a L^2(a-\lambda )I_{\phi }, \end{aligned}$$

(C30)

$$\begin{aligned} \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)}}. \end{aligned}$$

(C31)

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].

• $\underline{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,

  $$\begin{aligned} \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), \end{aligned}$$

  (C32)

  where i denotes the source s or observer o and,

  $$\begin{aligned} 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})}}. \end{aligned}$$

  (C33)

  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,

  $$\begin{aligned} \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})}}, \end{aligned}$$

  (C34)

  where,

  $$\begin{aligned} \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) . \end{aligned}$$

  (C35)

• $\underline{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,

  $$\begin{aligned} \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), \end{aligned}$$

  (C36)
  $$\begin{aligned} \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 -A r_{x}-B r_{x}}, \end{aligned}$$

  (C37)

  where the parameters $A, g_P,k_P$ and B have been defined as,

  $$\begin{aligned} 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}}{4A B}, \end{aligned}$$

  (C38)
  $$\begin{aligned} a_{1}&=\sqrt{-\frac{(r_{3}-r_{4})^{2}}{4}},\quad b_{1}=\frac{r_{3}+r_{4}}{2}, \end{aligned}$$

  (C39)

  we have also defined,

  $$\begin{aligned} 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{B r_{2}+A r_{1}-r_{x}(A+B)}{B r_{2}-r_{1}A+r_{x}(A+B)},\quad \alpha _{2}=\frac{A+B}{B-A}, \end{aligned}$$

  (C40)
  $$\begin{aligned} 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| }, \end{aligned}$$

  (C41)
  $$\begin{aligned} 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] . \end{aligned}$$

  (C42)

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.

1.2 2. Kerr de Sitter revisited

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

$$\begin{aligned} {\hat{I}}_{{\hat{\theta }}}&= \int _{{\hat{\theta }}_{s}}^{{\hat{\theta }}_{o}} \frac{\text {d}\theta }{\sqrt{{\hat{\Theta }}({\hat{\theta }})}}, \end{aligned}$$

(C43)

$$\begin{aligned} {\hat{I}}_r&=\int _{r_{s}}^{r_{o}} \frac{\text {dr}}{\sqrt{{\hat{R}}(r)}}, \end{aligned}$$

(C44)

$$\begin{aligned} {\hat{\phi }}_{o}-{\hat{\phi }}_{s}&= \int _{r_{s}}^{r_{o}} \left( \frac{a r^{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 }, \end{aligned}$$

(C45)

$$\begin{aligned} 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}, \end{aligned}$$

(C46)

for which we have defined,

$$\begin{aligned} {\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 )}}. \end{aligned}$$

(C47)

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

$$\begin{aligned} {\hat{\Theta }}_{{\hat{u}}}= \frac{a^{2}}{1-{\hat{u}}}({\hat{u}}_{+}-{\hat{u}})({\hat{u}}-{\hat{u}}_{-}). \end{aligned}$$

(C48)

The angular integrals still unpack as Eq. (D1). We obtain the general solutions as,

• 1.
  $$\begin{aligned} \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) , \end{aligned}$$

  (C49)

  where,

  $$\begin{aligned} {\hat{x}}_{i}= \arcsin \left( \sqrt{\frac{\cos ^{2}\theta _{i}}{{\hat{u}}_{+}}}\right) , {\hat{k}}=\frac{{\hat{u}}_{+}}{{\hat{u}}_{-}} . \end{aligned}$$

  (C50)
• 2.
  $$\begin{aligned} \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}}]. \end{aligned}$$

  (C51)
• 3.
  $$\begin{aligned} \left| \int _{\pi /2}^{\theta _{i}} \frac{\cos ^{2}\theta \text {d}\theta }{\sqrt{{\hat{\Theta }}(\theta )}}\right| = \dfrac{1}{\sqrt{-u_{-}a^{2}}}\int _{0}^{{\hat{x}}_{i}} \dfrac{{\hat{u}}_{+}\sin ^{2}x \text {dx}}{\sqrt{1-k\sin ^{2}x }} =\dfrac{1}{\sqrt{-u_{-}a^{2}}} \frac{{\hat{u}}_{+}(F({\hat{x}}_{i} |{\hat{k}}) -E({\hat{x}}_{i} |{\hat{k}}))}{{\hat{k}}}. \end{aligned}$$

  (C52)

Moreover, in order to obtain the time delay and azimuthal angle parameter, we employ a similar approach as in the case of KdS. This involves integrating over half orbits and treating the radial coordinate as a constant. Thus, we express Eqs. (C45) and (C46) as,

$$\begin{aligned} {\hat{\delta }}&= \left( \frac{a r^{2}+a^{3}-a {\hat{\Delta }}-a^{2}{\hat{\lambda }}}{{\hat{\Delta }}} \right) \left| \int _{\theta _1}^{\theta _4} \frac{\text {d}\theta }{\sqrt{{\hat{\Theta }}(\theta )}}\right| + {\hat{\lambda }} \left| \int _{\theta _1}^{\theta _{4}} \frac{\text {d}\theta }{\sin ^{2} \theta \sqrt{{\hat{\Theta }}(\theta )}}\right| , \end{aligned}$$

(C53)

$$\begin{aligned} {\hat{\zeta }}&= \left( \frac{(r^{2}+a^{2})(r^{2}+a^{2}-a \lambda )}{{\hat{\Delta }}} +a {\hat{\lambda }}-a^{2} \right) \left| \int _{\theta _1}^{\theta _4} \frac{\text {d}\theta }{\sqrt{{\hat{\Theta }}(\theta )}}\right| +a^{2} \left| \int _{\theta _1}^{\theta _{4}} \frac{\cos ^{2}\theta \text {d}\theta }{\sqrt{{\hat{\Theta }}(\theta )}}\right| , \end{aligned}$$

(C54)

where we have employed the equality relation given by Eq. (C6) to transform the Mino parameter into an angular integral. Making use of the general solutions Eqs. (C49–C52), we obtain,

$$\begin{aligned} \left| \int _{\theta _1}^{\theta _4} \frac{\text {d}\theta }{\sqrt{{\hat{\Theta }}(\theta )}}\right|&=\dfrac{2}{\sqrt{-{\hat{u}}_{-} a^{2}}} K\left( {\hat{k}}\right) ={\hat{\Gamma }}_\theta , \left| \int _{\theta _1}^{\theta _{4}} \frac{\cos ^{2}\theta \text {d}\theta }{\sqrt{{\hat{\Theta }}(\theta )}}\right| =\dfrac{2{\hat{u}}_{+}}{{\hat{k}}\sqrt{-u_{-}a^{2}}}\left( K({\hat{k}})-E({\hat{k}})\right) ={\hat{\Gamma }}_t, \end{aligned}$$

(C55)

$$\begin{aligned} \left| \int _{\theta _1}^{\theta _{4}} \frac{\text {d}\theta }{\sin ^{2} \theta \sqrt{{\hat{\Theta }}(\theta )}}\right|&=\dfrac{2}{\sqrt{-u_{-}a^{2}}} \Pi ({\hat{u}}_{+} |{\hat{k}} )={\hat{\Gamma }}_\phi . \end{aligned}$$

(C56)

Moreover, the Lyapunov exponent per half orbits is derived using the same approach in KdS through the relation,

$$\begin{aligned} {\hat{\gamma }}=\sqrt{\frac{\hat{R''}(r_b)}{2}} {\hat{\Gamma }}_\theta . \end{aligned}$$

(C57)

To obtain solutions for the radial integrals, firstly we express the radial potential and the parameter ${\hat{\delta }}_r$ in terms of their roots,

$$\begin{aligned} {\hat{R}}(r)=(r-{\hat{r}}_{1})(r-{\hat{r}}_{2})(r-{\hat{r}}_{3})(r-{\hat{r}}_{4}), {\hat{\Delta }}=-\frac{3}{\Lambda }(r-{\hat{r}}_{-c})(r-{\hat{r}}_{-})(r-{\hat{r}}_{+})(r-{\hat{r}}_{c}). \end{aligned}$$

(C58)

Using Eq. (C58) we rewrite Eqs. (C44–C46) as,

$$\begin{aligned} {\hat{I}}_{r}&=\int _{r_{s}}^{r_{o}} \frac{\text {dr}}{\pm \sqrt{(r-{\hat{r}}_{1})(r-{\hat{r}}_{2})(r-{\hat{r}}_{3})(r-{\hat{r}}_{4})}}, \end{aligned}$$

(C59)

$$\begin{aligned} {\hat{\phi }}_{o}-{\hat{\phi }}_{s}&=a \left( \frac{{\hat{\Omega }}_{-c}}{r_{c,-c}r_{-,-c}r_{+,-c}}{\hat{I}}_{-c}+\frac{{\hat{\Omega }}_{-}}{r_{c,-}r_{-c,-}r_{+,-}}{\hat{I}}_{-}+\frac{{\hat{\Omega }}_{+}}{r_{c,+}r_{-c,+}r_{-,+}}{\hat{I}}_{+}+\frac{{\hat{\Omega }}_{c}}{r_{-c,c}r_{-,c}r_{+,c}}{\hat{I}}_{c}-\tau \right) +{\hat{\lambda }} {\hat{I}}_{\phi }, \end{aligned}$$

(C60)

$$\begin{aligned} {\hat{t}}_{o}-{\hat{t}}_{s}&= \frac{{\hat{\Phi }}_{c}}{{\hat{r}}_{-c,c} {\hat{r}}_{-,c} {\hat{r}}_{+,c}}+\frac{{\hat{\Phi }}_{-c}}{{\hat{r}}_{c,-c} {\hat{r}}_{c,-c} {\hat{r}}_{+,-c}}+\frac{{\hat{\Phi }}_{-}}{{\hat{r}}_{c,-} {\hat{r}}_{-c,-} {\hat{r}}_{+,-}}+\frac{{\hat{\Phi }}_{+}}{{\hat{r}}_{c,+} {\hat{r}}_{-c,+} {\hat{r}}_{-,+}}+\tau \left( -a^2+a {\hat{\lambda }} -\frac{\Lambda }{3}\right) , \end{aligned}$$

(C61)

$$\begin{aligned} {\hat{\Phi }}_{x}&=\left( a^2+{\hat{r}}_{x}^2\right) {\hat{\Omega }}_{x} {\hat{I}}_{x}, {\hat{\Omega }}_{x}=\frac{3 }{\Lambda }\left( a^2-a {\hat{\lambda }} +{\hat{r}}_{x}^2\right) ,{\hat{I}}_{x}=\int _{r_{s}}^{r_{o}} \frac{\text {dr}}{(r-{\hat{r}}_{x}) \sqrt{{\hat{R}}(r)}}. \end{aligned}$$

(C62)

The subscripts x denotes the horizons $-c,-,+,c$, where ${\hat{r}}_{+}$, ${\hat{r}}_{-}$, ${\hat{r}}_{c}$, and ${\hat{r}}_{-c}$ represent the event horizon, Cauchy horizon, the cosmological horizon and its dual, respectively. Further, ${\hat{r}}_{a,b}={\hat{r}}_{a}-r_{b}$. Further, the subscript o and s denote the source and the observer respectively, where $r_o$ and $r_s$ need to be chosen within the domain of outer communication. The radial integrals ${\hat{I}}_r$ and ${\hat{I}}_x$ takes the same form as that of KdS. In the same way we obtain their solutions for the roots (${\hat{r}}_{1}<{\hat{r}}_{2}<{\hat{r}}_{3}<{\hat{r}}_{4}<r_s$) and (${\hat{r}}_{1}<{\hat{r}}_{2}<r_{s}, {\hat{r}}_{3}=\bar{{\hat{r}}}_{4}$).

• $\underline{{\hat{r}}_{1}<{\hat{r}}_{2}<{\hat{r}}_{3}<{\hat{r}}_{4}<r_s }$

  $$\begin{aligned}&\int _{r_{4}}^{r_{i}} \frac{\text {dr}}{\pm \sqrt{(r-{\hat{r}}_{1})(r-{\hat{r}}_{2})(r-{\hat{r}}_{3})(r-{\hat{r}}_{4})}}={\hat{g}}_{E}F({\hat{\varphi }}_{E,i}|{\hat{k}}_{E}), \end{aligned}$$

  (C63)
  $$\begin{aligned}&\int _{r_{4}}^{r_{i}} \frac{\text {dr}}{{\pm _r}(r-{\hat{r}}_x)\sqrt{(r-{\hat{r}}_{1})(r-{\hat{r}}_{2})(r-{\hat{r}}_{3})(r-{\hat{r}}_{4})}}=-\frac{2}{{\hat{r}}_{x,3}. \sqrt{{\hat{r}}_{4,2} {\hat{r}}_{3,1}}} \left( F\left( {\hat{\varphi }}_{E,i}|{\hat{k}}_{E}\right) +\frac{{\hat{r}}_{4,3}{\hat{\Pi }}_{i})}{{\hat{r}}_{x,4}} \right) , \end{aligned}$$

  (C64)

  where we have defined,

  $$\begin{aligned} {\hat{\varphi }}_{E,i}=\arcsin \sqrt{\frac{{\hat{r}}_{3,1}{\hat{r}}_{i,4}}{{\hat{r}}_{4,1}{\hat{r}}_{i,3}}}, {\hat{k}}_{E}=\frac{{\hat{r}}_{4,1} {\hat{r}}_{3,2}}{{\hat{r}}_{4,2} {\hat{r}}_{3,1}} , {\hat{g}}_{E}=\frac{2}{\sqrt{{\hat{r}}_{4,2}{\hat{r}}_{3,1}}}, {\hat{\Pi }}_{i}=\Pi \left( \frac{{\hat{r}}_{4,1} {\hat{r}}_{x,3}}{{\hat{r}}_{x,4} {\hat{r}}_{3,1}};{\hat{\varphi }}_{E,i}|{\hat{k}}_{E}\right) \end{aligned}$$

  (C65)

• $\underline{{\hat{r}}_{1}<{\hat{r}}_{2}<r_{s}, {\hat{r}}_{3}=\bar{{\hat{r}}}_{4}}$

  In the same way as Eq. (C36) and (C37) we obtain the solution for this case in RKdS as,

  $$\begin{aligned}&\int _{r_{2}}^{r_{i}} \frac{\text {dr}}{\pm \sqrt{(r-{\hat{r}}_{1})(r-{\hat{r}}_{2})(r-{\hat{r}}_{3})(r-{\hat{r}}_{4})}}={\hat{g}}_P F({\hat{\varphi }}_{P,i}\mid k_P), \end{aligned}$$

  (C66)
  $$\begin{aligned}&\int _{r_{2}}^{r_{i}} \frac{\text {dr}}{{\pm _r}(r-{\hat{r}}_x)\sqrt{(r-{\hat{r}}_{1})(r-{\hat{r}}_{2})(r-{\hat{r}}_{3})(r-{\hat{r}}_{4})}}= \frac{({\hat{B}}-{\hat{A}}) \left( ({\hat{\alpha }} -{\hat{\alpha }}_{2}) {\hat{R}}_{1,i}+{\hat{\alpha }}_{2}F({\hat{\varphi }}_{P,i} \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})}. \end{aligned}$$

  (C67)

  The parameters ${\hat{A}},{\hat{B}},{\hat{k}}_P$ have been defined as,

  $$\begin{aligned}&{\hat{A}}=\sqrt{\left( {\hat{r}}_{2}-\frac{{\hat{r}}_{3}+{\hat{r}}_{4}}{2}\right) ^{2}-\frac{({\hat{r}}_{3}-{\hat{r}}_{4})^{2}}{4}}, {\hat{B}}=\sqrt{\left( {\hat{r}}_{1}-\frac{{\hat{r}}_{3}+{\hat{r}}_{4}}{2} \right) ^{2}-\frac{({\hat{r}}_{3}-{\hat{r}}_{4})^{2}}{4}}, \end{aligned}$$

  (C68)
  $$\begin{aligned}&{\hat{k}}_{P}=\frac{({\hat{A}}+{\hat{B}})^{2}-({\hat{r}}_{2}-{\hat{r}}_{1})^{2}}{4{\hat{A}} {\hat{B}}}. \end{aligned}$$

  (C69)

  We have also defined,

  $$\begin{aligned}&{\hat{R}}_{1,i}=\frac{1}{1-{\hat{\alpha }}^2}\left( \Pi \left[ \frac{{\hat{\alpha }}^2}{{\hat{\alpha }}^2-1};{\hat{\varphi }}_{P,i} \mid {\hat{k}}_{P} \right] -{\hat{\alpha }} {\hat{f}}_{1,i}\right) ,\quad {\hat{\alpha }}=\frac{{\hat{B}} r_{2}+{\hat{A}}r_{1}-r_{x}(A+{\hat{B}})}{{\hat{B}} r_{2}-r_{1}{\hat{A}}+r_{x}({\hat{A}}+{\hat{B}})},{\hat{\alpha }}_{2}=\frac{{\hat{A}}+{\hat{B}}}{{\hat{B}}-{\hat{A}}}, \end{aligned}$$

  (C70)
  $$\begin{aligned}&{\hat{f}}_{1,i}=\sqrt{\frac{{\hat{\alpha }}^2-1}{{\hat{k}}_{P}+(1-{\hat{k}}_{P}){\hat{\alpha }}^2}}\ln {\left| \frac{\sqrt{{\hat{k}}_{P}+(1-{\hat{k}}_{P}){\hat{\alpha }}^2}dn({\hat{u}}_{i}\mid {\hat{k}}_{P})+\sqrt{{\hat{\alpha }}^2-1}sn({\hat{u}}_{i}\mid {\hat{k}}_{P})}{\sqrt{{\hat{k}}_{P}+(1-{\hat{k}}_{P}){\hat{\alpha }}^2}dn({\hat{u}}_{i}\mid {\hat{k}}_{P})-\sqrt{{\hat{\alpha }}^2-1}sn({\hat{u}}_{i}\mid {\hat{k}}_{P})} \right| }, \end{aligned}$$

  (C71)
  $$\begin{aligned}&{\hat{u}}_{i}=cn^{-1}\left[ \cos ({\hat{\varphi }}_{P,i}) \mid {\hat{k}}_{P} \right] ,{\hat{\varphi }}_{P,i}=\arccos \left[ \frac{({\hat{A}}-{\hat{B}})r_{i}+{\hat{r}}_{2} {\hat{B}}-{\hat{r}}_{1} A}{({\hat{A}}+{\hat{B}})r_{i}-{\hat{r}}_{2} {\hat{B}}-{\hat{r}}_{1}{\hat{A}}} \right] . \end{aligned}$$

  (C72)

Appendix D: Integrals in terms of angular turning points

For purposes of analytic ray tracing, it becomes imperative to express the solutions in terms of the turning points experienced by photons. In this section we will express the angular integrals in terms of the number of turning points that the photons encounter in the $\theta$ direction. Considering all possible configurations, Ref. [31] demonstrates that angular path integrals for ordinary motion unpack as,

$$\begin{aligned} \int _{\theta _s}^{\theta _o}=2m\left| \int _{\pi /2}^{\theta _\pm } \right| +\eta _s \left| \int _{\pi /2}^{\theta _s}\right| -\eta _o \left| \int _{\pi /2}^{\theta _o}\right| . \end{aligned}$$

(D1)

Where,

$$\begin{aligned} \eta _s=sign(p^{\theta }_{s})sign(\cos \theta _{s})=(-1)^m sign(p^{\theta }_{o})sign(\cos \theta _s), \eta _o=sign(p^{\theta }_{o})sign(\cos \theta _{o})=(-1)^m sign(p^{\theta }_{s})sign(\cos \theta _o), \end{aligned}$$

(D2)

and m is the number of turning points in the latitudinal direction. We have already obtained all the general solutions thus we will directly apply our solutions to Eq. (D1).

In KdS spacetime we obtain,

$$\begin{aligned}&{\mathbb {I}}_\theta = \dfrac{2 m}{\sqrt{-u_{-}\Xi }} K\left( k\right) + \dfrac{ (-1)^m sign(p^{\theta }_{o})}{\sqrt{-u_{-}\Xi }} F\left( x_{s}| k\right) -\dfrac{sign(p^{\theta }_{o})}{\sqrt{-u_{-}\Xi }} F\left( x_{o}| k\right) , \end{aligned}$$

(D3)

$$\begin{aligned}&{\mathbb {I}}_\phi =\dfrac{2 m}{\sqrt{-u_{-}\Xi }}\Pi \left( -u_{+} {\mathcal {Y}};|k\right) +\dfrac{(-1)^m sign(p^{\theta }_{o})}{\sqrt{-u_{-}\Xi }}\Pi \left( -u_{+} {\mathcal {Y}};x_s|k\right) -\dfrac{ sign(p^{\theta }_{o})}{\sqrt{-u_{-}\Xi }}\Pi \left( -u_{+} {\mathcal {Y}};x_o|k\right) , \end{aligned}$$

(D4)

$$\begin{aligned}&{\mathbb {I}}_t= \frac{2 m u_{+} J(-u_+ {\mathcal {Y}}|k)}{\sqrt{-u_{-}\Xi }}+\frac{(-1)^m sign(p^{\theta }_{o}) u_{+} J(-u_+ {\mathcal {Y}},\arcsin t_s|k)}{\sqrt{-u_{-}\Xi }}-\frac{ sign(p^{\theta }_{o}) u_{+} J(-u_+ {\mathcal {Y}},\arcsin t_o|k)}{\sqrt{-u_{-}\Xi }}, \end{aligned}$$

(D5)

$$\begin{aligned}&\bar{{\mathbb {I}}}_\phi =\dfrac{2 m}{\sqrt{-u_{-}\Xi }} \dfrac{1}{1+{\mathcal {Y}}}[\Pi (u_{+}|k)+{\mathcal {Y}} \Pi (-u_{+}{\mathcal {Y}}|k) ]+\dfrac{(-1)^m sign(p^{\theta }_{o})}{\sqrt{-u_{-}\Xi }} \dfrac{1}{1+{\mathcal {Y}}}[\Pi (u_{+};\arcsin t_s|k)+{\mathcal {Y}} \Pi (-u_{+}{\mathcal {Y}};\arcsin t_s|k) ] \nonumber \\&\quad -\dfrac{ sign(p^{\theta }_{o})}{\sqrt{-u_{-}\Xi }} \dfrac{1}{1+{\mathcal {Y}}}[\Pi (u_{+};\arcsin t_o|k)+{\mathcal {Y}} \Pi (-u_{+}{\mathcal {Y}};\arcsin t_o|k) ], \end{aligned}$$

(D6)

while in RKdS spacetime we obtain,

$$\begin{aligned}&\hat{{\mathbb {I}}}_\theta =\dfrac{2 m}{\sqrt{-{\hat{u}}_{-}a^2}} K\left( {\hat{k}}\right) + \dfrac{ (-1)^m sign(p^{\theta }_{o})}{\sqrt{-{\hat{u}}_{-}a^2}} F\left( {\hat{x}}_{s}| {\hat{k}}\right) -\dfrac{sign(p^{\theta }_{o})}{\sqrt{-{\hat{u}}_{-}a^2}} F\left( {\hat{x}}_{o}| {\hat{k}}\right) , \end{aligned}$$

(D7)

$$\begin{aligned}&\hat{{\mathbb {I}}}_\phi =\dfrac{2m}{\sqrt{-{\hat{u}}a^{2}}} \Pi ({\hat{u}}_{+} |{\hat{k}} )+\dfrac{(-1)^m sign(p^{\theta }_{o})}{\sqrt{-{\hat{u}}a^{2}}} \Pi [{\hat{u}}_{+};{\hat{x}}_{s} |{\hat{k}}]-\dfrac{sign(p^{\theta }_{o})}{\sqrt{-{\hat{u}}a^{2}}} \Pi [{\hat{u}}_{+};{\hat{x}}_{o} |{\hat{k}}], \end{aligned}$$

(D8)

$$\begin{aligned}&\hat{{\mathbb {I}}}_t=\dfrac{2m}{\sqrt{-{\hat{u}}a^{2}}} \frac{{\hat{u}}_{+}(F({\hat{k}}) -E({\hat{k}}))}{{\hat{k}}}+\dfrac{(-1)^m sign(p^{\theta }_{o})}{\sqrt{-{\hat{u}}a^{2}}} \frac{{\hat{u}}_{+}(F({\hat{x}}_{s} |{\hat{k}}) -E({\hat{x}}_{s} |{\hat{k}}))}{{\hat{k}}}\nonumber \\&\qquad -\dfrac{ sign(p^{\theta }_{o})}{\sqrt{-{\hat{u}}a^{2}}} \frac{{\hat{u}}_{+}(F({\hat{x}}_{o} |{\hat{k}}) -E({\hat{x}}_{o} |{\hat{k}}))}{{\hat{k}}}. \end{aligned}$$

(D9)

The $sign(p^{\theta }_{o})$ is equivalent to sign(y). The above equations are only valid for observers at latitude $\theta _o \in (0,\pi /2)$.

We postpone the case of $\theta _o=0,\pi /2$ to our future work.

Appendix E: Change in $\phi$ and t for subsequent images

In sect. 3, an analysis of direct images, lensing rings and photon rings has been done. In the same section we have slightly talked about subsequent images, ($(m+1)^{th}$,$m^{th}$), undergoing a rotation due to change in azimuthal angle and a delay in time of detection. In this part we explicitly derive these quantities. Giving the explicit forms of these quantities is important as we shall see that they are related to the critical parameters controlling the photon ring structure.

To obtain the KdS change in azimuthal angle between subsequent images $\Delta \phi _{m+1}-\Delta \phi _{m}$ we proceed as follows,

$$\begin{aligned} \Delta \phi _{m+1}-\Delta \phi _{m}&= \left( \int _{r_{o}}^{r_{s}} \frac{\left( a L^2\right) \left( a (a-\lambda )+\text {r}^2\right) \text {dr}}{\pm _{r} {\Delta _{r} } \sqrt{R(r)}} -L^{2}a {\mathbb {I}}_{\phi }+ L^{2}\lambda \bar{{\mathbb {I}}}_{\phi }\right) _{m+1}\nonumber \\&\qquad - \left( \int _{r_{s}}^{r_{o}} \frac{\left( a L^2\right) \left( a (a-\lambda )+\text {r}^2\right) \text {dr}}{\pm _{r} {\Delta _{r} } \sqrt{R(r)}} -L^{2}a {\mathbb {I}}_{\phi }+ L^{2}\lambda \bar{{\mathbb {I}}}_{\phi }\right) _{m}, \end{aligned}$$

(E1a)

$$\begin{aligned}&= \int _{r_{s}^{m+1}}^{r_{s}^{m}} \frac{\left( a L^2\right) \left( a (a-\lambda )+\text {r}^2\right) \text {dr}}{\pm _{r} {\Delta _{r} } \sqrt{R(r)}} -\dfrac{2 L^{2}a}{\sqrt{-u_{-}\Xi }}\Pi \left( -u_{+} {\mathcal {Y}};|k\right) + \dfrac{2 L^{2}\lambda }{(1+{\mathcal {Y}})\sqrt{-u_{-}\Xi }} [\Pi (u_{+}|k)+{\mathcal {Y}} \Pi (-u_{+}{\mathcal {Y}}|k) ], \end{aligned}$$

(E1b)

$$\begin{aligned}&= \frac{\left( a L^2\right) \left( a (a-\lambda )+\text {r}^2\right) }{{\Delta _{r} }} \Gamma _\theta -\dfrac{2 L^{2}a}{\sqrt{-u_{-}\Xi }}\Pi \left( -u_{+} {\mathcal {Y}};|k\right) + \dfrac{2 L^{2}\lambda }{(1+{\mathcal {Y}})\sqrt{-u_{-}\Xi }} [\Pi (u_{+}|k)+{\mathcal {Y}} \Pi (-u_{+}{\mathcal {Y}}|k) ]. \end{aligned}$$

(E1c)

In Eq. (E1b) we have performed a subtraction alongside inserting Eqs. (D4) and (D6). Subsequent images are just images of same sources that have undergone extreme lensing. Therefore, for each image, $r_{s}^{m+1}$ is simply the same as $r_{s}^{m}$. Thus on the radial integral, we have kept r to be constant. Additionally, $r_{s}^{m}$ to $r_{s}^{m+1}$ is equivalent to half an orbit. The Mino period for half an orbit is as given in Eq. (C18c). This is the approach we have used to arrive at Eq. (E1c).

Performing the same analysis for $\Delta t$ gives us,

$$\begin{aligned} \Delta t_{m+1}-\Delta t_{m}&= \dfrac{L^{2}((r^{2}+a^{2})^{2}-a\lambda (a^{2}+r^{2}))}{ \Delta _{r}} \Gamma _\theta +\frac{2 a^2 L^2 u_{+} J(-u_+ {\mathcal {Y}}|k)}{\sqrt{-u_{-}\Xi }}-\dfrac{2 a L^2(a-\lambda )}{\sqrt{-u_{-}\Xi }}\Pi \left( -u_{+} {\mathcal {Y}};|k\right) . \end{aligned}$$

(E2)

Likewise in RKdS spacetime we obtain these parameters as,

$$\begin{aligned} {\hat{\Delta }}\phi _{m+1}-{\hat{\Delta }}\phi _{m}&= \left( \frac{a r^{2}+a^{3}-a {\hat{\Delta }}-a^{2}{\hat{\lambda }}}{{\hat{\Delta }}} \right) {\hat{\Gamma }}_\theta +\dfrac{2 {\hat{\lambda }}}{\sqrt{-u_{-}a^{2}}} \Pi ({\hat{u}}_{+} |{\hat{k}} ), \end{aligned}$$

(E3)

$$\begin{aligned} {\hat{\Delta }}t_{m+1}-{\hat{\Delta }} t_{m}&= \left( \frac{(r^{2}+a^{2})(r^{2}+a^{2}-a \lambda )}{{\hat{\Delta }}} +a {\hat{\lambda }}-a^{2} \right) {\hat{\Gamma }}_\theta + \dfrac{2 a^{2}}{\sqrt{-u_{-}a^{2}}} \frac{{\hat{u}}_{+}(E({\hat{k}}) -E({\hat{k}}))}{{\hat{k}}}. \end{aligned}$$

(E4)