GYOTO: a new general relativistic ray-tracing code

We use the standard Boyer–Lindquist (x^α) = (t, r, θ, φ) coordinates to express the metric according to [27]:

$$\begin{eqnarray} \fl {\mathrm{d}}s^{2} &=& g_{\alpha \beta } \, {\mathrm{d}}x^{\alpha }{\mathrm{d}}x^{\beta } \nonumber \\ \fl &=& -\left(1 - \frac{2 \,M\, r}{\Sigma } \right) {\mathrm{d}}t^{2} - \frac{4 \,M \,a \,r \,\sin ^{2}\theta }{\Sigma } \, {\mathrm{d}}t\, {\mathrm{d}}\varphi + \frac{\Sigma }{\Delta }\, {\mathrm{d}}r^{2} + \Sigma \, {\mathrm{d}}\theta ^{2} \nonumber \\ \fl & & + \left( r^{2} + a^{2} + \frac{2\, M\, a^{2}\, r\, \sin \theta }{\Sigma } \right) \sin ^{2} \theta \, {\mathrm{d}}\varphi ^{2} , \end{eqnarray} \tag{ 1 }$$

where M is the black hole's mass, a is its spin parameter (aM being the total angular momentum), Σ ≡ r² + a² cos ²θ and Δ ≡ r² − 2 M r + a².

For a particle of mass μ and 4-momentum p_α, two constants of motion are provided by the spacetime symmetries: the energy 'at infinity' E ≡ −p_t and the axial component of the angular momentum L ≡ p_φ. The constant nature of E and L results, respectively, from the spacetime stationarity and the spacetime axisymmetry. In addition, the specific structure of Kerr spacetime gives birth to a third constant of motion, the Carter constant [28]:

$$\begin{equation} Q \equiv p_{\theta }^{2} + \cos ^{2} \theta \left[ a^{2}(\mu ^{2} - E^{2}) + \sin ^{-2}\theta \, L^{2}\right]. \end{equation} \tag{ 2 }$$

The problem is thus completely integrable by using the four constants (μ, E, L, Q).

Following [15], we use a Hamiltonian formulation for the equations of geodesics, which consists in using the variables (t, r, θ, φ, p_t, p_r, p_θ, p_φ) to express the equations of motion according to

$$\begin{eqnarray} \dot{t} &=& \frac{1}{2\, \Delta \Sigma } \frac{\partial }{\partial E} \left( R + \Delta \Theta \right) \end{eqnarray} \tag{ 3a }$$

$$\begin{eqnarray} \dot{r} &=& \frac{\Delta }{\Sigma }p_{r} \end{eqnarray} \tag{ 3b }$$

$$\begin{eqnarray} \dot{\theta } &=& \frac{1}{\Sigma }p_{\theta } \end{eqnarray} \tag{ 3c }$$

$$\begin{eqnarray} \dot{\varphi } &=& -\frac{1}{2\, \Delta \Sigma } \frac{\partial }{\partial L} \left( R + \Delta \Theta \right) \end{eqnarray} \tag{ 3d }$$

$$\begin{eqnarray} \dot{p_{t}} &=& 0 \end{eqnarray} \tag{ 3e }$$

$$\begin{eqnarray} \dot{p_{r}} &=& -\left( \frac{\Delta }{2\, \Sigma }\right)^{\prime } p_{r}^{2} - \left( \frac{1}{2\, \Sigma }\right)^{\prime } p_{\theta }^{2} + \left( \frac{R+\Delta \Theta }{2 \, \Delta \Sigma } \right)^{\prime } \end{eqnarray} \tag{ 3f }$$

$$\begin{eqnarray} \dot{p_{\theta }} &=& -\left( \frac{\Delta }{2\, \Sigma }\right)^{\theta } p_{r}^{2} - \left( \frac{1}{2\, \Sigma }\right)^{\theta } p_{\theta }^{2} + \left( \frac{R+\Delta \Theta }{2 \, \Delta \Sigma } \right)^{\theta } \end{eqnarray} \tag{ 3g }$$

$$\begin{eqnarray} \dot{p_{\varphi }} &=& 0, \end{eqnarray} \tag{ 3h }$$

where the superscripts ' and ^θ denote differentiation with respect to r and θ, respectively, a dot denoting differentiation with respect to proper time (for timelike geodesics) or to an affine parameter (for null geodesics). In addition, the following abbreviations have been used:

$$\begin{eqnarray} \Theta &\equiv &Q - \cos ^{2} \left[ a^{2} (\mu ^{2}-E^{2}) + \sin ^{-2}\theta L^{2}\right] \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} R&\equiv &\left[ E (r^{2} + a^{2}) - a L\right]^{2} - \Delta \left[ \mu ^{2} r^{2} + (L- a E)^{2} + Q \right]. \end{eqnarray} \tag{ 5 }$$

Let us consider the integration of a single null geodesic by GYOTO. The initial conditions are the position of the observer and the direction of incidence of the photon. These quantities allow us to determine the tangent vector to the photon's geodesic at the observer's position.

The integration is then performed backward in time, from the observer toward a distant astrophysical object. This backward integration is of particular use when the target object has a large angular size on the observer's sky: a large fraction of the computed geodesics will hit the object. If the integration was done from the object, most of the geodesics would not reach the observer, leading to a waste of computing time. However, owing to its modular nature, GYOTO could be fairly easily extended to support the forward ray-tracing paradigm.

The equations of motion are solved by means of a Runge–Kutta algorithm of fourth order (RK4), with an adaptive step (see e.g. [29]). The integration goes on until one of the following stop conditions is fulfilled:

• the photon reaches the emitting object (however, see section 2.2.1 for the particular treatment inside optically thin objects),
• the photon escapes too far from the target object (what 'far' means depends of the kind of object considered),
• the photon approaches too closely to the event horizon. This limit is typically defined as when the photon's radial coordinate becomes only a few percent larger than the radial coordinate of the event horizon. However, this limit must not be too large in order not to stop the integration of photons swirling close to the black hole before reaching the observer, thereby creating higher order images.

In order to ensure a proper integration, the conservation of the constants of motion is continuously checked, and the solution of the integration is modified to impose this conservation. This choice of modifying the computed coordinates to ensure the conservation of constants is a consequence of our choice of Hamiltonian formulation for the geodesic equations. The usual form of Kerr geodesic equations (see e.g. [27]), which directly enforces the conservation of constants, involves square roots that lead to sign indeterminations. The Hamiltonian formulation allows us to get rid of this problem [15], at the expense of slightly modifying the solution as explained below.

The constancy of E and L is directly enforced by equations (3e) and (3h). The constancy of Q is imposed by modifying the value of $\dot{\theta }$ given by the RK4 algorithm, according to (2) and (3c):

$$\begin{equation} \dot{\theta }= \pm \frac{1}{\Sigma } \sqrt{Q - \cos ^{2} \theta [ a^{2}(\mu ^{2} - E^{2}) + \sin ^{-2}\theta \, L^{2}] } , \end{equation} \tag{ 6 }$$

where the ± sign is chosen according to the sign of $\dot{\theta }$ resulting from the RK4 algorithm. Finally, the constancy of the 4-momentum's norm is ensured by modifying the $\dot{r}$ coordinate furnished by the RK4 algorithm. For instance, when considering a null geodesic, the value of $\dot{r}$ is chosen to enforce the relation

$$\begin{equation} g_{tt}\,\dot{t}^{2}+2\,g_{t\varphi }\,\dot{t}\,\dot{\varphi }+g_{rr}\,\dot{r}^{2}+g_{\theta \theta }\,\dot{\theta }^{2}+g_{\varphi \varphi }\,\dot{\varphi }^{2} = 0. \end{equation} \tag{ 7 }$$

Thus,

$$\begin{equation} \dot{r} = \pm \sqrt{-\frac{g_{tt}\,\dot{t}^{2}+2\,g_{t\varphi }\,\dot{t}\,\dot{\varphi }+g_{\theta \theta }\,\dot{\theta }^{2}+g_{\varphi \varphi }\,\dot{\varphi }^{2}}{g_{rr}}}, \end{equation} \tag{ 8 }$$

where the ± sign is chosen according to the sign of $\dot{r}$ resulting from the RK4 algorithm. However, in order not to depart too much from the RK4 output, the maximum modification of $\dot{\theta }$ and $\dot{r}$ is limited to 1% of their initial values.

The numerical accuracy of GYOTO is investigated in the appendix where a convergence test is given.

Once the null geodesic, integrated backward in time from the observer's screen, reaches the emitting object, the equation of radiative transfer must be integrated along the part of the geodesic that lies inside the emitter. In order to perform this computation, two basic quantities have to be known at each integration step: the emission coefficient j_ν and the absorption coefficient α_ν in the frame of a given observer (typically the observer comoving with the emitting matter). These coefficients are related to the specific intensity's increment dI_ν as along the geodesic according to

$$\begin{equation} {\mathrm{d}}I_{\nu } = - \alpha _{\nu }\,I_{\nu }\,{\mathrm{d}}s \qquad \mbox{and}\qquad {\mathrm{d}}I_{\nu } = j_{\nu }\,{\mathrm{d}}s , \end{equation} \tag{ 9 }$$

where ds is the element of length as measured by the observer. The impact of scattering is not taken into account by the code. However, it is possible, if needed, to include in the absorption and emission coefficients the effect of scattering. Moreover, the modular nature of GYOTO makes it possible to develop a specific scattering treatment that would be plugged into the existing code.

As shown in [30], the relativistic equation of radiative transfer is

$$\begin{equation} \frac{{\mathrm{d}}\mathcal {I}}{{\mathrm{d}}\lambda } = \mathcal {E} - \mathcal {A} \, \mathcal {I} , \end{equation} \tag{ 10 }$$

where λ is an affine parameter along the considered geodesic, and where the invariant specific intensity $\mathcal {I}$, the invariant emission coefficient $\mathcal {E}$ and the invariant absorption coefficient $\mathcal {A}$ have been used:

$$\begin{equation} \mathcal {I}\equiv \frac{I_{\nu }}{\nu ^{3}}, \qquad \mathcal {E}\equiv \frac{j_{\nu }}{\nu ^{2}}, \qquad \mathcal {A} \equiv \nu \, \alpha _{\nu }, \end{equation} \tag{ 11 }$$

ν being the frequency of the radiation. These quantities are invariant in the sense that they do not depend on the reference frame in which they are evaluated.

Let us now consider the reference frame comoving with the fluid emitting the radiation. The invariant equation (10) becomes in this frame

$$\begin{equation} \frac{{\mathrm{d}}I_{\nu _{\mathrm{em}}}}{{\mathrm{d}}\lambda } = \nu _{\mathrm{em}} [ j_{\nu _{\mathrm{em}}} - \alpha _{\nu _{\mathrm{em}}} \, I_{\nu _{\mathrm{em}}}], \end{equation} \tag{ 12 }$$

where ν_em is the emitted frequency of the radiation. It is straightforward to obtain the following relation:

$$\begin{equation} {\mathrm{d}}\lambda \, \nu _{\mathrm{em}} = {\mathrm{d}}s_{\mathrm{em}} , \end{equation} \tag{ 13 }$$

where ds_em is the amount of proper length as measured by the emitter. This expression leads to

$$\begin{equation} \frac{{\mathrm{d}}I_{\nu _{\mathrm{em}}}}{{\mathrm{d}}s_{\mathrm{em}}} = j_{\nu _{\mathrm{em}}} - \alpha _{\nu _{\mathrm{em}}} \, I_{\nu _{\mathrm{em}}}, \end{equation} \tag{ 14 }$$

where all quantities are computed in the emitter's frame: this is the standard form of the equation of radiative transfer in the emitter's frame.

Equation (14) can be immediately integrated between some value s₀ where the specific intensity is vanishing (the 'opposite' border of the emitting region) and some position s:

$$\begin{equation} I_{\nu } (s) = \int _{s_{0}}^{s} \mathrm{exp}\left( - \int _{s^{\prime }}^{s} \alpha _{\nu }(s^{\prime \prime })\, {\mathrm{d}}s^{\prime \prime }\right)\,j_{\nu }(s^{\prime })\, \,{\mathrm{d}}s^{\prime }. \end{equation} \tag{ 15 }$$

Provided that the absorption and emission coefficients are furnished at each position of spacetime, (15) allows one to compute the integrated specific intensity. This will be illustrated in section 2.3.

When considering an optically thick object, the computation of a spectrum is quite straightforward. Let us consider the spectrum emitted by a geometrically thin, optically thick accretion disk around a Schwarzschild black hole. Following [6], the specific intensity emitted at a given position of the disk is written as

$$\begin{equation} I_{\nu _{\mathrm{em}}} \propto \delta (\nu _{\mathrm{em}} - \nu _{\mathrm{line}})\,\epsilon (r), \end{equation} \tag{ 16 }$$

where δ is the Dirac distribution and where (r) follows a power law with parameter p:

$$\begin{equation} \epsilon (r) \propto r^{-p}. \end{equation} \tag{ 17 }$$

By means of the invariant intensity $\mathcal {I}=I_{\nu } / \nu ^{3}$, the specific intensity observed by a distant observer can be related to the emitted specific intensity according to

$$\begin{equation} I_{\nu _{\mathrm{obs}}} = g^{3}\,I_{\nu _{\mathrm{em}}} , \end{equation} \tag{ 18 }$$

where ν_obs is the observed frequency, and

$$\begin{equation} g \equiv \frac{\nu _{\mathrm{obs}}}{\nu _{\mathrm{em}}}. \end{equation} \tag{ 19 }$$

The observed flux F_ν is then related to the observed specific intensity according to

$$\begin{equation} {\mathrm{d}}F_{\nu _{\mathrm{obs}}} = I_{\nu _{\mathrm{obs}}}\,\cos \theta \,{\mathrm{d}}\Omega , \end{equation} \tag{ 20 }$$

where Ω is the solid angle under which the emitting element is seen, and θ is the angle between the normal to the observer's screen and the direction of incidence. The GYOTO screen is a very simplistic model for a physical telescope equipped with a spectro-imaging camera. It this model, the screen is assumed to be point-like. A pixel on this screen corresponds to a direction on the sky, just like a pixel on a camera detector. Our screen object is defined, among others, by the field-of-view it covers on the sky, the number of samples (or pixels) to cover this field-of-view and the spectral properties of the detector (number of spectral channels and their wavelength). The screen covers a solid angle Ω on the sky and each pixel screen corresponds to a small solid angle around a given direction of incidence of the photons. Hence,

$$\begin{equation} F_{\nu _{\mathrm{obs}}} = \sum _{\mathrm{pixels}} I_{\nu _{\mathrm{obs}},\mathrm{pixel}}\,\cos (\theta _{\mathrm{pixel}})\,\delta \Omega _{\mathrm{pixel}} , \end{equation} \tag{ 21 }$$

where $I_{\nu _{\mathrm{obs}},\mathrm{pixel}}$ is the specific intensity reaching the given pixel, θ_pixel is the angle between the normal to the screen and the direction of incidence corresponding to this pixel and δΩ_pixel is the element of solid angle covered by this pixel on the sky, defined as the total solid angle covered by the screen divided by the number of pixels:

$$\begin{equation} \delta \Omega _{\mathrm{pixel}} = \frac{2\,\pi \,(1-\cos f)}{N_{\mathrm{pixels}}}, \end{equation} \tag{ 22 }$$

where f is the angle between the normal to the screen and the most external incident direction of photons on the screen. Figure 1 shows the resulting emitted spectrum for a Schwarzschild black hole seen under an inclination of 45°.

Zoom In Zoom Out Reset image size

When considering an optically thin object, the spectrum can also be computed quite easily. It is possible to compute maps of specific intensities for a given panel of N observed frequencies ν_{obs, i}, 1 ⩽ i ⩽ N, by using (15), and to compute the observed flux at these frequencies by using (21). Illustrations of such computations will be given in section 2.3.

GYOTO can compute the image of a moving star, orbiting around a Kerr black hole. The model for the star is very simple, only the timelike geodesic of the star's center is computed, and the star is defined as the points whose Euclidian distance to the center is less than a given radius R, i.e. the points whose Cartesian-like coordinates (x ≡ r sin θ cos φ, y ≡ r sin θ sin φ, z ≡ r cos θ) satisfy

$$\begin{equation} (x-x_{c})^{2}+(y-y_{c})^{2}+({z}-{z}_{c})^{2} \le R^{2}, \end{equation} \tag{ 23 }$$

where (x_c, y_c, z_c) are the stellar center's coordinates. No internal physics nor tidal effects are taken into account.

Figure 2 shows the resulting image³ for a star orbiting on the innermost stable circular orbit (ISCO) of a Kerr black hole of spin parameter a = 0.9 M. Four particular rays are depicted on the right side of the figure, in order to show the trajectory followed by photons responsible for the primary, secondary, tertiary and quaternary images. The corresponding geodesics swirl around the black hole zero, one, two or three times before reaching the object.

Zoom In Zoom Out Reset image size

The capacity of GYOTO to compute the trajectories of stars around a Kerr black hole will allow us to develop relativistic orbit-fitting scripts, in particular by using the Yorick package presented in section 4.

A basic object implemented in GYOTO is the geometrically thin infinite accretion disk, with interior radius corresponding to the ISCO [31]. The expression of the emitted flux is given in [12], and is proportional to the emitted specific intensity provided the emission is isotropic, yielding

$$\begin{equation} I_{\nu } \propto \frac{1}{(\rho ^{2}-3)\,\rho ^{5}}\left\lbrace \rho - \frac{3}{\sqrt{2}}\,\ln \left[(3-2\sqrt{2})\,\frac{\rho +\sqrt{3}}{\rho +\sqrt{3}}\right] \right\rbrace , \end{equation} \tag{ 24 }$$

where $\rho = \sqrt{r}/M$, M being the central black hole's mass. The resulting image is depicted in figure 3 for the case of a Schwarzschild black hole and an inclination angle of 70°. It is in good agreement with the results of [12].

Zoom In Zoom Out Reset image size

In a geometrically thin disk, a Rossby wave instability (RWI) can be triggered if the density profile exhibits an extremum for some value of the radius r_RWI [32]. Density waves will be emitted away from the extremum region, with vortices appearing in the vicinity of r_RWI. Such a disk has been implemented in GYOTO, the density profile being computed by means of the VAC code (see [33]). The observed flux emitted by such a disk is computed according to the prescription in [34], as a function of the surface density and radius. The result is presented in figure 4, which can be compared with those obtained in [34].

Zoom In Zoom Out Reset image size

Ion tori are geometrically thick, optically thin accretion structures with constant specific angular momentum. They are the optically thin counterparts of the optically thick Polish doughnuts [35, 36].

For a torus made of a perfect fluid, it can be shown [35] that the energy–momentum conservation law along with the assumption of constant specific angular momentum leads to

$$\begin{equation} \frac{\nabla _{\mu } p}{p+\epsilon } = \nabla _{\mu } \mathcal {W}, \end{equation} \tag{ 25 }$$

where p is the pressure, is the proper energy density of the fluid and the potential $\mathcal {W}$ is related to the covariant component u_t of the fluid 4-velocity by $\mathcal {W} = -\ln |u_{t}|$. The isobaric surfaces thus coincide with the equipotential surfaces of $\mathcal {W}$.

The cross-sectional shape of the equipotential surfaces is given in figure 1 of [35]. It appears that one particular surface has a cusp at some r = r_crit: it crosses itself in the equatorial plane. Equipotential surfaces contained inside this critical surface are not connected to the central object; thus, matter cannot be accreted and swallowed by the black hole. It is thus assumed that the torus physical surface coincides with this critical surface. The central point, r = r_central, coincides with the innermost equipotential surface and to the point of maximum pressure.

It is convenient to introduce the dimensionless parameter

$$\begin{equation} w=\frac{\mathcal {W}-\mathcal {W}_{\mathrm{crit}}}{\mathcal {W}_{\mathrm{central}}-\mathcal {W}_{\mathrm{crit}}}. \end{equation} \tag{ 26 }$$

The potential (and the pressure) increases continuously between r_crit and r_central. Thus, w varies from 0 (cusp) to 1 (center) inside the torus. Outside, w can take any other value (positive or negative).

The physics of the radiative transfer used for the ion torus is based on [37] and is presented in detail in [38]. Figure 5 shows the image of an ion torus around a Kerr black hole with a spin parameter a = 0.5 M, with trivial radiative transfer included. This image can be compared with figure 5 of [39]. Figure 6 shows the image of such an ion torus at the Galactic center as seen by an observer on Earth, with radiative transfer including the emission of synchrotron radiation (see [38] for the details of this computation).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 7 shows the synchrotron emitted spectrum of the previous ion torus. This is an example of the computation of the spectrum of an optically thin object by GYOTO. Astrophysically interesting perspectives of this kind of computation will be presented in [38].

Zoom In Zoom Out Reset image size

From the 3+1 data (N, βⁱ, γ_ij), reconstruct the 4-metric g_αβ according to (27); then, compute the Christoffel symbols ⁴Γ^α_μν of g_αβ with respect to the coordinates (x^α) and integrate the geodesic equation in the standard four-dimensional form:

$$\begin{equation} \frac{{\mathrm{d}}^2 X^\alpha }{{\mathrm{d}}\lambda ^2} + {}^4 \Gamma ^\alpha _{\ \, \mu \nu } \frac{{\mathrm{d}}X^\mu }{{\mathrm{d}}\lambda } \frac{{\mathrm{d}}X^\nu }{{\mathrm{d}}\lambda } = 0 , \end{equation} \tag{ 28 }$$

where λ is either the proper time (timelike geodesics) or some affine parameter (null geodesics), the geodesic being defined by x^α = X^α(λ).

The geodesic is searched in the form

$$\begin{equation} x^i = X^i(t) , \end{equation} \tag{ 29 }$$

i.e. the geodesic is parametrized by the coordinate time t and not by λ. In [43], we have derived a system of differential equations for the functions Xⁱ(t), which involve only the 3+1 quantities (N, βⁱ, γ_ij, K_ij) and their spatial derivatives (denoted by ∂_i ≡ ∂/∂xⁱ) :

$$\begin{eqnarray} &&\fl \left\{\begin{array}{@{}l@{}} \dsty\frac{{\mathrm{d}}X^i}{{\mathrm{d}}t} = N V^i - \beta ^i \hspace{270pt} \,(30a)\nonumber\\ \dsty\frac{{\mathrm{d}}V^i}{{\mathrm{d}}t} = N[ V^i( V^j\partial _j\ln N - K_{jk} V^j V^k ) + 2 K^i_{\ \, j} V^j - {}^3\Gamma ^i_{\ jk} V^j V^k ]- \gamma ^{ij}\partial _j N - V^j \partial _j \beta ^i.\ \ \ \ \ \ \,(30b) \end{array} \right. \end{eqnarray} \tag{ 30 }$$

Here, the ³Γⁱ_jk are the Christoffel symbols of the 3-metric γ_ij with respect to the coordinates (xⁱ). The auxiliary variable Vⁱ(t) = N⁻¹(dXⁱ/dt + βⁱ) is the linear momentum of the considered particle divided by its energy, both quantities being measured by the Eulerian observer (i.e. the observer whose worldline is normal to the Σ_t hypersurfaces) [43].