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

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.


I. INTRODUCTION AND MOTIVATION
Since the foundational concepts of black holes were established by Schwarzschild [1] and Finkelstein [2], 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 [3,4] to the recent shadow images of M87* [5] and Sgr A* [6] 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 [7][8][9].Shadow images can provide information about the geometric structure near the black hole's horizon [10] and its physical characteristics [11].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 [12][13][14][15][16][17].
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) [18][19][20].As a result, f (R) theories of gravity have garnered significant interest and have been thoroughly scrutinized for their consistency [21][22][23][24][25][26][27] (see also the reviews in Refs.[28,29]).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 [30].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 [31,32].To address these issues properly, a suitable f (R) action model has been proposed in Refs.[33][34][35], which is consistent with both galactic and cosmological scales.In Ref. [36], 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. ).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 [73].Generally, this type of accretion is based on the Shakura-Sunyaev [74], Novikov-Thorne [75], and Page-Thorne [76] 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. [77].Since then, this method has been applied in various publications (e.g., Refs.[78][79][80][81][82][83][84][85][86][87][88]), 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 [89][90][91][92][93][94][95][96][97][98][99].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 [100][101][102][103][104], 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. [77], 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 in which, κ is a coupling constant, f (R) is a function of the Ricci scalar of the spacetime with the metric determinant g, and S m is a matter field action.Accordingly, the field equations are derived as by varying the action S with respect to the metric, where F (R) = df (R) dR , □ = ∇ λ ∇ λ and T µν is the energy-momentum tensor.Adopting a generic spherically symmetric line element in the usual Schwarzschild coordinates x µ = (t, r, θ, ϕ), the field equation ( 2) can be recast as [36] 2F (r) X ′ (r) X(r) for the vacuum space (i.e.T µν = 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) ≡ 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) −α [36], where α and d are the free parameters of the action.Here α is dimensionless whereas [d] = m and is of galactic size.Applying this ansatz, Eq. ( 4) results in the solution where X 0 represents an integration constant.When α = 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 [36] B in which β = α/d.Here, Λ represents the cosmological constant having the value |Λ| ≤ 10 −52 m −2 [105] 1 .It is worth noting that the above expression resembles the Mannheim-Kazanas vacuum solution in fourth-order Weyl conformal gravity [106].In that solution, the linear term β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 α 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 [107,108], dynamics of the accretion disk [109], thermodynamic geometry [110], and particle collision in the black hole's exterior [111].
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 ≪ d), the model is derived as [36] f where R 0 = 6α 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 [36] For cosmological scales (α ≪ 1), this model reduces to f (R) = R + Λ, 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. [36], the generic model has been proposed, which aligns with both of the above models under different scale criteria.In the regime of strong curvature, characterized by R ≫ Λ and R/R 0 ≫ 2/α, we obtain the model (8).Conversely, in cosmological scales, where R ≃ R 0 ≃ Λ and α ≪ 1, it simplifies to f (R) = R + Λ.Consequently, the small-valued parameter β 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 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 where The existence of real values for the above radii, however, depends on the sign of the polynomial's discriminant, i.e. ∆ = , which is always of positive values for β, Λ ≪ 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 A. Motion of mass-less particles The motion of test particles can be described by the Lagrangian where "dot" stands for differentiation with respect to the affine parameter τ of the geodesic curves.Enjoying the spherical symmetry of the spacetime, we confine the motion of particles to the equatorial plane (i.e.θ = π/2) without loss of generality.One can then define the conjugate momenta which provides the two constants of motion 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 ≡ 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 L = 0 which characterizes the null geodesics.Thus, by means of Eq. ( 17), the equations of motion are obtained as in which represents the effective potential for photons.Then, the turning points r t in the orbits correspond to ṙ = 0, which is encountered when V (r t ) = E 2 .This potential has a maximum at given that β 0 = √ 1 + 6β, 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 Λ, and decreases as β increases.One can also verify that 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 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 [112], and the cone of gravitational radiation capture by Zeldovich and Novikov [113].Bardeen, Chandrasekhar, and Luminet popularized the terms optical appearance of black holes and black hole image to describe the visual characteristics of black holes [73,[114][115][116].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 [117].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 [118,119], although in this study, the term holds a broader significance (see Ref. [120] 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 sh = b p ). Specifically, for the f (R) black hole being examined, the expression for R sh is given by Eq. ( 25).This way, the theoretical shadow diameter for this black hole is given by d theo sh = 2R sh = 2b p .In Fig. 1(a), the β-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 [121] which calculates the shadow diameter as observed by an observer positioned at a distance D (in parsecs) from the black hole, where γ is the mass ratio of the black hole and the Sun, which is γ = (6.5 ± 0.90) × 10 9 for M87* at the distance D = 16.8Mpc [5], and is γ = (4.3±0.013)×10 6for Sgr A* at D = 8.127 kpc [6].In Eq. ( 26), θ * is the angular diameter of the shadow, which has been measured as θ * = 42 ± 3 µas for M87*, and θ * = 48.7 ± 7 µas for Sgr A*.This way, one can calculate the shadow diameters as d M87 * sh = 11 ± 1.5 and d SgrA * sh = 9.5 ± 1.4.In Fig. 1(b), these values are displayed within the 1σ uncertainties for both black holes.However, these observed diameters do not correspond to the calculated theoretical shadow diameters d theo sh 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 A).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 β for the f (R) black hole.Instead, for the numerical studies in the rest of the paper, we consider the values of the β-parameter within a specific range, marked in the shaded region of the β-profile of the theoretical shadow diameter (as seen in Fig. 1(a)).This segment confines the parameter within the range 0 ≤ β < 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. In the left panel, five radial profiles of the effective potential have been plotted which correspond to different values of the β-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 (dr/dϕ) 2 = 0. Applying the Eqs.( 21) and (22), this results in which beside the trivial solution at r = 0, has one negative root and the two positive roots 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 turning points have been given for different values of β.As it is inferred from the table, increase in the β-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 β (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. [122].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.

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 .= r i /r (r i = r d , r f ), one can rewrite the differential equation ( 27) as dz dϕ A further change of variable u .= 1 2 (z/r i − 1/6) provides us with the Weierstrassian differential equation in which are known as the Weierstrass invariants.This leads to the integrals respectively, for the OFK and OSK, in which ϕ 0 is the initial azimuth angle, and Taking into account the applied changes of variables, the above integrals yield for the OFK and for the OSK, where ℘(x) ≡ ℘(x; g2 , g3 ) is the ℘-Weierstrassian elliptic function [123], and we have defined

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 Θ that an observer O at the radial position r O from the black hole (the lens) measures.Accordingly, we have [124] with the Weierstrass invariants given in Eqs.(33).Using the analytical expression for r d , we have plotted the behavior of Θ in terms of the impact parameter b in Fig. 3.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 β-parameter.However, in general, Θ decreases for a given b, when the β-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.

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 in < r ++ , and the critical orbits of the second kind (COSK) which correspond to photons approaching r p from r + < r in < r p .When λ −2 → λ −2 p = b −2 p + Λ, the point r = r p is a double root of P 4 (r) in Eq. (27).The differential equation of motion can then be factorized as by means of the method of synthetic division, in which 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 I (ϕ) for the COFK and r II (ϕ) for the COSK, which are given as where Φ = λ −1 p (ϕ − ϕ 0 ) r p (A + r p ) + B, and In Fig. 4, some examples of the possible photon orbits have been shown.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 [74], proposed by Novikov and Thorne in Ref. [75].To proceed with applying this model, let us return to the Lagrangian (17), which is now specified as 2L = −1, for massive particles of energy E and angular momentum L that constitute the accretion disk.This way, one can rewrite the equations of motion ( 20) and ( 21) as dr dϕ , and 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 V(r) which has been plotted for the adopted values of β.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 P 6 (r) = 0 = P ′ 6 (r).Note that, although this system of equations cannot be solved by means of common radicals, however, the ISCO radius obeys the following relation [125] 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 β-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 for the energy and angular momentum of particles residing in the ISCO, in which 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 β-parameter.It is observed that by increasing β, all the quantities increase which is a consequence of the relevant changes in the effective potential.Note that, on the ISCO, one can recast the characteristic polynomial as P 6 (r) = Λ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, A r , is constant all the time, in the way that in which √ −g = r2 , Σ is the surface density of the disk, and U 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 [76,126] where F(r) is the flux of the radiated energy from the disk, and is given by By considering the fact that Now taking the expressions in Eqs. ( 48), ( 49) and ( 50) up to the first-order in Λ, and employing them in the integrand of the above relation, one can obtain the analytical solution for the flux, for which 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 in which, σ 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 β-parameter.It can be observed that by increasing β, 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. 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 A, 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.

Direct emission, lensing rings and photon rings
Here, we follow the method introduced in Ref. [77], to characterize the light rings in the observer's sky, where the number of orbits is defined as in which, ϕ 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. [77], 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 β-parameter, with distinct colors indicating the domains of direct emissions, lensing rings, and photon rings.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 Λ = 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 β-parameter, the width of the lensing rings and photon rings shrinks.In Table II, 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 β.According to the data presented in this table, it can be checked that by an increase in the β-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 β, 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.

Transfer functions and the observed intensities
The radiation of the accretion disk is supposed to be isotropic in its rest frame.By I e (r), we denote the specific intensity of an emitted radiation of frequency ν e from the disk.From Liouville's theorem, we know that the quantity I e (r)/ν 3 e is conserved along the entire path of light propagation.Hence, the observed intensity I o of frequency ν o are related in terms of the relation I e (r)/ν 3 e = I o (r)/ν 3 o [127].Accordingly, we have which in our model h = B(r).Now by integrating over the range of all the observed frequencies, the total observed specific intensity is obtained as in which, the total emission intensity is given by I emit (r) = I e (r)dν 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 where r m (b) is the transfer function that relates the impact parameter of the light ray trajectories with the radial coordinate of the mth intersection of light rays with the accretion disk2 .Hence, the slope of the transfer function indicates its (de)magnification scale [77,128].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 β.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 ≥ 4 (black hole subrings) do not possess significant observational features, although they have been found to produce certain interferometric signatures [119].

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 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 • 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 • 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 Table III 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 β-parameter considered thus far.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.
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 β-parameter, except β = 0.031, the photon ring intensity is smaller than the direct emission, while the lensing ring intensity is always larger 3 .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 β 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 β, both peaks have significant ranges and approximately equal intensities.For most cases of β, the direct emission dominates the observed intensity, while the rings are strongly demagnified.However, in the exceptional case of β = 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 β = 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 β-parameter.Additionally, the brightness of the accretion disk increases with an increase in the β-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 β = 0.022).This difference arises from changes in both the β-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.

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 [129].In this model, the observed intensity is expressed as over the null geodesic congruence γ, in which R is the redshift factor, ν e is the frequency of emitted photons from the accretion disk, dI prop is the infinitesimal proper length, and is the permittivity per unit volume in the emitter's rest frame, in which ν f is the monochromatic rest-frame emission frequency, and δ is the delta function.In this construction, the redshift factor is given by where u o and u e are, respectively, the four-velocities associated with a distant static observer, and the infalling accreting matter.Accordingly, u µ o = (1, 0, 0, 0), and in the spacetime of the f (R) black hole we can write The Π 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 Π, which yields [129] Π in which the ± signs correspond, respectively, to whether the photons approach or recede from the black hole.One can therefore recast the redshift factor as and from here, the infinitesimal proper length is obtained as Accordingly, Eq. ( 64) takes the form 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 β-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 β-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 β, 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

B. 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.[128,130] 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 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 B(r) = 1 − q − 2M/r, corresponding to a Schwazschild black hole associated with cloud of strings [131].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 B(r) = 1 − qr − 2M/r, comparable to the lapse function (7) for Λ = 0.This specific case has been studied in Ref. [128], where only the quintessential dark energy affects the spacetime, and in Ref. [130], 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 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 β-parameter.Similar considerations have been explored in Ref. [132], 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.[130,133].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 Λ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 β, along with a cosmological constant.While the cosmological constant is meant to compensate for vacuum energy and the accelerated expansion of the universe, the β-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 β.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 µas for the assumed values of the β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. [77] 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 β-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 β and the radial position of the direct emission.Furthermore, the accretion disk's brightness is enhanced with an increase in the β-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 β values.Our analysis demonstrated that as the β-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.

FIG. 1 .
FIG. 1.The diagrams show (a) the β-profile of d theo sh and (b) the observed shadow diameters d sh for M87* (green region) and Sgr A* (red region), each with 1σ uncertainties.The shaded region in panel (a) represents the assumed range for the β-parameter, limiting it within the domain 0 ≤ β < 0.05.

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

β
FIG. 3. The plot of the deflection angle Θ (in µas) versus the changes in the impact parameter b, plotted for the different values of the βparameter, and r O = 10 5 .

FIG. 4 .
FIG. 4. Examples of (a) OFK, (b) OSK, (c) COFK and (d) COSK, plotted for β = 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 rp.

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

FIG. 6 .
FIG.6.Radial profiles of (a) Ec, (b) Lc and (c) ϖc, for the adopted values of the β-parameter and the same color-coding as in Fig.5.

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

FIG. 10 .
FIG. 10.The observational characteristics of the accretion disk around the f (R) black hole, depicted for the case of β = 0 (the Schwarzschildde 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.

5 . 69 FIG. 15 .
FIG.15.Blurring the shadow in Fig.12(a) for β = 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 bp = 4.74 for the case of β = 0.022.

FIG. 16 .
FIG. 16.(a) The b-profiles of the observed intensity of the infalling spherical accretion, which from bottom to top correspond to the cases of β = 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.

TABLE I .
Table.I, the The turning points of photonic trajectories together with their corresponding values of bp.

TABLE III .
The radii of the origins for the radiation emission models, provided for various values of the β-parameter.