Testing the validity of the raytracing code GYOTO
LESIA, Observatoire de Paris, CNRS, Université Pierre Marie Curie,
Université Paris Diderot,
5 place Jules Janssen,
92190
Meudon,
France
email:
marion.grould@obspm.fr
Received: 17 June 2015
Accepted: 22 April 2016
Context. In the next few years, the nearinfrared interferometer GRAVITY will be able to observe the Galactic center. Astrometric data will be obtained with an anticipated accuracy of 10 μas. To analyze these future data, we have developed a code called GYOTO to compute orbits and raytrace images.
Aims. We want to assess the validity and accuracy of GYOTO in a variety of contexts, in particular for stellar astrometry in the Galactic center. Furthermore, we want to tackle and complete a study made on the astrometric displacements that are due to lensing effects of a star of the central parsec with GYOTO.
Methods. We first validate GYOTO in the weakdeflection limit (WDL) by studying primary caustics and primary critical curves obtained for a Kerr black hole. We compare GYOTO results to available analytical approximations and estimate GYOTO errors using an intrinsic estimator. In the strongdeflection limit (SDL), we choose to compare null geodesics computed by GYOTO and the raytracing code named Geokerr. Finally, we use GYOTO to estimate the apparent astrometric displacements of a star for different angles from Sagittarius A* (Sgr A*).
Results. In the WDL, we find a good coherence between GYOTO results and approximations. The maximal difference is around 10^{5}μas. Our intrinsic estimator finds a conservative uncertainty estimate also around 10^{5}μas. In the SDL, both raytracing codes find the same photon’s coordinates with a maximal difference of about 10^{3}μas. The shift of a star located behind the plane of sky containing Sgr A* is consistent with the current study. In addition, the effect of lensing on any star in this plane of sky is a radial shift by 5 μas, independent of the distance from Sgr A* up to a very large distance.
Conclusions. We have demonstrated that GYOTO is accurate to a very high level, orders of magnitude better than the GRAVITY requirements. GYOTO is also valid in weak and strongdeflection regimes and for very long integrations. At the astrometric precision that GRAVITY is aiming for, lensing effects must always be taken into account when fitting stellar orbits in the central parsec of the Galaxy.
Key words: black hole physics / gravitational lensing: weak / gravitational lensing: strong / Galaxy: center
© ESO, 2016
1. Introduction
General relativitY OrbiT of Observatoire de Paris (GYOTO^{1}) is a raytracing code developed by Vincent et al. (2011). It integrates null and timelike geodesics in any analytical or numerical metrics. GYOTO can compute images and spectra for a variety of astrophysical objects, such as moving stars or accretion disks, around a Kerr black hole. The code’s ability to take numerical metrics into account allows it to compute images or trajectories of stars orbiting exotic objects such as a boson star (Grandclément et al. 2014; Vincent et al. 2016).
The main motivation for the development of GYOTO was to interpret the data to be obtained with the secondgeneration very long baseline telescope interferometry (VLTI) instrument GRAVITY (Eisenhauer et al. 2011). This instrument observes stars and flares orbiting Sgr A*. It probes spacetime near the central object with an expected astrometric accuracy of 10 μas. In a preliminary work, Vincent et al. (2014) used GYOTO to show that GRAVITY is capable of distinguishing an ejected blob from alternative flare models.
The stellar orbits measured by GRAVITY is affected by several effects such as periastron shift and LenseThirring effects (Will 2008; Merritt et al. 2010). In addition, the individual astrometric measurements are affected by relativistic effects: time delay and lensing (Bozza & Mancini 2012).
All these effects need to be considered in an apparent orbit model that will be fitted to the GRAVITY data, enabling the nature of Sgr A* to be constrained. Since the goal of GRAVITY is to deliver astrometry at an accuracy of 10 μas, models need to be more accurate than this value so as not to limit the accuracy of final results. We therefore aim for a model with an astrometric accuracy of 1 μas. In this paper, we study the accuracy of GYOTO to determine whether this tool can be used as a foundation for a future apparent orbit model to fit the GRAVITY data. Using the star images computed by GYOTO, it will be possible to get the apparent position of the star. However, the accuracy of this position will depend on the precision of the photon trajectories. Null geodesics need to be properly computed by the integrator implemented in GYOTO to take into account the correct bending effect. Besides, because of the 2′′ fieldofview of GRAVITY, a wide range of distances between stars and Sgr A* will be possible. GYOTO has never been used with this type of a configuration, so we need to ensure that geodesics are well computed.
We first focus on the Einstein ring formation, more precisely we study primary caustics and primary critical curves considering a Kerr black hole. The aims are both to compare our numerical results with the analytical study on primary caustics and primary critical curves performed by Sereno and De Luca in 2008 (Sereno & de Luca 2008), and to check whether the numerical error is sufficiently low, which means inferior or equal to 1 μas. The comparison between GYOTO and approximations is a validation of GYOTO in the weakdeflection limit (WDL), however we also have to check if this raytracing code is valid in the strongdeflection limit (SDL). To do so, we choose to compare null geodesics computed in GYOTO and with another code named Geokerr^{2} (Dexter & Agol 2009). Before starting the validation and test of GYOTO, we discuss the different integrators implemented in it to choose the most appropriate one. To do so, we consider a Schwarzschild black hole and compare GYOTO with one of the approximations developed by Sereno & de Luca (2008).
In this paper, we also investigate lensing effects on the apparent position of a star to compare the astrometric displacements found in Bozza & Mancini (2012). We demonstrate that the effect cannot be neglected even when the star is in front of the plane of the sky containing Sgr A*. We also discuss the impact of neglecting the lensing effects on the estimation of the orbital parameters and on the detection of other effects affecting the orbits (see Sect. 5).
The paper is organized as follows. In Sect. 2, we introduce our notations and define primary caustics and primary critical curves in the Schwarzschild and Kerr cases. We present the different analytical approximations we choose for our study. In Sect. 3, we discuss the validity of the different integrators implemented in GYOTO. Section 4 is devoted to the comparison of the GYOTO results with those obtained by other methods: analytical formulae of primary caustics and primary critical curves (see Sect. 4.1) and Geokerr (4.2). Methods implemented to determine the different parameters with GYOTO are explained in Appendix A. A discussion on lensing effects is presented in Sect. 5. Concluding remarks are given in Sect. 6.
2. The Einstein ring formation: focus on primary caustics and primary critical curves
Fig. 1 Lensing configuration: O, S, and L denote the observer, the source, and the lens, respectively. 

Open with DEXTER 
To understand how the Einstein ring is formed, we recall the basics of gravitational lensing using a Schwarzschild lens. Then, we focus on the Einstein ring obtained with a Kerr black hole. In both cases we consider a static observer. The geometry for lensing configuration is shown in Fig. 1. The spin axis coincides with the zaxis. Spherical coordinates of the observer and the source, relative to the lens L, are noted (r_{0},ϑ_{0},φ_{0}) and (r_{s},ϑ_{s},φ_{s}), respectively. Without loss of generality, we choose to work in the equatorial plane of the black hole: ϑ_{0} = π/ 2, φ_{0} = 0 and ϑ_{s} = π/ 2. This yields (x_{0},0,0) for the observer and (x_{s},y_{s},0) for the source. We note M the lens mass and a the spin of the black hole ranging from 0 (Schwarzschild black hole) to 1 (extremal Kerr black hole). From the point of view of the observer, the black hole rotates from the left to the right. In this paper, we use two different units for the distance: parsecs and geometrical units. This last unit is equal to GM/c^{2} with G the Newton’s constant and c the speed of light, but we will consider G = c = 1 and note it M.
2.1. Schwarzschild lens
Fig. 2 Spatial projection of a Schwarzschild lensing situation: S corresponds to the source, L to the lens and O to the observer. 

Open with DEXTER 
Using the notation of Fig. 2, we can write the lens equation (Schneider et al. 1992) (1)with β the unlensed angular position of the source, θ the lensed angular position of the source equal to ξ/r_{0} and the deflection angle depending on the impact parameter ξ. The latter angle is given by (2)with R_{S} = 2GM/c^{2} the Schwarzschild radius. We can rewrite the lens equation as (3)with (4)which corresponds to the Einstein angle. The magnification of the source in the lens plane is a function of the lensed and unlensed angles as (5)A is infinite when det(∂β/∂θ) = 0. In the source plane, these positions are called caustic points. For a Schwarzschild lens, the caustic is a line behind the lens starting from it and extending toward infinity (Rauch & Blandford 1994). If the source lies on the caustic line then β = 0. Thus, the solution of the lens equation is θ = α_{0}. A circle called a critical curve is formed in the lens plane with a radius of α_{0}. Considering the source as a star, the observer sees the wellknown Einstein ring. The radius of the ring corresponds to the critical curve radius so we get α_{0} = θ_{E} with θ_{E} the Einstein ring radius. If the star does not lie on the caustic, the observer will see two images called primary and secondary images. These images are formed by lensing effects. Light rays are deviated because of the curvature of spacetime by the black hole. At the caustic points, the lensed images merge into the Einstein ring.
Fig. 3 GYOTO image of a fixed star, with a radius of 1M, located at 6M in front of a Schwarzschild black hole. The observer is at 8 kpc from the black hole. The star lies on the secondorder caustic line: the biggest image corresponds to the primary image of the star and the ring to the critical curve formed thanks to the secondorder caustic. Axes are labeled in μas and the field of view of the image is equal to 70 μas. 

Open with DEXTER 
Several studies have been performed on caustics and critical curves with Schwarzschild black holes. Rauch & Blandford (1994) showed the existence of several orders of caustics. The caustic line we previously discussed is the firstorder (also known as primary) caustic. This caustic is generated by light rays, which do not wind around the black hole. Higher orders correspond to photons winding one or several times around the black hole. Bozza (2008) showed that these caustics are also lines but caustics of even order start from the black hole and extend to the observer. It is also possible to form an Einstein ring if the star lies on these caustics of even order. But it is harder to detect as the largest fraction of the flux is concentrated in the primary image of the star (see Fig. 3).
2.2. Kerr black hole lens
In case of the Kerr black hole, there is also a primary caustic and higherorder caustics but these are no longer lines (Rauch & Blandford 1994; Bozza 2008). Rauch & Blandford (1994) were the first to discover that the primary caustic is a tube with an astroid (fourcusped) crosssection (see the upper scheme in Fig. 4). At large distances the crosssection is symmetric but becomes distorted near the horizon. In addition, the closer to the BH the source is, the larger the tube shifts with respect to the Schwarzschild’s case. For sources near the horizon, the caustic winds around the black hole in the opposite sense with respect to its rotation. Very far from the black hole, the shift is still significant but the size of the caustic (distance between the right and the left cusp of the astroid crosssection) decreases and tends towards zero.
To form the critical curve the source must cover all of the astroid crosssections. If we consider a point source in the caustic, four images are formed in the observer’s sky. If the point source is on the caustic surface or outside the caustic, only two images are formed. An illustration is given in Fig. 4.
Fig. 4 Effect on the displacement of a point source, in the equatorial plane of the black hole, with respect to the caustic (top) on the formation of images in the observer plane (bottom). The black astroid represents the primary caustic and the black circle represents the primary critical curve. We define the z′ and y′ axis whose origin corresponds to the middle of the astroid crosssection. The α and δ axis correspond to the observer’s sky axis, the origin is centered in the middle of the screen (because of the rotation of the black hole, the critical curve is shifted from the center of the screen). 

Open with DEXTER 
Because of the rotation of the black hole, if the source lies on the right cusp one of the two images will be formed on the left side of the critical curve.
Analytical studies of primary caustics and primary critical curves has been made (Sereno & de Luca 2006, 2008). All of their analytical approximations were obtained in both weakdeflection and weakfield regimes. In the WDL, photons do not wind around the black hole which means that the minimum distance between the photon and the lens r_{minγ} must satisfy R_{S} ≪ r_{minγ}. In the weakfield regime, we verify the condition r_{s} ≫ R_{S}. In both regimes, the primary caustic is only shifted and keeps a symmetric shape. Because of the shift of the caustic, the critical curve is not centered on the black hole. In Sereno & de Luca (2008), three equations are developed through a Taylor expansion of the null geodesics in (6)where D = r_{s}/ (r_{s} + r_{0}) and θ_{E} is the Einstein ring radius. All of the equations in this paper are expressed in the equatorial plane of the black hole. One of the three equations gives the radius of the critical curve: (7)The left and the right radius of the critical curve are equal and depend on the spin in the thirdorder term in ε. The second equation gives the offset Δ of the center of the critical curve relative to the black hole: (8)Because the black hole is spinning from the left to the right, the critical curve is shifted to the right (y> 0). The last equation gives the position y_{C} of the primary caustic for a given x_{s} and z_{s}: (9)with (10)where Δ_{C} is the size of the caustic. The angle η ranges from 0 to 2π. In this study, we decide to focus on the right cusp of the primary caustic and to determine the angular position B_{C} of this cusp seen from the Earth. It means that we assume η = ± π and z_{s} = 0. This yields, for a given x_{s}(11)Analytical approximations in the SDL have also been derived (Bozza et al. 2005; Bozza & Scarpetta 2007) but not for the primary caustic since it is formed in the WDL. A numerical study was made by Bozza (2008) who study the full structure of caustics and critical curves and compared the results with available analytical formulas. In the case of primary caustics, the author compared the size and the position of the left cusp of the caustic and showed that the approximations of Sereno and De Luca only fail at very small distances from the horizon.
To validate GYOTO in the WDL, we estimate the three parameters Θ_{E}, Δ, and B_{C} with GYOTO and make a comparison with formulae (7), (8), and (11), respectively. To reproduce the observational conditions of GRAVITY, we consider an observer at r_{0} = 8 kpc from a black hole of mass M equal to 4.31 × 10^{6}M_{⊙}. We also consider a source far enough from the black hole to be compliant with the domain of validity of these approximations. In the next two sections, we estimate the three parameters as a function of the position x_{s} of the source. The methods used to measure these parameters using GYOTO are given in Appendix A. For each distance of the source, we estimate the error made on the different parameters. Since the goal of this paper is to determine whether the accuracy of GYOTO is better than 1 μas, we only consider the maximum error for each parameter (see for instance the end of Appendix A.1.1 for the maximal error estimation of the angular position of the caustic point). The next section is devoted to discussing the different integrators implemented in GYOTO and used to integrate null geodesics.
Fig. 5 Offset estimated with GYOTO in μas versus the distance  x_{s}  of the point source in parsec, in the Schwarzschild case. Each plot represents different values of the tuning parameters h and AbsTol: 10^{18} (left), 10^{16} (middle), and 10^{14} (right). The different types of curves denote the different integrators: Legacy Generic Integrator (solid), Runge Kutta Cash Karp 54 (dash), and Runge Kutta Fehlberg 78 (dot). 

Open with DEXTER 
3. Comparison of various integrators of GYOTO
Fig. 6 Comparison of the three integrators by measuring the offset of the Einstein ring in the Schwarzschild case. Left: computing time versus the tuning parameters. Middle: offset versus the tuning parameters. Right: computing time versus the offset. The different types of curves denote the different integrators: Legacy Generic Integrator (solid), Runge Kutta Cash Karp 54 (dash), and Runge Kutta Fehlberg 78 (dot). 

Open with DEXTER 
To compute null geodesics in GYOTO, different integrators were implemented to solve the equations of motion of the photon. Before estimating the three parameters, we need to do an intrinsic test of GYOTO to know which integrators are best. To test their validity, we estimate the offset of the critical curve in the Schwarzschild case, which means that Δ should tend to zero. In this section, we do not estimate the position of the caustic point, we directly put the source at (x_{s}, 0, 0). For more details on the estimation of the offset see Appendix A.2. Six different integrators are available in GYOTO, however we choose to compare the following three:

the Legacy Generic integrator. This is the first integrator that has been implemented in GYOTO. It is a fourthorder Runge Kutta integration controlled by several parameters. An important tuning parameter for this integrator, DeltaMaxOverR (h), depends on the maximal integration step (δ_{max}) and the current distance to the center of the coordinate system (r) as h = δ_{max}/r. Thus, the integration step cannot be larger than h × r.

the Runge Kutta Cash Karp 54 integrator. This is an integrator of the Odeint library (for numerically solving ordinary differential equations) belonging to the Boost C++ Libraries^{3}. Several Boost integrators have recently been added to GYOTO because they are well debugged thanks to their wide user base and enable the user to make runtime choices depending on the problem at hand. Several tuning parameters are used to optimize the integration. For instance, the integrator controls the error on the step integration by computing two steps with different orders, the fourth and the fifth. It estimates the difference between these two steps to evaluate the error. This error is compared against an error tolerance depending on two parameters. They are called AbsTol and RelTol and represent the absolute and the relative tolerance, respectively. RelTol is used for high values of solutions and AbsTol for solutions close to zero.

the Runge Kutta Fehlberg 78 integrator. This is another integrator from the Boost family. It computes the two steps with higher orders, the seventh and the eighth.
The offset is obtained taking into consideration these three integrators and different values for the tuning parameters h, AbsTol, and RelTol. We choose to work with RelTol=AbsTol but this is not a unique choice. We did not investigate independent variations of the two parameters because this is not relevant to our study. In what follows, the two parameters are set equal and we only use the name AbsTol.
We also note that there is a systematic error in GYOTO generated by the machine accuracy. More precisely, this bias corresponds to the position of the observer, which is not exactly at y_{s} = 0 but is shifted by approximately y_{s} = 10^{6}M owing to the choice of coordinate system. This shift leads to a (negligible) systematic effect around 10^{5}μas in the observer plane.
The offset of the critical curve for h^{3}/ 10^{8} and AbsTol equal to 10^{18}, 10^{16}, and 10^{14}, versus the absolute distance of the source  x_{s}  in parsec (the source is always behind the black hole thus x_{s}< 0) is presented in Fig. 5. We consider varying h^{3}/ 10^{8} and not h for the Legacy integrator because the tuning parameters AbsTol and h are not defined in the same way. The quantity h^{3}/ 10^{8} is obtained empirically to globally recover the same error behavior between Legacy and Boost integrators. In Fig. 5, the noise generated by GYOTO increases with the distance of the source to the black hole for all integrators. Boost integrators are more accurate than Legacy for tuning parameters equal to 10^{18}. However, the offset of the Legacy integrator remains small (< 10^{2}μas), even for large distances and for the three values of h. The offset obtained with the Runge Kutta Fehlberg 78 integrator significantly increases with the distance for AbsTol = 10^{16} and 10^{14}, and can exceeds 1 μas for very large  x_{s} . The other Boost integrator Runge Kutta Cash Karp 54 is accurate (Δ ≈ 10^{5}μas) for all the distances investigated and for the three values of AbsTol.
When the star is located between  x_{s}  ≈ 10^{4} parsec and  x_{s}  ≈ 10^{1} parsec, the offsets Δ, which are evaluated for each integrator and each value of tuning parameter, fluctuate around 10^{5}μas. These errors are dominated by the observer position shift that is due to machine accuracy. However, a second error appears for larger distances of the star and is due to the interpolation made in GYOTO. More precisely, this interpolation is made in a function called MinDistance in GYOTO. This function is used to evaluate the right and left Einstein radii needed to estimate the offset Δ (see Appendix A for the definition of the MinDistance function and an explanation on the method implemented to obtain Δ). This limitation is due to the step size in which the interpolation is made to estimate the MinDistance function. Indeed, when the step size increases, the interpolation is less efficient and the offset Δ increases. This second limitation is only apparent for the Legacy and Runge Kutta Fehlberg78 integrators because these two integrators generate less intermediate points in this context, and therefore incur more interpolation. This source of error is therefore not cumulative (it is reset at each point estimated by the integrator). We have looked at the statistics for this error, which are essentially Gaussian.
Fig. 7 Absolute difference between analytical approximations (7), (8), and (11), and GYOTO measurements for the three parameters B_{C}, Θ_{E}, and Δ. The types of line denote different values of the spin: 0.2 in solid, 0.5 in dotted, and 0.9 in dashed. 

Open with DEXTER 
All integrators seem to meet the requirements for the studies to follow. However, we also decided to evaluate the computing time needed by the various integrators to reach their best astrometric accuracy. To do so, we estimate the offset of the Einstein ring for different values of the tuning parameters (ranging from 10^{10} to 10^{18}), and the computing time needed by each integrator to evaluate the offsets. We used a 2.66 GHz Intel Core 2 Duo processor to compute null geodesics. We also consider a star at  x_{s}  = 2 parsecs.
Figure 6 shows that Runge Kutta Fehlberg 78 is the fastest integrator for all values of the tuning parameters. Its best accuracy (≈ 10^{5}μas) is reached for AbsTol in the range 10^{17}−10^{18}, corresponding to a computing time in the range 15−35 s. The Runge Kutta Cash Karp 54 integrator is always accurate (≈ 10^{5}μas), its computing time is between 15 s (AbsTol= 10^{10}) and 45 s (AbsTol= 10^{18}). The best accuracy for the last integrator (≈ 10^{5}−10^{4}μas) is reached for h^{3}/ 10^{8} between 10^{18} to 10^{12}, which means a computing time between 10 min to 10 s, respectively. The optimal compromise between precision and computing time is thus reached at AbsTol= 10^{17} for the Runge Kutta Fehlberg 78 integrator, AbsTol= 10^{10} for Runge Kutta Cash Karp54 and h^{3}/ 10^{8} = 10^{12} for Legacy, where the offset is around 10^{5}μas and the time needed is less than one minute.
The important time difference between Boost integrators and the Legacy integrator at AbsTol= 10^{18} is due to the method used to estimate the adaptive step during the integration. Indeed, considering the number of steps used to estimate one null geodesic, ≈90 000 steps are needed for the Legacy integrator, however, Boost integrators only need one order of magnitude less than Legacy. The accuracy limitation on Δ on this figure is also dominated by the two types of error discussed above.
This quick study shows that the three integrators are appropriate for the following studies. In addition to the optimal compromise mentioned above, other reasonable choices of numerical parameters can be made such as AbsTol= 10^{18} with Runge Kutta Fehlberg 78, where we still have an accuracy of around 10^{5}μas and a computingtime of less than one minute. Since the following studies do not use timeconsuming methods, we choose to consider AbsTol= 10^{18} and Runge Kutta Fehlberg 78. This is the choice that we made for the rest of this paper.
4. Results
4.1. Validation of GYOTO in WDL
In Fig. 7, we present the absolute difference between analytical approximations and GYOTO measurements obtained for three values of the spin, 0.2, 0.5, and 0.9. For the range of parameters investigated, the differences are always extremely small (≲10^{2}μas).
Fig. 8 Absolute difference between analytical approximations and GYOTO measurements considering three different equations for the approximation Θ_{E}: the solid line is obtained considering only the zeroth and first orders in ε in the Eq. (7), the dashed line considers the zeroth, first, and second orders in ε and the dotted line considers the full equation. Only a = 0.5 is plotted. 

Open with DEXTER 
On each plot two different regimes can be observed. For small  x_{s} , the curve is marked by a smooth, powerlaw decrease: GYOTO and the numerical approximation agrees better and better for larger and larger values of  x_{s} . After reaching a minimum, the curve raises again, with a much more noisy appearance. This is due to GYOTO being better than analytical approximations for small  x_{s} . The difference between the two traces the order of the approximation. On the other hand, for large values of  x_{s} , the approximations win over GYOTO and the difference is dominated by the numerical error of GYOTO. To demonstrate this interpretation, we compare GYOTO to three approximations of different order in Fig. 8. The first absolute difference is obtained considering the zero and the first orders in ε (ε^{0} and ε) in Eq. (7). For the second absolute difference, we consider one more order in ε (ε^{0}, ε and ε^{2}). The last one is obtained considering the full Eq. (7), which is the same plot shown in the middle of Fig. 7. We can see that, when smaller orders are considered, the differences decrease for all distances of the source. We also note that the noisy part appears earlier and earlier when considering accurate approximations. We thus confirm that the first part of the curves presented in Fig. 7 is due to analytical approximations.
The maximal errors on each parameter evaluated with GYOTO, and for each spin, are presented in Table 1. The method to estimate this error is presented in Appendix A.2. All the maximal errors in Table 1 are around 10^{5}μas. The limitation in the estimation of these errors is due to the systematic bias induced in the observer screen. The raytracing code is very accurate, even for sources far behind the black hole (e.g. δ_{ΘEgyoto} ≈ 10^{5}μas at 200 parsecs).
Maximal errors evaluated for each parameter and each spin in μas.
The requirement on accuracy (as) is largely met in the WDL. We note that our method (see Appendix A) uses photons both inside and outside of the equatorial plane of the black hole. Therefore we have demonstrated our conclusion for any photon in the WDL, not only in the equatorial plane. However, an equivalent test is necessary in the SDL.
4.2. Comparison of GYOTO null geodesics with Geokerr in SDL
The aim of this subsection is to check if null geodesics computed with GYOTO in the SDL are accurate enough. To do so, we decided to compare photon trajectories computed with GYOTO with those computed with the raytracing code Geokerr. Contrary to GYOTO, Geokerr computes photon coordinates semianalytically, reducing the equations of motion expressed in the Hamiltonian formulation to Carlson elliptic integrals.
The comparison is made using the same observer coordinates and black hole parameters as before. We evaluate null geodesics for three different values of the spin (0.2, 0.5 and 0.998) and we consider photons launched from the center of the observer screen (α = δ = 0, see the beginning of Sect. 1 of the Appendix A for explanation of GYOTO null geodesics integrations). We first compute the geodesics with Geokerr and get the time coordinate of each point of the photon trajectory. Then, to get the null geodesics with GYOTO, we interpolate the positions of photons with our raytracing code, taking these dates into consideration. The maximal errors found on the photon coordinates is presented in Table 2. Errors correspond to the difference between positions evaluated with GYOTO and those evaluated with Geokerr. The maximal errors are about 10^{3}μas which shows a very good consistency between the two raytracing codes. Even for a photon, which is not launched from the center of the screen (α = 5μas and δ = −5μas), with a = 0.998, we find a very small error (≈ 4 × 10^{3}μas). An example of null geodesics computed with both codes is shown in Fig. 9.
As mentioned, we use the null geodesic integrated in our raytracing code to estimate the geodesic coordinates of GYOTO at Geokerr time coordinates. We noted that when the step size is high, the interpolation is less well evaluated. Thus, considering one fixed position on the null geodesic, each Cartesian coordinate is affected by a random error due to the interpolation in GYOTO. If we now focus on the maximal error along each Cartesian coordinate, we found that this error is constant when varying the tolerance parameter. It thus confirms that the error is only linked to the interpolation made in GYOTO.
These results show that GYOTO is accurate to a high level, even for large distances and in both weak and strongdeflection regimes. Next, we extend the study of Bozza & Mancini (2012) on lensing effects in the central parsec of our Galaxy and discuss the impact of neglecting lensing effects in future stellarorbit models on the fitted orbital parameters and the detection of the other effects affecting the orbit.
Maximal errors on the photon coordinates evaluated for each spin in μas.
Fig. 9 Null geodesics computed with GYOTO (dotted lines) and Geokerr (dashed lines). We consider a spin of 0.998 and a photon launched from the center of the screen. b) is a zoom of geodesics plotted on a): δ_{x} and δ_{y} on b) is about 4 × 10^{4}M. 

Open with DEXTER 
Input and fitted parameters for the S2 and E0 stars.
5. Short studies on lensing effects in the central parsec
In this section, we want to complete the study performed by Bozza & Mancini (2012) using the raytracing code GYOTO. Considering the astrometric accuracy of GRAVITY, they showed that this instrument is sensitive to the lensing effects generated by a black hole of 4.3 × 10^{6}M_{⊙}. They divided the parameter space in three regions depending on the position of the star. However, they only focused on stars located behind the plane of the sky containing Sgr A*. We present the astrometric shift of the primary image in Fig. 10. We consider stars with r_{s} ranging from 10^{4} to 1 parsec and vary the angle γ corresponding to the angle between the observer, the lens and the source, from 20° (the star is behind the plane of the black hole) to 135° (the star is in front of the plane of the sky containing Sgr A*). The shift presented here is obtained in the weakdeflection and weakfield regimes since we satisfy the conditions r_{s} ≫ R_{S} and r_{minγ} ≫ R_{S}. Indeed, the smallest r_{s} and r_{minγ} we get are of about 250 R_{S} and 90 R_{S}, respectively. We consider the same mass of the black hole and distance of the observer (r_{0} = 8.3 kpc) as Bozza & Mancini (2012).
We find the same results as Bozza & Mancini (2012). Although the effect appears to be small for a star in front of the plane of Sgr A* (e.g. ≈ 2μas at 110° inclination), it is a systematic effect which always displaces the image of the star away from the black hole. For this reason, any attempt at orbitfitting must take the lensing effects into account. For instance, the major axis of an orbit in the plane of the sky will appear 5 μas larger than it really is.
To get an idea of the impact of neglecting lensing effects in stellarorbit models on fitted orbital parameters, we decide to fit these parameters to astrometric and radial velocity data with this type of a model. The orbital parameters we fit are: the period of the orbit T, the semimajor axis a_{sma}, the eccentricity e, the time of the pericenter passage t_{p}, the inclination i, the position angle of the line of nodes Ω, and the angle from ascending node to pericenter ω. To simulate the positions and radial velocities, we use a Keplerian model and simulate the lensing effects using analytical formulae from Sereno & de Luca (2006) (only applied when the star is behind the plane of sky containing Sgr A*). Although not as accurate as a fullfledged numerical simulation, this simplistic model is sufficient for our purpose. No noise is added to the mock observations and we consider 1000 dates to generate one complete orbit. We simulated two orbits: one based on the bestfit orbital parameters of the S2 star (Gillessen et al. 2009), the other one for a fictional star (hereafter E0), which is closer to Sgr A* than S2. The parameters of each orbit are reported in Table 3. The impact of the lensing effects on the astrometric data is visible in Fig. 11, where we only consider a portion of the orbits of the stars.
Fig. 10 Astrometric shift in μas of the primary image owing to lensing effects. This plot is valid for sources ranging from 10^{4} to 1 parsec. The shift is the difference between lensed and unlensed angular positions of the star. 

Open with DEXTER 
Fig. 11 Lensing effects affecting a portion of the S2 orbit (dashed blue line) and the E0 orbit (solid black). These simulated shifts in the astrometric observations are obtained considering a Keplerian model taking into account the lensing effects. This effect is added by using the analytical formulas developed in Sereno & de Luca (2006). 

Open with DEXTER 
We then fit the mock data using a second Keplerian model, but neglecting the lensing effects. We finally recomputed mock (lensed) data using the fitted (nonlensed) parameters. The difference of the fitted parameters to the input parameters (also reported in Table 3) seem small in relative value (the maximal relative difference is around ≃ 10^{5} for S2 and ≃ 10^{3} for E0), but they are significant in the sense that using the fitted parameters instead of the input parameters in the mock data yields a significant error in the astrometry and in the radial velocity (up to 20 μas and 0.8 km/s respectively for S2, and 60 μas and 10 km s^{1}, respectively, for E0).
The other point that should be discussed is the impact of using models neglecting lensing effects, for the detection of other effects affecting the orbit. To simulate these different effects, we consider the raytracing code GYOTO and the E0 star. We estimate the maximal astrometric influence of different effects such as Roemer, pericenteradvance, Shapiro and LenseThirring. We find that, in two orbital periods, the maximal effects reach: ≈ 220μas for Roemer, ≈ 240μas for the pericenteradvance, ≈ 90μas for gravitational lensing (see the solid black curve on Fig. 11), and ≈ 8μas for both Shapiro and LenseThirring. In two periods of monitoring, the Roemer and pericenteradvance effects dominate the lensing effects most of the time. However, if we only consider the first pericenter passage, we note that the impact of lensing is of the same order of magnitude as the pericenteradvance and the Roemer effects. The gravitational lensing will thus probably interfere with detecting these effects. Beside, we find that higherorder effects such as Shapiro and LenseThirring are smaller than gravitational lensing. The detection of these smaller effects will be difficult if we do not take into account lensing effects in stellarorbit models.
To summarize, we need to take into account the lensing effects in future stellarorbit models. As shown here, it can not be neglected in most cases. If these effects are not included in stellarorbit models, they will probably prevent the measurement of other relativistic effects and we will not be able to efficiently constrain all key parameters, such as the spin of the black hole.
6. Conclusion
The Galactic center is a unique laboratory to observe stars close enough to a compact object to test general relativity. Thanks to GRAVITY it will be possible to measure astrometric positions of stars orbiting Sgr A* with an expected astrometric precision of 10 μas. We have shown that GYOTO is extremely accurate even in complex configurations. For the purpose of the interpretation of the future astrometric positions observed by GRAVITY, GYOTO is accurate enough to model star trajectories and fit the GRAVITY data.
The various integrators implemented in GYOTO are all very accurate and fast. In most scenarios, the default parameters yield good results. On the other hand, when the astrometric error must be minimized or simply evaluated, or when computing time needs to be optimized, some effort must be put into choosing the optimal integrator and parameters. The scenario that we exhibit in this paper is particularly challenging: since each photon needs to be integrated over a long distance (several parsecs) before and after the interaction with the black hole, the deflection angle must be estimated to an extreme precision. One method to evaluate the numerical errors for a given integrator and parameters is to compare a typical geodesic that has been computed with these parameters to the same geodesic that has been computed with one of the most accurate setup: Runge Kutta Fehlberg 78 and AbsTol= RelTol = 10^{18}.
In this paper, we also showed that lensing effects need to be taken into account to reach a model with an accuracy of 1 μas, even for stars located inside or in front of the plane of sky containing Sgr A*. The shift of the primary image can reach 5 μas for stars in this plane. The half angle of aperture, in which the astrometric shift is less than 1 μas, is about 50°. Out
of this cone, the lensing effects must be taken into account to avoid a systematic error in each modeled astrometric position. It is especially true for the next generation of instruments such as MICADO on the EELT (Davies et al. 2010), which are expected to observe stellar orbits with a high accuracy.
To constrain the nature of Sgr A* with future accurate instruments, we need to take into account the lensing effects in stellarorbit models. Otherwise, fitting with models that neglect these effects will lead to inaccurate orbital parameter estimations and prevent the measurement of others effects, such as the LenseThirring effect.
Freely available at the URL http://GYOTO.obspm.fr
Freely available at the URL http://www.astro.washington.edu/users/agol/geokerr/
Acknowledgments
This work was supported by the French ANR POLCA project (Processing of pOLychromatic interferometriC data for Astrophysics, ANR10BLAN0511), and by the OPTICON project (EC FP7 grant agreement 312430) in “Image reconstruction in optical interferometry” WP4. We want to thank, Eric Gourgoulhon, Claire Somé, Jason Dexter and Stefan Gillessen for helpful discussions. We also would like to thank the referee and Frédéric Vincent for their advice, which enabled us to improve the quality of the paper.
References
 Bozza, V. 2008, Phys. Rev. D, 78, 063014 [NASA ADS] [CrossRef] [Google Scholar]
 Bozza, V., & Mancini, L. 2012, ApJ, 753, 56 [NASA ADS] [CrossRef] [Google Scholar]
 Bozza, V., & Scarpetta, G. 2007, Phys. Rev. D, 76, 083008 [NASA ADS] [CrossRef] [Google Scholar]
 Bozza, V., de Luca, F., Scarpetta, G., & Sereno, M. 2005, Phys. Rev. D, 72, 083003 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Davies, R., Ageorges, N., Barl, L., et al. 2010, in Groundbased and Airborne Instrumentation for Astronomy III, Proc. SPIE, 7735, 77352A [CrossRef] [Google Scholar]
 Dexter, J., & Agol, E. 2009, ApJ, 696, 1616 [NASA ADS] [CrossRef] [Google Scholar]
 Eisenhauer, F., Perrin, G., Brandner, W., et al. 2011, The Messenger, 143, 16 [NASA ADS] [Google Scholar]
 Gillessen, S., Eisenhauer, F., Fritz, T. K., et al. 2009, ApJ, 707, L114 [NASA ADS] [CrossRef] [Google Scholar]
 Grandclément, P., Somé, C., & Gourgoulhon, E. 2014, Phys. Rev. D, 90, 024068 [NASA ADS] [CrossRef] [Google Scholar]
 Merritt, D., Alexander, T., Mikkola, S., & Will, C. M. 2010, Phys. Rev. D, 81, 062002 [NASA ADS] [CrossRef] [Google Scholar]
 Rauch, K. P., & Blandford, R. D. 1994, ApJ, 421, 46 [NASA ADS] [CrossRef] [Google Scholar]
 Schneider, P., Ehlers, J., & Falco, E. E. 1992, Gravitational Lenses (Heidelberg: Springer) [Google Scholar]
 Sereno, M., & de Luca, F. 2006, Phys. Rev. D, 74, 123009 [NASA ADS] [CrossRef] [Google Scholar]
 Sereno, M., & de Luca, F. 2008, Phys. Rev. D, 78, 023008 [NASA ADS] [CrossRef] [Google Scholar]
 Vincent, F. H., Paumard, T., Gourgoulhon, E., & Perrin, G. 2011, Class. Quant. Grav., 28, 225011 [NASA ADS] [CrossRef] [Google Scholar]
 Vincent, F. H., Paumard, T., Perrin, G., et al. 2014, MNRAS, 441, 3477 [NASA ADS] [CrossRef] [Google Scholar]
 Vincent, F. H., Meliani, Z., Grandclément, P., Gourgoulhon, E., & Straub, O. 2016, Class. quant. grav., 33, 105015 [NASA ADS] [CrossRef] [Google Scholar]
 Will, C. M. 2008, ApJ, 674, L25 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: How to determine primary caustics and critical curves parameters with GYOTO?
GYOTO is a raytracing code integrating null geodesics backward in time. Each pixel on the observer’s screen corresponds to the initial direction of each photon. When a photon reaches the star, the pixel illuminates and we get the star image, as presented in Fig. 3. However, we consider another type of image to estimate the three parameters Θ_{E}, Δ, and B_{C}. We use a mathematical function named MinDistance in GYOTO. This represents the square minimum distance between the photon and the surface of the star. Zeroes of this function correspond to photons, which have reached the surface of the star. The advantage of this function is that it enables images of the object to be located, even if no geodesic that is actually computed emanated from the star. Finding images of the object corresponds to finding minima of MinDistance, and checking that the minimum reaches zero. To illustrate this function, we consider a point source behind a Schwarzschild black hole (see Fig. A.1). The following subsections explain the measurement methods of the different parameters with the MinDistance function.
Appendix A.1: Angular position B_{Cgyoto} of the caustic
Appendix A.1.1: Caustic point
Firstly, we are interested in measuring the angular position of caustic points. This means that we focus on caustics obtained in the Schwarzschild case or in the Kerr case, when the star is very far from the black hole. We note this angle B_{Cp} in this subsection. To get the angular position of the caustic point, we need to estimate its position relative to the black hole, y_{Cp}. At this position, a critical curve is formed in the observer’s sky. Thus, we have to search the position of the source y_{Cp} leading to the formation of this critical curve. If the source is not on the caustic point, minima of MinDistance will be located on the αaxis and will correspond to the primary and secondary images. If the source lies on the caustic point, the two images merge and form the critical curve, and minima will be distributed all along the curve (see Fig. A.1). Thus, to converge to y_{Cp} we search the best minimum of MinDistance along the positive δaxis (α = 0) varying the position of the source y_{s}. If we consider a Schwarzschild lens, it means that we look for the value of y_{s}, enabling us to reach the best minimum at the top of the critical curve. However, with a Kerr black hole, the critical curve is shifted from the δaxis (the center of the critical curve is no longer at α = 0) but we still search the best minimum of MinDistance along the positive δaxis. Even if we are not exactly on the axis passing through the top of the critical curve, the best minimum is reached when the critical curve is formed (because best minimums will be distributed all along the critical curve).
The method we selected to find y_{Cp} is a golden section search. This method is the same as the bisection method, but divides the interval by the golden number instead of 2. The convergence is faster than the bisection method. The search process stops when the MinDistance function cannot be well estimated for two close y_{s}. The error on y_{Cp} is given by the size of the last step used to divide the final interval. Using error propagation, we get the error on B_{Cp} given by (A.1)The final error we consider is the biggest error evaluated for all distances of the source x_{s} we selected.
The stop conditions and the method to estimate the maximal errors for the next parameters are the same as above.
Fig. A.1 Illustration of the square root of the MinDistance function with a point source behind a Schwarzschild black hole. The black hole and the source are located at the center of the image. The source is at 100M. The darker the pixels, the smaller the distance between the photon and the point source. The colorbar is labeled in M. The minimal and maximal values of the image is about 7 × 10^{4}M and 102M (the maximal value of the colorbar is cut at 100M), respectively. 

Open with DEXTER 
Appendix A.1.2: Four cusps astroid
The first step to measure the angular position of the right cusp, noted B_{Cr}, is also to find its position relative to the black hole, y_{Cr}. To estimate y_{Cr}, we consider the primary caustics and critical curves properties that we outlined in Sect. 2.2. We use the link between the position of the source relative to the astroid, and images formed in the observer’s sky (source outside or on the caustic lead to the formation of two images, source inside the caustic leads to the formation of four images). If we consider a source at the position (x_{s}, y_{Cr}, 0), we observe two images in the observer plane, with one of these images forming on the left side of the critical curve (see the left triangle on the critical curve on the Fig. 4). The aim is to use a bisection method to converge to the right cusp. Initially, we consider two different positions of the source:

One inside the caustic, y_{C1}. This is obtained using the same method presented in Sect. A.1.1. We use a golden section search to look for the position of the source leading to the best minimum of the MinDistance function along the positive δaxis. In this case, it means that we search the position of the source leading to the formation of one image along the positive δaxis. If we take Fig. 4 into consideration, this image corresponds to the upper cross inside the critical curve.

One outside the caustic, y_{C2}. This is given by y_{C2} = y_{C1} ± ζ with ζ superior to Δ_{C}/ 2. The upper and the lower signs hold for the right cusp and for the left cusp, respectively.
The next value of y is given by the bisection method: y_{C3} = (y_{C1} + y_{C2})/2. At each division of the interval we always satisfy the condition: one source outside (two images) and one source inside (four images). To estimate the maximal error of B_{Cr} we also use the propagation error.
Appendix A.2: Radius Θ_{Egyoto} and offset Δ_{gyoto} of the critical curve
These two parameters are estimated thanks to simple relations that depend on the right and the left radii of the critical curve relative to the black hole and measured along the α axis where r_{l} and r_{r} are the left and right radii, respectively. In the Schwarzschild case the radii are equal. To measure the right and the left radii of the critical curve, we need to find the position of the left and the right cusps, respectively. We perform the method explained before to obtain the two positions of the two cusps. Then, we estimate the radii using another golden section search. The final error on Θ_{E} and Δ is given by (A.4)
where δ_{rl} = δ_{rl1} + δ_{rl2} and δ_{rr} = δ_{rr1} + δ_{rr2} are the errors on the left and right radii, respectively. δ_{rl1} and δ_{rr1} are the errors obtained with the golden section search that was used to converge to the two radii. This is given by the last step of the interval. δ_{rl2} and δ_{rr2} correspond to the errors on the two radii owing to the errors made in the evaluation of the right y_{Cr} and the left y_{Cl} position of the two cusps, respectively. To evaluate δ_{rl2} (δ_{rr2}), we estimate the left (right) radii obtained considering a source located at y_{Cr} + δ_{yCr} (y_{Cl} + δ_{yCl}) and y_{Cr} − δ_{yCr} (y_{Cl} − δ_{yCl}). We note the first radius r_{l2+} (r_{r2+}) and the second radius r_{l2−} (r_{r2−}). Thus, we get (A.5)and (A.6)For caustic points, we directly estimate the right and left radii at y_{Cp}. To estimate the errors, we just need to obtain the right and left radii for two values of y_{s}: y_{Cp} + δ_{yCp} and y_{Cp} − δ_{yCp}, and use the same process explained in the previous paragraph.
All Tables
All Figures
Fig. 1 Lensing configuration: O, S, and L denote the observer, the source, and the lens, respectively. 

Open with DEXTER  
In the text 
Fig. 2 Spatial projection of a Schwarzschild lensing situation: S corresponds to the source, L to the lens and O to the observer. 

Open with DEXTER  
In the text 
Fig. 3 GYOTO image of a fixed star, with a radius of 1M, located at 6M in front of a Schwarzschild black hole. The observer is at 8 kpc from the black hole. The star lies on the secondorder caustic line: the biggest image corresponds to the primary image of the star and the ring to the critical curve formed thanks to the secondorder caustic. Axes are labeled in μas and the field of view of the image is equal to 70 μas. 

Open with DEXTER  
In the text 
Fig. 4 Effect on the displacement of a point source, in the equatorial plane of the black hole, with respect to the caustic (top) on the formation of images in the observer plane (bottom). The black astroid represents the primary caustic and the black circle represents the primary critical curve. We define the z′ and y′ axis whose origin corresponds to the middle of the astroid crosssection. The α and δ axis correspond to the observer’s sky axis, the origin is centered in the middle of the screen (because of the rotation of the black hole, the critical curve is shifted from the center of the screen). 

Open with DEXTER  
In the text 
Fig. 5 Offset estimated with GYOTO in μas versus the distance  x_{s}  of the point source in parsec, in the Schwarzschild case. Each plot represents different values of the tuning parameters h and AbsTol: 10^{18} (left), 10^{16} (middle), and 10^{14} (right). The different types of curves denote the different integrators: Legacy Generic Integrator (solid), Runge Kutta Cash Karp 54 (dash), and Runge Kutta Fehlberg 78 (dot). 

Open with DEXTER  
In the text 
Fig. 6 Comparison of the three integrators by measuring the offset of the Einstein ring in the Schwarzschild case. Left: computing time versus the tuning parameters. Middle: offset versus the tuning parameters. Right: computing time versus the offset. The different types of curves denote the different integrators: Legacy Generic Integrator (solid), Runge Kutta Cash Karp 54 (dash), and Runge Kutta Fehlberg 78 (dot). 

Open with DEXTER  
In the text 
Fig. 7 Absolute difference between analytical approximations (7), (8), and (11), and GYOTO measurements for the three parameters B_{C}, Θ_{E}, and Δ. The types of line denote different values of the spin: 0.2 in solid, 0.5 in dotted, and 0.9 in dashed. 

Open with DEXTER  
In the text 
Fig. 8 Absolute difference between analytical approximations and GYOTO measurements considering three different equations for the approximation Θ_{E}: the solid line is obtained considering only the zeroth and first orders in ε in the Eq. (7), the dashed line considers the zeroth, first, and second orders in ε and the dotted line considers the full equation. Only a = 0.5 is plotted. 

Open with DEXTER  
In the text 
Fig. 9 Null geodesics computed with GYOTO (dotted lines) and Geokerr (dashed lines). We consider a spin of 0.998 and a photon launched from the center of the screen. b) is a zoom of geodesics plotted on a): δ_{x} and δ_{y} on b) is about 4 × 10^{4}M. 

Open with DEXTER  
In the text 
Fig. 10 Astrometric shift in μas of the primary image owing to lensing effects. This plot is valid for sources ranging from 10^{4} to 1 parsec. The shift is the difference between lensed and unlensed angular positions of the star. 

Open with DEXTER  
In the text 
Fig. 11 Lensing effects affecting a portion of the S2 orbit (dashed blue line) and the E0 orbit (solid black). These simulated shifts in the astrometric observations are obtained considering a Keplerian model taking into account the lensing effects. This effect is added by using the analytical formulas developed in Sereno & de Luca (2006). 

Open with DEXTER  
In the text 
Fig. A.1 Illustration of the square root of the MinDistance function with a point source behind a Schwarzschild black hole. The black hole and the source are located at the center of the image. The source is at 100M. The darker the pixels, the smaller the distance between the photon and the point source. The colorbar is labeled in M. The minimal and maximal values of the image is about 7 × 10^{4}M and 102M (the maximal value of the colorbar is cut at 100M), respectively. 

Open with DEXTER  
In the text 