1 Schwarzschild geodesics in terms of elliptic 1 0 functions and the related red shift 2 n a J ∗ Gu¨nter Scharf 1 3 Institut fu¨r Theoretische Physik, Universit¨at Zu¨rich, ] A Winterthurerstr. 190 , CH-8057 Zu¨rich, Switzerland G . h p - o r t s a [ 2 v Abstract 7 0 Using Weierstrassian elliptic functions the exact geodesics in the 2 Schwarzschild metric are expressed in a simple and most transparent 1 form. The results are useful for analyticaland numericalapplications. . 1 Forexamplewecalculatetheperihelionprecessionandthelightdeflec- 0 tion in the post-Einsteinian approximation. The bounded orbits are 1 computedinthepost-Newtonianorder. Asatopicalapplicationwecal- 1 culatethegravitationalredshiftforastarmovingintheSchwarzschild : v field. i X r a ∗e-mail: [email protected] 1 1 Introduction Schwarzschild geodesicsareellipticfunctions,therefore,theyshouldbewrit- ten as such. For this purpose the Weierstrassian elliptic functions are most useful because they lead to simple expressions. The reason for this is that the solution of quartic or cubic equations can be avoided in this way. In a recent paper [1] an analytic solution for the geodesic in the weak- field approximation was given. As pointed out in that paper the progress in the astronomical observations call for better analytical methods. In this respect it is desirable to have the exact geodesics in a form most suited for applications. For the orbits in polar coordinates (next section) this goal can be achieved by using Weierstrass’ -function for which many analytical and P numerical methods are known [2]. Considering the motion in time (section 4) the related ζ- and σ-functions of Weierstrass appear. Jacobian elliptic functions have been used by Darwin [8] for the form of the orbits. After some transformation our result (2.13) agrees with his. But in his second paper he abandons the elliptic functions because they were “not so well adapted to a study of the time in those orbits”. Obviously the Weierstrass functions are better suited for the problem. Indeed, expressing them by Theta functions one gets the natural expansion of the geodesics in powers of the Schwarzschild radius, this expansion involves elementary functions only. The Weierstrass functions have also been used by Hagihara 1 [11] . But he has chosen the variables and constants of integration in a manner which leads to less explicit results. So it is difficult to derive the post-Newtonian corrections to the geodesics given here from his formulas. As a topical application we finally calculate the red shift for a star moving in the Schwarzschild field. The geodesics are also needed for the study of modifications of general relativity ([10], section 5.13). 2 The orbits in polar coordinates r = r(ϕ) 0 1 2 3 We take the coordinates x = ct, x = r, x = ϑ, x = ϕ and write the Schwarzschild metric in the form r r r 2 s 2 2 2 2 2 2 2 ds = − c dt dr r (dϑ +sin ϑdϕ ) (2.1) r − r r − s − 1I am indebted to C. L¨ammerzahl and P. Fiziev for bringing this reference to my attention 2 where r = 2GM is the Schwarzschild radius. We shall assume c = 1 in the s c2 following. The geodesic equation d2xα dxβ dxγ +Γα = 0 (2.2) ds2 βγ ds ds with the Christoffel Γα leads to the following three differential equations βγ 2 d t dt dr ′ +ν = 0 (2.3) ds2 dsds d2r ν′ dt 2 ν′ dr 2 dϕ 2 + e2ν reν = 0 (2.4) ds2 2 ds −2 ds − ds (cid:16) (cid:17) (cid:16) (cid:17) (cid:16) (cid:17) 2 d ϕ 2drdϕ + = 0. (2.5) ds2 rds ds Here we have used the standard representation r 1 s = eν (2.6) − r and have chosen ϑ = π/2 as the plain of motion. The Christoffel symbols can be taken from the Appendix of [3]. Multiplying (2.3) by expν we find ∂ dt eν = 0 ds ds (cid:16) (cid:17) so that dt eν =const. = E ds dt = Ee−ν. (2.7) ds 2 Next multiplying (2.5) by r we get dϕ 2 r = const. = L, ds hence dϕ L = . (2.8) ds r2 For the constants of integration we use the notation of Chandrasekhar [4]. 3 Finally,substituting(2.7)and(2.8)into(2.5)andmultiplyingby2exp( ν) − dr/ds we obtain × d dr 2 L2 e−ν E2e−ν + = 0. (2.9) ds ds − r2 h (cid:16) (cid:17) i Consequently, the square bracket is equal to another constant = b. Then the resulting differential equation can be written as dr 2 L2 = E2+eν b . (2.10) ds − r2 (cid:16) (cid:17) (cid:16) (cid:17) The constant b can be arbitrarily adjusted by rescaling the affine parameter 2 s. Below weshalltakeb = m wheremistherestmassofthetestparticle. − 2 This will enable us to include null geodesics (light rays) with m = 0. Each geodesic is characterized by two constants of the motion: energy E and angular momentum L. Taking the square root of (2.10) and dividing by (2.8) we get dr E2 m2 m2 = − r4+ r r3 r2+r r f(r). (2.11) dϕ s L2 L2 s − s ≡ q Now ϕ = ϕ(r) can be written as an elliptic integral. However, it is better to consider the inverse r = r(ϕ) in terms of elliptic function by usinga formula of Weierstrass ([5], p.452). Let the quartic f(r) be written as 4 3 2 f(r)= a0r +4a1r +6a2r +4a3r+a4, (2.12) and let r1 be a zero f(r1) =0, then a solution of (2.11) is given by ′ f (r1) r = r1+ 4 (ϕ;g2,g3) f′′(r1)/6. (2.13) P − Here (ϕ;g2,g3) is Weierstrass’ -function with invariants P P 2 g2 = a0a4 4a1a3+3a2 (2.14) − 3 2 2 g3 = a0a2a4+2a1a2a3 a2 a0a3 a1a4. (2.15) − − − In our case we have a4 = 0. For the convenience of the reader we reproduce the short proof in the Appendix. The result (2.13) is not yet the solution of our problem because it con- tains too many constants: the invariants g2,g3 and the derivatives of f 4 depend on E,L, but in addition the zero r1 appears. Of course one could calculate r1 as a function of E,L by solving the quartic equation f(r) = 0, but this gives complicated expressions. It is much better to use r1 and a second zero r2 as constants of integration instead of E,L. This is even de- sirable from the astronomers point of view because the zeros of derivative (2.11) are turning points of the geodesic, for example in case of a bounded orbit they can be identified with the perihelion and aphelion of the orbit. In order to express E,L by r1,r2 we write our quartic in the form f(r)= a0r(r r1)(r r2)(r r3) (2.16) − − − 3 2 and compare the coefficients of r ,r ,r with (2.12). This leads to 2 m 4a1 =−a0(r1+r2+r3)= L2rs 6a2 = a0(r1r2+r1r3+r2r3) = 1 (2.17) − 4a3 = a0r1r2r3 = rs. − Since 2 2 E m a0 = L−2 (2.18) we can solve for 2 m r1+r2+r3 r = , (2.19) L2 s r1r2+r1r3+r2r3 2 m r 2 2 s E m = . (2.20) − −r1+r2+r3 In addition we obtain the third zero r1r2 r3 = rs . (2.21) r1r2 r1rs r2rs − − Therelations (2.19-21) allow toexpresseverythingintermsofr1,r2. For the invariants we find 2 1 m 2 1 rs r1+r2+r3 g2 = 12 − 4L2rs = 12 − 4 r1r2+r1r3+r2r3 (2.22) 2 1 m 2 a0 2 g3 = 63 − 48L2rs − 16rs 2 = 1 1 (r1+r2+r3)rs + 1 rs . (2.23) 63 − 48r1r2+r1r3+r2r3 16r1r2+r1r3+r2r3 5 ′ ′′ Here r3 has to be substituted by (2.21). For the derivatives f (r1), f (r1) which appear in our solution (2.13) we obtain ′ r1(r1 r2)(r1 r3) f (r1)= − − (2.24) − r1r2+r1r3+r2r3 ′′ (r1 r2)r1+(r1 r3)r1+(r1 r2)(r1 r3) f (r1)= 2 − − − − . (2.25) − r1r2+r1r3+r2r3 With these substitutions the result (2.13) gives all possible geodesics in the form r = r(ϕ;r1,r2). This will be discussed in the next section. As a first check of the solution (2.13) we consider the Newtonian limit. Let the two zeros r1, r2 be real and large compared to the Schwarzschild radiusr inabsolutevalue. ThenneglectingO(r )intheinvariants(2.22-23) s s the -function becomes elementary ([2],p.652, equation 18.12.27): P 1 1 1 1 1 −3 (ϕ; ,6 )= + = + . (2.26) 2 P 12 −12 4sin ϕ/2 −12 2(1 cosϕ) − The leading order in the derivatives of f is given by ′ r1 f (r1) = (r1 r2) −r2 − ′′ r1 f (r1) = 2 3 2 . (2.27) − r2 − (cid:16) (cid:17) It is convenient to introduce the eccentricity ε by r1 1 ε = − . (2.28) r2 1+ε Using all this in (2.13) we find the wellknown conic (1+ε)r1 r = . (2.29) 1+εcosϕ Assuming both zeros r1,r2 positive and r1 < r2 we have ε < 1 and the orbit is an ellipse with perihelion r1 and aphelion r2. In the hyperbolic case ε > 1 we see from (2.28) that if r1 is positive r2 must be negative. Then there is only one physical turning point r1 which is the point of closest approach. Thelatteralwayscorrespondstoϕ= 0. Therelativisticcorrectionsto(2.29) are calculated in the following section. 6 3 Discussion of the solution Thesolution r = r(ϕ) (2.13) is an elliptic function of ϕ which implies that it is doubly-periodic ([2], p.629 or any book on elliptic functions). The values ′ of the two half-periods ω,ω depend on the three roots of the fundamental cubic equation 3 4e g2e g3 = 0. (3.1) − − Again it is not necessary to solve this equation because the solutions e ,j = j 1,2,3 can be easily obtained from the roots 0,r1,r2,r3 of our quartic f(r)= 0. To see this we transform f(r) to Weierstrass’ normal form as follows. First we set r = 1/x so that from (2.16) we get 1 3 2 f(r)= x4(4a3x +6a2x +4a1x+a0) Next we remove the quadratic term by introducing 1 1 a2 = x = e . (3.2) r a3 − 2 (cid:16) (cid:17) This gives the normal form of Weierstrass 2 a 3 3 f(r)= (e 1a2)4(4e −g2e−g3), (3.3) − 2 with the above invariants (2.14-15). That means roots of f(r) are simply related to roots of (3.1) by the transformation a3 a2 rs 1 e = + = . (3.4) j r 2 4r − 12 j j The cubic equation (3.1) with real coefficients has either three real roots or one real and two complex conjugated roots. The first case occurs if the discriminant 3 2 = g 27g (3.5) 2 3 △ − is positive, in the second case is negative. In terms of theroots is given △ △ by ([2], p.629, equation 18.1.8) 2 2 2 =16(e1 e2) (e2 e3) (e3 e1) . (3.6) △ − − − The physically interesting orbits correspond to the first case of real roots. ∗ If we have two complex conjugated zeros r2 = r1 then (2.28) implies that the eccentricity ε is imaginary. Such orbits have been discussed by Chan- drasekhar ([4], p.111). Now we discuss the various cases. 7 3.1 Bound orbits In this case we have two positive turning points r2 > r1 > 0, consequently there are three real roots e1 > 0> e2 > e3 given by rs 1 1 rsr1+r2 e1 = = (3.7) 4r3 − 12 6 − 4 r1r2 1 r 1 r s s e2 = + , e3 = + , (3.8) −12 4r1 −12 4r2 Our convention is chosen in agreement with [2]. The real half-period ω of the -function is given by ([2], p.549, equation 18.9.8) P ∞ 2 dt K(k ) ω = = (3.9) 4t3 g2t g3 √e1 e3 eZ1 − − − p 2 whereK(k )is thecomplete elliptic integral ofthefirstkindwithparameter 2 e2 e3 r2 r1 k = − = r − s e1 e3 r1r2 rs(2r1+r2) − − r 2ε 3 ε −1 s = 1 r − . (3.10) s r11+ε − 1+ε (cid:16) (cid:17) As a first application let us give the post-Einsteinian correction to the 2 orbital precession. If k (3.10) is small we can use the expansion ([2], p.591, equation 17.3.11) π 1 2 1 3 2 2 2 4 K(k ) = 1+ k + · k +... (3.11) 2 2 2 4 h (cid:16) (cid:17) (cid:16) · (cid:17) i From the roots e1,e3 we find 1 rs2r1+r2 3 2 2r1+r2 2 3 = 2 1+ + r +O(r ) . √e1 e3 2 r1r2 8 s r1r2 s − h (cid:16) (cid:17) i This finally leads to the half-period 2 ω = π 1+ 34rsr1r1+r2r2 + 8(r31rrs2)2 (2r1+r2)(2r2 +r1)+ 83(r2−r1)2 n h io 2 2 3r 1 3 r 18+ε = π 1+ s + s +O(r3) (3.12) 2r11+ε 18r12(1+ε)2 s n o 8 The perihelion precession is given by ϕ = 2(ω π). Then the order r s △ 2 − in (3.12) is Einstein’s result and the O(r ) gives the correction to it. The s accurate computation of the half-period is necessary to control the orbit in the large. To compute the relativistic corrections for r(ϕ) from (2.13) we express the -function by Theta functions ([5], p.464) P ϑ1(z,q) = 2q1/4(sinz q2sin3z+q6sin5z ...) − − ϑ2(z,q) = 2q1/4(cosz+q2cos3z+q6cos5z+...) 3 8 ϑ3(z,q) = 1+2q(cos2z+q cos4z+q cos6z+...) (3.13) 3 8 ϑ4(z,q) = 1 2q(cos2z q cos4z+q cos6z ...). − − − Here q is the so-called Nome ([2], eq. 17.3.21) k2 k2 2 q = +8 +... (3.14) 16 16 (cid:16) (cid:17) 2 These series are rapidly converging since k is small (3.10), they give the natural expansion in powers of the Schwarzschild radius r . Now the - s P function is given in terms of Theta functions by ([2], eq. 18.10.5) π2 ϑ′1(0)ϑ3(φ) 2 P(ϕ) = e2+ 4ω2 ϑ3(0)ϑ1(φ) (3.15) (cid:16) (cid:17) 2 1 r π 1 s 2 2 = + + 1+4q(cos φ 1)+O(q ) , −12 4r1 4ω2sin2φ − h i where π φ= ϕ. (3.16) 2ω Using 2 ε r 3ε+ε ′ s f (r1) =2r1 1+ε − r1(1+ε)2 (cid:16) (cid:17) 2 1 5ε r 1 4ε ε ′′ s f (r1) = −2 1−+ε +6r1 (−1+ε−)2 (3.17) this leads to ′′ 2 f (r1) 1+ε 2εsin φ rs 1 ε (ϕ) = − 1+ 3 (1 cosφ)+ 2 P − 24 4(1+ε)sin φ r11+εcos2φ − − 2 − h (cid:16) 9 3+ε r 2 2 s +2ε sin φ +O . (3.18) 1+ε r1 (cid:17)i (cid:16) (cid:17) Substituting this into (2.13) gives the desired orbit to O(r ) s 2 r(ϕ) 1+ε r 2sin φ 1 ε s = + ε 3+ (1 cosφ) r1 1+εcos2φ r1 1+εcos2φ 1+εcos2φ 2 − − n h 3+ε 3+ε 2 2ε sin φ . (3.19) − 1+ε −1+ε i o It is important to insertthe periodω in φ (3.16) according to (3.12) in order to describe the perihelion precession correctly. ′ If the two roots r1 = r2 coincide, it follows from (2.24) that f (r1) = 0. According to (2.13) we then have circular motion r = r1. If all three zeros coincide r1 = r2 = r3 then (2.21)gives r3 = 3rs which is the innermost circular orbit. 3.2 Unbound orbits Inthis casethereis onlyonephysicalpoint, thepointofclosest approachr1. The other root r2 is negative, therefore, it is better to use the eccentricity ε (2.27) as the second basic quantity. With r3 given by (2.21) we then have 1+ε r1 > r3 > 0 > r2 = r1, (3.20) 1 ε − becauseε 1. Theperiodicityof (2.13) inϕisnow realized byajumptoan ≥ unphysical branch with r < 0. In reality a comet moves on one branch only, but it is a tricky problem to decide on which one. This is due to the fact thattheperioddiffersalittlefrom2π asintheboundedcase. Consequently, neighboring physical branches r > 0 are a little rotated against each other and the distinction between them is not easy. The quantity of physical interest is the direction ϕ∞ of the asymptote. It follows from the original equation (2.11) by integrating the inverse over r from r1 to ∞ ∞ dr ϕ∞ = . (3.21) f(r) rZ1 p This is an elliptic integral which can be transformed to Legendre’s normal form Φ2 µ dΦ ϕ∞ = (3.22) √a0 Z0 1 k2sin2φ − q 10