Propagation of light in the presence of gravity generated by static and spherically symmetric curved space-times using Maxwell equations

In this manuscript, we present an alternative method for calculating null geodesics in General Static Isotropic Metrics in General Relativity and Extended Theories of Gravity. By applying a conformal transformation, we are able to consider an analogue gravity model, where curvature is encoded in the dielectric and magnetic properties of a medium. In other words, we pass from curved to flat space-times, where instead of the Einstein field equations, the Maxwell equations are solved. Within this geometrical background, the photon geodesics are calculated. Then, given different black hole and wormhole metrics, we apply this method obtaining an excellent agreement with respect to the exact solutions in the original gravity framework by committing angular deviations below $3^\circ$. Finally, we provide the image of a Schwarzschild black hole surrounded by a thin accretion disk, and the apparent image of a Morris&Thorne-like wormhole within an angular discrepancy below $4^\circ$.


Introduction
When an electromagnetic (EM) field interacts with a strong gravitational field, its characteristics (such as frequency, amplitude, and propagation direction) are affected by the curvature of the space-time. The great interest in studying the behavior of the EM waves' propagation in different gravitational fields motivated the efforts of some research groups towards the development of EM-gravity analog models conformally settled in flat space-times in order to investigate grava e-mail: legarcia@ing.uc3m.es (corresponding author) ity through an alternative point of view and with very handy tools, which can be reproducible and built up in laboratory. In order to inquiry the curvature properties of the underlying space-time we study the geodesic orbits followed by both photons and massive bodies when they are in free-fall and only gravity is acting on them, neglecting thus other minor perturbing effects.
The first idea of this conformal EM-gravity analogue pattern can be traced back to Eddington, who in 1920 thought to describe the light deflection effects from the Sun's gravitational field by treating it as it was a refractive medium [1]. This concept was better mathematically worked out in 1960 by Plebanski [2], who presented, for the first time, an EM model of the gravitational field based on a set of suitable dielectric permittivity and magnetic permeability tensors.
There are several international research groups aiming at reducing the computational times for calculating the photon trajectories starting from static and spherically symmetric space-times. These developments permit to inquire gravity and compact object physics, such as black holes (BHs) and wormholes (WHs), via numerical imaging, where theoretical results can be benchmarked with those furnished by the Event Horizon Telescope Collaboration [3,4,5,6,7].
In 1979 Luminet provided the first ever numerical image of the optical appearance of a Schwarzschild BH surrounded by a thin accretion disk [8]. Beloborodov [9], La Placa et al. [10], and Poutanen [11] presented different and very accurate approximated expressions of photon geodesics in the Schwarzschild geometry. Also De Falco and collaborators approximated photon geodesics in Schwarzschild BH, but introducing a mathematical arXiv:2303.06085v1 [gr-qc] 10 Mar 2023 2 method to obtain polynomials [12] and then tested different theories of gravity through ray-tracing procedures [13]. Bao [14] and Aldi [15] employed geodesics in the estimation of emission lines from the accretion disk around a BH. Johannsen et al. [16,17] introduced a framework to test the no-hair theorem via electromagnetic observables (i.e., spectrum emission lines from accretion disks, measurement of the Inner Most Stable Orbit (ISCO), and BH images), which they estimated by numerically integrating the geodesics. Bakala et al. [18] computed the optical appearance of a Schwarzschildde Sitter BH. Fernández-Núñez et al. [19] applied a conformal transformation to mimic the Schwarzschild BH space-time with metamaterials. In [20], the authors presented and validated a method for computing null geodesics in General Static Isotropic Metrics (GSIMs) by applying a Conformally Flat Space-time Transformation (CFST), based on the Snell's refraction law and isotropic coordinates, particularizing therefore the Plebanski's results to isotropic EM materials.
Regarding the last work, in this paper we aim to extend those goals by applying the aforementioned strategy to different BH and WH metrics by also discussing about the numerical approximation errors and calculating the apparent image of an Ellis WH. Our methodology might represent an interesting tool to understand and study gravity via feasible EM apparati. This article is organized as follows: in Sec. 2 the theoretical background and the proposed ray-tracing algorithm are presented; in Sec. 3, we show the outcomes of our approach; finally in Sec. 4, we discuss the obtained results and draw the conclusions.

The General Static Isotropic Metric
In General Relativity (GR), the features of a spacetime geometry are encoded in the metric tensor g αβ . The GSIM describes static and isotropic gravitational fields, namely, they are independent from the coordinate time x 0 = t and depend only on the radial coordinate [21]. Their line element, ds 2 = g αβ dx α dx β (where Einstein' summation convention is used), expressed in spherical coordinates (t, r, θ, φ), geometric units (c = G = 1, where c is the speed of light in vacuum and G is the gravitational constant), and signature (−, +, +, +) reads as ds 2 = −g tt (r)dt 2 + g rr (r)dr 2 + r 2 (dθ 2 + sin 2 θdφ 2 ), (1) where g tt (r), g rr (r) are two unknown positive functions (depending only on the radius, r), whose explicit ex-pressions allow to determine the background space-time geometry.

Photon geodesic equations in GSIM
In GR, the null geodesic trajectories represent the minimal path connecting two points followed by photons on the space-time background. They can be obtained via the Euler-Lagrange equations d dλ ∂L ∂ẋ α − ∂L ∂ x α = 0, where x α are the coordinates, 2L = g αβ dx α dλ dx β dλ is the Lagrangian function, and λ is the affine parameter along the null trajectories [22]. Due to the spherical symmetry of the GSIM, the photon energy E and angular momentum L (along any direction) are conserved along the photon orbit. The photon moves in a single invariant plane [12], characterized by the initial position and spatial velocity vector of the photon. Without loss of generality, we can settle the photon dynamics in the equatorial plane, i.e., θ = π/2. Its orbital equation in a GSIM, derived from the Euler-Lagrange equations, has the following expression where b = L/E is the photon impact parameter, which can be defined as follows with α being the emission angle. The so-called critical impact parameter, b c , is obtained by Eq. (3) substituting r = r ps and α = π/2, where r ps is the photon-sphere radius (obtained from the effective potential V Eff = −g tt (r)/r 2 by imposing dV Eff dr = 0) [13]. We can classify three possible families of null-geodesic trajectories depending on the values of b with respect to the critical photon impact parameter b c : (1) b < b c , which include all photons coming from infinity and ending their motion inside the compact object's boundary; (2) b = b c , corresponding to the "critical orbits", where the photons from infinity conclude their motions on a stable circular orbit around the compact body, namely the photon sphere r ps ; (3) b > b c , encompassing all photons coming from infinity, being deflected by the compact object's curvature, and sent again to infinity. As a final remark, not all GSIM-like metric may allow photons to move along critical orbits. 3

Conformally flat spacetime transformation
A CFST [21,19,20] maps the metric (1) into another being isotropic, whose explicit expression is (4) where H(ρ) = −g tt (r(ρ)), J(ρ) = r(ρ) 2 /ρ 2 , and the radial coordinate r changes into a new suitable radial coordinate ρ, whose transformation is where C o is a constant that can be computed knowing that J(ρ) = 1 for r → ∞, for the asymptotically flat condition. As a final remark, this conformal transformation applies only to GSIM space-times.

Analogous Electromagnetic-gravity Model
The EM approach for treating the GSIM gravitational fields is based on the analog-gravity model derived in [2]. Since GSIM-like metrics are diagonal in spherical coordinates, the equivalent dielectric permittivity and magnetic permeability tensors are related to the metric tensor via where g = det(g µν ). This formalism can be exploited when the EM field amplitudes are low enough such that their effect on the metric tensor is negligible [2]. When the metric is expressed in its isotropic form (4) thanks to the CFST transformation (5), the dielectric permittivity and magnetic permeability tensors assume the following form where δ ij is the Kronecker delta and n(ρ) being the isotropic refractive index, which allows to rewrite the metric (4) as [23] Equation (9) represents the EM-length line element for calculating the optical distance between two neighboring points for a Transversal-ElectroMagnetic (TEM) mode of propagation.

Ray-tracing technique in analogous electromagnetic model spacetimes
We compute the photon geodesics in GSIM-like metrics model by exploiting the conventional ray-tracing scheme via the Snell's refraction law. In Fig. 1, we report a visual sketch of the strategy adopted in this article. Fig. 1 Scheme of the method employed to discretize a null geodesic path via straight segments and the Snell's refraction law. The colored circular regions represent a generic GSIM equivalent spacetime described in an isotropic coordinate system through the EM-analog model used in this work. The figure shows how a null geodesic in the inhomogeneous media is approximated by using discrete steps when a photon moves from a region of low refractive index towards a region of high refractive index.  Fig. 2 Flowchart of the algorithm employed to compute the photon trajectories based on the Snell's refraction law.
In Fig. 2, we sketch the flowchart of the approximation method to calculate photon trajectories in GSIM space-times proposed in this work. Firstly, we set the initial position (r i , φ i ) and momentum κ of a photon emitted with an emission angle α in the standard spher-4 ical coordinate system or GSIM space-time (see Fig. 2). The momentum, as a function of the initial conditions, is given by Then, as shown in Fig. 2, we transform the initial conditions in spherical coordinates: the initial position (r i , φ i ), and the wave vector ( κ = k r e r +k φ e φ ) to isotropic coordinates: (ρ i , φ i ), and ( κ = k ρ e ρ + k φ e φ ), respectively, by using the Conformally Flat Space-time Transformation. The azimuthal angle φ i and the azimuthal component of the wave vector are not altered, because the coordinate transformation is conformal. Instead, the radial component of the wave vector transforms as As shown in Fig. 1, the ray-tracing algorithm can be adapted to any GSIM space-time (9) considering the movement of photons in an inhomogeneous medium and using the Snell's refraction law, n i−1 sin(φ) = n i sin(θ), where φ and θ are the incidence and refraction angles, respectively, and n i−1 and n i are the refractive indices at the point i − 1 and i, respectively (see inset in Fig.  1), to compute how the photon's trajectory is refracted because of the medium inhomogeneity. In the following, we shall explain the main characteristics of the solving-technique based on Snell's law. Let us consider a photon moving, such as in Fig. 1, from a point i − 1 to a point i. We divide the medium into three regions separated by radial distances ρ i−1 and ρ i , respectively, where, for simplicity, we assume that ρ i < ρ i−1 as reported in Fig. 1. We have that ρ i−1 is the initial radial position and ρ i = ρ i (∆µ i ) is the arrival radial position, where ∆µ i is the straight-length followed by the photon along its wave vector direction. We state that the refractive index inside the region ρ i < ρ < ρ i−1 is a constant defined by n(ρ) = n(ρ i−1 ) since we consider very small displacements. The photon's movement gives rise to the so-called electrical displacement (∆µ ele ), which provides information on the phase difference between two points, obtained by integrating Eq. (9) as follows: where λ is an affine parameter along the photon trajectory, useful for the integration process. Equation (12) is a line integral along the aforementioned photon straightline trajectory. We define ∆µ max as the maximum electrical displacement in order to obtain a good resolution and low approximation error (after having fixed an a-priori tolerance). Between two points, the smaller the electrical displacement is, the smaller the variation of the refractive index will be. This ensures that the solving technique based on Snell's law provides a reliable approximation. Therefore, at each step, we check whether the electrical displacement is larger than ∆µ max , otherwise ∆µ i (defining the photon's length-element along the wave vector direction) is reduced until the condition ∆µ ele < ∆µ max is met. After having determined the optimum ∆µ i , we can then compute the arrival radial distance ρ i and the related refractive index in that point, n(ρ) = n(ρ i ). We calculate the refracted wave vector to be used in the next photon displacement by invoking the Snell's refraction law. This process continues until the photon reaches a maximum number of displacements N or when the photon reaches the spatial infinity, namely the limits of the computational domain. The adaptability of the displacement distance ∆µ i at each propagation step is mandatory provided that the isotropic refractive index, n(ρ), is non-homogeneous and the optimum ∆µ i is not known at the beginning of the displacement where it is set to a default value. Once the approximation finishes, the Conformally Flat Space-time Transformation is inverted to come back to the standard (spherical) coordinate system.
This algorithm must be slightly revised for WHs, since the photons can traverse them. Indeed, care must be encompassed during the integration process on the WH throat (boundary between universe 1 and universe 2). In this region, we use as initial position and fourvelocity for the universe 2 the arrival coordinates obtained in the universe 1, together with the radial component of four-velocity inverted, while maintaining fixed the angular component. Physically, this corresponds to the angular momentum conservation. To compute the arrival coordinates, we allow the photon to move towards the origin of the universe 1 closer to ρ ≈ 1.01ρ min (we set this condition heuristically, and it can adjusted in order to improve the precision), which can dub this point the "pre-arrival" point. Then, we compute the arrival coordinates by extending the pre-arrival coordinates to the WH throat by intersecting a straight trajectory along the pre-arrival wave vector with the WH mouth. This approximation may give rise to high deviations when going from one universe to another if the refractive index has strong variation near the WH mouth. This way of approximating null geodesics works very well for all kinds of null geodesic in a GSIM (see Sec. 2.2). However, due to the discrete nature of the algorithm and instability of critical geodesics, particular attention must be used in the calculations of these orbits. 5 Y,X ′ X i n v a r i a n t p l a n e accretion disk intersection accretion disk-invariant plane Sketch of image generation of a thin accretion disk around a BH. We define r min ≡ rs = 2M as the Schwarzschild radius, r i and re are the inner and outer radii of the accretion disk, respectively. d plate is the distance between the compact object and the center of the photographic plate, ∆Lp is the pixel squared-size, r pixel is the distance between the pixel center and the origin of the observer's plate, Φap is the rotation angle of the invariant plane that passes through the pixel center around the Y axis, and Φ d is the rotation angle of the accretion disk around the X axis. The red trajectory is a null geodesic computed in the Schwarzschild metric impinging the accretion disk in a point, which is emitted over there isotropically. n 1 and n 2 are two refraction indices, φ 1 and φ 2 are the incident and the refraction angles, respectively.n is the normal vector to the circles located in the invariant plane. The computational domain is set in the celestial sphere. This figure is taken from Ref. [20].
As already mentioned in Sec. 2.2, critical orbits in GSIM-like metric spiral around the WH mouth or BH event horizon stabilizing their motion around the circular photon orbit. In this approach, the photon critical orbits are difficult to be approximated, because they can easily deviate from their trajectories when reach the photon sphere, entailing thus large errors. Therefore, in order to avoid such an effect, we implement an analytical extension to the algorithm, which approximates the final part of the orbit (i.e., when the photons start spiralling around the photon sphere) through a circular trajectory.

Image generation of compact objects in curved spacetimes
In this section, we show how to produce the image of a compact object described by a GSIM framed in an analogue EM-gravity model. In Fig. 3, we depict the geometrical setup for imaging a thin accretion disk formed around a static and spherically symmetric compact object.
We consider that each element of the accretion disk emits X-ray EM radiation isotropically in its rest frame. Normally, we should compute an infinite number of photon trajectories from all the emission points coming from each geometrical disk element and then find those that reach the observer's photographic plate located at infinity. We follow instead the inverse procedure, consisting in shooting the photons from the observer's location towards the emission point . We make use of a collision detection algorithm to know where a light ray is being emitted from. The photon trajectories are colored according to the element they come from.
The photographic plate is divided into squared pixels of lengths ∆L p r min , and it is placed at a distance d plate from the center of the compact object such that the gravitational effects are negligible 1 . In practice, this corresponds to have |g rr − 1| < 6% (i.e, asymptotic flatness condition set in our simulations). The protocol for launching the photons relies on shooting them towards the compact body from the center of each pixel with normal direction. Since we are in the invariant plane (X − Y ), as reported in Fig. 3, we must consider only the initial radius and azimuthal velocity components on the invariant plane under consideration given by In order to speed up the integration process we exploit the spherical symmetry of the problem and the fact that photons are launched with a direction perpendicular to the photographic plate. If we look at Fig.  3, we can see the point the photon is being launched from to follow the red-color geodesic high-listed as a blue dot. If we rotate this point on the photographic plate around the Y axis, we will find out that there are at least three more pixels on the photographic plate (high-listed in green), at which photons will be launched with the same initial conditions in the local coordinate system on the the corresponding invariant plane X , Y (i.e., φ i , r pixel , and α i ). The geodesics at green points can be obtained from the geodesic at the blue point just by applying a rotation of coordinate system.

Results
In Fig. 4, we display the geodesics computed in a spacetime generated by a Schwarzschild BH [24,19] with r s = 2. Specifically, Fig. 4(a) shows in an embedding diagram [25] (gray colored mesh) geodesics obtained following the presented approach (colored solid lines) and exact solutions calculated by numerically solving the orbital equation (dashed black lines) for different initial emission angles. In order to measure the convergence of our approach, we plot the angular deviation, ∆θ, measuring the absolute angular error between the final ends of the geodesics computed via our method and the exact solutions. For Fig. 4(a), it amounts to ∆θ ≤ 0.33 • .
Instead in Fig. 4(b), we plot the angular deviation when photons are launched in a Schwarzschild BH spacetime with emission angles in the range 120−240 • . It can be seen that around the critical emission angles α crit1,2 the angular deviation is lower than 2.3 • , and when the emission angle is outside the range α crit1,2 ± 2 • , ∆θ ≤ 1 • . We expected this result due to the discrete nature of the algorithm and the instability of the critical orbits and trajectories for photons emitted with impact parameter near to the critical impact parameter, gravitational source, it is normally considered infinite in the theory for easing the ensuing calculations.  b c . This behavior can be improved by increasing the discretization density. As it can be noted from Fig. 4(a), photons with emission angles larger than the critical one do not fall into the BH, whereas, for critical orbits (cyan solid line), photons are trapped on the photon sphere region. In Fig. 4(a), for critical orbits, the numerical approximation is analytically extended using a circular trajectory as explained in Sec. 2.5, and photons with emission angles lower than the critical one fall into the BH event horizon. Now we focus our attention on geodesics computed in the Ellis WH space-time [25] setting the units for the WH throat as b 0 = 1, see Fig. 5. In particular, Fig. 5(a) shows in the embedding diagram (gray grid) the geodesics computed with our approach (solid lines) and the exact solutions (black dashed lines), having ∆θ ≤ 0.63 • in this case. Instead in Fig. 5(b), we report the angular discrepancy distribution for the emission angle in the range 120 − 240 • . In this case, the maxi-  mum ∆θ is located around α 1 = 170 • and α 2 = 190 • with ∆θ ≈ 3 • . Outside the range α 1,2 ± 2 • , ∆θ ≤ 1 • . As shown in Fig. 5(a), photons emitted with an emission angle larger than the critical one propagates in the same region (universe 1), whereas photons with emission angle shorter than the critical one go through the WH mouth into the so-called universe 2. The photon sphere of the Ellis WH coincides exactly with the WH mouth. Photons emitted with a critical emission angle (cyan solid line) are trapped in the photon-sphere.
We then consider a Morris & Thorne-like WH solution framed in a metric theory of gravity (see Ref. [13] and references therein for more details) endowed with (b 0 , γ) = (0.5, 1.5), see Fig. 6. As done for the previous cases and adopting the same color and style conventions, we first show the geodesics obtained via our method and those obtained integrating exactly the orbital equation using ∆θ ≤ 12.8 • (see Fig. 6(a)). Then, we highlight the angular discrepancy profile for emission angles in the range 120 − 240 • (see Fig. 6(b)). In the last case, the maximum ∆θ is located around α 1 = 144 • and α 2 = 216 • with maximum deviation of ∆θ ≈ 308.51 • . Outside the range α 1,2 ± 12 • , ∆θ ≤ 4 • . We note that in this situation, the angular deviation is considerably larger than the previous cases. This anomaly is related to the fact that the WH metric admits a singular point near the WH mouth and there the refractive index blows up.
Finally, we compare the numerical simulation of an accretion disk around a BH provided by Igor Bogush [26] (see Fig. 7(a)) with the one produced employing our approach (see Fig. 7(b)). The colors of the accretion disk in Fig. 7(b) refer to the observed bolometric flux, as implemented in [8] and [27]. This numerical simulation shows the characteristic gravitational redshift effects, where matter moving away from the observer 8 is redshifted, whereas that approaching the observer is blueshifted. In addition, the BH shadow (black circular spot surrounded by the BH accretion disk), as seen from infinity, agrees with the theoretical prediction (red dashed circle with radius r sh = 3 √ 3r s /2) as reported in [28]. In Fig. 8, we concentrate on the apparent image of an Ellis WH with b 0 = 5 as seen from the universe 1, where the observer is located. We can see, outside the WH mouth (white dashed circle) how the space-time is distorted because of gravity, whereas inside the mouth, a view of the universe 2 appears.
(a) (b) Fig. 7 Apparent image of a Schwarzschild BH surrounded by a thin accretion disk. We set the BH mass M = 1, rs = 2, r skyph = 60rs (computational domain limit or infinite), ∆Lp = rs/10, d plate = 40rs, r i = 3rs (internal radius of the accretion disk), re = 15rs (external radius of the accretion disk), and θo = 80 • (inclination angle of the observer). The red and white dashed circles represent the photon sphere radius projected at infinity (BH shadow), and the BH event horizon rs, respectively. The lateral colored scale represents the normalized time averaged energy flux of matter present in the accretion disk. Panel (a): Apparent numerical image of a Schwarzschild BH provided by Igor Bogush [26], where Σ is a coupled scalar field. The vertical and horizontal axes: β, and α, represent the projection of the BH image on the observer's sky. Panel (b): Result obtained with our approach. Credits for textures used in this image to [29].

Conclusion
The null geodesic trajectories, calculated by exploiting the proposed algorithm shown in Figs. 4, 5, and 6, exhibit an excellent agreement with the exact solutions obtained by solving the orbital equation (2). The main and only limitation of our approach is related to critical and quasi-critical geodesics due to the instability of the photon trajectories owed to the discrete nature of our algorithm. We obtained an angular deviation lower than 3 • for both, the Schwarzschild BH and the Ellis WH under the simulation conditions of this work (see Figs. 4 and 5). For the WH framed in the metric theory of gravity, we obtained considerably larger errors for quasi-critical geodesics. In general, for geodesics with b > b c , the exact and approximated solutions exhibited a reasonable agreement by presenting low angular deviations. This implies that the algorithm may be exploited for the study of emission lines from accretion disks. The accuracy can be improved by either increasing the resolution of the geodesics (but this naturally entails also as a consequence an increment of the computational times) or by implementing alternative strategies for estimating the displacement ∆µ i at every step. For photons falling toward a BH singularity, the approximation of the intersection between the geodesics and the BH event horizon is not critical; instead for photons traversing a WH mouth, the way in which the intersection point between the photon orbit and the WH mouth in the universe 1 is computed might give rise to strong deviations in the final ends of the geodesic segment in the universe 2. The isotropic refractive index attains higher values around the WH mouth, therefore there more precision in the calculations should be employed in order to avoid these problematic issues. The use of isotropic coordinates is extremely important not only for calculation purposes, but they are also able to intuitively capture the main properties of curved space-times. In addition, this algorithm can be straightforwardly generalized to non-zero EM wavelength cases by taking into account the diffraction effects through full-wave numerical solvers based for example on the Finite Element Method (FEM) [30] or the Finite-Difference Time-Domain (FDTD) method [31].
This algorithm can be applied in the numerical image generation of compact objects such as BHs, neutron stars, and WHs. The results presented in this work are very promising if compared with other established and accepted numerical simulations reported in the literature. The advantages and future potentiality of this approach can be shortly summarized in the following points: (1) investigating gravity and gravitational effects through a model easily buildable in laboratory in order to project tests of gravity in small scales, as well as to simulate also and probe gravitational effects occurring around compact objects; (2) vice versa, gravity can be of inspiration for developing new EM technologies, which can be widely employed for engineering applications and to improve the everyday life of the society.