# Light propagation around a Kerr-like black hole immersed in an inhomogeneous
anisotropic plasma in Rastall gravity I: Analytical solutions to the equations
of motion

Mohsen Fathi<EMAIL_ADDRESS>Departamento de Física, Universidad de
Santiago de Chile, Avenida Víctor Jara 3493, Estación Central, 9170124,
Santiago, Chile  Marco Olivares<EMAIL_ADDRESS>Facultad de
Ingeniería y Ciencias, Universidad Diego Portales, Avenida Ejército Libertador
441, Santiago, Chile  J.R. Villanueva<EMAIL_ADDRESS>Instituto de
Física y Astronomía, Universidad de Valparaíso, Avenida Gran Bretaña 1111,
Valparaíso, Chile  Norman Cruz<EMAIL_ADDRESS>Departamento de Física,
Universidad de Santiago de Chile, Avenida Víctor Jara 3493, Estación Central,
9170124, Santiago, Chile

###### Abstract

In this paper, we explore the behavior of light ray trajectories in the
exterior geometry of a rotating black hole within the Rastall theory of
gravity, which is surrounded by an inhomogeneous anisotropic electronic cold
plasma. By specifying the plasma’s frequency profile, we derive fully
analytical solutions for the temporal evolution of spacetime coordinates using
elliptic integrals and Jacobi elliptic functions. These solutions illustrate
various possible orbits. Throughout the study, we compare results with the
vacuum case, emphasizing the influence of the plasma. Additionally, we utilize
the analytical solutions to establish the lens equation for the considered
spacetime. The investigation also addresses the significance of spherical
photon orbits in critical trajectories, by presenting several examples.

keywords: Photon motion, Rastall gravity, plasmic medium, gravitational
lensing

PACS numbers: 04.20.Fy, 04.20.Jb, 04.25.-g

##  I Introduction and Motivation

While the general theory of relativity remains highly regarded in astrophysics
and cosmology for its elegance and beauty, recent observations like Type Ia
supernovae, cosmic microwave background radiation anisotropies, flat galactic
rotation curves, and weak lensing have shifted scientific perspectives. These
observations suggest that only a small portion of the observable universe
consists of ordinary baryonic matter, with the remainder attributed to dark
matter and dark energy riess_type_2004 ; komatsu_seven-year_2011 . To date,
scientists have devised numerous theoretical models to explain the phenomena
of dark matter and dark energy. (for example, see Refs. turner_coherent_1983 ;
press_single_1990 ; dubinski_structure_1991 ; sin_late-time_1994 ;
navarro_structure_1996 ; navarro_universal_1997 ; hu_fuzzy_2000 ;
peebles_fluid_2000 ; peebles_cosmological_2003 ; kiselev2003quintessential ;
viel_constraining_2005 ; amendola_dark_2006 ; copeland_dynamics_2006 ;
de_la_macorra_dark_2006 ; rahaman_perfect_2010 ; rahaman_modeling_2011 ;
li_dark_2011 ; schive_cosmic_2014 ; marsh_axion_2016 ; hui_ultralight_2017 ;
wang_holographic_2017 ). Nevertheless, the enigmatic nature of dark matter and
dark energy persists as unresolved quandaries in modern cosmology.
Consequently, some scientists speculate that by adjusting the general theory
of relativity and reshaping our understanding of gravitational behavior, it
might be plausible to elucidate the anomalous phenomena currently attributed
to the dark side of the universe. Among the various modified gravity models
proposed, the Rastall theory of gravity stands as one such alternative
Rastall19723357 ; rastall_theory_1976 , which follows a broader energy-
momentum conservation law, directly connected to the Mach principle. Since its
introduction, this theory has found various applications in both gravitation
and cosmology (for example, see Refs. batista_rastall_2012 ;
fabris_rastall_2012 ; fabris_rastalls_2015 ; heydarzade_black_2017-1 ;
heydarzade_black_2017 ; darabi_einstein_2018 ; lobo_thermodynamics_2018 ).
Additionally, the Rastall theory produces interesting outcomes, particularly
in its black hole solutions that incorporate additional terms to accommodate
the universe’s dark components. The incorporation of dark matter/energy
effects into black hole spacetimes has been a longstanding focus of scientific
inquiry. Kottler, for instance, investigated the impact of the cosmological
constant on the Schwarzschild black hole spacetime kottler_uber_1918 , later
extended to its Kerr-like solution by Carter Carter:1973rla . Kiselev was
among the first to reconcile quintessential dark energy with the Schwarzschild
black hole spacetime kiselev_quintessence_2003 , a concept later extended to
Kerr and Kerr-Newman-anti de Sitter black holes toshmatov_rotating_2017 . A
similar solution has been proposed for Rastall gravity in Ref.
heydarzade_black_2017 , involving a Reissner-Nordström black hole spacetime
embedded within perfect fluid matter. This solution was further generalized in
Refs. kumar_rotating_2018 ; xu_kerrnewman-ads_2018  to describe a Kerr-Newman-
like black hole within Rastall gravity, examining certain characteristics of
the dark fluid within this spacetime context. However, in this study, we focus
on a similar black hole with zero net charge. This leads to a Kerr-like black
hole within a cosmological fluid background, allowing us to investigate the
behavior of light rays within its exterior spacetime in this work.

In fact, the study of light behavior has held significant importance in
relativistic astrophysics. From Eddington’s groundbreaking confirmation of
light deflection during the 1919 solar eclipse Eddington:1920 , validating
Einstein’s predictions Einstein:1911 , to recent revelations from the Event
Horizon Telescope (EHT) capturing the shadows of M87* Akiyama:2019  and Sgr A*
Akiyama:2022 , the role of light as our primary source of information from
celestial objects has become increasingly evident. Black holes, in particular,
remain enigmatic as they reside beyond an event horizon where light becomes
trapped, rendering their interiors inscrutable. Thus, our only means of
studying light behavior around black holes is beyond this event horizon. This
study focuses on the limits to which light can approach a Kerr-like Rastall
black hole, employing fully analytical methods to analyze photon behavior
within these specific regions.

Indeed, the study of geodesics in black hole spacetimes began soon after the
advent of general relativity. However, it was Hagihara’s 1930 derivation of
particle trajectories in Schwarzschild spacetime that marked a reliable
analytical breakthrough. Efforts since then have consistently pursued
gravitational lensing caused by black holes, yielding significant
formulations, including methods to constrain the black hole shadow
Bardeen:1972a ; Bardeen:1973a ; Bardeen:1973b ; Chandrasekhar:2002 .
Analytical treatment of particle geodesics aids in strict theoretical
predictions and their comparison with observational data, attracting interest
among scientists. Such analysis also helped revive the study of gravitational
lensing Virbhadra:2000 . However, non-static black hole spacetimes involve
complex equations requiring advanced mathematical methods, leading to numerous
investigations. In the realm of elliptic integrals’ application in calculating
light ray geodesics (relevant to this work), various criteria have been
discussed and resolved in Refs. kraniotis_precise_2004 ; beckwith_extreme_2005
; kraniotis_frame_2005 ; hackmann_complete_2008 ; hackmann_geodesic_2008 ;
bisnovatyi-kogan_strong_2008 ; kagramanova_analytic_2010 ;
hackmann_analytical_2010 ; hackmann_complete_2010 ; enolski_inversion_2011 ;
kraniotis_precise_2011 ; enolski_inversion_2012 ; gibbons_application_2012 ;
munoz_orbits_2014 ; kraniotis_gravitational_2014 ; de_falco_approximate_2016 ;
soroushfar_detailed_2016 ; barlow_asymptotically_2017 ; uniyal_null_2018 ;
villanueva_null_2018 ; chatterjee_analytic_2019 ; hsiao_equatorial_2020 ;
gralla_null_2020 ; hendi_simulation_2020 ; fathi_classical_2020 ;
kraniotis_gravitational_2021 ; fathi_study_2022 ; fathi_spherical_2023 . Most
recently, these methods have been employed to conduct analytical ray-tracing
of the shadow of rotating black holes in cosmological backgrounds, as
demonstrated in Refs. Omwoyo:2023 ; Wang:2023 . Furthermore, Synge’s
discussions in Ref. Synge:1960  enabled mathematical methods for determining
and constraining light propagation in theoretical (non-)static black hole
spacetimes within plasmic or dark fluid environments bisnovatyi-
kogan_gravitational_2009 ; bisnovatyi-kogan_gravitational_2010 ;
tsupko_gravitational_2013 ; morozova_gravitational_2013 ; bisnovatyi-
kogan_gravitational_2015 ; perlick_influence_2015 ; atamurotov_optical_2015 ;
abdujabbarov_shadow_2016 ; bisnovatyi-kogan_gravitational_2017 ;
perlick_light_2017 ; schulze-koops_sachs_2017 ; abdujabbarov_shadow_2017 ;
liu_effects_2017 ; haroon_shadow_2019 ; kimpson_spatial_2019 ;
babar_optical_2020 ; junior_spinning_2020 ; Badia:2021 ; fathi_analytical_2021
. However, the analytical treatment of light ray trajectories in non-vacuum
black hole environments demands greater attention, offering the potential to
identify additional orbits often overlooked in determining the photon spheres
and shadows of black holes. Hence, building upon the introductory context
outlined earlier, we utilize these mathematical methods to compute light ray
paths within a non-vacuum, non-static spacetime governed by Rastall gravity.
Our study delves into exploring the optical characteristics inherent in the
outer spacetime of a Kerr-like black hole within the framework of Rastall
gravity. Accordingly, we focus on the possible photon orbits around a Kerr-
like black hole immersed in an inhomogeneous and anisotropic plasmic medium
whose structural properties remain constant throughout.

To obviate this objective, we have organize this study as follows: In Sect. II
we provide a brief introduction to the Rastall theory of gravity, outlining
its static and Kerr-like black hole solutions, and emphasizing the parameter
$\zeta$ in the cosmic fluid. We explore three distinct cases for this
parameter, discussing the resulting causal structure of the black hole. We
argue that $\zeta=3$ serves as the foundational choice, forming the basis for
subsequent sections. In Sect. III, we delve into the Hamiltonian dynamics of
light propagation within a refractive medium, particularly focusing on a non-
magnetized plasma. We derive and present the general equations of motion
governing light propagation within a stationary spacetime filled with such
plasma. Section IV constitutes the core of our study, analyzing the evolution
of spacetime coordinates, beginning with polar motion. We classify feasible
photon orbits using the radial effective potential and provide analytical
solutions to the equations of motion for all cases. We express these
solutions, when necessary, in terms of elliptic integrals and Jacobian
elliptic functions. Additionally, this section addresses spherical photon
orbits and their variations concerning the photons’ impact parameter. Further
insights into spherical orbits are discussed in Sect. V, followed by a
comprehensive summary of our findings in Sect. VI. Throughout this work, we
adopt a geometrized unit system ($G=c=1$) and adhere to the sign convention
$(-,+,+,+)$.

##  II Rastall gravity and its black hole solution

The field equations of the Rastall theory of gravity is given by
Rastall19723357

$R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R=\kappa\left(T_{\mu\nu}-\lambda
g_{\mu\nu}R\right),$ (1)

in which $R_{\mu\nu}$, $R$, and $T_{\mu\nu}$ represent the Ricci tensor, Ricci
scalar, and energy-momentum tensor, respectively. The gravitational coupling
constant in Rastall gravity, denoted as $\kappa$, contributes to our
discussions via a new parameter, $\psi=\kappa\lambda$, recognized as the
Rastall coupling parameter KUMAR2021100881 . This parameter holds significant
importance in distinguishing between non-vacuum general relativistic solutions
and those derived in Rastall gravity. Rastall gravity is founded on the
Rastall hypothesis, which, instead of the conventional energy-momentum
conservation law, ${T^{\mu\nu}}_{;\mu}=0$, adheres to the modified
conservation law ${T^{\nu}}_{\mu;\nu}=\lambda R_{,\mu}$. Notably, this theory
reinstates the common conservation law in flat spacetime. The parameter
$\lambda$ delineates the degree to which energy-momentum conservation is
upheld. As a result, the theory aligns with the Mach principle, postulating
that the inertia of a local mass correlates with the total mass/energy
distribution in the universe majernik2006rastalls .

In fact, the quest for exact solutions in Rastall gravity remains a relatively
recent trend. However, certain aspects of spherically symmetric black hole
solutions have surfaced in recent discussions. Notably, Ref.
heydarzade_black_2017-1  explored aspects of spherically symmetric black hole
solutions within the context of Rastall gravity, introducing a solution
grounded in an anisotropic fluid background, initially reported in Ref.
heydarzade_black_2017 . This solution is formulated based on the line element

$\mathrm{d}s^{2}=-f(r)\mathrm{d}t^{2}+\frac{\mathrm{d}r^{2}}{f(r)}+r^{2}\left(\mathrm{d}\theta^{2}+\sin^{2}\theta\mathrm{d}\phi^{2}\right),$
(2)

in the usual Schwarzschild coordinates, where the lapse function is given by

$f(r)=1-\frac{2M}{r}-N_{s}r^{-\frac{1+3w_{s}-6\psi(1+w_{s})}{1-3\psi(1+w_{s})}},$
(3)

where $N_{s}$ represents the matter density for a perfect fluid characterized
by the state parameter $w_{s}$, enveloping a static black hole with the mass
$M$. Note that, the lapse function (3) is valid for
$\psi\neq\frac{1}{3}(1+w_{s})^{-1}$. This fluid is characterized by an
anisotropic energy-momentum tensor, as described in Ref. Visser_2020 . In
Refs. kumar_rotating_2018 ; xu_kerrnewman-ads_2018 , the above static
spacetime was extended to a Kerr-Newman-like spacetime, by means of a modified
Newman-Janis algorithm, proposed by Azreg-AÏnou Azreg:2014 ; azreg-
ainou_static_2014 . In the Boyer–Lindquist coordinates, this stationary black
hole spacetime is characterized by the line element xu_kerrnewman-ads_2018

$\mathrm{d}s^{2}=-\left(\frac{\Delta-a^{2}\sin^{2}\theta}{\rho^{2}}\right)\mathrm{d}t^{2}+\frac{\rho^{2}}{\Delta}\mathrm{d}r^{2}+\rho^{2}\mathrm{d}\theta^{2}-2a\sin^{2}\theta\left(\frac{r^{2}+a^{2}-\Delta}{\rho^{2}}\right)\mathrm{d}t\mathrm{d}\phi\\\
+\sin^{2}\theta\left[\frac{(r^{2}+a^{2})^{2}-\Delta
a^{2}\sin^{2}\theta}{\rho^{2}}\right]\mathrm{d}\phi^{2},$ (4)

where for the case of a black hole with zero electric charge, we have

$\displaystyle\Delta(r)=r^{2}+a^{2}-2Mr-N_{s}r^{\zeta},$ (5a)
$\displaystyle\rho^{2}(r,\theta)={r^{2}+a^{2}\cos^{2}\theta},$ (5b)

with

$\zeta={1-3w_{s}\over 1-3\psi(1+w_{s})}.$ (6)

This way, the line element (4) characterizes the exterior spacetime of a Kerr-
like Rastall black hole (KRBH) with mass $M$ and spin parameter $a$, existing
within a fluid defined by the parameter $N_{s}$. In instances where $N_{s}=0$,
the KRBH converges to a Kerr black hole. However, for $N_{s}\neq 0$ and
$-1<w_{s}<-1/3$, the configuration yields a Kerr black hole encompassed by a
quintessential fluid when $\psi=0$. When $N_{s}=a=0$, the spacetime reverts to
that of a Schwarzschild black hole. It is noteworthy that akin to Kerr black
holes, the KRBH spacetime exhibits two Killing vectors $\bm{\xi}_{t}$ and
$\bm{\xi}_{\phi}$, aligning with symmetries respecting the $t$ and
$\phi$-coordinates.

We proceed by examining the crucial hypersurfaces that delineate the causal
boundaries within the KRBH spacetime.

###  II.1 The causal structure of the KRBH

Based on the specific spacetime metric, it is evident that the formation of
horizons around the KRBH diverges from that in general relativistic
spacetimes. This divergence arises due to contributions from the structure
parameter $N_{s}$ and the Rastall coupling $\psi$. Equation (5a) indicates the
potential existence of up to three horizons within the spacetime. However, the
black hole’s ability to manifest these horizons hinges strictly upon the value
of $\zeta$ within the function $\Delta$, determined by both $w_{s}$ and
$\psi$. In the subsequent subsection, we explore various specific cases,
discussing the availability of horizons and furnishing analytical solutions
for the radii of the characteristic hypersurfaces within the spacetime.

####  II.1.1 The case of $\zeta=1$

The polynomial equation $\Delta=0$ identifies the black hole horizons, which,
when $\zeta=1$, leads to a quadratic equation. Solving this equation yields
the analytical solutions

$r_{\mp}={2M+N_{s}\over 2}\mp\sqrt{\left({2M+N_{s}\over 2}\right)^{2}-a^{2}},$
(7)

which corresponds, respectively, to the Cauchy and event horizons of the black
hole, each present distinctly when the condition
$4a^{2}<\left({2M+N_{s}}\right)^{2}$ holds true. When
$a=(2M+N_{s})/2=a_{\mathrm{ext}}$, an extremal KRBH forms, where
$r_{+}=r_{-}=M+N_{s}/2=r_{\mathrm{ext}}$. However, for
$4a^{2}>\left({2M+N_{s}}\right)^{2}$, a naked singularity emerges. Notably, in
this specific case, no cosmological horizons are observed. Furthermore, from
Eq. (6), specific values of $\zeta$ impose constraints on ranges for $w_{s}$
and $\psi$, determining the maximum number of solutions to the equation
$\Delta=0$. The left panel of Fig. 1 illustrates these domains for $w_{s}$ and
$\psi$ within a diagram considering $\zeta=1$.

(a)  (b)

Figure 1: The diagrams depict (a) the interrelation between $\psi$ and
$w_{s}$, and (b) the radial profile of $\Delta(r)$ for three distinct spin
parameters, specifically for $\zeta=1$, with $N_{s}=0.05$. Throughout this
paper, including these diagrams and all subsequent ones, the unit of length
along the axes is standardized to $M$.

According to this diagram, the allowed range of the Rastall coupling, within
the domain $-1<w_{s}<-1/3$, is $-\infty<\psi<-0.5$. The plot of the radial
profile of the function $\Delta$ in the right panel of Fig. 1 confirms this
anticipated outcome that the number of real solutions remains at either two
(black hole), one (extremal black hole), or zero (naked singularity). Thus,
despite the presence of a quintessential field around the black hole, a
cosmological horizon is not observed for this specific value of $\zeta$.

It is important to note that apart from horizons, stationary black holes also
feature other hypersurfaces where static observers, defined by the world-lines
$\bm{u}=(-g_{tt})^{-1/2}\bm{\xi}_{t}$, cannot exist. These particular
hypersurfaces are termed as the static limits, where $\bm{u}$ becomes null.
The radii of these static limits are determined by the equation $g_{tt}=0$.
For the KRBH, utilizing Eq. (4), the radii of the static limits are derived as

$r_{\mathrm{SL}_{\mp}}(\theta)={2M+N_{s}\over 2}\mp\sqrt{\left({2M+N_{s}\over
2}\right)^{2}-a^{2}\cos^{2}\theta},$ (8)

which are located, respectively, inside the Cauchy horizon and outside the
event horizon. The spaces between the static limits and the horizons are known
as the ergoregions, and correspond to the domains, in which, all observers are
in the state of corotation with the black hole. Considering the common
Cartesian coordinates $X=r\sin\theta\cos\phi$, $Y=r\sin\theta\sin\phi$, and
$Z=r\cos\theta$, in Fig. 2, we have demonstrated the ergoregions for the case
of $\zeta=1$.

(a)  (b)  (c)

Figure 2: The cross-section of the ergoregions of the KRBH in the $Z$-$X$
plane, for (a) black hole with $a=0.8$, (b) extremal black hole with
$a=a_{\mathrm{ext}}=1.025$, and (c) naked singularity with $a=1.2$, for the
case of $\zeta=1$. The red and green dashed circles correspond, respectively,
to $r_{\mp}$ and $r_{\mathrm{ext}}$.

Clearly, as the spin parameter escalates, the radial extent of the ergoregion
diminishes. In scenarios involving an extremal black hole with a solitary
horizon at $r_{\mathrm{ext}}$, it is noteworthy that
$r_{\mathrm{SL}_{-}}^{0}=r_{\mathrm{SL}_{+}}^{0}=r_{\mathrm{SL}_{\mp}}(0)|_{a=a_{\mathrm{ext}}}$.
Despite the horizon vanishing for the naked singularity, the ergoregion
persists, albeit bifurcating into two distinct branches.

####  II.1.2 The case of $\zeta=2$

In this case, the horizons are found at

$r_{\mp}={M\over 1-N_{s}}\mp\sqrt{\left({M\over
1-N_{s}}\right)^{2}-{a^{2}\over 1-N_{s}}},$ (9)

which implies that the existence of the KRBH corresponds to the condition
$a^{2}(1-N_{s})<M^{2}$, and hence, $a_{\mathrm{ext}}=M/\sqrt{1-N_{s}}$,
implying that $r_{\mathrm{ext}}=M/(1-N_{s})$. As demonstrated in Fig. 3 for a
quintessential background fluid, the above choice of $\zeta$ confines the
Rastall coupling inside the domain $-\infty<\psi<0$.

(a)  (b)

Figure 3: (a) The relationship between $\psi$ and $w_{s}$, and (b) the radial
structure of $\Delta(r)$ exhibited across three distinct spin parameters,
under the setting of $\zeta=2$ and $N_{s}=0.05$.

It is evident that there are no cosmological horizons in this particular
scenario. In this chosen $\zeta$ configuration, the ergoregions are delineated
by the radii of the static limits, expressed as

$r_{\mathrm{SL}_{\mp}}(\theta)={M\over 1-N_{s}}\mp\sqrt{\left({M\over
1-N_{s}}\right)^{2}-{a^{2}\cos^{2}\theta\over 1-N_{s}}}.$ (10)

(a)  (b)  (c)

Figure 4: The cross-section of the ergoregions of the KRBH for three different
spin parameters, corresponding to the choice of $\zeta=2$. In this particular
case, $a_{\mathrm{ext}}=1.026$.

In Fig. 4, the cross-section of the corresponding ergoregions has been plotted
in the $Z$-$X$ plane.

####  II.1.3 The case of $\zeta=3$

The cosmological (quintessential) horizon first emerges when $\zeta=3$. This
specific scenario is particularly noteworthy as it allows us to capture all
feasible real solutions of the cubic equation $\Delta=0$. This results in the
three solutions

$\displaystyle
r_{+}=\sqrt{\frac{g_{2}}{3}}\cos\left(\frac{1}{3}\arccos\left(3g_{3}\sqrt{\frac{3}{g_{2}^{3}}}\,\right)+\frac{4\pi}{3}\right)+\frac{1}{3N_{s}},$
(11) $\displaystyle
r_{++}=\sqrt{\frac{g_{2}}{3}}\cos\left(\frac{1}{3}\arccos\left(3g_{3}\sqrt{\frac{3}{g_{2}^{3}}}\,\right)\right)+\frac{1}{3N_{s}},$
(12) $\displaystyle
r_{-}=\sqrt{\frac{g_{2}}{3}}\cos\left(\frac{1}{3}\arccos\left(3\sqrt{\frac{3}{g_{2}^{3}}}\,\right)+\frac{2\pi}{3}\right)+\frac{1}{3N_{s}},$
(13)

where

$\displaystyle g_{2}=4\left(\frac{1}{3N_{s}^{2}}-\frac{2M}{N_{s}}\right),$
(14a) $\displaystyle
g_{3}=4\left(\frac{a^{2}}{N_{s}}-\frac{2M}{3N_{s}^{2}}+\frac{2}{27N_{s}^{3}}\right).$
(14b)

It is straightforward to verify that the discriminant of the cubic is always
non-negative, i.e. $\delta_{\Delta}=16\left(g_{2}^{3}-27g_{3}^{2}\right)\geq
0$. Consequently, all the aforementioned roots are real. Moreover, as depicted
in Fig. 5(a), within a quintessential fluid background, the Rastall coupling
is constrained within the domain $-\infty<\psi<1/6=0.167$. Within this range,
all three zeros of $\Delta(r)=0$ will be positive (see Fig. 5(b)).

(a)  (b)

Figure 5: (a) The correlation between $\psi$ and $w_{s}$, and (b) the radial
profile of $\Delta(r)$ for three distinct spin parameters, for the case of
$\zeta=3$ and $N_{s}=0.1$.

Accordingly, the radii (11)–(12) indicate, respectively, the event,
cosmological, and Cauchy horizons of the KRBH for the above choice of $\zeta$.
Hence, the extremal KRBH black hole with two horizons (a cosmological horizon
and an event horizon) can be identified in terms of the equation
$\delta_{\Delta}=0$. This equation provides

$a^{2}_{\mathrm{ext}}=\frac{2}{27}\left(\frac{\sqrt{(1-6MN_{s})^{3}}}{N_{s}^{2}}+\frac{9M}{N_{s}}-\frac{1}{N_{s}^{2}}\right),$
(15)

as the extremal spin parameter. Therefore, the KRBH with three horizons is
present when $a\geq a_{\mathrm{ext}}$, while the naked singularity with a
solitary cosmological horizon emerges when $a\leq a_{\mathrm{ext}}$ (refer to
Fig. 5(b)). Similar to previous cases, the static limit can also be
identified. However, for $\zeta=3$, an additional ergoregion forms within the
cosmological horizon. Building upon the findings of the prior two cases, it is
evident that the radii of the static limits are $r_{\mathrm{SL}_{++}}$,
$r_{\mathrm{SL}_{+}}$, and $r_{\mathrm{SL}_{-}}$. These radii share identical
expressions as $r_{++}$, $r_{+}$, and $r_{-}$, respectively, considering the
exchange $a^{2}\rightarrow a^{2}\cos^{2}\theta$ within the $g_{3}$ coefficient
in Eq. (14b). These expressions are employed in Fig. 6 to illustrate the
ergoregions for the three instances of the KRBH, the extremal KRBH, and the
naked singularity.

(a)  (b)  (c)

Figure 6: The cross-section of the ergoregions of the KRBH for three different
spin parameters, corresponding to the choice of $\zeta=3$. In panel (a), the
three horizons are available separately for the KRBH and have been indicated
by red dashed circles. Hence there will be three ergoregions. The extremal
case is shown in panel (b) for $a_{\mathrm{ext}}=1.065$, where $r_{-}=r_{+}$
and $r_{++}$ have been indicated by green dashed circles. There will be two
ergoregions in this case. Finally, in panel (c) a naked singularity is shown
with only $r_{++}$ being available, which is indicated by a purple dashed
circle. As in the previous cases, there is an ergoregion connected to the
singularity, but in this case, there will be an exterior ergoregion formed
together with the cosmological horizon.

It is evident that the choice of $\zeta=3$ provides more diversity in the
spacetime structure and it is of more completeness. Hence, in the rest of the
paper, we consider this choice as our reference and we base our analytical and
numerical studies on the metric function

$\Delta(r)=r^{2}+a^{2}-2Mr-
N_{s}r^{3}=N_{s}\left(r_{++}-r\right)\left(r-r_{+}\right)\left(r-r_{-}\right),$
(16)

in accordance with the characteristic radii given in Eqs. (11)–(13). Moving
forward, we turn our attention to the inclusion of a plasmic medium around the
KRBH and proceed to present the corresponding equations of motion in the
following section.

##  III Light propagation in a non-magnetized plasma

In this section, we introduce some salient aspects of plasma effects on the
light ray trajectories that pass a black hole spacetime. In fact, we regard
the plasma as a medium that is only discernible by its influence on the rays
so the wave optics descriptions are neglected here. It is well-known that
inside a plasma, light rays do not travel on null geodesics, and in fact, it
is the frequency perlick_calculating_2022

$\omega_{p}(\bm{x})=\frac{4\pi e^{2}}{m_{e}}N_{p}(\bm{x}),$ (17)

of an electron-ion plasma, that affects the propagation of light in such a
medium. In the above relation, $e$, $m_{e}$, and $N_{p}$ are, respectively,
the electron charge, electron mass $m_{e}$, and the electron number density
$N_{p}$, in the coordinate system $\bm{x}\equiv
x^{\mu}=(x^{0},x^{1},x^{2},x^{3})$. Here, we assume that the plasma is non-
magnetized and pressure-less so the only factor that affects the light rays is
the electron density. In order to formulate the equations of motion for light
rays that propagate inside this medium, we rely on the phase space dynamics
and introduce the Hamiltonian

$\mathcal{H}(\bm{x},\bm{p})=\frac{1}{2}\left[g^{\mu\nu}(\bm{x})p_{\mu}p_{\nu}+\omega_{p}^{2}(\bm{x})\right],$
(18)

defined on the cotangent bundle of the Riemannian manifold
$(\mathcal{M},\bm{g})$, with $\bm{p}\equiv p_{\mu}=(p_{0},p_{1},p_{2},p_{3})$
being the canonical momentum covector. In this regard, the propagation of the
light rays is identified as solution to the canonical equations

$\displaystyle{\mathrm{d}x^{\mu}\over\mathrm{d}\lambda}=\frac{\partial\mathcal{H}}{\partial
p_{\mu}},$ (19)
$\displaystyle{\mathrm{d}p_{\mu}\over\mathrm{d}\lambda}=-\frac{\partial\mathcal{H}}{\partial
x^{\mu}},$ (20) $\displaystyle\mathcal{H}=0,$ (21)

in which, $\lambda$ is the curve parameter along the trajectories. In fact,
Breuer and Ehler were the first to derive the Hamiltonian (18) from Maxwell’s
equations in a general relativistic spacetime, while using a model of the
plasma involving two charged fluids Breuer:1980 .

Since the light rays travel on time-like trajectories with respect to
spacetime manifold’s metric, it is advantageous to consider a comoving
observer with the four-velocity $\bm{u}$, who measures the light ray frequency
$\omega(\bm{x})$, and characterize the plasma by the refractive index

$n^{2}(\bm{x})=1+\frac{p_{\mu}p^{\mu}}{(p_{\nu}u^{\nu})^{2}}=1-\frac{\omega_{p}^{2}(\bm{x})}{\omega^{2}(\bm{x})},$
(22)

in writing which, we have used the identity
$p_{\nu}u^{\nu}=p_{0}/\sqrt{g_{00}}=-\omega$, corresponding to an observer
located on the $x^{0}$ curves, having the four-velocity
$u^{\mu}={\delta_{0}^{\mu}}/{\sqrt{-g_{00}}}$. This way, one can recast the
Hamiltonian as

$\mathcal{H}=\frac{1}{2}\left[g^{\mu\nu}p_{\mu}p_{\nu}+(n^{2}-1)(p_{\nu}u^{\nu})^{2}\right]=\frac{1}{2}\mathfrak{o}^{\mu\nu}p_{\mu}p_{\nu},$
(23)

where $\mathfrak{o}^{\mu\nu}$ is the optical metric associated with the the
refractive index $n\equiv n(\vec{x},\omega)$. Contributing the Jacobi action
$\mathcal{S}$, with the properties

$\displaystyle p_{\mu}=\frac{\partial\mathcal{S}}{\partial x^{\mu}},$ (24a)
$\displaystyle\mathcal{H}=\frac{\partial\mathcal{S}}{\partial\lambda},$ (24b)

one can rewrite the Hamilton-Jacobi equations as Atamurotov:2015 ;
Perlick:2017

$\displaystyle\frac{\partial\mathcal{S}}{\partial\lambda}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\left[g^{\mu\nu}\frac{\partial\mathcal{S}}{\partial
x^{\mu}}\frac{\partial\mathcal{S}}{\partial
x^{\nu}}-\left(n^{2}-1\right)\omega^{2}\right]$ (25) $\displaystyle=$
$\displaystyle-\frac{1}{2}\left[g^{\mu\nu}\frac{\partial\mathcal{S}}{\partial
x^{\mu}}\frac{\partial\mathcal{S}}{\partial x^{\nu}}+\omega_{p}^{2}\right].$

In fact, the study of light propagation in dispersive media with frequency-
dependent refractive indices within the framework of general relativity was
pioneered by Synge Synge:1960 , although he did not consider the case of the
plasma (see also Refs. Bicak:1975 ; PhysRevD.45.525 ). If the spacetime
exhibits spherical symmetry, the plasmic frequency will vary with respect to
the radial distance $\omega_{p}\equiv\omega_{p}(r)$, and hence, the
propagation of light can be determined through appropriate integrations
Muhleman:1966 ; Perlick_ray_2000 . The gravitational deflection of light in
homogeneous plasma has been studied in Refs. bisnovatyi-
kogan_gravitational_2009 ; bisnovatyi-kogan_gravitational_2010 ; Tsupko:2013 .
It has been shown in Ref. perlick_light_2017  that the separation of the
Hamilton-Jacobi equations in the Kerr spacetimes is possible, only if the
plasma frequency adopts the form

$\omega_{p}^{2}\equiv\omega_{p}^{2}(r,\theta)={f_{r}(r)+f_{\theta}(\theta)\over
r^{2}+a^{2}\cos^{2}\theta},$ (26)

in which $f_{r}$ and $f_{\theta}$ can have arbitrary expressions. As we will
see in what follows, this holds, as well, for the case of the KRBH. Such an
expression has been considered in Ref. fathi_analytical_2021  to analyze the
light ray trajectories around a Kerr black hole in the presence of an
inhomogeneous plasma. Note that, the form (26) is also of crucial importance
in the analytical determination of the black hole shadow boundary, as shown in
Refs. perlick_light_2017 ; Yan:2019 . On the other hand, when the separation
of the Hamilton-Jacobi equation is not at hand, one has to make numerical
implementation in order to determine the shadow huang_revisiting_2018 .
However, since the current discussion mainly deals with the analytical
description of the light ray trajectories rather than the shadow, we do not
continue this review and proceed with the determination of the equations of
motion.

###  III.1 Separation of the Hamilton-Jacobi equation and the equations of
motion

The method of the separation of the Hamilton-Jacobi equation is based on
rewriting the action $\mathcal{S}$ as Carter:1968 ; Chandrasekhar:2002

$\mathcal{S}=-Et+L\phi+\mathcal{S}_{r}(r)+\mathcal{S}_{\theta}(\theta)+\frac{1}{2}m^{2}\lambda,$
(27)

in which, the energy $E$ and the angular momentum $L$ are the constants of
motion specific to the Hamilton equations, and $m$ is the mass of the test
particles (here, $m=0$). In fact, applying the line element (4) in the
Hamiltonian (18), we get

$\mathcal{H}(\bm{x},\bm{p})={1\over
2\rho^{2}}\left[\Delta\,p_{r}^{2}+p_{\theta}^{2}+\left(ap_{t}\sin\theta+{p_{\phi}\over\sin\theta}\right)^{2}-{1\over\Delta}\Big{(}p_{t}\left(r^{2}+a^{2}\right)+ap_{\phi}\Big{)}^{2}+\rho^{2}\omega_{p}^{2}\right],$
(28)

in which, the components of the four-momentum are determined in alignment with
Eqs. (24a) and (27), taking into account the expressions

$\displaystyle p_{t}=\frac{\partial\mathcal{S}}{\partial t}=-E=-\omega_{0},$
(29a) $\displaystyle p_{\phi}=\frac{\partial\mathcal{S}}{\partial\phi}=L.$
(29b)

for the constants of motion. It is then clear that $p_{\phi}$ corresponds to
the axial symmetry of the spacetime, while $\omega_{0}$ could be regarded as
the light rays’ frequency. Hence, the constants of motion $E$ and $L$,
together with the eikonal equation $\mathcal{H}=0$, characterize the
propagation of light. Now from Eqs. (28) and (29), the eikonal equation takes
the form

$0=\Delta\,p_{r}^{2}+\,p_{\theta}^{2}+\left(\omega_{0}\,a\sin\theta-{L\over\sin\theta}\right)^{2}-{1\over\Delta}\left[\omega_{0}\left(r^{2}+a^{2}\right)-aL\right]^{2}+\rho^{2}\omega_{p}^{2},$
(30)

from which, it is evident that the separability to the $r$-dependent and
$\theta$-dependent segments, is only possible when the condition (26) is
respected. By exploiting the identity

$\left(\omega_{0}\,a\sin\theta-{L\over\sin\theta}\right)^{2}=\left(L^{2}\csc^{2}\theta-a^{2}\omega_{0}^{2}\right)\cos^{2}\theta\\\
+(L-a\omega_{0})^{2},$ (31)

together with Eqs. (30) and (26), we get the separated form of the Hamilton-
Jacobi equation, which reads as

$\displaystyle\mathscr{Q}$ $\displaystyle=$
$\displaystyle\,p_{\theta}^{2}+\left(L^{2}\csc^{2}\theta-a^{2}\omega_{0}^{2}\right)\cos^{2}\theta+f_{\theta}(\theta)$
(32) $\displaystyle=$
$\displaystyle-\Delta\,p_{r}^{2}+{1\over\Delta}\left[\omega_{0}\left(r^{2}+a^{2}\right)-aL\right]^{2}-(L-a\omega_{0})^{2}-f_{r}(r),$

where $\mathscr{Q}$ is the Carter’s constant. Accordingly, one can write

$\displaystyle
p_{\theta}^{2}=\mathscr{Q}-\left(L^{2}\csc^{2}\theta-a^{2}\omega_{0}^{2}\right)\cos^{2}\theta-
f_{\theta}(\theta),$ (33)
$\displaystyle\Delta\,p_{r}^{2}={1\over\Delta}\left[\omega_{0}\left(r^{2}+a^{2}\right)-aL\right]^{2}-\mathscr{Q}-(L-a\omega_{0})^{2}-f_{r}(r).$
(34)

Incorporating Eqs. (33) and (34) into the Hamiltonian (28), and then applying
Eq. (19), yields the first order differential equations of motion as follows:

$\displaystyle\rho^{2}\frac{\mathrm{d}r}{\mathrm{d}\lambda}=\sqrt{\mathcal{R}(r)},$
(35)
$\displaystyle\rho^{2}\frac{\mathrm{d}\theta}{\mathrm{d}\lambda}=\sqrt{\Theta(\theta)},$
(36)
$\displaystyle\rho^{2}\frac{\mathrm{d}\phi}{\mathrm{d}\lambda}={L\left(\Delta-a^{2}\sin^{2}\theta\right)\csc^{2}\theta+a\omega_{0}\left(r^{2}+a^{2}-\Delta\right)\over\Delta},$
(37)
$\displaystyle\rho^{2}\frac{\mathrm{d}t}{\mathrm{d}\lambda}={\omega_{0}\Sigma^{2}-aL\left(r^{2}+a^{2}-\Delta\right)\over\Delta},$
(38)

where

$\displaystyle\mathcal{R}(r)=\left[\omega_{0}\left(r^{2}+a^{2}\right)-aL\right]^{2}-\Delta\left[\mathscr{Q}+f_{r}(r)+(L-a\omega_{0})^{2}\right],$
(39a)
$\displaystyle\Theta(\theta)=\mathscr{Q}-f_{\theta}(\theta)-\cos^{2}\theta\left(L^{2}\csc^{2}\theta-a^{2}\omega_{0}^{2}\right),$
(39b)
$\displaystyle\Sigma^{2}=\left(r^{2}+a^{2}\right)^{2}-\Delta\,a^{2}\sin^{2}\theta.$
(39c)

However, in order to facilitate the mathematical manipulations in what
follows, we work with the dimension-less Mino time, $\gamma$, which is defined
in terms of the relation $\rho^{2}\mathrm{d}\gamma=M\mathrm{d}\lambda$
Mino:2003 . In this context, the equations of motion adopt the forms

$\displaystyle M{\mathrm{d}r\over\mathrm{d}\gamma}$ $\displaystyle=$
$\displaystyle\sqrt{\mathcal{R}(r)},$ (40) $\displaystyle
M{\mathrm{d}\theta\over\mathrm{d}\gamma}$ $\displaystyle=$
$\displaystyle\sqrt{\Theta(\theta)},$ (41) $\displaystyle
M{\mathrm{d}\phi\over\mathrm{d}\gamma}$ $\displaystyle=$
$\displaystyle{L\left(\Delta-a^{2}\sin^{2}\theta\right)\csc^{2}\theta+a\omega_{0}\left(r^{2}+a^{2}-\Delta\right)\over\Delta},$
(42) $\displaystyle M{\mathrm{d}t\over\mathrm{d}\gamma}$ $\displaystyle=$
$\displaystyle{\omega_{0}\Sigma^{2}-aL\left(r^{2}+a^{2}-\Delta\right)\over\Delta}.$
(43)

We extend our discussion in the next section by implementing a specific
criterion for the frequency $\omega_{p}$. In line with the form (26), we
construct an inhomogeneous anisotropic plasma. Subsequently, we analyze the
aforementioned equations of motion to obtain precise analytical solutions for
the evolution of the coordinates. These solutions determine the trajectories
of light rays around the KRBH, considering the case where $\zeta=3$ as
previously mentioned in Sect. II.

##  IV General analysis of the equations of motion

The separation condition stated in Eq. (26) defines the plasma and its
geometric distribution in the spacetime surrounding the black hole.
Consequently, the characteristic functions $f_{r}(r)$ and $f_{\theta}(\theta)$
assume significant roles in delineating the configuration of the plasma. In
this section, consistent with the approach outlined in Ref.
fathi_analytical_2021 , we constrain the light rays to traverse through an
inhomogeneous anisotropic plasma by making the specific choices

$\displaystyle f_{r}(r)\equiv f_{r}=\mathrm{constant},$ (44) $\displaystyle
f_{\theta}(\theta)\equiv f_{\theta}=\mathrm{constant.}$ (45)

Furthermore, within the text, we use, frequently, the conventions
Chandrasekhar:2002

$\displaystyle\xi={L\over\omega_{0}},$ (46)
$\displaystyle\eta={\mathscr{Q}\over\omega_{0}^{2}},$ (47)

in order to simplify the analysis. In what follows, the temporal evolution of
the spacetime coordinates are analyzed separately for this plasmic medium, and
the relevant exact solutions are presented.

###  IV.1 Latitudinal motion

To analyze the evolution of the polar coordinate $\theta$, we need to examine
the angular potential $\Theta(\theta)$. To proceed, we apply the change of
variable $z=\cos^{2}\theta$, so that the potential adopts the form

$\Theta_{z}(z)=\frac{\omega_{0}^{2}}{1-z}\left[-a^{2}z^{2}+\tilde{\chi}z+\tilde{\eta}\right],$
(48)

in which we have notated $\tilde{\eta}=\eta-f_{\theta}/\omega_{0}^{2}$ and
$\tilde{\chi}=a^{2}-\tilde{\eta}-\xi^{2}<0$. This potential is consequently
second order in $z$, and its roots are obtained as

$z_{\pm}=\frac{\tilde{\chi}\pm\sqrt{4a^{2}\tilde{\eta}+\tilde{\chi}^{2}}}{2a^{2}},$
(49)

providing $z_{+}z_{-}=-\tilde{\eta}/a^{2}$. This way, one can recast the
angular potential as

$\Theta_{z}(z)=\frac{\omega_{0}^{2}}{1-z}\left[a^{2}(z_{+}-z)(z-z_{-})\right].$
(50)

Accordingly, the evolution of the polar angle may encounter four turning
points, which are given by

$\displaystyle\theta_{1}=\arccos(\sqrt{z_{+}}),$ (51)
$\displaystyle\theta_{2}=\arccos(\sqrt{z_{-}}),$ (52)
$\displaystyle\theta_{3}=\arccos(-\sqrt{z_{-}}),$ (53)
$\displaystyle\theta_{4}=\arccos(-\sqrt{z_{+}}).$ (54)

The nature of the roots becomes evident through the positivity or negativity
of $z_{\pm}$ or the sign of $\tilde{\eta}$. Therefore, in the subsequent
discussion, we categorize the latitudinal motion into three cases.

  * (i)

$\tilde{\eta}>0$ (ordinary trajectories): In this case $z_{+}z_{-}<0$. On the
other hand, it is immediately inferred from Eq. (49) that $z_{+}>0$ and
$z_{-}<0$. Accordingly, $\theta_{2,3}\in\mathbb{C}$ and
$\theta_{1,4}\in\mathbb{R}$. By analyzing the latter case, we find that
$\theta_{1}<\pi/2<\theta_{4}$, and hence, the polar coordinate oscillates
between these two limits in the course of its evolution, and passes the
equatorial plane in each oscillation.

  * (ii)

$\tilde{\eta}=0$ (planar trajectories): It is straightforward to see that, in
this case, $z_{+}=0$ and $z_{-}=\tilde{\chi}/a^{2}<0$. Hence
$\theta_{1}=\theta_{4}=\pi/2$ and $\theta_{2,3}\in\mathbb{C}$. In this regard,
the orbits are confined in the equatorial plane.

  * (iii)

$\tilde{\eta}<0$ (vortical trajectories): Under this condition,
$z_{+}z_{-}>0$, which based on the fact that $\tilde{\chi}<0$, then it holds
$z_{\pm}<0$, and all the roots $\theta_{i}$ are imaginary. Hence, this case is
not physically valid and is ruled out from consideration.

Accordingly, to obtain the analytical solution to the evolution of the
$\theta$-coordinate, we consider the case of $\tilde{\eta}\geq 0$. Assuming a
source located at the angular position $\theta_{\mathrm{s}}$ and an observer
at $\theta_{\mathrm{o}}$, the Mino time can be obtained by means of the
equation of motion (41) and the angular potential (50), which after proper
inversion, results in the analytical evolution relation

$\theta_{\mathrm{o}}(\gamma)=\arccos\bigg{(}-\sqrt{z_{+}}\operatorname{\mathbf{sn}}\Big{(}\frac{z_{0}}{M}\gamma-\operatorname{\mathbf{F}}(\varphi_{\mathrm{s}}|\mathcal{k})\Big{|}\mathcal{k}\Big{)}\bigg{)},$
(55)

in which $z_{0}=a\omega_{0}\sqrt{-z_{-}}$, $\mathcal{k}=z_{+}/z_{-}$ is the
elliptic modulus, $\operatorname{\mathbf{sn}}(x|\mathcal{k})$ is the Jacobian
elliptic sine function with argument $x$,
$\operatorname{\mathbf{F}}(\varphi_{s}|\mathcal{k})$ is the incomplete
elliptic integral of the first kind with argument $\varphi_{s}$
byrd_handbook_1971 , and

$\varphi_{\mathrm{s}}=\arcsin\left(\frac{\cos{\theta_{\mathrm{s}}}}{\sqrt{z_{+}}}\right).$
(56)

Based on the periodicity of the elliptic sine function, as expected, the
$\theta$-coordinate will oscillate between its maximum and minimum values in
the course of full trajectories. This has been shown in Fig. 7 for some
specific examples.

(a)  (b)

Figure 7: The oscillating evolution of the $\theta$-coordinate in the Mino
time, by considering $a=0.8$, $L=1$ and $\omega_{0}=0.75$ (corresponding to
$\xi=1.33$), and $\theta_{\mathrm{s}}=\pi/2$. Panel (a) shows a table of
graphs for $\theta_{o}(\gamma)$ for the impact parameter range
$\tilde{\eta}\in[0,2]$ by $0.1$ steps. As it is evident, by decreasing
$\tilde{\eta}$, the oscillations approach the $\theta=\pi/2$ line until they
finally stop oscillating (planar orbits). In panel (b) a typical graph has
been plotted by taking into account $\mathscr{Q}=9$ and $f_{\theta}=1$
(corresponding to $\tilde{\eta}=17.77$). In this diagram, the roots of the
angular potential have been shown as the maximum and minimum values of the
$\theta$-coordinate oscillations.

###  IV.2 Radial motion

In the analysis of particle trajectories in curved spacetimes, a fundamental
aspect involves comprehending the approach and recession of particles
influenced by gravity. Traditionally, this exploration revolves around
computing the radial effective gravitational potential governing these
particles. Consequently, when scrutinizing the motion of light rays in our
specified context, our focus shifts to the radial equation of motion in Eq.
(40). To advance our study, this subsection concentrates on solving the
aforementioned evolution equation under different conditions, leading to
various types of motion. We will discuss, separately, the deflecting
trajectories (both bound and unbound), captured trajectories, and the critical
orbits (both bound and unbound). However before delving into the analytical
solution of the radial equation of motion, let us consider some aspects of
spherical photon orbits around the black hole.

####  IV.2.1 Orbits of constant radius around the KRBH

Around black holes with a non-zero spin parameter, a multitude of light rays
follow trajectories where the radial distance $r$, remains constant. These
rays essentially skim the surface of a spherical structure, and carry a
critical nature, given that they either escape from the black hole or plunge
into it (as further explained below). Consequently, they constitute a major
portion of a more general form of orbits, termed as the critical orbits (see
below). Within these orbits, the $\theta$-coordinate oscillates between the
turning points of the angular potential (as elaborated in the last
subsection). For Kerr black holes, these particular orbits have been
extensively examined in standard texts such as Refs. Chandrasekhar:2002 ;
Bardeen:1972a ; Bardeen:1973b . These photon orbits are inherently unstable,
and aligned with the number of half orbits around the black hole, they
collectively constitute an infinite number of photon rings in the observer’s
sky that converge to a critical curve confining the black hole shadow, as
indicated in Gralla:2019 . From a theoretical standpoint, determining
spherical orbits for static spacetimes (like Schwarzschild) is
straightforward, as they are planar and essentially circles on the equatorial
plane. However, this simplicity changes for stationary (rotating) black holes
due to the frame-dragging effect. This effect creates a region within a
specific radial range around the black hole, bound by the innermost and
outermost planar orbits that follow circular paths. Essentially, the remaining
orbits are non-planar. The region populated by spherical photon orbits is
appropriately termed the photon region, representing the set of points that
these orbits traverse, forming a closed boundary. For Kerr black holes and by
applying the geodesic equations for light rays, the radii of planar and polar
orbits were initially provided in Refs. Bardeen:1972a ; Bardeen:1973b , and
discussed further in the context of Kerr geodesics in Ref. Chandrasekhar:2002
. Ever since, numerous publications have been devoted to the study of
spherical photon orbits, photon regions and photon rings in Kerr and Kerr-like
black hole spacetimes (see for example Refs. stoghianidis_polar_1987 ;
cramer_using_1997 ; Teo:2003 ; Johannsen:2013 ; Grenzebach:2014 ; Perlick:2017
; charbulak_spherical_2018 ; Johnson_universal_2020 ; Himwich:2020 ;
Gelles:2021 ; Ayzenberg:2022 ; Das:2022 ; fathi_spherical_2023 ;
ANJUM2023101195 ; Chen:2023 ; andaru_spherical_2023 ).

In fact, the spherical photon orbits form under specific conditions which are
established by the equations $\mathcal{R}(r)=0=\mathcal{R}^{\prime}(r)$
Teo:2003 . Given that from Eq. (39a) we can rewrite

$\mathcal{R}(r)=\omega_{0}^{2}\Big{(}a^{2}+r^{2}-a\xi\Big{)}^{2}-\omega_{0}^{2}\Delta(r)\Big{[}(\xi-a)^{2}+\eta_{r}\Big{]},$
(57)

with $\eta_{r}=\eta+f_{r}/\omega_{0}^{2}$, these mutual conditions provide the
criteria

$\displaystyle\xi^{c}=\frac{a^{2}+r^{2}}{a},$ (58)
$\displaystyle\eta_{r}^{c}=-\frac{r^{4}}{a^{2}},$ (59)

or

$\displaystyle\xi^{c}=\frac{a^{2}[2M+r(3N_{s}r+2)]-r^{2}[6M+r(N_{s}r-2)]}{a[2M+r(3N_{s}r-2)]},$
(60)
$\displaystyle\eta_{r}^{c}=-\frac{r^{3}\left[8a^{2}\left(N_{s}r^{2}-2M\right)+r\left[6M+r(N_{s}r-2)\right]^{2}\right]}{a^{2}[2M+r(3N_{s}r-2)]^{2}},$
(61)

for the case of $\zeta=3$ in Eq. (5a). These expressions indicate the critical
impact parameters for light rays on $r$-constant orbits. For the case of Eq.
(59), it is clear that $\eta^{c}_{r}<0$, and in accordance with the statements
in the classification (iii) of the $\theta$-motion in the previous subsection,
confronts us with a non-physical problem. Hence, we disregard this case and
consider the second set of the critical impact parameters in Eqs. (60) and
(61).

Recalling that for planar orbits, the condition $\theta=0$ holds, which
together with the condition of spherical orbits, leads us to the condition
$\eta_{r}^{c}=(f_{r}+f_{\theta})/{\omega_{0}^{2}}$. This condition corresponds
to spherical planar orbits (i.e., circular orbits) for photons. Utilizing Eq.
(61), this leads to the octic equation.

$P_{8}(r)=\sum_{j=0}^{8}b_{j}r^{j}=0,$ (62)

where

$\displaystyle b_{0}=4a^{2}M^{2}\left(f_{r}+f_{\theta}\right),$ (63a)
$\displaystyle b_{1}=8a^{2}M\left(f_{\theta}-f_{r}\right),$ (63b)
$\displaystyle
b_{2}=4a^{2}\left(f_{r}+f_{\theta}\right)\left(1+3MN_{s}\right),$ (63c)
$\displaystyle
b_{3}=4a^{2}\left[3N_{s}\left(f_{\theta}-f_{r}\right)-4M\omega_{0}^{2}\right],$
(63d) $\displaystyle
b_{4}=9a^{2}N_{s}^{2}\left(f_{r}+f_{\theta}\right)+36M^{2}\omega_{0}^{2},$
(63e) $\displaystyle b_{5}=8\omega_{0}^{2}\left(a^{2}N_{s}-3M\right),$ (63f)
$\displaystyle b_{6}=4\omega_{0}^{2}\left(1+3MN_{s}\right),$ (63g)
$\displaystyle b_{7}=-4N_{s}\omega_{0}^{2},$ (63h) $\displaystyle
b_{8}=N_{s}^{2}\omega_{0}^{2}.$ (63i)

Although such polynomial equations cannot be solved in terms of simple
radicals, they however have been treated by means of approximate methods (see
for example Ref. Tavlayan:2020  for the case of spherical photon orbits in
Kerr spacetime). On the other hand, since this study is aimed at general
analytical solutions to the equations of motion, we do not deal with solutions
to the above octic, since its numerical solutions suffice for our purpose.
This way, one can verify that the octic (62) has only two positive solutions
$r_{\mathrm{ph}\mp}$ which present, respectively, the inner prograde, and the
outer retrograde radii of circular photon orbits. These two radii constitute
the boundary of the radii of spherical orbits, in the way that the photon
region is confined within the range $r_{\mathrm{ph}-}\leq r\leq
r_{\mathrm{ph}+}$. The photon region is indeed determined by the condition
$\Theta(\theta)\geq 0$ for spherical photon orbits, which is specified by
applying the critical impact parameters (60) and (61) in the angular potential
(39b). In Fig. 8, the photon region around the KRBH has been demonstrated for
three different spin parameters. Note that, photons with no initial angular
momentum are able to move through the whole polar angle range, and hence,
travel on polar orbits. The radius of such orbits can be determined by the
equation $\xi^{c}=0$, which by means of Eq. (60) results in the quartic
equation
$N_{s}r^{4}-2r^{3}+3r^{2}\left(2M-a^{2}N_{s}\right)-2a^{2}r-2a^{2}M=0$, whose
the only positive solution in the domain $r\in[r_{+},r_{++}]$ can be presented
regularly as $r_{\mathrm{po}}$ (see appendix A). One should, however, note
that the presence of plasma has no effect on the radial position of polar
orbits.

(a)  (b)  (c)

Figure 8: The photon region around the KRBH plotted for three different spin
parameters, considering $N_{s}=0.1$ and $f_{r}=f_{\theta}=1$. In these
diagrams, the blue dashed curves indicate $r_{\mathrm{ph}-}$ and
$r_{\mathrm{ph}+}$, the white dashed curve is the radius of polar orbits
$r_{\mathrm{po}}$, and the black disk shows the region $r\leq r_{+}$.

####  IV.2.2 Deflecting trajectories

The radial potential (57) can be recast as

$\mathcal{R}(r)=\mathcal{p}(r)\big{[}\omega_{0}-V_{-}(r)\big{]}\big{[}\omega_{0}-V_{+}(r)\big{]},$
(64)

in which $\mathcal{p}(r)=(r^{2}+a^{2})^{2}-a^{2}\,\Delta$, and

$V_{\mp}(r)=\frac{1}{\mathcal{p}(r)}\left[a^{3}L+aL\left(r^{2}-\Delta\right)\mp\sqrt{({\mathscr{Q}+f_{r}}){\Delta}\left[r^{4}\left(1+\frac{L^{2}}{\mathscr{Q}+f_{r}}\right)+a^{2}\left(2r^{2}-\Delta\right)+a^{4}\right]}\,\right],$
(65)

indicate the radial gravitational potentials. Since the negative branch has no
classical interpretations, we rule it out from consideration and assign
$V_{\mathrm{eff}}=V_{+}(r)\equiv V(r)$, as the effective potential for photons
affected by the gravitational field of the KRBH. In Fig. 9, the above relation
has been used to demonstrate this effective potential for certain values for
the constants of motion and the plasma.

(a)  (b)

Figure 9: The radial profiles of the effective potential plotted for
$\mathscr{Q}=9$, $f_{r}=1$, and $N_{s}=0.05$. The diagrams show (a) the
profiles for three different angular momenta and $a=0.85$, where the
retrograde profile is obtained by letting $a=-0.85$, and (b) three different
profiles for $L=1$ and different spin parameters.

As it can be observed from the diagrams, the effective potential possesses no
minimums, and hence, no stable photon orbits are expected. Furthermore,
effective potentials for retrograde orbits (i.e. for $a<0$) have lower
maximums in their energy levels, compared to those for prograde orbits. Hence,
photons traveling in the opposite direction of the black hole’s spin,
experience a remarkably smoother fall in their energies.

Now, in order to study the possible photon orbits at the vicinity of the KRBH,
let us consider the typical effective potential as demonstrated in Fig. 10.
The orbits are then categorized according to the position of the initial
frequency $\omega_{0}$ with respect to $\omega_{c}$, which is the frequency of
photons on critical orbits (CO).

Figure 10: The radial profile of a typical effective potential, plotted for
$a=0.8$, $N_{s}=0.1$, $L=1$, $\mathscr{Q}=9$, and $f_{r}=1$. This way, the
radius of critical orbits is $r_{c}=2.86$, which corresponds to the critical
frequency $\omega_{c}=0.41$, or the critical impact parameter $\xi^{c}=2.437$,
in accordance with Eq. (60). The deflecting trajectories occur when
$\omega_{0}<\omega_{c}$ (or $\xi>\xi^{c}$), at $r_{d}$ or $r_{f}$. The case of
$\omega_{0}>\omega_{c}$ (or $\xi<\xi^{c}$) corresponds to the capture zone.

The analytical study of the trajectories, however, requires direct analysis of
the equation of motion (40), and hence, the radial potential $\mathcal{R}(r)$.
In fact, using the form given in Eq. (57), one can recast the radial equation
of motion as

$M{\mathrm{d}r\over\mathrm{d}\gamma}=\omega_{0}\sqrt{\mathcal{P}(r)},$ (66)

given in terms of the characteristic polynomial

$\mathcal{P}(r)=r^{4}+\sum_{j=0}^{3}\mathrm{a}_{j}r^{j},$ (67)

where

$\displaystyle\mathrm{a}_{0}=a^{4}-2a^{3}\xi-a^{2}\left[\left(\xi-a\right)^{2}+\eta_{r}\right]+a^{2}\xi^{2},$
(68a)
$\displaystyle\mathrm{a}_{1}=2M\left[\left(\xi-a\right)^{2}+\eta_{r}\right],$
(68b) $\displaystyle\mathrm{a}_{2}=2a^{2}-2a\xi-(\xi-a)^{2}-\eta_{r},$ (68c)
$\displaystyle\mathrm{a}_{3}=N_{s}\left[\left(\xi-a\right)^{2}+\eta_{r}\right].$
(68d)

As is clear from Fig. 10, deflecting trajectories correspond to
$\omega_{0}<\omega_{c}$ (or $\xi>\xi^{c}$), for which the characteristic
polynomial possesses four real roots $r_{d}>r_{+}$, $r_{+}<r_{f}<r_{d}$,
$0<r_{3}<r_{-}$ and $r_{4}<0$ (see appendix A). In fact, for the
aforementioned initial frequency (or impact parameter), Photons that approach
from $r\geq r_{d}$ will deflect from the black hole, escape towards a distant
observer, and are detected. This kind of deflecting trajectory remains
dominant up to the energy level $\omega_{0}\lesssim\omega_{c}$ and the initial
radius $r_{d}\gtrsim r_{c}$, and is known as the deflection of the first kind
(DFK). Based on the fact that the solutions to the characteristic polynomial
respect the condition $r_{4}<r_{3}<r_{f}<r_{d}<r_{\mathrm{s}}<r_{\mathrm{o}}$,
the solution to the radial equation of motion (66) for the DFK, after the
integration and performing the proper inversion in accordance with the
observer at $r_{\mathrm{o}}$, is obtained as

$r_{\mathrm{o}}(\gamma)={r_{d}\left(r_{f}-r_{4}\right)-r_{f}\left(r_{d}-r_{4}\right)\operatorname{\mathbf{sn}}^{2}\Big{(}\tilde{\kappa}\gamma+\operatorname{\mathbf{F}}(\tilde{\varphi}_{\mathrm{s}}|\tilde{\mathcal{k}})\Big{|}\tilde{\mathcal{k}}\Big{)}\over\left(r_{f}-r_{4}\right)-\left(r_{d}-r_{4}\right)\operatorname{\mathbf{sn}}^{2}\Big{(}\tilde{\kappa}\gamma+\operatorname{\mathbf{F}}(\tilde{\varphi}_{\mathrm{s}}|\tilde{\mathcal{k}})\Big{|}\tilde{\mathcal{k}}\Big{)}},$
(69)

in which

$\displaystyle\tilde{\kappa}=\frac{\omega_{0}\sqrt{(r_{d}-r_{3})(r_{f}-r_{4})}}{2M},$
(70a)
$\displaystyle\tilde{\varphi}_{\mathrm{s}}=\arcsin\left(\sqrt{\frac{(r_{f}-r_{4})(r_{\mathrm{s}}-r_{d})}{(r_{d}-r_{4})(r_{\mathrm{s}}-r_{f})}}\right),$
(70b)
$\displaystyle\tilde{\mathcal{k}}=\frac{(r_{f}-r_{3})(r_{d}-r_{4})}{(r_{d}-r_{3})(r_{f}-r_{4})}.$
(70c)

On the other hand, photons with the same initial frequencies experience a
plunge into the black hole, if they approach the black hole from the radial
distances $r\leq r_{f}$. This way, they deflect into the event horizon and are
never detected by a distant observer. This kind of deflection is dominant in
the same energy levels as in the DFK and is radially valid for the initial
radius $r_{f}\lesssim r_{c}$, and is known as the deflection of the second
kind (DSK). In this case, the inequality
$r_{4}<r_{3}<r_{\mathrm{o}}<r_{f}<r_{d}<r_{\mathrm{s}}$ holds, and hence, the
evolution of the $r$-coordinate for the DSK is obtained as

$r_{\mathrm{o}}(\gamma)={r_{f}\left(r_{d}-r_{3}\right)-r_{d}\left(r_{f}-r_{3}\right)\operatorname{\mathbf{sn}}^{2}\Big{(}\tilde{\kappa}\gamma-\operatorname{\mathbf{F}}(\stackon[-1.3ex]{\stackon[-1ex]{\varphi}{\mathchar
12414\relax}}{\mathchar
12414\relax}_{\mathrm{s}}|\tilde{\mathcal{k}})\Big{|}\tilde{\mathcal{k}}\Big{)}\over\left(r_{d}-r_{3}\right)-\left(r_{f}-r_{3}\right)\operatorname{\mathbf{sn}}^{2}\Big{(}\tilde{\kappa}\gamma-\operatorname{\mathbf{F}}(\stackon[-1.3ex]{\stackon[-1ex]{\varphi}{\mathchar
12414\relax}}{\mathchar
12414\relax}_{\mathrm{s}}|\tilde{\mathcal{k}})\Big{|}\tilde{\mathcal{k}}\Big{)}},$
(71)

in which

$\stackon[-1.3ex]{\stackon[-1ex]{\varphi}{\mathchar 12414\relax}}{\mathchar
12414\relax}_{\mathrm{s}}=\arcsin\left(\sqrt{\frac{(r_{d}-r_{3})(r_{f}-r_{\mathrm{s}})}{(r_{f}-r_{3})(r_{d}-r_{s})}}\right).$
(72)

In Fig. 11, an example of the DFK and DSK is presented in the polar plane,
aligning with the data provided in Fig. 10.

(a)  (b)

Figure 11: Examples of (a) DFK, and (b) DSK, inside the cone of
$\theta$-confinement (confined between $\theta_{1}$ and $\theta_{4}$ in
accordance with Eqs. (51) and (54)), plotted for $a=0.8$, $N_{s}=0.1$,
$\mathscr{Q}=9$, $f_{r}=f_{\theta}=1$, $L=1$, $\omega_{0}=0.3$ (corresponding
to $\tilde{\eta}=88.89$, $\eta_{r}=111.11$ and $\xi=3.33$), $\theta_{s}=\pi/2$
and $r_{s}=5$. The minimum and maximum polar angles are
$\theta_{1}=19.90^{\circ}$ and $\theta_{4}=160.10^{\circ}$ from the $Z$-axis.
The red, blue, and green dashed circles correspond, respectively, to $r_{+}$,
$r_{d}$ and $r_{f}$.

####  IV.2.3 Critical trajectories

As stated above, photons on the orbits of constant radius, are indeed
traveling on CO with the initial frequency $\omega_{0}=\omega_{c}$. According
to the effective potential in Fig. 10, such orbits are indeed unstable.
However, based on the point from which the photons approach the black hole,
this instability may result in different fates; the photons on CO may finally
plunge into the black hole or escape from it. Hence, the critical orbit of the
first kind (COFK) occurs when the photons approach the black hole from $r\geq
r_{c}$, and after performing spherical orbits for a certain duration, they
finally deflect and escape from the black hole. On the other hand, photons
perform a critical orbit of the second kind (COSK), when they approach from
$r\leq r_{c}$, which finally leads to a plunge into the black hole after
performing spherical orbits. To obtain the analytical solutions for the
evolution of the $r$-coordinate for such orbits, let us highlight the fact
that the characteristic polynomial (67) in this case can be recast as
$\mathcal{P}(r)=(r-r_{c})^{2}(r-r_{3})(r-r_{4})$. This way, direct integration
from the radial equation of motion yields

$r_{\mathrm{o}}(\gamma)=r_{c}+{2\left(r_{c}-r_{3}\right)\left(r_{c}-r_{4}\right)\over\left(r_{3}-r_{4}\right)\cosh{\left(\tilde{\kappa}_{c}\gamma+\tilde{\varphi}_{c}\right)}-\left(2r_{c}-r_{3}-r_{4}\right)},$
(73)

for the COFk, and

$r_{\mathrm{o}}(\gamma)=r_{c}-{2\left(r_{c}-r_{3}\right)\left(r_{c}-r_{4}\right)\over\left(r_{3}-r_{4}\right)\cosh{\left(\tilde{\kappa}_{c}\gamma\right)}+\left(2r_{c}-r_{3}-r_{4}\right)},$
(74)

for the COSK, where

$\displaystyle\tilde{\kappa}_{c}=\frac{\omega_{c}\sqrt{\left(r_{c}-r_{3}\right)\left(r_{c}-r_{4}\right)}}{M},$
(75a)
$\displaystyle\cosh{(\tilde{\varphi}_{c})}=\frac{2r_{c}-r_{3}-r_{4}}{r_{3}-r_{4}}.$
(75b)

In Fig. 12, an example for each of the above kinds of orbits has been
demonstrated in the polar plane, in accordance with the data given in Fig. 10

(a)  (b)

Figure 12: Examples of (a) COFK, and (b) COSK, inside the cone of confinement
for the $\theta$-coordinate ($\theta_{1}=49.55^{\circ}$ and
$\theta_{4}=130.44^{\circ}$), plotted for $a=0.8$, $N_{s}=0.1$,
$\mathscr{Q}=9$, $f_{r}=f_{\theta}=1$, $L=1$ and $\omega_{0}=\omega_{c}=0.41$
(corresponding to $\tilde{\eta}=47.52$, $\eta_{r}=59.40$ and
$\xi=2.44=\xi^{c}$). The red and brown dashed circles correspond,
respectively, to $r_{+}$ and $r_{c}=2.86$.

####  IV.2.4 The capture zone

In the case that $\omega_{0}>\omega_{c}$ (or $\xi<\xi^{c}$), the photons do
not encounter any turning points, and hence, become completely unstable in the
exterior region of the black hole. Such photons do travel on particular bound
orbits which finally plunge into the black hole, and form the capture zone. To
obtain the analytical solution for the evolution of the radial coordinate for
this type of motion, we should bear in mind that in this case, the
characteristic polynomial takes the form
$\mathcal{P}(r)=(r-r_{1})(r-r_{1}^{*})(r-r_{3})(r-r_{4})$, in which
$r_{1}\in\mathbb{C}$ and is obtained by means of the relation $r_{1}=r_{d}$,
together with the extra condition $r_{1}^{*}=r_{f}$, in accordance with the
values given in Eqs. (115) and (116). According to the fact that in this case
the inequalities $r_{4}<r_{3}<r_{\mathrm{o}}<r_{\mathrm{s}}$ hold, the direct
integration of Eq. (66) results in

$r_{\mathrm{o}}(\gamma)=\frac{\tilde{B}r_{3}-\tilde{A}r_{4}+(\tilde{B}r_{3}+\tilde{A}r_{4})\operatorname{\mathbf{cn}}\Big{(}\tilde{\kappa}_{p}\gamma+\operatorname{\mathbf{F}}(\tilde{\varphi}_{p,\mathrm{s}}|\tilde{\mathcal{k}}_{p})\Big{|}\tilde{\mathcal{k}}_{p}\Big{)}}{(\tilde{B}-\tilde{A})+(\tilde{A}+\tilde{B})\operatorname{\mathbf{cn}}\Big{(}\tilde{\kappa}_{p}\gamma+\operatorname{\mathbf{F}}(\tilde{\varphi}_{p,\mathrm{s}}|\tilde{\mathcal{k}}_{p})\Big{|}\tilde{\mathcal{k}}_{p}\Big{)}},$
(76)

in which $\operatorname{\mathbf{cn}}(x|\tilde{\mathcal{k}}_{p})$ is the
Jacobian elliptic cosine function with argument $x$ byrd_handbook_1971 , and

$\displaystyle\tilde{A}=\sqrt{\left(r_{d}-r_{3}\right)\left(r_{f}-r_{3}\right)},$
(77a)
$\displaystyle\tilde{B}=\sqrt{\left(r_{d}-r_{4}\right)\left(r_{f}-r_{4}\right)},$
(77b)
$\displaystyle\tilde{\kappa}_{p}=\frac{\omega_{0}\sqrt{\tilde{A}\tilde{B}}}{M},$
(77c)
$\displaystyle\tilde{\mathcal{k}}_{p}=\frac{(\tilde{A}+\tilde{B})^{2}-\left(r_{3}-r_{4}\right)^{2}}{4\tilde{A}\tilde{B}},$
(77d)
$\displaystyle\tilde{\varphi}_{p,\mathrm{s}}=\arccos\left(\frac{r_{3}\tilde{B}-r_{4}\tilde{A}+(\tilde{A}-\tilde{B})r_{s}}{-r_{3}\tilde{B}-r_{4}\tilde{A}+(\tilde{A}-\tilde{B})r_{s}}\right).$
(77e)

In Fig. 13, two examples of this kind of orbit have been plotted, which
correspond to the capture zone of the KRBH.

Figure 13: two examples of captured orbits for $a=0.8$, $N_{s}=0.1$,
$\mathscr{Q}=9$, $f_{r}=f_{\theta}=1$, $L=1$,
$\omega_{0}=\omega_{c}+0.3=0.71$, $\theta_{s}=\pi/2$, and the initial radii
$r_{s}=5,10$. These orbits are confined inside the angular domain
$\theta_{1}=19.15$ and $\theta_{4}=160.85$. The red dashed circle corresponds
to $r_{+}$.

###  IV.3 Azimuth motion

In order to obtain the evolution of the $\phi$-coordinate, let us divide the
integral of motion (42) into the two parts

$\phi(\gamma)=\Phi_{\theta}(\gamma)+\Phi_{r}(\gamma),$ (78)

in which

$\displaystyle\Phi_{\theta}(\gamma)=\frac{L}{M}\int^{\pi/2}_{\theta_{i}(\gamma)}\frac{\mathrm{d}\theta}{\sin^{2}\theta\sqrt{\Theta(\theta)}},$
(79a)
$\displaystyle\Phi_{r}(\gamma)={\frac{1}{M}}\int^{r(\gamma)}_{r_{i}}\frac{\left(2Ma\omega_{0}r+N_{s}a\omega_{0}r^{3}-a^{2}L\right)\mathrm{d}r}{\Delta\sqrt{\mathcal{R}(r)}}.$
(79b)

Without loss of generality, we have taken the minimum value of the
$\theta$-coordinate to coincide with the initial point $r_{i}$, which can be
set as either of the turning points, in accordance with the considered type of
orbit. For $\tilde{\eta}>0$, direct integration of the integral (79a) results
in

$\Phi_{\theta}(\gamma)=\frac{L}{Ma\sqrt{-z_{-}}}\Big{[}\operatorname{\mathbf{\Pi}}(z_{+};\varphi_{\mathrm{s}}|\mathcal{k})-\operatorname{\mathbf{\Pi}}(z_{+};\varphi_{\mathrm{o}}(\gamma)|\mathcal{k})\Big{]},$
(80)

with

$\varphi_{\mathrm{o}}(\gamma)=\arcsin\left(\frac{\cos\left(\theta_{\mathrm{o}}(\gamma)\right)}{\sqrt{z_{+}}}\right),$
(81)

in which $\theta_{\mathrm{o}}(\gamma)$ has been derived in Eq. (55), and all
the other included parameters are the same as those in that equation.

Now, the $r$-dependent integral yields the following categories, corresponding
to the different possible types of orbits:

  * (i)

DFK: In this case, after proper manipulation of the integral, we get to the
solution

$\Phi_{r}(\gamma)=-\frac{a\omega_{0}}{M}g_{d}\left[\operatorname{\mathbf{F}}(\tilde{\varphi}_{\mathrm{o}}(\gamma)|\tilde{\mathcal{k}})-\operatorname{\mathbf{F}}(\tilde{\varphi}_{\mathrm{s}}|\tilde{\mathcal{k}})\right]+\frac{g_{d}}{MN_{s}}\sum_{j=1}^{3}A_{j}^{\mathrm{DFK}}I_{j}^{\mathrm{DFK}}(\gamma),$
(82)

where

$\displaystyle g_{d}=\frac{2}{\sqrt{(r_{d}-r_{3})(r_{f}-r_{4})}},$ (83a)
$\displaystyle\tilde{\varphi}_{\mathrm{o}}(\gamma)=\arcsin\left(\sqrt{\frac{(r_{f}-r_{4})(r_{\mathrm{o}}(\gamma)-r_{d})}{(r_{d}-r_{4})(r_{\mathrm{o}}(\gamma)-r_{f})}}\right),$
(83b)

in which $r_{\mathrm{o}}(\gamma)$ is identified in Eq. (69) for the DFK,
$\tilde{\varphi}_{\mathrm{s}}$ and $\tilde{\mathcal{k}}$ are given in Eqs.
(70b) and (70c), and the coefficients $A_{j}^{\mathrm{DFK}}$ and the elliptic
integrals $I_{j}^{\mathrm{DFK}}$ are presented in appendix B. With these
values considered, Fig. 14 illustrates two examples of the DFK in the context
of the effective potential shown in Fig. 10, assuming a thin disk in the black
hole’s equatorial plane. These diagrams highlight two significant scenarios.
Firstly, direct emission occurs when the initial frequency differs
significantly from $\omega_{c}$ (panel (a)). In such cases, photons emitted
from the disk directly deflect at $r_{d}$ and escape the black hole without
intersecting the disk. Conversely, as $\omega_{0}$ approaches $\omega_{c}$
(panel (b)), the emitted photons approach the turning point $r_{d}$ and then
curve away from the black hole after crossing the disk. As a result, observers
perceive a demagnified lensed image of the disk. This phenomenon generates an
image of higher order compared to that from direct emission. The diagrams also
display both initial frequency scenarios, but for the case where the KRBH
exists in a vacuum. In this instance, direct emission exhibits more variation
in the $\theta$-coordinate, whereas the additional lensed image yields a
demagnified representation of the disk’s backside.

(a)  (b)

Figure 14: Examples of the DFK plotted for the specific parameter values
$a=0.8$, $N_{s}=0.1$, $\mathscr{Q}=9$, $L=1$, initial radii $r_{\mathrm{s}}=4$
(indicated as the yellow point), $\theta_{s}=\pi/2$, and two distinct initial
frequencies, namely, (a) $\omega_{0}=0.3$ and (b) $\omega_{0}=0.375$. The
black sphere in the diagrams represents the black hole with a radius of
$r_{+}$, and the orange region indicates the thin emission disk situated in
the equatorial plane, illuminating the black hole. The red and blue curves in
the diagrams correspond, respectively, to the paths of light rays propagating
around the KRBH in a plasmic medium with $f_{r}=f_{\theta}=1$, and in a vacuum
with $f_{r}=f_{\theta}=0$. In diagram (a), the DFK for both cases results in
simple lensed images produced by the direct emission from the disk. In diagram
(b), both scenarios involve light rays crossing the disk before escaping the
black hole. This leads to the observer detecting demagnified lensed images of
points on the disk, stemming from the rays completing half orbits around the
black hole before their ultimate deflection. In the scenario involving the
KRBH surrounded by a plasmic medium, this mechanism yields higher-order images
originating from the front side of the disk. Conversely, for the KRBH in a
vacuum, these images correspond to the backside of the disk.

In fact, by having at hand the analytical solutions for the azimuth motion in
the course of the DFK, one can also calculate the deflection angle, as
detected by an observer at a finite radial distance $r_{\mathrm{o}}$.
Confining ourselves to the equatorial plane and exploiting Eqs. (37) and (40),
one yields the equation

$\phi_{i}=\int_{r_{d}}^{r_{i}}\frac{\left(2Ma\omega_{0}r+N_{s}a\omega_{0}r^{3}-a^{2}L\right)\mathrm{d}r}{\Delta\sqrt{\mathcal{R}(r)}}-L\int_{r_{d}}^{r_{i}}\frac{\mathrm{d}r}{\sqrt{\mathcal{R}(r)}},$
(84)

which provides the extent of azimuth angle covered as the light traverses the
distance between the turning point $r_{d}$ and the fixed radial position
$r_{i}$. Accordingly, the lens equation can be written as
fathi_analytical_2021 ; fathi_study_2022

$\hat{\vartheta}=(\phi_{\mathrm{o}}+\phi_{\mathrm{s}})-\pi.$ (85)

This relation calculates the deflection angle $\hat{\vartheta}$ of the light
rays around the KRBH, where

$\displaystyle\phi_{\mathrm{o}}=\left.-(a\omega_{0}+L)g_{d}\operatorname{\mathbf{F}}(\tilde{\varphi}_{i}|\tilde{\mathcal{k}})\right|_{r_{d}}^{r_{\mathrm{o}}}+\frac{g_{d}}{N_{s}}\sum_{j=1}^{3}A_{j}^{\mathrm{DFK}}I_{j}^{\mathrm{L}_{\mathrm{o}}},$
(86a)
$\displaystyle\phi_{\mathrm{s}}=\left.-(a\omega_{0}+L)g_{d}\operatorname{\mathbf{F}}(\tilde{\varphi}_{i}|\tilde{\mathcal{k}})\right|_{r_{d}}^{r_{\mathrm{s}}}+\frac{g_{d}}{N_{s}}\sum_{j=1}^{3}A_{j}^{\mathrm{DFK}}I_{j}^{\mathrm{L}_{\mathrm{s}}},$
(86b)

in which $I_{j}^{\mathrm{L}_{\mathrm{o}}}$ and
$I_{j}^{\mathrm{L}_{\mathrm{s}}}$ are obtained, respectively, by doing the
pair of exchanges $\left\\{r_{\mathrm{s}}\rightarrow
r_{d},r_{\mathrm{o}}(\gamma)\rightarrow
r_{\mathrm{o}}=\mathrm{const.}\right\\}$, and
$\left\\{r_{\mathrm{s}}\rightarrow r_{d},r_{\mathrm{o}}(\gamma)\rightarrow
r_{\mathrm{s}}=\mathrm{const.}\right\\}$ in $I^{\mathrm{DFK}}_{j}$. In Fig.
15, the deflection angle profile $\hat{\vartheta}$ has been plotted against
the impact parameter $\xi$ for an illustrative case. As observed in the
figure, the deflection angle diverges at the critical impact parameter
$\xi_{c}$, and subsequently decreases rapidly until it reaches a minimum
value. It then gradually increases until peaking at a specific impact
parameter, followed by a rapid decline once more. From these diagrams, it’s
evident that when light propagates around the KRBH in a vacuum, the deflection
angle exhibits a lower minimum but a higher maximum compared to its behavior
within a plasmic medium. Generally, the deflection angle values in a plasmic
medium tend to be larger than those in a vacuum, an anticipated outcome due to
the nonunity refractive index of the plasma. Furthermore, gravitational
lensing by the KRBH spans a broader range of impact parameters when the black
hole exists in a vacuum as opposed to being immersed in plasma. This concept
will be explored further in an upcoming study, delving into the extensive
gravitational lensing by the KRBH within a plasmic medium.

Figure 15: The $\xi$-profile of the deflection angle $\hat{\vartheta}$ plotted
for $a=0.8$, $N_{s}=0.1$, $\mathscr{Q}=9$, $r_{\mathrm{s}}=5$ and
$r_{\mathrm{o}}=100$. The solid curve corresponds to the plasmic medium with
$f_{r}=1$, whereas the dashed curve corresponds to the vacuum case.

  * (ii)

DSK: Based on the same methods, the direct integration of Eq. (79b) for the
case of the DSK, results in

$\Phi_{r}(\gamma)=-\frac{a\omega_{0}}{M}g_{d}\Big{[}\operatorname{\mathbf{F}}(\stackon[-1.3ex]{\stackon[-1ex]{\varphi}{\mathchar
12414\relax}}{\mathchar
12414\relax}_{\mathrm{o}}(\gamma)|\tilde{\mathcal{k}})-\operatorname{\mathbf{F}}(\stackon[-1.3ex]{\stackon[-1ex]{\varphi}{\mathchar
12414\relax}}{\mathchar
12414\relax}_{\mathrm{s}}|\tilde{\mathcal{k}})\Big{]}+\frac{g_{d}}{MN_{s}}\sum_{j=0}^{3}A^{\mathrm{DSK}}_{j}I_{j}^{\mathrm{DSK}},$
(87)

for which $A_{j}^{\mathrm{DSK}}=A_{j}^{\mathrm{DFK}}$ as given in appendix B,
$\stackon[-1.3ex]{\stackon[-1ex]{\varphi}{\mathchar 12414\relax}}{\mathchar
12414\relax}_{\mathrm{s}}$ is defined in Eq. (72), and

$\stackon[-1.3ex]{\stackon[-1ex]{\varphi}{\mathchar 12414\relax}}{\mathchar
12414\relax}_{\mathrm{o}}(\gamma)=\arcsin\left(\sqrt{\frac{(r_{d}-r_{3})(r_{f}-r_{\mathrm{o}}(\gamma))}{(r_{f}-r_{3})(r_{d}-r_{\mathrm{o}}(\gamma))}}\right),$
(88)

where $r_{\mathrm{o}}(\gamma)$ has been expressed in Eq. (71). The relevant
expressions for the integrals $I_{j}^{\mathrm{DSK}}$ have been given in Eqs.
(127) and (128) in appendix B. Based on the fact that the DSK occurs in the
vicinity of the event horizon, it does not have any particular latitudinal
behavior, as illustrated in Fig. 11(b). We have therefore omitted any diagrams
representing this type of orbit.

  * (iii)

COFK: Based on the discussion provided in Subsect. IV.2.3, the integral (79b)
for this type of orbit results in

$\Phi_{r}(\gamma)=\frac{2}{MN_{s}}\sum_{j=1}^{4}A_{j}^{\mathrm{COFK}}I_{j}^{\mathrm{COFK}}(\gamma),$
(89)

for which, the coefficients $A_{j}^{\mathrm{COFK}}$ and the hyperbolic
expressions $I_{j}^{\mathrm{COFK}}(\gamma)$ have been given in Eqs. (129) and
(130). Some examples of such orbits are depicted in Fig. 16. In each scenario,
the photon emission initiates at the radial distance $r_{c}$, where the
photons execute unstable spherical orbits before escaping the black hole. As
observed from the diagrams, a higher initial angular momentum results in a
narrower range of oscillations of the $\theta$-coordinate. Additionally, as a
general trend, the number of times the orbits traverse the equatorial plane
(i.e., the number of nodes) is greater in a vacuum environment compared to
that within a plasmic medium.

(a)  (b)  (c)  (d)

Figure 16: Examples of the COFK for specific parameter values: $a=0.8$,
$N_{s}=0.1$, $\mathscr{Q}=9$, and $L=1$. The red and blue curves represent
light propagation in plasma with $f_{r}=f_{\theta}=1$ and in a vacuum,
respectively. The diagrams showcase the following cases: (a) $L=1$:
$r_{c}=2.86$ and $\omega_{c}=0.41$ in plasma, and $r_{c}=2.85$ and
$\omega_{c}=0.39$ in vacuum, (b) $L=5$: $r_{c}=2.51$, $\omega_{c}=1.02$ in
plasma, and $r_{c}=2.50$, $\omega_{c}=1.01$ in vacuum, (c) $L=10$:
$r_{c}=2.458$, $\omega_{c}=1.914$ in plasma, and $r_{c}=2.456$,
$\omega_{c}=1.909$ in vacuum, (d) $L=0$: $r_{c}=3.12$, $\omega_{c}=0.32$ in
plasma, and $r_{c}=3.12$, $\omega_{c}=0.30$ in vacuum. In each case, photons
start from the peak of the effective potential, perform spherical orbits at
the radial distance $r_{c}$, and eventually escape the black hole.

Due to the similarity in the evolution of the radial coordinates between COFK
and COSK, as seen in Eqs. (73) and (74), their analytical azimuth motion
solutions show minimal distinctions. Consequently, the shapes of the orbits
remain fundamentally unchanged. As a result, we will omit this scenario and
move forward with the plunge orbit analysis

  * (iv)

The capture zone: For plunge orbits discussed in Subsect. IV.2.4, the integral
equation (79b) provides the analytical solution

$\Phi_{r}(\gamma)=-\frac{a\omega_{0}}{M}\Big{[}\operatorname{\mathbf{F}}(\tilde{\varphi}_{p,\mathrm{o}}|\tilde{\mathcal{k}}_{p})-\operatorname{\mathbf{F}}(\tilde{\varphi}_{p,\mathrm{s}}|\tilde{\mathcal{k}}_{p})\Big{]}+\frac{1}{MN_{s}}\sum_{j=1}^{3}A_{j}^{\mathrm{PO}}I_{j}^{\mathrm{PO}}(\gamma),$
(90)

for which the coefficients $A_{j}^{\mathrm{PO}}$ and the elliptic integrals
$I_{j}^{\mathrm{PO}}$ have been given in Eqs. (131)–(137). Some examples of
this kind of orbit have been demonstrated in Fig. 17.

(a)  (b)

Figure 17: Two examples of plunge orbits, plotted by considering $a=0.8$,
$N_{s}=0.1$, $\mathscr{Q}=9$, $L=1$ and $\theta_{s}=\pi/2$. The red and blue
curves correspond respectively to the propagation of light in plasma with
$f_{r}=f_{\theta}=1$ and $\omega_{0}=\omega_{c}+0.02=0.43$, and in vacuum with
$\omega_{0}=0.41$. The diagrams correspond to (a) $r_{s}=5$, and (b)
$r_{s}=6.5$.

As evident from the diagrams, in both scenarios of light propagation through
plasma and in a vacuum, the light rays initiate their trajectory towards the
black hole from the rear side of the emission disk. This occurrence is solely
a result of the spacetime geometry, which dictates this behavior in the rays
for plunge orbits (see also Fig. 13).

###  IV.4 The $t$-coordinate evolution

By employing the same methodology utilized in computing integrals for the
$\phi$-motion, and considering Eq. (43) alongside Eqs. (40) and (41), the
integral governing the evolution of the $t$-coordinate can be isolated as

$t(\gamma)=t_{\theta}(\gamma)+t_{r}(\gamma),$ (91)

with

$\displaystyle
t_{\theta}(\gamma)=-\omega_{0}a^{2}\int^{\pi/2}_{\theta_{i}(\gamma)}\frac{\sin^{2}\theta\,\mathrm{d}\theta}{\sqrt{\Theta(\theta)}},$
(92a) $\displaystyle
t_{r}(\gamma)=\frac{1}{M}\int^{r(\gamma)}_{r_{i}}\frac{\left[\omega_{0}(r^{2}+a^{2})^{2}-2MaLr-
N_{s}aLr^{3}\right]\mathrm{d}r}{\Delta\sqrt{\mathcal{R}(r)}}.$ (92b)

Direct integration of the angular part in Eq. (92a), then results in the
analytical solution

$t_{\theta}(\gamma)=-\frac{a}{\sqrt{-z_{-}}}\left[\left(1-{z_{-}\over\mathcal{k}}\right)\Big{(}\operatorname{\mathbf{F}}(\varphi_{\mathrm{s}}|\mathcal{k})-\operatorname{\mathbf{F}}(\varphi_{\mathrm{o}}(\gamma)|\mathcal{k})\Big{)}+{z_{-}\over\mathcal{k}}\Big{(}\operatorname{\mathbf{E}}(\varphi_{\mathrm{s}}|\mathcal{k})-\operatorname{\mathbf{E}}(\varphi_{\mathrm{o}}(\gamma)|\mathcal{k})\Big{)}\right],$
(93)

in which $\varphi_{\mathrm{o}}(\gamma)$ has been given in Eq. (81), and
$\operatorname{\mathbf{E}}(x|\mathcal{k})$ is the elliptic integral of the
second kind with argument $x$ and modulus $\mathcal{k}$ byrd_handbook_1971 .
Moreover, when considering the most general case for the radial integral
(92b), which corresponds to the DFK, the solution is given by

$t_{r}(\gamma)=-\frac{aLN_{s}+\omega_{0}(r_{++}+r_{+}+r_{-})}{MN_{s}}g_{d}\left[\operatorname{\mathbf{F}}(\tilde{\varphi}_{\mathrm{o}}(\gamma)|\tilde{\mathcal{k}})-\operatorname{\mathbf{F}}(\tilde{\varphi}_{\mathrm{s}}|\tilde{\mathcal{k}})\right]+\sum_{j=1}^{3}B_{j}^{t}I^{\mathrm{DFK}}_{j}(\gamma)+\frac{\omega_{0}u_{f}\sqrt{-u_{d}u_{f}u_{3}u_{4}}}{u_{d}^{2}}g_{t}I^{t}(\gamma),$
(94)

where $\tilde{\varphi}_{\mathrm{o}}(\gamma)$ and
$\tilde{\varphi}_{\mathrm{s}}$ are provided in Eqs. (83b) and (70b),
respectively, while $I_{j}^{\mathrm{DFK}}(\gamma)$ is defined in Eqs. (125),
and the coefficients $B_{j}^{t}$ are outlined in Eqs. (138). Additionally,
we’ve introduced $u_{i}=1/r_{i}$, along with

$g_{t}=\frac{2}{\sqrt{(u_{f}-u_{d})(u_{f}-u_{4})}},$ (95)

and

$I^{t}(\gamma)=\left.-\frac{(u_{d}-u_{4})(u_{f}-u_{d})}{u_{d}(u_{f}-u_{4})}\operatorname{\mathbf{F}}(\tilde{\varphi}^{t}_{i}|\tilde{\mathcal{k}}_{t})\right|_{\tilde{\varphi}^{t}_{\mathrm{o}}(\gamma)}^{\tilde{\varphi}^{t}_{\mathrm{s}}}\left.-\frac{u_{f}-u_{d}}{u_{f}}\operatorname{\mathbf{\Pi}}\left(\frac{u_{d}-u_{4}}{u_{f}-u_{4}},\tilde{\varphi}_{i}^{t}|\tilde{\mathcal{k}}_{t}\right)\right|_{\tilde{\varphi}^{t}_{\mathrm{o}}(\gamma)}^{\tilde{\varphi}^{t}_{\mathrm{s}}},$
(96)

in which

$\tilde{\mathcal{k}}_{t}=\sqrt{\frac{(u_{3}-u_{f})(u_{d}-u_{4})}{(u_{3}-u_{d})(u_{f}-u_{4})}},$
(97)

and

$\displaystyle\tilde{\varphi}^{t}_{\mathrm{s}}=\arcsin\left(\sqrt{\frac{(u_{f}-u_{4})(u_{d}-u_{\mathrm{s}})}{(u_{d}-u_{4})(u_{f}-u_{\mathrm{s}})}}\,\right),$
(98a)
$\displaystyle\tilde{\varphi}^{t}_{\mathrm{o}}(\gamma)=\arcsin\left(\sqrt{\frac{(u_{f}-u_{4})\left(u_{d}-u_{\mathrm{o}}(\gamma)\right)}{(u_{d}-u_{4})\left(u_{f}-u_{\mathrm{o}}(\gamma)\right)}}\,\right),$
(98b)

where $u_{\mathrm{o}}(\gamma)=1/r_{\mathrm{o}}(\gamma)$, and
$r_{\mathrm{o}}(\gamma)$ is outlined in Eq. (69). It is worth noting that our
orbit simulations exclusively utilized the evolution of the three spatial
coordinates. Hence, the inclusion of the $t$-coordinate motion analysis is
only for the sake of completeness in our study.

##  V More on the spherical orbits

As highlighted in point (iii) of Subsect. IV.3, spherical orbits are pivotal
within photon trajectories navigating critical paths. These orbits notably
contribute to forming the black hole shadow, representing the closest orbits
that escape the black hole and become observable to the observer. However,
photons traversing these orbits possess two crucial characteristics requiring
careful consideration. First, photons take a finite time to complete each half
orbit around the black hole, thereby inducing a time delay for these photons
contributing to the total time for subsequent images reaching the observer.
Second, the azimuth angle experiences distinct changes during each half orbit,
influencing the rotation of the photon ring and ensuing subrings. These
parameters significantly impact the determination of higher-order images
generated by the black hole. Hence, this section focuses on calculating the
azimuth angle changes and time delays occurring during critical orbits at
$r_{c}$. Given these orbits remain at a constant radius from the black hole,
the radial part of the integrals of motion remains constant. The determination
of each parameter involves integration over half orbits, between the angular
turning points of the angular potential, $\theta_{1}$ and $\theta_{4}$, as
discussed in Subsect. IV.1.

We begin by computing the changes in the azimuth angle,
$\phi_{\mathrm{o}}-\phi_{\mathrm{s}}$, throughout the orbits. By incorporating
this constant radius into Eq. (42), we arrive at the integral

$\delta\phi_{c}\equiv\phi_{\mathrm{o}}-\phi_{\mathrm{s}}=\frac{1}{M}\left[\frac{a\omega_{0}\left(r_{c}^{2}+a^{2}-\Delta_{c}\right)-La^{2}}{\Delta_{c}}\right]\int_{\theta_{1}}^{\theta_{4}}\frac{\mathrm{d}\theta}{\sqrt{\Theta(\theta)}}+\frac{L}{M}\int_{\theta_{1}}^{\theta_{4}}\frac{\mathrm{d}\theta}{\sin^{2}\theta\sqrt{\Theta(\theta)}},$
(99)

in which, we have considered the identity
$\gamma=M\int\frac{\mathrm{d}r}{\sqrt{\mathcal{R}(r)}}=M\int\frac{\mathrm{d}\theta}{\sqrt{\Theta(\theta)}}$,
and $\Delta_{c}\equiv\Delta(r_{c})$. Direct integration of Eq. (99) results in

$\delta\phi_{c}=\frac{2}{z_{0}}\left[\frac{a\omega_{0}\left(r_{c}^{2}+a^{2}-\Delta_{c}\right)-La^{2}}{\Delta_{c}}\right]\operatorname{\mathbf{K}}(\mathcal{k})+\frac{2L\omega_{0}}{Mz_{0}}\operatorname{\mathbf{\Pi}}(z_{+}|\mathcal{k}),$
(100)

where
$\operatorname{\mathbf{K}}(\mathcal{k})=\operatorname{\mathbf{F}}(\pi/2|\mathcal{k})$
is the complete elliptic integral of the first kind, and
$\operatorname{\mathbf{\Pi}}(z_{+}|\mathcal{k})=\operatorname{\mathbf{\Pi}}(z_{+};\pi/2|\mathcal{k})$
is the complete elliptic integral of the third kind byrd_handbook_1971 .
Similarly, to determine the time delay during these orbits, we analyze Eq.
(43) for $r=r_{c}$, which yields the relationship

$\delta t_{c}\equiv
t_{\mathrm{o}}-t_{\mathrm{s}}=\frac{1}{M}\left[\frac{\omega_{0}\left(r_{c}^{2}+a^{2}\right)^{2}-a\left(r_{c}^{2}+a^{2}-\Delta_{c}\right)}{\Delta_{c}}\right]\int_{\theta_{1}}^{\theta_{4}}\frac{\mathrm{d}\theta}{\sqrt{\Theta(\theta)}}-a^{2}\omega_{0}\int_{\theta_{1}}^{\theta_{4}}\frac{\sin^{2}\theta\,\mathrm{d}\theta}{\sqrt{\Theta(\theta)}},$
(101)

and has the solution

$\delta
t_{c}=\frac{2}{z_{0}}\left[\frac{\omega_{0}\left(r_{c}^{2}+a^{2}\right)^{2}-a\left(r_{c}^{2}+a^{2}-\Delta_{c}\right)}{\Delta_{c}}\right]\operatorname{\mathbf{K}}(\mathcal{k})-\frac{2a^{2}\omega_{0}^{2}}{z_{0}}\left[\operatorname{\mathbf{K}}(\mathcal{k})-\frac{4z_{+}}{\mathcal{k}}\operatorname{\mathbf{{J}}}(\mathcal{k})\right],$
(102)

where we have defined

$\operatorname{\mathbf{{J}}}(z|\mathcal{k})=\operatorname{\mathbf{F}}(z|\mathcal{k})-\operatorname{\mathbf{E}}(z|\mathcal{k}).$
(103)

Hence,
$\operatorname{\mathbf{{J}}}(\mathcal{k})=\operatorname{\mathbf{{J}}}(\pi/2|\mathcal{k})$
represents the complete form of this function. These findings hold promise for
computing alterations in crucial parameters during light propagation around
black holes. However, for the KRBH in a plasmic medium, the application of
these outcomes to examine photon rings and subrings as images of the black
hole’s accretion disk remains a task for future research. Consequently, we
conclude our current discussion here and provide a summary of our results in
the next section.

##  VI Summary and conclusions

This study delved into the behavior of light near a rotating black hole within
Rastall gravity, considering a structure filled with cold, non-magnetized
electronic plasma. Our investigation commenced with an overview of Rastall
gravity, introducing its static spherically symmetric black hole solution
involving an anisotropic cosmic fluid. We extended this solution to the
rotating scenario (KRBH) using a modified Newman-Janis algorithm and
investigated diverse values of the parameter $\zeta$ in the cosmological term.
Throughout, we emphasized the ergoregion and identified $\zeta=3$ as a
favorable choice due to the absence of imaginary horizons. Shifting focus to
Hamiltonian dynamics, we constructed a framework for light propagation within
a non-magnetized electronic plasma. Assuming constant plasma parameters, we
derived analytical solutions for the temporal evolution of spacetime
coordinates in terms of the Mino time. Our attention turned to motion along
the latitude, revealing two turning points for the polar angle and expressing
the solution using elliptic integrals and Jacobi elliptic functions. We
highlighted the solution’s oscillatory nature within allowable polar angle
ranges. Analyzing the radial motion, we concentrated on spherical orbits,
pivotal for black hole imaging. Calculations involved critical impact
parameters and photon regions, and we discerned between deflecting
trajectories, distinguishing orbits easily escaping the black hole (DFK) from
those drawn into it (DSK). We derived analytical solutions for both cases and
discussed critical trajectories and plunge orbits. Simulating two-dimensional
light ray trajectories around the KRBH in plasma, our study encompassed all
orbit types. We calculated azimuth angle evolution, expressing analytical
solutions using elliptic integrals and highlighting plasma versus vacuum
differences in light propagation. Illustrative examples with an equatorial
disk as the illumination source showcased the impact of plasma on light.
Exploring the role of DFK in gravitational lensing, we presented the lens
equation and calculated deflection angles for light passing the black hole,
observing larger deflection angles in plasma compared to vacuum. Additionally,
we explored the COFK, focusing on spherical orbits and escape, demonstrating
how higher initial impact parameters led to smaller polar ranges. Analytical
solutions for plunge orbits were also provided. In conclusion, we revisited
spherical orbits, offering expressions for azimuth angle change and time delay
during each half orbit. These findings and preceding solutions lay the
foundation for future studies. This paper marks the initial step in
forthcoming research exploring strong lensing and accretion disk imaging for
KRBH in plasmic media, to be detailed in future papers.

## Acknowledgements

M.F. acknowledges financial support from Vicerrectoría de Investigación,
Desarrollo e Innovación - Universidad de Santiago de Chile (USACH), Proyecto
DICYT, Código 042331CM$\\_$Postdoc. J.R.V. is partially supported by the
Centro de Astrofísica de Valparaíso (CAV). M.F. would like to thank the
hospitality of the faculty of Instituto de Física y Astronomía (IFA) of
Universidad de Valparaíso during his visit to this institute, where this work
was accomplished.

##  Appendix A Derivation of the roots of a quartic equation

Here, we provide the method of solving the quartic equation $\xi^{c}=0$ that
generates the radii of polar orbits $r_{\mathrm{po}}$. This method is also
used to find the turning points in the radial motion. In fact, the equation
$\xi^{c}=0$ can be rewritten as

$P_{4}(r)=r^{4}-\frac{2}{N_{s}}r^{3}+\frac{3}{N_{s}}\left(M-a^{2}N_{s}\right)r^{2}-\frac{2a^{2}}{N_{s}}r-\frac{2a^{2}}{M}=0.$
(104)

Applying the change of variable $r=u-\frac{1}{2N_{s}}$, the above polynomial
equation reduces to

$P_{4}(u)=u^{4}+\mathrm{p}_{2}u^{2}+\mathrm{p}_{1}u+\mathrm{p}_{0},$ (105)

where

$\displaystyle\mathrm{p}_{2}=\frac{3(4MN_{s}-1)}{2N_{s}^{2}}-3a^{2},$ (106a)
$\displaystyle\mathrm{p}_{1}=\frac{6MN_{s}-5a^{2}N_{s}^{2}-1}{N_{s}^{3}},$
(106b)
$\displaystyle\mathrm{p}_{0}=-\frac{3+4N_{s}\left[7a^{2}N_{s}+2M\left(4a^{2}N_{s}^{2}-3\right)\right]}{16N_{s}^{4}}.$
(106c)

The solutions to this quartic are then obtained as (see the appendices of
Refs. fathi_gravitational_2021 ; fathi_spherical_2023 )

$\displaystyle
u_{1}=\bar{\mathrm{a}}+\sqrt{\bar{\mathrm{a}}^{2}-\bar{\mathrm{b}}},$ (107)
$\displaystyle
u_{2}=\bar{\mathrm{a}}-\sqrt{\bar{\mathrm{a}}^{2}-\bar{\mathrm{b}}},$ (108)
$\displaystyle
u_{3}=-\bar{\mathrm{a}}+\sqrt{\bar{\mathrm{a}}^{2}-\bar{\mathrm{c}}},$ (109)
$\displaystyle
u_{4}=-\bar{\mathrm{a}}-\sqrt{\bar{\mathrm{a}}^{2}-\bar{\mathrm{c}}},$ (110)

where

$\displaystyle\bar{\mathrm{a}}=\sqrt{\bar{\mathrm{u}}-\frac{\mathrm{p}_{2}}{6}},$
(111a)
$\displaystyle\bar{\mathrm{b}}=2\bar{\mathrm{a}}^{2}+\frac{\mathrm{p}_{2}}{2}+\frac{\mathrm{p}_{1}}{4\bar{\mathrm{a}}},$
(111b)
$\displaystyle\bar{\mathrm{c}}=2\bar{\mathrm{a}}^{2}+\frac{\mathrm{p}_{2}}{2}-\frac{\mathrm{p}_{1}}{4\bar{\mathrm{a}}},$
(111c)

and

$\bar{\mathrm{u}}=\sqrt{\frac{\bar{\epsilon}_{2}}{3}}\cosh\left(\frac{1}{3}\operatorname{arccosh}\left(3\bar{\epsilon}_{3}\sqrt{\frac{3}{\bar{\epsilon}_{2}^{3}}}\right)\right),$
(112)

with

$\displaystyle\bar{\epsilon}_{2}=\frac{\mathrm{p}_{2}^{2}}{12}+\mathrm{p}_{0},$
(113a)
$\displaystyle\bar{\mathrm{\epsilon}}_{3}=\frac{\mathrm{p}_{2}^{3}}{216}-\frac{\mathrm{p}_{2}\mathrm{p}_{0}}{6}+\frac{\mathrm{p}_{1}^{2}}{16}.$
(113b)

The positive radii of polar photon orbits outside the event horizon are then
given in terms of the values in Eqs. (107) and (108). However, the only
reliable solution inside the domain $r\in[r_{+},r_{++}]$ is based on Eq.
(108), and is obtained as

$r_{\mathrm{po}}=u_{2}-\frac{1}{2N_{s}}.$ (114)

Employing the same methods, the roots of the characteristic polynomial (67)
are obtained as

$\displaystyle r_{d}$ $\displaystyle=$
$\displaystyle\tilde{\alpha}+\sqrt{\tilde{\alpha}^{2}-\tilde{\beta}}-{\bar{A}\over
4},$ (115) $\displaystyle r_{f}$ $\displaystyle=$
$\displaystyle\tilde{\alpha}-\sqrt{\tilde{\alpha}^{2}-\tilde{\beta}}-{\bar{A}\over
4},$ (116) $\displaystyle r_{3}$ $\displaystyle=$
$\displaystyle-\tilde{\alpha}+\sqrt{\tilde{\alpha}^{2}-\tilde{\gamma}}-{\bar{A}\over
4},$ (117) $\displaystyle r_{4}$ $\displaystyle=$
$\displaystyle-\tilde{\alpha}-\sqrt{\tilde{\alpha}^{2}-\tilde{\gamma}}-{\bar{A}\over
4},$ (118)

where $\bar{A}=N_{s}\left[\left(\xi-a\right)^{2}+\eta_{r}\right]$, and

$\displaystyle\tilde{\alpha}=\sqrt{\Xi-{\tilde{p}\over 3}},$ (119a)
$\displaystyle\tilde{\beta}=2\tilde{\alpha}^{2}+{\mathcal{A}\over
2}+{\mathcal{B}\over 4\tilde{\alpha}},$ (119b)
$\displaystyle\tilde{\gamma}=2\tilde{\alpha}^{2}+{\mathcal{A}\over
2}-{\mathcal{B}\over 4\tilde{\alpha}},$ (119c)

with

$\Xi=\sqrt{{\chi_{2}\over
3}}\cosh\left(\frac{1}{3}\operatorname{arccosh}\left(3\chi_{3}\sqrt{{3\over\chi_{2}^{3}}}\right)\right),$
(120)

in which

$\displaystyle\chi_{2}=4\left({\tilde{p}^{2}\over 3}-\tilde{q}\right),$ (121a)
$\displaystyle\chi_{3}=4\left({2\tilde{p}^{3}\over
27}+{\tilde{p}\,\tilde{q}\over 3}-\tilde{r}\right),$ (121b)

where

$\displaystyle\tilde{r}=-{\mathcal{B}^{2}\over 64},$ (122a)
$\displaystyle\tilde{q}={1\over 4}\left({\mathcal{A}^{2}\over
4}-\mathcal{C}\right),$ (122b) $\displaystyle\tilde{p}={\mathcal{A}\over 2},$
(122c)

and

$\displaystyle\mathcal{A}=\mathrm{a}_{2}-{3\mathrm{a}_{3}^{2}\over 8},$ (123a)
$\displaystyle\mathcal{B}=\mathrm{a}_{1}+{\mathrm{a}_{3}^{3}\over
8}-{\mathrm{a}_{3}\,\mathrm{a}_{2}\over 2},$ (123b)
$\displaystyle\mathcal{C}=\mathrm{a}_{0}-{3\mathrm{a}_{3}^{4}\over
256}+{\mathrm{a}_{3}^{2}\,\mathrm{a}_{2}\over
16}-{\mathrm{a}_{3}\,\bar{C}\over 4},$ (123c)

with $\bar{C}=2M\left[\left(\xi-a\right)^{2}+\eta_{r}\right].$

##  Appendix B Coefficients and the analytical relations for the integrals of
$\Phi_{r}(\gamma)$ and $t_{r}(\gamma)$

For the DFK, we have defined in Eq. (82) that

$\displaystyle
A_{1}^{\mathrm{DFK}}=-\frac{La^{2}-2Ma\omega_{0}r_{-}-N_{s}a\omega_{0}r_{-}^{3}}{(r_{+}-r_{-})(r_{++}-r_{-})(r_{d}-r_{-})(r_{f}-r_{-})},$
(124a) $\displaystyle
A_{2}^{\mathrm{DFK}}=\frac{La^{2}-2Ma\omega_{0}r_{+}-N_{s}a\omega_{0}r_{+}^{3}}{(r_{+}-r_{-})(r_{++}-r_{+})(r_{d}-r_{+})(r_{f}-r_{+})},$
(124b) $\displaystyle
A_{3}^{\mathrm{DFK}}=\frac{-La^{2}+2Ma\omega_{0}r_{++}+N_{s}a\omega_{0}r_{++}^{3}}{(r_{++}-r_{+})(r_{++}-r_{-})(r_{++}-r_{d})(r_{++}-r_{f})}.$
(124c)

Furthermore, the elliptic integrals in the summation in Eq. (82) are expressed
as

$\displaystyle
I_{1}^{\mathrm{DFK}}(\gamma)=\Big{[}(r_{d}-r_{f})\operatorname{\mathbf{\Pi}}_{i}^{-}+(r_{d}-r_{-})\operatorname{\mathbf{F}}_{i}\Big{]}\Big{|}_{r_{\mathrm{s}}}^{r_{\mathrm{o}}(\gamma)},$
(125a) $\displaystyle
I_{2}^{\mathrm{DFK}}(\gamma)=\Big{[}(r_{d}-r_{f})\operatorname{\mathbf{\Pi}}_{i}^{+}+(r_{d}-r_{+})\operatorname{\mathbf{F}}_{i}\Big{]}\Big{|}_{r_{\mathrm{s}}}^{r_{\mathrm{o}}(\gamma)},$
(125b) $\displaystyle
I_{3}^{\mathrm{DFK}}(\gamma)=\Big{[}(r_{d}-r_{f})\operatorname{\mathbf{\Pi}}_{i}^{++}-(r_{++}-r_{d})\operatorname{\mathbf{F}}_{i}\Big{]}\Big{|}_{r_{\mathrm{s}}}^{r_{\mathrm{o}}(\gamma)},$
(125c)

with

$\displaystyle\operatorname{\mathbf{\Pi}}_{i}^{-}=\operatorname{\mathbf{\Pi}}\left(\frac{(r_{d}-r_{4})(r_{f}-r_{-})}{(r_{d}-r_{-})(r_{f}-r_{4})};\tilde{\varphi}_{i}\Big{|}\tilde{\mathcal{k}}\right),$
(126a)
$\displaystyle\operatorname{\mathbf{\Pi}}_{i}^{+}=\operatorname{\mathbf{\Pi}}\left(\frac{(r_{d}-r_{4})(r_{f}-r_{+})}{(r_{d}-r_{+})(r_{f}-r_{4})};\tilde{\varphi}_{i}\Big{|}\tilde{\mathcal{k}}\right),$
(126b)
$\displaystyle\operatorname{\mathbf{\Pi}}_{i}^{++}=\operatorname{\mathbf{\Pi}}\left(\frac{(r_{d}-r_{4})(r_{++}-r_{f})}{(r_{++}-r_{d})(r_{f}-r_{4})};\tilde{\varphi}_{i}\Big{|}\tilde{\mathcal{k}}\right),$
(126c)
$\displaystyle\operatorname{\mathbf{F}}_{i}=\operatorname{\mathbf{F}}(\tilde{\varphi}_{i}|\tilde{\mathcal{k}}),$
(126d)

in which $\tilde{\varphi}_{i}$ changes from $\tilde{\varphi}_{\mathrm{s}}$ to
$\tilde{\varphi}_{\mathrm{o}}(\gamma)$. For the DSK, the integrals
$I_{j}^{\mathrm{DSK}}$ are defined as

$\displaystyle
I_{1}^{\mathrm{DSK}}(\gamma)=\Big{[}(r_{d}-r_{f})\tilde{\operatorname{\mathbf{\Pi}}}_{i}^{-}+(r_{f}-r_{-})\tilde{\operatorname{\mathbf{F}}}_{i}\Big{]}\Big{|}_{r_{\mathrm{s}}}^{r_{\mathrm{o}}(\gamma)},$
(127a) $\displaystyle
I_{2}^{\mathrm{DSK}}(\gamma)=\Big{[}(r_{d}-r_{f})\tilde{\operatorname{\mathbf{\Pi}}}_{i}^{+}+(r_{f}-r_{+})\tilde{\operatorname{\mathbf{F}}}_{i}\Big{]}\Big{|}_{r_{\mathrm{s}}}^{r_{\mathrm{o}}(\gamma)},$
(127b) $\displaystyle
I_{3}^{\mathrm{DSK}}(\gamma)=\Big{[}(r_{d}-r_{f})\tilde{\operatorname{\mathbf{\Pi}}}_{i}^{++}-(r_{++}-r_{f})\tilde{\operatorname{\mathbf{F}}}_{i}\Big{]}\Big{|}_{r_{\mathrm{s}}}^{r_{\mathrm{o}}(\gamma)},$
(127c)

in which

$\displaystyle\tilde{\operatorname{\mathbf{\Pi}}}_{i}^{-}=\operatorname{\mathbf{\Pi}}\left(\frac{(r_{f}-r_{3})(r_{d}-r_{-})}{(r_{d}-r_{3})(r_{f}-r_{-})};\stackon[-1.3ex]{\stackon[-1ex]{\varphi}{\mathchar
12414\relax}}{\mathchar 12414\relax}_{i}\Big{|}\tilde{\mathcal{k}}\right),$
(128a)
$\displaystyle\tilde{\operatorname{\mathbf{\Pi}}}_{i}^{+}=\operatorname{\mathbf{\Pi}}\left(\frac{(r_{f}-r_{3})(r_{d}-r_{+})}{(r_{d}-r_{3})(r_{f}-r_{+})};\stackon[-1.3ex]{\stackon[-1ex]{\varphi}{\mathchar
12414\relax}}{\mathchar 12414\relax}_{i}\Big{|}\tilde{\mathcal{k}}\right),$
(128b)
$\displaystyle\tilde{\operatorname{\mathbf{\Pi}}}_{i}^{++}=\operatorname{\mathbf{\Pi}}\left(\frac{(r_{f}-r_{3})(r_{++}-r_{d})}{(r_{d}-r_{3})(r_{++}-r_{f})};\stackon[-1.3ex]{\stackon[-1ex]{\varphi}{\mathchar
12414\relax}}{\mathchar 12414\relax}_{i}\Big{|}\tilde{\mathcal{k}}\right),$
(128c)
$\displaystyle\tilde{\operatorname{\mathbf{F}}}_{i}=\operatorname{\mathbf{F}}(\stackon[-1.3ex]{\stackon[-1ex]{\varphi}{\mathchar
12414\relax}}{\mathchar 12414\relax}_{i}|\tilde{\mathcal{k}}),$ (128d)

where $\stackon[-1.3ex]{\stackon[-1ex]{\varphi}{\mathchar
12414\relax}}{\mathchar 12414\relax}_{i}$ changes from
$\stackon[-1.3ex]{\stackon[-1ex]{\varphi}{\mathchar 12414\relax}}{\mathchar
12414\relax}_{\mathrm{s}}$ to
$\stackon[-1.3ex]{\stackon[-1ex]{\varphi}{\mathchar 12414\relax}}{\mathchar
12414\relax}_{\mathrm{o}}(\gamma)$.

The parameters involved in the $r$-dependent expression for the COFK, as
provided in Eq. (89), are

$\displaystyle
A_{1}^{\mathrm{COFK}}=\frac{La^{2}-2Ma\omega_{0}r_{c}-N_{s}a\omega_{0}r_{c}^{3}}{(r_{c}-r_{-})(r_{c}-r_{+})(r_{++}-r_{c})\sqrt{(r_{c}-r_{3})(r_{c}-r_{4})}},$
(129a) $\displaystyle
A_{2}^{\mathrm{COFK}}=\frac{La^{2}-2Ma\omega_{0}r_{-}-N_{s}a\omega_{0}r_{-}^{3}}{(r_{c}-r_{-})(r_{+}-r_{-})(r_{++}-r_{-})\sqrt{(r_{-}-r_{3})(r_{-}-r_{4})}},$
(129b) $\displaystyle
A_{3}^{\mathrm{COFK}}=\frac{-La^{2}+2Ma\omega_{0}r_{+}+N_{s}a\omega_{0}r_{+}^{3}}{(r_{c}-r_{+})(r_{+}-r_{-})(r_{++}-r_{+})\sqrt{(r_{+}-r_{3})(r_{+}-r_{4})}},$
(129c) $\displaystyle
A_{4}^{\mathrm{COFK}}=\frac{-La^{2}+2Ma\omega_{0}r_{++}+N_{s}a\omega_{0}r_{++}^{3}}{(r_{++}-r_{c})(r_{++}-r_{-})(r_{++}-r_{+})\sqrt{(r_{++}-r_{3})(r_{++}-r_{4})}},$
(129d)

together with the inverse hyperbolic functions

$\displaystyle
I_{1}^{\mathrm{COFK}}(\gamma)=\left.\operatorname{arccosh}\left(\sqrt{\frac{(r_{c}-r_{4})(r-r_{3})}{(r_{3}-r_{4})(r-r_{c})}}\,\right)\right|_{r_{\mathrm{s}}}^{r_{\mathrm{o}}(\gamma)},$
(130a) $\displaystyle
I_{2}^{\mathrm{COFK}}(\gamma)=\left.\operatorname{arccosh}\left(\sqrt{\frac{(r_{-}-r_{4})(r-r_{3})}{(r_{3}-r_{4})(r-r_{-})}}\,\right)\right|_{r_{\mathrm{s}}}^{r_{\mathrm{o}}(\gamma)},$
(130b) $\displaystyle
I_{3}^{\mathrm{COFK}}(\gamma)=\left.\operatorname{arccosh}\left(\sqrt{\frac{(r_{+}-r_{4})(r-r_{3})}{(r_{3}-r_{4})(r-r_{+})}}\,\right)\right|_{r_{\mathrm{s}}}^{r_{\mathrm{o}}(\gamma)},$
(130c) $\displaystyle
I_{4}^{\mathrm{COFK}}(\gamma)=\left.\operatorname{arccosh}\left(\sqrt{\frac{(r_{++}-r_{4})(r-r_{3})}{(r_{3}-r_{4})(r_{++}-r)}}\,\right)\right|_{r_{\mathrm{s}}}^{r_{\mathrm{o}}(\gamma)},$
(130d)

in which, $r_{\mathrm{o}}(\gamma)$ corresponds to the analytical solution (73)
for the COFK.

In the solution (90) we have defined

$\displaystyle g_{p}=\frac{1}{\sqrt{\tilde{A}\tilde{B}}},$ (131a)
$\displaystyle\tilde{\varphi}_{p,\mathrm{o}}=\arccos\left(\frac{r_{3}\tilde{B}-r_{4}\tilde{A}+(\tilde{A}-\tilde{B})r_{\mathrm{o}}(\gamma)}{-r_{3}\tilde{B}-r_{4}\tilde{A}+(\tilde{A}-\tilde{B})r_{\mathrm{o}}(\gamma)}\right),$
(131b)

and the expressions for $\tilde{A}$, $\tilde{B}$, $\tilde{\mathcal{k}}_{p}$
and $\tilde{\varphi}_{p,\mathrm{s}}$ have been given in Eqs. (77), and
$r_{\mathrm{o}}(\gamma)$ corresponds to the analytical solution (76) for
plunge orbits. Furthermore, the coefficients $A_{j}^{\mathrm{PO}}$ in Eq. (90)
are defined as

$\displaystyle
A_{1}^{\mathrm{PO}}=\frac{La^{2}-2Ma\omega_{0}r_{-}-N_{s}a\omega_{0}r_{-}^{3}}{(r_{+}-r_{-})(r_{++}-r_{-})},$
(132a) $\displaystyle
A_{2}^{\mathrm{PO}}=\frac{-La^{2}+2Ma\omega_{0}r_{+}+N_{s}a\omega_{0}r_{+}^{3}}{(r_{+}-r_{-})(r_{++}-r_{+})},$
(132b) $\displaystyle
A_{3}^{\mathrm{PO}}=\frac{-La^{2}+2Ma\omega_{0}r_{++}+N_{s}a\omega_{0}r_{++}^{3}}{(r_{++}-r_{-})(r_{++}-r_{+})}.$
(132c)

The included elliptic integrals in this solution are given by

$\displaystyle
I_{1}^{\mathrm{PO}}(\gamma)=\left.\frac{g_{p}(\tilde{B}-\tilde{A})\left[(\accentset{\circ}{\alpha}_{-}-\alpha_{p})\accentset{\circ}{\mathbf{R}}_{i}^{-}+\alpha_{p}\accentset{\circ}{\mathbf{F}}_{i}\right]}{r_{3}\tilde{B}-r_{4}\tilde{A}-(\tilde{A}+\tilde{B})r_{-}}\right|_{r_{\mathrm{s}}}^{r_{\mathrm{o}}(\gamma)},$
(133a) $\displaystyle
I_{2}^{\mathrm{PO}}(\gamma)=\left.\frac{g_{p}(\tilde{B}-\tilde{A})\left[(\accentset{\circ}{\alpha}_{+}-\alpha_{p})\accentset{\circ}{\mathbf{R}}_{i}^{+}+\alpha_{p}\accentset{\circ}{\mathbf{F}}_{i}\right]}{r_{3}\tilde{B}-r_{4}\tilde{A}-(\tilde{A}+\tilde{B})r_{+}}\right|_{r_{\mathrm{s}}}^{r_{\mathrm{o}}(\gamma)},$
(133b) $\displaystyle
I_{3}^{\mathrm{PO}}(\gamma)=\left.\frac{g_{p}(\tilde{B}-\tilde{A})\left[(\accentset{\circ}{\alpha}_{++}-\alpha_{p})\accentset{\circ}{\mathbf{R}}_{i}^{++}+\alpha_{p}\accentset{\circ}{\mathbf{F}}_{i}\right]}{r_{3}\tilde{B}-r_{4}\tilde{A}-(\tilde{A}+\tilde{B})r_{++}}\right|_{r_{\mathrm{s}}}^{r_{\mathrm{o}}(\gamma)},$
(133c)

in which
$\accentset{\circ}{\mathbf{F}}i=\operatorname{\mathbf{F}}(\tilde{\varphi}_{p,i}|\tilde{\mathcal{k}})$,
$\alpha_{p}=(\tilde{A}+\tilde{B})/(\tilde{B}-\tilde{A})$, and

$\displaystyle\accentset{\circ}{\alpha}_{-}=\frac{\tilde{B}r_{3}+\tilde{A}r_{4}-(\tilde{A}+\tilde{B})r_{-}}{\tilde{B}r_{3}-\tilde{A}r_{4}+(\tilde{A}+\tilde{B})r_{-}},$
(134a)
$\displaystyle\accentset{\circ}{\alpha}_{+}=\frac{\tilde{B}r_{3}+\tilde{A}r_{4}-(\tilde{A}+\tilde{B})r_{+}}{\tilde{B}r_{3}-\tilde{A}r_{4}+(\tilde{A}+\tilde{B})r_{+}},$
(134b)
$\displaystyle\accentset{\circ}{\alpha}_{++}=\frac{\tilde{B}r_{3}+\tilde{A}r_{4}-(\tilde{A}+\tilde{B})r_{++}}{\tilde{B}r_{3}-\tilde{A}r_{4}+(\tilde{A}+\tilde{B})r_{++}},$
(134c)

so that

$\displaystyle\accentset{\circ}{\mathbf{R}}_{i}^{-}=\frac{1}{1-\accentset{\circ}{\alpha}_{-}^{2}}\left[\accentset{\circ}{\mathbf{\Pi}}_{i}^{-}-\accentset{\circ}{\alpha}_{-}\accentset{\circ}{\mathbf{f}}_{i}^{-}\right],$
(135a)
$\displaystyle\accentset{\circ}{\mathbf{R}}_{i}^{+}=\frac{1}{1-\accentset{\circ}{\alpha}_{+}^{2}}\left[\accentset{\circ}{\mathbf{\Pi}}_{i}^{+}-\accentset{\circ}{\alpha}_{+}\accentset{\circ}{\mathbf{f}}_{i}^{+}\right],$
(135b)
$\displaystyle\accentset{\circ}{\mathbf{R}}_{i}^{++}=\frac{1}{1-\accentset{\circ}{\alpha}_{++}^{2}}\left[\accentset{\circ}{\mathbf{\Pi}}_{i}^{++}-\accentset{\circ}{\alpha}_{++}\accentset{\circ}{\mathbf{f}}_{i}^{++}\right],$
(135c)

where

$\displaystyle\accentset{\circ}{\mathbf{\Pi}}_{i}^{-}=\operatorname{\mathbf{\Pi}}\left(\frac{\accentset{\circ}{\alpha}_{-}^{2}}{\accentset{\circ}{\alpha}_{-}^{2}-1};\tilde{\varphi}_{p,i}\Big{|}\tilde{\mathcal{k}}_{p}\right),$
(136a)
$\displaystyle\accentset{\circ}{\mathbf{\Pi}}_{i}^{+}=\operatorname{\mathbf{\Pi}}\left(\frac{\accentset{\circ}{\alpha}_{+}^{2}}{\accentset{\circ}{\alpha}_{+}^{2}-1};\tilde{\varphi}_{p,i}\Big{|}\tilde{\mathcal{k}}_{p}\right),$
(136b)
$\displaystyle\accentset{\circ}{\mathbf{\Pi}}_{i}^{++}=\operatorname{\mathbf{\Pi}}\left(\frac{\accentset{\circ}{\alpha}_{++}^{2}}{\accentset{\circ}{\alpha}_{++}^{2}-1};\tilde{\varphi}_{p,i}\Big{|}\tilde{\mathcal{k}}_{p}\right),$
(136c)

with

$\displaystyle\accentset{\circ}{\mathbf{f}}_{i}^{-}=\sqrt{\frac{\accentset{\circ}{\alpha}_{-}^{2}-1}{\tilde{\mathcal{k}}_{p}+(1-\tilde{\mathcal{k}}_{p})\accentset{\circ}{\alpha}_{-}^{2}}}\,\ln\left|\frac{\sqrt{\tilde{\mathcal{k}}_{p}+(1-\tilde{\mathcal{k}}_{p})\accentset{\circ}{\alpha}_{-}^{2}}\,\operatorname{\mathbf{dn}}(\tilde{\mathbf{u}}_{p,i}|\tilde{\mathcal{k}}_{p})+\sqrt{\accentset{\circ}{\alpha}_{-}^{2}-1}\,\operatorname{\mathbf{sn}}(\tilde{\mathbf{u}}_{p,i}|\tilde{\mathcal{k}}_{p})}{\sqrt{\tilde{\mathcal{k}}_{p}+(1-\tilde{\mathcal{k}}_{p})\accentset{\circ}{\alpha}_{-}^{2}}\,\operatorname{\mathbf{dn}}(\tilde{\mathbf{u}}_{p,i}|\tilde{\mathcal{k}}_{p})-\sqrt{\accentset{\circ}{\alpha}_{-}^{2}-1}\,\operatorname{\mathbf{sn}}(\tilde{\mathbf{u}}_{p,i}|\tilde{\mathcal{k}}_{p})}\right|,$
(137a)
$\displaystyle\accentset{\circ}{\mathbf{f}}_{i}^{+}=\sqrt{\frac{\accentset{\circ}{\alpha}_{+}^{2}-1}{\tilde{\mathcal{k}}_{p}+(1-\tilde{\mathcal{k}}_{p})\accentset{\circ}{\alpha}_{+}^{2}}}\,\ln\left|\frac{\sqrt{\tilde{\mathcal{k}}_{p}+(1-\tilde{\mathcal{k}}_{p})\accentset{\circ}{\alpha}_{+}^{2}}\,\operatorname{\mathbf{dn}}(\tilde{\mathbf{u}}_{p,i}|\tilde{\mathcal{k}}_{p})+\sqrt{\accentset{\circ}{\alpha}_{+}^{2}-1}\,\operatorname{\mathbf{sn}}(\tilde{\mathbf{u}}_{p,i}|\tilde{\mathcal{k}}_{p})}{\sqrt{\tilde{\mathcal{k}}_{p}+(1-\tilde{\mathcal{k}}_{p})\accentset{\circ}{\alpha}_{+}^{2}}\,\operatorname{\mathbf{dn}}(\tilde{\mathbf{u}}_{p,i}|\tilde{\mathcal{k}}_{p})-\sqrt{\accentset{\circ}{\alpha}_{+}^{2}-1}\,\operatorname{\mathbf{sn}}(\tilde{\mathbf{u}}_{p,i}|\tilde{\mathcal{k}}_{p})}\right|,$
(137b)
$\displaystyle\accentset{\circ}{\mathbf{f}}_{i}^{++}=\sqrt{\frac{\accentset{\circ}{\alpha}_{++}^{2}-1}{\tilde{\mathcal{k}}_{p}+(1-\tilde{\mathcal{k}}_{p})\accentset{\circ}{\alpha}_{++}^{2}}}\,\ln\left|\frac{\sqrt{\tilde{\mathcal{k}}_{p}+(1-\tilde{\mathcal{k}}_{p})\accentset{\circ}{\alpha}_{++}^{2}}\,\operatorname{\mathbf{dn}}(\tilde{\mathbf{u}}_{p,i}|\tilde{\mathcal{k}}_{p})+\sqrt{\accentset{\circ}{\alpha}_{++}^{2}-1}\,\operatorname{\mathbf{sn}}(\tilde{\mathbf{u}}_{p,i}|\tilde{\mathcal{k}}_{p})}{\sqrt{\tilde{\mathcal{k}}_{p}+(1-\tilde{\mathcal{k}}_{p})\accentset{\circ}{\alpha}_{++}^{2}}\,\operatorname{\mathbf{dn}}(\tilde{\mathbf{u}}_{p,i}|\tilde{\mathcal{k}}_{p})-\sqrt{\accentset{\circ}{\alpha}_{++}^{2}-1}\,\operatorname{\mathbf{sn}}(\tilde{\mathbf{u}}_{p,i}|\tilde{\mathcal{k}}_{p})}\right|,$
(137c)

where $\operatorname{\mathbf{dn}}(x|\tilde{\mathcal{k}}_{p})$ is the Jacobian
elliptic delta amplitude of the argument $x$, and
$\tilde{\mathbf{u}}_{p,i}=\operatorname{\mathbf{cn}}^{-1}(\cos(\tilde{\varphi}_{p,i})|\tilde{\mathcal{k}}_{p})$,
in which $\operatorname{\mathbf{cn}}^{-1}(x|\tilde{\mathcal{k}}_{p})$ is the
inverse Jacobian elliptic cosine function of the argument $x$
byrd_handbook_1971 .

The coefficients $B_{j}^{t}$ in the solution (94) are given as

$\displaystyle B_{1}^{t}=\frac{-a^{4}\omega_{0}+r_{-}\left[2aLM-
r_{-}\left(2a^{2}\omega_{0}-aLN_{s}r_{-}+\omega_{0}r_{-}^{2}\right)\right]}{MN_{s}(r_{++}-r_{-})(r_{+}-r_{-})(r_{d}-r_{-})(r_{f}-r_{-})\sqrt{(r_{d}-r_{3})(r_{f}-r_{4})}},$
(138a) $\displaystyle
B_{2}^{t}=\frac{a^{4}\omega_{0}+r_{+}\left[-2aLM+r_{+}\left(2a^{2}\omega_{0}-aLN_{s}r_{+}+\omega_{0}r_{+}^{2}\right)\right]}{MN_{s}(r_{++}-r_{+})(r_{+}-r_{-})(r_{d}-r_{+})(r_{f}-r_{+})\sqrt{(r_{d}-r_{3})(r_{f}-r_{4})}},$
(138b) $\displaystyle
B_{3}^{t}=-\frac{a^{4}\omega_{0}+r_{++}\left[-2aLM+r_{++}\left(2a^{2}\omega_{0}-aLN_{s}r_{++}+\omega_{0}r_{++}^{2}\right)\right]}{MN_{s}(r_{++}-r_{+})(r_{++}-r_{-})(r_{++}-r_{d})(r_{++}-r_{f})\sqrt{(r_{d}-r_{3})(r_{f}-r_{4})}}.$
(138c)

## References

  * (1) A. G. Riess, L. Strolger, J. Tonry, S. Casertano, H. C. Ferguson, B. Mobasher, P. Challis, A. V. Filippenko, S. Jha, W. Li, R. Chornock, R. P. Kirshner, B. Leibundgut, M. Dickinson, M. Livio, M. Giavalisco, C. C. Steidel, T. Benitez, and Z. Tsvetanov, “Type Ia Supernova Discoveries at z > 1 from the Hubble Space Telescope : Evidence for Past Deceleration and Constraints on Dark Energy Evolution,” The Astrophysical Journal, vol. 607, pp. 665–687, June 2004\.
  * (2) E. Komatsu, K. M. Smith, J. Dunkley, C. L. Bennett, B. Gold, G. Hinshaw, N. Jarosik, D. Larson, M. R. Nolta, L. Page, D. N. Spergel, M. Halpern, R. S. Hill, A. Kogut, M. Limon, S. S. Meyer, N. Odegard, G. S. Tucker, J. L. Weiland, E. Wollack, and E. L. Wright, “SEVEN-YEAR WILKINSON MICROWAVE ANISOTROPY PROBE ( WMAP ) OBSERVATIONS: COSMOLOGICAL INTERPRETATION,” The Astrophysical Journal Supplement Series, vol. 192, p. 18, Feb. 2011.
  * (3) M. S. Turner, “Coherent scalar-field oscillations in an expanding universe,” Physical Review D, vol. 28, pp. 1243–1247, Sept. 1983.
  * (4) W. H. Press, B. S. Ryden, and D. N. Spergel, “Single mechanism for generating large-scale structure and providing dark missing matter,” Physical Review Letters, vol. 64, pp. 1084–1087, Mar. 1990.
  * (5) J. Dubinski and R. G. Carlberg, “The structure of cold dark matter halos,” The Astrophysical Journal, vol. 378, p. 496, Sept. 1991.
  * (6) S.-J. Sin, “Late-time phase transition and the galactic halo as a Bose liquid,” Physical Review D, vol. 50, pp. 3650–3654, Sept. 1994.
  * (7) J. F. Navarro, C. S. Frenk, and S. D. M. White, “The Structure of Cold Dark Matter Halos,” The Astrophysical Journal, vol. 462, p. 563, May 1996.
  * (8) J. F. Navarro, C. S. Frenk, and S. D. M. White, “A Universal Density Profile from Hierarchical Clustering,” The Astrophysical Journal, vol. 490, pp. 493–508, Dec. 1997.
  * (9) W. Hu, R. Barkana, and A. Gruzinov, “Fuzzy Cold Dark Matter: The Wave Properties of Ultralight Particles,” Physical Review Letters, vol. 85, pp. 1158–1161, Aug. 2000.
  * (10) P. J. E. Peebles, “Fluid Dark Matter,” The Astrophysical Journal, vol. 534, pp. L127–L129, May 2000.
  * (11) P. J. E. Peebles and B. Ratra, “The cosmological constant and dark energy,” Reviews of Modern Physics, vol. 75, pp. 559–606, Apr. 2003.
  * (12) V. V. Kiselev, “Quintessential solution of dark matter rotation curves and its simulation by extra dimensions,” arXiv e-prints, pp. gr–qc/0303031, Mar. 2003.
  * (13) M. Viel, J. Lesgourgues, M. G. Haehnelt, S. Matarrese, and A. Riotto, “Constraining warm dark matter candidates including sterile neutrinos and light gravitinos with WMAP and the Lyman-$\alpha$ forest,” Physical Review D, vol. 71, p. 063534, Mar. 2005.
  * (14) L. Amendola and R. Barbieri, “Dark matter from an ultra-light pseudo-Goldsone-boson,” Physics Letters B, vol. 642, pp. 192–196, Nov. 2006.
  * (15) E. J. Copeland, M. Sami, and S. Tsujikawa, “DYNAMICS OF DARK ENERGY,” International Journal of Modern Physics D, vol. 15, pp. 1753–1935, Nov. 2006.
  * (16) A. De La Macorra and T. Matos, “Dark Energy and Dark Matter,” in AIP Conference Proceedings, vol. 857, (Morelia, Michoacan (Mexico)), pp. 191–204, AIP, 2006.
  * (17) F. Rahaman, K. Nandi, A. Bhadra, M. Kalam, and K. Chakraborty, “Perfect fluid dark matter,” Physics Letters B, vol. 694, pp. 10–15, Oct. 2010.
  * (18) F. Rahaman, P. K. F. Kuhfittig, K. Chakraborty, M. Kalam, and D. Hossain, “Modeling Galactic Halos with Predominantly Quintessential Matter,” International Journal of Theoretical Physics, vol. 50, pp. 2655–2665, Sept. 2011.
  * (19) M. Li, X.-D. Li, S. Wang, and Y. Wang, “Dark Energy,” Communications in Theoretical Physics, vol. 56, pp. 525–604, Sept. 2011.
  * (20) H.-Y. Schive, T. Chiueh, and T. Broadhurst, “Cosmic structure as the quantum interference of a coherent dark wave,” Nature Physics, vol. 10, pp. 496–499, July 2014.
  * (21) D. J. Marsh, “Axion cosmology,” Physics Reports, vol. 643, pp. 1–79, July 2016.
  * (22) L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten, “Ultralight scalars as cosmological dark matter,” Physical Review D, vol. 95, p. 043541, Feb. 2017\.
  * (23) S. Wang, Y. Wang, and M. Li, “Holographic dark energy,” Physics Reports, vol. 696, pp. 1–57, June 2017.
  * (24) P. Rastall, “Generalization of the einstein theory,” Physical Review D, vol. 6, no. 12, p. 3357 – 3359, 1972.  Cited by: 291.
  * (25) P. Rastall, “A theory of gravity,” Canadian Journal of Physics, vol. 54, pp. 66–75, Jan. 1976.
  * (26) C. E. M. Batista, M. H. Daouda, J. C. Fabris, O. F. Piattella, and D. C. Rodrigues, “Rastall cosmology and the $\lambda$CDM model,” Physical Review D, vol. 85, p. 084008, Apr. 2012.
  * (27) J. C. Fabris, O. F. Piattella, D. C. Rodrigues, C. E. M. Batista, and M. H. Daouda, “RASTALL COSMOLOGY,” International Journal of Modern Physics: Conference Series, vol. 18, pp. 67–76, Jan. 2012.
  * (28) J. C. Fabris, O. F. Piattella, D. C. Rodrigues, and M. H. Daouda, “Rastall’s cosmology and its observational constraints,” (Valparaiso, Chile), pp. 50–53, 2015.
  * (29) Y. Heydarzade, H. Moradpour, and F. Darabi, “Black hole solutions in Rastall theory,” Canadian Journal of Physics, vol. 95, pp. 1253–1256, Dec. 2017\.
  * (30) Y. Heydarzade and F. Darabi, “Black hole solutions surrounded by perfect fluid in Rastall theory,” Physics Letters B, vol. 771, pp. 365–373, Aug. 2017\.
  * (31) F. Darabi, K. Atazadeh, and Y. Heydarzade, “Einstein static universe in the Rastall theory of gravity,” The European Physical Journal Plus, vol. 133, p. 249, July 2018.
  * (32) I. P. Lobo, H. Moradpour, J. P. Morais Graça, and I. G. Salako, “Thermodynamics of black holes in Rastall gravity,” International Journal of Modern Physics D, vol. 27, p. 1850069, May 2018.
  * (33) F. Kottler, “Über die physikalischen Grundlagen der Einsteinschen Gravitationstheorie,” Annalen der Physik, vol. 361, pp. 401–462, Jan. 1918.
  * (34) B. Carter, “Black holes equilibrium states,” in Les Houches Summer School of Theoretical Physics: Black Holes, pp. 57–214, 1973.
  * (35) V. V. Kiselev, “Quintessence and black holes,” Classical and Quantum Gravity, vol. 20, pp. 1187–1197, Mar. 2003.
  * (36) B. Toshmatov, Z. Stuchlík, and B. Ahmedov, “Rotating black hole solutions with quintessential energy,” The European Physical Journal Plus, vol. 132, p. 98, Feb. 2017.
  * (37) R. Kumar and S. G. Ghosh, “Rotating black hole in Rastall theory,” The European Physical Journal C, vol. 78, p. 750, Sept. 2018.
  * (38) Z. Xu, X. Hou, X. Gong, and J. Wang, “Kerr–Newman-AdS black hole surrounded by perfect fluid matter in Rastall gravity,” The European Physical Journal C, vol. 78, p. 513, June 2018.
  * (39) F. W. Dyson, A. S. Eddington, and C. Davidson, “Ix. a determination of the deflection of light by the sun’s gravitational field, from observations made at the total eclipse of may 29, 1919,” Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character, vol. 220, no. 571-581, pp. 291–333, 1920.
  * (40) A. Einstein, “Über den einfluß der schwerkraft auf die ausbreitung des lichtes,” Annalen der Physik, vol. 340, no. 10, pp. 898–908, 1911.
  * (41) K. Akiyama et al., “First M87 Event Horizon Telescope Results. IV. Imaging the Central Supermassive Black Hole,” Astrophys. J. Lett., vol. 875, no. 1, p. L4, 2019.
  * (42) K. Akiyama et al., “First Sagittarius A∗ Event Horizon Telescope results. I. The shadow of the supermassive black hole in the center of the Milky Way.,” Astrophys. J. Lett., vol. 930, no. 2, p. L12, 2022.
  * (43) J. M. Bardeen, W. H. Press, and S. A. Teukolsky, “Rotating Black Holes: Locally Nonrotating Frames, Energy Extraction, and Scalar Synchrotron Radiation,” Astrophysical Journal, vol. 178, pp. 347–370, Dec. 1972.
  * (44) J. Bardeen, “Timelike and null geodesics in the Kerr metric,” in Les Houches Summer School of Theoretical Physics: Black Holes, pp. 215–240, 1973\.
  * (45) J. Bardeen, “Timelike and null geodesics in the Kerr metric,” in Les Houches Summer School of Theoretical Physics: Black Holes, pp. 215–240, 1973\.
  * (46) S. Chandrasekhar, The mathematical theory of black holes.  Oxford classic texts in the physical sciences, Oxford: Oxford Univ. Press, 2002.
  * (47) K. S. Virbhadra and G. F. R. Ellis, “Schwarzschild black hole lensing,” Physical Review D, vol. 62, p. 084003, Sep 2000.
  * (48) G. V. Kraniotis, “Precise relativistic orbits in Kerr and Kerr–(anti) de Sitter spacetimes,” Classical and Quantum Gravity, vol. 21, pp. 4743–4769, Oct. 2004.
  * (49) K. Beckwith and C. Done, “Extreme gravitational lensing near rotating black holes,” Monthly Notices of the Royal Astronomical Society, vol. 359, pp. 1217–1228, June 2005.
  * (50) G. V. Kraniotis, “Frame dragging and bending of light in Kerr and Kerr–(anti) de Sitter spacetimes,” Classical and Quantum Gravity, vol. 22, pp. 4391–4424, Nov. 2005.
  * (51) E. Hackmann and C. Lämmerzahl, “Complete Analytic Solution of the Geodesic Equation in Schwarzschild–(Anti-)de Sitter Spacetimes,” Physical Review Letters, vol. 100, p. 171101, May 2008.
  * (52) E. Hackmann and C. Lämmerzahl, “Geodesic equation in Schwarzschild-(anti-)de Sitter space-times: Analytical solutions and applications,” Physical Review D, vol. 78, p. 024035, July 2008.
  * (53) G. S. Bisnovatyi-Kogan and O. Y. Tsupko, “Strong gravitational lensing by Schwarzschild black holes,” Astrophysics, vol. 51, pp. 99–111, Jan. 2008\.
  * (54) V. Kagramanova, J. Kunz, E. Hackmann, and C. Lämmerzahl, “Analytic treatment of complete and incomplete geodesics in Taub-NUT space-times,” Physical Review D, vol. 81, p. 124044, June 2010.
  * (55) E. Hackmann, C. Lämmerzahl, V. Kagramanova, and J. Kunz, “Analytical solution of the geodesic equation in Kerr-(anti-) de Sitter space-times,” Physical Review D, vol. 81, p. 044020, Feb. 2010.
  * (56) E. Hackmann, B. Hartmann, C. Lämmerzahl, and P. Sirimachan, “Complete set of solutions of the geodesic equation in the space-time of a Schwarzschild black hole pierced by a cosmic string,” Physical Review D, vol. 81, p. 064016, Mar. 2010.
  * (57) V. Enolski, E. Hackmann, V. Kagramanova, J. Kunz, and C. Lämmerzahl, “Inversion of hyperelliptic integrals of arbitrary genus with application to particle motion in general relativity,” Journal of Geometry and Physics, vol. 61, pp. 899–921, May 2011.
  * (58) G. V. Kraniotis, “Precise analytic treatment of Kerr and Kerr-(anti) de Sitter black holes as gravitational lenses,” Classical and Quantum Gravity, vol. 28, p. 085021, Apr. 2011.
  * (59) V. Enolski, B. Hartmann, V. Kagramanova, J. Kunz, C. Lämmerzahl, and P. Sirimachan, “Inversion of a general hyperelliptic integral and particle motion in Hořava–Lifshitz black hole space-times,” Journal of Mathematical Physics, vol. 53, p. 012504, Jan. 2012.
  * (60) G. W. Gibbons and M. Vyska, “The application of Weierstrass elliptic functions to Schwarzschild null geodesics,” Classical and Quantum Gravity, vol. 29, p. 065016, Mar. 2012.
  * (61) G. Muñoz, “Orbits of massless particles in the Schwarzschild metric: Exact solutions,” American Journal of Physics, vol. 82, pp. 564–573, June 2014.
  * (62) G. V. Kraniotis, “Gravitational lensing and frame dragging of light in the Kerr–Newman and the Kerr–Newman (anti) de Sitter black hole spacetimes,” General Relativity and Gravitation, vol. 46, p. 1818, Nov. 2014.
  * (63) V. De Falco, M. Falanga, and L. Stella, “Approximate analytical calculations of photon geodesics in the Schwarzschild metric,” Astronomy & Astrophysics, vol. 595, p. A38, Nov. 2016.
  * (64) S. Soroushfar, R. Saffari, S. Kazempour, S. Grunau, and J. Kunz, “Detailed study of geodesics in the Kerr-Newman-(A)dS spacetime and the rotating charged black hole spacetime in f ( R ) gravity,” Physical Review D, vol. 94, p. 024052, July 2016.
  * (65) N. S. Barlow, S. J. Weinstein, and J. A. Faber, “An asymptotically consistent approximant for the equatorial bending angle of light due to Kerr black holes,” Classical and Quantum Gravity, vol. 34, p. 135017, July 2017.
  * (66) R. Uniyal, H. Nandan, and K. D. Purohit, “Null geodesics and observables around the Kerr–Sen black hole,” Classical and Quantum Gravity, vol. 35, p. 025003, Jan. 2018.
  * (67) J. R. Villanueva, F. Tapia, M. Molina, and M. Olivares, “Null paths on a toroidal topological black hole in conformal Weyl gravity,” The European Physical Journal C, vol. 78, p. 853, Oct. 2018.
  * (68) A. K. Chatterjee, K. Flathmann, H. Nandan, and A. Rudra, “Analytic solutions of the geodesic equation for Reissner-Nordström–(anti–)de Sitter black holes surrounded by different kinds of regular and exotic matter fields,” Physical Review D, vol. 100, p. 024044, July 2019.
  * (69) Y.-W. Hsiao, D.-S. Lee, and C.-Y. Lin, “Equatorial light bending around Kerr-Newman black holes,” Physical Review D, vol. 101, p. 064070, Mar. 2020.
  * (70) S. E. Gralla and A. Lupsasca, “Null geodesics of the Kerr exterior,” Physical Review D, vol. 101, p. 044032, Feb. 2020.
  * (71) S. H. Hendi, A. M. Tavakkoli, S. Panahiyan, B. E. Panah, and E. Hackmann, “Simulation of geodesic trajectory of charged BTZ black holes in massive gravity,” The European Physical Journal C, vol. 80, p. 524, June 2020.
  * (72) M. Fathi, M. Olivares, and J. R. Villanueva, “Classical tests on a charged Weyl black hole: bending of light, Shapiro delay and Sagnac effect,” The European Physical Journal C, vol. 80, p. 51, Jan. 2020.
  * (73) G. V. Kraniotis, “Gravitational redshift/blueshift of light emitted by geodesic test particles, frame-dragging and pericentre-shift effects, in the Kerr–Newman–de Sitter and Kerr–Newman black hole geometries,” The European Physical Journal C, vol. 81, p. 147, Feb. 2021\.
  * (74) M. Fathi, M. Olivares, and J. R. Villanueva, “Study of null and time-like geodesics in the exterior of a Schwarzschild black hole with quintessence and cloud of strings,” The European Physical Journal C, vol. 82, p. 629, July 2022.
  * (75) M. Fathi, M. Olivares, and J. R. Villanueva, “Spherical photon orbits around a rotating black hole with quintessence and cloud of strings,” The European Physical Journal Plus, vol. 138, p. 7, Jan. 2023.
  * (76) E. Omwoyo, H. Belich, J. C. Fabris, and H. Velten, “Black hole lensing in Kerr-de Sitter spacetimes,” arXiv e-prints, p. arXiv:2212.10492, Dec. 2022.
  * (77) K. Wang, C.-J. Feng, and T. Wang, “Image of Kerr-de Sitter black holes: An additional avenue for testing the cosmological constant,” arXiv e-prints, p. arXiv:2309.16944, Sept. 2023.
  * (78) J. Synge, Relativity: The General Theory.  No. v. 1 in North-Holland series in physics, North-Holland Publishing Company, 1960.
  * (79) G. S. Bisnovatyi-Kogan and O. Y. Tsupko, “Gravitational radiospectrometer,” Gravitation and Cosmology, vol. 15, pp. 20–27, Jan. 2009.
  * (80) G. S. Bisnovatyi-Kogan and O. Y. Tsupko, “Gravitational lensing in a non-uniform plasma: Gravitational lensing in a non-uniform plasma,” Monthly Notices of the Royal Astronomical Society, pp. no–no, May 2010.
  * (81) O. Y. Tsupko and G. S. Bisnovatyi-Kogan, “Gravitational lensing in plasma: Relativistic images at homogeneous plasma,” Physical Review D, vol. 87, p. 124009, June 2013.
  * (82) V. S. Morozova, B. J. Ahmedov, and A. A. Tursunov, “Gravitational lensing by a rotating massive object in a plasma,” Astrophysics and Space Science, vol. 346, pp. 513–520, Aug. 2013.
  * (83) G. S. Bisnovatyi-Kogan and O. Y. Tsupko, “Gravitational lensing in plasmic medium,” Plasma Physics Reports, vol. 41, pp. 562–581, July 2015.
  * (84) V. Perlick, O. Y. Tsupko, and G. S. Bisnovatyi-Kogan, “Influence of a plasma on the shadow of a spherically symmetric black hole,” Physical Review D, vol. 92, p. 104031, Nov. 2015.
  * (85) F. Atamurotov, B. Ahmedov, and A. Abdujabbarov, “Optical properties of black holes in the presence of a plasma: The shadow,” Physical Review D, vol. 92, p. 084005, Oct. 2015.
  * (86) A. Abdujabbarov, M. Amir, B. Ahmedov, and S. G. Ghosh, “Shadow of rotating regular black holes,” Physical Review D, vol. 93, p. 104004, May 2016.
  * (87) G. Bisnovatyi-Kogan and O. Tsupko, “Gravitational Lensing in Presence of Plasma: Strong Lens Systems, Black Hole Lensing and Shadow,” Universe, vol. 3, p. 57, July 2017.
  * (88) V. Perlick and O. Y. Tsupko, “Light propagation in a plasma on Kerr spacetime: Separation of the Hamilton-Jacobi equation and calculation of the shadow,” Physical Review D, vol. 95, p. 104003, May 2017.
  * (89) K. Schulze-Koops, V. Perlick, and D. J. Schwarz, “Sachs equations for light bundles in a cold plasma,” Classical and Quantum Gravity, vol. 34, p. 215006, Nov. 2017.
  * (90) A. Abdujabbarov, B. Toshmatov, Z. Stuchlík, and B. Ahmedov, “Shadow of the rotating black hole with quintessential energy in the presence of plasma,” International Journal of Modern Physics D, vol. 26, p. 1750051, May 2017\.
  * (91) C.-Q. Liu, C.-K. Ding, and J.-L. Jing, “Effects of Homogeneous Plasma on Strong Gravitational Lensing of Kerr Black Holes,” Chinese Physics Letters, vol. 34, p. 090401, Aug. 2017.
  * (92) S. Haroon, M. Jamil, K. Jusufi, K. Lin, and R. B. Mann, “Shadow and deflection angle of rotating black holes in perfect fluid dark matter with a cosmological constant,” Physical Review D, vol. 99, p. 044015, Feb. 2019\.
  * (93) T. Kimpson, K. Wu, and S. Zane, “Spatial dispersion of light rays propagating through a plasma in Kerr space–time,” Monthly Notices of the Royal Astronomical Society, vol. 484, pp. 2411–2419, Apr. 2019.
  * (94) G. Z. Babar, A. Z. Babar, and F. Atamurotov, “Optical properties of Kerr–Newman spacetime in the presence of plasma,” The European Physical Journal C, vol. 80, p. 761, Aug. 2020.
  * (95) H. C. D. L. Junior, L. C. B. Crispino, P. V. P. Cunha, and C. A. R. Herdeiro, “Spinning black holes with a separable Hamilton–Jacobi equation from a modified Newman–Janis algorithm,” The European Physical Journal C, vol. 80, p. 1036, Nov. 2020.
  * (96) J. Badía and E. F. Eiroa, “Shadow of axisymmetric, stationary and asymptotically flat black holes in the presence of plasma,” 6 2021.
  * (97) M. Fathi, M. Olivares, and J. R. Villanueva, “Analytical study of light ray trajectories in Kerr spacetime in the presence of an inhomogeneous anisotropic plasma,” The European Physical Journal C, vol. 81, p. 987, Nov. 2021.
  * (98) R. Kumar, B. P. Singh, M. S. Ali, and S. G. Ghosh, “Shadows of black hole surrounded by anisotropic fluid in rastall theory,” Physics of the Dark Universe, vol. 34, p. 100881, 2021.
  * (99) V. Majernik and L. Richterek, “Rastall’s gravity equations and Mach’s Principle,” arXiv e-prints, pp. gr–qc/0610070, Oct. 2006.
  * (100) M. Visser, “The kiselev black hole is neither perfect fluid, nor is it quintessence,” Classical and Quantum Gravity, vol. 37, p. 045001, jan 2020\.
  * (101) M. Azreg-Aïnou, “Generating rotating regular black hole solutions without complexification,” Phys. Rev. D, vol. 90, p. 064041, Sep 2014.
  * (102) M. Azreg-Aïnou, “From static to rotating to conformal static solutions: rotating imperfect fluid wormholes with(out) electric or magnetic field,” The European Physical Journal C, vol. 74, p. 2865, May 2014.
  * (103) V. Perlick and O. Y. Tsupko, “Calculating black hole shadows: Review of analytical studies,” Physics Reports, vol. 947, pp. 1–39, Feb. 2022.
  * (104) R. A. Breuer and J. Ehlers, “Propagation of High-Frequency Electromagnetic Waves Through a Magnetized Plasma in Curved Space-Time. I,” Proceedings of the Royal Society of London Series A, vol. 370, pp. 389–406, Mar. 1980.
  * (105) A. Abdujabbarov, M. Amir, B. Ahmedov, and S. G. Ghosh, “Shadow of rotating regular black holes,” Phys. Rev. D, vol. 93, p. 104004, May 2016.
  * (106) V. Perlick and O. Y. Tsupko, “Light propagation in a plasma on kerr spacetime: Separation of the hamilton-jacobi equation and calculation of the shadow,” Phys. Rev. D, vol. 95, p. 104003, May 2017.
  * (107) J. Bicak and P. Hadrava, “General-relativistic radiative transfer theory in refractive and dispersive media.,” Astron. Astrophys., vol. 44, pp. 389–399, Nov. 1975.
  * (108) R. Kulsrud and A. Loeb, “Dynamics and gravitational interaction of waves in nonuniform media,” Phys. Rev. D, vol. 45, pp. 525–531, Jan 1992.
  * (109) D. O. Muhleman and I. D. Johnston, “Radio propagation in the solar gravitational field,” Phys. Rev. Lett., vol. 17, pp. 455–458, Aug 1966\.
  * (110) V. Perlick, Ray Optics, Fermat’s Principle, and Applications to General Relatively, vol. 61 of Lecture Notes in Physics.  Berlin, Heidelberg: Springer Berlin Heidelberg, 2000.
  * (111) O. Y. Tsupko and G. S. Bisnovatyi-Kogan, “Gravitational lensing in plasma: Relativistic images at homogeneous plasma,” Phys. Rev. D, vol. 87, p. 124009, Jun 2013.
  * (112) H. Yan, “Influence of a plasma on the observational signature of a high-spin kerr black hole,” Phys. Rev. D, vol. 99, p. 084050, Apr 2019.
  * (113) Y. Huang, Y.-P. Dong, and D.-J. Liu, “Revisiting the shadow of a black hole in the presence of a plasma,” International Journal of Modern Physics D, vol. 27, p. 1850114, Sept. 2018.
  * (114) B. Carter, “Global structure of the kerr family of gravitational fields,” Phys. Rev., vol. 174, pp. 1559–1571, Oct 1968.
  * (115) Y. Mino, “Perturbative approach to an orbital evolution around a supermassive black hole,” Physical Review D, vol. 67, p. 084027, Apr 2003.
  * (116) P. F. Byrd and M. D. Friedman, Handbook of Elliptic Integrals for Engineers and Scientists.  Berlin, Heidelberg: Springer Berlin Heidelberg, 1971.
  * (117) S. E. Gralla, D. E. Holz, and R. M. Wald, “Black hole shadows, photon rings, and lensing rings,” Phys. Rev. D, vol. 100, p. 024018, Jul 2019.
  * (118) E. Stoghianidis and D. Tsoubelis, “Polar orbits in the Kerr space-time,” General Relativity and Gravitation, vol. 19, pp. 1235–1249, Dec. 1987.