{'fraction_non_alphanumeric': 0.038800822444783445, 'fraction_numerical': 0.01949990051071168, 'mean_word_length': 4.109454422229753, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Optics in the Schwarzschild space-time

Realistic modelling of radiation transfer in and from variable accretion disks around black holes requires the solution of the problem: find the constants of motion and equation of motion of a light-like geodesic connecting two arbitrary points in space. Here we give the complete solution of this problem in the Schwarzschild space-time.


Introduction
The first light detected from a relativistic region about a black hole was discovered by the ASCA satellite (Tanaka et al. 1995, Yaqoob et al. 1998. The now accepted theoretical model describing the broad X-ray emission lines is that of an accretion disk around either a Kerr or a Schwarzschild black hole (Reynolds et al. 1995, Fanton et al. 1997, Bromley et al. 1997, Dabrowski & Lasenby 2001,Čadež et al. 1998, Martocchia et al. 2000. This discovery increased the interest for phenomena occurring in the vicinity of black holes. We now know that other interesting high energy phenomena, such as X-ray flares (Baganoff et al. 2001, Goldwurm et al. 2003 or quasi-periodic oscillations (Strohmayer 2001, Miller et al. 2002 are occuring in the environment of black holes. To further the understanding of such phenomena from the theoretical point of view, it is important to develop tools to model the phenomena themselves as well as to model radiation transfer in and from these strongly curved regions of space-time. Most current accretion disk models are very simple from the radiation transfer point of view. They consider disks as geometrically thin and optically thick. Light that reaches a far observer, comes directly (without beeing scattered) from a very definite point on the disk. Therefore, line profiles can be calculated by aiming light-like goedesic from a point on the disk to the observer or vice versa, aiming from the observer to points on the disk. In such ray-tracing procedures geodesic equations are usually solved by direct numerical integration (Laor 1991, Pariev et al. 2001, Reynolds & Begelman 1997. However, when modelling transient phenomena produced by small debris around black holes, or by other transient phenomena in the disk (moving hot spots, varying external illumination, waves), it becomes necessary to solve a more difficult radiation transfer problem, the problem of following a single photon through its more than one scattering and/or more than one possible path from the source to the eye of the observer (Schnittman 2005, Schnittman & Bertschinger 2004. The problem is further complicated by noting that photon arrival times from the same initial source to the observer may (and often will) be markedly different for photons reaching the observer along different possible paths. It is obvious, that a multiple scattering path cannot be effectively constructed by aiming geodesics between succesive scattering points as the number of succesive iterations required would soon blow up. The solution to the radiation transfer problem thus requires one to be able to follow a light ray from one to the other scattering point along its path. Thus, one would like to find the quickest way to determine all constants of motion of a geodesic that connects the given initial and final point. In this article we give analytic expressions which completely solve this problem in the Schwarzschild space-time. We also turn this analytic tool into a numerical code and demonstrate that it is much faster, more accurate and more transparent than aiming and integrating geodesic equations. This tool thus opens the possibility to solve complex radiation transfer problems in curved space-time using similar Monte Carlo techniques that are used in solving radiation problems in flat space-time.
Our work starts with the results of Chandrasekhar (1992) and Rauch & Blandford (1994), who expressed the solutions to geodesic equations in terms of elliptic integrals. By inverting their expressions into Jacobi elliptic functions (Čadež et al. 1998), we obtain simple solutions for all the three types of orbit equations that occur in Schwarzscild space-time. These solutions no longer contain branch ambiguities. Since these orbits are essentialy planar, their equations depend only on two nontrivial constants of motion: the angular momentum and the longitude of the periastron. Expressing the orbit equation at the initial and final point on the geodesic, one obtains two nonlinear equations for the two nontrivial constants of motion. However, since the longitude of the periastron occurs only linearly as the argument of elliptic functions, it is possible to use the elliptic functions addition theorem to elliminate the longitude of the periastron and obtain a single nonlinear equation for the angular momentum as function of initial and final coordinates. Here we derive these equations for all the three types of orbits and discuss their properties and solutions. We also write down all the other constants of motion in terms of final and initial point coordinates and give analytic expressions for travel times.

Connecting two points with a light-like geodesic
In the Schwarzschild space-time it is customary to introduce Schwarzscild coordinates t, r, θ and ϕ. In these coordinates geodesics are governed by the Hamiltonian: which admits 8 constants of motion: the value of the Hamiltonian (H) and the value of the Lagrangian (L) (the ratio of the two defining the relation between time and proper time), the energy E = p t , the angular momentum ( l), the longitude of the periastron (ω) and the time of periastron passage. The angular momentum is expressed as l = l ·n, wheren is the unit vector along l, which is defined by the inclination of the orbit (ε) and the longitude of the ascending node (Ω).
Consider a light-like geodesic connecting the initial point P i and the final point P f . Five constants H = 0, L = 0, E, Ω and ε are determined readily. So one can use the true anomaly λ as a parameter along the geodesic. From the initial to the final point λ increases by: and the winding number k = ... − 1, 0, 1... tells how many times a geodesic winds around the black hole. 1 The two constants Ω and ε are obtained from angular coordinates of the initial and final point (see Figure 1) with the help of basic spherical trigonometry: tan Since the variable λ is conjugate to the orbital angular momentum, it obeys the Poisson bracket relation [λ, l] = 1, so that the angles θ and ϕ are expressed with λ as follows (Čadež & Gomboc 1996, eqs. (6 -16)): and 2 the differential equation for orbits of light-like geodesics becomes: Here u = 2M/r and a = 2ME/l. The solutions to this equation, called orbit equations, depend only on the parameter a and are of three types ( Figure 2): type A: scattering orbits with both endpoints at infinity; their angular momentum parameter is on the interval 0 < a < 2 3 √ 3 . Scattering orbits can never extend below r = 3M.
type B: plunging orbits with one end at infinity and the other behind the horizon; a > 2 3 √ 3 , type C : near orbits with both ends behind the horizon of the black hole; their angular momentum parameter is on the interval 0 < a < 2 3 √ 3 . Near orbits can never reach beyond r = 3M.
For completeness, we list the solutions of the orbit equation (8) in terms of a few auxiliary parameters.
Type A: Introduce a = 2 3 √ 3 sin ψ 2 (0 < ψ < π) and the following auxiliary parameters, functions of ψ only: where u i = 2M/r i and similarly u f = 2M/r f . In terms of these, the type A differential equation for the orbit becomes (Figure 2, left): with ∆λ = λ−λ i and ∆λ/n ∈ Figure 4 left. Here F and K are the incomplete and complete elliptic integrals of the first kind (Wolfram 1996).
Critical C : A special limiting case of solution (13) for µ = 0 (or (15) for ψ = π 2 ) (a = 2 3 √ 3 ) and u i > 2/3 (r i < 3M) is: with ∆λ = λ − λ i and ∆λ ∈ 2Arcoth 4/3 − 2Arcoth u i + 1/3 , ∞ . The two constants of motion ω and a can now be determined in two steps. First we determine the type of the orbit following reasoning illustrated in Figure 3 and then we solve the two equations obtained from the orbit equation at the initial and the final point (Brajnik 1999, Gomboc 2001, Cadež et al. 2003, Kostić 2003. We note that the parameter ω appears in a contrived form as the starting point of the true anomaly only in the argument of the Jacobi functions, therefore, it can be eliminated from the two equations by using the Jacobi elliptic functions addition theorem. Let v = F (χ i |m) + ∆λ/n be the argument of elliptic functions at the final point of the orbit and z = F (χ i |m) be the argument of elliptic functions at the initial point of the orbit. Then one can use orbit equations to write:  : For a known initial point P i , the two critical orbits through this point divide the orbital plane into three regions. Left: initial point at r > 3M . If the final point P f is in region 3, then all orbits from P i are of type A. If P f is in region 1, then all orbits from P i are of type B. If P f is in region 2, then orbits leading to P f are of type A if the trajectory winds around the black hole, and of type B if it goes directly from P i to P f . Right: initial point at r < 3M . If the final point P f is in region 3, then all orbits from P i are of type B. If P f is in region 1, then orbits from P i are of type C. If P f is in region 2, then orbits leading to P f are of type C if the trajectory winds around the black hole, and of type B if it goes directly from P i to P f .  Figure 4: Functions cn(x|m) and sn(x|m) along orbits: left type A, middle type B, right type C. x min and x max are at the endpoints of the orbit and are designated by P − and P + if they are at spatial infinity, and P h− and P h+ if they are at the horizon of the black hole. Interval definitions for the endpoints follow after equations (10), (13) and (15). x i = z and x f = v are at the initial and the final point of the orbit.
This is a non-linear equation for a (i.e. for ψ or µ with respect to the type of the orbit). Equations (17a -18f, except 18a and 18d) only give squares, so functions cn and sn are determined only up to the sign; dn is always positive. In order to determine those signs we plot in Figure 4 the functions sn and cn along orbits together with the interval where the orbit is defined. Figure 4 left shows that for type A orbits, only the function cn changes sign at the periastron and the function sn is always positive along the orbit. Thus, the sign of cn is positive if the orbit did not pass the periaston and negative otherwise. For orbits of type B (Figure 4 middle) sn is always positive, while the sign of cn is determined by (18a) and (18d), thus no sign ambiguity arises. For orbits of type C the sign of cn is again unambiguous ((18a) and (18d)) while the function sn changes sign from negative before apastron to positive after apastron (Figure 4 right). We conclude that the right hand side of equation (19) has two branches if it contains a sign ambiguity. The cases are:  Figure 5: Determining the angular momentum parameter for a type A orbit. Left: the orbit from P i to P f does not pass periastron. Right: the orbit from P i to P ′ f passes periastron. The designations right 1 and right 2 refer to the branches (20) and (21) respectively. type A "right 1 " and "right 2 " and type C "right 1 " and "right 3 ", where right 1 = cn(v|m)cn(z|m) + sn(v|m)sn(z|m)dn(v|m)dn(z|m) 1 − msn 2 (v|m)sn 2 (z|m) (20) right 3 = cn(v|m)cn(z|m) − sn(v|m)sn(z|m)dn(v|m)dn(z|m) 1 − msn 2 (v|m)sn 2 (z|m) .
In the above the sign of cn or sn is positive when calculated from a square root. Some examples of determining the value of a are shown in figures 5-9. They show the left hand side and all the appropriate branches of the right hand side of (19) as a function of ψ or µ for different cases of r i , r f and ∆λ f . The solution of equation (19), which is the point where "right" crosses "left", is found numerically by using Brent method (Press et al. 1988). 3 With the now known parameter a it is straightforward to calculate the value of ω (equations 10,13,15).

Calculation of the time of flight
The time of flight for a photon, is given by the integral (Chandrasekhar 1992): .
The indefinite integral in the above formula can be expressed in analytic form using elliptic integrals. After considerable algebraic manipulations we obtain the analytic forms for the three cases of type A, B and C orbits as follows.  Figure 7: Determining the angular momentum parameter for a type B orbit. Left: a direct orbit from P i to P f . Right: a winding orbit from P i to P f (winding number k = 1).  Figure 8: Determining the angular momentum parameter for a type C orbit. Left: the orbit from P i to P f does not pass apastron. Right: the orbit from P i to P ′ f passes apastron. The designations right 1 and right 3 refer to the branches (20) and (22) Figure 9: Connecting P i and P f with a type B orbit and a type C winding orbit (winding number k = 1); no k = 0 type C orbit exists in this case. Left: the orbit from P i to P f of type C . Right: the orbit from P i to P f of type B .
The type C orbit starts at the event horizon at t → −∞, passes the apastron at t = M and again approaches the event horizon as t → ∞. In the above example a = 2 3 √ 3 − 10 −3 for type A and C orbits and a = 2 3 √ 3 + 5 · 10 −4 for type B orbit.
For type B orbits the parameters u 1 , m and n are those from (12a -12c), while for type C orbits these parameters are from (14a -14e).
In figure 10 we show examples of r(t) for all three types of orbits. We compared the effectivness of the analytic method and the pure numerical method to calculate the times of flight. It was found that the algorithm based on the analytic method and implemented with the Carlson's algorithm (Carlson 1979) is always more accurate and 3 to 6 times faster than fourth-order Runge-Kutta integration with adaptive stepsize control. 4 We also note (Kostić 2003) that it is possible to calculate the time of flight numerically by expanding the integrand (23) into piecewise rapidly convergent series of analytically integrable functions. Such series easily give results accurate to 10 −8 M for subcritical orbits and are some 6 times faster than the analytic method. The most effective method to calculate the time of flight would thus be a combination of the analytic method for the close to critical orbits and the series solution for the rest.

Conclusions
In this work we presented the complete solution of the ray-tracing problem in the Schwarzschild space-time. All the algorithms presented here have been tested for accuracy and for speed of execution and were found to be more accurate and considerably faster than commonly used direct integration methods. We hope that algorithms presented here will be used as a useful tool in solving more complex ray tracing problems that will elucidate the physics governing complicated transient phenomena in the vicinity of black holes. We would like to remind the community that a similar ray-tracing problem in the Kerr space-time still remains to be solved.