YNOGK: A NEW PUBLIC CODE FOR CALCULATING NULL GEODESICS IN THE KERR SPACETIME

From the equations of motion, we know that both R and Θ_θ must be nonnegative. This restriction divides the coordinate space into allowed (where R ⩾ 0 and Θ_θ ⩾ 0) and forbidden (where R < 0 or Θ_θ < 0) regions for the motion of a photon. The boundary points of these regions are the so-called turning points, and their coordinates r_tp and θ_tp satisfy equations R(r_tp) = 0 and Θ_θ(θ_tp) = 0. For a photon emitted at r_ini and θ_ini, its motion is confined between two turning points $r_{{\rm tp_1}}$ and $r_{{\rm tp_2}}$ for a radial coordinate, and $\theta _{{\rm tp_1}}$ and $\theta _{{\rm tp_2}}$ for a poloidal coordinate. If we assume $r_{{\rm tp_1}}\le r_{{\rm tp_2}}$ and $\theta _{{\rm tp_1}}\le \theta _{{\rm tp_2}}$, then we have $r_{{\rm ini}}\in [r_{{\rm tp_1}}, r_{{\rm tp_2}}]$ and $\theta _{{\rm ini}}\in [\theta _{{\rm tp_1}}, \theta _{{\rm tp_2}}]$. Because $p_r=\pm \sqrt{R}/\Delta$ and $p_\theta =\pm \sqrt{\Theta _\theta }$, if p_r = 0 (or p_θ = 0) at the initial position, we have R(r_ini) = 0 (or Θ_θ(θ_ini) = 0); therefore, r_ini (or θ_ini) must be a turning point and equal to one of $r_{{\rm tp_1}}$ and $r_{{\rm tp_2}}$ (or $\theta _{{\rm tp_1}}$ and $\theta _{{\rm tp_2}}$).

The radial motion of a photon can be unbounded, meaning that the photon can go to infinity or fall into the black hole. These cases usually correspond to the equation R(r) = 0 which has no real roots or $r_{{\rm tp_1}}$ is less than or equal to the radius of the event horizon. We regard infinity and the event horizon of the black hole as two special turning points in the radial motion. A photon asymptotically approaches them but never returns from them. Thus, $r_{{\rm tp_2}}$ can be infinity and $r_{{\rm tp_1}}$ can be less than or equal to r_h (here r_h is the radius of the event horizon).

For the poloidal motion, there are also two special positions: θ = 0 and θ = π, i.e., the spin axis of the black hole. A photon with λ = 0 goes through the spin axis due to zero angular momentum and changes the sign of its angular velocity dθ/dσ instantaneously, and its azimuthal coordinate jumps from ϕ to ϕ ± π (Shakura 1987), implying that the spin axis is not a turning point. From Equation (13), we also know that 0 and π are not the roots of the equation Θ_θ(θ) = 0.

First, we use a new variable μ to replace cos θ, and Equation (23) can be rewritten as

$$\begin{equation} p=\pm \int ^r\frac{dr}{\sqrt{R}}=\pm \int ^\mu \frac{d\mu }{\sqrt{\Theta _\mu }}, \end{equation} \tag{ 24 }$$

where

$$\begin{equation} \Theta _\mu =q-(q+\lambda ^2-a^2)\mu ^2-a^2\mu ^4. \end{equation} \tag{ 25 }$$

Both R and Θ_μ are quartic, but the polynomial of Weierstrass's standard elliptic integral is cubic. We need a variable transformation to make R and Θ_μ cubic. We define the following constants for poloidal motion:

$$\begin{eqnarray} &&b_0 = -4a^2\mu _{{\rm tp_1}}^3-2(q+\lambda ^2-a^2)\mu _{{\rm tp_1}}, \end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray} &&b_1 = -2a^2\mu _{{\rm tp_1}}^2-\frac{1}{3}(q+\lambda ^2-a^2), \end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray} &&b_2 = -\frac{4}{3}a^2\mu _{{\rm tp_1}}, \end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray} &&b_3 = -a^2, \end{eqnarray} \tag{ 29 }$$

where $\mu _{{\rm tp_1}}=\cos \theta _{{\rm tp_1}}$ and introduce a new variable t,

$$\begin{equation} t=\frac{b_0}{4}\frac{1}{(\mu -\mu _{{\rm tp_1}})}+\frac{b_1}{4}. \end{equation} \tag{ 30 }$$

Making the transformation from μ to t, the μ part of Equation (24) can be reduced to

$$\begin{equation} p = \pm \int ^{t(\mu)}\frac{dt}{\sqrt{4t^3-g_2t-g_3}}, \end{equation} \tag{ 31 }$$

where $g_2 = ({3}/{4})(b_1^2-b_0b_2)$, $g_3 = ({1}/{16})(3b_0b_1b_2-2b_1^3-b^2_0b_3).$ Using the definition of Weierstrass's elliptic function ℘(z; g₂, g₃) (Abramowitz & Stegun 1965) from Equation (31), we have t = ℘(p ± Π_μ; g₂, g₃). Solving Equation (30) for μ, we can express μ as the function of p:

$$\begin{equation} \mu (p)=\frac{b_0}{4\wp (p\pm \Pi _\mu;g_2,g_3)-b_1}+\mu _{{\rm tp_1}}, \end{equation} \tag{ 32 }$$

where

$$\begin{equation} \Pi _\mu =|\wp ^{-1}[t(\mu _{{\rm ini}});g_2,g_3]|. \end{equation} \tag{ 33 }$$

The sign in front of Π_μ depends on the initial value of p_θ, which is the θ component of the four-momentum of a photon, and

$$\begin{eqnarray} \left\lbrace \begin{array}{@{}ll}p_\theta >0,&\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \;\; +,\\ p_\theta =0,&\left\lbrace \begin{array}{@{}lll}\theta _{{\rm ini}} = \theta _{{\rm tp_1}},& \Pi _\mu =n\omega,&+,-,\\ \theta _{{\rm ini}} = \theta _{{\rm tp_2}},& \Pi _\mu =\left(\frac{1}{2}+n\right)\omega,&+,-,\\ \end{array} \right. \cr p_\theta <0,&\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad -, \end{array} \right. \end{eqnarray} \tag{ 34 }$$

where ω is the period of ℘(z; g₂, g₃) and n = 0, 1, 2, .... The sign in front of Π_μ can be "+" or "−" when p_θ = 0.

From the above discussion, we know that one root of equation Θ_μ = 0 is needed in the variable transformation, namely, $\mu _{{\rm tp_1}}$. To avoid the complexity caused by introducing complex numbers, we always use the real root. Luckily, the equation Θ_μ = 0 always has real roots, but the same is not true for the equation R(r) = 0. For cases in which the equation R(r) = 0 has no real roots, we use Jacobi's elliptic functions sn(z|k²), cn(z|k²) to express r.

If the equation R(r) = 0 has real roots, then $r_{{\rm tp_1}}$ exists. We can define the following constants using $r_{{\rm tp_1}}$:

$$\begin{eqnarray} &&b_0 = 4r_{{\rm tp_1}}^3-2(q+\lambda ^2-a^2)r_{{\rm tp_1}}+2[q+(\lambda -a)^2], \end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray} &&b_1 = 2r_{{\rm tp_1}}^2-\frac{1}{3}(q+\lambda ^2-a^2), \end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray} &&b_2 = \frac{4}{3}r_{{\rm tp_1}}, \end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray} &&b_3 = 1, \end{eqnarray} \tag{ 38 }$$

and introducing a new variable t,

$$\begin{equation} t=\frac{b_0}{4}\frac{1}{(r-r_{{\rm tp_1}})}+\frac{b_1}{4}. \end{equation} \tag{ 39 }$$

It is similar to μ, using t as an independent variable. We can reduce the r part of Equation (24) to the standard form of Weierstrass's elliptical integral

$$\begin{equation} p = \pm \int ^{t(r)}\frac{dt}{\sqrt{4t^3-g_2t-g_3}}, \end{equation} \tag{ 40 }$$

where $g_2 = ({3}/{4})(b_1^2-b_0b_2)$, $g_3 = ({1}/{16})(3b_0b_1b_2-2b_1^3-b^2_0b_3).$ Taking the inverse of the above equation, we get t = ℘(p ± Π_r; g₂, g₃). Solving Equation (39) for r, we have

$$\begin{equation} r(p)=\frac{b_0}{4\wp (p\pm \Pi _r;g_2,g_3)-b_1}+r_{{\rm tp_1}}, \end{equation} \tag{ 41 }$$

where

$$\begin{equation} \Pi _r=|\wp ^{-1}[t(r_{{\rm ini}});g_2,g_3]|.\quad \end{equation} \tag{ 42 }$$

The sign in front of Π_r also depends on the initial value of p_r, which is the r component of the four-momentum of a photon, and

$$\begin{eqnarray} \left\lbrace \begin{array}{@{}ll}p_r>0,&\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad +,\\ p_r=0,&\left\lbrace \begin{array}{@{}lll} r_{{\rm ini}} = r_{{\rm tp_1}},& \Pi _r=n\omega,&+,-,\\ r_{{\rm ini}} = r_{{\rm tp_2}},& \Pi _r=\left(\frac{1}{2}+n\right)\omega,&+,-,\\ \end{array} \right. \cr p_r<0,&\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad -, \end{array} \right. \end{eqnarray} \tag{ 43 }$$

where ω is the period of ℘(z; g₂, g₃) and n = 0, 1, 2, .... The sign in front of Π_r can be "+" or "−" when p_r = 0.

If the equation R(r) = 0 has no real roots, we use Jacobi's elliptic functions to express r. Since the coefficient of r³ is zero, the roots of equation R(r) = 0 satisfy r₁ + r₂ + r₃ + r₄ = 0. Therefore, the roots r₁, r₂, r₃, and r₄ can be written as

$$\begin{eqnarray} \nonumber &&r_1=u-iw, \,\,\, r_2=u+iw,\\ &&r_3=-u-iv, r_4=-u+iv. \end{eqnarray} \tag{ 44 }$$

Introducing two constants λ₁ and λ₂

$$\begin{equation} \lambda _{1,2}=\frac{1}{2w^2}[4u^2+v^2+w^2\pm \sqrt{(4u^2+w^2+v^2)^2-4w^2v^2}], \end{equation} \tag{ 45 }$$

which satisfy λ₁ > 1 > λ₂ > 0, and a new variable t,

$$\begin{equation} t=\sqrt{\frac{\lambda _1-1}{(\lambda _1-\lambda _2)[(r-u)^2+w^2]}}\left(r-u\frac{\lambda _1+1}{\lambda _1-1}\right), \end{equation} \tag{ 46 }$$

we can reduce the r part of Equation (24) to Legendre's standard elliptic integral

$$\begin{equation} p=\int ^t\frac{dt}{w\sqrt{\lambda _1}\sqrt{(1-t^2)(1-m_2t^2)}}, \end{equation} \tag{ 47 }$$

where

$$\begin{equation} m_2=\frac{\lambda _1-\lambda _2}{\lambda _1}. \end{equation} \tag{ 48 }$$

Using the definition of Jacobi's elliptic function sn(z|k²) (Abramowitz & Stegun 1965), from Equation (47) we obtain $t=\mathrm{sn}(pw\sqrt{\lambda _1}\pm \Pi _0|m_2)$. Solving Equation (46) for r, we get the expression for r as a function of p

$$\begin{equation} r_{\pm }(p)=u+\frac{-2u\pm w(\lambda _1-\lambda _2)\mathrm{sn}(pw\sqrt{\lambda _1}\pm \Pi _0|m_2)|\mathrm{cn}(pw\sqrt{\lambda _1}\pm \Pi _0|m_2)|}{(\lambda _1-\lambda _2)\mathrm{sn}^2(pw\sqrt{\lambda _1}\pm \Pi _0|m_2)-(\lambda _1-1)}, \end{equation} \tag{ 49 }$$

where

$$\begin{equation} \Pi _0=|{\rm sn}^{-1}\left[t(r_{{\rm ini}})|m_2\right]|. \end{equation} \tag{ 50 }$$

When the initial value of p_r > 0, we have r = r₋, and when p_r < 0, we have r = r₊. p_r cannot be zero, otherwise the initial radial coordinate r_ini of the photon becomes one root of the equation R(r) = 0, which is the case that has been discussed above.

In this section, we express the coordinates t, ϕ and the affine parameter σ as the numerical functions of the parameter p. All of these variables have been expressed as the integrals of r and θ in Equations (9)–(11) and (17)–(22). We achieve our goal if we can compute all of these integrals along a geodesic for a specified p. Making transformations from r and μ to a new variable t (defined by Equations (30), (39), and (46)), we compute these integrals under the new variable t. For simplicity, we use F_r(t) and F_θ(t) to denote the complicated integrands (see below) for r and θ, respectively.

First, we discuss the integral path, which starts from the initial position and terminates at the photon. If the photon encounters the turning points along the geodesic, then the whole integral path is not monotonic, as shown in Figure 1 for poloidal motion (radial motion is similar). In this figure, the projected poloidal motion of a photon in the r–θ plane is illustrated. The motion is confined between two turning points: $\mu _{{\rm tp_1}}$ and $\mu _{{\rm tp_2}}$. The photon encounters the turning points three times. Obviously any section of the path that contains one or more than one turning point is not monotonic, such as paths CDE, EFP, etc. The path between any two neighboring turning points has the maximum monotonic length and the total integrals should be computed along each of them and summed.

Zoom In Zoom Out Reset image size

There are four important points involved in the limits of these integrals, i.e., $\mu _{{\rm tp_1}}$, $\mu _{{\rm tp_2}}$, μ_ini, and μ_p. These points are the μ coordinates of turning points, the initial point, and the photon position for a given p, respectively. The values of these points corresponding to the new variable t are $t_{{\rm tp_1}}$, $t_{{\rm tp_2}}$, t_ini, which can be calculated from Equation (30), and t_p, which can be calculated from t = ℘(p ± Π_μ; g₂, g₃) with a given p. Because the function t = t(μ) expressed by Equation (30) is monotonically increasing, we have $t_{{\rm tp_2}}\le t_{{\rm ini}}\le t_{{\rm tp_1}}$ and $t_{{\rm tp_2}}\le t_{p}\le t_{{\rm tp_1}}$.

We use Nt₁ and Nt₂ to denote the number of times a photon meets the turning points $\mu _{{\rm tp_1}}$ and $\mu _{{\rm tp_2}}$ for a given p, respectively, and define the following integrals (cf. Figure 1):

$$\begin{eqnarray} \nonumber &&I_0=\int ^{t_p}_{t_{{\rm ini}}}F_\theta (t) dt,\quad I_1=\int ^{t_{{\rm tp_1}}}_{t_{p}}F_\theta (t) dt,\quad I_2=\int ^{t_p}_{t_{{\rm tp_2}}}F_\theta (t) dt,\\ &&I_{01}=I_0+I_1=\int ^{t_{{\rm tp_1}}}_{t_{{\rm ini}}}F_\theta (t) dt,\quad I_{02}=I_2-I_0=\int ^{t_{{\rm ini}}}_{t_{{\rm tp_2}}}F_\theta (t) dt. \end{eqnarray} \tag{ 51 }$$

The integrals of θ in σ, t, and ϕ can be written as (cf. Figure 1)

$$\begin{eqnarray} \nonumber \sigma _\theta \,\,\,(\,\,\phi _\theta,\,\,\,t_\theta \,\,)&=& -\mathrm{sign}(p_\theta)I_0+2Nt_1I_1+2Nt_2I_2,\\ &=&-[\mathrm{sign}(p_\theta)+2Nt_1-2Nt_2]I_0+2Nt_1I_{01}+2Nt_2I_{02}, \end{eqnarray} \tag{ 52 }$$

where p_θ is the θ component of the initial four-momentum of a photon. In order to evaluate the above expression, we need to know Nt₁ and Nt₂ for a given p. Similarly, if we define five integrals from Equation (31) as follows:

$$\begin{eqnarray} && p_0=\int ^{t_p}_{t_{\mathrm{ini}}}\frac{dt}{\sqrt{W(t)}},\quad p_1=\int ^{t_{{\rm tp_1}}}_{t_p}\frac{dt}{\sqrt{W(t)}},\quad p_2=\int ^{t_p}_{t_{{\rm tp_2}}}\frac{dt}{\sqrt{W(t)}}, \nonumber\\ \ms && p_{01}=p_0+p_1=\int ^{t_{{\rm tp_1}}}_{t_{\mathrm{ini}}}\frac{dt}{\sqrt{W(t)}},\quad p_{02} = p_2-p_0=\int ^{t_{\mathrm{ini}}}_{t_{{\rm tp_2}}}\frac{dt}{\sqrt{W(t)}}, \end{eqnarray} \tag{ 53 }$$

where W(t) = 4t³ − g₂t − g₃, we get the following identity:

$$\begin{eqnarray} \nonumber p&=&-\mathrm{sign}(p_\theta)p_0+2Nt_1p_1+2Nt_2p_2,\\ &=&-[\mathrm{sign}(p_\theta)+2Nt_1-2Nt_2]p_0+2Nt_1p_{01}+2Nt_2p_{02}. \end{eqnarray} \tag{ 54 }$$

Note that Nt₁ and Nt₂ are not arbitrary and are instead related to the initial direction of the photon in poloidal motion. For p_θ > 0 (or p_θ = 0 and $\theta _{{\rm ini}}=\theta _{{\rm tp_1}}$), they increase as

$$\begin{eqnarray*} &&Nt_1= 0\quad 0\quad 1\quad 1\quad 2\quad 2\ldots,\\ &&Nt_2= 0\quad 1\quad 1\quad 2\quad 2\quad 3\ldots. \end{eqnarray*}$$

For p_θ < 0 (or p_θ = 0 and $\theta _{{\rm ini}}=\theta _{{\rm tp_2}}$), they increase as

$$\begin{eqnarray*} &&Nt_1= 0\quad 1\quad 1\quad 2\quad 2\quad 3\ldots,\\ &&Nt_2= 0\quad 0\quad 1\quad 1\quad 2\quad 2\ldots. \end{eqnarray*}$$

For a given p, we find that one pair of Nt₁ and Nt₂ values that satisfies Equation (54) always exists. These values are the number of photons meeting the turning points. With Nt₁ and Nt₂, Equation (52) now can be evaluated readily.

For the r coordinate, the process is similar to that above. Nt₁ and Nt₂ also represent the number of times a photon meets the turning points $r_{{\rm tp_1}}$ and $r_{{\rm tp_2}}$, respectively. Five integrals are defined as

$$\begin{eqnarray} &&\nonumber I_0=\int _{t_{{\rm ini}}}^{t_p}F_r(t) dt,\quad I_1=\int ^{t_{{\rm tp_1}}}_{t_p}F_r(t) dt, \quad I_2 = \int ^{t_{p}}_{t_{{\rm tp_2}}} F_r(t) dt,\\ \ms && I_{01}=I_0+I_1=\int ^{t_{{\rm tp_1}}}_{t_{{\rm ini}}}F_r(t)dt, \quad I_{02}=I_2-I_0=\int ^{t_{{\rm ini}}}_{t_{{\rm tp_2}}}F_r(t)dt, \end{eqnarray} \tag{ 55 }$$

where t_p = ℘(p ± Π_r; g₂, g₃) (or $t_p=\mathrm{sn}(pw\sqrt{\lambda _1}\pm \Pi _0|m_2)$ when the equation R(r) = 0 has no real roots), and $t_{p_1}$, $t_{p_2}$, and t_ini are calculated from Equation (39) or Equation (46). Then the integrals of r in σ, t, ϕ can be written as

$$\begin{eqnarray} \nonumber \sigma _r\,\,\,(\,\,t_r,\,\, \phi _r\,\,)&=&-\mathrm{sign}(p_r)I_0+2Nt_1I_1+2Nt_2I_2,\\ &=&-[\mathrm{sign}(p_r)+2Nt_1-2Nt_2]I_0+2Nt_1I_{01}+2Nt_2I_{02}. \end{eqnarray} \tag{ 56 }$$

To get Nt₁ and Nt₂, we define p₀, p₁, and p₂ from Equation (40) as

$$\begin{eqnarray} &&\nonumber p_0=\int ^{t_p}_{t_{{\rm ini}}}\frac{dt}{\sqrt{W(t)}},\quad p_1=\int ^{t_{p_1}}_{t_p}\frac{dt}{\sqrt{W(t)}},\quad p_2=\int ^{t_{p}}_{t_{p_2}}\frac{dt}{\sqrt{W(t)}},\\ \ms &&p_{01}=p_0+p_1=\int ^{t_{{\rm tp_1}}}_{t_{{\rm ini}}}\frac{dt}{\sqrt{W(t)}},\quad p_{02}=p_2-p_0=\int ^{t_{{\rm ini}}}_{t_{{\rm tp_2}}}\frac{dt}{\sqrt{W(t)}}. \end{eqnarray} \tag{ 57 }$$

For a given p, we have

$$\begin{eqnarray} \nonumber p&=&-\mathrm{sign}(p_r)p_0+2Nt_1p_1+2Nt_2p_2,\\ &=&-[\mathrm{sign}(p_r)+2Nt_1-2Nt_2]p_0+2Nt_1p_{01}+2Nt_2p_{02}. \end{eqnarray} \tag{ 58 }$$

To get Nt₁ and Nt₂ from the above equation, one just needs to note that when p_r > 0 (or p_r = 0 and $r_{{\rm ini}}=r_{{\rm tp_1}}$), they increase as

$$\begin{eqnarray*} &&Nt_1= 0\quad 0\quad 1\quad 1\quad 2\quad 2\ldots,\\ &&Nt_2= 0\quad 1\quad 1\quad 2\quad 2\quad 3\ldots, \end{eqnarray*}$$

when p_r < 0 (or p_r = 0 and $r_{{\rm ini}}=r_{{\rm tp_2}}$), they increase as

$$\begin{eqnarray*} &&Nt_1= 0\quad 1\quad 1\quad 2\quad 2\quad 3\ldots,\\ &&Nt_2= 0\quad 0\quad 1\quad 1\quad 2\quad 2\ldots. \end{eqnarray*}$$

Similarly, for a given p, there is one pair of Nt₁ and Nt₂ values that satisfies Equation (58). With Nt₁ and Nt₂ values, Equation (56) can now be evaluated. In many cases, the number of times a photon meets the turning point in r is fewer than 2. Specifically, the number becomes zero when the equation R(r) = 0 has no real roots.

In previous sections, the four coordinates r, θ, ϕ, t and the affine parameter σ have been expressed as functions of p; many elliptic integrals need to be calculated using this routine. In this section, we reduce these integrals to standard forms and then evaluate them using Carlson's method as Dexter & Agol (2009) did.

First, we introduce two notations J_k(h) and I_k(h) with the following definitions:

$$\begin{eqnarray} &&J_k(h) = \int ^{t_2}_{t_1}\frac{dt}{(t-h)^k\sqrt{4t^3-g_2t-g_3}}, \end{eqnarray} \tag{ 59 }$$

$$\begin{eqnarray} &&I_k(h) = \int ^{r_2}_{r_1}\frac{dr}{(r-h)^k\sqrt{(r^2-2ur+u^2+w^2)(r^2+2ur+u^2+v^2)}}, \end{eqnarray} \tag{ 60 }$$

where k is an integer. From Equation (8), we get one of the standard forms as

$$\begin{equation} J_0=\int ^{t_2}_{t_1}\frac{dt}{\sqrt{4t^3-g_2t-g_3}}. \end{equation} \tag{ 61 }$$

After having been reduced to J₀, the forms of the integrals of r and θ in Equation (8) are exactly the same. Noting the definition of the parameter p, we have J₀ = p. The radial integrals in Equation (9) are reduced to

$$\begin{eqnarray} \sigma _r&=&\frac{b_0^2}{16}J_2\left(\frac{b_1}{4}\right)+\frac{b_0r_{{\rm tp_1}}}{2}J_1\left(\frac{b_1}{4}\right)+r^2_{{\rm tp_1}}p, \end{eqnarray} \tag{ 62 }$$

where J₀ has been replaced by p. The radial integrals in Equation (10) can be reduced to

$$\begin{eqnarray} &&t_r=\sigma _r+\frac{b_0}{2}J_1\left(\frac{b_1}{4}\right)+(2r_{{\rm tp_1}}+4+A_{t+}-A_{t-})p-B_{t+}J_1(t_+)+B_{t-}J_1(t_-), \end{eqnarray} \tag{ 63 }$$

where

$$\begin{eqnarray} \nonumber r_\pm &=&1\pm \sqrt{1-a^2},\\ \nonumber A_{t\pm }&=&\frac{2[r_\pm (4-a\lambda)-2a^2]}{(r_+-r_-)(r_{{\rm tp_1}}-r_{\pm })},\\ \nonumber B_{t\pm }&=&\frac{[r_\pm (4-a\lambda)-2a^2]b_0}{2(r_+-r_-)(r_{{\rm tp_1}}-r_{\pm })^2},\\ t_\pm &=&\frac{b_1}{4}+\frac{b_0}{4(r_\pm -r_{{\rm tp_1}})}. \end{eqnarray} \tag{ 64 }$$

Similarly, the radial integrals in Equation (11) have the following form:

$$\begin{equation} \phi _r=a[(A_{\phi +}-A_{\phi -})p-B_{\phi +}J_1(t_+)+B_{\phi -}J_1(t_-)], \end{equation} \tag{ 65 }$$

where

$$\begin{eqnarray} \nonumber A_{\phi \pm }&=&\frac{2r_\pm -a\lambda }{(r_+-r_-)(r_{{\rm tp_1}}-r_\pm)},\\ B_{\phi \pm }&=&\frac{(2r_\pm -a\lambda)b_0}{4(r_+-r_-)(r_{{\rm tp_1}}-r_\pm)^2}. \end{eqnarray} \tag{ 66 }$$

When equation R(r) = 0 has no real roots, the integrals of r in σ, t, and ϕ can be written as

$$\begin{eqnarray} &&\sigma _r = I_{-2}(0), \end{eqnarray} \tag{ 67 }$$

$$\begin{eqnarray} &&t_r = \sigma _r+4p+2I_{-1}(0)+C_{t+}I_{1}(r_+)-C_{t-}I_1(r_-), \end{eqnarray} \tag{ 68 }$$

$$\begin{eqnarray} &&\phi _r = a[C_{\phi +}I_1(r_+)-C_{\phi -}I_1(r_-)], \end{eqnarray} \tag{ 69 }$$

where

$$\begin{eqnarray} &&C_{t\pm } = \frac{2[r_\pm (4-a\lambda)-2a^2]}{r_+-r_-},\quad C_{\phi \pm } = \frac{2r_\pm -a\lambda }{r_+-r_-}. \end{eqnarray} \tag{ 70 }$$

The integrals concerning μ in Equation (9) are reduced to

$$\begin{equation} \sigma _\mu =a^2\left[\frac{b_0^2}{16}J_2\left(\frac{b_1}{4}\right)+\frac{b_0\mu _{{\rm tp_1}}}{2}J_1\left(\frac{b_1}{4}\right)+\mu ^2_{{\rm tp_1}}p\right], \end{equation} \tag{ 71 }$$

and t_μ = σ_μ. The μ integrals in Equation (11) can be reduced to

$$\begin{equation} \phi _\mu =\lambda \left[\frac{p}{1-\mu _{{\rm tp_1}}^2}+W_{\mu _+}J_1(t_+)-W_{\mu _-}J_1(t_-)\right], \end{equation} \tag{ 72 }$$

where

$$\begin{eqnarray} \nonumber W_{\mu _\pm }&=&\frac{b_0}{8(-1\pm \mu _{{\rm tp_1}})^2},\\\ms t_\pm &=&\frac{b_1}{4}+\frac{b_0}{4(\pm 1-\mu _{{\rm tp_1}})}. \end{eqnarray} \tag{ 73 }$$

Finally, we have

$$\begin{eqnarray} \nonumber \sigma &= &\sigma _r+\sigma _\mu,\\ \nonumber t&=&t_r+t_\mu,\\ \phi &=&\phi _r+\phi _\mu. \end{eqnarray} \tag{ 74 }$$

From Equations (61)–(72) the integrals that need to be calculated are J₀, J₁, J₂, and I₁, I₋₁, I₋₂. Now, we use Carlson's method to evaluate them. When equation 4t³ − g₂t − g₃ = 0 has three real roots denoted by e₁, e₂, and e₃, J_k(h) can be written as (Carlson 1988)

$$\begin{eqnarray} \nonumber J_k(h)&=&s_h\frac{1}{2}\int ^x_y\frac{dt}{\sqrt{(t-e_1)(t-e_2)(t-e_3)(t-h)^{2k}}},\\ &=&s_h\frac{1}{2}[-1,-1,-1,-2k], \end{eqnarray} \tag{ 75 }$$

where s_h = sign[(y − h)^k]. When equation 4t³ − g₂t − g₃ = 0 has one real root e₁ and one pair of complex conjugate roots u+iv and u−iv, J_k(h) can be written as (Carlson 1991)

$$\begin{eqnarray} \nonumber J_k(h)&=&s_h\frac{1}{2}\int ^x_y\frac{dt}{\sqrt{(t-e_1)(t^2-2ut+u^2+v^2)(t-h)^{2k}}},\\ &=&s_h\frac{1}{2}[-1,-1,-1,-2k]. \end{eqnarray} \tag{ 76 }$$

From Equations (30) and (39), we know that when $r = r_{{\rm tp_1}}$ or $\mu = \mu _{{\rm tp_1}}$, t becomes ∞; thus, one limit of these integrals can be infinity.

The integrals I_k(h) can be reduced to Carlson's integrals directly (Carlson 1992)

$$\begin{eqnarray} \nonumber I_{k}(h) &=&s_h\int ^x_y\frac{dr}{\sqrt{(r^2-2ur+u^2+w^2)(r^2+2ur+u^2+v^2)(r-h)^{2k}}},\\ &=&s_h[-1,-1,-1,-1,-2k], \end{eqnarray} \tag{ 77 }$$

where [p₁, ..., p_k] is a symbol used by Carlson (1992) to denote the elliptic integrals and has the following definition:

$$\begin{eqnarray} &&[p_1,\ldots,p_k] = \int ^x_y \prod _{i = 1}^{k}(a_i+b_it)^{p_i/2}dt, \end{eqnarray} \tag{ 78 }$$

which can be evaluated by the formulae provided in Carlson (1988, 1989, 1991, 1992).

In previous sections, we have expressed all coordinates as functions of a parameter p and discussed how to calculate them using Carlson's method, but before the calculation, one needs to specify the constants of motion and p_r, p_θ, which determine the signs in front of Π_μ, Π_r, Π₀ and also how the number of turning points increases. In this section, we discuss how to compute λ and q and p_r, p_θ from $\bar{p}_{(a)}$, which are the components of the four-momentum measured in the locally nonrotating reference frame (LNRF) and have been specified by the user.

First, following Bardeen et al. (1972), we introduce the LNRF observers or the zero angular momentum observer; the basis vectors of the orthonormal tetrad of the observers are given by

$$\begin{equation} \mathbf {e}_{(a)}(\mathrm{LNRF})=e_{{(a)}}^{\nu }\partial _{\nu }, \end{equation} \tag{ 79 }$$

where

$$\begin{eqnarray} e_{{(a)}}^{\nu }=\left(\begin{array}{@{}cccc@{}}e^{-\nu }&0 &0 & \omega e^{-\nu }\\ 0 &e^{-\mu _1} &0 & 0\\ 0&0& e^{-\mu _2} &0\\ 0 & 0 & 0 & e^{-\psi } \end{array}\right). \end{eqnarray} \tag{ 80 }$$

The covariant components of the four-momentum of a photon in B–L coordinates can be expressed as

$$\begin{eqnarray} &&p_\mu =E\left(-1,s_r\frac{\sqrt{R}}{\Delta },s_\theta \sqrt{\Theta _{\theta }},\lambda \right), \end{eqnarray} \tag{ 81 }$$

where s_r and s_θ are the signs of the r and θ components. One can easily show that $\bar{p}_{(a)} = e_{(a)}^\mu p_\mu$; namely (Shakura 1987)

$$\begin{eqnarray} \bar{p}_{(t)}&=&-Ee^{-\nu }(1-\lambda \omega), \end{eqnarray} \tag{ 82 }$$

$$\begin{eqnarray} \bar{p}_{(r)}&=&s_r Ee^{-\mu _1}\frac{\sqrt{R}}{\Delta }, \end{eqnarray} \tag{ 83 }$$

$$\begin{eqnarray} \bar{p}_{(\theta)}&=&s_\theta Ee^{-\mu _2} \sqrt{\Theta _\theta }, \end{eqnarray} \tag{ 84 }$$

$$\begin{eqnarray} \bar{p}_{(\phi)}&=&E\lambda e^{-\psi }. \end{eqnarray} \tag{ 85 }$$

From Equations (83) and (84) we have $s_r=\mathrm{sign}(\bar{p}_{(r)})$ and $s_\theta =\mathrm{sign}(\bar{p}_{(\theta)})$, which determine the initial direction of the photon in the B–L system (cf. Equations (34) and (43)), thus determining how the number of turning points increases.

Solving Equations (82) and (85) simultaneously for λ, one obtains

$$\begin{eqnarray} &&\lambda = \frac{\sin \theta \bar{p}_{(\phi)}/\bar{p}_{(t)}}{-\sqrt{\Delta }\Sigma /A+\omega \sin \theta \bar{p}_{(\phi)}/\bar{p}_{(t)}}. \end{eqnarray} \tag{ 86 }$$

Using λ and Equation (82), one obtains $E=\bar{p}_{(t)}e^\nu /(1-\lambda \omega)$. Using λ and E from Equation (84), one obtains the formula for calculating the motion constant q,

$$\begin{eqnarray} &&q=\left[\left(\frac{\bar{p}_{(\phi)}/\bar{p}_{(t)}}{-\sqrt{\Delta }\Sigma /A+\omega \sin \theta \bar{p}_{(\phi)}/\bar{p}_{(t)}}\right)^2-a^2\right]\cos ^2\theta +\left[\frac{\bar{p}_{(\theta)}}{\bar{p}_{(t)}}(1-\lambda \omega)\right]^2\frac{A}{\Delta }. \end{eqnarray} \tag{ 87 }$$

Thus, we have obtained the basic Equations (86) and (87) connecting λ, q and the components of the four-momentum $\bar{p}_{(a)}$ measured in the LNRF.

To prescribe $\bar{p}_{(a)}$, one should note that it satisfies the following equation:

$$\begin{eqnarray} &&-\bar{p}_{(t)}^2+\bar{p}_{(r)}^2+\bar{p}_{(\theta)}^2+\bar{p}_{(\phi)}^2=0. \end{eqnarray} \tag{ 88 }$$

Thus there are only three independent components. Obviously the user can specify the four-momentum $\bar{p}_{(a)}$ directly in the LNRF or equivalently specify $p^{\prime }_{(a)}$ in any other reference frame of his/her own choice and then transform it into the LNRF using a Lorentz transformation, i.e., $\bar{p}_{(a)}=\alpha _{a}^{(b)}p^{\prime }_{(b)}$, where $\alpha _{a}^{(b)}$ is the transformation matrix. From Equations (86) and (87), we know that what one needs is just $\bar{p}_{(i)}/\bar{p}_{(t)}$ and

$$\begin{eqnarray} &&\frac{\bar{p}_{(i)}}{\bar{p}_{(t)}}=\frac{\alpha _i^{(t)}+\alpha _{i}^{(j)}p^{\prime }_{(j)}/p^{\prime }_{(t)}}{\alpha _t^{(t)}+\alpha _{t}^{(j)}p^{\prime }_{(j)}/p^{\prime }_{(t)}}. \end{eqnarray} \tag{ 89 }$$

$\alpha _{a}^{(b)}$ should be specified by the user according to his/her own needs.⁴

As an example, in Figure 2, we show a group of null geodesics emitted isotropically from a particle moving around a black hole in a marginally stable circular orbit (r_ms) with a = 0.9375. The physical velocities of the particle with respect to the LNRF are $\upsilon _r=\upsilon _\theta =0,\upsilon _\phi =e^{\psi -\nu }(\Omega -\omega)|_{r=r_{{\rm ms}}}$. The four-momentum $p^{\prime }_{(a)}$ are specified isotropically in the reference of the particle and then transformed into the LNRF by the Lorentz transformation expressed in Equation (90), i.e., $\bar{p}_{(a)}=\alpha _{a}^{(b)}p^{\prime }_{(b)}$. With $\bar{p}_{(a)}$, the constants of motion are computed readily. The light bending and beaming effects are illustrated clearly in this figure.

Zoom In Zoom Out Reset image size

In the next section, we discuss how to compute $\bar{p}_{(i)}/\bar{p}_{(t)}$ from impact parameters, which play a key role in imaging. For simplicity, we use the transformation expressed by Equations (90) and (91) if the observer has motion.

From the works of Cunningham & Bardeen (1973) and Cunningham (1975), we know that λ and q can be calculated from impact parameters, usually denoted α and β, which are the coordinates of the impact position of a photon on the photographic plate of the observer. The formulae provided by these authors read as follows (Cunningham & Bardeen 1973):

$$\begin{eqnarray} \lambda &=&-\alpha \sin \theta _{{\rm obs}}, \end{eqnarray} \tag{ 92 }$$

$$\begin{eqnarray} q&=&\beta ^2+(\alpha ^2-a^2)\cos ^2\theta _{{\rm obs}}. \end{eqnarray} \tag{ 93 }$$

The above equations are valid only when the distance between the observer and the emitter is infinite and the observer is stationary. Practically speaking, the distance is not infinite, otherwise the integrals of coordinate t become divergent. When the distance is finite, the above formulae should be modified. We extend those formulae to general situations in which both the finite distance and the motion state of the observer are considered.

Considering the finite distance is very easy. One just needs to substitute the coordinates r_obs, θ_obs of the observer into Equations (86) and (87) in the calculation. Considering the motion state is more complicated, however. Obviously we can divide the motion states of the observer into two types. In the first one, the observer is stationary and in the second one the observer has physical velocities $\upsilon _r,\upsilon _\theta,\upsilon _\phi$ with respect to the LNRF.

In the first instance, the observer is just an LNRF observer whose orthonormal tetrad is given by Equation (79), namely, e_(a)(obs) = e_(a)(LNRF) while in the second instance, the tetrad of the observer can be created by a Lorentz transformation, i.e., $\mathbf {e}_{(a)}(\mathrm{obs})=\widetilde{\alpha }_{(a)}^b \mathbf {e}_{(b)}(\mathrm{LNRF})$. Here, $\widetilde{\alpha }_{(a)}^b$ is given by Equation (91).

As shown in Figure 3, we can plot the image of a photon hitting the photographic plate. The plate is located in the plane determined by the basis vectors e_(θ)(obs) and e_(ϕ)(obs), and in which an orthonormal coordinate system α, β has been established. The basis vectors e_α and e_β of the system are aligned with e_(ϕ)(obs) and e_(θ)(obs), respectively. All photons go through the center of the lens before hitting the plate. From this figure, one obtains the relationships between the impact parameters and $p^{\prime }_{(a)}$ as follows:

$$\begin{eqnarray} &&\alpha =\left. r {scal} \frac{p^{\prime }_{(\phi)}}{p^{\prime }_{(r)}}\right|_{r=r_{{\rm obs}},\,\theta =\theta _{{\rm obs}}}, \end{eqnarray} \tag{ 94 }$$

$$\begin{eqnarray} &&\beta =\left. r {scal} \frac{p^{\prime }_{(\theta)}}{p^{\prime }_{(r)}}\right|_{r=r_{{\rm obs}},\,\theta =\theta _{{\rm obs}}}. \end{eqnarray} \tag{ 95 }$$

Clearly, one can read off $p^{\prime }_{(r)}> 0$, $\alpha p^{\prime }_{(\phi)}\ge 0$, and $\beta p^{\prime }_{(\theta)} \ge 0$.

Zoom In Zoom Out Reset image size

Two dimensionless factors, r and scal, have been multiplied to amplify the size of the image; otherwise it will be infinitely small since the distance D between the central compact object and the observer and the size of the target object L satisfy D ≫ L.

In the rest frame of the observer, spacetime is locally flat, and we still have

$$\begin{equation} -{p^{\prime }}_{(t)}^2+{p^{\prime }}_{(r)}^2+{p^{\prime }}_{(\theta)}^2+{p^{\prime }}_{(\phi)}^2=0. \end{equation} \tag{ 96 }$$

Using Equations (94), (95), and (96) and noting the signs of $p^{\prime }_{(a)}, \alpha, \beta$, we obtain

$$\begin{eqnarray} &&\frac{p^{\prime }_{(r)}}{p^{\prime }_{(t)}}=-\frac{1}{\sqrt{1+(\alpha /r\, {scal})^2+\left(\beta /r\, {scal}\right)^2}}, \end{eqnarray} \tag{ 97 }$$

$$\begin{eqnarray} &&\frac{p^{\prime }_{(\theta)}}{p^{\prime }_{(t)}}=-\frac{\beta /r {scal}}{\sqrt{1+(\alpha /r\, {scal})^2+\left(\beta /r\, {scal}\right)^2}}, \end{eqnarray} \tag{ 98 }$$

$$\begin{eqnarray} &&\frac{p^{\prime }_{(\phi)}}{p^{\prime }_{(t)}}=-\frac{\alpha /r {scal}}{\sqrt{1+(\alpha /r\, {scal})^2+\left(\beta /r\, {scal}\right)^2}}. \end{eqnarray} \tag{ 99 }$$

Substituting Equations (97)–(99) into Equation (89), we get the functions $\bar{p}_{(i)}(\alpha,\beta)/\bar{p}_{(t)}(\alpha,\beta)$. We can then calculate λ and q from the impact parameters by using Equations (86) and (87).

One can verify directly that when r_obs → ∞, scal = 1, and $\upsilon _r= \upsilon _\theta = \upsilon _\phi =0$, Equations (86) and (87) reduce to Equations (92) and (93) immediately.

If the observer has motion, the image on the plate has a displacement in comparison with the image when the observer is stationary. The displacement is proportional to the observer's velocity and can be described by α_c and β_c, which are the coordinates of the image point of the origin of the B–L coordinate system on the photographic plate. Obviously, α_c and β_c satisfy the following equations:

$$\begin{eqnarray} &&\bar{p}_{(\theta)}(\alpha _c,\beta _c)=0, \end{eqnarray} \tag{ 100 }$$

$$\begin{eqnarray} &&\bar{p}_{(\phi)}(\alpha _c,\beta _c)=0. \end{eqnarray} \tag{ 101 }$$

Actually, $\bar{p}_{(\phi)}(\alpha, \beta)=0$ represents the projection of the spin axis of the black hole onto the plate. When $\upsilon _r=\upsilon _\theta =\upsilon _\phi =0$, the above equations become α_c = 0 and β_c = 0. The region of the image on the plate therefore is [ − ΔL + α_c, ΔL + α_c] and [ − ΔL + β_c, ΔL + β_c], where ΔL is the half-length of the image.

The redshift g of a photon is defined by g = E_obs/E_em. From the above discussion, we know that $E_{{\rm obs}}=-p^{\prime }_{(t)}$ and $E_{{\rm em}}=-p_{\mu }u^{\mu }_{{\rm em}}$, where $u^{\mu }_{{\rm em}}$ is the four-velocity of the emitter. We define

$$\begin{equation} f_t = \alpha _t^{(t)}+\alpha _t^{(i)}\frac{p^{\prime }_{(i)}}{p^{\prime }_{(t)}}. \end{equation} \tag{ 102 }$$

From Equation (82), we have $-p^{\prime }_t=Ee^{-\nu }(1-\lambda \omega)/f_t$. Using Equation (81), g can be expressed as

$$\begin{equation} g=\frac{[e^{-\nu }(1-\lambda \omega)/f_t]_{{\rm obs}}}{[u^t (1+s_r\dot{r}\sqrt{R}/\Delta +s_{\theta }\dot{\theta }\sqrt{\Theta _\theta }-\lambda \Omega)]_{{\rm em}}}, \end{equation} \tag{ 103 }$$

where $\dot{r}=u^r/u^t,\dot{\theta }=u^\theta /u^t, \Omega =\dot{\phi }=u^\phi /u^t$ are the coordinate velocities.

With $\dot{r}, \dot{\theta }, \Omega$, the physical velocities of the emitter with respect to the LNRF $\upsilon _r, \upsilon _\theta, \upsilon _\phi$ can be written as (Bardeen et al. 1972)

$$\begin{equation} \upsilon _r=e^{\mu _1-\nu }\dot{r},\quad \upsilon _\theta =e^{\mu _2-\nu }\dot{\theta },\quad \upsilon _\phi =e^{\psi -\nu }(\Omega -\omega), \end{equation} \tag{ 104 }$$

with which the four-velocity of the emitter can be expressed as

$$\begin{equation} u^\mu _{{\rm em}}=\gamma (e^{-\nu },\upsilon _re^{-\mu _1},\upsilon _\theta e^{-\mu _2},\Omega e^{-\nu }). \end{equation} \tag{ 105 }$$

Then g can be rewritten as (Müller & Camenzind 2004)

$$\begin{equation} g=\frac{[e^{-\nu }(1-\lambda \omega)/f_t]_{{\rm obs}}}{[\gamma e^{-\nu }(1+s_r e^{\nu } \upsilon _{r}\sqrt{R}/\sqrt{\Sigma \Delta }+s_{\theta }e^{\nu }\upsilon _{\theta }\sqrt{\Theta _\theta }/\sqrt{\Sigma }-\lambda \Omega)]_{{\rm em}}}. \end{equation} \tag{ 106 }$$

For an emitter moving in a Keplerian orbit, the formula for g reduces to

$$\begin{equation} g=\frac{[e^{-\nu }(1-\lambda \omega)/f_t]_{{\rm obs}}}{[\gamma e^{-\nu }\left(1-\lambda \Omega \right)]_{{\rm em}}}. \end{equation} \tag{ 107 }$$

We have expressed the four coordinates r, μ, ϕ, t and the affine parameters σ as functions of p. We denote them as follows:

$$\begin{equation} r(p),\mu (p),\phi (p),t(p),\sigma (p). \end{equation} \tag{ 108 }$$

In practical applications, we are interested in determining the original position where the photon was emitted or the paths traveled by the photon. To make the calculations effectively, all photons are traced backward from the observer to the emitter along geodesics, but not all photons that start from the observer go through the emission region one is interested in, and the tracing process is terminated when these photons either reach infinity or fall into the event horizon of a black hole.

Now we discuss how to determine the intersection of a geodesic with the surface of an optically thick emission region, assuming that the optical depth is so high that a sharp emission surface exits. The surface is smooth and continuous and can be described by an algebraic equation,

$$\begin{equation} F(r,\theta,\phi)=F_0,\quad {\rm or} \quad F(x,y,z)=F_0, \end{equation} \tag{ 109 }$$

where x, y, and z are the pseudo-Cartesian coordinates defined by

$$\begin{equation} x=\sqrt{r^2+a^2}\sin \theta \cos \phi,\,\, y=\sqrt{r^2+a^2}\sin \theta \sin \phi,\,\, z=r\cos \theta. \end{equation} \tag{ 110 }$$

In some special cases, the surface one considers may not remain stationary, and the surface equation becomes a function of time t, i.e., F(r, θ, ϕ, t) = F₀. We introduce a function f(p) defined by

$$\begin{equation} f(p)=F[r(p),\theta (p),\phi (p)]-F_0. \end{equation} \tag{ 111 }$$

Then the roots of the equation f(p) = 0 correspond to the intersections of the geodesic with the target surface. Therefore, if the equation f(p) = 0 has no roots, the geodesic will never intersect with the surface. To solve this equation efficiently, we classify the geodesics into four classes denoted A, B, C, and D, according to their relationships with respect to a shell, as shown in Figure 4. The shell completely includes the emission region, and its inner and outer radii are r_in and r_out. We remind the reader that $r_{{\rm tp_1}}$ and $r_{{\rm tp_2}}$ are the turning points between which the radial motion of a photon is confined. Geodesics in the four classes satisfy the conditions—A: $r_{{\rm tp_1}}>r_{{\rm out}}$; B: $r_{{\rm in}}\le r_{{\rm tp_1}}\le r_{{\rm out}}$; C: $r_h< r_{{\rm tp_1}}< r_{{\rm in}}$; and D: $r_{{\rm tp_1}}\le r_h$, respectively, where r_h is the radius of the event horizon.

Zoom In Zoom Out Reset image size

The values of the parameter p corresponding to the intersections of the geodesic with the shell are denoted p₁, p₂, p₃, and p₄, which satisfy p₁ < p₂ < p₃ < p₄. Obviously, the roots of the equation f(p) = 0 may exist on intervals [p₁, p₂] and [p₃, p₄]. We use the Bisection or the Newton–Raphson method to search for the roots (Press et al. 2007).

Solving the radiative transfer equation in optically thin or thick media, one needs to evaluate integrations along the geodesic taking the affine parameter σ as an independent variable. Since we have taken p as an independent variable, we can replace σ by p to evaluate these integrals. From the definition of p, i.e., Equation (23), one has

$$\begin{equation} dp = \pm \frac{dr}{\sqrt{R(r)}}=\pm \frac{d\theta }{\sqrt{\Theta _{\theta }}}, \end{equation} \tag{ 112 }$$

and from Equations (4) and (5) one gets

$$\begin{equation} d\sigma = \pm \Sigma \frac{dr}{\sqrt{R(r)}}=\pm \Sigma \frac{d\theta }{\sqrt{\Theta _\theta }}. \end{equation} \tag{ 113 }$$

From the above equations one immediately obtains

$$\begin{equation} d\sigma = \Sigma dp, \end{equation} \tag{ 114 }$$

which converts the independent variable from σ to p in the radiative transfer applications (Yuan et al. 2009).

Finally, we discuss briefly the determination of a geodesic connecting the emitter and the observer (Viergutz 1993; Beckwith & Done 2005). We use α_em and β_em to represent the impact parameters of the geodesic connecting the observer and the emitter and p_em to indicate the position of the emitter on the geodesic, in which the coordinates of the emitter are r_em, μ_em, and ϕ_em. We have the following set of equations:

$$\begin{eqnarray} r(p_{{\rm em}},\alpha _{{\rm em}},\beta _{{\rm em}})&=&r_{{\rm em}}, \end{eqnarray} \tag{ 115 }$$

$$\begin{eqnarray} \mu (p_{{\rm em}},\alpha _{{\rm em}},\beta _{{\rm em}})&=&\mu _{{\rm em}}, \end{eqnarray} \tag{ 116 }$$

$$\begin{eqnarray} \phi (p_{{\rm em}},\alpha _{{\rm em}},\beta _{{\rm em}})&=&\phi _{{\rm em}}. \end{eqnarray} \tag{ 117 }$$

In principle, if we can solve this set of nonlinear equations simultaneously for p_em, α_em, and β_em, the geodesic is determined uniquely. Therefore, the observer–emitter problem also becomes a root-finding problem. In our code, we use the Newton–Raphson method (Press et al. 2007) to solve these equations.

In this section we give a brief introduction to the code; a more detailed introduction is given in the README⁵ file. The code is named YNOGK (Yun-Nan Observatory Geodesics Kerr) and it was written in Fortran 95 using the object-oriented method. The code is composed of a couple of modules. For each module, a special function has been implemented and one can use all supporting functions and subroutines in that module by the command "use module-name." By adding corresponding modules into one's own code, one can easily develop new modules to handle special and more sophisticated applications.

Two modules named ell-function and B–L coordinate are the most important ones in YNOGK; the former includes supporting functions and subroutines for calculating Carlson's elliptic integrals and the R-functions. Many routines in this module come from geokerr (Dexter & Agol 2009) and Numerical Recipes (Press et al. 2007). The latter includes routines for computing all coordinates and the affine parameter functions: t(p), r(p), θ(p), ϕ(p), and σ(p). To call these routines, constants of motion and the components of the four-momentum $\bar{p}_{(r)}$, $\bar{p}_{(\theta)}$, $\bar{p}_{(\phi)}$ measured in an LNRF must be provided. These values can be computed by two subroutines named lambdaq and initial-direction in YNOGK. The former routine computes the constants of motion from the impact parameters, while the latter computes λ and q from the initial $p^{\prime }_{(a)}$ given in a reference K', which has physical velocities $\upsilon _r, \upsilon _\theta, \upsilon _\phi$ with respect to the LNRF. Of course, one can compute them by one's own subroutines according to one's needs. As discussed in Section 5.1, we present a module named pem-finding to search for the minimum root p_em of the equation f(p) = 0; f(p) as an external function should be given by the user. In the module obs-emitter, we present routines to find the root of Equations (115)–(117) using the Newton–Raphson algorithm (Press et al. 2007). In the testing section of the code, this module has been used to determine the geodesic connecting the observer and the central point of a hot spot, which moves in the innermost stable circular orbit (ISCO). With the geodesic, the motion of the spot can be described easily. The results agree very well with those of previous works, in which a very different method was used to determine the motion of the spot, i.e., by tabulating the motion according to the time of the observer over one period (Schnittman & Bertschinger 2004; Dexter & Agol 2009).

All routines for computing Carlson's elliptic integrals have been extensively checked by the NIntegrate function in Mathematica. The original code for these routines comes from geokerr (Dexter & Agol 2009) and has been modified to adapt to our code. The original code for the computation of the R-functions comes from Numerical Recipes (Press et al. 2007). The same check has also been done for functions t(p), r(p), θ(p), ϕ(p), σ(p). When |α| and |β| ≲ 10⁻⁷, we let them be zero, since Carlson's integrals cannot maintain their accuracy. The treatment is the same for any other parameters if they take offending values. For some critical cases, special treatments also have been implemented.

Our code has many points in common with the geokerr code of Dexter & Agol (2009). Both use Carlson's method to compute the elliptic integrals and use the elliptic functions to express all coordinates as functions of an independent variable, but the elliptic function we used is mainly Weierstrass's elliptic function ℘(z; g₂, g₃), which has a cubic polynomial, leading to a simpler root distribution and a reduction in the number of integrals. In our code, the four B–L coordinates r, θ, ϕ, t and the affine parameter σ are expressed as analytical or numerical functions of a parameter p, which corresponds to I_u or I_μ in Dexter & Agol (2009). With this treatment, one can compute the geodesics directly without providing any information about the turning points in advance, enabling one to track emission from a more sophisticated surface instead of being limited to a standard thin accretion disk. In the code testing section, we show the images of a warped disk, which has a curved surface.

Our strategy, i.e., expressing coordinates as functions of p semi-analytically, can be extended to compute the timelike geodesics directly, almost without any modifications. As mentioned in Dexter & Agol (2009), the calculations of the timelike geodesics involve many more cases. The main challenge is to specify the number of radial turning points for bounded orbits in advance, but our strategy does not require the specification of the number of turning points both in radial and poloidal coordinates in advance. Therefore, our strategy can be used naturally and effectively in the calculations of timelike geodesics even in a Kerr–Newmann spacetime.

In our code, we give the orthonormal tetrad of the emitter or the observer analytically, provided that the physical velocities with respect to the LNRF are specified. They may be useful in a Monte-Carlo-type radiative transfer code, which requires one to make frequent transformations between the references frames of the emitter and the B–L coordinate system (Dolence et al. 2009). As illustrated in Figure 2, emission in the reference frame of the emitter is isotropic, but from the perspective of the B–L coordinate system the emission is anisotropic due to the Doppler beaming effect.

Our testing results for various toy problems agree well with those of Dexter & Agol (2009). In Figure 5, we illustrate the projection of a uniform orthonormal grid from the photographic plate of the observer onto the equatorial plane of a black hole. The solid and dotted lines represent the results from our code and geokerr, respectively. The results agree with each other very well.

Zoom In Zoom Out Reset image size

The basic strategies used in YNOGK to compute elliptic integrals and functions are very similar to those used in geokerr. For example, we make it possible for YNOGK to compute the minimum number of R-functions and share them between routines. This strategy greatly improves the speed of our code, but there are still some differences between the two codes. First, we assemble r(p) and μ(p) into routines for computing t(p), ϕ(p) and σ(p); thus, the repeated calculations for the same integrals among those functions can be avoided. We provide a routine named ynogk to compute the four B–L coordinates and affine parameter simultaneously. We also provide two independent routines named radius and mucos to compute r(p) and μ(p), respectively. Secondly, YNOGK can save the values of variables used in calculations for the same geodesic but for a different p. These values can be used repeatedly.

YNOGK is almost as fast as geokerr in tracing radiations from an optically and geometrically thin disk because only one point, i.e., the intersection of the ray with the disk surface, needs to be calculated. For the radiative transfer calculations, in which many points are needed along each geodesic, YNOGK is a little slower than geokerr. The speed tests of a code are dependent not only on the applications but also on the environment of the running code. From the testing results of YNOGK, we expect that its speed is almost the same as that of geokerr in many other applications.

For more detailed introductions, we refer the reader to the README file. In the next section, we show the testing results of our code for toy problems.

As the first test of our code, we provide the image of a black hole shadow. We trace all photons backward from the photographic plate to the black hole along geodesics. The intensities of the image are taken to be the affine parameter σ evaluated from the observer to the terminations on the geodesic—either when the photon intersects the event horizon of the black hole or reaches a turning point and returns to the starting radius. The evaluations of the affine parameter outside the shadow are multiplied by a factor of 1/2. In Figure 6, we show the image from an edge-on view. We take the spin a to be 0.998, and the distance of the observer to be 10⁶ r_g, where r_g is the gravitational radius.

Zoom In Zoom Out Reset image size

To evaluate the affine parameter σ from the function σ(p), we need p_h, which is the value of parameter p corresponding to the event horizon and also the root of the equation r(p) = r_h. We can get p_h by evaluating the integral of r in the definition of p and do not need to solve this equation directly. We provide a routine named r2p to complete this evaluation.

Next, we present images of the accretion disks around a Kerr black hole, including the standard thin, thick, and warped disks. The imaging of the disks is usually the first step in calculating the line profiles of the Fe Kα line and the spectrum (Li et al. 2005). Usually, the pseudo colors of the image represent the redshift g or the observed flux intensity I_ν of emission that comes from the disk.

In Figure 7, we show the image of a standard thin disk, the inner and outer radii of which are r_ms and 22 r_g, respectively. The black hole spin a is 0.998 and the inclination angle θ_obs is 86°. The distance of the observer is 40 r_g. The shape of the image is quite different from that observed from infinitely far away. We also illustrate the higher order images of the disk in this figure. Due to light bending and focusing, one can see the part behind the black hole and the bottom side of the disk. The color intensities represent the redshift g of emission that comes from the disk.

Zoom In Zoom Out Reset image size

In order to determine the intersection of the geodesics with the disk, we need to know the minimum root of the equation μ(p) = 0. We provide two routines named pemdisk and pemdisk-all to compute the root by evaluating the integral of μ in the definition of p. Using pemdisk, one can draw the direct image; while using pemdisk-all, one can draw the direct and higher-order images.

The surface of the thick disk has a constant inclination angle δ with respect to the equatorial plane (Wu & Wang 2007). To trace the thick disk, we need to solve the equation μ(p) = cos (π/2 − δ) to get p for the upper surface and μ(p) = cos (π/2 + δ) for the bottom surface; the roots of these two equations can also be computed with pemdisk and pemdisk-all. Since the surface particles of the disk no longer remain in the equatorial plane, they exhibit sub-Keplerian motion with an angular velocity given by (Ruszkowski & Fabian 2000)

$$\begin{equation} \Omega =\left(\frac{\theta }{\pi /2}\right)^{\frac{1}{n}}\Omega _K+ \left[1-\left(\frac{\theta }{\pi /2}\right)^{\frac{1}{n}}\right]\omega, \end{equation} \tag{ 118 }$$

where Ω_K = 1/(r^3/2 + a) is the Keplerian velocity and the parameter n is taken to be 3 here. Using Equation (107), we can calculate the redshift of the emission that comes from the disk. The images are shown in Figure 8 and they agree very well with Figure 10 of Wu & Wang (2007).

Zoom In Zoom Out Reset image size

Warped accretion disks are also very interesting objects in astrophysics (Bardeen & Petterson 1975; Wu & Wang 2007; Wang & Li 2012). Here, we discuss a very simple model of a warped disk, in which the disk is assumed to be optically thick and its surface can be described by (Wang & Li 2012)

$$\begin{equation} \cot \theta =-\tan \beta \cos (\phi -\gamma), \end{equation} \tag{ 119 }$$

where the parameters γ and β are defined by

$$\begin{eqnarray} \gamma (p)&=&\gamma _0+n_1\exp \left[n_2\frac{r_{{\rm in}}-r(p)}{r_{{\rm out}}-r_{{\rm in}}}\right], \end{eqnarray} \tag{ 120 }$$

$$\begin{eqnarray} \beta (p)&=&n_3\sin \left[\frac{\pi }{2}\frac{r(p)-r_{{\rm in}}}{r_{{\rm out}}-r_{{\rm in}}}\right], \end{eqnarray} \tag{ 121 }$$

where r_in and r_out are the inner and outer radii of the disk, and n₁, n₂, and n₃ are the warping parameters. With the above equations, we get f(p) as follows:

$$\begin{equation} f(p) = \tan \beta (p) \cos [\phi (p)-\gamma (p)]+\frac{\mu (p)}{\sqrt{1-\mu ^2(p)}}. \end{equation} \tag{ 122 }$$

With the minimum root of the equation f(p) = 0, we can image the warped disk. If the poloidal velocity $\dot{\theta }$ of the particle is nonzero, Equations (103) or (106) are used to calculate the redshift g (cf. Wang & Li 2012). The images of the warped disk are shown in Figure 9, in which the warping parameters n₁ and n₂ are nonzero, leading to disk warps along the azimuthal direction. For comparison, one can see Figure 3 of Wang & Li (2012), in which n₁ is taken to be zero for simplicity. The shape of the disk from Wang & Li (2012) is quite different from that illustrated here.

Zoom In Zoom Out Reset image size

In this section, we show the images of a rotationally supported torus. For simplicity, we give a brief introduction to the torus model here; for more detailed discussions we recommended papers by Fuerst & Wu (2004) and Younsi et al. (2012). The torus is assumed to be stationary and axisymmetric. Due to the balance of the centrifugal force, gravity, and the pressure force, the structure of the torus is stratified and the isobaric surfaces can be described by a set of differential equations (Younsi et al. 2012),

$$\begin{eqnarray} &&\frac{dr}{d\zeta } = \frac{\psi _2}{\sqrt{\psi _2^2+\Delta \psi _1^2}}, \end{eqnarray} \tag{ 123 }$$

$$\begin{eqnarray} &&\frac{d\theta }{d\zeta } = \frac{-\psi _1}{\sqrt{\psi _2^2+\Delta \psi _1^2}}, \end{eqnarray} \tag{ 124 }$$

where

$$\begin{eqnarray} \psi _1 & =& M\left(\frac{\Sigma -2r^2}{\Sigma ^2}\right)(\Omega ^{-1}-a\sin \theta)^2+r\sin ^2\theta, \end{eqnarray} \tag{ 125 }$$

$$\begin{eqnarray} \psi _2 &=& \sin 2\theta \left(\frac{Mr}{\Sigma ^2}[a\Omega ^{-1}-(r^2+a^2)]^2+\frac{\Delta }{2}\right), \end{eqnarray} \tag{ 126 }$$

$$\begin{eqnarray} \Omega &=& \frac{\sqrt{M}}{(r\sin \theta)^{3/2}+a\sqrt{M}}\left(\frac{r_k}{r\sin \theta }\right)^{n}, \end{eqnarray} \tag{ 127 }$$

and Ω is the angular velocity. r_k represents the radius at which the particle orbits with a Keplerian velocity. The index parameter n is crucial for regulating the angular velocity profile and adjusting the geometrical aspect ratio of the torus. ζ is an auxiliary parameter. In order to obtain the outer surface of the torus, one needs to specify the innermost radius of the torus. This radius is usually regarded as the intersection of the isobaric surface with the ISCO. Taking the innermost radius in the equatorial plane to be the initial condition, differential Equations (123) and (124) can now readily to be integrated. In Figure 10, we illustrate the images of the torus; the torus has the same parameters as in Figure 3 of Younsi et al. (2012) and the results agree with each other very well.

Zoom In Zoom Out Reset image size

When the trajectory of a photon is close to a compact object, the photon's path will be bent or focused and multiple images will be observed. This effect is the so-called gravitational lensing effect due to a strong gravitational field. Here, this phenomenon is illustrated by a simple example in which a ball moves around a near extremal black hole (a = 0.998) in a Keplerian orbit. If the radius of the orbit is R₀, then the angular velocity of the ball is $\Omega = 1/(R_0^{3/2}+a)$. The coordinates of the center of the ball become

$$\begin{eqnarray} &&x_0(p) = \sqrt{R_0^2+a^2}\cos [\Omega t(p)], \end{eqnarray} \tag{ 128 }$$

$$\begin{eqnarray} &&y_0(p) = \sqrt{R_0^2+a^2}\sin [\Omega t(p)], \end{eqnarray} \tag{ 129 }$$

$$\begin{eqnarray} &&z_0(p) = 0. \end{eqnarray} \tag{ 130 }$$

Then, the function of the surface of the ball can be expressed as follows:

$$\begin{equation} f(p)=\sqrt{[x(p)-x_0(p)]^2+[y(p)-y_0(p)]^2+z(p)^2}-r_1, \end{equation} \tag{ 131 }$$

where r₁ is the radius of the ball. The images observed from an edge-on view are illustrated in Figure 11. For different positions of the ball in its orbit, the image changes significantly and even becomes a ring as the ball moves to the back of the event horizon.

Zoom In Zoom Out Reset image size

The calculation of line profiles is very easy, provided that the structure of the disk is specified (Bromley et al. 1997). For simplicity, we assume that the particles of the accretion flow have Keplerian motion and that the disk is geometrically thin and optically thick. The inner and outer radii of the disk are located at r_ms and 15 r_g, respectively. The emission is monochromatic and its profile can be described by Dirac's δ function in the local rest frame of the flow

$$\begin{equation} I_{{\rm em}}(\nu)=\frac{\epsilon _0}{4\pi r^n}\delta (\nu -\nu _{{\rm em}}), \end{equation} \tag{ 132 }$$

where n is the index of emissivity and is assumed to be 3. Since I_ν/ν³ is invariant along a geodesic (Misner et al. 1973), we obtain the observed intensity I_obs = g³I_em, where g = ν_obs/ν_em is the redshift. Then, the observed flux density F_ν at a frequency ν can be computed by integrating I_obs over the whole plate as in the following expression:

$$\begin{equation} F_{\nu }=\int \frac{\epsilon _0}{4\pi r^n}\delta (\nu -\nu _{{\rm em}})g^3 d\alpha d\beta. \end{equation} \tag{ 133 }$$

The observed intensities have been normalized in the computation. The results are shown in Figure 12, which agrees very well with Figure 3 of Čadež et al. (1998). From this figure, one can see that larger black hole spins lead to broadening at low frequencies, for the ISCO is closer to the event horizon, and gravitational redshift effects are more significant. Higher inclination angles lead to broadening at high frequencies for the Doppler beaming effect.

Zoom In Zoom Out Reset image size

In this section, we compute the spectrum of a Keplerian disk around a Kerr black hole to illustrate the effects of black hole spin and the observer's inclination angle on the observed profile of the spectrum (Li et al. 2005). Similarly, the disk is assumed to be geometrically thin and optically thick, and the radiation spectrum of the disk in its local rest frame is assumed to be an isotropic blackbody spectrum. We denote the effective temperature of the disk by T_eff. Then, the radiation intensity at frequency ν can be written as

$$\begin{equation} I_{{\rm em}}(\nu)=\frac{h\nu ^3}{\exp (h\nu /k_{\mathrm{B}}T_{\mathrm{eff}})-1}, \end{equation} \tag{ 134 }$$

where h and k_B are the Planck and Boltzmann constants, respectively. For blackbody radiation, the effective temperature is simply

$$\begin{eqnarray} &&T_{\mathrm{eff}}=\left[\frac{F(r)}{\sigma _{{\rm SB}}}\right]^{1/4}, \end{eqnarray} \tag{ 135 }$$

where σ_SB is the Stefan–Boltzmann constant. Here, we do not consider the effect of the returning radiation of the disk on the spectrum; therefore F(r) is just the energy flux emitted from the disk's surface measured by a locally corotating observer. For the Keplerian accretion disk around a Kerr black hole, Page & Thorne (1974) presented the analytical expression for F(r), i.e.,

$$\begin{eqnarray} &&F(r)=\frac{\dot{M}}{4\pi r}f, \end{eqnarray} \tag{ 136 }$$

where $\dot{M}$ is the mass accretion rate, f is a function of r, a, r_ms, the seminal expression of which is given by Equations (15d) and (15n) of Page & Thorne (1974). With F(r), the effective temperature of the disk can be computed readily. Using the invariance I_ν/ν³, one can get the observed intensity I_obs = g³I_em. The total observed flux density at frequency ν therefore is the integration of I_obs over the whole plate:

$$\begin{equation} F(\nu _{{\rm obs}})=\int \frac{h\nu _{{\rm obs}}^3 }{\exp (h\nu _{{\rm obs}}/gk_\mathrm{B}T_{\mathrm{eff}})-1}d\alpha d\beta. \end{equation} \tag{ 137 }$$

Then the photon number flux density is N_obs = F(ν_obs)/E_obs.

The results are plotted in Figure 13. Compared with Figure 5 of Li et al. (2005), we find that the basic features of the two figures are in agreement. For example, as shown in the top panel, we see that as the spin of the black hole increases the spectrum becomes harder. Physically, this effect is due to the fact that as the spin a increases, the system of the accretion disk has a higher radiation efficiency and a higher temperature. In the bottom panel of Figure 13, we can see that at the low-energy end, the flux density decreases as θ_obs increases. As explained by Li et al. (2005), this is caused by the projection effect. While at the high-energy end, the flux density increases as θ_obs increases. As pointed out by Li et al. (2005), this effect results from both Doppler beaming and gravitational focusing.

Zoom In Zoom Out Reset image size

In the top panel, there is a noticeable effect: even though we do not consider the returning radiation, the flux density increases as the spin increases at the low-energy end. Li et al. (2005) suggested that this effect is caused by the returning radiation. We propose that this explanation may not be correct and the effect is just caused by the simple fact that a higher spin leads to higher radiation efficiency and temperature.

In order to test the validity of the function t(p) and illustrate the time delay effect in the Kerr spacetime, we image a hot spot that moves around a black hole in a prograde marginally stable circular orbit for various spins. We also compute the observed time-dependent light curve and spectra for a spot in retrograde ISCO around a Schwarzschild black hole for various inclination angles. The radius of the spot is R_spot = 0.5 r_g. The emissivity of the spot is taken to be Gaussian in its rest frame (Schnittman & Bertschinger 2004), i.e.,

$$\begin{equation} j(\mathbf {x})\propto \exp \left[-\frac{|\mathbf {x}-\mathbf {x}_{\mathrm{spot}}(t)|^2}{2 R_{\mathrm{spot}}}\right]. \end{equation} \tag{ 138 }$$

For the motion of the spot, one must consider the time delay effect and the azimuthal position of the spot when imaging the spot and calculating its spectra. In order to compute the time delay Δt for each geodesic starting from the photographic plate, a reference time t_obs, which is taken to be the time required for a photon to travel from the central point of the spot to the observer, needs to be specified. Meanwhile, the position of the spot can be determined by its central coordinates (r_ms, μ = 0, ϕ_spot). Then, with the method discussed in Section 5 (i.e., the method for determining a geodesic connecting the observer and emitter with the given coordinates), we can determine the geodesic connecting the central point of the spot and the observer. With this geodesic, the reference time t_obs can be calculated readily. Using t_obs, we can easily calculate the time delay Δt for each geodesic, and Δt = t_geo − t_obs, where t_geo is the time required for a photon to travel from the observer to the disk following the geodesic. With Δt and the position of the spot, we can compute the distance between the intersection of the geodesic with the disk and the center of the spot, i.e., |x − x_spot|. Thus, the emissivity can be computed readily.

In Figure 14, we illustrate the images of the spot with different black hole spins. As the spin increases, the marginally stable circular orbit is closer to the event horizon of the black hole, and the time delay effect becomes significant. The image of the spot is strongly warped, especially when the spot moves to the back of the event horizon.

Zoom In Zoom Out Reset image size

When an image is obtained, the redshift and Gaussian emissivity of all points on the spot can be computed. Repeating this procedure over one period of the motion gives a time-dependent spectrum. Integrating the spectrum over frequency, or equivalently over the impact parameters, gives the light curve. The spectrum and light curve are shown in Figure 15, and agree well with the results shown in Figures 6 and 7 of Dexter & Agol (2009).

Zoom In Zoom Out Reset image size

6.8.1. The Radiative Transfer Formulation

In this section, we briefly discuss radiative transfer in Kerr spacetime. One can find more detailed discussions from Fuerst & Wu (2004) and Younsi et al. (2012). It is well known that $\mathcal {I} = I_\nu /\nu ^3$, χ = να_ν, and η = j_ν/ν² are Lorentz invariants, where I_ν is the specific intensity of the radiation and α_ν and j_ν are the absorption and emission coefficients at frequency ν, respectively. The radiative transfer equation reads (Younsi et al. 2012)

$$\begin{equation} \frac{d\mathcal {I}}{d\tau _\nu } = -\mathcal {I}+\frac{\eta }{\chi }, \end{equation} \tag{ 139 }$$

where τ_ν is the optical depth at frequency ν and is defined by dτ_ν = α_νds and ds = −p_μu^μdσ, in which ds is the differential distance element of a photon traveling in the rest frame of the medium, σ is the affine parameter, p_μ is the four-momentum of the photon, and u^μ is the four-velocity of the medium. Then, the radiative transfer equation can be rewritten as (Younsi et al. 2012)

$$\begin{equation} \frac{d\mathcal {I}}{d\sigma } = -p_\mu u^\mu |_\sigma \left(-\alpha _\nu \mathcal {I}+\frac{j_\nu }{\nu ^3}\right). \end{equation} \tag{ 140 }$$

The solution of the above equation is (Younsi et al. 2012)

$$\begin{equation} \mathcal {I}(\sigma) = \mathcal {I}(\sigma _0)e^{-\tau _\nu (\sigma)}-\int ^\sigma _{\sigma _0}\frac{j_\nu (\sigma ^{\prime \prime })}{\nu ^3} \exp \left(-\int ^\sigma _{\sigma ^{\prime \prime }}{\alpha _\nu (\sigma ^{\prime }) |p_\mu u^\mu |_{\sigma ^{\prime }}|d\sigma ^{\prime }}\right) p_\mu u^\mu |_{\sigma ^{\prime \prime }}d\sigma,^{\prime \prime } \end{equation} \tag{ 141 }$$

where the optical depth is

$$\begin{equation} \tau _\nu (\sigma) = -\int ^\sigma _{\sigma _0}\alpha _\nu (\sigma ^{\prime }) p_\mu u^\mu |_{\sigma ^{\prime }} d\sigma ^{\prime }. \end{equation} \tag{ 142 }$$

As discussed in Section 5.1, we can convert the independent variable from the affine parameter σ into the parameter p. Using σ = σ(p) and dσ = Σdp, we can rewrite the solution as the integration of the parameter p (Yuan et al. 2009)

$$\begin{equation} \mathcal {I}(p) = \mathcal {I}(p_0)e^{-\tau _\nu (p)}-\int ^p_{p_0}\frac{j_\nu (p^{\prime \prime })}{\nu ^3} \exp \left(-\int ^p_{p^{\prime \prime }}{\alpha _\nu (p^{\prime }) |p_\mu u^\mu |_{p^{\prime }}|\Sigma ^{\prime } dp^{\prime }}\right) p_\mu u^\mu |_{p^{\prime \prime }}\Sigma ^{\prime \prime } dp,^{\prime \prime } \end{equation} \tag{ 143 }$$

where

$$\begin{equation} \tau _\nu (p) = -\int ^p_{p_0}\alpha _\nu (p^{\prime }) p_\mu u^\mu |_{p^{\prime }}\Sigma ^{\prime } dp^{\prime }. \end{equation} \tag{ 144 }$$

With the above formulae one can deal with radiative transfer problems without considering the scattering contributions to the absorption and emission coefficients as was done by Yuan et al. (2009) and Younsi et al. (2012).

6.8.2. Radiative Transfer in Pressure-Supported Torus

In Section 6.3, we have discussed a rotationally supported torus and showed its images. When the torus is optically thick, only the emission that comes from the boundary surface is considered. When the torus is optically thin, all parts of the torus contribute to the observed emission. We need to consider the radiative transfer procedure along a ray inside the torus. To get the absorption and emission coefficients, we need to know the structure model of the tours, which determines the distributions of temperature, mass density, pressure, etc.

Firstly, we construct a model of the torus, in which the torus is a perfect fluid and its energy-momentum tensor is given by (Younsi et al. 2012)

$$\begin{equation} T^{\alpha \beta }=(\rho +P+\epsilon)u^\alpha u^\beta +P g^{\alpha \beta }, \end{equation} \tag{ 145 }$$

where ρ is the mass density, P is the pressure, is the internal energy, u^α is the four-velocity of the fluid, and g^αβ are the contravariant components of the Kerr metric. From the conservation law, namely, $T^{\alpha \beta }_{\quad;\beta } = 0$, we get the equation of motion of the fluid as follows (Abramowicz et al. 1978):

$$\begin{equation} \frac{\partial _\alpha P}{\rho +P+\epsilon }=-u_{\alpha {;\beta }}u^\beta, \end{equation} \tag{ 146 }$$

where the semicolon ";" represents the covariant derivative and u_{α; β}u^β = a_α is the four acceleration of the fluid. For a stationary and axisymmetric torus, we have a_t = 0, a_ϕ = 0, and a_r, a_θ given by (Younsi et al. 2012)

$$\begin{eqnarray} a_r & =& -\dot{t}^2\left[M\left(\frac{\Sigma -2r^2}{\Sigma ^2}\right)\left(1-a\sin \theta \Omega \right)^2+r\sin ^2\theta \Omega ^2\right], \end{eqnarray} \tag{ 147 }$$

$$\begin{eqnarray} a_\theta &=& -\dot{t}^2\sin 2\theta \left(\frac{Mr}{\Sigma ^2}[a-(r^2+a^2)\Omega ]^2+\frac{\Delta \Omega ^2}{2}\right), \end{eqnarray} \tag{ 148 }$$

where $\dot{t}=u^t$ is the time component of the four-velocity and Ω is the angular velocity and takes the form of Equation (127). They satisfy the following equation

$$\begin{eqnarray} &&u^t=\frac{1}{\sqrt{-(g_{tt}+2g_{t\phi }\Omega +g_{\phi \phi }\Omega ^2)}}. \end{eqnarray} \tag{ 149 }$$

Since the torus is assumed to be radiation dominated, the pressure P can be regarded as the sum of the gas pressure P_gas and the radiation pressure P_rad, and

$$\begin{eqnarray} P_{\mathrm{gas}} &=& \frac{\rho k_\mathrm{B} T}{\mu m_\mathrm{H}}=\beta P, \end{eqnarray} \tag{ 150 }$$

$$\begin{eqnarray} P_{\mathrm{rad}} &=& \frac{\sigma T^4}{3} = (1-\beta)P, \end{eqnarray} \tag{ 151 }$$

where k_B is the Boltzmann constant, μ is the mean molecular weight, m_H is the mass of hydrogen, β is the ratio of gas pressure to the total pressure, and σ = π²k⁴/15ℏ³c³ is the blackbody emission constant. From the above equations, one finally obtains

$$\begin{eqnarray} &&P = \hbar c\left[\frac{45(1-\beta)}{\pi ^2(\mu m_\mathrm{H}\beta)^4}\right]^{1/3}\rho ^{4/3}, \end{eqnarray} \tag{ 152 }$$

$$\begin{eqnarray} &&k T =\hbar c\left[\frac{45(1-\beta)}{\pi ^2\mu m_\mathrm{H}\beta }\right]^{1/3}\rho ^{1/3}. \end{eqnarray} \tag{ 153 }$$

Thus $P = \kappa \rho ^\Gamma$, which implies that the state equation of the fluid is polytropic, therefore its internal energy is proportional to the pressure = P/(Γ − 1), and the equation of motion of the fluid (146) becomes

$$\begin{eqnarray} &&\left({\rho + \frac{\Gamma }{\Gamma -1}P}\right)a_\alpha =- \partial _\alpha P. \end{eqnarray} \tag{ 154 }$$

Substituting ∂_αP = κΓρ^{Γ − 1}∂_αρ and $P = \kappa \rho ^\Gamma$ into the above equation, one obtains

$$\begin{eqnarray} &&\partial _\alpha \rho = -a_\alpha \left(\frac{\rho ^{2-\Gamma }}{\kappa \Gamma }+\frac{\rho }{\Gamma -1}\right). \end{eqnarray} \tag{ 155 }$$

Introducing a new variable ξ defined by ξ = ln (Γ − 1 + κΓρ^{Γ − 1}), the above equation is simplified as

$$\begin{eqnarray} &&\partial _\alpha \xi = -a_\alpha, \end{eqnarray} \tag{ 156 }$$

which implies that the vector n = (a_r, a_θ) in the r–θ plane can be regarded as the normal vector of the contours of the density ρ. Thus, if we use t = (dr, dθ) to denote the tangent vector of the contours, we have n · t = 0, or equivalently

$$\begin{eqnarray} &&a_rdr+a_\theta d\theta = 0. \end{eqnarray} \tag{ 157 }$$

If we use ds to denote the differential proper length of the tangent vector, we have

$$\begin{eqnarray} &&ds^2= \mathbf {t}\cdot \mathbf {t}=g_{rr}dr^2+g_{\theta \theta }d\theta ^2, \end{eqnarray} \tag{ 158 }$$

where g_rr and g_θθ are the components of the Kerr metric and g_rr = Σ/Δ, g_θθ = Σ. Solving Equations (157) and (158) simultaneously, we get a set of differential equations to describe the contours of density ρ

$$\begin{eqnarray} &&\frac{dr}{ds}=\sqrt{\frac{\Delta }{\Sigma }}\frac{|a_\theta |}{\sqrt{a_\theta ^2+\Delta a_r^2}}, \end{eqnarray} \tag{ 159 }$$

$$\begin{eqnarray} &&\frac{d\theta }{ds}=-\sqrt{\frac{\Delta }{\Sigma }}\frac{|a_r|}{\sqrt{a_\theta ^2+\Delta a_r^2}}. \end{eqnarray} \tag{ 160 }$$

If we introduce an auxiliary variable ζ defined by $d\zeta =\sqrt{\Delta /\Sigma }ds$ and substitute Equations (147) and (148) into the above equations, we get

$$\begin{eqnarray} &&\frac{dr}{d\zeta } = \frac{\psi _2}{\sqrt{\psi _2^2+\Delta \psi _1^2}}, \end{eqnarray} \tag{ 161 }$$

$$\begin{eqnarray} &&\frac{d\theta }{d\zeta } = \frac{-\psi _1}{\sqrt{\psi _2^2+\Delta \psi _1^2}}, \end{eqnarray} \tag{ 162 }$$

which have exactly the same forms as Equations (123) and (124), where ψ₁ and ψ₂ are given by Equations (125) and (126). With these equations, the distributions of the mass density ρ of the torus now can be computed by evaluating the integral of ξ from the torus center r = r_k, ρ = ρ_c to the location (r, θ) along a path C which is orthogonal to the density contours everywhere. The integral of ξ is

$$\begin{eqnarray} &&\xi = -\int _{C}a_r dr+a_\theta d\theta. \end{eqnarray} \tag{ 163 }$$

From Equations (152) and (153), one can get the total pressure and temperature distributions immediately with the given density ρ.

Knowing the structure model of the torus, the absorption and emission coefficients are now readily specified. We use these coefficients to discuss the radiative transfer process inside the torus. Using the above torus model, we give two examples of radiative transfer applications.

Firstly, we consider a rather simple case in which the torus is optically thin. The emissivity is taken to be proportional to the mass density ρ, i.e., j_em∝ρ, and is independent of the frequency ν. The absorption coefficient α_ν is simply assumed to be zero. The torus parameters are n = 0.21 and r_k = 12 r_g. The black hole spin a is 0.998. The ratio of gas pressure to total pressure β is 2.87 × 10⁻⁸. In Figure 16, we show images of the torus, which is optically thin and radiation pressure dominated. We see that the emission mainly comes from the central region of the torus where the density is higher. As the inclination angle of the observer increases, the frequency shift of the emission from the Doppler boosting becomes larger. In this figure, the false color represents the observed intensities of the emission, showing that the approaching side of the torus is brighter than the receding side, especially at higher inclination angles.

Zoom In Zoom Out Reset image size

Secondly, we mimic a more realistic case, namely, thermal free–free emission and absorption, in which the torus is semi-opaque. The emission and absorption coefficients of the torus for a photon with energy E₀ are given by (Younsi et al. 2012)

$$\begin{eqnarray} j(E_0)& =& \mathcal {K}\left(\frac{n_e}{{\rm cm^{-3}}}\right)^2\left(\frac{E_0}{{\rm keV}}\right)^{-1}\left(\frac{\Theta }{{\rm keV}}\right)^{-1/2}e^{-E_0/\Theta }, \end{eqnarray} \tag{ 164 }$$

$$\begin{eqnarray} \alpha (E_0)& = &B_1\left(\frac{n_e}{{\rm cm^{-3}}}\right)^2 \frac{\sigma _\mathrm{T}}{E_0^2}, \end{eqnarray} \tag{ 165 }$$

where Θ = k_BT, $\mathcal {K}$ and B₁ are the normalization constants, n_e is the electron number density and n_e = ρ/μm_H, σ_T is the Thompson cross-section. The observed intensity images of the optically thick and semi-opaque torus are plotted in Figure 17. These images are quite different from those of an optically thin torus. The emissivity now depends on the temperature, which decreases toward the outer surface of the torus, leading to the limb darkening phenomenon. When the rays are nearly tangential to the layers of the torus, they travel a longer distance and go through the outer (colder) layers. On the other hand, when the rays are perpendicular to the layers of the torus, they travel a shorter distance and go through the inner (hotter) layers. Consequently, the observed intensity at lower inclination angles becomes much brighter than that at higher inclination angles (Younsi et al. 2012).

Zoom In Zoom Out Reset image size