Optical appearance of numerical black hole solutions in higher derivative gravity

The optical appearance of the numerically black hole solutions within the higher derivative gravity illuminated by an accretion disk context is discussed. We obtain solutions for non-Schwarzschild black holes with r0 = 1, r0 = 2, and r0 = 3. Further analysis of spacetime trajectories reveals properties similar to Schwarzschild black holes, while the r0 = 2 black hole exhibits significant differences. The results reveal the presence of a repulsive potential barrier for the black hole, allowing only particles with energies exceeding a certain threshold to approach it, providing a unique gravitational scenario for non-Schwarzschild black holes. Additionally, the optical images are derived through numerical simulations by discussing the trajectories of photons intheblackholespacetime.The distribution of radiation flux and the effects of gravitational redshift and Doppler shift on the observed radiation flux are considered. Interestingly, previous analyses of the optical appearance of black holes were conducted within the framework of analytic solutions, whereas the analysis of numerical black hole solutions first appears in our analysis.


I. INTRODUCTION
Einstein's general relativity (GR) regarded as the standard theory of gravity, but it encounters conflicts with quantum physics in some details.One of the more important issues in theoretical physics is to construct a theory of gravity that is renormalizable.The GR can be understood as an effective theory at low energies.As energy scales increase, higher-order corrections become increasingly significant [1].One possible solution involves incorporating higher-order curvature terms into the standard Einstein-Hilbert action.By including all possible quadratic curvature terms, Stelle constructed a theory amenable to renormalization, which is the higher derivative gravity [2,3].The inclusion of a quadratic curvature term comes with the drawback of introducing a non-physical ghost field [4].Furthermore, in vacuum quadratic gravity within higher derivative gravity, Einstein gravity is equated with some form of stress energy, rendering the Birkhoff theorem invalid [5].Consequently, theoretically, there exist spherically symmetric solutions beyond the Schwarzschild solutions.Several static spherically symmetric black hole (BH) solutions have been discovered in higher derivative gravity, with discussions on the stability of non-Schwarzschild BHs [6][7][8][9][10].Additionally, investigations into charged BH solutions and high-precision analytical approximations of numerical solutions have been undertaken [11,12].
In addition to theoretical investigations, the astronomical detection of real images of black holes holds tremen-dous interest.The Event Horizon Telescope (EHT) Collaboration published an image of the supermassive object at the center of the Messier(M) 87* galaxy, marking the first direct visual evidence for the existence of BHs [13].Subsequently, the EHT collaboration released an image of the BH at the center of the Sagittarius (Sgr) A* [14].The polarized image of M87* BH revealed the structure of magnetic fields and the plasma properties near the BH [15].These works provide enlightening insights into understanding the matter distribution and the electromagnetic interactions around BHs.
The BH shadow, caused by light being captured by the BH, stands as one of the most crucial elements in BH images [16].The shape and size of these shadows are intimately linked to the geometry of the BH spacetime, offering a method for identifying BHs observed through astronomy [17][18][19].Grenzebach et al. derived an analytical formula for the shadow of a Kerr-Newman-NUT-(anti-)de Sitter BH, facilitating the convenient visualization of photon regions and shadows [20,21].The EHT collaboration's tests on the radiation transmission code of the BH images revealed that the error in the observed radiation flux was merely one percent, suggesting the feasibility of parameter estimation using BH observations [22].Numerous detailed studies have delved into gravitational lensing and BH shadows within various gravitational spacetime contexts [23][24][25][26][27][28][29][30][31][32][33][34][35][36].
In real astrophysical observation, observable BHs are usually encircled by an accretion disk.Gravitational energy released from material within this disk emits a characteristic spectrum of electromagnetic radiation, serving as a light source against the background of the BH [37].If the accretion rate surpasses that of other objects in the Milky Way, it can manifest as an active galactic nucleus (AGN) theough Eddington accretion.The redshift of iron lines in the spectra of AGN finds solid explanation through BH accretion disks [38].Research on accre-tion disk models traces back to the 1970s, where Shakura et al. proposed a geometrical thin, optically thick accretion disk model, also known as the standard model [39].Thorne et al. extended this standard model to relativity, establishing the Novikov-Thorne model [40].Buding upon these studies, Luminet obtained direct and secondary images of the disk of a Schwarzschild BH with a thin accretion disk, calculating the accretion disk's brightness using the radiation flux analytical formula derived in [41].Utilizing Luminet's approach, images of BHs and naked singularities within various modified gravitis have been investigated [42][43][44][45][46].
In addition to optically thick, geometrically thin accretion disks, discussions also encopass other types of accretion disks.Gralla et al. provided a precise description of shadows, photon rings, and lensing rings for optically thin disk [47].Pedro et al. discussed the lensing and shadow of a BH surrounded by a heavy accretion disk [48].The cornerstone of these works lies in ray-tracing techniques.However, Luminet's ray-tracing method is only applicable to spherically symmetric spacetime.For axially symmetric BHs, the Newman-Janis algorithm can be utilized to reduce the second-order geodesic equations to a set of first-order geodesic equations by means of a constant of motion [49].However, in cases where the dynamical system of photons is nonintegrable [50][51][52][53][54], resulting in equations of motion lacking a carter constant, the motion of the photons is chaotic, and the disk image can only be obtained through a numerical ray-tracing code [53].
Research on the optical appearance of BHs in the past has been conducted in the context of analytical solutions to BH, which typically apply to simplified BH models.In real astronomical observations, the considered physical backgrounds are complex, involving asymmetry, rotation, or the combination of gravitational fields with other matter, which severely restricts the applicability of analytical solutions.The introduction of numerical BH solutions helps fill this gap, allowing for a more accurate test of the validity of GR.Therefore, the main goal of this article is to discuss the observable features of BH shadows under the backdrop of numerical solutions.The structure of this paper is as follows: in Sec.II, we briefly review the non-Schwarzschild BH numerical solutions in higher derivative gravity.Sec.III analyzes the time-like geodesics and null geodesics of non-Schwarzschild BH.In Sec.IV, we explore the optical appearance of this intriguing numerical BH solution within the framework of a thin accretion disk model, while considering the influences of redshift and observed inclination.Finally, we draw our conclusions in Sec.V.

II. NUMERICAL BLACK HOLE SOLUTION
The most general Einstein-Hilbert action with additional quadratic curvature terms can be expressed as [6] where C µνρσ represents the Weyl tensor, and α, β, and γ denote coupling constants.Since the resulting tensor in the field equations arising from the Weyl tensor is traceless, the equations of motion should exclude the contribution of the R 2 term .Therefore, it is feasible to assign β = 0. Assuming that γ = 1 to streamline the calculation, the equations of motion can be simplified as [6] where the expression B µν = (∇ ρ ∇ σ + 1 2 R ρσ )C µρνσ defines the Bach tensor.It is notable that the Schwarzschild solution satisfies Eq. (2).A numerical solution for a non-Schwarzschild BH with static spherical symmetry was presented in [6].The metric for a static spherically symmetric spacetime can be represented as: Substituting Eq. (3) into Eq.( 2), one obtains the field equations: Obtaining an analytical solution of Eq. ( 4) and Eq. ( 5) is challenging, necessitating a numerical approach for solution.We postulate the presence of a BH event horizon at r = r 0 > 0, where h(r) = f (r) = 0. Near this horizon, one can derive Taylor expansions for f (r) and h(r): (7) in which c can be absorbed into a rescaling of the time coordinate, and f i and h i denote the constant coefficients characterizing h(r) and f (r) in the vicinity of the event horizon.We designate h 1 = f 1 to accommodate the arbitrary rescaling of the time coordinate.Setting α = 1 2 without loss of generality allows for the calculation of all coefficients h i and f i with i > 2 from two nontrivial free parameters, namely r 0 and f 1 .The parameter values are presented in Table I.Utilizing numerical routines in MATHEMATICA, we integrate the equations.It's crucial to emphasize that r 0 needs to be determined  first, and subsequently, a specific f 1 is obtained through a shooting method to ensure the integration remains free of singularities over a large radius.Figure 1 illustrates the behavior of h(r) and f (r) for various values of r 0 .In the analysis by [6], a critical event horizon radius r c ≈ 1.143 is identified.For r 0 < r c , both h(r) and f (r) are increasing functions, while for r 0 > r c , h(r) and f (r) display extremal points.In both two cases, the metric function f (r) will approach 1 when r → +∞, so the non-Schwarzschild solutions are asymptotically flat.Furthermore, the mass of the non-Schwartzschild black hole decreases as r 0 increases, and becoming negative when r 0 > r c [6].We will examine the class of r 0 < r c and r 0 > r c by taking r 0 as 1 and 2 respectively.

III. GEODESICS OF HIGHER DERIVATIVE GRAVITY BLACK HOLE A. Time-like Geodesic
By leveraging the spherical symmetry of spacetime, we confine the particle motion to the equatorial plane (θ = π 2 ) without sacrificing generality.As the metric lacks explicit dependence on the coordinates t and φ, two conserved quantities governing the motion emerge: where τ represents the proper time of the particle, while E and L denote the energy and angular momentum, respectively.The time-like geodesic is governed by ds 2 = −dτ 2 , By Substituting Eq. ( 8) and Eq. ( 9) into the geodesic equation, one can derive the equations governing radial motion and orbital dynamics for particles: Taking the derivative of both sides of Eq. ( 10) with respect to τ , one obtains the expression for radial acceleration: Considering the orbit of the particle on the accretion disk as a circular orbit, we initially examine the scenario of such orbit.For particles in a circular orbit, both radial velocity and acceleration are zero ( ṙ = r = 0).Hence, one can derive expressions for E and L in terms of r: As energy and angular momentum must be real numbers, the conditions for a circular orbit are 2h(r) − rh ′ (r) > 0 and h ′ (r) > 0. Utilizing the numerical values of BH calculated in Sec.II, one can determine the range of circular orbits for non-Schwarzschild BH, i.e, r r0=1 > 1.437, 2.666 < r r0=2 < 3.996.(14) Moreover, the stability of a circular orbit necessitates V ′′ (r) ≤ 0. Our calculations for non-Schwarzschild BH reveal that the radius of the inner stable circular orbit is approximately r isco ≈ 3.503.However, in cases where r 0 > r c , our calculations demonstrate that V ′′ (r) > 0 within the range of conditions for circular orbit.This indicates the absence of unstable circular orbits for BH with r 0 > r c .For general particle orbits, to simplify the description of radial motion, we introduce the particle's effective potential V ef f , defined as the energy with zero radial velocity.This potential can be derived from Eq. ( 10) Figure 2 illustrates the radial distribution of the effective potential for two types of BHs with r 0 = 1 and r 0 = 2.In the case of r 0 = 1, there is minimal distinction between the Schwarzschild BH and the non-Schwarzschild BH counterpart.For the Schwarzschild BH, extreme points in the effective potential emerge only when L > √ 3.This implies that particles with low angular momentum moving towards the BH are inevitably captured, while those with high angular momentum must overcome a potential barrier to approach the BH.However, in the case of r 0 = 2, the non-Schwarzschild BH exhibits a distinct behavior.Even when angular momentum is negligible, the potential function still features a unique extreme value V max .Consequently, particle with any angular momentum outside the corresponding radius of this extreme point must traverse a potential barrier to approach the BH.
To illustrate the contrast between the orbits of two types of BHs with r 0 = 2, Fig. 3 depicts the trajectories of unbound particles with varying energies around the BHs, assuming L = 6.The plot reveals the repulsive effect exerted by the non-Schwarzschild BH on distant particles.In particular, if the radial velocity of the particle is 0, Eq. ( 12) can be reduced to The sign on the right side of Eq. ( 16) denotes the direction of the BH force.It's evident that r > 0 holds true whenever h ′ (r) < 0, indicating the effective repulsive nature of the BH on particles.For a particle in a circular orbit, its radial velocity is zero, hence its acceleration can be described by Eq. ( 16).This elucidates why there exists an upper limit for the orbit radius of the BHs with r 0 = 2 and r 0 = 3 in Eq. (14).Beyond this limit, the force from the BH is repulsive, thereby preventing the formation of circular orbits.

B. Null Geodesic
In our earlier examination, we determined that a non-Schwarzschild BH with r 0 > r c lacks a stable circular orbit, thereby precluding the formation of a stable accretion disk.Consequently, our discussion regarding the accretion disk's image solely focuses on the scenario where r 0 = 1.Analogous to the scrutiny of time-like geodesics, we delve into the trajectory of photons within the equatorial plane: where λ represents the affine parameter, E denotes the energy, and L signifies the angular momentum of the photon.The null geodesic adheres to the condition ds 2 = 0. Following the approach employed in defining the effective potential for particles in the preceding section, we define the effective potential for photons as: For photons, the angular momentum L solely influences the magnitude of the effective potential, without altering its behavior.The radius of the photon sphere is determined at points where V ef f = 0 and V ′ ef f = 0. Hence, One can derive the radius of the photon sphere sphere and the corresponding impact parameters by: where b = L E represents the impact parameter for photons at infinity.Subsequently, incorporating Eq. ( 17) and Eq. ( 18) into ds 2 = 0, and eliminating τ from the equation, yields the orbital equation: To obtain the deflection angle of light, integrate r in the orbital equation (22): In this equation, γ 1 denotes the change in the φcoordinate of a photon from the radius r sourse to infinity.It's worth noting that for the photon passing the perihelion, the expression for the deflection angle is: where r p represents the radius at perihelion.

IV. OBSERVABLE FEATURES OF THE THIN ACCRETION DISK
In this section, we consider a scenario where the BH is encircled by a geometrically thin, optically thick accretion disk.We utilize the ray-tracing method outlined in [41] to explore the direct and secondary images of the BH.The schematic diagram of the ray-tracing method is depicted in Fig. 5, with detailed geometric relationships elaborated in [54].Here, the BH is situated at point o, the accretion disk lies in the xoy plane, the photographic plate at infinity is positioned in the x ′′ o ′ y ′′ plane with an inclination angle θ 0 , and light emanates from point M (r, φ) on the accretion disk, arriving at point m on the observation plane in the progressive direction oo ′ .Our primary focus lies on the direct image (b (d) , α) generated by photons facing the photographic plate and the secondary image (b (s) , α + π) generated by photons facing away from the photographic plate.

A. Ray-tracing
As depicted in Fig. 5, the photon traverses the oo ′ M plane, and the deflection angle γ precisely corresponds to ∠M oo ′ .By applying the sine theorem of spherical triangles to △M yy ′ and △M xx ′ , one can derive: In Section III, we utilized numerical integration to illustrate the relationship between the deflection angle and the impact parameter b as outlined in Eq. ( 23) and Eq.(24).In Luminet's analysis, the deflection angle is expressed through elliptic integration.However, our calculations concerning Schwarzschild BHs demonstrate that both numerical and elliptic integration yield nearly identical results within a highly precise range.Hence, we employ the numerical integration method to compute the deflection angle for non-Schwarzschild BHs.For photons corresponding to the direct image and higher-order images, one can derive from [42]: in which the parameter n signifies the order of the image.Given that higher-order images closely resemble secondary images, we solely focus on secondary images where n = 1 [42].
In Section II, we calculate the non-Schwarzschild solution for α = 1 2 .Despite α and r 0 are not independent [6], we can still obtain a solution with r 0 = 1 in a small neighborhood around α = 1 2 .Figure 6 shows the direct and secondary images of non-Schwarzschild BH with α = 7 16 , α = 1 2 and α = 9  16 , We can see that as α increases, the secondary image gradually becomes larger, while the direct image changes slightly.For the sake of convenience in calculations, we still consider the case of α = 1 2 to investigate the differences between non-Schwarzschild BHs and Schwarzschild BHs.
Figure 7 illustrates the direct and secondary images of two types of BHs.It's evident that the observed viewing angle profoundly influences both the direct and secondary images.When the viewing angle is small, the secondary image resembles the direct image, forming a structure akin to the "photon ring".As the viewing angle increases, the direct image assumes a "hat" shape, and the secondary image gradually separate.We observed that while the direct images of the two types of BHs exhibit minimal disparity, the secondary images of non-Schwarzschild BH are notably less separated compared to Schwarzschild BH.

B. Observed Flux
In this subsection, we will employ the Novikov-Thorne model to analyze the radiation flux of the accretion disk.The expression for the radiation flux is provided in [40]: Where Ṁ denotes the mass accretion rate, g represents the metric determinant, and r in signifies the radius of the inner edge of the accretion disk.E, L, and Ω respectively stand for the energy, angular momentum, and angular velocity of the particle in circular orbit.
It's important to note that while the Novikov-Thorne model is originally formulated for Kerr BHs, several studies have adapted this model for spherically symmetric black holes [41][42][43][44][45].Only when g tφ = 0 for spherically symmetric black holes, modifications are required for the expressions of E, L, and Ω.Under the metric of spherical symmetry ds 2 = g tt dt 2 + g rr dr 2 + g θθ dθ 2 + g φφ dφ 2 , the energy, angular momentum, and angular velocity of particles on the equatorial plane are respectively: To simplify calculations, we utilize the maximum radiation emitted by the Schwarzschild BH as the unit of radiation flux.Fig. 8 illustrates the radiation flux distribution of two types of BHs at different viewing angles.It's observable that the radiation from the accretion disk approaches zero at the inner edge of the disk, and the radiation flux initially rises and then decreases with increasing radius.Although the radiation distribution patterns of the non-Schwarzschild BH align with those of the Schwarzschild BH, the corresponding radiation values are smaller for the non-Schwarzschild BH.
For observers positioned at infinity, both gravitational redshift and Doppler shift effects influence the observed radiation flux.As per the findings presented in [56], the observed radiation flux corresponds to the ratio of the radiation flux to the fourth power of the redshift factor (z + 1).
The redshift factor is defined by the ratio of the energy E em at photon emission to the energy E obs observed at infinity.The energy at photon emission can be expressed in terms of the momentum of the photon p and the velocity of the emitting particle u: For a distant observer, the ratio p t /p φ precisely represents the impact parameter of the photon at infinity with respect to the z axis.By combining Eq. ( 29) and Eq. ( 33), one can derive the expression for the redshift factor: Figure 9 displays the redshift distribution at various viewing angles.It's notable that, besides the symmetric gravitational redshift, the overall redshift distribution exhibits asymmetry due to the rotation of the accretion disk.On the side rotating towards the observer, the redshift is diminished, and may even shift towards blueshift, whereas on the other side, the redshift effect is amplified.Additionally, the viewing angle influences the redshift distribution, with a more pronounced redshift effect observed as the viewing angle increases.The redshift distribution patterns of the two types of black holes are fundamentally similar, although the effect is slightly weaker for non-Schwarzschild BHs.
Considering both the radiation flux and redshift of the accretion disk, we can deduce the observed flux distribution of the accretion disk (refer to Fig. 10 and Fig. 11).Subsequently, we simulate the images of the two types of BHs based on these observed flux distributions (Fig. 12).It's evident that for both BHs with r 0 = 1, the images exhibit strong asymmetry, and the brightness distribution patterns are similar.However, notable differences arise: the non-Schwarzschild BHs appear weaker, and their secondary images are smaller compared to those of the Schwarzschild BHs.These distinctions could po-tentially serve as observable features for distinguishing between the two types of BHs.

V. CONCLUSIONS
In this analysis, we explore the optical characteristics of non-Schwarzschild BH surrounded by thin disks within higher derivative gravity, identifying observable features that distinguish them from Schwarzschild BHs.We replicate the investigations conducted in [6], obtaining solutions for non-Schwarzschild BHs with r 0 = 1, r 0 = 2, and r 0 = 3. Subsequent analysis of time-like geodesics reveals that the geodesic properties of non-Schwarzschild BHs with r 0 = 1 closely resemble those of Schwarzschild BHs, whereas non-Schwarzschild BHs with r 0 = 2 exhibit distinct characteristics.Notably, particles with any angular momentum, even zero angular momentum, must traverse a potential barrier to approach the BH from a distance, suggesting a effective repulsive force exerted by the BH in this scenario.To illustrate this property effectively, we present the unbound orbit diagram of particles with varying energies for the case r 0 = 2 (Fig. 3), revealing that only particles with energies greater than the maximum potential energy V max can reach the BH.Furthermore, numerical results regarding the stability of circular orbits indicate that BHs with r 0 > r c lack stable time-like circular orbits, precluding the formation of accretion disks around the BH.Consequently, observing such BHs may prove challenging.
For the case of r 0 = 1, we conduct further analysis on null geodesics within the BH spacetime, calcu-   lating the photon sphere radius and the corresponding critical impact parameters of non-Schwarzschild BHs using the effective potential.Additionally, we employ numerical integration to determine the deflection angle of photons around Schwarzschild BHs, and plot the trajectories of photons around both Schwarzschild and non-Schwarzschild BHs (Fig. 4).To visualize the BH images, we utilize the ray-tracing method proposed by Luminet [41] to obtain the direct and secondary images of the black hole accretion disk (Fig. 7).Our findings indicate that while there's no significant difference between the direct images of the two types of BHs, the degree of separation of secondary images for non-Schwarzschild BHs diminishes compared to Schwarzschild BHs as the radius increases.Furthermore, in line with previous research, the primary image of non-Schwarzschild BHs also exhibits a "hat" shape with increasing viewing angles, and the separation of secondary images becomes more apparent.
Furthermore, we take into account the radiation emitted by the accretion disk and calculate the distribution of the radiation flux using the Novikov-Thorne accretion disk formula from [40] (Fig. 8).Our results indicate that while the distribution patterns of the radiation flux for both types of black holes are essentially the same, the radiation intensity of non-Schwarzschild BHs is compar-atively weaker.Additionally, gravitational redshift and the Doppler effect influence the observed radiation flux.By combining the radiation flux distribution with the redshift distribution of the accretion disk, we calculate the observed radiation flux distribution (Fig. 10 and Fig. 11).Finally, we conduct numerical simulations of the observed radiation flux to generate optical images of the two types of BHs (Fig. 12).The results also reveal that the brightness distribution patterns of the two types of BHs are similar.However, due to the redshift effect, both types of images exhibit asymmetry, with the degree of asymmetry increasing with the viewing angle.Notably, we observe that non-Schwarzschild BHs appear dimmer, and the secondary images of non-Schwarzschild BHs are less separated compared to those of Schwarzschild BHs.These differences may serve as observable features to distinguish between the two types of BHs.

FIG. 1 .
FIG.1.Numerical solutions for f (r), h(r) in the case of r0 = 1, r0 = 2 and r0 = 3.The dashed black line represents the unity.In each plot, We set c =3  4 to separate h(r) from f (r).

FIG. 3 .
FIG.3.The particle trajectories of two types of BHs with r0 = 2.The BHs are shown as the black disks, and the dashed blue lines represent the circular orbit.

FIG. 4 .FIG. 5 .
FIG. 4. The light trajectories of two types of BHs with r0 = 1, the BHs are shown as the black disks, and the yellow dashed circle represent the circular orbit.

Fig. 4 15 ¡ = 9 / 16 , 2 .
illustrates the photon orbit of two types of BHs.For b > b c , the green line corresponds to direct image, the orange line corresponds to lensed ring and the red line corresponds to photon ring.For b > b c , The cyan line, blue line, purple line represent the direct image, lens ring, and photon ring, respectively.[55] It's observable that light exhibiting the same deflection around non-Schwarzschild BHs possesses slightly smaller impact parameters.¢ = 80°F IG. 6.The direct (solid line) and secondary (dashed line) of non-Schwarzschild BHs with the viewing angles θ0 = 17 • , 53 • and 80 • .Left panel: α = 7 16 .Middle panel: α = 1 Right panel: α = 9 16 .

TABLE I .
The numerical value of fi and hi (i ≤ 3).