© 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 (r_s = 2GM/c²). 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 (1)where G = c = 1, and r_s = 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): (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 (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).

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 (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 b_c into Eq. (3) (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 (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 p³−b²p + b²r_s = 0. The polynomial in p has three real solutions (because b ≥ b_c): one is negative, one is lower than 3M, and we consider only the solution satisfying p ≥ r_c, where r_c = 3M is the critical radius associated to b_c (see, e.g., Luminet 1979). The time delay is composed of the time delay Δt_S from point α_S, as determined by the Eq. (5), plus the time delay between [α_S,α_p], Δt_p−S, and [α_p,α_E], Δt_E−p. Since the integrand is symmetric with respect to α_p, the latter two time delays are equal (Δt_E−p = Δt_p−S), the equation can be written (see Fig. 2) (6)

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, dΩ, 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 b ≈ D·ψ, where D is the distance from the emission point to the observer, (7)In the emission reference frame, Eq. (7) becomes (8)where we considered the following dependencies ϕ = ϕ(ϕ′) and b = b(r,ψ). The Jacobian of the transformation is always independent of the value of , 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. 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)¹(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 (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 small² and aim at obtaining an integrable polynomial function (11)P contains even powers of g(z), since f(r,M,g(z)) ∝ g(z)². This condition is given by substituting 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 , 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 ∗, r_f = r_f ∗ and r_i = r_i ∗ that permits solving the integral I easily and obtain (12)The other parameter can be determined through the initial condition . 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 (13)Setting and neglecting all the terms up to the second order in z, Eq. (13) becomes (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 ψ(15)where we consider the terms to the second order in z. If we choose A = −(B/ 2)² , we obtain a simple linear approximation, 1−cosψ ≈ Bz/ (2(1−u)), in which z² coefficients vanish.

We now solve Eq. (2) for the special values u = 0,R = 1 and obtain (16)Using the same values (u = 0,R = 1) for the approximated polynomial equation, 1−cosψ ≈ Bz/ (2(1−u)), we obtain (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)(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 = 3r_s, while for R = 5r_s the error does not exceed 1%. We note that R = 3r_s 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 = 2r_s the equation is not anymore applicable after α = π/ 2. In Fig. 3 we also show the exact light bending curve for R = 1.55r_s; 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).

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 (19)we set again and neglect all terms up to the third order in z, so that (20)To determine (A,B) we compare the original Eq. (5) with Eq. (20) both evaluated for u = 0 and R = 1; we find³(21)where on the left and right hand sides are the results of Eqs. (5) and (20), respectively. By imposing A + B²/ 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) (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 = 2r_s, while after R = 3r_s 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.

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 (23)where (24)P₁ 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 (25)with The function P₂/I is not yet a polynomial function since it contains a ratio of polynomials and square root functions in P₂. For these reasons we expand P₂/I in a Taylor series around z = 0 and neglect all the terms up to third order in z(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 (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 (30)where (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 = 3r_s the error is ~ 5% and after R = 5r_s 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 = 2r_s because the approximate formula Eq. (30) does not give adequately accurate results.

We note that Eq. (A3) in Beloborodov (2002) is obtained by approximating the derivative 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 = 3r_s and 0 ≤ cosα ≤ 0.3.

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 (32)with θ = π/ 2, ϕ = ω_kt and t = 0 when the emission point is closest to the observer. The photon arrival time, T_obs, is the sum of the emission time, T_orb = ϕ/ω_k, plus the photon propagation delay, ΔT(b), from the emission point to the observer (see Eq. (5)).

The observed flux is , where I_νobs is the specific intensity at the photon frequency ν_obs. We use the Lorentz invariant ratio (see, e.g., Misner et al. 1973), where 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, . 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 (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/R³)^1/2. The geometry is shown in Fig. 6. For simplicity we assume . 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).

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

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 R_NS = 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.5r_s, 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).

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< 3r_s, 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.5r_s, ψ_max our solution approaches asymptotically 270°. Instead, for R ≥ 3r_s, 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 ≥ 3r_s, since for R− → ∞, ψ_max approaches 180°.

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.

---

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

All Figures