OSIRIS: A New Code for Ray Tracing Around Compact Objects

The radiation observed in quasars and active galactic nuclei is mainly produced by a relativistic plasma orbiting close to the black hole event horizon, where strong gravitational effects are relevant. The observational data of such systems can be compared with theoretical models to infer the black hole and plasma properties. In the comparison process, ray tracing algorithms are essential to computing the trajectories followed by the photons from the source to our telescopes. In this paper, we present OSIRIS: a new stable FORTRAN code capable of efficiently computing null geodesics around compact objects, including general relativistic effects such as gravitational lensing, redshift, and relativistic boosting. The algorithm is based on the Hamiltonian formulation and uses different integration schemes to evolve null geodesics while tracking the error in the Hamiltonian constrain to ensure physical results. We found from an error analysis that the integration schemes are all stable, and the best one maintains an error below $10^{-11}$. Particularly, to test the robustness and ability of the code to evolve geodesics in curved spacetime, we compute the shadow and Einstein rings of a Kerr black hole with different rotation parameters and obtain the image of a thin Keplerian accretion disk around a Schwarzschild black hole. Although OSIRIS is parallelized neither with MPI nor with CUDA, the computation times are of the same order as those reported by other codes with these types of parallel computing platforms.


Introduction
The study of the null geodesics in the curved space-time generated by a black hole or any other compact object serves to analyze the strong gravity effects that are expected near to them as well as to understand the dynamics of the accretion processes. In 2019, the Event Horizon Telescope collaboration published the first images-ever of the shadow produced by the supermassive black hole at the center of the M87 galaxy [8]. M87 * is one of the largest known supermassive black holes and is located at the center of the gargantuan elliptical galaxy, 53 million light-years away. The estimated mass for the black hole is M BH = (6.5 ± 0.7) × 10 9 M and was obtained by comparing the images with an extensive library of ray-traced general-relativistic MHD simulations.
Through numerical simulations, it is possible to process, analyze and compare the observational data. That is one of the reasons why a large number of ray-tracing codes in general relativistic space-times have been developed. From these codes, it has been possible to simulate the motion of photons around compact objects and thus get a better understanding of the gravitational field exerted over them. Solving the equations of the null geodesics analytically allows studying spherical orbits of photons that constitute the wellknown photon ring [2,20,24,38], which conform to the apparent edge of the event horizon and delimit the black hole shadow. However, obtaining an analytical solution for an arbitrary space time is not always possible. Therefore, it is necessary to construct numerical codes that allow finding solutions to the geodesic equations for a large number of photons. Based on this, several numerical studies, related to raytracing, have been carried out [3,15,18,25,26,29,39,43], distinguishing in the literature codes such as GeoKerr [13], Gyoto [40], Ray [34], GRay [7], GeoVis [28], Pyhole [9], Oddisey [35], among others. Each of these stands out for its different attractions, be it for the speed of its algorithms, the ability to solve the geodesic equation in general space-times (stationary or dynamic) in different formalisms, having a user-friendly interface design, or the powerful architecture in which they are written.
Even though nowadays the codes are extremely advanced, including polarized radiative transfer in general relativity [12,19,27,31], in in this work we present and certify our code OSIRIS (Orbits and Shadows In RelativIstic Spacetimes), which in its first version focuses on the solution of null geodesics equations, using the Hamiltonian formulation. OSIRIS uses the Backward ray-tracing method to efficiently calculate the shadow, the Einstein ring of a black hole, and different general relativistic effects, such as gravitational lensing, redshift, and relativistic boosting. Particularly, OSIRIS has been used to study the shadow produced by a naked singularity described by the q-metric, which is the simplest static and axially symmetric solution of Einstein equations with a non-vanishing quadrupole moment [1]. The validation of the code has been tested through different applications like the Black hole shadow in Kerr spacetime and the shadow produced by a thin accretion disk.
This article is organized as follows. In Section 2, we briefly describe the Hamiltonian formulation and the motion equations used to compute the null geodesics. Additionally, we show the backward ray-tracing method to visualize the effects of the strong gravitational field around a compact object. Moreover, to obtain the shadow produced by a black hole, we present the impact parameters, which correspond to all points of the image plane where the observer is located. In Section 3, we describe the numerical scheme to integrate in time. Particularly, we carried out the simulations with four different adaptive time step time integrators and show the norm of the Hamiltonian constraint for each one of them. Besides, as an example, we simulate the shadow produced by the Kerr black hole. In particular, we compare the numerical solution obtained with OSIRIS with the respective analytical solution for each time integrator. Here we show the capability of the code to deal with several initial conditions. Subsequently, to compute the gravitational lensing, it is necessary to classify the orbits to elucidate the effect of gravity over the light. For this reason, in Section 4, we describe this phenomenon and its implementation in OSIRIS. In Section 5, we show the shadow produced by the Schwarzschild black hole surrounded by a Thin accretion disk. Finally, in Section 6, we present our main conclusions and discussions.
Throughout this paper, we employ the signature (−, +, +, +) and geometrized units, in which the gravitational constant and the speed of light are equal to unity. Additionally, we set the mass of the black hole M BH = 1. 2 Shadow formulation

Motion equations
The main goal of OSIRIS is to evolve null geodesics in asymptotically flat and stationary axisymmetric space-times, which implies that there are two constants of motion: the energy E and the azimuthal angular momentum L. In this way, given a metric tensor being g αβ the contravariant components of the tensor and dx α the base of 1-forms, we can obtain the motion equations from the hamiltonian formulatioṅ where the overdot implies differentiation with respect to an affine parameter λ . The Hamiltonian is defined in terms of the contravariant components of the metric tensor, g αβ , and for null particles, it satisfies that where p α are the components of the four-momentum. For simplicity and without loss of generality, we will choose the labels {x µ } = {t, r, θ , φ } for the coordinate system to be spherical-like. It is worth mentioning that OSIRIS also evolves time-like geodesics. In Appendix A, we show trajectories for test particles around a compact object with a nonvanishing quadrupole moment. Thus, from the expressions in 2 and 3, it is possible to obtain the coordinates and momentum of the photons' trajectories as followṡ t = g tt p t + g tφ p φ ,ṙ = g rr p r ,θ = g θ θ p θ ,φ = g tφ p t + g φ φ p φ , p t = 0, p φ = 0, those are the expressions that we solved with OSIRIS.

Backward ray-tracing
To visualize the effects of the strong gravitational field produced by a black hole, it is necessary to analyze the trajectories followed by photons that orbit in its vicinity and translate the information of those that escape and manage to be detected. In essence, this is the main algorithm of OSIRIS, and to understand its operation, we propose two scenarios: Scenario 1: a bright source emits photons in all directions. Among all, it is possible to find those that orbit around the black hole and escape to infinity, others that fall into the event horizon, or even those that do not approach the black hole. The photons that escape and arrive at the observer represent a minimal fraction of all the radiation originally emitted by the source, whereby calculating numerically the motion equations in this scenario translates into a waste of computational time. In other words, this scenario is not optimal and therefore is discarded.

Scenario 2:
The observer is the origin of the null geodesics that evolve back in time. In this way, it is possible to reconstruct the bright source from which they come. This method is known as backward ray tracing, and it is a standard in the simulation of numerical shadows [40]. In figure 1, we illustrate this method, where the gravitational source is a black hole, and the trajectories followed by photons that fall into the event horizon are classified with black color. Whereby, those are physically undetected geodesics whose starting point is the observer.

Impact parameters
To obtain the numerical shadow produced by a black hole, it is necessary to know the information encoded in the photons that orbit in its vicinity. It is then defined as impact parameters to all points (x, y) of the image plane where an observer, located at infinity, makes the measurements [22], as we illustrate in figure 2.
Since the observer is far from the source, it is possible to make a small angles approximation and Fig. 1 Graphic representation of the backward ray tracing: the Earth is the observer from which the photons are emitted, and the black sphere plays the role of a black hole. The lines represent null geodesics, among which are those that are trapped by gravitational attraction (black lines) and those that escape to infinity (colored lines).

Fig. 2
Illustration of the image plane associated with an inertial observer at a distance r 0 from the black hole. The final position of the photons is characterized by the spatial part of four-momentum P measured by the observer (the temporal component is irrelevant since it does not influence the trajectory of the photons), whose origin is taken in the center of the black hole. The impact parameters are related to α and β , where α is the angle between P and its projection in the plane xz, and β corresponds to the angle between this projection and r 0 . Also, the observer is located such that the radial direction coincides with the z-axis, where it is clear thatê where P α correspond to the components of the four-momentum in the observer's frame. These are determined from the canonical momentums through the transformation where η α µ are the components of the metric tensor in the Minkowski space-time, and Λ µ ν are the components of a base change matrix that relates the momentums measured in both reference systems, which is given by Explicitly, the momentums can be written as follows Additionally, the magnitude of the four-momentum satisfies the Hamiltonian constrain Therefore, the expressions in 5 take the form and replacing 9 in 4, we obtain The equations in 10 are general expressions for the impact parameters in an axisymmetric background, which only depends on the metric tensor components [11]. Thus, setting set P t = E = 1 without loss of generality, and combining equations 10 and 11, initial momentums can be computed in terms of each point on the image plane as follows This image plane setup has been used to numerically calculate shadows generated by compact objects in various spacetimes, as well as to study the impact of the intense gavitational field on phenomena such as the so-called gravitational lensing [9-11, 39, 41, 42]. Furthermore, thanks to this configuration, in some particular cases in which the system of equations is variable-separable, it is possible to find analytical solutions to determine the rim of the shadow through the calculation of the spherical orbits of photons around certain compact bodies. Next, we show an example in Kerr spacetime.  [5]. The Runge-Kutta methods are the commonly staged methods used for the approximation of solutions of ordinary differential equations. On the other hand, the BS algorithm combines three different ideas: Richardson extrapolation that considers the final answer of a numerical calculation as an analytic function of an adjustable parameter like the stepsize. Once there is enough information about the function, it is fit to some analytic form through a rational function extrapolation, to finally carry out the integration by the modified midpoint method [32].
Preserving the Hamiltonian constraint for a photon that orbits around a Kerr black hole with dimensionless spin parameter a = 0.98 and initial conditions (t 0 , r 0 , θ 0 , φ 0 ), we considered two scenarios: In figures 4 and 5 we plot the hamiltonian norm for two particular orbits: an escape one and a fall one, respectively. Our numerical results suggest that in both cases, the RK Dormand-Prince preserves the hamiltonian norm better than one part in 10 −10 , being the best option, with the BS algorithm the second one. It is worth mentioning that Dormand-Prince method has seven-time steps but only uses six func-tion evaluations per step. This method chooses the coefficients to minimize the error of the fifth-order solution. It is the main difference with the Fehlberg Runge-Kutta time integrator, which was constructed so that the fourth-order solution has a small error. It is clear that the error grows as the photon approaches the event horizon at r 0 = 100. The oscillations in the curves with each integrator occur due to the step refinement in regions where the gravitational field is stronger, i.e., near to the event horizon.
where the points (x, y) are parameterized by the radii r of the spherical orbits through ξ and η, which are given by the expressions In figure (6), we show the analytical and numerical shadow of a black hole with a spin parameter a = 0.98. Comparing both results, we obtained that the simulation is a little smaller than the analytical. It is because this solution is calculated at infinity, but it is not possible to simulate an observer with these characteristics; in addition, initializing the motion too much larger radii means that the error accumulates more and more between iterations. On the other hand, to carry out a global analysis, we plot the average of the hamiltonian norm in a square scattering region bounded by −8 ≤ x ≤ 8 and −8 ≤ y ≤ 8, with a dense uniform grid of 625 × 625 initial conditions. Using a logarithmic color palette, we displayed the final value of the hamiltonian norm for the entire phase space ( Figure 7). As we have seen before, the RK Dormand-Prince preserves the constraint better than a part in 10 −11 (in the worst scenario e.g., near the black hole, better than a part in 10 −10 ). In contrast, the remaining methods "fail" to preserve the constraint, yielding, in the best case, to values of the order of 10 −3 for the RK methods, and of the order of 10 −8 in the case of the BS algorithm. Finally, Figure 8 exposes simulation time machine as a function of the mesh resolution. We can see that the differences between the times required by each integrator increase with the resolution. Thus, at higher resolution, the RK Dormand-Prince and the BS algorithm require more computational time than the RK Cash-Karp and the RK-Fehlberg. However, due to the better numerical behavior of the former when maintaining the norm, and the almost negligible increase in the time machine, the RK-Dormand-Prince emerges as  the best option for our purpose. Based on the above, we compare our result with another codes. In [7], it is possible to appreciate a time comparison between three different codes: GRay [7], GeoKerr [13] and Ray [34], obtaining that in the process of integration of 10 6 geodesics the time orders in seconds employed for each code respectively were 10 3 , 10 4 and 10 5 , making of GRay the fastest of the three. OSIRIS's algorithm evolves 1024 2 geodesics in the same order of time as GRay, see Figure 8. Furthermore, we show that OSIRIS, besides being fast, preserves both the the hamiltonian constrain and the hamiltonian norm below 10 −11 . Nevertheless, comparing our results with Odyssey [35], it evolves 1024 2 geodesics in two orders of magnitude of time less than OSIRIS. It is necessary to mention that both GRay and Odyssey employ a powerful programming architecture based in GPU usage and parallelization in CUDA, in which each code evolves millions of geodesics at the same time. With the error analysis discussed in this section, we concluded that the calibration of the code is optimal. Then, we proceed to study the influence of the gravitational field on the trajectory of the photons.

Celestial sphere
Now we are going to classify the orbits to elucidate the effect of gravity on null geodesics. For this purpose, following [3], we define a celestial sphere as a bright source from which the light rays will be emitted. This sphere is concentric with the black hole, surrounds the observer, and is divided in four colors: blue, green, yellow, and red. Furthermore, the notion of curvature is provided by a black mesh with constant latitude and longitude lines, separated by 6 • ( figure 9). In this way, the classification will be as follows: if the photon's radial position, r, satisfies the condition r < r H + δ r , where r H is the event horizon radius and δ r is a little buffer, then it corresponds to a black point in the image plane. On the contrary, if the photon escapes from the strong gravitational attraction and r > r cs , where r cs is the radius of the celestial sphere, the integration process is stopped and a color is assigned depending on the sector of the sphere where the photon strikes. Fig. 9 Scheme of the celestial sphere where we remove a section of this to see inside the observer (Earth) and the compact object. The sphere is divided into four quadrants as follows. Top, for 0 ≤ θ < π/2: green quadrant, if 0 ≤ φ < π; red quadrant, π ≤ φ < 2π. Bottom, for π/2 ≤ θ < π: blue quadrant, if 0 ≤ φ < π; yellow quadrant, if π ≤ φ < 2π. The white line represents the direction of observation that coincide with radial direction.

Gravitational lensing
Due to the deflection of light in the presence of a compact object, the optical perception that we could obtain from a bright source can provide interesting results at the time of observation. This phenomenon is known as gravitational lensing. To have a better idea about it, first of all, we will consider the Minkowski space-time, where the rays travel in a straight line ( figure 10). Whereby, we obtain four quadrants that correspond to the real image of the celestial sphere. On the other hand, when a compact object interposes between the observer and the source, the strong gravitational field that it generates deflects the path of the photons. In figure 11, we show a pictorial representation of it, where a region of the celestial sphere may appear to have been emitted from a different area due to the curved trajectories of photons.

Gravitational lensing in Kerr space−time
As the first application of our code, we present numerical simulations of the gravitational lens produced by a Kerr black hole for different values of the dimensionless rotation parameter a, seen from the equatorial plane in contrast with flat space-time. The line element that describes this geometry, in Boyer-Lindquist coordinates {t, r, θ , φ }, is given by Fig. 10 Graphic representation of the trajectory of photons in the Minkowski space-time. The observer is looking directly at the celestial sphere, and in the absence of any gravitational field source, the path of the photons is a straight line. Therefore, in the observer's plane, we obtain an image that corresponds with a section of the celestial sphere. Fig. 11 Graphic representation of the deflection in the path of the photons in the presence of a compact object. Due to the strong gravitational field, some photons emitted from a quadrant can be detected in the image plane as if they had a different point of origin.
where a = J/M relates the angular momentum of rotation J and the mass M of the gravitational source, and the expression for the outer event horizon is given by In the upper left row of figure 12, we show the gravitational lensing produced by the Minkowski space-time, where the black lines of the mesh are not straight due to the curvature of the sphere itself. The upper right quadrant is the non-rotating, a = 0, corresponding to a Schwarzschild black hole. An appreciable effect of the gravitational lens is the Einstein ring, a phenomenon that appears when the observer, the compact object, and the source are aligned. Due, the light source looks like a concentric ring around the black hole. Inside this ring, the deflection angle of the photons's trajectories is big enough to generate an inversion of the colors [3]. Additionally, a second ring can be seen near the edge of the shadow, where the images are inverted again. The gravitational lens produced by a black Kerr hole with rotation a = 0.5 and a = 0.98 is shown at the bottom row of figure  12, where the Einstein ring appears as in the Schwarzschild case. However, as a result of the gravitational drag due to the black hole rotation, an asymmetry occurs in the lens. In the lower-left image, the blue and green quadrants spread along the inner edge of the ring, which generates an "ear" in yellow and red quadrants. This ear extends along with the shadow's silhouette. These effects are more noticeable as the rotation increases. In the lower-right quadrant of figure 12, the ear spreads over the edge of the shadow and envelops it. Additionally, a consecutive succession of partial Einstein rings appears in the vicinity of the flat edge of the shadow as a consequence of extreme rotation. Furthermore, a characteristic in common is the deformation suffered by the lines of the time-space mesh, providing the notion of curvature due to the strong gravitational field.

Thin accretion disk
Accretion disks are a great mechanism to determine the properties of a black hole, like its mass or spin [21,30]. Furthermore, the emission spectrum and the radiative flux from these structures can depend on space-time generated by the compact object [37], emitting from the radio band to x-rays [23,33]. The numerical simulations are useful to compare and interpret the observational results with different theoretical models. Next, we consider a thin accretion disk around a Kerr black hole in the equatorial plane, assuming an ideal non-self-gravitating disk moving in circular orbits (u r = u θ = 0, with u α the components of the 4-velocity for time-like particles moving on the disk). Let Ω the angular velocity and l 0 the specific angular momentum of time-like particles on the disk, then Based on our ray-tracing, the integration process for null geodesics stops when photons reach the disk, then an observed intensity is assigned to each point (x, y) on the image plane. Through the Lorentz invariant I/ν 3 , the observed intensity I obs , can be expressed in terms of the emitted intensity I obs and the respective frequencies as follows with the redshift factor due to the gravitational spectral shift and the Doppler effect (recalling that calligraph letters means physical magnitudes measured by the observer). Now, employing the normalization condition the equation 20 can be written as As result of our numerical simulation, in the figure 13, we show the image of a Kerr black hole surrounded by a thin accretion disk for different dimensionless spin parameter (a = 0, 0.5, 0.95) with a map for the observed intensity using the disk model proposed in [30], where it is evident that high rotation increases the intensity measured by a distant observer. We can appreciate a bright spot located at the left of the disk as an effect of the redshift, which decreases its size as rotation increases, concentrating the maximum observed intensity in a small region. On the other hand, the bend of light due to the extreme gravitational field allows us to see the back of the disk, which should be hidden by the black hole (a deeper discussion and more details are carried out by Luminet [25]). Additionally, when considering the disk as a bright source, it is possible to appreciate how the edge of the shadow is delimited by a thin layer of light coming from the disk.

Conlusions
In this paper, we present OSIRIS, new stable ray-tracing FORTRAN code capable of efficiently compute null-geodesics in stationary and axisymmetric space−times. OSIRIS incorporates general expressions for the impact parameters relating to the measurements of an observer located in the vicinity of the gravitational source and another at infinity. As a first application, we simulate the image of a rotating black hole seen by a distant observer and study (qualitatively) the effect of the gravitational dragging in the shadow and gravitational lensing. The image of a compact object is obtained by solving the motion equations for photons in each point of the observer's screen. We implement four different stable integration schemes with adaptive step: Runge-Kutta Fehlberg 45, Runge-Kutta Cash-Karp 45, Runge-Kutta Dormand-Prince 45, and a Bulirsch-Stoer algorithm. We analyze the error in the Hamiltonian constraint for null particles (H = 0) and find that the RK Dormand-Prince preserves the constraint better than a part in 10 −11 , while the other schemes produce considerably larger errors. Although this method requires more computational time than the RK Fehlberg and RK Cash-Karp, especially when dealing with high-resolution simulations, the time differences between the schemes are negligible for reasonable resolutions. Therefore, the RK Dormand-Prince emerges as the best option for the main goal of OSIRIS (photons dynamics around compact objects). Additionally, it should be noted that even though OSIRIS is parallelized neither with MPI nor with CUDA, the computation times are similar to those obtained by other codes programmed with these parallel computing platforms.
OSIRIS can simulate accretion disks around black holes, from which it is clear that the apparent form of the black hole and the observed intensity depend on the space-time that is considered. For this reason, develop theoretical models for accretion disk surrounding black holes is of great importance when comparing observational results with the computational simulations, since from these it is possible to obtain information about the space-time itself. Fig. 13 Simulation of a thin accretion disk for a black hole with different spin parameter (a = 0, 0.5, 0.998). The observed intensity show a bright dot at the left of the disk produced by the redshift, whose location is due to the direction of rotation of the disk. It is clear that the observed intensity increases as the spin parameter increases, being an order of magnitude higher for an extreme rotation black hole (a = 0.95) compared to the non-rotating case (a = 0) and the intermediate rotation case (a = 0.5). The simulation parameters are: observer location at r 0 = 1000, θ 0 = 85 • , φ 0 = 0 • and t 0 = 0; the limits of the disk are r in = r isco (for every case of rotation) and r out = 20. The specific angular momentum, l 0 , is 2. corresponds to naked singularities without event horizons [4]. Besides, for q = 0, this metric reduces to Schwarzschild space-time. In figure 14, we show time-like geodesics in the equatorial plane for different values of q. In the first row, we appreciate how the influence of q transforms unbounded Schwarzschild trajectories into bounded ones for particles with non-initial radial velocity (ṙ = 0). In the second row, we show a Schwarzschild bounded orbit and the geodesics withṙ = 0 around a prolate (q < 0) and oblate (q > 0) source. We can see that the quadrupole only affects the structure of the trajectory. Finally, in the third row, we present particles withṙ = 0. In this case, all of them scape toward infinity, and the parameter q only defines the scape direction.