# Observational signatures of a static $f(R)$ black hole with thin accretion
disk

Mohsen Fathi<EMAIL_ADDRESS>Departamento de Física, Universidad de
Santiago de Chile, Avenida Víctor Jara 3493, Estación Central, 9170124,
Santiago, 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 study, we focus on a static spherically symmetric $f(R)$ black hole
spacetime characterized by a linear dark matter-related parameter. Our
investigation delves into understanding the influence of different assumed
values of this parameter on the observable characteristics of the black hole.
To fulfill this task, we investigate the light deflection angles, which are
inferred from direct analytical calculations of null geodesics. To examine the
black hole’s properties further, we assume an optically thin accretion disk
and explore various emission profiles. Additionally, we investigate the shadow
cast by the illuminated black hole when affected by the disk. Furthermore, we
simulate the brightness of an infalling spherical accretion in the context of
silhouette imaging for the black hole. Our findings indicate that, except for
some specific cases, the observed brightness of the accretion disk
predominantly arises from direct emission, rather than lensing and photon
rings. Moreover, we reveal that the linear dark parameter of the black hole
significantly influences the shadow size and brightness. Our discussion covers
both analytical and numerical approaches, and we utilize ray-tracing methods
to produce accurate visualizations.

keywords: Black holes, $f(R)$ gravity, thin accretion, shadow

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

##  I Introduction and Motivation

Since the foundational concepts of black holes were established by
Schwarzschild Schwarzschild:1916  and Finkelstein Finkelstein:1958 , the
search for identifying these enigmatic objects has been on an exhilarating
journey. From the initial observational evidence gathered for Cygnus X-1 in
1971 webster_cygnus_1972 ; bolton_identification_1972  to the recent shadow
images of M87* Akiyama:2019  and Sgr A* Akiyama:2022  captured by the Event
Horizon Telescope (EHT), the pursuit of understanding black holes has
relentlessly advanced. By comparing theoretical predictions with observed
shadows, valuable insights into the behavior of light in extremely gravitating
systems can be obtained. The EHT results also unveiled the presence of a
magnetic field around M87*, possibly connected to the formation of jets
emanating from the black hole 2021ApJ…910L..12E ; 2021ApJ…910L..13E ;
2021PhRvD.103j4047K . Shadow images can provide information about the
geometric structure near the black hole’s horizon PhysRevD.97.084020  and its
physical characteristics Takahashi_2004 . However, despite these compelling
evidences that support the general theory of relativity, there are limitations
in the cosmological context where general relativity falls short. These
limitations include issues with flat galactic rotation curves, anti-lensing,
the accelerated expansion of the universe, observed anisotropies in the cosmic
microwave background radiation, and the coincidence problem Rubin1980 ;
Massey2010 ; Bolejko2013 ; Riess:1998 ; Perlmutter:1999 ; Astier:2012 .

Many scientists believe that the phenomena described above are attributed to
the dark side of the universe, which remains inadequately explained. For
instance, introducing the cosmological constant term to Einstein’s field
equations as a nonzero vacuum energy can account for the acceleration of the
universe. However, the reason behind the small value of the cosmological
constant remains elusive. To address unresolved cosmological problems, such as
the late-time acceleration of the universe, some researchers turn to modified
theories of gravity to mimic the effects of dark matter and dark energy and
provide an effective time-varying equation of state. In these models, the
Einstein-Hilbert action is generalized or extended to explain the dynamics of
the universe on cosmic, galactic, or astrophysical scales. One of the most
intuitive extensions of general relativity is the $f(R)$ theory, where the
Einstein-Hilbert action is replaced with a generic function $f(R)$
RevModPhys.82.451 ; de_felice_fr_2010 ; CAPOZZIELLO2011167 . As a result,
$f(R)$ theories of gravity have garnered significant interest and have been
thoroughly scrutinized for their consistency Navarro_2007 ; NOJIRI2007238 ;
JE-AN:2011 ; PhysRevD.78.063504 ; Capozziello_2008 ; cognola_class_2008 ;
nojiri_modified_2008  (see also the reviews in Refs. nojiri_unified_2011 ;
nojiri_modified_2017 ). Within the context of $f(R)$ gravity, researchers are
also interested in black hole solutions. One primary solution derived from the
Palatini formalism of $f(R)$ theory is the Schwarzschild-(anti-)de Sitter
metric, which includes an effective cosmological constant. However, this
solution faces incompatibilities with primary general relativistic tests since
the cosmological constant appears to have no significant role in solar system
scales KAGRAMANOVA2006465 . Nonetheless, this problem can be circumvented by
appropriately manipulating the action, rendering the cosmological constant
effective on cosmological scales while negligible on solar system scales
starobinsky_disappearing_2007 ; PhysRevD.76.064004 . To address these issues
properly, a suitable $f(R)$ action model has been proposed in Refs. refId0 ;
refId0_1 ; RAHVAR:2008 , which is consistent with both galactic and
cosmological scales. In Ref. Saffari:2008 , this model was further elaborated
using a generic function in the gravitational action to ensure its
compatibility with solar system tests, galactic rotation curves, and the late-
time acceleration of the universe. In this context, the authors also propose a
static spherically symmetric black hole solution, which is of interest in this
paper, particularly concerning light propagation in its geometry and its
shadow.

Indeed, the theoretical constraining of black hole shadows based on
observational data has been a subject of special interest to scientists,
leading to numerous publications dedicated to this area (see, for example,
Refs. Abdujabbarov:2012bn ; Yumoto:2012kz ; Atamurotov:2013sca ;
Zakharov:2014lqa ; Papnoi:2014aaa ; Johannsen:2015hib ; Johannsen:2015mdd ;
Moffat:2015kva ; Giddings:2016btb ; Cunha:2016bjh ; Tsukamoto:2017fxq ;
Hennigar:2018hza ; Cunha:2018acu ; Allahyari:2019jqz ; Abdikamalov:2019ztb ;
Ovgun:2019jdo ; Shaikh:2019fpu ; bambi_testing_2019 ; vagnozzi_hunting_2019 ;
Kumar:2020hgm ; Li2020 ; Ovgun:2020gjz ; khodadi_black_2020 ; Zhong:2021mty ;
Zuluaga:2021vjc ; Stashko:2021lad ; Rahaman:2021kge ; Ovgun:2021ttv ;
Pantig:2022whj ; Bisnovatyi-Kogan:2022ujt ; Kazempour:2022asl ; Pantig:2022ely
; roy_superradiance_2022 ; chen_superradiant_2022 ; vagnozzi_horizon-
scale_2022 ; omwoyo_remarks_2022 ). However, the recent silhouette imaging by
the EHT has added even greater importance to the need for reliable methods to
visualize black holes with accretion disks as their illumination sources. This
interest was initially ignited by Luminet in 1979 when he calculated the
radiation emitted from a thin accretion disk surrounding a Schwarzschild black
hole and proposed a ray-traced image of the disk Luminet:1979nyg . Generally,
this type of accretion is based on the Shakura-Sunyaev Shakura:1973 , Novikov-
Thorne Novikov:1973 , and Page-Thorne page_disk-accretion_1974  models, where
the disk is assumed to be thin, geometrically and optically. These
assumptions, along with the growing interest in black hole imaging, led to the
development of a new method for simulating higher-order light rings for black
holes with thin accretion disks, proposed in Ref. Gralla:2019 . Since then,
this method has been applied in various publications (e.g., Refs.
Guerrero_2021 ; li_observational_2021 ; okyay_nonlinear_2022 ;
PhysRevD.105.084057 ; hu_observational_2022 ; PhysRevD.106.044070 ; Guo_2022 ;
WANG2022116026 ; PhysRevD.105.064031 ; uniyal_probing_2023 ;
uniyal_nonlinearly_2023 ), and it also plays a significant role in our paper.

We structure our discussion as follows: In Sec. II, we provide a concise
overview of the $f(R)$ black hole solution and introduce its cosmological
parameters. Moving on to Sec. III, we initiate our investigation by studying
the causal structure of the spacetime. We then apply a Lagrangian formalism to
derive the equations of motion for massless particles (light rays). By
calculating the critical impact parameter of photon trajectories, we determine
the points at which the orbits become unstable. As a result, we are able to
calculate the theoretical size of the black hole shadow, which directly
correlates with the critical impact parameter. In the same section, we proceed
to find the turning points of light ray trajectories as they approach the
black hole. We obtain exact analytical solutions for the angular equation of
motion for deflecting and critical trajectories. In fact, gravitational
lensing serves as a remarkable tool in examining black hole solutions in the
strong-field regime Virbhadra:1998 ; Virbhadra:2000 ; Bozza:2001 ;
Virbhadra:2002 ; Bozza:2002 ; Hasse:2002 ; Virbhadra:2008 ; Virbhadra:2009 ;
He:2020 ; Adler:2022 ; Virbhadra:2022 . Additionally, weak lensing is of great
importance to astrophysicists and cosmologists, allowing them to estimate
matter distribution profiles within galaxies and other observable regions of
the universe Hoekstra:2003 ; Sheldon:2004 ; Mandelbaum:2006 ; Gavazzi:2007 ;
Parker:2007 , thus providing insights into dark matter and dark energy
properties. Hence, we use the above analytical solutions as instruments in
finding the lens equation and the deflection angles, which we calculate
analytically using the Weierstrassian elliptic function. In Sec. IV, we
construct a thin accretion disk for the black hole based on the Novikov-Thorne
model. We calculate the dynamical characteristics of accreting particles in
stable orbits and derive the radial profiles of the disk’s radiation flux and
temperature. Applying the method from Ref. Gralla:2019 , we visualize the
light rings and accretion disk of the $f(R)$ black hole for three different
disk emission profiles. Additionally, we calculate the thickness of the rings,
which is inferred from the observed effective intensity profiles. Finally, we
assume a spherically symmetric infalling accretion for the black hole and
calculate the observed disk emission, concluding this section by simulating
the black hole shadow under this condition. Our paper concludes with Sec. V,
where we summarize our findings and discuss their implications. Throughout
this work, we use the signature convention $(-++\,+)$, and primes on functions
denote differentiation with respect to the radial coordinate. We apply the
geometrized unit system, where $G=c=1$.

##  II A particular model of $f(R)$ gravity and its black hole solution

The gravitational action of the theory can be written in its most generic form
as

$\mathcal{S}=\frac{1}{2\kappa}\int\mathrm{d}x^{4}\sqrt{-g}~{}f(R)+\mathcal{S}_{m},$
(1)

in which, $\kappa$ is a coupling constant, $f(R)$ is a function of the Ricci
scalar of the spacetime with the metric determinant $g$, and $\mathcal{S}_{m}$
is a matter field action. Accordingly, the field equations are derived as

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

by varying the action $\mathcal{S}$ with respect to the metric, where
$F(R)=\frac{\mathrm{d}f(R)}{\mathrm{d}R}$,
$\square=\nabla_{\lambda}\nabla^{\lambda}$ and $T_{\mu\nu}$ is the energy-
momentum tensor. Adopting a generic spherically symmetric line element

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

in the usual Schwarzschild coordinates $x^{\mu}=(t,r,\theta,\phi)$, the field
equation (2) can be recast as Saffari:2008

$\displaystyle
2F(r)\frac{X^{\prime}(r)}{X(r)}+rF^{\prime}(r)\frac{X^{\prime}(r)}{X(r)}-2rF^{\prime\prime}(r)=0,$
(4) $\displaystyle
B^{\prime\prime}(r)+\left[\frac{F^{\prime}(r)}{F(r)}-\frac{1}{2}\frac{X^{\prime}(r)}{X(r)}\right]B^{\prime}(r)-\frac{2}{r}\left[\frac{F^{\prime}(r)}{F(r)}-\frac{1}{2}\frac{X^{\prime}(r)}{X(r)}\right]B(r)-\frac{2}{r^{2}}B(r)+\frac{2}{r^{2}}X(r)=0,$
(5)

for the vacuum space (i.e. $T_{\mu\nu}=0$), in which $X(r)=B(r)A(r)$, and
based on the $r$-dependence of the Ricci scalar, we have taken into account
$F=F(R)\equiv F(r)$. Obviously, general relativity is recovered for $F(R)=1$,
for which Eq. (4) results in $X(r)=1$ (or $B(r)=A(r)^{-1}$), and from Eq. (5),
the Schwarzschild solution is obtained. We adopt the ansatz
$F(r)=(1+r/d)^{-\alpha}$ Saffari:2008 , where $\alpha$ and $d$ are the free
parameters of the action. Here $\alpha$ is dimensionless whereas
$[d]=\mathrm{m}$ and is of galactic size. Applying this ansatz, Eq. (4)
results in the solution

$X(r)=X_{0}\left(1+\frac{r}{d}\right)^{-2(1+\alpha)}\left(1+\frac{2-\alpha}{2}\frac{r}{d}\right)^{{4(1+\alpha)}/{(2-\alpha)}},$
(6)

where $X_{0}$ represents an integration constant. When $\alpha=0$, we obtain
$X_{0}=1$, which corresponds to the Schwarzschild solution. Solving Eq. (5)
allows us to determine $B(r)$. By substituting the expression from Eq. (6)
into Eq. (5), and considering terms up to the first-order of the free
parameters, the lapse function can be obtained as Saffari:2008

$B(r)=1-\frac{2M}{r}+\beta r-\frac{1}{3}\Lambda r^{2},$ (7)

in which $\beta={\alpha}/{d}$. Here, $\Lambda$ represents the cosmological
constant having the value $|\Lambda|\leq 10^{-52}~{}\mathrm{m}^{-2}$
Carmeli:2001 111Unless otherwise specified, this value of $\Lambda$ is assumed
throughout the subsequent sections of the paper.. It is worth noting that the
above expression resembles the Mannheim-Kazanas vacuum solution in fourth-
order Weyl conformal gravity mannheim_exact_1989 . In that solution, the
linear term $\beta r$ acts as an additional potential that compensates for the
flat galactic rotation curves. Similarly, the model presented in Eq. (7)
suggests that small values of $\alpha$ can yield flat galactic rotation curves
for galaxies in a similar manner. The same solution has been further explored
in various aspects, including investigations on the quasinormal modes
Aragon:2020xtm ; Malekmakan:2022kfd , dynamics of the accretion disk
Soroushfar:2020kgb , thermodynamic geometry Upadhyay:2018bqy , and particle
collision in the black hole’s exterior Majeed:2017txa . The lapse function (7)
can be seen as a linear combination of solutions from two distinct $f(R)$
gravity models, each pertinent to a different length scale. In the solar
system scale ($r\ll d$), the model is derived as Saffari:2008

$f(R)=R+R_{0}\ln\frac{R}{R_{c}},$ (8)

where $R_{0}={6\alpha^{2}}/{d^{2}}$, and $R_{c}$ is an integration constant.
The solution to the above model covers the first three terms in the lapse
function (7). On the other hand, in the galactic scales ($r>d$), the $f(R)$
theory is based on the model Saffari:2008

$f(R)=R_{c}^{-{\alpha}/{2}}\left(R+\Lambda\right)^{1+{\alpha}/{2}}.$ (9)

For cosmological scales ($\alpha\ll 1$), this model reduces to
$f(R)=R+\Lambda$, which is equivalent to the Einstein-Hilbert action with a
cosmological constant. Therefore, in the context of the exact solution, this
model contributes the final term in the lapse function (7). It is worth noting
that in Ref. Saffari:2008 , the generic model

$f(R)=R+\Lambda+\frac{R+\Lambda}{{R}/{R_{0}}+{2}/{\alpha}}\ln\frac{R+\Lambda}{R_{c}},$
(10)

has been proposed, which aligns with both of the above models under different
scale criteria. In the regime of strong curvature, characterized by
$R\gg\Lambda$ and ${R}/{R_{0}}\gg{2}/{\alpha}$, we obtain the model (8).
Conversely, in cosmological scales, where $R\simeq R_{0}\simeq\Lambda$ and
$\alpha\ll 1$, it simplifies to $f(R)=R+\Lambda$. Consequently, the small-
valued parameter $\beta$ in Eq. (7) can effectively capture both large and
small-scale phenomena in the universe. In the following section, we will
assign presumed values to this model parameter within a favorable range. These
values will serve as the basis for our subsequent analysis and simulations in
the discussion.

We begin by studying the black hole exterior geometry and its causal
structure. For the sake of convenience in the calculations and demonstrations,
we adimensionalize the parameters by introducing the quantities

$\tilde{r}\rightarrow\frac{r}{M},\quad\tilde{\beta}\rightarrow\beta
M,\quad\tilde{\Lambda}\rightarrow\frac{1}{3}\Lambda M^{2}.$ (11)

In the forthcoming sections, however, we remove the ”tilde” overscript from
the dimensionless parameters in Eq. (11), which is equivalent to letting
$M=1$.

##  III Propagation of light and unstable photon orbits

The spacetime’s causal structure, characterized by the lapse function (7), can
be understood by studying the hypersurfaces where $B(r)=0$, corresponding to
the black hole horizons, resulting in a cubic equation with the solutions

$\displaystyle
r_{1}=\frac{\beta}{3\Lambda}-\frac{4}{\Lambda}\sqrt{\frac{g_{2}}{3}}\cos\left(\frac{1}{3}\arccos\left(\frac{3g_{3}}{g_{2}}\sqrt{\frac{3}{g_{2}}}\right)-\frac{4\pi}{3}\right),$
(12) $\displaystyle
r_{2}=\frac{\beta}{3\Lambda}-\frac{4}{\Lambda}\sqrt{\frac{g_{2}}{3}}\cos\left(\frac{1}{3}\arccos\left(\frac{3g_{3}}{g_{2}}\sqrt{\frac{3}{g_{2}}}\right)-\frac{2\pi}{3}\right),$
(13) $\displaystyle
r_{3}=\frac{\beta}{3\Lambda}-\frac{4}{\Lambda}\sqrt{\frac{g_{2}}{3}}\cos\left(\frac{1}{3}\arccos\left(\frac{3g_{3}}{g_{2}}\sqrt{\frac{3}{g_{2}}}\right)\right),$
(14)

where

$\displaystyle g_{2}=\frac{1}{12}\left(\beta^{2}+3\Lambda\right),$ (15a)
$\displaystyle
g_{3}=-\frac{1}{16}\left(\frac{2\beta^{3}}{27}+\frac{\beta\Lambda}{3}-2\Lambda^{2}\right).$
(15b)

The existence of real values for the above radii, however, depends on the sign
of the polynomial’s discriminant, i.e.
$\Delta=g_{2}^{3}-27g_{3}^{2}=\frac{\Lambda^{2}}{256}[\beta^{2}(1+8\beta)+4(1+9\beta)-108\Lambda^{2}]$,
which is always of positive values for $\beta,\Lambda\ll 1$. So, we infer that
all the solutions (12)–(14) are real-valued. It is straightforward to verify
that $r_{1}>r_{2}>0$ and $r_{3}<0$. Hence we identify $r_{++}=r_{1}$, as the
cosmological horizon of the black hole where the infinite blueshift happens,
and $r_{+}=r_{2}$, as its event horizon, where the infinite redshift happens.
This way, the lapse function can be rewritten as

$B(r)=\frac{\Lambda}{r}\left(r_{++}-r\right)\left(r-r_{+}\right)\left(r-r_{3}\right).$
(16)

###  III.1 Motion of mass-less particles

The motion of test particles can be described by the Lagrangian

$\displaystyle 2\mathscr{L}$ $\displaystyle=$ $\displaystyle
g_{\mu\nu}\dot{x}^{\mu}\dot{x}^{\nu}$ (17) $\displaystyle=$
$\displaystyle-B(r)\dot{t}^{2}+\frac{\dot{r}^{2}}{B(r)}+r^{2}\dot{\theta}^{2}+r^{2}\sin^{2}\theta\dot{\phi}^{2},$

where ”dot” stands for differentiation with respect to the affine parameter
$\tau$ of the geodesic curves. Enjoying the spherical symmetry of the
spacetime, we confine the motion of particles to the equatorial plane (i.e.
$\theta={\pi}/{2}$) without loss of generality. One can then define the
conjugate momenta

$\Pi_{\mu}=\frac{\partial\mathscr{L}}{\partial{\dot{x}^{\mu}}},$ (18)

which provides the two constants of motion

$\displaystyle\Pi_{t}=-B(r)\dot{t}=-E,$ (19a)
$\displaystyle\Pi_{\phi}=r^{2}\dot{\phi}=L,$ (19b)

in accordance with the Killing symmetries of the spacetime, and we name them,
respectively, as the energy and the angular momentum of the test particles.
These two quantities allow us to define the impact parameter $b\equiv{L}/{E}$.
This parameter corresponds to the vertical distance between the tangent to the
geodesic curves and the line passing the black hole singularity and is of
importance in the identification of possible photon trajectories. In fact, the
motion of photons can be described by the equation $\mathscr{L}=0$ which
characterizes the null geodesics. Thus, by means of Eq. (17), the equations of
motion are obtained as

$\displaystyle\dot{r}^{2}=E^{2}-V(r),$ (20)
$\displaystyle\left(\frac{\mathrm{d}r}{\mathrm{d}\phi}\right)^{2}=\frac{r^{4}}{b^{2}}\left[1-\frac{V(r)}{E^{2}}\right],$
(21)

in which

$V(r)=L^{2}\frac{B(r)}{r^{2}},$ (22)

represents the effective potential for photons. Then, the turning points
$r_{t}$ in the orbits correspond to $\dot{r}=0$, which is encountered when
$V(r_{t})=E^{2}$. This potential has a maximum at

$r_{p}=\frac{\beta_{0}-1}{\beta},$ (23)

given that $\beta_{0}=\sqrt{1+6\beta}$, which is where the photon orbits
become unstable. Hence, this maximum corresponds to the radius of the photon
sphere. As it is observed, this radius is independent of $\Lambda$, and
decreases as $\beta$ increases. One can also verify that

$\lim_{\beta\rightarrow 0}r_{p}=3,$ (24)

which is the radius of unstable photon orbits for a Schwarzschild-de Sitter
black hole. It is straightforward to calculate the critical value of the
impact parameter which is obtained as

$b_{p}=\frac{\sqrt{3}\left(\beta_{0}-1\right)}{\sqrt{\beta^{2}\left(2\beta_{0}-1\right)-18\beta\Lambda+6\left(\beta_{0}-1\right)\Lambda}},$
(25)

which is crucial in determining the size of the black hole’s apparent shape in
the observer’s sky. Historically, various definitions have been proposed in
the literature to describe the visual appearance of black holes. These include
the escape cone by Synge Synge:1966 , and the cone of gravitational radiation
capture by Zeldovich and Novikov zeldovich_relativistic_1966 . Bardeen,
Chandrasekhar, and Luminet popularized the terms optical appearance of black
holes and black hole image to describe the visual characteristics of black
holes Bardeen:1973b ; Chandrasekhar:1998 ; Luminet:1979nyg ;
luminet_seeing_2018 . These terms are widely used in the field. In addition,
it is worth noting that the term ”black hole shadow,” which is also utilized
in this study, was introduced by Falcke, Melia, and Agol Falcke_2000 . This
term refers to the dark region observed within the apparent boundary of the
black hole. The apparent boundary is occasionally referred to as the photon
ring johannsen_testing_2010 ; johnson_universal_2020 , although in this study,
the term holds a broader significance (see Ref. perlick_calculating_2022  for
further historical insights on this matter). In this context, the shadow of
black holes is created by the presence of light rings, which represent the
lensed images of their luminous background. The radius of these rings can be
determined by the critical impact parameter (i.e. $R_{\mathrm{sh}}=b_{p}$).
Specifically, for the $f(R)$ black hole being examined, the expression for
$R_{\mathrm{sh}}$ is given by Eq. (25). This way, the theoretical shadow
diameter for this black hole is given by
$d_{\mathrm{sh}}^{\mathrm{theo}}=2R_{\mathrm{sh}}=2b_{p}$. In Fig. 1(a), the
$\beta$-profile of this theoretical diameter has been plotted. It is worth
noting that one can calculate the diameter of the observed shadows in the
recent EHT images of M87* and Sgr A*. To do so, let us consider the relation
Bambi:2019tjh

$d_{\text{sh}}=\frac{D\theta_{*}}{\gamma},$ (26)

which calculates the shadow diameter as observed by an observer positioned at
a distance $D$ (in parsecs) from the black hole, where $\gamma$ is the mass
ratio of the black hole and the Sun, which is $\gamma=(6.5\pm 0.90)\times
10^{9}$ for M87* at the distance $D=16.8\,\mathrm{Mpc}$ Akiyama:2019 , and is
$\gamma=(4.3\pm 0.013)\times 10^{6}$ for Sgr A* at $D=8.127\,\mathrm{kpc}$
Akiyama:2022 . In Eq. (26), $\theta_{*}$ is the angular diameter of the
shadow, which has been measured as $\theta_{*}=42\pm 3\,\mathrm{\mu as}$ for
M87*, and $\theta_{*}=48.7\pm 7\,\mathrm{\mu as}$ for Sgr A*. This way, one
can calculate the shadow diameters as
$d_{\mathrm{sh}}^{\mathrm{M87^{*}}}=11\pm 1.5$ and
$d_{\mathrm{sh}}^{\mathrm{SgrA^{*}}}=9.5\pm 1.4$. In Fig. 1(b), these values
are displayed within the $1\sigma$ uncertainties for both black holes.
However, these observed diameters do not correspond to the calculated
theoretical shadow diameters $d_{\mathrm{sh}}^{\mathrm{theo}}$ in the same
physical sense. The theoretical diameter relates to the photon ring and, in
general, the photon ring’s brightness is predominantly influenced by direct
disk emission (further elaborated in Subsec. IV.1). Additionally,
uncertainties stemming from the physical characteristics of the accretion disk
and its emission profiles introduce further complexities. Consequently,
attempting to constrain black hole parameters by comparing the theoretical
shadow diameter to observations from the EHT is not viable. Therefore, such
analysis does not offer a means to constrain the model parameter $\beta$ for
the $f(R)$ black hole.

(a)  (b)

Figure 1: The diagrams show (a) the $\beta$-profile of
$d_{\mathrm{sh}}^{\mathrm{theo}}$ and (b) the observed shadow diameters
$d_{\mathrm{sh}}$ for M87* (green region) and Sgr A* (red region), each with
$1\sigma$ uncertainties. The shaded region in panel (a) represents the assumed
range for the $\beta$-parameter, limiting it within the domain
$0\leq\beta<0.05$.

Instead, for the numerical studies in the rest of the paper, we consider the
values of the $\beta$-parameter within a specific range, marked in the shaded
region of the $\beta$-profile of the theoretical shadow diameter (as seen in
Fig. 1(a)). This segment confines the parameter within the range
$0\leq\beta<0.05$, which we adhere to in the subsequent discussions in this
paper.

Hence, it is now possible to visualize the behavior of the effective potential
(22), which is shown in Fig. 2.

(a)  (b)

Figure 2: In panel (a), the radial profile of the effective potential is shown
for various values of $\beta$, within the assumed range. In panel (b), by
adopting $\beta=0.022$, a typical effective potential has been shown together
with the turning points and their corresponding impact parameters.

In the left panel, five radial profiles of the effective potential have been
plotted which correspond to different values of the $\beta$-parameter within
the assumed range. As it is inferred from the diagram, the potential possesses
one maximum, at which, unstable orbits can occur. This is shown in more detail
in the right panel of Fig. 2, where the orbits are categorized in accordance
with the values of the impact parameter $b$. When $b>b_{p}$, the photons may
approach from either of the turning points $r_{d}$ (where they are recessively
deflected by the black hole) or $r_{f}$ (where they are deflected towards the
event horizon). In fact, the turning points can be obtained analytically for
the spacetime of the $f(R)$ black hole, by solving the equation
$({\mathrm{d}r}/{\mathrm{d}\phi})^{2}=0$. Applying the Eqs. (21) and (22),
this results in

$\left(\frac{\mathrm{d}r}{\mathrm{d}\phi}\right)^{2}=\mathcal{P}_{4}(r)\equiv
r\left(\frac{r^{3}}{\lambda^{2}}-\beta r^{2}-r+2\right)=0,$ (27)

which beside the trivial solution at $r=0$, has one negative root and the two
positive roots $r_{d}=x_{d}^{-1}$ and $r_{f}=x_{f}^{-1}$, where

$\displaystyle
x_{f}=\frac{1}{6}-2\sqrt{\frac{\bar{g}_{2}}{3}}\sin\left(\frac{1}{3}\arcsin\left(\frac{3\bar{g}_{3}}{\bar{g}_{2}}\sqrt{\frac{3}{\bar{g}_{2}}}\right)-\frac{2\pi}{3}\right),$
(28) $\displaystyle
x_{d}=\frac{1}{6}-2\sqrt{\frac{\bar{g}_{2}}{3}}\sin\left(\frac{1}{3}\arcsin\left(\frac{3\bar{g}_{3}}{\bar{g}_{2}}\sqrt{\frac{3}{\bar{g}_{2}}}\right)\right),$
(29)

with ${\lambda^{-2}}={b^{-2}}+\Lambda$, and

$\displaystyle\bar{g}_{2}=\frac{1}{4}\left(\frac{1}{3}+2\beta\right),$ (30a)
$\displaystyle\bar{g}_{3}=-\frac{1}{4}\left(\frac{1}{\lambda^{2}}-\frac{\beta}{6}-\frac{1}{54}\right).$
(30b)

For the case of $b=b_{p}$, the photons encounter the turning point $r_{p}$
(see the right panel of Fig. 2), which has been determined in Eq. (23). At
this stage, the photons travel on unstable (or critical) orbits, which form
the black hole shadow. In Table. 1, the turning points have been given for
different values of $\beta$.

$\beta$ | $0.0$ | 0.011 | 0.022 | 0.031 | 0.041
---|---|---|---|---|---
$r_{+}$ | 2.0001 | 1.9563 | 1.9196 | 1.8880 | 1.8583
$r_{p}$ | 3.0000 | 2.9503 | 2.9077 | 2.8704 | 2.8350
$b_{p}$ | 5.1968 | 4.9468 | 4.7445 | 4.5762 | 4.4228
$r_{d}$ | 6.3692 | 6.8091 | 7.2156 | 7.5957 | 7.9813
$r_{f}$ | 2.1797 | 2.1076 | 2.0544 | 2.0103 | 1.9699
Table 1: The turning points of photonic trajectories together with their
corresponding values of $b_{p}$.

As it is inferred from the table, increase in the $\beta$-parameter leads to a
smaller black hole (decrease in $r_{+}$), a wider effective potential
(increase in the distance between $r_{d}$ and $r_{f}$), and a lower potential
maximum. So unstable orbits are less likely to happen for larger $\beta$
(decrease in $b_{p}$), and the black hole shadow decreases in size. In such
cases, the light rays detected by a distant observer are mostly dominated by
the direct emission, which is a simply lensed image of the black hole’s
emitting disk or its luminous background. This will be discussed in more
details in the forthcoming sections.

Now that the turning points have been obtained and analyzed, we proceed with
the determination of the exact solutions for the aforementioned possible
photon orbits around the $f(R)$ black hole. In fact, the null geodesics for
this black hole have been studied in their most general form in Ref.
soroushfar_analytical_2015 . In what follows, however, we base our exact
solutions on the analytically known turning points, which helps us
investigate, separately, the deflecting and critical trajectories that are of
importance for the purpose of this paper.

####  III.1.1 Deflecting trajectories

As mentioned above, photonic trajectories become deflected at the turning
points $r_{d}$ and $r_{f}$ which leads to different fates for the photons.
Accordingly, the possible deflecting trajectories can be ramified as the orbit
of the first kind (OFK) at $r_{d}$, and the orbit of the second kind (OSK) at
$r_{f}$. Since both of these turning points have been identified analytically,
hence, by applying the change of variable $z\doteq{r_{i}}/{r}$
($r_{i}=r_{d},r_{f}$), one can rewrite the differential equation (27) as

$\left(\frac{\mathrm{d}z}{\mathrm{d}\phi}\right)^{2}=\mathcal{P}_{3}(z)\equiv\frac{2}{r_{i}}z^{3}-z^{2}-r_{i}\beta
z+\frac{r_{i}^{2}}{\lambda^{2}}.$ (31)

A further change of variable $u\doteq\frac{1}{2}({z}/{r_{i}}-{1}/{6})$
provides us with the Weierstrassian differential equation

$\left(\frac{\mathrm{d}u}{\mathrm{d}\phi}\right)^{2}=\tilde{\mathcal{P}}_{3}(u)\equiv
4u^{3}-\tilde{g}_{2}u-\tilde{g}_{3},$ (32)

in which

$\displaystyle\tilde{g}_{2}=\frac{1}{12}(1+6\beta),$ (33a)
$\displaystyle\tilde{g}_{3}=-\frac{1}{216}\left(\frac{54}{\lambda^{2}}-9\beta-1\right),$
(33b)

are known as the Weierstrass invariants. This leads to the integrals

$\displaystyle\phi-\phi_{0}=\int_{u_{d}}^{u}\frac{\mathrm{d}u^{\prime}}{\sqrt{\tilde{\mathcal{P}}_{3}(u^{\prime})}}~{}~{}(\mathrm{with}~{}u_{d}<u),$
(34)
$\displaystyle\phi-\phi_{0}=\int_{u}^{u_{f}}\frac{\mathrm{d}u^{\prime}}{\sqrt{\tilde{\mathcal{P}}_{3}(u^{\prime})}}~{}~{}(\mathrm{with}~{}u_{f}>u),$
(35)

respectively, for the OFK and OSK, in which $\phi_{0}$ is the initial azimuth
angle, and $u_{i}=\frac{1}{2}({1}/{r_{i}}-{1}/{6})$. Taking into account the
applied changes of variables, the above integrals yield

$r(\phi)=\frac{6}{1+12\wp\left(\omega_{d}-(\phi-\phi_{0})\right)},$ (36)

for the OFK and

$r(\phi)=\frac{6}{1+12\wp\left(\omega_{f}+(\phi-\phi_{0})\right)},$ (37)

for the OSK, where $\wp(x)\equiv\wp(x;\tilde{g}_{2},\tilde{g}_{3})$ is the
$\wp$-Weierstrassian elliptic function byrd_handbook_1971 , and we have
defined

$\omega_{i}=\wp^{-1}\left(\frac{1}{2r_{i}}-\frac{1}{12}\right).$ (38)

####  III.1.2 Deflection angle

The OFK is in fact related to the gravitational lensing that is caused by the
black hole and is, in part, responsible for the formation of the black hole
shadow. Hence, by having at hand the integral equation (34), one can calculate
the deflection angle $\hat{\Theta}$ that an observer $\mathbb{O}$ at the
radial position $r_{\mathbb{O}}$ from the black hole (the lens) measures.
Accordingly, we have misner_gravitation_2017

$\displaystyle\hat{\Theta}$ $\displaystyle=$ $\displaystyle 2\Delta\phi-\pi$
(39) $\displaystyle=$ $\displaystyle
2\int_{u_{\mathbb{O}}}^{u_{d}}\frac{\mathrm{d}u^{\prime}}{\sqrt{\tilde{\mathcal{P}}_{3}(u)}}-\pi=2\left[\wp^{-1}(u_{\mathbb{O}})-\wp^{-1}(u_{d})\right]-\pi,$

with the Weierstrass invariants given in Eqs. (33). Using the analytical
expression for $r_{d}$, we have plotted the behavior of $\hat{\Theta}$ in
terms of the impact parameter $b$ in Fig. 3.

Figure 3: The plot of the deflection angle $\hat{\Theta}$ (in $\mathrm{\mu
as}$) versus the changes in the impact parameter $b$, plotted for the
different values of the $\beta$-parameter, and $r_{\mathbb{O}}=10^{5}$.

As it can be inferred from the diagram, there is no significant sensitivity in
the behavior of the deflection angle, given the small changes in the
$\beta$-parameter. However, in general, $\hat{\Theta}$ decreases for a given
$b$, when the $\beta$-parameter increases from 0 to its maximum value. In this
sense, the Schwarzschild-de Sitter spacetime causes the highest deflection
angle. Naturally, by approaching the black hole (smaller $b$) the deflection
angle increases until it diverges at $b_{p}$ for each of the cases. So, strong
lensing occurs in the near-horizon regions whereas by receding the black hole,
the light deflection process will correspond to weak lensing.

####  III.1.3 Critical trajectories

Similar to the deflecting trajectories, unstable orbits can also lead to
different fates, which we name after as the critical orbits of the first kind
(COFK) which happen when photons approach $r_{p}$ from an initial distance
$r_{p}<r_{\mathrm{in}}<r_{++}$, and the critical orbits of the second kind
(COSK) which correspond to photons approaching $r_{p}$ from
$r_{+}<r_{\mathrm{in}}<r_{p}$. When
${\lambda^{-2}}\rightarrow{\lambda_{p}^{-2}}={b_{p}^{-2}}+\Lambda$, the point
$r=r_{p}$ is a double root of $\mathcal{P}_{4}(r)$ in Eq. (27). The
differential equation of motion can then be factorized as

$\frac{\mathrm{d}r}{\mathrm{d}\phi}\equiv\mathcal{P}_{4}^{p}(r)=\left|r-r_{p}\right|\sqrt{\frac{r^{2}}{\lambda_{p}^{2}}+\left(\frac{r_{p}}{\lambda_{p}^{2}}+\chi_{1}\right)r+\frac{r_{p}^{2}}{\lambda_{p}^{2}}+\chi_{1}r_{p}+\chi_{0}},$
(40)

by means of the method of synthetic division, in which

$\displaystyle\chi_{0}=r_{p}\left(\frac{r_{p}}{\lambda_{p}^{2}}-\beta\right)-1,$
(41a) $\displaystyle\chi_{1}=\frac{r_{p}}{\lambda_{p}^{2}}-\beta.$ (41b)

One can therefore obtain the exact solutions to Eq. (40), by means of direct
integration and applying the inversion. This yields the two solutions
$r_{\mathrm{I}}(\phi)$ for the COFK and $r_{\mathrm{II}}(\phi)$ for the COSK,
which are given as

$r_{\begin{subarray}{c}\mathrm{I}\\\
\mathrm{II}\end{subarray}}(\phi)=\frac{1}{\mathcal{A}^{2}+8r_{p}^{2}+(\mathcal{A}^{2}-4\mathcal{B})\cosh\Phi}\Big{[}2r_{p}\left(\mathcal{A}^{2}-4\mathcal{B}\right)\cosh^{2}\Phi\pm\left[r_{p}(r_{p}+\mathcal{A})+\mathcal{B}\right]\\\
\times\sqrt{(\mathcal{A}^{2}-4\mathcal{B})\operatorname{sech}^{2}\Phi\tanh^{2}\Phi}\big{(}\cosh(2\Phi)\mp
2\mathcal{A}\big{)}\Big{]},$ (42)

where
$\Phi=\lambda_{p}^{-1}(\phi-\phi_{0})\sqrt{r_{p}(\mathcal{A}+r_{p})+\mathcal{B}}$,
and

$\displaystyle\mathcal{A}=r_{p}+\chi_{1}\lambda_{p}^{2},$ (43a)
$\displaystyle\mathcal{B}=r_{p}^{2}+\lambda_{p}^{2}\left(\chi_{0}+\chi_{1}r_{p}\right).$
(43b)

In Fig. 4, some examples of the possible photon orbits have been shown.

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

Figure 4: Examples of (a) OFK, (b) OSK, (c) COFK and (d) COSK, plotted for
$\beta=0.022.$ In panels (a) and (b), the red circle at the center denotes
$r_{+}$, while the blue and purple dashed circles correspond to $r_{f}$ and
$r_{d}$. In panels (c) and (d), the exterior dashed circles indicate $r_{p}$.

As can be inferred from the diagrams, the boundary of the black hole shadow is
based on these four types of orbits. Together, these photon orbits are able to
produce the bright ring surrounding the dark shadow in the observer’s sky,
which is produced as a result of the strong gravitational lensing in the near-
horizon regions.

As previously mentioned, the lensing phenomenon and the existence of critical
photon orbits play a crucial role in shaping the black hole shadow and forming
the light rings. However, the formation of these features is highly dependent
on how the black hole is illuminated. In the case of astrophysical black
holes, illumination is often provided by an accretion disk, which is also
taken into account in this study. In the following section, we begin by
developing a thin accretion model and subsequently delve into both analytical
and numerical examinations of the emission process from this disk, which
contributes to the formation of the rings.

##  IV Thin accretion disk model and emission from the black hole

In this section, we study the observational signatures of the black hole in
the case of the existence of a thin accretion disk. We assume that the
accretion process is explained by a generalized version of the well-known
Shakura-Sunyaev model Shakura:1973 , proposed by Novikov and Thorne in Ref.
Novikov:1973 . To proceed with applying this model, let us return to the
Lagrangian (17), which is now specified as $2\mathscr{L}=-1$, for massive
particles of energy $\mathcal{E}$ and angular momentum $\mathcal{L}$ that
constitute the accretion disk. This way, one can rewrite the equations of
motion (20) and (21) as

$\displaystyle\dot{r}=\mathcal{E}^{2}-\mathcal{V}(r),$ (44)
$\displaystyle\left(\frac{\mathrm{d}r}{\mathrm{d}\phi}\right)^{2}=\frac{\mathscr{P}_{6}(r)}{\mathcal{L}^{2}},$
(45)

in which $\mathscr{P}_{6}(r)=r\left[\Lambda r^{5}-\beta
r^{4}-(1-\mathcal{E}^{2}-\mathcal{L}^{2}\Lambda)r^{3}+(2-\mathcal{L}^{2}\beta)r^{2}-\mathcal{L}^{2}r-2\mathcal{L}^{2}\right]$,
and

$\mathcal{V}(r)=B(r)\left(1+\frac{\mathcal{L}^{2}}{r^{2}}\right),$ (46)

is the effective potential for massive particles orbiting the black hole in
the equatorial plan. The left panel of Fig. 5 shows a typical radial profile
of $\mathcal{V}(r)$ which has been plotted for the adopted values of $\beta$.

(a)  (b)

Figure 5: (a) The radial profile of the effective potential $\mathcal{V}(r)$
plotted for different values of the $\beta$-parameter and $\mathcal{L}=10$.
(b) The position of ISCO (shown by a point on each of the curves) for the same
values of $\beta$. From bottom to top, the corresponding values of the angular
momentum are $\mathcal{L}=3.46,3.64,3.74,3.81$ and $3.86$.

According to the diagram, the effective potential possesses a minimum which
allows for stable circular orbits for the particles. The latter is a necessary
condition for the formation of accretion disks in the context of the innermost
stable circular orbit (ISCO), whose corresponding radius, $r_{c}$, can be
obtained by the conditions $\mathscr{P}_{6}(r)=0=\mathscr{P}_{6}^{\prime}(r)$.
Note that, although this system of equations cannot be solved by means of
common radicals, however, the ISCO radius obeys the following relation
PhysRevD.105.023024

$r_{c}=\frac{3B(r_{c})B^{\prime}(r_{c})}{2B^{\prime}(r_{c})^{2}-B(r_{c})B^{\prime\prime}(r_{c})},$
(47)

in the spacetime geometry of the $f(R)$ black hole. This radius corresponds to
the inner edge of the accretion disk and as we move away from the black hole,
particles appear to moving on Keplerian bound orbits. In the right panel of
Fig. 5, the position of the ISCO has been indicated for each of the cases. As
it can be inferred, an increase in the $\beta$-parameter decreases $r_{c}$,
and hence, affects the structure of the accretion disk. Furthermore, the above
conditions make it possible to obtain the analytical expressions

$\displaystyle\mathcal{E}_{c}(r)=\frac{B(r)}{\sqrt{B(r)-r^{2}\varpi_{c}(r)^{2}}}=\frac{\sqrt{2\Lambda}(r_{++}-r)(r-r_{+})(r-r_{3})}{r\sqrt{\frac{3r_{++}r_{+}r_{3}}{r}+r(r_{++}+r_{+}+r_{3})-2\left[r_{++}(r_{+}+r_{3})+r_{+}r_{3}\right]}},$
(48)
$\displaystyle\mathcal{L}_{c}(r)=\frac{r^{2}\varpi_{c}(r)}{\sqrt{B(r)-r^{2}\varpi_{c}(r)^{2}}}=\frac{r^{\frac{3}{2}}\sqrt{\left(-\frac{r_{++}r_{+}r_{3}}{r^{2}}-2r+r_{++}+r_{+}+r_{3}\right)}}{\sqrt{\frac{3r_{++}r_{+}r_{3}}{r}+r(r_{++}+r_{+}+r_{3})-2\left[r_{++}(r_{+}+r_{3})+r_{+}r_{3}\right]}},$
(49)

for the energy and angular momentum of particles residing in the ISCO, in
which

$\varpi_{c}(r)=\frac{\mathrm{d}\phi}{\mathrm{d}t}=\sqrt{\frac{B^{\prime}(r)}{2r}}=\sqrt{\frac{\Lambda}{2}}\left(\frac{r_{++}+r_{+}+r_{3}}{r}-\frac{r_{++}r_{+}r_{3}}{r^{3}}-2\right)^{\frac{1}{2}},$
(50)

is the angular velocity of orbiting particles, and we have used the expression
given in Eq. (16). In Fig. 6, the radial profile of the above quantities has
been plotted for given values of the $\beta$-parameter. It is observed that by
increasing $\beta$, all the quantities increase which is a consequence of the
relevant changes in the effective potential.

(a)  (b)  (c)

Figure 6: Radial profiles of (a) $\mathcal{E}_{c}$, (b) $\mathcal{L}_{c}$ and
(c) $\varpi_{c}$, for the adopted values of the $\beta$-parameter and the same
color-coding as in Fig. 5.

Note that, on the ISCO, one can recast the characteristic polynomial as
$\mathscr{P}_{6}(r)=\Lambda r(r-r_{c})^{3}(r-r_{4})(r-r_{5})$, where $r_{4}>0$
and $r_{5}<0$ are the remaining two real roots of the characteristic
polynomial, which can be expressed in terms of $r_{c}$ by means of the method
of synthetic division.

For an accretion disk to be thin, its radius must be large compared to its
thickness. In addition, the disk is considered to be in local hydrodynamical
equilibrium at each point, which implies low pressure and vertical gradients
within the disk. We assume that the cooling process in the disk is fast enough
to prevent heat buildup due to particle friction. To ensure the stability of
the disk, we assume that the accretion rate along the radial axis,
$\mathscr{A}^{r}$, is constant all the time, in the way that

$\mathscr{A}^{r}=-2\pi\sqrt{-g}\,\Sigma\,\mathcal{U}^{r}=\mathrm{const.},$
(51)

in which $\sqrt{-g}=r^{2}$, $\Sigma$ is the surface density of the disk, and
$\mathcal{U}^{r}=\dot{r}$ is the radial component of the four-velocity of the
accreting particles. From the conservation of energy and angular momentum, one
can obtain the differential of the luminosity as page_disk-accretion_1974 ;
Joshi:2014

$\frac{\mathrm{d}\ell}{\mathrm{d}\ln r}=4\pi
r\sqrt{-g}\,\mathcal{E}_{c}(r)\mathcal{F}(r),$ (52)

where $\mathcal{F}(r)$ is the flux of the radiated energy from the disk, and
is given by

$\mathcal{F}(r)=-\frac{\mathscr{A}^{r}}{4\pi\sqrt{-g}}\frac{\varpi^{\prime}_{c}(r)}{\Big{[}\mathcal{E}_{c}(r)-\varpi_{c}(r)\mathcal{L}_{c}(r)\Big{]}^{2}}\int_{r_{c}}^{r}\Big{[}\mathcal{E}_{c}(r)-\varpi_{c}(r)\mathcal{L}_{c}(r)\Big{]}\mathcal{L}_{c}^{\prime}(r)\,\mathrm{d}r.$
(53)

By considering the fact that
$\mathcal{E}_{c}^{\prime}(r)=\varpi_{c}(r)\mathcal{L}_{c}^{\prime}(r)$, one
can write page_disk-accretion_1974

$\int_{r_{c}}^{r}\Big{[}\mathcal{E}_{c}(r)-\varpi_{c}(r)\mathcal{L}_{c}(r)\Big{]}\mathcal{L}_{c}^{\prime}(r)\,\mathrm{d}r=\mathcal{E}_{c}(r)\mathcal{L}_{c}(r)-\mathcal{E}_{c}(r_{c})\mathcal{L}_{c}(r_{c})-2\int_{r_{c}}^{r}\mathcal{L}_{c}(r)\mathcal{E}^{\prime}_{c}(r)\,\mathrm{d}r.$
(54)

Now taking the expressions in Eqs. (48), (49) and (50) up to the first-order
in $\Lambda$, and employing them in the integrand of the above relation, one
can obtain the analytical solution

$\mathcal{F}(r)=-\frac{\mathscr{A}^{r}}{4\pi\sqrt{-g}}\frac{\varpi_{c}^{\prime}(r)}{\Big{[}\mathcal{E}_{c}(r)-\varpi_{c}(r)\mathcal{L}_{c}(r)\Big{]}^{2}}\Big{[}\mathcal{E}_{c}(r)\mathcal{L}_{c}(r)-\mathcal{E}_{c}(r_{c})\mathcal{L}_{c}(r_{c})-2\mathscr{J}(r)+2\mathscr{J}(r_{c})\Big{]},$
(55)

for the flux, for which $\mathscr{J}(r)$ has been given in appendix A. Since
the disk is thin, we can assume that the emission follows the radiation of a
black body whose temperature profile is given by

$\mathcal{T}(r)^{4}=\frac{\mathcal{F}(r)}{\sigma},$ (56)

in which, $\sigma$ is the Stefan-Boltzmann constant. In Fig. 7, the above
relations have been employed to plot the radial profiles of the flux,
temperature, and differential luminosity for different values of the
$\beta$-parameter. It can be observed that by increasing $\beta$, all of the
above quantities will increase, which implies that the more the black hole
differs from the Schwarzschild-de Sitter solution, the more intense the
radiation of the accretion disk is, and the disk’s temperature becomes higher.

(a)  (b)  (c)

Figure 7: Radial profiles of (a) flux, (b) differential luminosity, and (c)
temperature, for the same values of the $\beta$-parameter and color-coding as
in Fig. 5.

###  IV.1 Shadow and rings of the black hole with thin accretion

To a distant observer, a real black hole appears as a dark shaded region
surrounded by an illuminated area. This illuminated area is formed by light
rays originating from various parts of the accretion disk, which have the
potential to escape from the black hole. As discussed in Subsec. III.1,
certain photons with specific impact parameters can escape the black hole
either through direct deflection or by following critical orbits (OFK and
COFK) as shown in diagrams (a) and (c) of Fig.4. Consequently, incoming
photons from the accretion disk may undergo different numbers of orbits around
the black hole before exiting, giving rise to several light rings that confine
the shadow.

####  IV.1.1 Direct emission, lensing rings and photon rings

Here, we follow the method introduced in Ref. Gralla:2019 , to characterize
the light rings in the observer’s sky, where the number of orbits is defined
as

$n=\frac{\phi}{2\pi},$ (57)

in which, $\phi$ represents the final azimuth angle of photons just before
they escape the black hole. The parameter $n$ corresponds to the number of
times that the light ray geodesics cross the plane of the accretion disk. In
Ref. Gralla:2019 , these rays were classified into three cases: $0.25<n<0.75$,
where the light rays intersect the accretion disk only once, forming the
direct emission; $0.75<n<1.25$, in which the rays cross the accretion disk
twice, generating the lensed (lensing) ring; and $n>1.25$, corresponding to
the formation of the photon ring, where the rays intersect the accretion disk
more than twice. Fig. 8 illustrates the behavior of $n$ with respect to the
impact parameter $b$ for different values of the $\beta$-parameter, with
distinct colors indicating the domains of direct emissions, lensing rings, and
photon rings.

(a)  (b)  (c)  (d)  (e)  (f)  (g)  (h)  (i)  (j)

Figure 8: The $b$-profile of the total number of photon orbits $n$, together with the behavior of the null geodesics in the near-horizon regions. In these diagrams, the direct, lensing ring and photon ring emissions have been color-coded appropriately. The black disk indicates the event horizon of the black hole whereas the green dashed circle denotes the radius of unstable (critical) orbits, $r_{p}$. The diagrams correspond to the cases of (a,d) $\beta=0$, (b,e) $\beta=0.011$, (c,f) $\beta=0.022$, (g,i) $\beta=0.031$, and (h,j) $\beta=0.041$. $\beta$ | Direct emission $(0.25>n>0.75)$ | Lensing ring $(0.75<n<1.25)$ |  Photon ring $(n>1.25)$
---|---|---|---
0.0 |  $b<5.01685$; $b>6.14685$ |  $5.01685<b<b_{p}$; $5.23685<b<6.14685$ |  $b_{p}<b<5.23685$
0.011 |  $b<4.80681$; $b>5.71683$ |  $4.80681<b<b_{p}$; $4.98682<b<5.71683$ |  $b_{p}<b<4.98682$
0.022 |  $b<4.62449$; $b>5.38449$ |  $4.62449<b<b_{p}$; $4.78449<b<5.38449$ |  $b_{p}<b<4.78449$
0.031 |  $b<4.47621$; $b>5.13621$ |  $4.47621<b<b_{p}$; $4.61621<b<5.13621$ |  $b_{p}<b<4.61621$
0.041 |  $b<4.33281$; $b>4.91281$ |  $4.33281<b<b_{p}$; $4.46281<b<4.91281$ |  $b_{p}<b<4.46281$
Table 2: The impact parameter domains corresponding to the direct emission,
lensing rings and photon rings of the black hole given for different values of
the $\beta$-parameter.

In this diagram, a large number of geodesics have been simulated for each
case, including the OFK, OSK, COFK, and COSK. Notably, for the remainder of
this section, we assume the value of $\Lambda=10^{-8}$, which does not
significantly alter the photon orbit properties or the characteristic
distances of the spacetime but is essential for proper functioning of our ray-
tracing codes. As seen from the diagrams, increasing $b$ leads to an increase
in the total number of orbits within the domain $b<b_{p}$ until it reaches a
narrow peak, after which it decreases within the domain $b>b_{p}$. On the
other hand, for larger values of the $\beta$-parameter, the width of the
lensing rings and photon rings shrinks. In Table 2, this has also been shown
numerically by writing down the range of $b$ for the direct, lensing ring, and
photon ring emissions, for different values of $\beta$. According to the data
presented in this table, it can be checked that by an increase in the
$\beta$-parameter, the range of $b$ for all emission types is shrunk.
Therefore, the thickness of the photon and lensing rings is decreased in this
sense. Accordingly, the angular size of the shadow is also decreased for
larger $\beta$, and hence, the contribution to the brightness of the rings is
reduced.

We continue by studying the observed emission intensity from the thin
accretion disk in the framework of the $f(R)$ black hole model.

####  IV.1.2 Transfer functions and the observed intensities

The radiation of the accretion disk is supposed to be isotropic in its rest
frame. By $I_{\mathrm{e}}(r)$, we denote the specific intensity of an emitted
radiation of frequency $\nu_{\mathrm{e}}$ from the disk. From Liouville’s
theorem, we know that the quantity
${I_{\mathrm{e}}(r)}/{\nu^{3}_{\mathrm{e}}}$ is conserved along the entire
path of light propagation. Hence, the observed intensity $I_{\mathrm{o}}$ of
frequency $\nu_{\mathrm{o}}$ are related in terms of the relation
${I_{\mathrm{e}}(r)}/{\nu^{3}_{\mathrm{e}}}={I_{\mathrm{o}}(r)}/{\nu^{3}_{\mathrm{o}}}$
bromley_line_1997 . Accordingly, we have

$I_{\mathrm{o}}(r)=\mathfrak{h}^{3}I_{\mathrm{e}}(r),$ (58)

which in our model $\mathfrak{h}=\sqrt{B(r)}$. Now by integrating over the
range of all the observed frequencies, the total observed specific intensity
is obtained as

$I_{\mathrm{o}}^{t}(r)=\int
I_{\mathrm{o}}(r)\,\mathrm{d}\nu_{\mathrm{o}}=\mathfrak{h}^{4}I_{\mathrm{emit}}(r),$
(59)

in which, the total emission intensity is given by $I_{\mathrm{emit}}(r)=\int
I_{\mathrm{e}}(r)\mathrm{d}\nu_{\mathrm{e}}$. Note that, since each
intersection of the light rays with the accretion disk generates an additional
brightness, the reliable total observed intensity of the direct emission and
the rings is given by

$I_{\mathrm{obs}}(r)=\sum_{m}I_{\mathrm{o}}^{t}(r)|_{r=r_{m}(b)},$ (60)

where $r_{m}(b)$ is the transfer function that relates the impact parameter of
the light ray trajectories with the radial coordinate of the $m$th
intersection of light rays with the accretion disk222The number of
intersections $m$ and the number of orbits $n$, are related as
$n={m}/{2}-{1}/{4}$ hu_observational_2022 .. Hence, the slope of the transfer
function indicates its (de)magnification scale Gralla:2019 ;
zeng_influence_2020 . Therefore, this slope is called the (de)magnification
factor. In Fig. 9, we have demonstrated the $b$-profile of the transfer
function for different values of $\beta$.

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

Figure 9: The transfer function $r_{m}(b)$ plotted for different values of
$\beta$. The panels (a)–(e) correspond respectively to the cases of
$\beta=0.011,0.022,0.031$ and $0.041$. The color coding is the same as that in
Fig. 8.

In these diagrams, the black line with an approximately constant slope
corresponds to the case of $m=1$ and represents direct emission. The nearly
unit slope indicates a redshifted source. For the lensing ring with $m=2$, the
impact parameter $b_{p}$ is approached, but the slope increases significantly
from one, indicating demagnification of the back side image of the accretion
disk. In the case of $m=3$ and the formation of the photon ring, the slope
tends to infinity, demonstrating an extreme demagnification of the front side
image of the accretion disk. Based on these observations, one can infer that
the contribution of the lensing ring and photon ring in the observed intensity
is negligible compared to direct emission. It is noteworthy that higher-order
rings with $n\geq 4$ (black hole subrings) do not possess significant
observational features, although they have been found to produce certain
interferometric signatures johnson_universal_2020 .

####  IV.1.3 Observational signatures of emissions from the accretion disk

In this part of the paper, we employ a ray-tracing procedure to generate the
black hole shadow along with its accretion disk image. We adopt a face-on
view, which provides greater generality and informative insights into the
silhouette imaging of black holes.

From the perspective of a distant observer, the accretion disk serves as the
primary light source that illuminates the black hole. The brightness of this
source solely depends on the radial coordinate $r$, and as discussed earlier,
it can be represented in terms of the emitted intensity $I_{\mathrm{emit}}$.
To explore the observational features of the $f(R)$ black hole further, we
consider three toy models for the intensity profile of the thin accretion
disk, as described below:

  * •

Model 1: In this model, the emission comes from the ISCO, and the intensity
profile is given by the decaying function

$\displaystyle
I_{\mathrm{emit}}(r)=\begin{cases}{\left[r-(r_{c}-1)\right]^{-2}}&\text{for
$r>r_{c}$}\\\ 0&\text{for $r\leq r_{c}$}\end{cases}.$ (61)

  * •

Model 2: We assume that the radiation is originated from the photon sphere of
the radius $r_{p}$, and the emission intensity profile is expressed as

$\displaystyle
I_{\mathrm{emit}}(r)=\begin{cases}{\left[r-(r_{p}-1)\right]^{-3}}&\text{for
$r>r_{p}$}\\\ 0&\text{for $r\leq r_{p}$}\end{cases}.$ (62)

  * •

Model 3: For the case that the emission starts from the event horizon radius
$r_{+}$, we consider an emission profile in which, the decay is more moderate
compared with the last two models, and is given by

$\displaystyle
I_{\mathrm{emit}}(r)=\begin{cases}\left[\frac{\pi}{2}-\arctan\left(r-[r_{c}-1]\right)\right]\left[{\frac{\pi}{2}-\arctan(r_{p})}\right]^{-1}&\text{for
$r>r_{+}$}\\\ 0&\text{for $r\leq r_{+}$}\end{cases}.$ (63)

Table 3 presents the radii of the event horizon, the ISCO, and the photon
sphere, serving as the origins for the aforementioned emission models, based
on different values of the $\beta$-parameter considered thus far.

$\beta$ | $r_{+}$ | $r_{c}$ | $r_{p}$
---|---|---|---
0.0 | $2.00008$ | $6.01981$ | $3.00000$
0.011 | $1.95638$ | $5.31895$ | $2.95031$
0.022 | $1.91959$ | $5.01730$ | $2.90767$
0.031 | $1.88795$ | $4.82752$ | $2.87042$
0.041 | $1.85826$ | $4.68106$ | $2.83497$
Table 3: The radii of the origins for the radiation emission models, provided
for various values of the $\beta$-parameter.

Note that, the above models have their own specific properties respecting the
black hole shadow, and the second model emission profile shows the largest
decay. These models, despite being rather idealized, however, can provide
useful insights into the light propagation in the exterior of black holes. In
Fig. 10–14, the observational appearance of the accretion disk around the
$f(R)$ black hole has been shown for each of the above models, together with
the plots of the emitted and observed intensities for each of the cases.

(a)  (b)  (c)

Figure 10: The observational characteristics of the accretion disk around the
$f(R)$ black hole, depicted for the case of $\beta=0$ (the Schwarzschild-de
Sitter black hole). The panels, from top to bottom, correspond to (a) model 1,
(b) model 2, and (c) model 3 emission profiles. The left and middle panels in
each row display the $b$-profiles of the emitted and observed intensities,
respectively. The right panels present two-dimensional face-on ray-traced
shadow images for each of the models.

(a)  (b)  (c)

Figure 11: The case of $\beta=0.011$.

(a)  (b)  (c)

Figure 12: The case of $\beta=0.022$.

(a)  (b)  (c)

Figure 13: The case of $\beta=0.031$.

(a)  (b)  (c)

Figure 14: The case of $\beta=0.041$.

Figure 15: Blurring the shadow in Fig. 12(a) for $\beta=0.022$ using a
Gaussian filter, to emulate the EHT nominal resolutions for the images of M87*
and Sgr A*. In the left panel, the starting radius of the direct emission has
been shown to be about 5.69, which forms the boundary of the shadow. After
applying the Gaussian filter, the lensing and photon rings disappear in the
right panel, and hence, the radius of the black hole shadow is estimated as
5.69. This value is much larger than the radius of the photon ring, which is
$b_{p}=4.74$ for the case of $\beta=0.022$.

The figures in each row represent, respectively, the emitted and observed
intensity, and the shadow of the $f(R)$ black hole for models 1, 2, and 3. In
model 1, the emitted intensity shows an asymptotic behavior near $b_{p}$,
decreasing as the radial distance increases and approaching zero. Here,
spherical photon orbits occur inside the disk’s emission region. The observed
intensity exhibits two independent peaks within the lensing ring and photon
ring domains. For all values of the $\beta$-parameter, except $\beta=0.031$,
the photon ring intensity is smaller than the direct emission, while the
lensing ring intensity is always larger333In the two-dimensional ray-traced
shadow images, zooming in reveals the photon ring as a thin circle inside the
thicker lensing ring.. Both peaks have remarkably narrow observational ranges,
indicating that, at long distances where the observer is located, the
contribution of the lensing and photon rings to the observed intensity is
dominated by direct emission. Thus, the observed emission from the black hole
in this model is mainly direct emission, as evident from the shadow images.
Furthermore, comparing the diagrams for different $\beta$ values reveals that
the shadow size decreases with an increase in this parameter, making the
Schwarzschild-de Sitter black hole have the largest shadow for this model. In
model 2, the emitted intensity peaks at $r_{p}$ and then sharply drops with
increasing radial distance. The observed intensity has two peaks, with the
first corresponding to direct emission and the second to a combination of
lensing and photon rings. In some cases of $\beta$, both peaks have
significant ranges and approximately equal intensities. For most cases of
$\beta$, the direct emission dominates the observed intensity, while the rings
are strongly demagnified. However, in the exceptional case of $\beta=0.031$,
both direct emission and the rings contribute relatively equally to the
observed intensity, as evident from the corresponding shadow image. In model
3, the emitted intensity peaks at the event horizon $r_{+}$ and decreases with
radial distance. In this case, direct emission, lensing ring, and photon ring
merge and occupy a significant range in the observational domain. Outside the
event horizon, there is a smooth uplift in the observed intensity profiles,
where direct emission dominates, followed by an intense peak corresponding to
the lensing and photon rings. In the case of $\beta=0$ (the Schwarzschild-de
Sitter black hole), the observed intensity exhibits a narrow but intense peak
for the photon ring, followed by a smaller peak where both rings contribute,
forming a wide and bright ring. This bright ring is observed in the intensity
profiles and shadow images for all cases of the $\beta$-parameter.
Additionally, the brightness of the accretion disk increases with an increase
in the $\beta$-parameter. Notably, a thin accretion disk significantly
influences the size of the observed black hole shadow. For instance, in Fig.
15, we have applied a Gaussian filter to simulate the angular resolution
generated by the EHT on the shadow image presented in Fig. 12(a) as a
reference. Calculating the maximums of the observed emission in Eq. (60) for
this case yields an estimated radius of $5.69$ for the direct emission peak,
which is significantly larger than the theoretical value ($b_{p}=4.7445$ for
$\beta=0.022$). This difference arises from changes in both the
$\beta$-parameter and the disk’s emission profile. Therefore, directly
inferring the value of $b_{p}$ from the black hole shadow size makes it
challenging to test general relativity using the EHT results.

####  IV.1.4 Observational signatures of infalling spherical accretion

Here, we investigate the shadow cast of the $f(R)$ black hole while it
accretes spherically the radiative gas, constituting its thin emission disk
bambi_can_2013 . In this model, the observed intensity is expressed as

$I_{\mathrm{obs}}=\int_{\bm{\gamma}}\mathscr{R}^{3}\mathcal{J}(\nu_{\mathrm{e}})\,\mathrm{d}I_{\mathrm{prop}},$
(64)

over the null geodesic congruence $\bm{\gamma}$, in which $\mathscr{R}$ is the
redshift factor, $\nu_{\mathrm{e}}$ is the frequency of emitted photons from
the accretion disk, $\mathrm{d}I_{\mathrm{prop}}$ is the infinitesimal proper
length, and

$\mathcal{J}(\nu_{\mathrm{e}})\propto\frac{\delta(\nu_{\mathrm{e}}-\nu_{\mathrm{f}})}{r^{2}},$
(65)

is the permittivity per unit volume in the emitter’s rest frame, in which
$\nu_{\mathrm{f}}$ is the monochromatic rest-frame emission frequency, and
$\delta$ is the delta function. In this construction, the redshift factor is
given by

$\mathscr{R}=\frac{\Pi_{\mu}u_{\mathrm{o}}^{\mu}}{\Pi_{\nu}u_{\mathrm{e}}^{\nu}},$
(66)

where $\bm{u}_{\mathrm{o}}$ and $\bm{u}_{\mathrm{e}}$ are, respectively, the
four-velocities associated with a distant static observer, and the infalling
accreting matter. Accordingly, $u_{\mathrm{o}}^{\mu}=(1,0,0,0)$, and in the
spacetime of the $f(R)$ black hole we can write

$u_{\mathrm{e}}^{\mu}=\left(\frac{1}{B(r)},-\sqrt{1-B(r)},0,0\right).$ (67)

The $\bm{\Pi}$ covector in Eq. (66) is the four-momentum of the emitted
photons from the accretion disk and has the same definition as in Eq. (18).
Since the accretion is supposed to be only in the radial direction, it is then
sufficient to recalculate the fraction of the temporal and radial components
of $\bm{\Pi}$, which yields bambi_can_2013

$\frac{\Pi_{r}}{\Pi_{t}}=\pm\frac{1}{B(r)}\sqrt{1-B(r)\frac{b^{2}}{r^{2}}},$
(68)

in which the $\pm$ signs correspond, respectively, to whether the photons
approach or recede from the black hole. One can therefore recast the redshift
factor as

$\displaystyle\mathscr{R}$ $\displaystyle=$
$\displaystyle\left(u_{\mathrm{e}}^{t}+\frac{\Pi_{r}}{\Pi_{t}}u_{\mathrm{e}}^{r}\right)^{-1}$
(69) $\displaystyle=$
$\displaystyle\left[\frac{1}{B(r)}\pm\sqrt{\left(\frac{1}{B(r)}-1\right)\left(\frac{1}{B(r)}-\frac{b^{2}}{r^{2}}\right)}\,\right]^{-1},$

and from here, the infinitesimal proper length is obtained as

$\mathrm{d}I_{\mathrm{prop}}=\Pi_{\mu}u_{\mathrm{e}}^{\mu}\mathrm{d}\tau=\frac{\Pi_{t}}{\mathscr{R}|\Pi_{r}|}\mathrm{d}r.$
(70)

Accordingly, Eq. (64) takes the form

$I_{\mathrm{obs}}\propto\int_{\bm{\gamma}}\frac{\mathscr{R}^{3}}{r^{2}}\frac{\Pi_{t}}{|\Pi_{r}|}\mathrm{d}r,$
(71)

and the observed intensity is, therefore, obtained by doing the above
integration over all the frequencies. Using this method, we can study the
brightness of infalling accretion and the shadow of the $f(R)$ black hole. As
shown in Fig. 16, for all cases of the $\beta$-parameter, by the increase in
the impact parameter and by moving away from the origin, the specific observed
intensity increases until it reaches a peak around $b_{p}$, and falls of
remarkably in the region $b>b_{p}$ and goes to zero. On the other hand, by
altering the $\beta$-parameter from zero, the peak becomes slightly higher and
the bottom line of the profile is lifted by the same value, whereas its width
is decreased relevantly. This means that for larger $\beta$, the accretion
disk appears brighter to the observer, and the silhouette becomes less dark
but smaller in size. Hence, to a distant observer, the Schwarzschild-de Sitter
black hole has the darkest and the largest silhouette, and the least bright
accretion disk. In Fig. 17, these effects have been visualized for the $f(R)$
black hole.

(a)  (b)

Figure 16: (a) The $b$-profiles of the observed intensity of the infalling
spherical accretion, which from bottom to top correspond to the cases of
$\beta=0,0.011,0.022,0.031$ and $0.041$. (b) The same as panel (a), but
showing only a part of the $b$ domain within the positive values.

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

Figure 17: The images of the disk and the silhouette of the black hole with
infalling spherical accretion, given for (a) $\beta=0$, (b) $\beta=0.011$, (c)
$\beta=0.022$, (d) $\beta=0.031$ and (e) $\beta=0.041$.

###  IV.2 Discussion

As mentioned in Section II, the first and second-order linear terms present in
the lapse function (7) can be interpreted as mimicking the flat galactic
rotation curves and the accelerated expansion of the universe, and hence, they
can be associated with the dark matter and dark energy components of the
$f(R)$ black hole spacetime. In recent years, similar investigations to the
one presented in this paper have been conducted on spherically symmetric black
holes with linear or nonlinear extra terms, aiming to account for the
properties of the dark side of the universe. Notably, the lapse function (7)
exhibits similarities with quintessential dark fields, which have also been
explored in Refs. zeng_influence_2020 ; he_shadow_2022  concerning the shadow
and observational signatures of infalling spherical accretions. When
quintessential dark energy is employed in a black hole spacetime model, the
lapse function takes on the form

$\bar{\mathrm{B}}(r)=1-\frac{2M}{r}-\frac{q}{r^{3w+1}},$ (72)

where the quintessential and state parameters, are denoted by $q$ and $w$,
respectively. For $w=-1$, the resulting lapse function resembles the
Schwarzschild-de Sitter black hole spacetime, which has not been addressed in
the previously mentioned references. However, in this paper, we consider this
case as a baseline to facilitate reliable comparisons between the impacts of
different values for the $f(R)$ black hole parameters (see Figs. 10 and 17).
When $w=-1/3$, the lapse function becomes $\bar{\mathrm{B}}(r)=1-q-2M/r$,
corresponding to a Schwazschild black hole associated with cloud of strings
Letelier:1979 . Therefore, the quintessential dark energy field is obtained
within the range $-1<w<-1/3$. Another important case is achieved by setting
$w=-2/3$, resulting in $\bar{\mathrm{B}}(r)=1-qr-2M/r$, comparable to the
lapse function (7) for $\Lambda=0$. This specific case has been studied in
Ref. zeng_influence_2020 , where only the quintessential dark energy affects
the spacetime, and in Ref. he_shadow_2022 , when the cloud of strings is also
present. However, the resultant spacetimes for these cases cannot compensate
for the flat rotation curves in galactic scales, as this requires a positive
contribution of the linear term (i.e., for $q<0$ in $\bar{\mathrm{B}}(r)$),
which has not been addressed in the aforementioned references. On the other
hand, this aspect is fully considered in our paper, as all studied examples
for the $f(R)$ black hole involve a positive $\beta$-parameter. Similar
considerations have been explored in Ref. Li:2021ypw , where a global monopole
black hole in $f(R)$ gravity exhibits a lapse function similar to the
Schwarzschild black hole associated with quintessence and cloud of strings, as
discussed in Refs. he_shadow_2022 ; Fathi:2022a . However, the aforementioned
$f(R)$ black hole contains a negative linear term, and thus cannot compensate
for dark matter effects in galactic scales. In contrast, our study benefits
from the positive first-order term, allowing for compensation at both large
and galactic scales. Furthermore, our investigation also accounts for
gravitational effects at large distances, based on the existence of the square
term $\Lambda r^{2}$ in the lapse function. Therefore, this investigation
offers some important new features that are potentially interesting regarding
current astrophysical observations.

Having studied the analytical structure of light propagation around the $f(R)$
black hole and demonstrated its observational features across various
criteria, we can now summarize our results in the next section.

##  V Summary and Conclusion

The theoretical study of black holes with accretion disks provides a more
realistic framework for making reliable comparisons with observational
effects. In this work, we focused on a static spherically symmetric black hole
derived from a special $f(R)$ theory of gravity, which is compatible with both
small and large-scale structure tests. This black hole features a linear
first-order term with a coefficient $\beta$, along with a cosmological
constant. While the cosmological constant is meant to compensate for vacuum
energy and the accelerated expansion of the universe, the $\beta$-parameter is
responsible for mimicking flat galactic rotation curves. We initially derived
the theoretical diameter of the black hole’s photon ring and demonstrated its
behavior concerning variations in the parameter $\beta$. However, we clarified
that this parameter cannot be constrained by comparing the above theoretical
diameter with the observed shadow diameters of M87* and Sgr A* since the
photon ring is not observed in the images generated by the EHT. Consequently,
for the numerical computations throughout the paper, we adopted a specific
range for the parameter, setting it between $0$ and $0.05$. This range allows
us to incorporate the Schwarzschild-de Sitter black hole scenario and funds a
more realistic investigation. Next, we delved into solving the angular
equations of motion for deflecting and critical trajectories, providing exact
analytical solutions for each case. These types of orbits play a crucial role
in shaping the shadow of black holes when illuminated by an accretion disk. By
employing the obtained analytical solutions for deflecting trajectories, we
derived the lens equation and calculated the deflection angle, which turned
out to be around 10 $\mathrm{\mu as}$ for the assumed values of the
$\beta$-parameter. In the latter part of our investigation, we constructed a
geometrically and optically thin accretion disk around the black hole using
the Novikov-Thorne model and computed the disk’s characteristics. We then
applied the method introduced in Ref. Gralla:2019  to classify the light rings
and different types of accretion emission profiles. Based on the number of
half orbits $n$ completed by light rays around the black hole, the rings were
categorized into lensing rings and photon rings, both of which experienced
demagnification due to extreme gravitational lensing. We also determined the
impact parameter ranges for each ring, providing their respective thicknesses.
By considering three types of disk emission profiles, we observed that an
increase in the $\beta$-parameter leads to a decrease in the size of the
shadow compared to that of the Schwarzschild-de Sitter black hole.
Additionally, the brightness of the rings may vary depending on the value of
$\beta$ and the radial position of the direct emission. Furthermore, the
accretion disk’s brightness is enhanced with an increase in the
$\beta$-parameter. Moreover, we applied a Gaussian filter to a particular case
for simulating EHT images, further demonstrating that inferring the black hole
size directly from observations remains inconclusive, as we have detailed
earlier. Lastly, we explored a spherically asymmetric infalling accretion
scenario for the black hole and obtained the observed intensity profiles for
different $\beta$ values. Our analysis demonstrated that as the
$\beta$-parameter increases, the black hole becomes brighter in terms of the
accretion disk, while the silhouette gradually shrinks. These findings were
visually confirmed through appropriate ray-tracing methods. An interesting
topic for future exploration could involve examining a plasmic medium
characterized by a specific index profile encircling the black hole. This
approach allows us to delve into how plasma components might impact the
observable characteristics of black holes, serving as more reliable
constituents within the accretion disks.

## Acknowledgements

We thank the referee for valuable comments, which greatly helped us improve
the paper’s discussion. We would like to thank Mert Okyay and Ali Övgün for
providing some Mathematica codes they have used in Ref. okyay_nonlinear_2022 .
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.

##  Appendix A The full expression of $\mathscr{J}(r)$

Direct integration of the integral in Eq. (54) results in

$\displaystyle\int\mathcal{L}_{c}(r)\mathcal{E}^{\prime}_{c}(r)\,\mathrm{d}r\equiv\mathscr{J}(r)$
$\displaystyle=$
$\displaystyle\frac{1}{30\sqrt{2}\beta^{3}\sqrt{\frac{r}{\beta
r^{2}+2}}}\left\\{\frac{2\beta r\left(-30\Lambda-24\beta^{2}\Lambda
r^{3}+2\beta
r^{2}\left(5\beta^{2}+\Lambda\right)+r\left(5\beta^{2}+9\beta\Lambda+10\Lambda\right)\right)}{\beta
r^{2}+2r-6}\right.$
$\displaystyle+\left.\frac{{\beta^{\frac{3}{4}}}}{\sqrt{(6\beta+1)}\left(\beta
r^{2}+2\right)}\left[-60\sqrt{6\beta+1}\beta^{\frac{7}{4}}-30\sqrt{6\beta+1}\beta^{\frac{3}{4}}r^{2}-120\sqrt{\frac{(6\beta+1)}{\sqrt{\beta}}}\Lambda\right.\right.$
$\displaystyle\left.\left.-324\sqrt{\frac{(6\beta+1)}{\sqrt{\beta}}}\beta\Lambda-162\sqrt{6\beta+1}\beta^{\frac{7}{4}}\Lambda
r^{2}-60\beta\sqrt{\frac{(6\beta+1)}{\sqrt{\beta}}}\Lambda
r^{2}\right.\right.$ $\displaystyle\left.\left.+60\
2^{\frac{3}{4}}\Lambda\sqrt{(6\beta+1)r\left(\beta
r^{2}+2\right)}\,\mathbf{E}\left(\left.\operatorname{arcsinh}\left(\frac{2^{\frac{3}{4}}}{\sqrt{r}{\beta^{\frac{3}{4}}}}\right)\right|-1\right)\right.\right.$
$\displaystyle\left.\left.+6\
{2^{\frac{3}{4}}}\beta(5\beta+27\Lambda)\sqrt{(6\beta+1)r\left(\beta
r^{2}+2\right)}\,\mathbf{E}\left(\left.\operatorname{arcsinh}\left(\frac{1}{{2^{\frac{3}{4}}}\sqrt{r}{\beta^{\frac{3}{4}}}}\right)\right|-1\right)\right.\right.$
$\displaystyle\left.\left.-30\ 2^{\frac{3}{4}}\sqrt{r(6\beta+1)\left(\beta
r^{2}+2\right)}\,\mathbf{F}\left(\left.\operatorname{arcsinh}\left(\frac{1}{\sqrt[4]{2\beta}\sqrt{r}}\right)\right|-1\right)\beta^{2}\right.\right.$
$\displaystyle\left.\left.-162\ 2^{\frac{3}{4}}\sqrt{r(6\beta+1)\left(\beta
r^{2}+2\right)}\Lambda\,\mathbf{F}\left(\left.\operatorname{arcsinh}\left(\frac{1}{\sqrt[4]{2\beta}\sqrt{r}}\right)\right|-1\right)\beta\right.\right.$
$\displaystyle\left.\left.-10\
2^{\frac{3}{4}}\sqrt{r\beta(6\beta+1)\left(\beta
r^{2}+2\right)}\Lambda\,\mathbf{F}\left(\left.\operatorname{arcsinh}\left(\frac{1}{\sqrt[4]{2\beta}\sqrt{r}}\right)\right|-1\right)\right.\right.$
$\displaystyle\left.\left.-20\
2^{\frac{3}{4}}\sqrt{r\beta^{5}(6\beta+1)\left(\beta
r^{2}+2\right)}\,\mathbf{F}\left(\left.\operatorname{arcsinh}\left(\frac{1}{\sqrt[4]{2\beta}\sqrt{r}}\right)\right|-1\right)\right.\right.$
$\displaystyle\left.\left.-30\
2^{\frac{3}{4}}\sqrt{r^{3}+\frac{2r}{\beta}}\Lambda\,\mathbf{\Pi}\left(\frac{3\sqrt{2\beta}}{\sqrt{6\beta+1}+1};\left.\operatorname{arcsinh}\left(\frac{1}{\sqrt[4]{2\beta}\sqrt{r}}\right)\right|-1\right)\right.\right.$
$\displaystyle\left.\left.-225\ 2^{\frac{3}{4}}\sqrt{r\beta\left(\beta
r^{2}+2\right)}\Lambda\,\mathbf{\Pi}\left(\frac{3\sqrt{2\beta}}{\sqrt{6\beta+1}+1};\left.\operatorname{arcsinh}\left(\frac{1}{\sqrt[4]{2\beta}\sqrt{r}}\right)\right|-1\right)\right.\right.$
$\displaystyle\left.\left.-270\ 2^{\frac{3}{4}}\sqrt{r\beta^{3}\left(\beta
r^{2}+2\right)}\Lambda\,\mathbf{\Pi}\left(\frac{3\sqrt{2\beta}}{\sqrt{6\beta+1}+1};\left.\operatorname{arcsinh}\left(\frac{1}{\sqrt[4]{2\beta}\sqrt{r}}\right)\right|-1\right)\right.\right.$
$\displaystyle\left.\left.-15\ 2^{\frac{3}{4}}\sqrt{r\beta^{3}\left(\beta
r^{2}+2\right)}\,\mathbf{\Pi}\left(\frac{3\sqrt{2\beta}}{\sqrt{6\beta+1}+1};\left.\operatorname{arcsinh}\left(\frac{1}{\sqrt[4]{2\beta}\sqrt{r}}\right)\right|-1\right)\right.\right.$
$\displaystyle\left.\left.-90\ 2^{\frac{3}{4}}\sqrt{r\beta^{5}\left(\beta
r^{2}+2\right)}\,\mathbf{\Pi}\left(\frac{3\sqrt{2\beta}}{\sqrt{6\beta+1}+1};\left.\operatorname{arcsinh}\left(\frac{1}{\sqrt[4]{2\beta}\sqrt{r}}\right)\right|-1\right)\right.\right.$
$\displaystyle\left.\left.+30\
2^{\frac{3}{4}}\Lambda\,\mathbf{\Pi}\left(-\frac{3\sqrt{2\beta}}{\sqrt{6\beta+1}-1};\left.\operatorname{arcsinh}\left(\frac{1}{\sqrt[4]{2\beta}\sqrt{r}}\right)\right|-1\right)\sqrt{r^{3}+\frac{2r}{\beta}}\right.\right.$
$\displaystyle\left.\left.+30\
2^{\frac{3}{4}}\Lambda\,\mathbf{\Pi}\left(-\frac{3\sqrt{2\beta}}{\sqrt{6\beta+1}-1};\left.\operatorname{arcsinh}\left(\frac{1}{\sqrt[4]{2\beta}\sqrt{r}}\right)\right|-1\right)\sqrt{\frac{r(6\beta+1)\left(\beta
r^{2}+2\right)}{\beta}}\right.\right.$ $\displaystyle\left.\left.+30\
2^{\frac{3}{4}}\Lambda\,\mathbf{\Pi}\left(\frac{3\sqrt{2\beta}}{\sqrt{6\beta+1}+1};\left.\operatorname{arcsinh}\left(\frac{1}{\sqrt[4]{2\beta}\sqrt{r}}\right)\right|-1\right)\sqrt{\frac{r(6\beta+1)\left(\beta
r^{2}+2\right)}{\beta}}\right.\right.$ $\displaystyle\left.\left.+225\
2^{\frac{3}{4}}\Lambda\,\mathbf{\Pi}\left(-\frac{3\sqrt{2\beta}}{\sqrt{6\beta+1}-1};\left.\operatorname{arcsinh}\left(\frac{1}{\sqrt[4]{2\beta}\sqrt{r}}\right)\right|-1\right)\sqrt{r\beta\left(\beta
r^{2}+2\right)}\right.\right.$ $\displaystyle\left.\left.+15\
2^{\frac{3}{4}}\,\mathbf{\Pi}\left(-\frac{3\sqrt{2\beta}}{\sqrt{6\beta+1}-1};\left.\operatorname{arcsinh}\left(\frac{1}{\sqrt[4]{2\beta}\sqrt{r}}\right)\right|-1\right)\sqrt{r\beta^{3}\left(\beta
r^{2}+2\right)}\right.\right.$ $\displaystyle\left.\left.+270\
2^{\frac{3}{4}}\Lambda\,\mathbf{\Pi}\left(-\frac{3\sqrt{2\beta}}{\sqrt{6\beta+1}-1};\left.\operatorname{arcsinh}\left(\frac{1}{\sqrt[4]{2\beta}\sqrt{r}}\right)\right|-1\right)\sqrt{r\beta^{3}\left(\beta
r^{2}+2\right)}\right.\right.$ $\displaystyle\left.\left.+90\
2^{\frac{3}{4}}\,\mathbf{\Pi}\left(-\frac{3\sqrt{2\beta}}{\sqrt{6\beta+1}-1};\left.\operatorname{arcsinh}\left(\frac{1}{\sqrt[4]{2\beta}\sqrt{r}}\right)\right|-1\right)\sqrt{r\beta^{5}\left(\beta
r^{2}+2\right)}\right.\right.$ $\displaystyle\left.\left.+60\
2^{\frac{1}{4}}\Lambda\,\mathbf{F}\left(\left.\operatorname{arcsinh}\left(\frac{{2^{\frac{1}{4}}}}{\sqrt{r}\beta^{\frac{1}{4}}}\right)\right|-1\right)\sqrt{r(6\beta+1)\left(\beta
r^{2}+2\right)}\right.\right.$ $\displaystyle\left.\left.+135\
2^{\frac{3}{4}}\Lambda\,\mathbf{\Pi}\left(-\frac{3\sqrt{2\beta}}{\sqrt{6\beta+1}-1};\left.\operatorname{arcsinh}\left(\frac{1}{\sqrt[4]{2\beta}\sqrt{r}}\right)\right|-1\right)\sqrt{r\beta(6\beta+1)\left(\beta
r^{2}+2\right)}\right.\right.$ $\displaystyle\left.\left.+135\
2^{\frac{3}{4}}\Lambda\,\mathbf{\Pi}\left(\frac{3\sqrt{2\beta}}{\sqrt{6\beta+1}+1};\left.\operatorname{arcsinh}\left(\frac{1}{\sqrt[4]{2\beta}\sqrt{r}}\right)\right|-1\right)\sqrt{r\beta(6\beta+1)\left(\beta
r^{2}+2\right)}\right.\right.$ $\displaystyle\left.\left.+15\
2^{\frac{3}{4}}\,\mathbf{\Pi}\left(-\frac{3\sqrt{2\beta}}{\sqrt{6\beta+1}-1};\left.\operatorname{arcsinh}\left(\frac{1}{\sqrt[4]{2\beta}\sqrt{r}}\right)\right|-1\right)\sqrt{r\beta^{3}(6\beta+1)\left(\beta
r^{2}+2\right)}\right.\right.$ $\displaystyle\left.\left.+15\
2^{\frac{3}{4}}\,\mathbf{\Pi}\left(\frac{3\sqrt{2\beta}}{\sqrt{6\beta+1}+1};\left.\operatorname{arcsinh}\left(\frac{1}{\sqrt[4]{2\beta}\sqrt{r}}\right)\right|-1\right)\sqrt{r\beta^{3}(6\beta+1)\left(\beta
r^{2}+2\right)}\right.\right.$ $\displaystyle\left.\left.+60\
2^{\frac{3}{4}}\,\mathbf{\Pi}\left(-\frac{3\sqrt{2\beta}}{\sqrt{6\beta+1}-1};\left.\operatorname{arcsinh}\left(\frac{1}{\sqrt[4]{2\beta}\sqrt{r}}\right)\right|-1\right)\sqrt{r\beta^{5}(6\beta+1)\left(\beta
r^{2}+2\right)}\right.\right.$ $\displaystyle\left.\left.+60\
2^{\frac{3}{4}}\,\mathbf{\Pi}\left(\frac{3\sqrt{2\beta}}{\sqrt{6\beta+1}+1};\left.\operatorname{arcsinh}\left(\frac{1}{\sqrt[4]{2\beta}\sqrt{r}}\right)\right|-1\right)\sqrt{r\beta^{5}(6\beta+1)\left(\beta
r^{2}+2\right)}\right]\right\\},$ (73)

where $\mathbf{F}(\varphi|\mathfrak{m})$, $\mathbf{E}(\varphi|\mathfrak{m})$
and $\mathbf{\Pi}(\mathfrak{n};\varphi|\mathfrak{m})$ are, respectively, the
incomplete elliptic integrals of the first, second, and third kind of argument
$\varphi$, modulus $\mathfrak{m}$ and characteristic $\mathfrak{n}$
byrd_handbook_1971 . Note that the above expression does not cover the case of
$\beta=0$, so the corresponding profile has to be obtained by doing numerical
integration of Eq. (54).

## References

  * (1) K. Schwarzschild, “Über das gravitationsfeld eines massenpunktes nach der einsteinschen theorie,” Sitzungsberichte der Königlich Preußischen Akademie der Wissenschaften (Berlin), pp. 189–196, Jan. 1916.
  * (2) D. Finkelstein, “Past-future asymmetry of the gravitational field of a point particle,” Phys. Rev., vol. 110, pp. 965–967, May 1958.
  * (3) B. L. Webster and P. Murdin, “Cygnus X-1—a Spectroscopic Binary with a Heavy Companion ?,” Nature, vol. 235, pp. 37–38, Jan. 1972.
  * (4) C. T. Bolton, “Identification of Cygnus X-1 with HDE 226868,” Nature, vol. 235, pp. 271–273, Feb. 1972.
  * (5) 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.
  * (6) 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.
  * (7) K. Akiyama et al., “First M87 Event Horizon Telescope Results. VII. Polarization of the Ring,” Astrophys. J. Lett., vol. 910, p. L12, Mar. 2021.
  * (8) K. Akiyama et al., “First M87 Event Horizon Telescope Results. VIII. Magnetic Field Structure near The Event Horizon,” Astrophys. J. Lett., vol. 910, p. L13, Mar. 2021.
  * (9) P. Kocherlakota et al. Phys. Rev. D, vol. 103, p. 104047, May 2021.
  * (10) P. V. P. Cunha, C. A. R. Herdeiro, and M. J. Rodriguez, “Does the black hole shadow probe the event horizon geometry?,” Phys. Rev. D, vol. 97, p. 084020, Apr 2018.
  * (11) R. Takahashi, “Shapes and positions of black hole shadows in accretion disks and spin parameters of black holes,” The Astrophysical Journal, vol. 611, p. 996, aug 2004.
  * (12) V. C. Rubin, W. K. Ford, Jr., and N. Thonnard, “Rotational properties of 21 sc galaxies with a large range of luminosities and radii, from ngc 4605 /r = 4kpc/ to ugc 2885 /r = 122 kpc/,” Astrophys. J, vol. 238, pp. 471–487, jun 1980.
  * (13) R. Massey, T. Kitching, and J. Richard, “The dark matter of gravitational lensing,” Rept. Prog. Phys., vol. 73, p. 086901, 2010.
  * (14) K. Bolejko, C. Clarkson, R. Maartens, D. Bacon, N. Meures, and E. Beynon, “Antilensing: The bright side of voids,” Phys. Rev. Lett., vol. 110, p. 021302, Jan 2013.
  * (15) A. G. Riess et al., “Observational evidence from supernovae for an accelerating universe and a cosmological constant,” Astron. J., vol. 116, pp. 1009–1038, 1998.
  * (16) S. Perlmutter et al., “Measurements of omega and lambda from 42 high redshift supernovae,” Astrophys. J., vol. 517, pp. 565–586, 1999.
  * (17) P. Astier, “The expansion of the universe observed with supernovae,” Rept. Prog. Phys., vol. 75, p. 116901, 2012.
  * (18) T. P. Sotiriou and V. Faraoni, “f(R) theories of gravity,” Rev. Mod. Phys., vol. 82, pp. 451–497, Mar 2010.
  * (19) A. De Felice and S. Tsujikawa, “f(R) Theories,” Living Reviews in Relativity, vol. 13, p. 3, Dec. 2010.
  * (20) S. Capozziello and M. De Laurentis, “Extended theories of gravity,” Physics Reports, vol. 509, no. 4, pp. 167–321, 2011.
  * (21) I. Navarro and K. V. Acoleyen, “f(R) actions, cosmic acceleration and local tests of gravity,” Journal of Cosmology and Astroparticle Physics, vol. 2007, p. 022, feb 2007.
  * (22) S. Nojiri and S. D. Odintsov, “Unifying inflation with $\Lambda$CDM epoch in modified f(R) gravity consistent with solar system tests,” Physics Letters B, vol. 657, no. 4, pp. 238–245, 2007.
  * (23) J.-A. Gu, “Cosmological and solar-system tests of f(R) modified gravity,” International Journal of Modern Physics D, vol. 20, no. 08, pp. 1357–1362, 2011.
  * (24) S. Capozziello, V. F. Cardone, and V. Salzano, “Cosmography of f(R) gravity,” Phys. Rev. D, vol. 78, p. 063504, Sep 2008.
  * (25) S. Capozziello and A. D. Felice, “f(R) cosmology from noether’s symmetry,” Journal of Cosmology and Astroparticle Physics, vol. 2008, p. 016, aug 2008.
  * (26) G. Cognola, E. Elizalde, S. Nojiri, S. D. Odintsov, L. Sebastiani, and S. Zerbini, “Class of viable modified $f({R})$ gravities describing inflation and the onset of accelerated expansion,” Physical Review D, vol. 77, p. 046009, Feb. 2008.
  * (27) S. Nojiri and S. D. Odintsov, “Modified $f(r)$ gravity unifying $r^{m}$ inflation with the $\lambda$cdm epoch,” Physical Review D, vol. 77, p. 026007, Jan. 2008.
  * (28) S. Nojiri and S. D. Odintsov, “Unified cosmic history in modified gravity: From ${F}({R})$ theory to Lorentz non-invariant models,” Physics Reports, vol. 505, pp. 59–144, Aug. 2011.
  * (29) S. Nojiri, S. Odintsov, and V. Oikonomou, “Modified gravity theories on a nutshell: Inflation, bounce and late-time evolution,” Physics Reports, vol. 692, pp. 1–104, June 2017.
  * (30) V. Kagramanova, J. Kunz, and C. Lämmerzahl, “Solar system effects in schwarzschild–de sitter space–time,” Physics Letters B, vol. 634, no. 5, pp. 465–470, 2006.
  * (31) A. A. Starobinsky, “Disappearing cosmological constant in f(R) gravity,” JETP Letters, vol. 86, pp. 157–163, Oct. 2007.
  * (32) W. Hu and I. Sawicki, “Models of f(R) cosmic acceleration that evade solar system tests,” Phys. Rev. D, vol. 76, p. 064004, Sep 2007.
  * (33) Sobouti, Y., “An f(R) gravitation for galactic environments,” A&A, vol. 464, no. 3, pp. 921–925, 2007.
  * (34) Saffari, R. and Sobouti, Y., “An f(R) gravitation for galactic environments,” A&A, vol. 472, no. 3, p. 833, 2007.
  * (35) S. Rahvar and Y. Sobouti, “An inverse f(R) gravitation for cosmic speed up, and dark energy equivalent,” Modern Physics Letters A, vol. 23, no. 23, pp. 1929–1937, 2008.
  * (36) R. Saffari and S. Rahvar, “f(R) gravity: From the pioneer anomaly to cosmic acceleration,” Phys. Rev. D, vol. 77, p. 104028, May 2008.
  * (37) A. Abdujabbarov, F. Atamurotov, Y. Kucukakca, B. Ahmedov, and U. Camci, “Shadow of Kerr-Taub-NUT black hole,” Astrophys. Space Sci., vol. 344, pp. 429–435, 2013.
  * (38) A. Yumoto, D. Nitta, T. Chiba, and N. Sugiyama, “Shadows of Multi-Black Holes: Analytic Exploration,” Phys. Rev. D, vol. 86, p. 103001, 2012.
  * (39) F. Atamurotov, A. Abdujabbarov, and B. Ahmedov, “Shadow of rotating non-Kerr black hole,” Phys. Rev. D, vol. 88, no. 6, p. 064004, 2013.
  * (40) A. F. Zakharov, “Constraints on a charge in the Reissner-Nordström metric for the black hole at the Galactic Center,” Phys. Rev. D, vol. 90, no. 6, p. 062007, 2014.
  * (41) U. Papnoi, F. Atamurotov, S. G. Ghosh, and B. Ahmedov, “Shadow of five-dimensional rotating Myers-Perry black hole,” Phys. Rev. D, vol. 90, no. 2, p. 024073, 2014.
  * (42) T. Johannsen, A. E. Broderick, P. M. Plewa, S. Chatzopoulos, S. S. Doeleman, F. Eisenhauer, V. L. Fish, R. Genzel, O. Gerhard, and M. D. Johnson, “Testing General Relativity with the Shadow Size of Sgr A*,” Phys. Rev. Lett., vol. 116, no. 3, p. 031101, 2016.
  * (43) T. Johannsen, “Sgr A* and General Relativity,” Class. Quant. Grav., vol. 33, no. 11, p. 113001, 2016.
  * (44) J. W. Moffat, “Modified Gravity Black Holes and their Observable Shadows,” Eur. Phys. J. C, vol. 75, no. 3, p. 130, 2015.
  * (45) S. B. Giddings and D. Psaltis, “Event Horizon Telescope Observations as Probes for Quantum Structure of Astrophysical Black Holes,” Phys. Rev. D, vol. 97, no. 8, p. 084035, 2018.
  * (46) P. V. P. Cunha, J. Grover, C. Herdeiro, E. Radu, H. Runarsson, and A. Wittig, “Chaotic lensing around boson stars and Kerr black holes with scalar hair,” Phys. Rev. D, vol. 94, no. 10, p. 104023, 2016.
  * (47) N. Tsukamoto, “Black hole shadow in an asymptotically-flat, stationary, and axisymmetric spacetime: The Kerr-Newman and rotating regular black holes,” Phys. Rev. D, vol. 97, no. 6, p. 064021, 2018.
  * (48) R. A. Hennigar, M. B. J. Poshteh, and R. B. Mann, “Shadows, Signals, and Stability in Einsteinian Cubic Gravity,” Phys. Rev. D, vol. 97, no. 6, p. 064041, 2018.
  * (49) P. V. P. Cunha and C. A. R. Herdeiro, “Shadows and strong gravitational lensing: a brief review,” Gen. Rel. Grav., vol. 50, no. 4, p. 42, 2018.
  * (50) A. Allahyari, M. Khodadi, S. Vagnozzi, and D. F. Mota, “Magnetically charged black holes from non-linear electrodynamics and the Event Horizon Telescope,” JCAP, vol. 02, p. 003, 2020.
  * (51) A. B. Abdikamalov, A. A. Abdujabbarov, D. Ayzenberg, D. Malafarina, C. Bambi, and B. Ahmedov, “Black hole mimicker hiding in the shadow: Optical properties of the $\gamma$ metric,” Phys. Rev. D, vol. 100, no. 2, p. 024014, 2019.
  * (52) A. Övgün, I. Sakallı, J. Saavedra, and C. Leiva, “Shadow cast of noncommutative black holes in Rastall gravity,” Mod. Phys. Lett. A, vol. 35, no. 20, p. 2050163, 2020.
  * (53) R. Shaikh, “Black hole shadow in a general rotating spacetime obtained through Newman-Janis algorithm,” Phys. Rev. D, vol. 100, no. 2, p. 024028, 2019.
  * (54) C. Bambi, K. Freese, S. Vagnozzi, and L. Visinelli, “Testing the rotational nature of the supermassive object M87* from the circularity and size of its first image,” Physical Review D, vol. 100, p. 044057, Aug. 2019.
  * (55) S. Vagnozzi and L. Visinelli, “Hunting for extra dimensions in the shadow of M87*,” Physical Review D, vol. 100, p. 024020, July 2019.
  * (56) R. Kumar, S. G. Ghosh, and A. Wang, “Gravitational deflection of light and shadow cast by rotating Kalb-Ramond black holes,” Phys. Rev. D, vol. 101, no. 10, p. 104001, 2020.
  * (57) P.-C. Li, M. Guo, and B. Chen, “Shadow of a spinning black hole in an expanding universe,” Phys. Rev. D, vol. 101, no. 8, pp. 1–26, 2020.
  * (58) A. Övgün and I. Sakallı, “Testing generalized Einstein–Cartan–Kibble–Sciama gravity using weak deflection angle and shadow cast,” Class. Quant. Grav., vol. 37, no. 22, p. 225003, 2020.
  * (59) M. Khodadi, A. Allahyari, S. Vagnozzi, and D. F. Mota, “Black holes with scalar hair in light of the Event Horizon Telescope,” Journal of Cosmology and Astroparticle Physics, vol. 2020, pp. 026–026, Sept. 2020.
  * (60) Z. Zhong, Z. Hu, H. Yan, M. Guo, and B. Chen, “QED effects on Kerr black hole shadows immersed in uniform magnetic fields,” Phys. Rev. D, vol. 104, no. 10, p. 104028, 2021.
  * (61) F. H. Zuluaga and L. A. Sánchez, “Accretion disk around a Schwarzschild black hole in asymptotic safety,” Eur. Phys. J. C, vol. 81, no. 9, p. 840, 2021.
  * (62) O. S. Stashko, V. I. Zhdanov, and A. N. Alexandrov, “Thin accretion discs around spherically symmetric configurations with nonlinear scalar fields,” Phys. Rev. D, vol. 104, no. 10, p. 104055, 2021.
  * (63) F. Rahaman, T. Manna, R. Shaikh, S. Aktar, M. Mondal, and B. Samanta, “Thin accretion disks around traversable wormholes,” Nucl. Phys. B, vol. 972, p. 115548, 2021.
  * (64) A. Övgün, “Black hole with confining electric potential in scalar-tensor description of regularized 4-dimensional Einstein-Gauss-Bonnet gravity,” Phys. Lett. B, vol. 820, p. 136517, 2021.
  * (65) R. C. Pantig and A. Övgün, “Dehnen halo effect on a black hole in an ultra-faint dwarf galaxy,” JCAP, vol. 08, no. 08, p. 056, 2022.
  * (66) G. S. Bisnovatyi-Kogan and O. Y. Tsupko, “Analytical study of higher-order ring images of the accretion disk around a black hole,” Phys. Rev. D, vol. 105, no. 6, p. 064040, 2022.
  * (67) S. Kazempour, Y.-C. Zou, and A. R. Akbarieh, “Analysis of accretion disk around a black hole in dRGT massive gravity,” Eur. Phys. J. C, vol. 82, no. 3, p. 190, 2022.
  * (68) R. C. Pantig and A. Övgün, “Testing dynamical torsion effects on the charged black hole’s shadow, deflection angle and greybody with M87* and Sgr. A* from EHT,” Annals Phys., vol. 448, p. 169197, 2023.
  * (69) R. Roy, S. Vagnozzi, and L. Visinelli, “Superradiance evolution of black hole shadows revisited,” Physical Review D, vol. 105, p. 083002, Apr. 2022.
  * (70) Y. Chen, R. Roy, S. Vagnozzi, and L. Visinelli, “Superradiant evolution of the shadow and photon ring of sgr a*,” Physical Review D, vol. 106, p. 043021, Aug. 2022.
  * (71) S. Vagnozzi, R. Roy, Y.-D. Tsai, L. Visinelli, M. Afrin, A. Allahyari, P. Bambhaniya, D. Dey, S. G. Ghosh, P. S. Joshi, K. Jusufi, M. Khodadi, R. K. Walia, A. Övg”un, and C. Bambi, “Horizon-scale tests of gravity theories and fundamental physics from the Event Horizon Telescope image of Sagittarius A*,” Classical and Quantum Gravity, vol. 40, p. 165007, Aug. 2023.
  * (72) E. Omwoyo, H. Belich, J. C. Fabris, and H. Velten, “Remarks on the black hole shadows in Kerr-de Sitter space times,” The European Physical Journal C, vol. 82, p. 395, May 2022.
  * (73) J. P. Luminet, “Image of a spherical black hole with thin accretion disk,” Astron. Astrophys., vol. 75, pp. 228–235, 1979.
  * (74) N. I. Shakura and R. A. Sunyaev, “Black holes in binary systems. Observational appearance.,” A&A, vol. 24, pp. 337–355, Jan. 1973.
  * (75) I. D. Novikov and K. S. Thorne, “Astrophysics of black holes.,” in Black holes. Les Astres occlus. Edited by C. DeWitt and B. S. DeWitt, pp. 343–450, Jan. 1973.
  * (76) D. N. Page and K. S. Thorne, “Disk-Accretion onto a Black Hole. Time-Averaged Structure of Accretion Disk,” The Astrophysical Journal, vol. 191, p. 499, July 1974.
  * (77) 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.
  * (78) M. Guerrero, G. J. Olmo, D. Rubiera-Garcia, and D. S.-C. Gómez, “Shadows and optical appearance of black bounces illuminated by a thin accretion disk,” Journal of Cosmology and Astroparticle Physics, vol. 2021, p. 036, aug 2021.
  * (79) G.-P. Li and K.-J. He, “Observational appearances of a f(R) global monopole black hole illuminated by various accretions,” The European Physical Journal C, vol. 81, p. 1018, Nov. 2021.
  * (80) M. Okyay and A. Övgün, “Nonlinear electrodynamics effects on the black hole shadow, deflection angle, quasinormal modes and greybody factors,” Journal of Cosmology and Astroparticle Physics, vol. 2022, p. 009, Jan. 2022.
  * (81) M. Guerrero, G. J. Olmo, D. Rubiera-Garcia, and D. S.-C. Gómez, “Light ring images of double photon spheres in black hole and wormhole spacetimes,” Phys. Rev. D, vol. 105, p. 084057, Apr 2022.
  * (82) S. Hu, C. Deng, D. Li, X. Wu, and E. Liang, “Observational signatures of Schwarzschild-MOG black holes in scalar-tensor-vector gravity: shadows and rings with different accretions,” The European Physical Journal C, vol. 82, p. 885, Oct. 2022.
  * (83) M. Guerrero, G. J. Olmo, D. Rubiera-Garcia, and D. S.-C. Gómez, “Multiring images of thin accretion disk of a regular naked compact object,” Phys. Rev. D, vol. 106, p. 044070, Aug 2022.
  * (84) S. Guo, G.-R. Li, and E.-W. Liang, “Observable characteristics of the charged black hole surrounded by thin disk accretion in rastall gravity,” Classical and Quantum Gravity, vol. 39, p. 135004, jun 2022.
  * (85) H.-M. Wang, Z.-C. Lin, and S.-W. Wei, “Optical appearance of einstein-Æther black hole surrounded by thin disk,” Nuclear Physics B, vol. 985, p. 116026, 2022.
  * (86) L. Chakhchi, H. El Moumni, and K. Masmar, “Shadows and optical appearance of a power-yang-mills black hole surrounded by different accretion disk profiles,” Phys. Rev. D, vol. 105, p. 064031, Mar 2022.
  * (87) A. Uniyal, R. C. Pantig, and A. Övgün, “Probing a non-linear electrodynamics black hole with thin accretion disk, shadow, and deflection angle with M87* and Sgr A* from EHT,” Physics of the Dark Universe, vol. 40, p. 101178, May 2023.
  * (88) A. Uniyal, S. Chakrabarti, R. C. Pantig, and A. Övgün, “Nonlinearly charged black holes: Shadow and Thin-accretion disk,” arXiv e-prints, p. arXiv:2303.07174, Mar. 2023.
  * (89) K. S. Virbhadra, D. Narasimha, and S. M. Chitre, “Role of the scalar field in gravitational lensing,” A&A, vol. 337, pp. 1–8, Sept. 1998.
  * (90) K. S. Virbhadra and G. F. R. Ellis, “Schwarzschild black hole lensing,” Phys. Rev. D, vol. 62, p. 084003, Sep 2000.
  * (91) V. Bozza, S. Capozziello, G. Iovane, and G. Scarpetta, “Strong Field Limit of Black Hole Gravitational Lensing,” General Relativity and Gravitation, vol. 33, pp. 1535–1548, Sept. 2001.
  * (92) K. S. Virbhadra and G. F. R. Ellis, “Gravitational lensing by naked singularities,” Physical Review D, vol. 65, p. 103004, May 2002.
  * (93) V. Bozza, “Gravitational lensing in the strong field limit,” Phys. Rev. D, vol. 66, p. 103001, Nov 2002.
  * (94) W. Hasse and V. Perlick, “Gravitational Lensing in Spherically Symmetric Static Spacetimes with Centrifugal Force Reversal,” General Relativity and Gravitation, vol. 34, pp. 415–433, Mar. 2002.
  * (95) K. S. Virbhadra and C. R. Keeton, “Time delay and magnification centroid due to gravitational lensing by black holes and naked singularities,” Phys. Rev. D, vol. 77, p. 124014, Jun 2008.
  * (96) K. S. Virbhadra, “Relativistic images of schwarzschild black hole lensing,” Phys. Rev. D, vol. 79, p. 083004, Apr 2009.
  * (97) G. He, X. Zhou, Z. Feng, X. Mu, H. Wang, W. Li, C. Pan, and W. Lin, “Gravitational deflection of massive particles in Schwarzschild-de Sitter spacetime,” The European Physical Journal C, vol. 80, p. 835, Sept. 2020.
  * (98) S. L. Adler and K. S. Virbhadra, “Cosmological constant corrections to the photon sphere and black hole shadow radii,” General Relativity and Gravitation, vol. 54, p. 93, Aug. 2022.
  * (99) K. S. Virbhadra, “Distortions of images of schwarzschild lensing,” Phys. Rev. D, vol. 106, p. 064038, Sep 2022.
  * (100) H. Hoekstra, M. Franx, K. Kuijken, R. G. Carlberg, and H. K. C. Yee, “Lensing by galaxies in CNOC2 fields,” Monthly Notices of the Royal Astronomical Society, vol. 340, pp. 609–622, 04 2003.
  * (101) E. S. Sheldon, D. E. Johnston, J. A. Frieman, R. Scranton, T. A. McKay, A. J. Connolly, T. Budavári, I. Zehavi, N. A. Bahcall, J. Brinkmann, and M. Fukugita, “The galaxy-mass correlation function measured from weak lensing in the sloan digital sky survey,” The Astronomical Journal, vol. 127, p. 2544, may 2004.
  * (102) R. Mandelbaum, U. Seljak, G. Kauffmann, C. M. Hirata, and J. Brinkmann, “Galaxy halo masses and satellite fractions from galaxy–galaxy lensing in the Sloan Digital Sky Survey: stellar mass, luminosity, morphology and environment dependencies,” Monthly Notices of the Royal Astronomical Society, vol. 368, pp. 715–731, 04 2006.
  * (103) R. Gavazzi, T. Treu, J. D. Rhodes, L. V. E. Koopmans, A. S. Bolton, S. Burles, R. J. Massey, and L. A. Moustakas, “The sloan lens acs survey. iv. the mass density profile of early-type galaxies out to 100 effective radii,” The Astrophysical Journal, vol. 667, p. 176, sep 2007.
  * (104) L. C. Parker, H. Hoekstra, M. J. Hudson, L. van Waerbeke, and Y. Mellier, “The masses and shapes of dark matter halos from galaxy-galaxy lensing in the cfht legacy survey,” The Astrophysical Journal, vol. 669, p. 21, nov 2007.
  * (105) M. Carmeli and T. Kuzmenko, “Value of the cosmological constant: Theory versus experiment,” AIP Conference Proceedings, vol. 586, no. 1, pp. 316–318, 2001.
  * (106) P. D. Mannheim and D. Kazanas, “Exact vacuum solution to conformal Weyl gravity and galactic rotation curves,” The Astrophysical Journal, vol. 342, p. 635, July 1989.
  * (107) A. Aragón, P. A. González, E. Papantonopoulos, and Y. Vásquez, “Quasinormal modes and their anomalous behavior for black holes in $f(R)$ gravity,” Eur. Phys. J. C, vol. 81, no. 5, p. 407, 2021.
  * (108) K. Malekmakan, R. Saffari, S. Soroushfar, and B. Pourhassan, “Gravitation field perturbation quasinormal modes of a black hole in F(R) gravity,” Eur. Phys. J. Plus, vol. 137, no. 11, p. 1255, 2022.
  * (109) S. Soroushfar and S. Upadhyay, “Accretion disks around a static black hole in $f(R)$ gravity,” Eur. Phys. J. Plus, vol. 135, no. 3, p. 338, 2020.
  * (110) S. Upadhyay, S. Soroushfar, and R. Saffari, “Perturbed thermodynamics and thermodynamic geometry of a static black hole in f(R) gravity,” Mod. Phys. Lett. A, vol. 36, no. 29, p. 2150212, 2021.
  * (111) B. Majeed and M. Jamil, “Dynamics and center of mass energy of colliding particles around black hole in $f(R)$ gravity,” Int. J. Mod. Phys. D, vol. 26, no. 5, p. 1741017, 2017.
  * (112) J. L. Synge, “The Escape of Photons from Gravitationally Intense Stars,” Monthly Notices of the Royal Astronomical Society, vol. 131, pp. 463–466, 02 1966.
  * (113) Y. B. Zel’dovich and I. D. Novikov, “RELATIVISTIC ASTROPHYSICS. II,” Soviet Physics Uspekhi, vol. 8, pp. 522–577, Apr. 1966.
  * (114) J. Bardeen, “Timelike and null geodesics in the Kerr metric,” in Les Houches Summer School of Theoretical Physics: Black Holes, pp. 215–240, 1973.
  * (115) S. Chandrasekhar, The mathematical theory of black holes.  Oxford classic texts in the physical sciences, Oxford University Press, 1998.
  * (116) J.-P. Luminet, “Seeing Black Holes: From the Computer to the Telescope,” Universe, vol. 4, p. 86, Aug. 2018.
  * (117) H. Falcke, F. Melia, and E. Agol, “Viewing the shadow of the black hole at the galactic center,” The Astrophysical Journal, vol. 528, p. L13, dec 1999.
  * (118) T. Johannsen and D. Psaltis, “TESTING THE NO-HAIR THEOREM WITH OBSERVATIONS IN THE ELECTROMAGNETIC SPECTRUM. II. BLACK HOLE IMAGES,” The Astrophysical Journal, vol. 718, pp. 446–454, July 2010.
  * (119) M. D. Johnson, A. Lupsasca, A. Strominger, G. N. Wong, S. Hadar, D. Kapec, R. Narayan, A. Chael, C. F. Gammie, P. Galison, D. C. M. Palumbo, S. S. Doeleman, L. Blackburn, M. Wielgus, D. W. Pesce, J. R. Farah, and J. M. Moran, “Universal interferometric signatures of a black hole’s photon ring,” Science Advances, vol. 6, p. eaaz1310, Mar. 2020.
  * (120) V. Perlick and O. Y. Tsupko, “Calculating black hole shadows: Review of analytical studies,” Physics Reports, vol. 947, pp. 1–39, Feb. 2022.
  * (121) C. Bambi, K. Freese, S. Vagnozzi, and L. Visinelli, “Testing the rotational nature of the supermassive object M87* from the circularity and size of its first image,” Phys. Rev. D, vol. 100, no. 4, p. 044057, 2019.
  * (122) S. Soroushfar, R. Saffari, J. Kunz, and C. Lämmerzahl, “Analytical solutions of the geodesic equation in the spacetime of a black hole in f(R) gravity,” Physical Review D, vol. 92, p. 044010, Aug. 2015.
  * (123) P. F. Byrd and M. D. Friedman, Handbook of Elliptic Integrals for Engineers and Scientists.  Berlin, Heidelberg: Springer Berlin Heidelberg, 1971.
  * (124) C. W. Misner, K. S. Thorne, J. A. Wheeler, and D. Kaiser, Gravitation.  Princeton, N.J: Princeton University Press, 2017.  OCLC: on1006427790.
  * (125) S. Guo, G.-R. Li, and E.-W. Liang, “Influence of accretion flow and magnetic charge on the observed shadows and rings of the hayward black hole,” Phys. Rev. D, vol. 105, p. 023024, Jan 2022.
  * (126) P. S. Joshi, D. Malafarina, and R. Narayan, “Distinguishing black holes from naked singularities through their accretion disc properties,” Classical and Quantum Gravity, vol. 31, p. 015002, nov 2013.
  * (127) B. C. Bromley, K. Chen, and W. A. Miller, “Line Emission from an Accretion Disk around a Rotating Black Hole: Toward a Measurement of Frame Dragging,” The Astrophysical Journal, vol. 475, pp. 57–64, Jan. 1997.
  * (128) X.-X. Zeng and H.-Q. Zhang, “Influence of quintessence dark energy on the shadow of black hole,” The European Physical Journal C, vol. 80, p. 1058, Nov. 2020.
  * (129) C. Bambi, “Can the supermassive objects at the centers of galaxies be traversable wormholes? The first test of strong gravity for mm/sub-mm very long baseline interferometry facilities,” Physical Review D, vol. 87, p. 107501, May 2013.
  * (130) A. He, J. Tao, Y. Xue, and L. Zhang, “Shadow and photon sphere of black hole in clouds of strings and quintessence *,” Chinese Physics C, vol. 46, p. 065102, June 2022.
  * (131) P. S. Letelier, “Clouds of strings in general relativity,” Phys. Rev. D, vol. 20, pp. 1294–1302, Sep 1979.
  * (132) G.-P. Li and K.-J. He, “Observational appearances of a f(R) global monopole black hole illuminated by various accretions,” Eur. Phys. J. C, vol. 81, no. 11, p. 1018, 2021.
  * (133) 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.