© ESO, 2016

1. Introduction

General relativitY OrbiT of Observatoire de Paris (GYOTO¹) is a ray-tracing 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 second-generation 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 Lense-Thirring 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′′ field-of-view 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 weak-deflection limit (WDL), however we also have to check if this ray-tracing code is valid in the strong-deflection limit (SDL). To do so, we choose to compare null geodesics computed in GYOTO and with another code named Geokerr² (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.

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 z-axis. Spherical coordinates of the observer and the source, relative to the lens L, are noted (r₀,ϑ₀,φ₀) and (r_s,ϑ_s,φ_s), respectively. Without loss of generality, we choose to work in the equatorial plane of the black hole: ϑ₀ = π/ 2, φ₀ = 0 and ϑ_s = π/ 2. This yields (x₀,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² 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.

Using the notation of Fig. 2, we can write the lens equation (Schneider et al. 1992) $\begin{array}{ccc}{\beta }\mathrm{=}{\theta }\mathrm{-}\frac{{\mathit{r}}_{\mathrm{s}}}{\mathrm{(}{\mathit{r}}_{\mathrm{s}}\mathrm{+}{\mathit{r}}_{\mathrm{0}}\mathrm{)}}{\alpha ̂}{\left(}{\xi }{\right)}\mathit{,}& & \end{array}$(1)with β the unlensed angular position of the source, θ the lensed angular position of the source equal to ξ/r₀ and the deflection angle depending on the impact parameter ξ. The latter angle is given by $\begin{array}{ccc}{\alpha ̂}{\left(}{\xi }{\right)}\mathrm{=}\frac{\mathrm{2}{\mathit{R}}_{\mathrm{S}}}{{\xi }}\mathit{,}& & \end{array}$(2)with R_S = 2GM/c² the Schwarzschild radius. We can rewrite the lens equation as $\begin{array}{ccc}{{\theta }}^{\mathrm{2}}\mathrm{-}{\beta \theta }\mathrm{-}{\mathit{\alpha }}_{\mathrm{0}}^{\mathrm{2}}\mathrm{=}\mathrm{0}\mathit{,}& & \end{array}$(3)with $\begin{array}{ccc}{\mathit{\alpha }}_{\mathrm{0}}\mathrm{=}\sqrt{\mathrm{2}{\mathit{R}}_{\mathrm{S}}\frac{{\mathit{r}}_{\mathrm{s}}}{{\mathit{r}}_{\mathrm{0}}\mathrm{(}{\mathit{r}}_{\mathrm{s}}\mathrm{+}{\mathit{r}}_{\mathrm{0}}\mathrm{)}}}\mathit{,}& & \end{array}$(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 $\mathit{A}\mathrm{=}|\mathit{det}\frac{\mathit{\partial }{\beta }}{\mathit{\partial }{\theta }}{|}^{-1}\mathrm{·}$(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 θ = α₀. A circle called a critical curve is formed in the lens plane with a radius of α₀. Considering the source as a star, the observer sees the well-known Einstein ring. The radius of the ring corresponds to the critical curve radius so we get α₀ = θ_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 second-order caustic line: the biggest image corresponds to the primary image of the star and the ring to the critical curve formed thanks to the second-order caustic. Axes are labeled in μas and the field of view of the image is equal to 70 μas.

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 first-order (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 higher-order 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 (four-cusped) cross-section (see the upper scheme in Fig. 4). At large distances the cross-section 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 cross-section) decreases and tends towards zero.

To form the critical curve the source must cover all of the astroid cross-sections. 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 cross-section. 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).

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 weak-deflection and weak-field 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 weak-field 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 $\mathit{\epsilon }\mathrm{=}\frac{{\mathit{\theta }}_{\mathrm{E}}}{\mathrm{4}\mathit{D}}\mathit{,}$(6)where D = r_s/ (r_s + r₀) 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: $\begin{array}{ccc}{\mathrm{\Theta }}_{\mathrm{E}}& \mathrm{\approx }& {\mathit{\theta }}_{\mathrm{E}}\{\mathrm{1}\mathrm{+}\frac{\mathrm{15}\mathit{\pi }}{\mathrm{32}}\mathit{\epsilon }\mathrm{+}\left[\mathrm{4}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{D}}^{\mathrm{2}}\mathrm{)}\mathrm{-}\frac{\mathrm{675}{\mathit{\pi }}^{\mathrm{2}}}{\mathrm{2048}}\right]{\mathit{\epsilon }}^{\mathrm{2}}\\ & & \end{array}$(7)The left and the right radius of the critical curve are equal and depend on the spin in the third-order term in ε. The second equation gives the offset Δ of the center of the critical curve relative to the black hole: $\mathrm{\Delta }\mathrm{\approx }\mathrm{4}\mathit{D}{\mathit{\epsilon }}^{\mathrm{2}}\mathit{a}\{\mathrm{1}\mathrm{+}\frac{\mathrm{15}\mathit{\pi }}{\mathrm{32}}\mathit{\epsilon }\mathrm{+}\left[\mathrm{4}\mathrm{(}\mathrm{2}\mathrm{-}\mathit{D}\mathrm{+}\mathrm{4}{\mathit{D}}^{\mathrm{2}}\mathrm{)}\mathrm{-}\frac{\mathrm{225}{\mathit{\pi }}^{\mathrm{2}}}{\mathrm{256}}\right]{\mathit{\epsilon }}^{\mathrm{2}}\}\mathit{.}$(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: $\begin{array}{ccc}{\mathit{y}}_{\mathrm{C}}\mathrm{\approx }\frac{\mathit{a}}{\mathrm{1}\mathrm{-}\mathit{D}}\left[\mathrm{1}\mathrm{+}\frac{\mathrm{5}\mathit{\pi }}{\mathrm{16}}\mathit{\epsilon }\mathrm{+}\left(\mathrm{4}\mathrm{-}\frac{\mathrm{225}{\mathit{\pi }}^{\mathrm{2}}}{\mathrm{512}}{\mathit{\epsilon }}^{\mathrm{2}}\right)\right]\mathrm{+}{\mathrm{\Delta }}_{\mathrm{C}}\cos \mathrm{3}\mathit{\eta }\mathit{,}& & \end{array}$(9)with $\begin{array}{ccc}{\mathrm{\Delta }}_{\mathrm{C}}\mathrm{\approx }\frac{\mathrm{15}\mathit{\pi }}{\mathrm{256}}\frac{{\mathit{a}}^{\mathrm{2}}}{\mathit{D}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{D}\mathrm{)}{\mathit{r}}_{\mathrm{0}}}\mathit{,}& & \end{array}$(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${\mathit{B}}_{\mathrm{C}}\mathrm{\approx }\arctan \frac{{\mathit{y}}_{\mathrm{C}}}{{\mathit{x}}_{\mathrm{0}}\mathrm{+}\mathrm{\left|}{\mathit{x}}_{\mathrm{s}}\mathrm{\right|}}\mathrm{·}$(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₀ = 8 kpc from a black hole of mass M equal to 4.31 × 10⁶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 | xs | 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).

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).

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 fourth-order 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³. 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 run-time 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^-6M 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³/ 10⁸ 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³/ 10⁸ and not h for the Legacy integrator because the tuning parameters AbsTol and h are not defined in the same way. The quantity h³/ 10⁸ 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 BC, Θ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.

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³/ 10⁸ 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³/ 10⁸ = 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 computing-time of less than one minute. Since the following studies do not use time-consuming 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.

On each plot two different regimes can be observed. For small | x_s |, the curve is marked by a smooth, power-law 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 ε (ε⁰ and ε) in Eq. (7). For the second absolute difference, we consider one more order in ε (ε⁰, ε and ε²). 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 ray-tracing code is very accurate, even for sources far behind the black hole (e.g. δ_ΘEgyoto ≈ 10^-5μas at 200 parsecs).

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 ray-tracing code Geokerr. Contrary to GYOTO, Geokerr computes photon coordinates semi-analytically, 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 ray-tracing 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 ray-tracing 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 ray-tracing 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 strong-deflection 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 stellar-orbit models on the fitted orbital parameters and the detection of the other effects affecting the orbit.

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-4M.

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 ray-tracing 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⁶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 weak-deflection and weak-field 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₀ = 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 orbit-fitting 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 stellar-orbit 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 semi-major 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 full-fledged 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 best-fit 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.

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).

We then fit the mock data using a second Keplerian model, but neglecting the lensing effects. We finally re-computed mock (lensed) data using the fitted (non-lensed) 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 ray-tracing code GYOTO and the E0 star. We estimate the maximal astrometric influence of different effects such as Roemer, pericenter-advance, Shapiro and Lense-Thirring. We find that, in two orbital periods, the maximal effects reach: ≈ 220μas for Roemer, ≈ 240μas for the pericenter-advance, ≈ 90μas for gravitational lensing (see the solid black curve on Fig. 11), and ≈ 8μas for both Shapiro and Lense-Thirring. In two periods of monitoring, the Roemer and pericenter-advance 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 pericenter-advance and the Roemer effects. The gravitational lensing will thus probably interfere with detecting these effects. Beside, we find that higher-order effects such as Shapiro and Lense-Thirring are smaller than gravitational lensing. The detection of these smaller effects will be difficult if we do not take into account lensing effects in stellar-orbit models.

To summarize, we need to take into account the lensing effects in future stellar-orbit models. As shown here, it can not be neglected in most cases. If these effects are not included in stellar-orbit 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 E-ELT (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 stellar-orbit 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 Lense-Thirring effect.

Acknowledgments

This work was supported by the French ANR POLCA project (Processing of pOLychromatic interferometriC data for Astrophysics, ANR-10-BLAN-0511), 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

Appendix A: How to determine primary caustics and critical curves parameters with GYOTO?

GYOTO is a ray-tracing 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 ${\mathit{\delta }}_{{\mathit{B}}_{\mathrm{Cp}}}\mathrm{=}{\mathit{\delta }}_{{\mathit{y}}_{\mathrm{Cp}}}\left|\frac{\mathrm{1}}{{\mathit{x}}_{\mathrm{0}}\mathrm{+}\mathrm{\left|}{\mathit{x}}_{\mathrm{s}}\mathrm{\right|}\mathrm{+}{\mathit{y}}_{\mathrm{Cp}}}\right|\mathrm{·}$(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 color-bar is labeled in M. The minimal and maximal values of the image is about 7 × 10-4M and 102M (the maximal value of the color-bar is cut at 100M), respectively.

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 $\begin{array}{ccc}& & {\mathrm{\Theta }}_{\mathrm{Egyoto}}\mathrm{=}\frac{{\mathit{r}}_{\mathrm{l}}\mathrm{+}{\mathit{r}}_{\mathrm{r}}}{\mathrm{2}}\mathit{,}\\ & & {\mathrm{\Delta }}_{\mathrm{gyoto}}\mathrm{=}\mathrm{|}{\mathrm{\Theta }}_{\mathrm{E}}\mathrm{-}{\mathit{r}}_{\mathrm{r}}\mathrm{|}\mathit{,}\end{array}$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 ${\mathit{\delta }}_{{\mathrm{\Theta }}_{\mathrm{Egyoto}}}\mathrm{=}{\mathit{\delta }}_{{\mathrm{\Delta }}_{\mathrm{gyoto}}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{\mathit{\delta }}_{{\mathit{r}}_{\mathrm{l}}}\mathrm{+}{\mathit{\delta }}_{{\mathit{r}}_{\mathrm{r}}}\mathrm{)}\mathit{,}$(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 ${\mathit{\delta }}_{{\mathit{r}}_{\mathrm{l}\mathrm{2}}}\mathrm{=}\mathrm{|}{\mathit{r}}_{\mathrm{l}{\mathrm{2}}_{\mathrm{+}}}\mathrm{-}{\mathit{r}}_{\mathrm{l}{\mathrm{2}}_{\mathrm{-}}}\mathrm{|}\mathit{,}$(A.5)and ${\mathit{\delta }}_{{\mathit{r}}_{\mathrm{r}\mathrm{2}}}\mathrm{=}\mathrm{|}{\mathit{r}}_{\mathrm{r}{\mathrm{2}}_{\mathrm{+}}}\mathrm{-}{\mathit{r}}_{\mathrm{r}{\mathrm{2}}_{\mathrm{-}}}\mathrm{|}\mathit{,}$(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.
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.
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 second-order caustic line: the biggest image corresponds to the primary image of the star and the ring to the critical curve formed thanks to the second-order caustic. Axes are labeled in μas and the field of view of the image is equal to 70 μas.

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 cross-section. 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).

In the text

Fig. 5 Offset estimated with GYOTO in μas versus the distance | xs | 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).

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).

In the text

	Fig. 7 Absolute difference between analytical approximations (7), (8), and (11), and GYOTO measurements for the three parameters BC, Θ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.
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.
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-4M.
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.
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).
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 color-bar is labeled in M. The minimal and maximal values of the image is about 7 × 10-4M and 102M (the maximal value of the color-bar is cut at 100M), respectively.

In the text