Magnetized discs and photon rings around Yukawa-like black holes
Alejandro Cruz-Osorio, Sergio Gimeno-Soler, José A. Font, Mariafelicia De Laurentis, Sergio Mendoza
MMagnetized discs and photon rings around Yukawa-like black holes
Alejandro Cruz-Osorio, Sergio Gimeno-Soler, Jos´e A. Font,
2, 3
Mariafelicia De Laurentis,
4, 5 and Sergio Mendoza Institut fur Theoretische Physik, Goethe Universit¨at Frankfurt,Max-von-Laue-Str.1, 60438 Frankfurt am Main, Germany Departamento de Astronom´ıa y Astrof´ısica, Universitat de Val`encia,C/ Dr. Moliner 50, 46100, Burjassot (Val`encia), Spain Observatori Astron`omic, Universitat de Val`encia,C/ Catedr´atico Jos´e Beltr´an 2, 46980, Paterna (Val`encia), Spain Dipartimento di Fisica, Universit´a di Napoli “Federico II”,Compl. Univ. di Monte S. Angelo, Edificio G, Via Cinthia, I-80126, Napoli, Italy INFN Sezione di Napoli, Compl. Univ. di Monte S. Angelo, Edificio G, Via Cinthia, I-80126, Napoli, Italy Instituto de Astronom´ıa, Universidad Nacional Aut´onoma de M´exico, AP 70-264, Ciudad de M´exico 04510, M´exico
We present stationary solutions of geometrically thick discs (or tori) endowed with a self-consistenttoroidal magnetic field distribution surrounding a non-rotating black hole in an analytical, static,spherically-symmetric f ( R )-gravity background. These f ( R )-gravity models introduce a Yukawa-like modification to the Newtonian potential, encoded in a single parameter δ which controls thestrength of the modified potential and whose specific values affect the disc configurations whencompared to the general relativistic case. Our models span different magnetic field strengths, frompurely hydrodynamical discs to highly magnetized tori. The characteristics of the solutions areidentified by analysing the central density, mass, geometrical size, angular size, and the black holemetric deviations from the Schwarzschild space-time. In the general relativistic limit ( δ = 0) ourmodels reproduce previous results for a Schwarzschild black hole. For small values of the δ parameter,corresponding to ∼
10% deviations from general relativity, we find variations of ∼
2% in the eventhorizon size, a ∼
5% shift in the location of the inner edge and center of the disc, while the outeredge increases by ∼ | δ | > .
1, however, reveals notable changes in the blackhole space-time solution which have a major impact in the morphological and thermodynamicalproperties of the discs. The comparison with general relativity is further investigated by computingthe size of the photon ring produced by a source located at infinity. This allows us to place constraintson the parameters of the f ( R )-gravity model based on the Event Horizon Telescope observations ofthe size of the light ring in M87 and SgrA ∗ . I. INTRODUCTION
Astrophysical systems comprising a rotating black holesurrounded by an accretion thick disc of plasma are rec-ognized as natural end results of highly dynamical eventsinvolving compact objects in a general-relativistic regime.Stellar-origin systems are produced in mergers of com-pact binaries comprising either a black hole and a neutronstar or two neutron stars, as well as in the gravitationalcollapse of massive stars (“failed” supernovae) [1, 2].Mergers of compact binaries have been dramatically dis-closed in recent times thanks to the LIGO-Virgo obser-vations of gravitational waves from GW170817 and thescores of multi-wavelength electromagnetic observationsthat followed [3, 4]. Additionally, black hole-disc systemsare used to explain astrophysical phenomenology of su-permassive black holes in active galactic nuclei [5, 6]. Ma-jor observational advances of the strong-gravity region ofsuch systems have recently been accomplished throughthe ground-breaking first image of the M87 black hole bythe Event Horizon Telescope [7].Theoretical models describing the morphology of sta-tionary thick discs around black holes in general relativ-ity were first developed in the seminal papers of Fish-bone and Moncrief [8] and Kozlowski et al. [9] for isen-tropic and barotropic discs respectively, assuming a con-stant distribution of angular momentum. The modelling has gradually improved through the elapsing decades [seee.g. 10, for a review]. Proposals to construct the ini-tial data for magnetized discs with weak magnetic fieldsexploring different configurations were put forward andevolved by various authors, e.g. advection-dominated ac-cretion flows (ADAF) [11, 12], standard and normal evo-lution (SANE) flows with poloidal magnetic field [13, 14]– for a comparison between codes for SANE evolutionssee [15] –, and magnetically arrested dominated (MAD)flows [16–19]. Self-consistent solutions for magnetizedthick discs with toroidal distributions of the magneticfield were obtained by Komissarov [20] for constant angu-lar momentum discs. This solution was extended to thenon-constant angular momentum case by Montero et al. [21] (see also [22]) and by Pimentel et al. [23, 24] who in-corporated magnetic polarization. Recently, equilibriumsolutions of self-gravitating magnetized discs in generalrelativity have been reported by Mach et al. [25], buildingon a procedure introduced by Shibata [26] for unmagne-tized tori (see also Stergioulas [27]). Numerical evolutionof those solutions, mainly in the purely hydrodynami-cal case, have been used to study the development ofpossible dynamical instabilities in the discs. These stud-ies include the runaway instability [28], the Papaloizou-Pringle instability [29]), the magneto-rotational instabil-ity [30] as well as the formation of jets and outflows [seee.g. 17, 19, 31–47]). a r X i v : . [ a s t r o - ph . H E ] F e b TABLE I. Summary of space-time and disc properties- in geometrized units - for our three values of the length scale λ andsome representative values of δ . From left to right the columns report the radii of the event horizon r EH , of the marginallybound orbit r mb , of the cusp of the gravitational potential r cusp , of the inner edge of the disc r in , of the center of the disc r c ,and of the outer edge of the disc r out , as well as the specific angular momentum l mb at r mb , the gap of the potential ∆ W , theangular velocity Ω, and the orbital period t orb at the center of the disc. δ r EH r mb r cusp r in r c r out l mb ∆ W Ω t orb λ = 10 +0 .
995 1.83 3.98 3.97 4.70 9.63 120.97 3.39 0.022 0.0313 32.0+0 .
500 1.97 3.98 3.98 4.74 10.32 102.83 3.88 0.039 0.0306 33.6+0 .
100 1.89 3.96 3.96 4.67 9.92 99.30 3.58 0.028 0.0298 32.7+0 .
000 2.00 4.00 4.00 4.67 10.47 108.52 4.00 0.043 0.0295 33.9 − .
100 2.04 4.03 4.03 4.80 10.63 111.28 4.15 0.048 0.0293 34.1 − .
300 2.17 4.16 4.15 4.95 11.08 116.95 4.61 0.062 0.0288 34.7 − .
500 2.43 4.50 4.51 5.35 11.89 107.87 5.57 0.081 0.0282 35.5 − .
800 5.28 9.91 9.90 11.59 24.03 269.63 15.50 0.099 0.0165 60.7 − .
995 400.0 800.0 800.0 952.0 2094.4 21704.8 800.0 0.043 0.0001 6777.7 λ = 60 +0 .
995 1.97 4.00 4.01 4.67 10.04 69.45 3.87 0.036 0.0314 31.9+0 .
500 1.98 4.00 4.01 4.67 10.16 76.27 3.91 0.038 0.0308 32.4+0 .
100 1.99 4.00 3.99 4.72 10.43 111.80 3.98 0.043 0.0297 33.7+0 .
000 2.00 4.00 4.00 4.67 10.47 108.52 4.00 0.043 0.0295 33.9 − .
100 2.01 4.00 3.99 4.67 10.57 120.04 4.03 0.045 0.0291 34.3 − .
300 2.03 4.01 4.01 4.77 10.78 130.68 4.11 0.048 0.0283 35.3 − .
500 2.07 4.02 4.02 4.78 11.26 182.45 4.27 0.057 0.0267 37.5 − .
800 2.30 4.22 4.21 5.06 13.28 262.41 5.20 0.099 0.0215 46.4 − .
995 399.0 800.0 799.9 952.0 2094.4 21708.8 800.0 0.043 0.0001 6777.7 λ = 1000 +0 .
995 2.00 4.00 4.00 4.76 10.42 98.48 3.99 0.042 0.0297 33.6+0 .
500 2.00 4.00 4.01 4.76 10.41 94.07 3.99 0.042 0.0298 33.6+0 .
100 2.00 4.00 3.99 4.76 10.48 111.68 4.00 0.043 0.0295 33.9+0 .
000 2.00 4.00 4.00 4.76 10.47 108.52 4.00 0.043 0.0295 33.9 − .
100 2.00 4.00 4.00 4.76 10.46 104.90 4.00 0.043 0.0296 33.8 − .
300 2.00 4.00 3.99 4.76 10.53 123.53 4.01 0.044 0.0293 34.2 − .
500 2.00 4.00 3.99 4.76 10.57 135.08 4.02 0.045 0.0291 34.4 − .
800 2.02 4.00 4.01 4.76 10.67 147.41 4.06 0.046 0.0287 34.8 − .
995 3.32 5.50 5.51 6.60 25.16 3671.75 9.80 0.277 0.0082 122.5All of those studies have been performed within theframework of general relativity in which astrophysicalblack holes are described by the Schwarzschild or Kerrsolutions. However, other types of black hole solutionshave been obtained in extended theories of gravity. Theobservational capabilities offered by the Event HorizonTelescope, targeted to measure black hole shadows andassociated strong-field lensing patterns from accretiondiscs around black holes, allows to test the validity of theblack hole solutions of general relativity. As a recent ex-ample Mizuno et al. [48] used the parameterized Einstein-Maxwell-dilaton-axion gravity solutions of Garc´ıa et al. [49], Konoplya et al. [50] to compare the shadows from aKerr black hole and a dilaton one, offering a proof of con-cept for the feasibility of such tests. More recently, theshadow of a boson star -surfaceless black hole mimicker-was also studied in Olivares et al. [51], where consideringrealistic astronomical observing conditions shows that ispossible to distinguish between Kerr black holes and non- rotating boson stars. Motivated by those works here weexplore the consequences of a simple extended model ofgravity that at “zeroth” perturbation order reproducesthe standard Schwarzschild space-time geometry of gen-eral relativity, but greatly differs from it when additionalterms are taken into account. This is done using a puremetric f ( R ) theory of gravity [52–56] by the introductionof a Yukawa-like potential following the proposal of DeMartino et al. [57], De Laurentis et al. [58]. As an as-trophysical application, in this article we explore two di-rections: (1) building stationary solutions of magnetizedthick discs with a self-consistent toroidal magnetic fieldaround a Yukawa-like black hole and (2) computing thephoton ring size in both general relativity and in ourextended theory of choice. The two solutions are com-pared with the aim of exploring whether the extended f ( R ) theory involving a Yukawa-like potential can stillbe valid with current observations. Previous attemptsto construct accretion discs around black holes in f ( R )theories of gravity can be found in P´erez et al. [59] andAlipour et al. [60] for thin and thick discs, respectively.In general relativity, the precise balance of the matter-energy content with the curvature reaction that this pro-duces on space-time is given by the well known Einsteinequations [see e.g 61]. When certain gravitational anoma-lies are observed and cannot be accounted for with a grav-itational theory, there are two paths to follow : either(a) the equations describing gravitation are modified andextended or (b) hypothetical dark matter and/or energycomponents are postulated. Despite the theoretical suc-cess of postulating such dark components the non-directdetection of dark matter and the less reliable constitu-tion of what dark energy might be, opens up the ideathat perhaps the Einstein equations might require exten-sions. For systems with mass-to-length ratios closely re-sembling that of the solar system, Newtonian gravity andgeneral relativity represent the correct theory of gravity[see e.g. 63]. When this ratio is smaller, as it occurs e.g.at the outskirts of many pressure-supported astronomi-cal systems like globular clusters, ellipical galaxies andclusters of galaxies, and for rotationally-supported sys-tems like spiral galaxies [see e.g. 64–66, and referencestherein] and in the universe as a whole at the presentepoch [cf. 67–69], an extension to the Einstein equationscan be postulated. At the other extreme, when gravi-tation is very strong, i.e. the mass-to-length ratio of asystem is much greater than that of the solar system, cur-rent gravitational-wave observations favour the Einsteinfield equations of general relativity [70, 71]. However, itcould still be possible that certain extended theories ofgravity may well describe current observations.This paper is organized as follows: In Section II wesummarize the problem setup, describing the black holespace-time in f ( R ) gravity and the procedure to buildthe disc solution. Stationary models of thick discs vary-ing the space-time parameters and the disc magnetizationare presented in Section III. This section also discussesthe dependence of the photon ring size on the space-timeparameters Finally Section IV summarizes our conclu-sions. Unless stated otherwise we use geometrized unitsin which the light speed, Newton’s constant, and themass of the black hole are equal to one, c = G = M = 1,the Kerr metric has the signature ( − , + , + , +), and the1 / π factor in the MHD equations is assumed to be one. These ideas were first conceived by Bouvard [62] after noticinganomalies in the orbital motion of the planet Uranus about thesun. Since then, every time an anomaly is seen in an astrophys-ical environment for which gravity has to be accounted for, amodification or extension of the gravitational theory is postu-lated or the existence of a hypothetical matter and/or energyentity is assumed.
II. SETUP
The immediate generalization of the Einstein equationsis done by allowing the Ricci scalar R in the gravita-tional action to be a general analytical function f ( R )[see e.g. 52–56, 72], with the matter action written inits usual form [see e.g. 73, 74]. The field equations inthis pure metric construction are then obtained by thenull variations of the whole action, i.e. the sum of thegravitational and matter actions, with respect to thespace-time metric. The obtained field equations turn outto be fourth-order differential equations for the metricand therefore, finding solutions of a well-posed partic-ular problem constitute a much harder task. By con-struction, when f ( R ) = R , the Einstein field equationsare recovered and the differential field equations are ofsecond order in the metric. A. Spherically symmetric black hole space-time ina pure metric static f ( R ) model In this article we consider a static spherically-symmetric space-time. The field equations are obtainedby the a specific choice of an f ( R ) function. In order toprovide a general scenario, Capozziello et al. [75] showedthat it was possible to find a weak-field limit solution thatcan be satisfied for all analytic f ( R ) functions (see alsoCapozziello and de Laurentis [54] and references therein).The idea is to expand in a Taylor series the function f ( R )and keep terms up to order 1 /c in the field equations.The resulting field equations at that order of approxima-tion have a Yukawa-like potential solution and, as shownby De Laurentis et al. [58] at that perturbation order,the general solution can be written as [see also 57, 76]: ds = − [1 + Φ( r )] dt + [1 − Φ( r )] dr + r d ˜Ω , (1)with Φ( r ) = − M (cid:0) δe − rλ + 1 (cid:1) r ( δ + 1) . (2)This Yukawa-like black hole (YBH) solution consti-tutes a generalization of the Schwarzschild space-timefor all analytic f ( R ) functions in the post-Newtonianlimit [see e.g. 63]. In the previous two equations d ˜Ω ≡ dθ + sin θdφ is the angular displacement and Φ( r ) isa Yukawa-like potential. Mathematically, the constantparameters λ and δ are related to the coefficients ofthe Taylor expansion of the function f ( R ) about a fixed R . In fact λ := (cid:112) − f (cid:48)(cid:48) /f (cid:48) and δ := f (cid:48) −
1, where[ ] (cid:48) := d[ ] / d R . As such, when f ( R ) = R , the grav-itational action becomes the Hilbert action of generalrelativity and so δ = 0 which leads to the Newtonianpotential φ ( r ) = − M/r for a point mass source and equa-tion (1) converges to the Schwarzschild exterior solution.
TABLE II. Maximum values of the rest-mass density ρ max and of the baryon mass of the disc M disc as a function of themagnetization of the plasma β . The quantities are shown for all values of the YBH space-time parameters δ and λ . We assumethat the baryon mass is M disc = 0 . M in the purely hydrodynamical case ( β = 10 ). δ β = 10 β = 10 β = 10 β = 10 β = 10 − β = 10 − β = 10 − ρ max λ = 10 +0 .
995 3 . × − . × − . × − . × − . × − . × − . × − +0 .
500 2 . × − . × − . × − . × − . × − . × − . × − +0 .
100 2 . × − . × − . × − . × − . × − . × − . × − +0 .
000 2 . × − . × − . × − . × − . × − . × − . × − − .
100 1 . × − . × − . × − . × − . × − . × − . × − − .
500 1 . × − . × − . × − . × − . × − . × − . × − − .
995 2 . × − . × − . × − . × − . × − . × − . × − λ = 60 +0 .
995 4 . × − . × − . × − . × − . × − . × − . × − +0 .
500 3 . × − . × − . × − . × − . × − . × − . × − +0 .
100 2 . × − . × − . × − . × − . × − . × − . × − +0 .
000 2 . × − . × − . × − . × − . × − . × − . × − − .
100 1 . × − . × − . × − . × − . × − . × − . × − − .
500 6 . × − . × − . × − . × − . × − . × − . × − − .
995 2 . × − . × − . × − . × − . × − . × − . × − λ = 1000 +0 .
995 2 . × − . × − . × − . × − . × − . × − . × − +0 .
500 2 . × − . × − . × − . × − . × − . × − . × − +0 .
100 1 . × − . × − . × − . × − . × − . × − . × − +0 .
000 2 . × − . × − . × − . × − . × − . × − . × − − .
100 2 . × − . × − . × − . × − . × − . × − . × − − .
500 1 . × − . × − . × − . × − . × − . × − . × − − .
995 5 . × − . × − . × − . × − . × − . × − . × − M disc = (cid:82) √ γW ρd xλ = 10 +0 .
995 1 . × − . × − . × − . × − . × − . × − . × − +0 .
500 1 . × − . × − . × − . × − . × − . × − . × − +0 .
100 1 . × − . × − . × − . × − . × − . × − . × − +0 .
000 1 . × − . × − . × − . × − . × − . × − . × − − .
100 1 . × − . × − . × − . × − . × − . × − . × − − .
500 1 . × − . × − . × − . × − . × − . × − . × − − .
995 1 . × − . × − . × − . × − . × − . × − . × − λ = 60 +0 .
995 1 . × − . × − . × − . × − . × − . × − . × − +0 .
500 1 . × − . × − . × − . × − . × − . × − . × − +0 .
100 1 . × − . × − . × − . × − . × − . × − . × − +0 .
000 1 . × − . × − . × − . × − . × − . × − . × − − .
100 1 . × − . × − . × − . × − . × − . × − . × − − .
500 1 . × − . × − . × − . × − . × − . × − . × − − .
995 1 . × − . × − . × − . × − . × − . × − . × − λ = 1000 +0 .
995 1 . × − . × − . × − . × − . × − . × − . × − +0 .
500 1 . × − . × − . × − . × − . × − . × − . × − +0 .
100 1 . × − . × − . × − . × − . × − . × − . × − +0 .
000 1 . × − . × − . × − . × − . × − . × − . × − − .
100 1 . × − . × − . × − . × − . × − . × − . × − − .
500 1 . × − . × − . × − . × − . × − . × − . × − − .
995 1 . × − . × − . × − . × − . × − . × − . × − The parameter λ is a length scale which can in prin-ciple be adjusted depending on the spatial scale of theparticular astrophysical system (see below). Moreover, δ is the parameter of the theory and governs the strength of the Yukawa-like potential (general relativity and there-fore Newton’s potential is recovered when δ = 0). Theevent horizon of the YBH is computed in the same wayas for the Schwarzschild black hole, solving the condi-tion g tt (r) = 0, where the surface of infinite redshift andthe event horizon coincide. Since Yukawa’s potential hasa nonlinear dependence on the radial coordinate we ob-tain a transcendental equation for r EH . Hence, we use aNewton-Raphson root-finder to compute the event hori-zon.For our YBH solution the angular velocity of Kepleriancircular orbits around the black hole readsΩ ( r ) = M δe − r/λ r λ ( δ + 1) − Φ2 r . (3)For circular orbits, we can write the angular velocity andthe specific angular momentum in terms of the nonzerocomponents of the 4-velocity u µ , namely Ω = u φ /u t and l = − u φ /u t . Using this the expression for the Keplerianspecific angular momentum can be written as l K = ± r (cid:115) M δe − r/λ r λ ( δ + 1) − Φ2 r . (4)We also define the specific bound angular momentumfunction l b ( r ) which corresponds to the specific angu-lar momentum of a marginally bound orbit at a certainradius r and can be written in our case as l ( r ) = − r Φ1 + Φ . (5)If we consider prograde (retrograde) motion, finding theminimum (maximum) of the function l b ( r ) gives the loca-tion of the innermost marginally bound circular orbit r mb and the value of the specific angular momentum there( l b ( r mb ) = l mb ). It is also worth noticing that at saidpoint, the Keplerian angular momentum is equal to l mb as well (i.e. l K ( r mb ) = l mb ). B. Procedure to build equilibrium magnetizedthick discs
We build sequences of equilibrium thick discs endowedwith a toroidal magnetic field in YBH space-time follow-ing the procedure first presented in [20] and generalizedby [21, 22]. For simplicity we assume that the plasma inthe disc obeys a constant distribution of specific angularmomentum l = l K ( r mb ) given by the Keplerian angularmomentum equation (4), evaluated at r mb . The funda-mental equation to describe a non-self-gravitating equi-librium torus around a black hole is obtained by applyingthe projection tensor h αβ = δ αβ + u α u β to the conservationlaw of the energy-momentum tensor [22]. This equationreads ∂ i (ln | u t | ) − Ω ∂ i l − l Ω + ∂ i pρh + ∂ i (cid:2) L b (cid:3) L ρh = 0 , (6)where i = r, θ . To obtain the previous equation we haveassumed that the thermodynamical relationship between the rest-mass density ρ and the thermal pressure p isgiven by a barotropic equation of state (EoS), ρ = ρ ( p ).In particular, we choose a polytropic EoS such as p = Kρ Γ and an EoS for the magnetic pressure p m ≡ b / p m = K m L q − ( ρh ) q , where K , K m , q and Γ areconstants and L ≡ g tφ − g tt g φφ . Moreover, h and b inEq. (6) are the enthalpy and the modulus (squared) ofthe magnetic field 4-vector. Using this relations we canrewrite Eq. (6) as W − W in + ln (cid:18) K ΓΓ − ρ Γ − (cid:19) + qK m q − (cid:20) L (cid:18) ρ + K Γ ρ Γ Γ − (cid:19)(cid:21) q − = 0 , (7)where W = ln | u t | is the (gravitational plus centrifugal)potential. To solve equation (7) we fix q = Γ = 4 /
3, thedensity at the center of the disc ρ c = 1, the specific an-gular momentum l = l mb and we fill 80% of the potentialgap ∆ W ≡ W in − W cusp , where subindices ‘in’ and ‘cusp’indicate that the potential is calculated at the inner edgeof the disc or at the cusp (see below). Therefore, in ourmodels the discs will always be inside their correspondingRoche lobes (∆ W < r, θ ) grid in a domain r ∈ [ r EH , r out ], whose specificvalues are reported in Table I. The number of zones inour base grid is 252 ×
256 in r and θ , respectively. III. RESULTSA. YBH parameters
Before constructing the YBH-disc solutions we ex-plore suitable values of the freely specifiable parame-ters of the theory, λ and δ . This serves the purposeof understanding the intrinsic properties of the space-time and offers the possibility of comparing our findingswith the analysis of [57, 58]. We consider three lengthscales, λ = 10 , , and 1000 in geometrized units. Inphysical units and for the case of M87 the first casecorresponds to the scale of the black hole photon ringshadow, λ phys (10) ∼ . × − pc (23 µ as), the secondone to the size of the inner core of the jet, λ phys (60) ∼ . × − pc (230 µ as), and the third case to the largescale jet of M87, λ phys (10 ) ∼ . × − pc (3 . λ phys (10) ∼ . × − pc (50 µ as), λ phys (60) ∼ . × − pc (300 µ as), and λ phys (10 ) ∼ . × − pc (5mas) . For each value of λ we use 14 To estimate the length scales in µ as we assume the followingblack hole masses and distances to the source: M M87 = (6 . ± . × M (cid:12) and D = 16 . ± . SgrA ∗ =(4 . ± . × M (cid:12) and D = 8 .
178 Mpc for the galacticcenter SgrA* [77]. values of δ . We focus our attention in the case δ < δ and λ and the strength of the toroidal magnetic field at thecenter of the tori, β c .Note that the choices δ < M < r ) < B. Tori geometry
Table I reports the values of selected geometrical quan-tities of the discs for a subset of representative models(parameterized by δ ) for all three values of λ . We notethat those quantities are only related to the space-timeand hence do not depend on the magnetization param-eter of the discs. Varying parameter δ ∈ [ − ,
1] yieldsYBHs with different event horizon radii r EH . We find thatthe event horizon size increases for negative values of δ and small values of λ reaching r EH = 400 M for λ = 10and 60. On the other hand, for λ = 1000 the variationof r EH with negative δ is not too pronounced, stayingat r EH = 2 for most models and reaching r EH = 3 . δ = − . δ values, the values of r EH we obtain are comparable to the event horizon size ofa slowly rotating black hole with r EH = 1 .
83 and 1 . λ = 10 and 60) corresponding to Kerr black holeswith spin parameters a = 0 . , and 0 . δ and λ = 1000 no changes areobserved with respect to the event horizon radius of anon-rotating black hole in general relativity.Table I also reports the radius of the marginally boundorbit r mb and its corresponding specific angular momen-tum l mb . For all λ , these two quantities show a weakdependence on δ except for extreme values very close to δ = −
1. Additional disc radii reported in Table I are thelocation of the cusp, r cusp , the center of the disc, r c , andthe inner and outer edges of the disc, r in , r out , respec-tively. The latter are computed assuming the discs fill80% of the gap of the potential ∆ W (i.e. of their Rochelobes). The center of the disc is defined as the minimumof W . The orbital period of the disc, t orb , reported in thelast column of Table I is measured at the center of thedisc. It is found that all characteristic quantities defin-ing the torus size, r cusp , r c , r in , and r out , only show astrong dependence on the YBH parameters as δ → − δ goes from positive to nega-tive values, except for δ = − .
995 where the incrementis significantly larger.The gap of the potential, ∆ W , also reported in TableI, defines the regions where the equilibrium plasma is lo- cated around the black hole. In particular, the surface ofthe magnetized disc is defined as the equipotential surfacewith W = W in and the solution of the thermodynamicalquantities (see Eq. (7)) also depends on this gap. In gen-eral relativity ( δ = 0) this value is ∆ W = 0 . λ . The largest deviations found are ∆ W = 0 . δ = − . , .
176 (for δ = − . , and 0 .
277 (for δ = − . λ = 60 and for different valuesof δ and considering a magnetization parameter at thecenter of the disc of β c = 10 − . Therefore, the examplesshown in Fig. 1 correspond to highly magnetized tori. Forclarity in the comparison we rescale the spatial domain( r cos θ, r sin θ ) by the event horizon size of each YBH(the specific values are reported in Table I) showing adomain [ − r EH , r EH ] in the x − z plane.As δ increases from 0 towards δ = 1, the size of the discdecreases while r EH is kept essentially constant, slightlychanging from 2.0 to 1.97 (see Table I). The most no-ticeable modifications are visible, however, only for thelargest values. For δ = 0 .
995 the disc is about 40%smaller than in general relativity. On the other hand, as δ decreases from 0 towards δ = −
1, the size of the discsbecomes gradually larger. The largest value in Fig. 1 cor-responds to the δ = − .
995 case. Note that the appar-ent smaller size of this model as compared e.g. with the δ = − . r EH which is 2.3 for the latterand about 400 for the former. Despite these significantmodifications in the geometrical size, the distribution ofthe density and of the magnetization in the torus seemonly weakly affected by the changes in δ . For the mostnegative values of δ further extended low-density layersare obtained as well as high-density regions along thesymmetry axis of the black hole. Similarly, the angularthickness of the disc also increases.Fig. 1 shows that the inner edge of each torus is at thesame distance of the YBH horizon since we are using asunit of distance the event horizon size. The same occursfor the location of the cusp of the potential, as confirmedby the values reported in Table I. Nevertheless, we noticedifferences in the location of the center of the discs r c and of their outer edge r out . The top plot in Figure 2displays the position of the center of the tori in unitsof the event horizon radii of the corresponding YBH asfunction of δ and for all three values of λ . Since the centerof the disc is defined as the location of the minimum ofthe potential, it does not depend on the magnetic fieldstrength. We observe noticeable changes in the positionof the disc center for negative values of δ . For the case λ = 60, r c increases monotonically with r EH as δ → − δ = − . r in ≈ . r EH . For δ = − . r c decreases to about 5 . r EH . Note, however, that in terms z / r E H x / r EH z / r E H x / r EH x / r EH x / r EH = 0.995 = 0.10 = 0.00 = 0.10= 0.30 = 0.50 = 0.80 = 0.9956 5 4 3 2 1 0 1 log ( ) log ( ) FIG. 1. Logarithm of the rest-mass density (left half portion of each panel) and magnetization parameter (right half) for anillustrative sample of tori around YBHs with λ = 60 and different values of δ . All tori are built assuming high magnetization,using a central magnetization parameter of β c = 10 − . Each box covers a spatial domain [ − r EH , r EH ] in the x − z plane.Note that each solution has a different length scale if expressed in conventional black hole mass M units, due to the dependenceof r EH with δ (see Table I for details). of the YBH mass, the radial position of the center of thedisc always increases monotonically (see Table I).The bottom plot of Figure 2 shows the correspondingangular thickness - expressed in radians - between thesurface of the tori, defined by the equipotential surface,and the x -axis as a function of δ . In the Schwarzschildcase, r c ∼ . r EH and the angular thickness is θ disc ∼ π/
5. Therefore, for negative values of δ more elon-gated discs with higher angular thickness are obtained,as shown in the bottom plot of Fig. 2.The dependence we have just described is affected bythe value of the length scale parameter λ . The corre-sponding results for λ = 10 and 1000 are also plottedin Fig. 2. For λ = 10, r c decreases monotonically with r EH as δ → −
1, up to δ ≈ − .
8. For even more neg-ative values of δ , r c increases. The discs are in generalsmaller than in the λ = 60 case and their angular sizes,which depend weakly with δ are also smaller. Finally,for λ = 1000 both r c and θ disc increase monotonicallywith r EH as δ → −
1. The largest discs with the highestangular thicknesss are found for this value of λ .The range of variation of r out with δ in our modelsis also significant, as reported in Table I. For λ = 10, r out ∼ − r EH , for λ = 60, r out ∼ − r EH and for λ = 1000, r out ∼ − r EH . This effect is aconsequence of the particular value of the gravitational- centrifugal potential gap for each model (which is a non-linear function of the Yukawa-like potential), which in-creases as we move from positive to negative values of δ , modifying the Roche lobes of each YBH solutions andthe morphology and thermodynamics of the equilibriumtorus. Such behavior depends directly on the Yukawa-likepotential; for small λ ’s the exponential function in equa-tion (1) has an important contribution on the space-time.On the other hand, for large values of λ the exponen-tial part goes to zero and the potential mostly dependson 1 /δ , consistent with the geometry and the thermody-namics of the magnetized discs (see Figs. 2, 4, 5). C. Tori thermodynamics
The differences found in the geometry of the discsare accompanied by quantitative differences in the phys-ical magnitudes characterizing the matter content of ourmodels and their thermodynamics. Figure 3 shows ra-dial profiles along the equatorial plane ( θ = π/
2) of therest-mass density in logarithmic scale for two values ofthe magnetization parameter and for all values of pa-rameters λ and δ . The profiles are shown for both anunmagnetized disc ( β = 10 ) and a highly magnetizedone ( β = 10 − ). For YBH space-time, the maximum r c [ r E H ] = 10 = 60 = 1000 d i c s [] = 10 = 60 = 1000 FIG. 2. Dependence on δ of the position of the center of thedisc r c in units of r EH (top) and of the angular thickness ofthe disc measured with respect to the x -axis at r c (bottom).Note that the dependence of both quantities on δ does notchange with the magnetic field strength. of the rest-mass density increases for high magnetizeddiscs and its location shifts towards the inner edge ofthe torus. This is in agreement with what is found forthe case of Schwarzschild space-time [79]. However, thedeviations are not large and they only become signifi-cant for δ = − . λ = 10 , δ = − .
995 and β = 10 ).Table II reports, for all of our models, the values ofthe maximum of the rest-mass density and of the baryonmass of the disc, defined as M disc = (cid:90) √ γW ρd x , (8)where W is the Lorentz factor. This table allows to quan-tify the effects on these two quantities of the parameters δ and λ that characterize the YBH spectime and to findthe dependence on the magnetic-field strength, from un-magnetized to highly magnetized discs. We assume that the mass of the disc is 10% of the mass of the black holefor the unmagnetized case, M HD = 0 . M , fixing the rest-mass density, ρ HD . The highest value of M disc is attainedfor the unmagnetized model ( β = 10 ) and the mass ofthe tori decreases monotonically as the magnetization in-creases.Figures 4 and 5 show the maximum of the rest-massdensity and of the baryon mass of the disc as a function of δ , normalized by ρ HD and M HD , respectively. Both quan-tities do not show a strong dependence on δ , irrespectiveof the value of λ , except for values of δ close to -1, wherethe largest deviations are found. The increase in ρ max ismonotonic with the increase of the disc magnetization,and it is similar for all three values of λ . Except for δ → − . ρ HD which is the value found for toriaround a Schwarzschild black hole by [79]. Only when δ → − ∼ . , . , . ρ HD for λ = 10 , , ∼ M HD / , M HD /
10 and M HD / D. Constraining the YBH parameters with thephoton ring size
The Event Horizon Telescope observations of the blackhole shadow in M87 [7, 81] and the forthcoming ob-servations of SgrA* provide a laboratory to test gen-eral relativity and modified theories of gravity by usingthe shadow properties. We turn now to analyze howthe shadow size of our YBH solution differs from thatof a Schwarzschild black hole, constraining the parame-ters of the space-time with the size of the photon ring.To this aim we solve the geodesic equations for pho-ton (null) trajectories in the YBH space-time given byEq. (1) following the procedure developed by Hioki andMaeda [82], Hou et al. [83]. We neglect any contributionfrom the emissivities and absorptivities of the magnetizedtorus (i.e. we do not consider the spectral content of thelight) because for our stationary models we assume thatthe emission from the photon ring is brighter than in thedisc, but the disk is sufficiently illuminated so that a pho-ton ring is visible. A detailed assessment of the shadowsbased on numerical evolutions of our YBH-torus systemswill be presented elsewhere. The apparent shape of theYBH shadow is computed by introducing the celestialcoordinates ( α, β ) assuming a source located at infinityand a viewing angle i o = π/
2. Such coordinates measureapparent angular distances of the image on the celestialsphere. The shadow of the black hole is defined by abright ring at the radius of the lensed photon sphere orphoton ring [84]. For a Kerr black hole at different view-ing angles such radius varies between r sh ∼ √ g ± √ r g is the photon ring of the Schwarzschildblack hole, and r g is the gravitational radius. The ± = 10= 10 = 60= 10 = 10 = 10 r [ r c ] = 10= 10 r [ r c ] = 60= 10 r [ r c ] = 10 = 10 = + 0.995 = + 0.500 = + 0.100 = + 0.000 = 0.100 = 0.500 = 0.995 FIG. 3. Radial profiles of the rest-mass density of the magnetized discs in logarithmic scale at the equatorial plane. The leftpanels correspond to the YBH space-time with λ = 10, the middle panels to λ = 60 and the right panels to λ = 10 . The toppanels show low magnetized discs and the bottom ones display highly magnetized cases. In each plot seven values for δ areshown. The radial coordinate is in units of the center of the disc to facilitate the comparison. m a x / H D = 10 = 10 = 10 = 10 = 10 = 10 FIG. 4. Maximum of the rest-mass density normalized by the corresponding value for a non-magnetized disc ρ HD (see Table I),for λ = 10 (left), λ = 60 (middle), and λ = 10 (right). The maximum density increases with the magnetization, in agreementwith the Schwarzschild black hole case. The largest increments are found for δ < black hole shadow correspond to the maximally rotat-ing Kerr black hole cases (with Kerr dimensionless spinparameter a = ± o = π/ λ = 10 ,
60, and 10 . In the middle panels we show thecomplete shape of the photon rings in celestial coordi- nates, normalized by the event horizon of the black hole,for representative values of the δ parameter. In the toppanels we show a close-up of the upper region of thering while the bottom panels display the radius of thephoton ring r sh as a function of δ . The orange regionsin these plots are the shadow sizes estimated from gen-eral relativity for Kerr black holes. Almost all of ourmodels are consistent with Event Horizon Telescope ob-0 M d i s c / M H D = 10 = 10 = 10 = 10 = 10 = 10 = 10 FIG. 5. Disc mass as a function of δ normalized by the mass of an unmagnetized disc, M HD (see Table II), for differentmagnetizations. From left to right the panels correspond to λ = 10, λ = 60, and λ = 10 , respectively. The disc mass decreaseswhen the magnitude of the magnetic field increases. The least massive torus is obtained for λ = 10 and δ → − / r E H = 10 / r E H = 60 = +0.995= +0.500= +0.100 = +0.000= 0.100= 0.300 = 0.500= 0.800= 0.995 / r E H = 1000 / r EH / r E H / r EH / r E H / r EH / r E H r s h / r E H r s h / r E H r s h / r E H FIG. 6. Photon ring size (shadow size) of different YBH space-times with varying δ parameter and length scales λ . The fullrings are displayed in celestial coordinates ( α, β ) normalized by the event horizon size in the middle panels, while a close-up isshown in the top panels. The shadow radius as a function of δ is shown in the bottom panels (blue curves), where the orangeregions indicate the range of the shadow size for maximally rotating Kerr black holes. The dashed gray lines correspond to theSchwarzschild black hole shadow radius, r sh = 3 √ g [80] and the black line is the radius of the event horizon. An observerview angle i o = π/ servations [84]. The largest deviations are obtained forthe most negative values of δ , irrespective of λ , albeitfor λ = 60 the variations are the smallest. On the otherhand, for positive values of δ the photon ring size onlyincreases slowly with δ for the two lowest values of λ we consider.To obtain a better understanding of the dependence ofthe shadow size on δ and λ we finish our analysis by ex-ploring 10 space-time models varying both parameters.The results are displayed in Figure 7 which shows the1shadow radius normalized by the event horizon radius asa function of δ and λ . The blue isocontour correspondsto a family of YBH with the same photon ring size asthe Schwarzschild black hole. The red iscontours are thebounds for 3 √ r g ± δ > λ <
50 generate small photon rings, whereas largephoton rings are produced for δ < λ .Irrespective of the sign of δ we observe nonlinear cor-relations between the shadow size and the δ parameter.In particular, at small scales, for λ <
20, the effects ofthe Yukawa-like gravitational potential are strong, lead-ing to a high variability of the photon ring size - fromminimum to maximum - when the δ parameter decreasesfrom δ ∼ . δ ∼ − .
25. This is expected since lightbending is more noticeable in stronger gravity regimes.Correspondingly, for asymptotic values λ → ∞ the YBHphoton rings tend asymptotically to the Schwarzschildphoton ring.Applying the δ parameter constriction to the super-massive black holes of M87 and the galactic center ( λ ∼
60) we find that astrophysically accepted values for δ would fall in the range − . < δ < .
0. For our sta-tionary accretion disc solutions this implies that we canneglect very thick torus with large angular thicknesss, θ disc > . π , for which the location of the center isr c > EH and have large event horizon, r EH ∼ M .The maximum densities in YBH-torus systems compati-ble with the constrained range of δ are in agreement withthe values found in general relativity [79] with small de-viations ∼ ± . IV. SUMMARY
We have presented stationary solutions of geometri-cally thick discs (or tori) with constant angular momen-tum and endowed with a self-consistent toroidal mag-netic field distribution surrounding a non-rotating blackhole in f ( R )-gravity. The particular f ( R )-gravity modelwe have employed introduces a Yukawa-like modificationto the Newtonian potential, encoded in a single param-eter δ and whose specific values affect the disc config-urations compared to the general relativistic case. Wehave built models for different magnetic field strengths,from low magnetized discs (essentially hydrodynamic) tohighly magnetized tori. This has been achieved by ad-justing the magnetization parameter β , i.e., the ratioof thermal pressure to magnetic pressure, in the rangelog β ∈ {− , } . Our stationary solutions have beenobtained numerically, employing the approach discussedin detail in [79].The characteristics of our solutions have been quanti-fied by analyzing the central density of the discs, theirbaryonic mass, their geometrical size and angular thick- Shadow Radius [r EH ]1.0 0.5 0.0 0.5 1.010 FIG. 7. Shadow radius of the YBH as function of δ and λ , nor-malized by the event horizon. The blue contour correspondsto shadow size for Schwarzschild black hole r sh ∼ √
3, in r g radius units and red contours to r sh ± ness, as well as the effects of the deviations of the YBHmetric from the Schwarzschild metric. We have foundthat in the general relativistic limit ( δ = 0) our modelsreproduce our previous results for a Schwarzschild blackhole [22, 79]. For small values of the δ parameter, cor-responding to ∼
10% deviations from general relativity,we have found small geometrical variations in the mod-els with respect to the general relativistic results, namely ∼
2% in the event horizon size, a ∼
5% shift in the loca-tion of the inner edge and center of the disc, and ∼ | δ | > . δ →−
1. Those modifications of the gravitational potentialhave a large direct impact in the torus solution. We havefound that the influence of the magnetic field in the discproperties becomes stronger in this case. In particular wehave observed an increment of the YBH event horizon ofabout four orders of magnitude with respect to the eventhorizon of a Schwarzschild black hole, a ∼
10% increaseof r c ( r EH ) and three orders of magnitude increase of thelocations of the outer edge of the disc, r out .The impact of the strength of the toroidal magneticfield in the morphology of the discs follows the sametrend of previous analysis for non-rotating black holesin general relativity [22, 79]. The maximum density forour most highly magnetized disc is ∼ . ρ HD and thetotal mass of the disc is ∼ . M HD . The variationin these two quantities is less than 10% for small devia-2tions from general relativity. For | δ | > . δ influence the angular size of the discs, whichbecome more elongated along the z -axis. This may af-fect the angular size of outflows and jets that might formwhen evolving these magnetized discs in the alternativetheory of gravity discussed here.Finally, we have analysed the differences between theYBH space-time and the Schwarzschild space-time bycomputing the size of the photon ring produced by asource located at infinity. This has allowed us to con-strain the parameters δ and λ of the YBH space-timewhen applying our approach to the supermassive blackholes of M87 and SgrA ∗ . The variations in the photonring reported in this work, even for small deviations fromgeneral relativity, might be measurable in upcoming ob-servations of the galactic center by the Event HorizonTelescope Collaboration.Even though the variations we have found in our YBH-disc models are less than 10% for realistic deviations ofgeneral relativity, namely δ = ± .
1, non-linear time evo-lution of these initial data might still disclose potentialdifferences in observable quantities like jet power or massaccretion rates, as well as in the time-dependent morphol- ogy of the photon ring produced by a turbulent, illumi-nating disc. Those results will be presented elsewhere.
ACKNOWLEDGEMENTS
We thank H´ector Olivares, Oliver Porth, and ZiriYounsi for numerous helpful discussions and comments.ACO gratefully acknowledges support from the COSTAction CA16214 “PHAROS”, the LOEWE-Programin HIC for FAIR, and the EU Horizon 2020 ResearchERC Synergy Grant “Black-HoleCam: Imaging theEvent Horizon of Black Holes” (grant No. 610058).JAF is supported by the Spanish Agencia Estatal deInvestigaci´on (PGC2018-095984-B-I00) and by the Gen-eralitat Valenciana (PROMETEO/2019/071). M.D.L.acknowledges INFN sez. di Napoli, Iniziative SpecificheQGSKY and TEONGRAV. SM acknowledges supportfrom DGAPA-UNAM (IN112019), Consejo Nacional deCiencia y Teconolog´ıa (CONACyT), M´exico (CB-2014-01 No. 240512) and CONACyT 26344. This work hasfurther been supported by the European Union’s Hori-zon 2020 Research and Innovation (RISE) programmeH2020-MSCA-RISE-2017 Grant No. FunFiCO-777740.The simulations were performed on the LOEWE clusterin CSC in Frankfurt. [1] S. E. Woosley, Astrophys. J. , 273 (1993).[2] L. Baiotti and L. Rezzolla, Reports on Progress inPhysics , 096901 (2017), arXiv:1607.03540 [gr-qc].[3] B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese,K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Ad-hikari, V. B. Adya, and et al., Physical Review Letters , 161101 (2017), arXiv:1710.05832 [gr-qc].[4] B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese,K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Ad-hikari, V. B. Adya, and et al., Astrophysical JournalLetters , L12 (2017), arXiv:1710.05833 [astro-ph.HE].[5] N. I. Shakura and R. A. Sunyaev, Astron. Astrophys. , 33 (1973).[6] M. J. Rees, Annual Review of Astronomy and Astro-physics , 471 (1984).[7] Event Horizon Telescope Collaboration, K. Akiyama,A. Alberdi, W. Alef, K. Asada, R. Azulay, A.-K. Baczko,D. Ball, M. Balokovi´c, J. Barrett, et al. , Astrophys. J.Lett. , L1 (2019).[8] L. G. Fishbone and V. Moncrief, Astrophys. J. , 962(1976).[9] M. Kozlowski, M. Jaroszynski, and M. A. Abramowicz,Astron. Astrophys. , 209 (1978).[10] M. A. Abramowicz and P. C. Fragile, Living Reviews inRelativity , 1 (2013), arXiv:1104.5499 [astro-ph.HE].[11] R. Narayan and I. Yi, Astrophys. J. Lett. , L13(1994), arXiv:astro-ph/9403052.[12] F. Yuan and R. Narayan, Annual Review of Astronomyand Astrophysics , 529 (2014), arXiv:1401.0586 [astro-ph.HE]. [13] R. Narayan, A. S¸adowski, R. F. Penna, and A. K.Kulkarni, Mon. Not. R. Astron. Soc. , 3241 (2012),arXiv:1206.1213 [astro-ph.HE].[14] A. S¸adowski, R. Narayan, A. Tchekhovskoy, andY. Zhu, Mon. Not. R. Astron. Soc. , 3533 (2013),arXiv:1212.5050 [astro-ph.HE].[15] O. Porth, K. Chatterjee, R. Narayan, C. F. Gammie,Y. Mizuno, P. Anninos, J. G. Baker, M. Bugli, C.-k.Chan, J. Davelaar, L. Del Zanna, Z. B. Etienne, P. C.Fragile, B. J. Kelly, M. Liska, S. Markoff, J. C. McK-inney, B. Mishra, S. C. Noble, H. Olivares, B. Prather,L. Rezzolla, B. R. Ryan, J. M. Stone, N. Tomei, C. J.White, Z. Younsi, and The Event Horizon TelescopeCollaboration, arXiv e-prints , arXiv:1904.04923 (2019),arXiv:1904.04923 [astro-ph.HE].[16] R. Narayan, I. V. Igumenshchev, and M. A. Abramow-icz, Publications of the ASJ , L69 (2003), astro-ph/0305029.[17] J.-P. De Villiers and J. F. Hawley, Astrophys. J. ,1060 (2003), arXiv:astro-ph/0303241 [astro-ph].[18] A. Tchekhovskoy, R. Narayan, and J. C. McKin-ney, Mon. Not. R. Astron. Soc. , L79 (2011),arXiv:1108.0412 [astro-ph.HE].[19] J. C. McKinney, A. Tchekhovskoy, and R. D. Bland-ford, Mon. Not. R. Astron. Soc. , 3083 (2012),arXiv:1201.4163 [astro-ph.HE].[20] S. S. Komissarov, Mon. Not. R. Astron. Soc. , 993(2006), astro-ph/0601678.[21] P. J. Montero, O. Zanotti, J. A. Font, and L. Rezzolla,Mon. Not. R. Astron. Soc. , 1101 (2007), arXiv:astro- ph/0702485.[22] S. Gimeno-Soler and J. A. Font, Astron. Astrophys. ,A68 (2017).[23] O. M. Pimentel, F. D. Lora-Clavijo, and G. A. Gonzalez,Astron. Astrophys. , A57 (2018), arXiv:1808.07400[astro-ph.HE].[24] O. M. Pimentel, F. D. Lora-Clavijo, and G. A. Gonz´alez,Astrophys. J. , 115 (2018), arXiv:1806.02266 [astro-ph.HE].[25] P. Mach, S. Gimeno-Soler, J. A. Font, A. Odrzywo(cid:32)lek,and M. Pir´og, Phys. Rev. D , 104063 (2019).[26] M. Shibata, Phys. Rev. D , 064035 (2007).[27] N. Stergioulas, in Journal of Physics Conference Series ,Journal of Physics Conference Series, Vol. 283 (2011) p.012036.[28] M. A. Abramowicz, M. Calvani, and L. Nobili, Nature , 597 (1983).[29] J. C. B. Papaloizou and J. E. Pringle, Mon. Not. R. As-tron. Soc. , 721 (1984).[30] S. A. Balbus and J. F. Hawley, Astrophys. J. , 214(1991).[31] J. F. Hawley, Astrophys. J. , 496 (1991).[32] J. A. Font and F. Daigne, Mon. Not. R. Astron. Soc. ,383 (2002), arXiv:astro-ph/0203403 [astro-ph].[33] C. F. Gammie, J. C. McKinney, and G. T´oth, Astrophys.J. , 444 (2003), arXiv:astro-ph/0301509 [astro-ph].[34] L. Rezzolla, O. Zanotti, and J. A. Font, Astron. Astro-phys. , 603 (2003), arXiv:gr-qc/0310045 [gr-qc].[35] F. Daigne and J. A. Font, Mon. Not. R. Astron. Soc. ,841 (2004), arXiv:astro-ph/0311618.[36] O. Zanotti, J. A. Font, L. Rezzolla, and P. J. Montero,Mon. Not. R. Astron. Soc. , 1371 (2005), arXiv:astro-ph/0411116 [astro-ph].[37] P. C. Fragile, O. M. Blaes, P. Anninos, andJ. D. Salmonson, Astrophys. J. , 417 (2007),arXiv:0706.4303 [astro-ph].[38] P. J. Montero, J. A. Font, and M. Shibata, Phys. Rev.Lett. , 191101 (2010), arXiv:1004.3102 [gr-qc].[39] K. Kiuchi, M. Shibata, P. J. Montero, and J. A. Font,Phys. Rev. Lett. , 251102 (2011), arXiv:1105.5035[astro-ph.HE].[40] O. Korobkin, E. B. Abdikamalov, E. Schnetter, N. Ster-gioulas, and B. Zink, Phys. Rev. D , 043007 (2011),arXiv:1011.3010 [astro-ph.HE].[41] O. Korobkin, E. Abdikamalov, N. Stergioulas, E. Schnet-ter, B. Zink, S. Rosswog, and C. D. Ott, Mon. Not.R. Astron. Soc. , 349 (2013), arXiv:1210.1214 [astro-ph.HE].[42] M. Wielgus, P. C. Fragile, Z. Wang, and J. Wilson, Mon.Not. R. Astron. Soc. , 3593 (2015).[43] V. Mewes, J. A. Font, F. Galeazzi, P. J. Montero,and N. Stergioulas, Phys. Rev. D , 064055 (2016),arXiv:1506.04056 [gr-qc].[44] P. C. Fragile and A. S¸adowski, Mon. Not. R. Astron. Soc. , 1838 (2017), arXiv:1701.01159 [astro-ph.HE].[45] M. Bugli, J. Guilet, E. M¨uller, L. Del Zanna, N. Buc-ciantini, and P. J. Montero, Mon. Not. R. Astron. Soc. , 108 (2018).[46] V. Witzany and P. Jefremov, Astron. Astrophys. ,A75 (2018), arXiv:1711.09241 [astro-ph.HE].[47] A. Cruz-Osorio, S. Gimeno-Soler, and J. A. Font, Mon.Not. R. Astron. Soc. , 5730 (2020), arXiv:2001.09669[gr-qc]. [48] Y. Mizuno, Z. Younsi, C. M. Fromm, O. Porth,M. De Laurentis, H. Olivares, H. Falcke, M. Kramer,and L. Rezzolla, Nature Astronomy , 585 (2018),arXiv:1804.05812 [astro-ph.GA].[49] A. Garc´ıa, D. Galtsov, and O. Kechkin, Phys. Rev. Lett. , 1276 (1995).[50] R. Konoplya, L. Rezzolla, and A. Zhidenko, Phys. Rev.D , 064015 (2016), arXiv:1602.02378 [gr-qc].[51] H. Olivares, Z. Younsi, C. M. Fromm, M. De Lauren-tis, O. Porth, Y. Mizuno, H. Falcke, M. Kramer, andL. Rezzolla, Mon. Not. R. Astron. Soc. , 521 (2020),arXiv:1809.08682 [gr-qc].[52] S. Capozziello and V. Faraoni, Beyond Einstein Gravity:A Survey of Gravitational Theories for Cosmology andAstrophysics , Fundamental Theories of Physics (SpringerNetherlands, 2010).[53] S. Capozziello and M. De Laurentis,
Invariance Princi-ples and Extended Gravity: Theory and Probes , PhysicsResearch and Technology (Nova Science Publishers,2010).[54] S. Capozziello and M. de Laurentis, Phys. Rep. , 167(2011), arXiv:1108.6266 [gr-qc].[55] S. Nojiri, S. D. Odintsov, and V. K. Oikonomou, Phys.Rep. , 1 (2017), arXiv:1705.11098 [gr-qc].[56] T. Harko and F. Lobo,
Extensions of f(R) Gravity:Curvature-Matter Couplings and Hybrid Metric-PalatiniTheory , Cambridge Monographs on Mathem (CambridgeUniversity Press, 2018).[57] I. De Martino, R. Lazkoz, and M. De Laurentis, Phys.Rev. D , 104067 (2018), arXiv:1801.08135 [gr-qc].[58] M. De Laurentis, I. De Martino, and R. Lazkoz, Phys.Rev. D , 104068 (2018), arXiv:1801.08136 [gr-qc].[59] D. P´erez, G. E. Romero, and S. E. Perez Bergliaffa, As-tron. Astrophys. , A4 (2013), arXiv:1212.2640 [astro-ph.CO].[60] N. Alipour, A. R. Khesali, and K. Nozari, Astrophysicsand Space Science , 240 (2016).[61] C. W. Misner, K. S. Thorne, and J. A. Wheeler, Gravi-tation , Misner73 (W. H. Freeman, San Francisco, 1973).[62] A. Bouvard,
Tables Astronomiques Publiees Par Le Bu-reau Des Longitudes De France, Contenant Les Ta-bles De Jupiter, De Saturne Et D’Uranus, ConstruitesD’Apres La Theorie De La Mecanique Celeste (Bachelieret Huzard, 1821).[63] C. Will,
Theory and Experiment in Gravitational Physics (Cambridge University Press, 2018).[64] B. Famaey and S. S. McGaugh, Living Reviews in Rela-tivity , 10 (2012), arXiv:1112.3960 [astro-ph.CO].[65] S. Mendoza, X. Hernandez, J. C. Hidalgo, andT. Bernal, Mon. Not. R. Astron. Soc. , 226 (2011),arXiv:1006.5037 [astro-ph.GA].[66] S. Mendoza, Canadian Journal of Physics , 217 (2015).[67] T. Bernal, S. Capozziello, G. Cristofano, and M. deLaurentis, Modern Physics Letters A , 2677 (2011),arXiv:1110.2580 [gr-qc].[68] D. A. Carranza, S. Mendoza, and L. A. Torres, EuropeanPhysical Journal C , 2282 (2013), arXiv:1208.2502[astro-ph.CO].[69] E. Barrientos, T. Bernal, and S. Mendoza, arXiv e-prints, arXiv:2008.01800 (2020), arXiv:2008.01800 [gr-qc].[70] B. P. Abbott, et al. (LIGO Scientific Collaboration, andV. Collaboration), Phys. Rev. D , 104036 (2019),arXiv:1903.04467 [gr-qc]. [71] R. Abbott, et al. (LIGO Scientific Collaboration, andV. Collaboration), arXiv e-prints , arXiv:2010.14529(2020), arXiv:2010.14529 [gr-qc].[72] T. P. Sotiriou and V. Faraoni, Reviews of Modern Physics , 451 (2010), arXiv:0805.1726 [gr-qc].[73] L. Landau and E. Lifshitz, The Classical Theory of Fields (Elsevier Science, 2013).[74] S. Mendoza and S. Silva, arXiv e-prints ,arXiv:2011.04175 (2020), arXiv:2011.04175 [gr-qc].[75] S. Capozziello, A. Stabile, and A. Troisi, Phys. Rev. D , 104019 (2007), arXiv:0708.0723 [gr-qc].[76] I. De Martino, M. De Laurentis, F. Atrio-Barandela,and S. Capozziello, Mon. Not. R. Astron. Soc. , 921(2014), arXiv:1310.0693 [astro-ph.CO].[77] Gravity Collaboration, R. Abuter, A. Amorim,M. Baub¨ock, J. P. Berger, H. Bonnet, W. Brandner, Y. Cl´enet, V. Coud´e Du Foresto, P. T. de Zeeuw,J. Dexter, G. Duvert, A. Eckart, F. Eisenhauer, N. M.F¨orster Schreiber, P. Garcia, F. Gao, E. Gendron,R. Genzel, O. Gerhard, S. Gillessen, M. Habibi,X. Haubois, T. Henning, S. Hippler, M. Horrobin,A. Jim´enez-Rosales, L. Jocou, P. Kervella, S. Lacour,V. Lapeyr`ere, J. B. Le Bouquin, P. L´ena, T. Ott,T. Paumard, K. Perraut, G. Perrin, O. Pfuhl, S. Ra- bien, G. Rodriguez Coira, G. Rousset, S. Scheithauer,A. Sternberg, O. Straub, C. Straubmeier, E. Sturm,L. J. Tacconi, F. Vincent, S. von Fellenberg, I. Waisberg,F. Widmann, E. Wieprecht, E. Wiezorrek, J. Woillez,and S. Yazici, Astron. Astrophys. , L10 (2019),arXiv:1904.05721 [astro-ph.GA].[78] P. K. Townsend, arXiv e-prints , gr-qc/9707012 (1997),arXiv:gr-qc/9707012 [gr-qc].[79] S. Gimeno-Soler, J. A. Font, C. Herdeiro, and E. Radu,Phys. Rev. D , 043002 (2019), arXiv:1811.11492 [gr-qc].[80] T. Johannsen and D. Psaltis, Astrophys. J. , 446(2010), arXiv:1005.1931 [astro-ph.HE].[81] Event Horizon Telescope Collaboration, K. Akiyama,A. Alberdi, W. Alef, K. Asada, R. Azulay, and thers, As-trophys. J. Lett. , L6 (2019), arXiv:1906.11243 [astro-ph.GA].[82] K. Hioki and K.-I. Maeda, Phys. Rev. D , 024042(2009), arXiv:0904.3575 [astro-ph.HE].[83] X. Hou, Z. Xu, M. Zhou, and J. Wang, Journal ofCosmology and Astroparticle Physics , 015 (2018),arXiv:1804.08110 [gr-qc].[84] Event Horizon Telescope Collaboration, K. Akiyama,A. Alberdi, W. Alef, K. Asada, et al. , Astrophys. J. Lett.875