Free Access
Issue
A&A
Volume 595, November 2016
Article Number A38
Number of page(s) 8
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/201629075
Published online 26 October 2016

© ESO, 2016

1. Introduction

Early studies (Luminet 1979; Pechenick et al. 1983) began a great interest in photons emitted by matter in a strong gravitational field, especially in relation to high-energy astrophysics. The relevant computations are carried out with ray-tracing techniques that are based on the photon geodesics in general relativistic spacetimes. Effects to be considered are (i) light bending; (ii) travel time delay; and (iii) gravitational lensing (Misner et al. 1973). The basic equations for the Schwarzschild metric are expressed through elliptic integrals that can be solved numerically. A powerful analytical approximation was introduced by Beloborodov (2002), who derived an approximate linear equation to describe the gravitational light bending of photons emitted at radius \hbox{$r \geqslant r_{\rm s}$} (rs = 2GM/c2). In the same vein, Poutanen & Beloborodov (2006) derived an approximate polynomial equation for photon travel time delays. These two analytical approximations were obtained by introducing an ad hoc parametrization of the photon emission angle (see Beloborodov 2002; Poutanen & Beloborodov 2006, for more details). Nevertheless, the equation for gravitational lensing, also known as solid angle equation, was still solved numerically by these authors.

In this paper we present a mathematical method through which the approximate polynomial equations for light bending and travel time delay in a Schwarzschild spacetime are derived without any ad hoc assumption. We then apply the same method to derive for the first time an approximate polynomial equation for gravitational lensing. High-accuracy approximate equations for photon geodesics translate into high-speed ray-tracing codes for different astrophysical applications in the strong gravitational field of Schwarzschild black holes (BHs). As examples we apply our approximate equations to calculate the light curve from a hot spot on the surface of a rotating neutron star (NS) and a clump in a circular orbit around BH. Moreover, we calculate the fluorescent iron Kα line profile from an accretion disk around a BH (e.g., Fabian et al. 1989).

2. Photons in the Schwarzschild spacetime

In this section we introduce the elliptical integrals of photon trajectories, travel time delay, and gravitational lensing in the Schwarzschild metric.

2.1. Schwarzschild metric

For static and spherically symmetric BHs of mass, M, the Schwarzschild metric in spherical coordinates (t,r,ϕ,ψ) is ds2=(1rsr)dt2(1rsr)-1dr2r2(dϕ2+sin2ϕdψ2),\begin{equation} {\rm d}s^2=\left (1-\frac{r_{\rm s}}{r} \right ) {\rm d}t^2 - \left (1-\frac{r_{\rm s}}{r} \right )^{-1}{\rm d}r^2-r^2\left({\rm d}\varphi^2+ \sin^2\varphi\, {\rm d}\psi^2\right), \label{eq:sm} \end{equation}(1)where G = c = 1, and rs = 2M is the Schwarzschild radius. In this standard system, the coordinate variables are time t, radius r, polar angle ϕ, and azimuthal angle ψ.

2.2. Gravitational light bending

Because of spherical symmetry, it is customary to use the equatorial plane at ϕ = π/ 2 to calculate geodesics in the Schwarzschild metric that are representative of all photon trajectories. A photon geodesic starting at radius R is described by the following elliptical integral (Chandrasekhar 1992; Misner et al. 1973): ψ=Rdrr2[1b21r2(1rsr)]12,\begin{equation} \label{eq:libe} \psi=\int_R^\infty \frac{{\rm d}r}{r^2}\left[\frac{1}{b^2}-\frac{1}{r^2}\left(1-\frac{r_{\rm s}}{r} \right) \right]^{-\frac{1}{2}}, \end{equation}(2)parametrized by the ratio of the angular momentum, L, and energy, E, of the photon, b = L/E. The impact parameter b represents the distance between the observer and the photon trajectory at infinity and is related to the photon emission angle, α by b=Rsinα1rs/R·\begin{equation} \label{eq:impact} b = \frac{R\sin\alpha}{\sqrt{1-r_{\rm s}/R}}\cdot \end{equation}(3)Equation (2) is strictly valid up to α = π/ 2, since the sine function is symmetric with respect to α = π/ 2. The photon deflection angle, ψ, can be directly determined in terms of the emission angle α through Eq. (3).

thumbnail Fig. 1

Two photon trajectories emitted at different radii, r1 and r2, and emission angles, α1 and α2, with their corresponding impact parameters, b1 and b2. Trajectory 1 is for a direct photon, while trajectory 2 has a turning point, i.e., passes through periastron p, the minimum distance between the trajectory and the BH. The observer is at infinity, and photons geodesics lie in a single invariant plane.

We must distinguish between direct photons, which have trajectories with an emission angle between 0 ≤ απ/ 2, and photons with a turning point, whose trajectories have an emission angle ranging between π/ 2 ≤ ααmax (see Fig. 1). Photon trajectories with a turning point can reach infinity only if their b is greater than the critical impact parameter bc=33M\hbox{$b_{\rm c}=3\sqrt{3}M$} (see, e.g., Luminet 1979). Since we are interested only in photons that are not captured by the BH, the maximum possible emission angle is obtained by substituting bc into Eq. (3) αmax=πarcsin[323(1rsR)rsR]·\begin{equation} \label{alphamax} \alpha_{\rm max}=\pi-\arcsin\left[\frac{3}{2}\sqrt{3\left(1-\frac{r_{\rm s}}{R}\right)}\frac{r_{\rm s}}{R}\right]\cdot \end{equation}(4)Photons emitted between π/ 2 ≤ ααmax follow trajectories with a turning point; therefore a periastron distance, p, is defined at an angle αp = π/ 2, which determines the minimum distance between the compact object and the photon trajectory. The emission point of a photon at ψE that passes through the turning point is symmetric with respect to the periastron angle, ψp, to the point ψS, (with an emission angle απ/ 2) along the same trajectory, as they have the same impact parameter at infinity. Based on this symmetry, we determine ψS = 2ψpψE, where αS = παE.

2.3. Travel time delay

A photon following its geodesic from an emission point, E, to an observer at infinity has an infinite travel time, Δτ, value. To have a finite quantity, we calculate the relative travel time delay between a photon emitted at a distance, R, following its geodesic and the photon emitted radially with b = 0, that is, Δt(b) = Δτ(b)−Δτ(b = 0) (Pechenick et al. 1983). In the Schwarzschild metric we have Δt=Rdr1rsr[1b2r2(1rsr)]121.\begin{equation} \label{eq:tide} \Delta t =\int_R^\infty \frac{{\rm d}r}{1-\frac{r_{\rm s}}{r}} \left \{ \left [1-\frac{b^2}{r^2}\left(1-\frac{r_{\rm s}}{r} \right) \right]^{-\frac{1}{2}} -1 \right \}. \end{equation}(5)To calculate the time delay for photons with a turning point, we need to calculate the periastron distance, p. For a given b we therefore consider the largest real solution of the following equation p3b2p + b2rs = 0. The polynomial in p has three real solutions (because bbc): one is negative, one is lower than 3M, and we consider only the solution satisfying prc, where rc = 3M is the critical radius associated to bc (see, e.g., Luminet 1979). The time delay is composed of the time delay ΔtS from point αS, as determined by the Eq. (5), plus the time delay between [αS,αp], Δtp−S, and [αp,αE], ΔtE−p. Since the integrand is symmetric with respect to αp, the latter two time delays are equal (ΔtE−p = Δtp−S), the equation can be written (see Fig. 2) Δt=ΔtS+2ΔtpS=ΔtS+2Rpdr1rsr[1b2r2(1rsr)]12·\begin{eqnarray} \label{eq:tidetp} \Delta t &=&\Delta t_{\rm S}+2\Delta t_{\rm p-S}\nonumber \\ &=&\Delta t_{\rm S}+2\int_R^p \frac{{\rm d}r}{1-\frac{r_{\rm s}}{r}} \left \{ \left [1-\frac{b^2}{r^2}\left(1-\frac{r_{\rm s}}{r} \right) \right]^{-\frac{1}{2}}\right \}\cdot \end{eqnarray}(6)

thumbnail Fig. 2

Calculation of travel time delay for trajectories with turning points. The photon is emitted at E with radius R and deflection angle ψE. The photon trajectory passes through point, S, which is symmetric to E with respect to periastron p, having the same impact parameter, b, and a deflection angle ψS = 2ψpψE.

2.4. Solid angle

We consider the emission reference frame of coordinates (x,y,z) and the observer reference frame of coordinates (x′,y′,z′), where the two systems are rotated with an angle, i, around y = y. The solid angle, , in the observer reference frame reads as dΩ = sinψ dψ dϕ. This equation can be expressed in terms of the impact parameter, b, by its first-order approximation for infinitesimally small ψ as bD·ψ, where D is the distance from the emission point to the observer, =bdbdϕD2·\begin{equation} \label{CASF} {\rm d}\Omega=\frac{b\ {\rm d}b\ {\rm d}\varphi'}{D^2}\cdot \end{equation}(7)In the emission reference frame, Eq. (7) becomes dΩ=bD2ϕ∂ϕ∂b∂rdrdϕ,\begin{equation} \label{NCASF} d\Omega=\frac{b}{D^2}\frac{\partial \varphi'}{\partial \varphi}\frac{\partial b}{\partial r} {\rm d}r {\rm d}\varphi, \end{equation}(8)where we considered the following dependencies ϕ = ϕ(ϕ′) and b = b(r,ψ). The Jacobian of the transformation is always ϕ∂ϕ∂b∂r\hbox{$\frac{\partial \varphi'}{\partial \varphi}\frac{\partial b}{\partial r}$} independent of the value of ∂b∂ψ\hbox{$\frac{\partial b}{\partial \psi}$}, since the photon moves in an invariant plane. Therefore, Eq. (8) is valid for any emission point. To calculate the Jacobian, we use the following coordinates transformation cosψ = sinicosϕ that relates the angles in the observer and emission reference frames. ∂b∂r=∂b∂ψ∂ψ∂r\hbox{$\frac{\partial b}{\partial r}=-\frac{\partial b}{\partial \psi}\frac{\partial \psi}{\partial r}$} is calculated using the light bending Eq. (2). The solid angle equation in the Schwarzschild metric is thus (see, e.g., Bao et al. 1994)1=cosiD2R2sin2ψb2cosαRdrr2[1b2r2(1rsr)]32drdϕ.\begin{equation} \label{EFSA} {\rm d}\Omega=\frac{\frac{\cos i}{D^2\ R^2\ \sin^2\psi}\frac{b^2}{\cos\alpha}} {\int_R^\infty \frac{{\rm d}r}{r^2}\left[1-\frac{b^2}{r^2}\left(1-\frac{r_{\rm s}}{r} \right) \right]^{-\frac{3}{2}}}\ {\rm d}r\ {\rm d}\varphi. \end{equation}(9)This equation contains an integral with the same functional form as those of light bending Eq. (2) and time delay Eq. (5), except for the −3/2 exponent and factors depending on the impact parameter b (or emission angle α).

3. Analytical approximations

In this section we present the general mathematical method used to approximate the elliptical equations in polynomials of light bending Eq. (2), time delay Eq. (5), and solid angle Eq. (9).

3.1. Mathematical method

Let f be an integrable function of radius, r, mass, M, and sine of the emission angle, sinα, that is, f = f(r,M,sinα) and I the following elliptic integral I=rirf1f(r,M,sinα)dr.\begin{equation} \label{INT} I=\int_{r_i}^{r_f} \frac{1}{\sqrt{f(r,M,\sin\alpha)}}\ {\rm d}r. \end{equation}(10)We are interested in deriving a polynomial approximation of the elliptic integral I. We first define sinα = g(z), where g(z) is a generic function of z(α). To expand Eq. (10) in Taylor series we assume that α is very small2 and aim at obtaining an integrable polynomial function I=rirf1f(r,M,g(z))drP(rf,ri,M,g(z)).\begin{equation} \label{POL} I=\int_{r_i}^{r_f} \frac{1}{\sqrt{f(r,M,g(z))}}\ {\rm d}r\approx P(r_f,r_i,M,g(z)). \end{equation}(11)P contains even powers of g(z), since f(r,M,g(z)) ∝ g(z)2. This condition is given by substituting b=(rsinα)/(1rs/r)\hbox{$b = (r\sin\alpha)/(\!\sqrt{1-r_{\rm s}/r})$} in the equations of the light bending Eq. (2), time delay Eq. (5), and solid angle Eq. (9). For an exact polynomial approximation, we therefore define g(z)=Az2+Bz\hbox{$g(z)=\sqrt{Az^2+Bz}$}, where A and B are general parameters. One of the two parameters (A,B) is determined by comparing Eq. (11) with the original integral I for special values of M = M, rf = rf and ri = ri that permits solving the integral I easily and obtain I(rf,ri,M,sinα)=P(rf,ri,M,Az2+Bz).\begin{equation} \label{EQL} I(r_f*,r_i*,M*,\sin\alpha)=P(r_f*,r_i*,M*,\sqrt{Az^2+Bz}). \end{equation}(12)The other parameter can be determined through the initial condition sinα=Az2+Bz\hbox{$\sin\alpha=\sqrt{Az^2+Bz}$}. We note that the polynomial approximation is valid for any emission angle α (not only for low values) since the parameters A,B are gauged on the whole range of I.

3.2. Light bending

For the light bending we Taylor-expand Eq. (2) up to the third order and defining u = 2M/R and sinα = g(z) we obtain ψbR[1+g2(z)6(1u)g2(z)u8(1u)+3g4(z)40(1u)2+3g4(z)u256(1u)2g4(z)8(1u)2u+5g6(z)112(1u)3g6(z)u332(1u)315g6(z)128(1u)3u+5g6(z)48(1u)3u2]·\begin{eqnarray} \label{eq:lb} \psi&\approx &\frac{b}{R} \left[1+\frac{g^2(z)}{6(1-u)}-\frac{g^2(z)u}{8(1-u)}+\frac{3g^4(z)}{40(1-u)^2} \right. \nonumber\\ &&\quad +\frac{3g^4(z)u^2}{56(1-u)^2}-\frac{g^4(z)}{8(1-u)^2}u+\frac{5g^6(z)}{112(1-u)^3} \nonumber\\ &&\quad \left. -\frac{g^6(z)u^3}{32(1-u)^3}-\frac{15g^6(z)}{128(1-u)^3}u+\frac{5g^6(z)}{48(1-u)^3}u^2\right]\cdot \end{eqnarray}(13)Setting g(z)=Az2+Bz\hbox{$g(z)=\sqrt{Az^2+Bz}$} and neglecting all the terms up to the second order in z, Eq. (13) becomes ψAz2+Bz1u[1+(B6(1u)Bu8(1u))z]·\begin{equation} \label{TEF} \psi\approx\sqrt{\frac{Az^2+Bz}{1-u}}\left[1+\left(\frac{B}{6(1-u)}-\frac{Bu}{8(1-u)}\right)z \right] \cdot \end{equation}(14)To approximate this equation with a polynomial, we introduce an even trigonometric function of ψ to remove the square root. The simplest choice is a cosine function expanded to the fourth order in ψ1cosψψ22ψ424Bz2(1u)+[B26(1u)2B2u8(1u)2+A2(1u)B224(1u)2]z2,\begin{eqnarray} \label{COSESP} 1-\cos\psi&\approx& \frac{\psi^2}{2}-\frac{\psi^4}{24}\nonumber\\ &\approx& \frac{Bz}{2(1-u)}+\left [\frac{B^2}{6(1-u)^2}-\frac{B^2u}{8(1-u)^2} \right .\nonumber\\ &&\quad\left. +\frac{A}{2(1-u)}-\frac{B^2}{24(1-u)^2} \right]z^2, \end{eqnarray}(15)where we consider the terms to the second order in z. If we choose A = −(B/ 2)2 , we obtain a simple linear approximation, 1−cosψBz/ (2(1−u)), in which z2 coefficients vanish.

We now solve Eq. (2) for the special values u = 0,R = 1 and obtain ψ=b1drr2[1sin2αr2]12=α.\begin{equation} \label{TELB} \psi=b\int_1^\infty \frac{{\rm d}r}{r^2}\left[1-\frac{\sin^2\alpha}{r^2} \right]^{-\frac{1}{2}}=\alpha. \end{equation}(16)Using the same values (u = 0,R = 1) for the approximated polynomial equation, 1−cosψBz/ (2(1−u)), we obtain 1cosα=Bz2·\begin{equation} 1-\cos\alpha=\frac{Bz}{2}\cdot \end{equation}(17)In this case, by defining B = 2 (implying A = −1), we find z = 1−cosα, which, when replaced in Eq. (15), gives the approximate light bending equation originally found by Beloborodov (2002)1cosψ=(1cosα)(1u)·\begin{equation} \label{AFLB} 1-\cos\psi=\frac{(1-\cos\alpha)}{(1-u)}\cdot \end{equation}(18)In Fig. (3) we show a comparison between the exact light bending curves for different emission radii, and curves obtained from the approximate equation. The accuracy of the latter between 0 ≤ ααmax is better than 3% for R = 3rs, while for R = 5rs the error does not exceed 1%. We note that R = 3rs corresponds to the innermost stable circular orbit (ISCO) for matter orbiting a Schwarzschild BH and is also close represent to a typical NS radius size of ~ 12 km for mass of 1.4 M. For values below R = 2rs the equation is not anymore applicable after α = π/ 2. In Fig. 3 we also show the exact light bending curve for R = 1.55rs; after a given minimum the photons are highly bent by strong-field effects. The largest error is at α = π/ 2 and then it tends to decrease until at αmax because of the symmetrization process around α = π/ 2 configuring as the maximum reachable angle (see Sect. 2.2). For more details about the accuracy between 0 ≤ απ/ 2 we refer to Beloborodov (2002).

thumbnail Fig. 3

Light bending curves from the exact Eq. (2) (solid lines), compared to those from the approximate Eq. (18) (dashed red lines) for R = 2rs, R = 3rs, and R = 5rs. The dotted blue line represents the threshold from trajectories without a turning point (0 ≤ απ/ 2) to trajectories with a turning point (π/ 2 ≤ ααmax). The exact light bending curve for R = 1.55rs is also plotted (dotted-dashed orange line) to show strong-field effects. The lower panels show the difference between the curves from the original and approximate equations.

3.3. Time delay

We now apply our method for deriving the approximate equation for the time delay. By expanding the integrand in Eq. (5) up to the third order Δt=R{g2(z)2(1u)+g4(z)8(1u)23g4(z)32(1u)2u+g6(z)16(1u)35g6(z)48(1u)3u+5g6(z)112(1u)3u2},\begin{eqnarray} \label{TDGE2} \Delta t &=& R \left \{\frac{g^2(z)}{2(1-u)}+\frac{g^4(z)}{8(1-u)^2}-\frac{3g^4(z)}{32(1-u)^2}u \right. \nonumber\\ &&\quad \left. +\frac{g^6(z)}{16(1-u)^3}-\frac{5g^6(z)}{48(1-u)^3}u+\frac{5g^6(z)}{112(1-u)^3}u^2\right \}, \end{eqnarray}(19)we set again g(z)=Az2+Bz\hbox{$g(z)=\sqrt{Az^2+Bz}$} and neglect all terms up to the third order in z, so that ΔtR{Az2+Bz2(1u)+B2z2+2ABz38(1u)23u(B2z2+2ABz3)32(1u)2+B3z316(1u)35uB3z348(1u)3+5u2B3z3112(1u)3}·\begin{eqnarray} \label{TDGE3} \frac{\Delta t}{R} &\approx & \left \{\frac{Az^2+Bz}{2(1-u)}+\frac{B^2z^2+2ABz^3}{8(1-u)^2}-\frac{3u(B^2z^2+2ABz^3)}{32(1-u)^2} \right. \nonumber\\ &&\quad \left. +\frac{B^3z^3}{16(1-u)^3}-\frac{5uB^3z^3}{48(1-u)^3}+\frac{5u^2B^3z^3}{112(1-u)^3}\right \}\cdot \end{eqnarray}(20)To determine (A,B) we compare the original Eq. (5) with Eq. (20) both evaluated for u = 0 and R = 1; we find31cosα=B2z+12(A+B24)z2+B4(A+B24)z3,\begin{equation} \label{AF} 1-\cos\alpha=\frac{B}{2}z+\frac{1}{2}\left(A+\frac{B^2}{4}\right)z^2+\frac{B}{4}\left(A+\frac{B^2}{4}\right)z^3, \end{equation}(21)where on the left and right hand sides are the results of Eqs. (5) and (20), respectively. By imposing A + B2/ 4 = 0 the coefficients of the second and third order in z vanish. Like in the light bending case, Eq. (21) reduces to 1−cosα = Bz/ 2; defining again B = 2 (implying A = −1) substituting in Eq. (20), we derive the approximate travel time delay equation (see for further details Poutanen & Beloborodov 2006) ΔtR=y[1+uy8+uy224u2y2112],\begin{equation} \label{FATIDE} \frac{\Delta t}{R}=y\left [1+\frac{uy}{8}+\frac{uy^2}{24}-\frac{u^2y^2}{112} \right], \end{equation}(22)where y = (1−cosψ).

In Fig. 4 we compare for different emission radii the exact travel time delay curves with the polynomial approximated equations. We here also extend the validity of the approximation to αmax-values accounting for turning points. The accuracy settles ~ 35% for R = 2rs, while after R = 3rs it is lower than 20%, according to the same symmetry argument explained in the Sect. 3.2. However, we refer to Poutanen & Beloborodov (2006) for the error estimation between 0 ≤ απ/ 2.

thumbnail Fig. 4

Continuous black curves are obtained from the original time delay Eq. (5), while the dashed red curves are from the polynomial approximate Eq. (22). Different panels show R = 2rs, R = 3rs and R = 5rs. The dotted blue line helps distinguishing trajectories without a turning point (i.e., 0 ≤ cosα ≤ 1) from those with a turning point (i.e., π/ 2 ≤ ααmax). The lower panels show the difference between the curves from the original and approximate equations.

3.4. Solid angle

We now apply the same method to derive for the first time a polynomial approximation to the solid angle Eq. (9). We note, at variance of light bending and time delay equations, that the solid angle equation has the integral in the denominator, and moreover, the emission angle, α, is also outside the integral. We first rewrite Eq. (9) as =P1P2Idrdϕ,\begin{equation} \label{INI} {\rm d}\Omega=\frac{P_{\rm 1}\ P_{\rm 2}}{I}\ {\rm d}r\ {\rm d}\varphi, \end{equation}(23)where P1=cosiD2sin2ψ(1u),P2=sin2αcosα,I=Rdrr2[1R2sin2αr2(1u)(1uRr)]32.\begin{eqnarray} \label{eq:soan} &&P_{\rm 1} =\frac{\cos i}{D^2\sin^2\psi\ (1-u)},\qquad P_{\rm 2}=\frac{\sin^2\alpha}{\cos\alpha}, \nonumber\\ && I=\int_R^\infty \frac{{\rm d}r}{r^2}\left[1-\frac{R^2\sin^2\alpha}{r^2(1-u)}\left(1-\frac{uR}{r} \right) \right]^{-\frac{3}{2}}. \end{eqnarray}(24)P1 is a constant because ψ is a function of the azimuthal angle, ϕ, the inclination angle, i, and the polar coordinate, θ, (for further details see Sect. 4). As a first step, we expand the integrand of I in a Taylor series up to the third order in z. We derive I1+Cz+Dz2R,\begin{equation} \label{AI} I\approx \frac{1+Cz+Dz^2}{R}, \end{equation}(25)with C=B2(1u)3Bu8(1u),D=A2(1u)3Au8(1u)+3B28(1u)2+15B2u216(1u)25B2u8(1u)2·\begin{eqnarray} C&=& \frac{B}{2(1-u)}-\frac{3Bu}{8(1-u)},\\ D&=& \frac{A}{2(1-u)}-\frac{3Au}{8(1-u)}+\frac{3B^2}{8(1-u)^2} \\ &&\quad + \frac{15B^2u^2}{16(1-u)^2}-\frac{5B^2u}{8(1-u)^2}\cdot \notag \end{eqnarray}The function P2/I is not yet a polynomial function since it contains a ratio of polynomials and square root functions in P2. For these reasons we expand P2/I in a Taylor series around z = 0 and neglect all the terms up to third order in zP2IAz2+Bz1Az2BzR1+Cz+Dz2R[Bz+(B22+ACB)z2+(AB+3B28CB22CA+BC2BD)z3].\begin{eqnarray} \label{AIA} \frac{P_{\rm 2}}{I}&\approx& \frac{Az^2+Bz}{\sqrt{1-Az^2-Bz}}\frac{R}{1+Cz+Dz^2}\nonumber\\ &\approx& R\left [Bz+\left(\frac{B^2}{2}+A-CB \right)z^2 \right.\nonumber\\ &&\quad\left. +\left(AB+\frac{3B^2}{8}-\frac{CB^2}{2}-CA+BC^2-BD \right)z^3 \right]. \end{eqnarray}(28)To determine (A,B) we compare the original solid angle Eq. (24) with the above approximate equation, evaluating both equations for u = 0 and R = 1; we find sin2α=Bz+Az2.\begin{equation} \sin^2 \alpha=Bz+Az^2. \end{equation}(29)The left- and right-hand sides are the result of original Eq. (24) and the polynomial Eq. (28), respectively. We can freely define the value of A and B because there are no particular constraints to impose. We set, as in the previous cases, A = −1 and B = 2, deriving again z = 1−cosα. The final approximate equation for the solid angle is cosiD2sin2ψ(1u)R[2z+(12C)z2+(1C+2C22D)z3]drdϕ,\begin{eqnarray} \label{AFSA} {\rm d}\Omega&\approx& \frac{\cos i}{D^2\sin^2\psi\ (1-u)}\ R \left[2z+\left(1-2C\right)z^2 \right. \nonumber\\[3mm] &&\quad \left. +\left(1-C+2C^2-2D\right)z^3\right]\ {\rm d}r\ {\rm d}\varphi, \end{eqnarray}(30)where C=43u4(1u),D=39u291u+5656(1u)2·\begin{equation} C=\frac{4-3u}{4(1-u)},\qquad D=\frac{39u^2-91u+56}{56(1-u)^2}\cdot \end{equation}(31)As for the previous two cases, in Fig. 5 we compare the exact solid angle curves with the polynomial approximated curves for different radii and inclination angles i. The comparison extends to αmax-values and thus accounts for trajectories with turning points in this case as well. For R = 3rs the error is ~ 5% and after R = 5rs it is lower than 1%. We note that for i = 30° the curves are fairly flat because the relativistic effects are small. Instead, passing from i = 60° to i = 80° , the curves become gradually steeper as general relativistic effects increase. Unlike the previous cases, we do not show here the case R = 2rs because the approximate formula Eq. (30) does not give adequately accurate results.

thumbnail Fig. 5

Continuous black curves are obtained from the original solid angle Eq. (9), while the dashed red curves are from the polynomial approximation in Eq. (30). Different panels are for R = 3rs and R = 5rs and three different inclination angles, i = 30°, i = 60°, and i = 80°. The dotted blue line helps distinguishing trajectories without a turning point (i.e., 0 ≤ cosα ≤ 1) from those with a turning point (i.e., π/ 2 ≤ ααmax). The lower panels show the difference between the curves from the original and approximate equations. The dotted-dashed lines represent the difference between the curves from the original and Beloborodov (2002) equations.

We note that Eq. (A3) in Beloborodov (2002) is obtained by approximating the derivative dcosψdcosα\hbox{$\frac{{\rm d}\!\cos\psi}{{\rm d}\!\cos\alpha}$} with the linear Eq. (18), while our Eq. (30) is a third-order polynomial that approximates the integral I and all the terms depending on the emission angle α. For example, our approximation is more accurate by a factor of ~3 to 10 for R = 3rs and 0 ≤ cosα ≤ 0.3.

thumbnail Fig. 6

Geometries adopted in the examples. Left: emission from a disk or clump orbiting a Schwarzschild BH. Right: emission from two hot opposite spots on an NS surface.

4. Examples of astrophysical applications

In this section we present three simple examples of astrophysical applications of the approximate equations. We consider the emission point at coordinates (r,ϕ,θ). The observer is located at infinity along the z-axis with a viewing angle, i, with respect to the z-axis; the observer polar coordinates are (r′,ϕ′,θ). Photons emitted from a point are deflected by an angle, ψ, and reach the observer with impact parameter, b. The plane containing the photon trajectory rotates around the line of sight as the emission point moves around the compact object. Two unit vectors are attached to the photon emission point, E: u is tangential to the photon trajectory, and n points in the same direction as the radius, R. The photon deflection angle, ψ, varies as cosψ=sinisinθcosϕ+cosicosθ,\begin{equation} \cos\psi=\sin i\sin\theta \cos\varphi+\cos i \cos\theta, \end{equation}(32)with θ = π/ 2, ϕ = ωkt and t = 0 when the emission point is closest to the observer. The photon arrival time, Tobs, is the sum of the emission time, Torb = ϕ/ωk, plus the photon propagation delay, ΔT(b), from the emission point to the observer (see Eq. (5)).

The observed flux is F=νobsΩIνobsdνobs\hbox{$F=\int_{\nu_{\rm obs}}\int_\Omega I_{\nu_{\rm obs}}{\rm d}\Omega\ {\rm d}\nu_{\rm obs}$}, where Iνobs is the specific intensity at the photon frequency νobs. We use the Lorentz invariant ratio Iνobs/νobs3=Iνem/νem3\hbox{$I_{\nu_{\rm obs}}/\nu^3_{\rm obs}=I_{\nu _{\rm em}}/\nu^3 _{\rm em}$} (see, e.g., Misner et al. 1973), where Iνem(ξ)=ϵ0ξq4πδ(νemνobs)\hbox{$I_{\nu _{\rm em}}(\xi)=\frac{\epsilon_0 \xi^{q}}{4\pi}\delta(\nu _{\rm em}-\nu_{\rm obs})$} is the specific intensity at the emission point E given by the product of the surface emissivity, varying as a power law of ξ = R/M with index q, and the delta function peaked at νem. Therefore, integrating over all the frequencies, we obtain the observed flux at frequency νem, Fνem=Ωϵ0ξq4π(1+z)-4\hbox{$F_{\nu _{\rm em}}=\int_{\Omega} \frac{\epsilon_0 \xi^{q}}{4\pi}\,(1+z)^{-4}\ {\rm d}\Omega$}. The redshift is defined as the ratio between the observed and the emitted energy, (1 + z)-1 = νobs/νem (Misner et al. 1973) and for matter orbiting in circular orbits around a compact object or for a spot on a NS surface reads as (1+z)-1=(1rsRω2R2sin2θ)1/2(1+sinisinϕsinθsinψ)-1·\begin{equation} \label{Redshift2} (1+z)^{-1}=\left(1-\frac{r_{\rm s}}{R}-\omega^{2}R^{2}\sin^{2}\theta\right)^{1/2}\left(1+b\omega \frac{\sin i\,\sin\varphi\sin\theta}{\sin\psi} \right)^{-1}\cdot \end{equation}(33)For ω = ωk we consider matter orbiting with Keplerian velocity around a BH, and for ω = ωspin we consider spots rotating with the NS spin frequency. The relevant geometry is shown in Fig. 6.

4.1. Light curve from an emitting clump orbiting a black hole

We first consider a clump defined as a small sphere radiating isotropically in its own rest frame, orbiting a Schwarzschild BH in a circular orbit with angular velocity ωk = (M/R3)1/2. The geometry is shown in Fig. 6. For simplicity we assume ϵ0ξq4π=1\hbox{$\frac{\epsilon_0 \xi^{q}}{4\pi}=1$}. Figure 7 shows the modulation of the Doppler factor (1 + z)-4, solid angle dΩ, and flux from the orbiting clump as a function of phase, ϕ(t), including light travel time delays. When the clump is behind the BH, gravitational lensing magnifies the solid angle from which the clump is seen by observer; the Doppler factor is greatest when the projected velocity along the photon trajectory reaching the observer is highest. The gravitational effects are stronger for larger inclination angles, and the observed peak flux is not at ϕ = π, but is significantly shifted especially for large inclination angles due to the travel time delays. The errors between the approximated and the original flux depend only on the emission radius, since the inclination angle figures as a constant. However, it is evident that the main errors derive from the approximated time delay equation (as shown in the Sect. 3.3).

thumbnail Fig. 7

Modulated flux (normalized to the maximum), Doppler factor (1 + z)-4 and solid angle (arbitrary units) in the rest coordinate frame of an emitting clump in a circular orbit around a Schwarzschild BH for different radii and inclinations angles. The continuous black lines are calculated with the exact equations, while the dashed red lines are calculated with the approximate equations. All quantities are plotted as a function of the arrival phase at the observer. In the left panels a self-eclipse of the spot is apparent.

4.2. Emission line profile from an accretion disk around a black hole

In Fig. 8 we calculate the steady relativistically broadened emission line profile from an accretion disk around a Schwarzschild BH (e.g., Fabian et al. 1989; Beckwith & Done 2004, and references therein). Fe lines at ~ 6−7 keV from a number of accreting stellar mass BHs and NSs in X-ray binaries, as well as supermassive BHs in the nuclei of active galaxies are interpreted based on this model (e.g., Tomsick et al. 2014). We integrate over the disk surface from an inner to an outer disk radius and ignore light propagation delays, as we consider a steady disk. The approximate equations reproduce very accurately the profiles obtained with the exact equations. A high accuracy is also retained for large inclination angles, even if larger inclination angles enhance the relativistic effects (see Sect. 4.1).

thumbnail Fig. 8

Line profile for isotropic radiation from Rin = 3 rs to Rout = 50 rs assuming surface emissivity q = −3. The continuous black lines represent the original equations and the dashed red lines are the polynomial approximate equations.

4.3. Light curve from a hot spot on the surface of a rotating neutron star

We calculate here the pulse profile generated by a point-like hot spot located on the surface of a NS, which emits like an isotropic blackbody. Calculations of this type have been carried out extensively to model the periodic signals of accreting millisecond pulsars (see, e.g., Pechenick et al. 1983; Poutanen & Beloborodov 2006; Leahy et al. 2011; Bauböck et al. 2015) as well as the so-called burst oscillations during Type I thermonuclear bursts in NS low-mass X-ray binaries (e.g., Nath et al. 2002; Miller & Lamb 2015); some of these calculations also include the angular size of the hot spot, the star oblateness, and the spacetime modifications induced by fast rotation. We use here a canonical NS mass of 1.4 M and radius RNS = 12 km, together different inclination angles, i, and colatitudes, θ of the spot. The NS spin frequency is chosen to be νs = 600 Hz. In Fig. 9 we report the corresponding pulse profiles; as expected, the case with higher values of i and θ displays larger departures from a sinusoidal shape. In this type of applications the value of α is always limited to π/ 2, as no turning points are involved. Therefore our approximate equations retain very high accuracy as long as the NS radius is ≥ 2.5rs, a range that encompasses a number of NS models for different equations of state, excluding only the upper end of the mass-radius branches. We conclude that our approximate equations can be usefully employed in calculations of the pulse profile of fast spinning NSs over a range of (but not all) models to be tested against the observation that the Neutron Star Interior Composition ExploreR (NICER), and other large-area X-ray missions of the future, such as Athena or LOFT, will obtain (see Watts et al. 2016, and references therein).

thumbnail Fig. 9

Modulation from a hot spot on an NS as a function of rotational phase for different inclination angles and hot spot colatitude. Light travel time delays are included. The continuous black lines represent the results from a numerical integration of the original equations; the dashed red lines are obtained from the polynomial approximate equations. The dashed-dotted orange line does not include light travel time delays.

4.4. Applicability regions

In Fig. 10 we plot ψmax as a function of the emission radius to investigate the applicability regions of the approximate equations. If we consider trajectories with turning points for radii R< 3rs, that is, smaller than the ISCO, then ψmax ≥ 180° and a polynomial treatment is no longer accurate because of strong field effects (see also Fig. 3). We note that for R− → 1.5rs, ψmax our solution approaches asymptotically 270°. Instead, for R ≥ 3rs, when the observer is located edge on (i.e., i = 90°), ψmax = 180° is attained; otherwise, for slightly smaller but still extreme inclination angles, for example, 87°, photon trajectories always remain below the critical bending angle, which guarantees a high accuracy of our polynomial approximations. This argument is valid for all the emission radii R ≥ 3rs, since for R− → ∞, ψmax approaches 180°.

thumbnail Fig. 10

Largest bending angle ψmax, vs. the emission radius (continuous black line). For inclination angles below i = 87° (dashed red line) the approximate equations provide a high accuracy, since they are below the ψmax-value. R = 3rs (dotted blue line) separates the applicability region from the strong-field regime (R< 3rs).

5. Conclusions

We developed an analytical method to approximate the elliptic integrals that describe gravitational light bending and light travel time delays of photon geodesics in the Schwarzschild metric. Based on this, we derived for the first time an approximate polynomial equation also for the solid angle. We discussed the accuracy and range of applicability of the approximate Eqs. (18), (22), and (30); adopting them can considerably speed up calculations related to a variety astrophysical problems, which normally require time-consuming numerical integrations. We also presented a few simple applications as examples. We will extend our treatment to the parallel transport of polarization vectors in a future work.


1

Equation (9) is equivalent to the formula (A3) in Beloborodov (2002).

2

Therefore, g(z) is small as well.

3

For Eq. (5) we used the following limit: limx+(x2a)12x=0\hbox{$\lim_{x\to +\infty}(x^2-a)^{\frac{1}{2}}-x=0$}.

Acknowledgments

This research was financed by the Swiss National Science Foundation project 200021_149865. V.d.F. and M.F. acknowledge the Department of Physics at the University of Basel and especially Friedrich-K. Thielemann. We also thank the International Space Science Institute in Bern for their support. V.d.F. is grateful to the International Space Science Institute in Beijing for the hospitality to carry out part of this work. L.S. acknowledges partial support under contract ASI INAF I/004/11/1.

References

  1. Bao, G., Hadrava, P., & Ostgaard, E. 1994, ApJ , 435, 55 [NASA ADS] [CrossRef] [Google Scholar]
  2. Bauböck, M., Psaltis, D., & Özel, F. 2015, ApJ , 811, 144 [NASA ADS] [CrossRef] [Google Scholar]
  3. Beckwith, K., & Done, C. 2004, MNRAS , 352, 353 [NASA ADS] [CrossRef] [Google Scholar]
  4. Beloborodov, A. M. 2002, ApJ, 566, L85 [NASA ADS] [CrossRef] [Google Scholar]
  5. Chandrasekhar, S. 1992, The mathematical theory of black holes (New York: Oxford University Press) [Google Scholar]
  6. Fabian, A. C., Rees, M. J., Stella, L., & White, N. E. 1989, MNRAS , 238, 729 [NASA ADS] [CrossRef] [Google Scholar]
  7. Leahy, D. A., Morsink, S. M., & Chou, Y. 2011, ApJ , 742, 17 [NASA ADS] [CrossRef] [Google Scholar]
  8. Luminet, J.-P. 1979, A&A , 75, 228 [NASA ADS] [Google Scholar]
  9. Miller, M. C., & Lamb, F. K. 2015, ApJ , 808, 31 [NASA ADS] [CrossRef] [Google Scholar]
  10. Misner, C. W., Thorne, K. S., & Wheeler, J. A. 1973, Gravitation (San Francisco: W.H. Freeman and Co) [Google Scholar]
  11. Nath, N. R., Strohmayer, T. E., & Swank, J. H. 2002, ApJ , 564, 353 [NASA ADS] [CrossRef] [Google Scholar]
  12. Pechenick, K. R., Ftaclas, C., & Cohen, J. M. 1983, ApJ , 274, 846 [NASA ADS] [CrossRef] [Google Scholar]
  13. Poutanen, J., & Beloborodov, A. M. 2006, MNRAS , 373, 836 [NASA ADS] [CrossRef] [Google Scholar]
  14. Tomsick, J. A., Nowak, M. A., Parker, M., et al. 2014, ApJ , 780, 78 [NASA ADS] [CrossRef] [Google Scholar]
  15. Watts, A. L., Andersson, N., Chakrabarty, D., et al. 2016, Rev. Mod. Phys., 88, 021001 [NASA ADS] [CrossRef] [Google Scholar]

All Figures

thumbnail Fig. 1

Two photon trajectories emitted at different radii, r1 and r2, and emission angles, α1 and α2, with their corresponding impact parameters, b1 and b2. Trajectory 1 is for a direct photon, while trajectory 2 has a turning point, i.e., passes through periastron p, the minimum distance between the trajectory and the BH. The observer is at infinity, and photons geodesics lie in a single invariant plane.

In the text
thumbnail Fig. 2

Calculation of travel time delay for trajectories with turning points. The photon is emitted at E with radius R and deflection angle ψE. The photon trajectory passes through point, S, which is symmetric to E with respect to periastron p, having the same impact parameter, b, and a deflection angle ψS = 2ψpψE.

In the text
thumbnail Fig. 3

Light bending curves from the exact Eq. (2) (solid lines), compared to those from the approximate Eq. (18) (dashed red lines) for R = 2rs, R = 3rs, and R = 5rs. The dotted blue line represents the threshold from trajectories without a turning point (0 ≤ απ/ 2) to trajectories with a turning point (π/ 2 ≤ ααmax). The exact light bending curve for R = 1.55rs is also plotted (dotted-dashed orange line) to show strong-field effects. The lower panels show the difference between the curves from the original and approximate equations.

In the text
thumbnail Fig. 4

Continuous black curves are obtained from the original time delay Eq. (5), while the dashed red curves are from the polynomial approximate Eq. (22). Different panels show R = 2rs, R = 3rs and R = 5rs. The dotted blue line helps distinguishing trajectories without a turning point (i.e., 0 ≤ cosα ≤ 1) from those with a turning point (i.e., π/ 2 ≤ ααmax). The lower panels show the difference between the curves from the original and approximate equations.

In the text
thumbnail Fig. 5

Continuous black curves are obtained from the original solid angle Eq. (9), while the dashed red curves are from the polynomial approximation in Eq. (30). Different panels are for R = 3rs and R = 5rs and three different inclination angles, i = 30°, i = 60°, and i = 80°. The dotted blue line helps distinguishing trajectories without a turning point (i.e., 0 ≤ cosα ≤ 1) from those with a turning point (i.e., π/ 2 ≤ ααmax). The lower panels show the difference between the curves from the original and approximate equations. The dotted-dashed lines represent the difference between the curves from the original and Beloborodov (2002) equations.

In the text
thumbnail Fig. 6

Geometries adopted in the examples. Left: emission from a disk or clump orbiting a Schwarzschild BH. Right: emission from two hot opposite spots on an NS surface.

In the text
thumbnail Fig. 7

Modulated flux (normalized to the maximum), Doppler factor (1 + z)-4 and solid angle (arbitrary units) in the rest coordinate frame of an emitting clump in a circular orbit around a Schwarzschild BH for different radii and inclinations angles. The continuous black lines are calculated with the exact equations, while the dashed red lines are calculated with the approximate equations. All quantities are plotted as a function of the arrival phase at the observer. In the left panels a self-eclipse of the spot is apparent.

In the text
thumbnail Fig. 8

Line profile for isotropic radiation from Rin = 3 rs to Rout = 50 rs assuming surface emissivity q = −3. The continuous black lines represent the original equations and the dashed red lines are the polynomial approximate equations.

In the text
thumbnail Fig. 9

Modulation from a hot spot on an NS as a function of rotational phase for different inclination angles and hot spot colatitude. Light travel time delays are included. The continuous black lines represent the results from a numerical integration of the original equations; the dashed red lines are obtained from the polynomial approximate equations. The dashed-dotted orange line does not include light travel time delays.

In the text
thumbnail Fig. 10

Largest bending angle ψmax, vs. the emission radius (continuous black line). For inclination angles below i = 87° (dashed red line) the approximate equations provide a high accuracy, since they are below the ψmax-value. R = 3rs (dotted blue line) separates the applicability region from the strong-field regime (R< 3rs).

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.