Exploring departures from Schwarzschild black hole in f(R) gravity

Different astrophysical methods can be combined to detect possible deviations from General Relativity. In this work, we consider a class of f(R) gravity models selected by the existence of Noether symmetries. In this framework, it is possible to determine a set of static and spherically symmetric black hole solutions, encompassing small departures from the Schwarzschild geometry. In particular, when gravity is the only dominating interaction, we exploit the ray-tracing technique to reconstruct the image of a black hole, the epicyclic frequencies, and the black hole shadow profile. Moreover, when matter dynamics is also affected by an electromagnetic radiation force, we take into account the general relativistic Poynting–Robertson effect. In light of the obtained results, the proposed strategy results to be robust and efficient: on the one hand, it allows to investigate gravity from strong to weak field regimes; on the other hand, it is capable of detecting small departures from General Relativity, depending on the current observational sensitivity.


Introduction
The recent breakthrough discoveries represented by the direct detection of gravitational waves [1,2], as well as the capacity to map the matter motion in the close vicinity of supermassive black holes (BHs) [3][4][5] have provided a large amount of complementary and sensitive observational data to inquire more deeply selfgravitating systems and gravitational interaction.General Relativity (GR) successfully passed several observational tests (see e.g., Refs [6][7][8]), still confirming as the most reliable and probed theory of gravity so far available.Nevertheless, from both theoretical and observational points of view, there are critical open questions that GR is not able to address, such as the dark energy problem [9][10][11][12][13][14][15][16], the dark matter issue [17,18], the unification of gravity with the other fundamental interactions [19] towards a full quantum description of gravity [20].To alleviate such shortcomings, alternatives or extensions of Einstein's theory have been proposed with the aim to relax some main assumptions of GR or to extend it in some way [21][22][23][24][25][26][27][28].However, none of the extensions or alternative models so far proposed has been proven to address all the aforementioned issues at once.It is widely believed that the strong gravitational field in the close vicinity of a BH, which has been not yet thoroughly investigated, may contain the key to uncover new physics and definitely confirm GR, or disregard it in favour of some new proposal [5,8,29].Therefore, it becomes of utmost importance to improve our insight for comparing theoretical results with observational data.Along this line of thought, some research programs resort to general BH parametrizations, where a theory-agnostic procedure is pursued [30][31][32].
This approach takes into account generic metrics capable of reproducing BH geometries belonging to different gravity theories and being functions of a restricted set of parameters.This method aims at detecting departures from GR metrics, which are usually constrained through the fit of the observational data [33][34][35].It consists of two important steps: p1q determining, through some astrophysical phenomena, whether there are significant deviations from GR; p2q if the first point is accomplished, a strategy has to be developed in order to reconstruct, from the observational data, the most suitable BH solutions, which, in turn, may provide indirect tests of gravity.
In this paper, we combine different astrophysical methods to accurately determine deviations from the Schwarzschild metric.For this purpose, we consider f pRq gravity as a generalization of the Hilbert-Einstein action including a generic function of the Ricci scalar, R [36][37][38][39][40][41][42][43].In particular, we exploit the Noether symmetry approach [44][45][46][47][48][49][50] to select f pRq power-law models such as f pRq " R k , with k P R, so that GR is recovered as soon as k Ñ 1.The same model has been considered in Ref. [51], within the cosmological context, whereas here we extend the analysis to astrophysical scales.Therefore, we search for static and spherically symmetric BH solutions to be probed by astrophysical methods.
As a first approximation, non-gravitational effects around compact objects can be considered quite small with respect to the gravitational field [52,53].In this scenario, our strategy relies on the ray-tracing technique to reconstruct the BH image in the observer plane [3,5,54], the epicyclic frequencies [55][56][57], and the BH shadow profile [58,59].
On the other hand, when high-energy electromagnetic processes occur around a BH (e.g. the emission from a hot corona [60]), we exploit the relativistic Poynting-Robertson (PR) effect [61][62][63][64].In this context, three forces come into the game: the gravitational pull directed inwards the BH, the radiation pressure exerted outwards, and the PR radiation drag force [61,62].When radiation is absorbed by matter, treated as a test particle, the latter re-emits isotropically in its own rest frame.This produces a back-reaction force, which efficiently removes energy and angular momentum [61,62,65,66].The relativistic PR mechanism has been usually invoked to solve some puzzling issues on the dynamics of accretion disks when invested by type-I Xray bursts [67][68][69].In addition, such an effect has been extensively studied in GR both in two and three dimensions [63,64,[70][71][72][73][74][75], featuring also chaotic behaviours under particular parameter configurations [76,77].Depending on the intensity of the electromagnetic field, if the radiation pressure dominates over the attracting forces, the test particles escape to infinity; otherwise, the PR drag force drives the test particle towards a critical hypersurface, which is located outside the BH event horizon [78].Such a hypersurface is a region where all forces balance and it is characterized by stable motions of test particles [65].This particular feature may represent a fundamental tool to inquiry the properties of the underlying geometry [79,80].
The present paper is organized as follows.In Sect.2, we introduce f pRq gravity, from which we select a class of models via the Noether symmetries.Assuming a static and spherically symmetric spacetime, we single out specific BH solutions (see [50] for details).In Sect.3, we analyse deviations from Schwarzschild geometry when only gravitational interaction is present.Moreover, in Sect.4, departures from GR are studied when the electromagnetic radiation force interacts with the surrounding matter.Finally, in Sect.5, we discuss the obtained results and outline future perspectives.

Static and spherically symmetric solutions in f pRq gravity
In this section, we select particular f pRq gravity models via the Noether symmetry approach.We look for static and spherically symmetric metrics and determine a particular BH solution via astrophysical considerations.

f pRq gravity
In f pRq gravity, the Hilbert-Einstein action is generalized as follows1 [38,81]: with g being the determinant of the metric tensor g µν and L m the matter Lagrangian density.By varying the action (1) with respect to g µν , one obtains the following field equations where f R " Bf BR .Here G µν " R µν ´1 2 g µν R is the Einstein tensor, ∇ µ is the covariant derivative, l " g µν ∇ µ ∇ ν is the d'Alembert operator, and T µν the energymomentum tensor of matter fields, satisfying the conservation law ∇ ν T µν " 0. It is worth noticing that, for f pRq " R, GR is fully recovered.

The Noether symmetry approach
Among all f pRq extensions of GR, viable models can be selected using the so-called Noether symmetry approach [45-47, 82, 83].This is a criterion, based on the Noether theorem, which allows to identify models endowed with symmetries.The method consists of assuming the existence of a point transformation, which leaves the point-like Lagrangian invariant and, afterwards, in selecting the related generator X , depending on the configuration variables q i : This can be also written in terms of the Noether vector X as with η i being unknown functions of q i .The first prolongation of X, including the first derivative transformation, is where ξ µ is the infinitesimal generator related to the spacetime coordinate transformations.According to the Noether theorem, the condition implies that the point transformation is a symmetry for the Lagrangian and the symmetry generator X automatically leads to the following integral of motion: In Eqs. ( 7) and ( 8), L is the Lagrangian density depending on the spacetime coordinates x µ , on the generalized variables q i and their derivatives B µ q i , namely L " L pq i , B µ q i , x µ q, while g µ " g µ px µ , q i q are gauge functions.
Considering the f pRq gravity Lagrangian, let us search for spherically symmetric solutions, whose line element, expressed in geometric units and spherical coordinates x µ " pt, r, θ, φq, reads where dΩ 2 " dθ 2 `sin 2 θdφ 2 , MprqdΩ 2 is a general 2sphere, while νprq and λprq are unknown smooth functions.In particular, we consider the case where the Ricci scalar depends only on the radius, i.e.R " Rprq, and the Birkhoff theorem holds.These assumptions thus lead to the static metric (9) (see [81] for details).Therefore, in f pRq gravity, the metric (9) leads to the point-like Lagrangian [84,85] L " e ν´λ 2 where the prime denotes the derivative with respect to the radial coordinate r.Notice that, though the configuration space is four-dimensional, the tangent space T Q contains only 7 variables: Consequently, from the Euler-Lagrange equation with respect to ν, we obtain the following relation: The other two Euler-Lagrange equations, calculated with respect to λ and M, provide the two remaining field equations: It is worth pointing out that the introduction of the function Mprq allows obtaining a non-dissipative Lagrangian so that the energy condition can be used to recover the tt component of the Einstein field equations, namely the Euler-Lagrange equation with respect to ν.More precisely, the energy condition is a requirement of zero energy that can be applied to non-dissipative Lagrangians in order to make the Lagrange multiplier method equivalent to the variational approach.It generally reads with q i being variables in the configuration space and with the "dot" denoting the derivative with respect to an affine parameter.
The Noether symmetry existence condition yields a system of 6 partial differential equations, which can be reduced after neglecting linear combinations.In our case, the Noether theorem provides a unique solution for the function f , namely f pRq " f 0 R k with k P R [85][86][87].As a consequence, the symmetry generator X and the conserved quantity j turn out to be A Schwarzschild-de Sitter solution is recovered for constant R. It is worth noticing that, for power-law models as f pRq " R k , dynamics can be reduced and solved both at cosmological [36,42,88,89] and astrophysical [90][91][92] scales.

Black hole solutions
The integral of motion j in Eq. ( 15) can be adopted to determine two relations among ν, M, and R, which can be distinguished for k ‰ 2, -, (16) and for k " 2, e ν " ´j 12α 0 r2 p4 `r2 Rq RR 1 . ( Specific subcases of the above general solutions are discussed in Refs.[51,85,90,93], where the first-order approximation of the field equations around k " 1 is taken into account.Hence, setting k " 1 `ϵ with |ϵ| ! 1, the metric components can be written as "" p1 ´ϵq 2 p1 ´2ϵ `4ϵ 2 q r1 ´2ϵp1 `ϵqs Let us consider the generic static and spherically symmetric geometry (9), where we choose Mprq " r 2 , corresponding to standard 2-spheres.This class of metrics has been extensively studied in Refs.[85,90,93].
In order to select a specific BH solution, we need to determine the parameters ϵ and C 1 such that Eq. ( 9) verifies the BH features, namely: piq the existence of a physical singularity, piiq the presence of an event horizon, piiiq the asymptotic flatness.
To avoid divergences when r Ñ 8, the exponents of r in Eqs.(18a) and (18b) give rise to the following system of inequalities 2ϵp2ϵ `1q from which we eventually obtain The case ϵ " 0 corresponds to the Schwarzschild geometry.Taking into account the interval (20), we can choose ϵ " ´0.01, for the purpose to investigate small departures from the Schwarzschild solution.Moreover, we can set C 1 " ´r1.01 0 , where r 0 is a real constant with dimension of length.Then, the metric components (18a) and (18b) become e νprq " 1 r 0.02 ´r1.01 0 r 1.03 , e λprq " 1.02 The parameter r 0 may be determined by imposing the photonsphere to be located at r ps " 3.00M , as in the Schwarzschild case, with M being the BH mass.In this case, we ensure our solution mimics the Schwarzschild properties with the following benefits: piq computations can be easily performed; piiq physically viable results can be provided.Considering the null geodesic equation, we obtain (see Ref. [94], for details) where E p and L p are the conserved energy and angular momentum along the photon trajectory, respectively.Imposing the condition for the orbit stability, dVpprq dr " 0, where V p prq is the photon geodesic potential and substituting r " r ps into Eq.( 22), we find r 0 " 2.01M .Therefore, we get where we have set M " 1 for simplicity.In this case, the distance between the BH and the observer is identified by the parameter r obs , which we choose to be r obs " 10 10 M as the distance to Sgr A ‹ , the BH at the center of our Galaxy 2 [5], which gives lim rÑr obs e νprq " 0.64, lim The BH event horizon r h is obtained by g tt " 0 and 1{g rr " 0, which leads to r h " r 0 " 2.01M .The metric (23) presents also a physical singularity located at r " 0.
3 Probing the BH solution in presence of gravity only Let us investigate the above static and spherically symmetric BH solution ( 23) via the comparison with the Schwarzschild spacetime.In both cases, gravity is the only acting force.The astrophysical techniques we are going to take into account are the following: A. the ray-tracing method for reconstructing the BH image in the observer screen; B. the epicyclic frequencies; C. the BH shadow profile.

Black hole imaging
Let us first consider the photon four-momentum k µ in the equatorial plane: Defining b :" L p {E p , as the photon impact parameter, we introduce the light bending angle ψpr, bq as [94] ψpr, bq :" where kµ " k µ {E p and R is the emission radius.The above formula describes how a photon geodesic bends along the path connecting the emission region, nearby the BH, with the observer, located at a very large distance from the gravitational source.Therefore, it is possible to ray-trace the BH environment with respect to a static observer (SO) located at distance r obs and inclined with respect to the axis orthogonal to the BH equatorial plane of an angle i3 [54,95].To do this, we introduce the emission angle α, which corresponds to the angle between the photon four-momentum vector k and the radial direction r on the equatorial plane.It is possible to relate the photon impact parameter to the photon emission angle α and radius R via We also define the critical photon impact parameter as b c " r ps ?
As such, for b ą b c , the light rays reach the observer, whereas, for b ď b c , the photons are swallowed by the BH.Another useful quantity is the maximum emission angle α max , defined (for R ą r ps ) as It is worth noticing that Eq. ( 26) is defined for α P r0, α p " π{2s.On the other hand, for α P rα p , α max s, turning points occur and a symmetrization process is needed.Thus, we first consider α S " π ´α and apply Eq. ( 26) to calculate ψ S .To correctly associate the right bending angle ψ, we then compute the periastron p, which is the minimum distance of the light ray from the BH, whose value is obtained by solving the following algebraic equation [94]: which depends on b and therefore on the emission angle α.We then calculate the periastron bending angle ψ p as where b p " p{ ?e νppq .Finally, the light bending angle ψ, related to α, is given by ψ " 2ψ p ´ψS [95].It is worth observing that, in our case, the upper bounds of integration in Eqs. ( 26) and ( 31) are replaced by r obs .
In order to reconstruct the image of matter moving around a BH in the observer plane, we adopt the celestial coordinates px obs , y obs q, defined as where φ P r0, 2πs is the azimuthal angle.For each φ, we associate the related ψ using cos ψ " sin i cos φ and finally we calculate α via interpolation.We consider the presence of an accretion disk around the BH and assume, for the sake of simplicity, that matter moves on Keplerian orbits at different radii.Therefore, the matter velocity in the disk frame is given by where is the Keplerian angular velocity and the prime denotes the derivative with respect to R. In our case, the innermost stable circular orbit (ISCO) radius associated to test particles is r isco " 6.23M4 .
We then define the gravitational redshift through the relation [94,95] p1 `zq ´1 :" We consider an accretion disk around a BH, which extends from r in " r isco " 6.23M to r out " 100M .In Fig. 1, we provide a visual scheme of what described above.
In Fig. 2, we show the BH image related to the Schwarzschild (ϵ " 0) and our BH (ϵ " ´0.01) solution.We note that the two images appear as identical at the first sight.In fact, a more detailed analysis must be performed in order to spot the differences Fig. 2 Accretion disk dynamics around a BH in coordinates px obs , y obs q for a distant observer with an inclination angle i " 80 ˝.The disk extends from r in " 6.23M to rout " 100M .Matter moves with Keplerian angular velocity Ω K .For each point, we report the gravitational redshift (see lateral bars).
between the two spacetimes not only in terms of geometric structure, but also through the gravitational redshift values.To this end, we develop a further investigation in Fig. 3.We first study the shape of the disk, which is closely related to the background spacetime geometry.In such a way, we can visually distinguish between the regions, where departures from the Schwarzschild metric are present.The simulations have been performed using the same number of points and the same radii to shape the whole disk.Therefore, to quantify deviations from the Schwarzschild geometry, we can calculate the relative error on each disk-point 5 .From this analysis, we can see that strong departures can be detected moving closer to r out .This is due to the fact that our BH solution does not properly go to infinity.The minimum relative error, " 1.5%, occurs at the ISCO radius, whereas the maximum, " 4.5%, is in correspondence of r out .The analysis of the gravitational redshift leads to the same conclusions, since the relative error increases when approaching the outer edge of the disk.The outcomes of our analysis can be compared to the Event Horizon Telescope (EHT) data, whose angular resolution is pλ{Dq " p20 ´25q µas [96].Indeed, if we consider, for example, the Sgr A ‹ source with mass M " 4 ˆ10 6 M d and distance to the Earth d " 8 kpc, we find that the ISCO radius is 30.16µas, being therefore in the field of view of EHT.However, the aforementioned departures from GR cannot be currently appreciated [97] as shorter wavelengths would be needed.However, these upgrades may be attained with one of the next releases.
In this astrophysical context, the observational Xray data, consisting of the persistent flux (or equivalently luminosity) emitted by the accretion disk, are provided by the present INTEGRAL, XMM-Newton, Swift instruments [98], and by near-future space missions, like Athena [99] and e-XTP [100].Modulo the emission properties of the disk, the flux generally depends on the gravitational redshift and the solid angle, which takes into account the light bending [95].In our approach, we separately examine the fundamental contributions (light bending and gravitational redshift) that, generally, should be inferred from the data reduction and analysis.
Alternative approaches aim to fit observational data (related to almost static or very slowly rotating compact objects) by making use of the Schwarzschild metric.If any deviations is detected, other metrics may be employed in the analysis [101][102][103][104][105] along with suitable reconstruction methodologies [57,80,106,107].Further alternative strategies involve the use of general BH parametrizations to be fitted with data, in order to trace back the underlying geometry [30-32, 108, 109]. 5The relative error used in our calculations is defined as ˇˇf Sch ´fsol fSch ˇˇ, where f Sch and f sol are quantities related to the Schwarzschild and our f pRq BH solution, respectively.

Epicyclic frequencies
Another approach for investigating compact objects is represented by the epicyclic frequencies.These generally encode the gravitational perturbations on the classical circular motion along the three fundamental directions ν r , ν θ , and ν φ .However, in a spherically symmetric spacetime, we have only ν r and ν φ , since ν θ " ν φ (see Fig. 1).These quantities exhibit a series of advantageous properties [94,110], such as: strong relation with the background spacetime geometry; production of observational effects in extreme gravitational field regimes; direct measurement via quasi-periodic oscillations (QPOs) through the present (e.g., XMM-Newton [98], LOFT [111]) and near future (e.g., e-XTP [100], IXPE [112]) observations; great availability of already existing observational data, extracted by applying the most suited QPO models to the gravitational system under investigation.
A common procedure makes use of the related epicyclic angular velocities tΩ r , Ω φ u, whose expressions are [57] Ω φ :" Ω r :" where Ω φ " Ω k (cf.Eq. ( 34)).It is possible to calculate r isco by determining the r-value such that Ω r is zero, in agreement with the calculations performed by resorting to the timelike geodesic equations [94].In Fig. 4, we display the two epicyclic angular frequencies both for the Schwarzschild case (black solid lines) and for our BH solution (red dashed lines), with the related relative errors.We can clearly notice a relevant difference along the radial direction, more precisely close to the ISCO radius up to r " 7M , and also moving further from the gravitational source.Furthermore, the radial epicyclic frequency is defined up to r « 49.31M , where it becomes zero.After this value, the radius becomes imaginary.The error bands of our analysis are within the current observational sensitivity (see e.g., Ref. [113], for details).

Shadow profile
The last method we focus on takes into account the BH shadow profile.The presence of a BH can be inferred from the interaction with the surrounding matter through the formation of accretion structures.The BH x obs (M)  gravitational potential energy is converted into kinetic energy and matter increases its temperature due to dissipation effects occurring in the disk plasma, causing electromagnetic emission mainly in the X-ray energy band.The strong gravitational field affects the light bending, so it is possible to obtain a two-dimensional dark zone in the observer plane, known as BH shadow [58], through the ray-tracing technique.This system is constructed by considering not only first-order direct photon images, but also higher-order terms, arising from photon trajectories undergoing one or more loops in the close vicinity of the BH [54,114].Therefore, this dark region encircled by a bright emission ring, only depends upon the BH spacetime geometry (see Fig. 1).In general, the shadow profile of a non-rotating BH is a standard circle with radius b c [54], whereas, for a rotating BH, the shadow will be elongated along the rotation direction.Regarding astrophysical BHs, their size and shape closely depend on the BH mass and spin, as well as on the inclination of the observer with respect to the BH equatorial plane.Such a shadow is fundamental for extracting information on gravity in extreme field regimes.
In Fig. 5, we display the shadow profile of our BH solution (red lines) together with the Schwarzschild metric (black lines), as well as the ISCO radius and the BH event horizon.When the inner disk radius coincides with the ISCO, it is hard to spot differences both from the shadow and ISCO radius.However, if we consider the inner disk radius to be located at the proper ISCO value for each metric, some differences may emerge.Nevertheless, difficulties in detecting relevant departures may also occur in the latter case, when the observer inclination angle is i À 30 ˝. Instead, discrepancies may be measured as the observer moves closer to the equatorial plane, due to the enhanced relativistic effects [54,95].
The observational data related to the shadow profile are essentially those provided by EHT.The possibility of detecting metric departures from this analysis follows the same discussions reported in Sect.3.1.Our case sets below the actual wavelength of EHT, but near future upgrades will most likely reach the requested precision.

Discussion of the results
Let us now discuss results obtained so far.For a more quantitative comparison, in Table 1, we report the outcomes of our solution (ϵ " ´0.01) compared with the Schwarzschild metric (ϵ " 0).The first significant difference relies on the value of the metric time component at r obs , as it can be clearly seen in Fig. 3.
Furthermore, the ray-tracing process in both spacetimes could reveal metric departures, as one can immediately note from the values of α max , ψ max , ψ p , and b.Therefore, we can conclude that a careful analysis of the accretion disk image could show whether there are departures from the Schwarzschild geometry.Nevertheless, this strategy is mainly based on the time component of the metric g tt , rather than g rr .Indeed, as an example, if two metrics differ only for the g rr component, they cannot be clearly distinguished by this method [80].
The aforementioned investigation must be complemented with further analyses.As shown in Fig. 4, the radial epicyclic frequency reveals some discrepancies in the close vicinity of the ISCO radius, up to " 7M , and of course far from the gravitational source.Besides being straightforward, this procedure involves simultaneously g tt and g rr , being therefore sensible in detecting departures in both metric components.
To conclude, we draw the BH shadow profile in Fig. 5, which simply reduces to a circle with radius b c .Compared to the techniques adopted in the previous sections, however, this approach has the advantage to inquiry gravity in stronger gravitational field regimes.In the present study, the radii of the two shadows exhibit a relative gap of the order of 1%, which thus does not provide significant information.On the other hand, determining the ISCO position turns out to be extremely important, as it allows to place constraints on the metric tensor, with an associated error of the order of 5.5%.Furthermore, by increasing the inclination angle of the observer, the relativistic effects are magnified, thus permitting to viably disclose the BH metric properties.
Therefore, the present analysis, based on three complementary astrophysical methods, is very solid for inquiring about the characteristics of the spacetime geometry around a static and spherically symmetric BH.From the observational point of view, we have seen that the current data on epicyclic frequencies are already sufficiently precise for detecting metric departures, whereas the required accuracy of EHT data (involving the image of a BH and the shadow profile) needs further sensitivity upgrades.

Probing the solutions in presence of gravity and radiation
The present section is devoted to the study of metric deviations occurring in presence of both the BH gravitational field and external forces generated by electromagnetic radiation processes.If the accreting matter is assumed to be made of uncharged test particles, the relativistic PR effect can be suitably exploited.
First, we shall introduce the gravitational setup to investigate the main features of the radiation stressenergy tensor.Then, we describe the test particle motion and its interaction with the radiation field.Based on these premises, we thus write the relativistic PR equations of motion.We eventually show that the trajectories driven by the relativistic PR effect can be used to distinguish the Schwarzschild metric from other similar geometries.

Gravitational setup
We consider a static BH, whose external spacetime geometry is described by the metric (21).Setting θ " π{2, we consider the orthonormal frame for both a SO at infinity pB B B ν q µ :" δ µ ν and for a SO located around the BH.In the latter case, the future-pointing unit vector orthogonal to the spatial hypersurfaces is given by n " x obs (M) Fig. 5 Inner disk radius, shadow profile, and BH event horizon as seen by a distant SO in its own coordinate frame px obs , y obs q.
The black and red lines refer to the Schwarzschild metric and our BH solution, respectively.Left panel.The inner disk radius is set to the same value for both metrics, namely r in " 6.23M , corresponding to the ISCO radius of our BH solution.Right panel.The inner disk regions are set to the proper ISCO radii for both metrics.
Table 1 Summary of the results from our BH solution (ϵ " ´0.01) compared with the Schwarzschild metric's parameters.
Parameter Units Our BH solution p p pϵ " " " ´0.01q q q Schwarschild p p pϵ " " " 0q q q gttpr obs q -0.64 1.00 Hereafter, all quantities (such as, v α and T αβ ) associated to the SO frame will be labeled by a hat, while all the scalar quantities measured in the SO frame (such as f ) will be followed by pnq (e.g., f pnq).In the kinematic decomposition of the SO congruence, the nonzero quantities are the acceleration apnq " ∇ n n " apnq r pnqe r and the relative Lie curvature k pLieq pnq " k pLieq pnq r e r , whose general definitions are apnq r pnq " 1 ?g rr B r lnpN q, (38a) The above expressions are reported in Table 2, together with some other useful quantities.´0.02

Stress-energy tensor for radiation
The radiation field is treated as a coherent flux of photons travelling along null geodesics on the background spacetime (21) and continuously hitting the test particle.The related stress-energy tensor is where I is the radiation field intensity and k is the photon four-momentum field.By decomposing k with respect to the SO frame, we obtain νpk, nq :" sin β e r `cos β e φ, with νpk, nq being the photon spatial unit relative velocity with respect to the SO, and β is the angle measured in the SO frame in the azimuthal direction.The radiation field is governed by the impact parameter b, associated to the emission angle β.The photon energy Epnq and the photon angular momentum along the zaxis (orthogonal to the equatorial plane) L ẑ pnq, can be expressed in the rest-frame of a distant observer as [63,64] Epnq :" where E :" ´kt ą 0, L z :" k φ , and b :" L z {E.The photon impact parameter b can be associated to its relative angle β in the SO frame as follows (cf.Eq. (41b)): For b " 0 (respectively, b ‰ 0) we model a radial (respectively, general) radiation field.The photon impact parameter b can be written in terms of physical quantities, such as the radius R ‹ and the angular velocity Ω ‹ of the emitting surface, namely [73] b with the maximum angular velocity given by Notice that the condition Ω ą Ω max is not physically viable, as it leads to superluminal velocities.From the conservation equation ∇ µ pI 2 ?´gq " 0, it is possible to determine the intensity parameter as [63,64]

Test particle motion
Let us now consider a test particle moving in the equatorial plane with timelike four-velocity U and spatial velocity νpU, nq with respect to the SO frame given by, respectively, ν :" νpsin α e r `cos α e φq.
In the above equations, γ " γpU, nq " 1{ a 1 ´||νpU, nq|| 2 is the Lorentz factor, ν " ||νpU, nq|| is the magnitude of the test particle spatial velocity νpU, nq, and α is the azimuthal angle of the vector νpU, nq measured in the SO frame clockwise from the positive direction in the r ´φ tangent plane (see Fig. 1).The explicit expression for the test particle velocity components with respect to the SO are where τ is the proper time parameter along U.

Radiation force field
We assume that the radiation-test particle interaction occurs through an elastic Thomson-like scattering, characterized by a constant momentum-transfer cross section σ, which is assumed to be independent of the direction and frequency of the radiation field.Therefore, the radiation force is [63,64] F pradq pU q α :" where P pU q α β " δ α β `U α U β projects a given vector in the spatial hypersurface orthogonal to U .Decomposing the photon four-momentum k with respect to the test particle four-velocity U and considering the SO frame, we obtain the relation which, substituted into Eq.( 49), leads to The equations of motion for test particles are mapU q " F pradq pU q, with m being the test particle mass.By defining σ " σ{m, we obtain apU q " σΦ 2 EpU q 2 Vpk, U q.
The components of Vpk, U q " Vt n `V r e r `V φ e φ are Vt " γν V r " sin β γr1 ´ν cospα ´βqs ´γν sin α, The normalized luminosity of the radiation field is given by A{M " σI 2 0 E 2 " L 8 {L Edd and lies within the range r0, 1s.Here, L 8 and L Edd :" 4πM m{σ are the luminosity evaluated at infinity and the Eddington luminosity, respectively.Using Eq. ( 53), the term σrIEpU qs 2 can be recast as

Critical surfaces
The set of Eqs.(55a)-(55d) admits a region on which gravitational attraction, radiation pressure, and PR drag force balance.In the two-dimensional case, this region gives rise to a circle dubbed critical surface, only characterized by its radius r pcritq .Therefore, let us consider a test particle in radial equilibrium, moving in a purely circular motion (i.e., α " 0 and ν " ν pcritq ).Then, setting dν dτ " 0, Eq. (55a) becomes from which we obtain ν pcritq " cos β.Thus, the velocity of the orbiting test particle is equal to the azimuthal velocity of the radiation field photons.It is worth noticing that, if the photons are radially emitted (that is b " 0), then ν pcritq " 0 and the test particle will stop on a point on the critical surface.On the other hand, for b ‰ 0, the test particle will move on the critical surface with constant velocity ν pcritq , being continuously pushed in the orbiting direction by the "non-radial" photons of the radiation field.If we set α " π in Eq. (55a), we obtain the same velocity ν pcritq , but with a negative sign.In the SO frame, the test particle's velocity is tangential to the critical surface (i.e., dα dτ " 0), so that Eq. (55b) becomes apnq r `kpLieq pnq r ν 2 pcritq " Therefore, once the values of the parameters tA, R ‹ , Ω ‹ u are assigned, the size of the critical surface can be obtained by solving Eq. ( 58) in the r-coordinate.

Trajectories of test particles
Let us now take into account the test particle trajectories with viable values for the underlying relativistic PR variables.In order to reduce the parameter space, we fix initial conditions on the test particle motion by means of physical considerations.By examining the outcomes, we can infer information on how to detect deviations from the Schwarzschild geometry.An important feature of the relativistic PR effect, which will be extensively used, relies on its sensitive dependence on the initial conditions [76].
To develop the graphs in Fig. 6, test particles in all simulations are assumed to start from the same position, namely the ISCO radius of our BH solution and to move with Keplerian velocity.In other words, the initial conditions are fixed as pν 0 , α 0 , r 0 , φ 0 q " pν K , 0, r ISCO , 0q, with ν K " r ISCO Ω K pr ISCO q.It is worth noticing that, albeit the test particles start from the same radius, they have different velocities, because the expression of Ω K is different in the two metrics.In Fig. 7, the test particles start from the corresponding ISCO radius for both the metrics, but they always move with Keplerian velocity.Therefore, defining νK " rISCO Ω K pr ISCO q, the initial conditions are fixed as pν 0 , α 0 , r 0 , φ 0 q " pν K , 0, rISCO , 0q, with rISCO " 6M for Schwarzschild and rISCO " r ISCO for our BH solution.
Regarding the radiation field variables, we choose three values for the luminosity parameter, i.e., A " p0.1, 0.4, 0.6q, corresponding to low, low-intermediate, and high-intermediate luminosities, respectively.Luminosities higher than A " 0.7 are not included, for avoiding the radiation pressure dominating over the attracting forces, and pushing the test particle to infinity.
Concerning the emitting surface, we fix its radius7 to R ‹ " 2.5M for both metrics and let its angular velocity vary as Ω ‹ " p0, 0.09, 0.16q M ´1.The maximum angular velocities are Ω max p2.5M q " 0.179M

´1
and Ω max p2.5M q " 0.178M ´1, corresponding to the Schwarzschild and to our BH solution, respectively.Therefore, we investigate the cases of static, medium, and high rotational spins.The parameters A, R ‹ , and Ω ‹ can be inferred from X-ray observational data (see Refs. [60,98], and references therein for details).As the photon impact parameter (43) depends on the metric, we provide its expression both in the Schwarzschild geometry (b Sch ) and in terms of our BH solution (b sol ).In our setup, if the critical region is located within the emission zone, then the test particle trajectory ends its motion when meeting the outermost external surface.The discrepancies among test particle trajectories are strongly enhanced by the relativistic PR effect due to its sensitive dependence on the initial data and the underlying geometry.Therefore, such a behaviour permits to single out particular features that can be surely detected through the data provided either by EHT or GRAVITY collaborations.

Discussion on the results
We provided two examples of trajectories with different initial conditions on the test particle motion, namely: (1) same radius, but different velocity (see Fig. 6); (2) different radius and velocity (see Fig. 7).We can immediately see that both simulations exhibit akin trajectories, although, in the second set, differences appear more pronounced.For low luminosities, we note more similarities and this makes it difficult to appreciate possible deviations.In the case of low-intermediate luminosities, discrepancies are more visible for medium angular velocity.Finally, for high-intermediate luminosities, relevant differences emerge.
It is important to note that motion of test particles can be mainly inferred by resorting to the ray-tracing technique.Depending on the matter type involved in the relativistic PR effect, we have: piq accreting plasma, captured by EHT data [3,5]; piiq astrophysical objects similar to stars orbiting around Sgr A ‹ tracked by Gemini, Subaru, and VLT data provided by the GRAV-ITY collaboration (see e.g., Ref. [115], and references therein, for details).
From our simulations, intermediate luminosities, medium and high rotational angular velocities show that the considered BH geometry is different from the Schwarzschild Fig. 6 Black trajectories refer to the Schwarzschild metric, whereas red dashed curves to our BH solution.The BH is placed at the center of the coordinate axes and its black boundary accounts for the Schwarzschild radius (r h " 2M ), while the red boundary for our BH solution (r h " 2.01M ).Due to the adopted scales, such boundaries are almost indistinguishable.The green and orange lines identify the critical region and the emitting surface, respectively.In each plot, we also report the values of the luminosity parameter A and the photon impact parameters b Sch and b sol .The initial conditions on the test particle are pν 0 , α 0 , r 0 , φ 0 q " p0.40, 0, 6.23M, 0q in the Schwarzschild case and pν 0 , α 0 , r 0 , φ 0 q " p0.39, 0, 6.23M, 0q for our BH solution.
one.For instance, setting A " 0.6 and Ω ‹ " 0.16, the test particle reaches the critical surface in the case of Schwarzschild metric, while it escapes to infinity for our BH solution.

Conclusions and perspectives
In this paper, we considered f pRq theories of gravity in a static and spherically symmetric spacetime.By means of Noether symmetries, we selected viable power-law models of the form f pRq " R k .We found that the metric solutions to the field equations can be parameterized by two real constants (ϵ, r 0 ), where r 0 is a constant and ϵ accounts for deviations from the Schwarzschild spacetime, being k " 1 `ϵ, with |ϵ| ! 1. Physical considerations constrain the parameter ϵ in the range ´1{2 ď ϵ ď 0. In particular, to detect small deviations from the Schwarzschild BH, we chose ϵ " ´0.01, as a value that could lie under the current observational sensitivity.
We proposed a combination of different and known astrophysical techniques capable of detecting possible small metric departures from the Schwarzschild spacetime.In particular, we considered two astrophysical scenarios: (1) gravity is the only force acting on the surrounding matter; (2) gravity and electromagnetic radiation processes alter the matter's motion.The first case is analysed by exploiting: piq the BH image reconstruction via the ray-tracing technique; piiq the epicyclic frequencies; piiiq the BH shadow profile.Our results demonstrate that the BH image can provide a deep and complete analysis of the underlying geometry.This is due to the extension of the image from the ISCO radius to large distances from the gravitational source.
A more complete and simple approach is represented by epicyclic frequencies.In this regard, the radial frequency is capable of highlighting differences with respect to the Schwarzschild metric, from the ISCO radius up to r " 7M and far from the gravitational source.Instead, the BH shadow profile cannot provide significant Fig. 7 The notation is the same as in Fig. 6.In this case, the test particle's initial conditions in the Schwarzschild metric are pν 0 , α 0 , r 0 , φ 0 q " p0.41, 0, 6M, 0q, whereas for our BH solution they are pν 0 , α 0 , r 0 , φ 0 q " p0.39, 0, 6.23M, 0q.
information, since it depends only on the critical impact parameter that, in our case, is extremely close to the Schwarzschild geometry.However, among all methods considered here, the last approach is the only one inquiring about gravity in strong regimes.
To summarize, these three astrophysical techniques result to be complementary and efficient in investigating BH physics and gravity.Indeed, they embrace a wide class of gravitational scenarios, ranging from strong (as the photon ring radius) to weak (as the outer edge of an accretion disk) regimes.In addition, such techniques also involve g tt and g rr , thus being very sensitive to small departures from GR in both metric components.The current observational data can partially detect these metric departures, whereas the EHT data have to be upgraded in wavelength sensitivity.Furthermore, we took into account the case where electromagnetic radiation processes are modelled via the relativistic PR effect.Since this phenomenon is very sensitive to the initial conditions [76], we first considered test particles with equal initial radii but different initial velocities, and then, with different initial radii and velocities.Even in this case, it is possible to in-quire about gravity in various regimes and configurations.In particular, we found that the physical systems exhibiting more pronounced deviations from GR are those with intermediate luminosities and an emitting surface located close to the BH event horizon and rotating with either medium or high angular velocities.Exploring test particle trajectories in terms of the underlying relativistic PR parameters, it is possible to find particular configurations identifying the presence of metric changes, independently of the sensitivity of current observational EHT or GRAVITY data.Thus, our study sheds light on the possibility to investigate metric departures via the synergy among various astrophysical techniques.Such methods could be applied to already existing data or very sensitive data from the near future.The primary goal of this work was to underline the predictive power of the above astrophysical techniques, which must be contextualized into the gravitational system under study for the sensibility and availability of the related observations.Certainly, additional and alternative strategies may be used to further constrain (and thus considerably reduce the degeneracy on) deviations from the Schwarzschild metric, such as gravitational waves, quasi-normal modes, and so forth.The methodology described here could be extended to axially symmetric and stationary spacetimes around uniformly spinning BHs.In the latter situation, the BH shadow profile can be useful to provide a more rich phenomenology, thus encoding valuable information on the gravitational background in the strong field regimes.

Fig. 1
Fig.1Sketch of different astrophysical techniques exploited to detect departures from the Schwarzschild geometry.They can be divided in two astrophysical situations: (1) only gravity is acting on the surrounding matter outside the BH event horizon; (2) also X-ray radiation processes are involved, besides gravity.The first case includes the image of a BH, the epicyclic frequencies, and the BH shadow profile, whereas the second one accounts for the relativistic PR effect.

Fig. 3 Fig. 4
Fig.3In the upper plots, black lines refer to the Schwarzschild metric whereas red ones to our BH solution.Left upper panel.The two overlapped disks are individually composed of 101 rings, each one consisting of 500 points, for a total of 50500 points.Left lower panel.Relative errors between the corresponding disk-points in the two metrics.Right upper panel.Gravitational redshift versus the disk-points.Right lower panel.Relative errors related to the gravitational redshift.

1 N
B B B t , where N " ?´gtt is the lapse function.A possible set of tetrad fields adapted to the SO frame

Table 2
Physical quantities related to the BH solution (ϵ " ´0.01) compared with Schwarzschild metric's parameters.The BH mass is set to unity.