Table Of ContentGeneral Relativity and Gravitation manuscript No.
(will be inserted by the editor)
Analytical time-like geodesics in Schwarzschild
space-time
Uroˇs Kosti´c
Received:29August2011/Accepted:15January2012
2 Abstract Time-like orbits in Schwarzschild space-time are presented and classi-
1 fied in a very transparent and straightforward way into four types. The analytical
0 solutions to orbit, time, and proper time equations are given for all orbit types in
2
the form r=r(λ), t=t(χ), and τ =τ(χ), where λ is the true anomaly and χ is a
n parameter along the orbit. A very simple relation between λ and χ is also shown.
a These solutions are very useful for modelling temporal evolution of transient phe-
J
nomena near black holes since they are expressed with Jacobi elliptic functions
6 and elliptic integrals, which can be calculated very efficiently and accurately.
2
Keywords Schwarzschild space-time · analytical solutions · time-like geodesics
]
c
q
- 1 Introduction
r
g
[ When modelling physical phenomena occurring in strong gravitational field of
black holes, it is common to work in the Schwarzschild space-time [1,2,3,4,5,6,7,
1
8,9,10] . The same applies also for determination of orbital parameters of objects
v
1 orbiting so close to a black hole that the orbits are affected by relativistic effects,
1 e.g. highly eccentric S stars at the Galactic Centre [11,12]. For this purpose, an
6 efficient and accurate way for calculating time-like orbits in the Schwarzschild
5 space-time is required. Moreover, if we are interested in time-dependence of these
.
1 phenomena, a method for solving the time equation is also required to calculate
0 the temporal evolution and the dynamics.
2 Since in such models, the number of calculations can rapidly increase either
1
because of increasing the number of points, or extending the time of the simu-
:
v lation, or reducing the time-step size, it is very desirable to have a very efficient
i and accurate method for solving these equations. For example, in a model which
X
includes time-dependant gravitational lensing, it is easy to miss the moment of
r
the strongest lensing when the Einstein ring appears, if the time-step is too large.
a
UroˇsKosti´c
FacultyofMathematicsandPhysics,UniversityofLjubljana
Jadranska19,1000Ljubljana,Slovenia
E-mail:[email protected]
2 UroˇsKosti´c
Consequently, the calculated signal, as received by a distant observer, lacks this
distinctive characteristic.
Although the numerical integration or post-Newtonian approximation yield
useful results in specific cases, analytical solutions of the orbit and time equation
are simpler, more efficient regardless of the accuracy required (as shown by Delva
[13]), and can be used in all cases (weak field limit, strong field limit). The well
knownworkofChandrasekhar[14]andRauch[15],wherethesolutionstogeodesic
equations are expressed in terms of elliptic integrals, has been followed by Cˇadeˇz
[16,17] and Gomboc [18] who inverted the expressions of Chandrasekhar [14] and
Rauch[15]intoJacobiellipticfunctions,whichnolongercontainthebranchambi-
guity. For light-like geodesics, Cˇadeˇz and Kosti´c [17] presented a very simple and
straightforwardwayofcharacterizingtheorbitswhichdependonlyononeparam-
eter, as well as giving analytical solutions to the time equation and a method of
determining a light-like geodesic between two arbitrary points (and thus facilitat-
ing ray-tracing used in numerical modelling of dynamical phenomena near black
holes).
Cruzetal.[19]haveclassifiedthelight-likeandtime-likegeodesicsaccordingto
the effective potential and found similar analytical solutions as [17], however they
givesolutionstotimeequationonlyforradialandcircularorbits.HioeandKuebel
[20]presentanalyticalsolutionstoorbitequations,classifythemaccordingtotwo
parameters, and show extensive tables of different values of these parameters for
corresponding orbits. They, however, do not give any solution to time equation.
To complement previous work on light-like orbits [17], the complete analytical
solutions of the time-like geodesics and the time equation for all orbit types are
presentedinthispaper:intheformr=r(λ)(whereλisthetrueanomaly)forthe
radial coordinate r, and in the form t=t(χ) and τ =τ(χ) for time t and proper
time τ, respectively, with a very simple relation between λ and χ.
2 Schwarzschild space-time
InSchwarzshildspace-timeweuseSchwarzshildcoordinatest,r,θ,ϕ.TheHamil-
tonian, from which geodesic equations are derived is
(cid:34) (cid:35)
1 1 (cid:16) 2M(cid:17) 1 (cid:16) 1 (cid:17)
H = − p2+ 1− p2+ p2+ p2 , (1)
2 1− 2M t r r r2 θ sin2θ ϕ
r
where p are canonical momenta and natural units c = G = 1 are used. The
µ
constants of motion are: value of Hamiltonian (H) and Lagrangian (L), energy
E = p , three components of angular momentum (l), longitude of periapsis (ω),
t
andtimeofperiapsispassage(t ).Fortime-likegeodesics,thevalueofHamiltonian
p
is H =−1/2.
In order to describe the position along the orbit, as well as the orientation of
the orbit, we introduce another local inertial (right-handed) orthonormal tetrad
nˆ, eˆ and eˆ as shown in Fig.1. The vector nˆ is a constant unit vector pointing
1 2
in the direction of angular momentum (l=lnˆ). The two unit vectors eˆ and eˆ in
1 2
theorbitalplaneareorientedsothateˆ pointsinthedirectionofinitialperiapsis,
1
apoapsis or toward the infinity (The choice depends on the orbit type and will be
Analyticaltime-likegeodesics 3
Fig.1 Theorbitalplaneinequatorialcoor-
dinates:nˆ unitnormal,ιinclination,Ω lon-
gitudeoftheascendingnode,ωlongitudeof
periapsisandλtrueanomaly.
explained further in the text.). The components of these vectors with respect to
the local Cartesian coordinate basis are expressed as in [21]:
cosωcosΩ−cosιsinωsinΩ
eˆ1 =cosωsinΩ+cosιsinωcosΩ (2a)
sinιsinω
−sinωcosΩ−cosιcosωsinΩ
eˆ2 =−sinωsinΩ+cosιcosωcosΩ (2b)
sinιcosω
nˆ =(sinιsinΩ,−sinιcosΩ,cosι) (2c)
whereΩ isthelongitudeoftheascendingnodeandιistheinclinationoftheorbit
with respect to the X−Y plane (see Fig.1).
By introducing a dimensionless variable
2M
u= (3)
r
and two dimensionless constants of motion related to orbital energy and orbital
angular momentum [18]:
2ME (cid:18)2M(cid:19)2
a= and b=2H (4)
l l
one can derive the differential orbit equation:
du (cid:112)
=± a2−u2(1−u)+b(1−u) , (5)
dλ
where λ is the true anomaly. As functions of u, time and proper time obey the
following differential equations:
dt 2Ma
= (6a)
du u2(1−u)(cid:112)a2−u2(1−u)+b(1−u)
dτ 2Ma 1
= . (6b)
du E u2(cid:112)a2−u2(1−u)+b(1−u)
4 UroˇsKosti´c
D>0 a2 D£0 a2
uL uL
PH a2+b PH
a2+b B 0 A C 0
0 1 0 u u 1
2 1
u u
D>0 a2 D£0 a2
C 0 PuHL PuHL
D C
a2+b a2+b 0
0 u1 1 0u3 u2 u1 1
u u
Fig.2 ThepolynomialP(u)anddistributionofitsrootsintheintervalu∈[0,1].Orbitsexist
onlywhereP(u)≥0,shownincolours.CorrespondingorbittypesaremarkedwithlettersA,
B, C or D and the roots are marked with u1, u2 and u3. The sign of the discriminant D is
alsonoted.
After (5) is solved for u as a function of λ, the orbit equation is written in vector
form as
2M
r(λ)= (eˆ cosλ+eˆ sinλ) . (7)
1 2
u(λ)
Solutions depend on the type of orbit, e.g. closed, scattering or plunging, and
in the following section we present them for all types of time-like geodesics.1
2.1 Types of orbits
Markingthepolynomialin(5)withP(u)=a2−u2(1−u)+b(1−u),thesolutions
to (5)–(6b) exist only on intervals where P(u) ≥ 0. This polynomial has three
roots, while the discriminant D, which is defined as:
27
α=1−9b− a2 (8)
2
β =−1−3b (9)
D=α2+β3 , (10)
determinesthenatureoftheseroots(i.e.thenumberofreal/complexroots).Since
orbits extend at most from u = 0 to u = 1, only roots on this interval are of
interest. In Fig.2, the polynomial P(u) is plotted for all the four possible orbit
types (according to the number of roots in the interval u∈(0,1)).
Theclassificationoforbitsismoreintuitivewhenitisdonewithrespecttothe
effective potential V defined as [22]
(cid:113)
V = (1−u)(1+˜l2u2) , (11)
1 The differential equations (5) and (6a) are formally the same for light-like geodesics [17]
(forlight-likegeodesicstakeb=0).
Analyticaltime-likegeodesics 5
Fig. 3 The effective potential V for
time-like geodesics (Eq.11) for ˜l = 2.2.
B By choosing appropriate value of E, or-
1.05 bits of any type can be constructed: or-
C A bits with E > Vmax (blue area) are
of type B (plunging orbits), orbits with
E 1.00 C D VHŽl=2.2L 1of≤typEe≤A(Vsmcaatxte(ryineglloowrbiatrse)ao)raorfetyeipteheCr
C (nearorbits),orbitswithVmin≤E<1
0.95 (greenarea)areeitheroftypeD(bound
orbits) or of type C, and orbits with
E < Vmin (red area) are only of type
0.90 C. Note that for Vmin ≤E ≤Vmax the
2 3 45 7 1100 20 30 50 100 200 600 discriminantisD≤0,andD>0other-
r@MD wise.
where ˜l=l/2M is the reduced angular momentum. Unlike in Keplerian case, the
effective potential gains a maximum V
max
(cid:32)(cid:32) ˜l2 (cid:33)(cid:32) 1 (cid:33)(cid:33)1/2
V = 1+ 1− (12)
max (cid:112) (cid:112)
(˜l2−˜l ˜l2−3)2 ˜l2−˜l ˜l2−3
at radius r 2
max
(cid:18) (cid:113) (cid:19)
r =2M˜l ˜l− ˜l2−3 . (13)
max
√
Clearly, the maximum exists only for ˜l > 3. The existence of V greatly
max
affects the nature of orbits,3 especially if the orbital energy E is E ∼V : such
max
orbits can wind around the black hole at r several times before continuing
max
eitherawayfromortowardstheblackhole,anddonotexistincaseofNewtonian
√
potential.If˜l= 3,themaximum(andtheminimum)ofthepotentialdisappears
at r =6M, which is the radius of the last stable circular orbit.
max
The effective potential and corresponding orbit types are shown in Fig.3. The
four types of orbits for massive particles have the following properties:
- type A: scattering orbits with both endpoints at infinity. Scattering orbits can
never extend below r=3M.
- type B: plunging orbits with one end at infinity and the other behind the
horizon,
- type C: near orbits with both ends behind the horizon of the black hole.
- type D: bound orbits. Highly eccentric orbits can never reach below r = 4M
while circular orbits can never reach below r=6M.4
SometypicalexamplesofalltypesareshowninFig.4.Notethat,onlyifE ≈V ,
max
then r corresponds to the radius of periapsis for type A and D orbits, and
max
apoapsis for type C orbits.
2 Note that rmax is not the maximal radius an orbit can extend to, but the radius where
Vmax=V(rmax).
3 Obviously,theorbitscanexistonlyforE≥V.
4 Highly eccentric orbits are orbits with energy E (cid:46) 1 (which makes the orbits almost
parabolic).Fromequations(12)and(13)itfollowsthatfortypeDorbits,rmaxisthesmallest
if Vmax ≈E ≈1, which happens for˜l=2 at rmax ≈4M. In this case, rmax corresponds to
theperiapsisdistance.
6 UroˇsKosti´c
C
A B D
Fig. 4 Time-like geodesics with˜l=2.32379 (for types A, B, and C) and˜l=3 (for ype D).
Fromlefttoright:OrbitsoftypeAwithE∈{1.0001,1.035,1.06,1.083},orbitsoftypeB with
E∈{1.0887,1.2,1.6,2.5},orbitsoftypeC withE∈{0.7,0.97372899,1.05,1.086},anorbitof
typeD withE=0.988.TheradiusoftheblackcircleistheSchwarzschildradius.
Fig. 5 The effective potential V for
1.02 time-like geodesics (Eq.11) for ˜l = 1.9.
B
BychoosingappropriatevalueofE,only
1.00
Ž orbits of type B, C, and D can be con-
0.98 C VHl=1.9L structed: orbits with E > 1 (blue area)
are of type B (plunging orbits), orbits
C D
E 0.96 with Vmax ≤ E ≤ 1 (top red area)
are of type C (near orbits), orbits with
C
0.94 Vmin ≤ E < Vmax (green area) are ei-
theroftypeD(boundorbits)oroftype
0.92
C, and orbits with E < Vmin (bottom
0.90 red area) are of type C. Note that for
2 3 45 7 1100 20 30 50 100 200 600 Vmin ≤ E ≤ Vmax the discriminant is
r@MD D≤0,andD>0otherwise.
Fig.6 OrbitoftypeCfor˜l=1.9andE=0.99.Theradius
oftheblackcircleistheSchwarzschildradius.
If, however, ˜l < 2 then V < 1 (see (12)) and consequently, orbits of type
max
A no longer exist, as shown in Fig.5. Moreover, while orbits of type C still have
both endpoints behind the horizon of the black hole, they can extend to infinity
for E → 1. An example of such extended type C orbit for ˜l = 1.9 is in Fig.6.
√
Furthermore,if˜l isloweredbelow 3,alsotypeDorbitsnolongerexist,andonly
orbits of type B and C remain.
Radial and circular orbits can be considered as special cases of type B and D
orbits, respectively. The corresponding equations and parameters for radial orbits
are:t˙=E/(1−u),r˙ =(cid:112)−(1−u)+E2 withzeroangularmomentum˜l=0,while
forcircularorbitstheyare:r˙ =0,λ˙ =l/r2 =const.,t˙=E/(1−u)=const.,with
energyE =V (˜l)whereV (˜l)istheminimumoftheeffectivepotential(11).
min min
2.2 Analytical solutions
The solutions of the equations (5)–(6b) are the following.
Analyticaltime-likegeodesics 7
2.2.1 Types A and D
In this case, the polynomial P(u) has either two (u and u ) or three (u , u ,
1 2 1 2
and u ) real roots on the interval (0,1) which can be elegantly expressed with
3
the constants α, β, and D using Cardan’s formula [23] by introducing two more
intermediary constants D and ψ [18]:
|D|=(cid:112)−β (14)
(cid:32) √ (cid:33)
−D
ψ =2arctan . (15)
(cid:112)
α+ −β3
In terms of these, the roots can be written in the trigonometric form:
1(cid:16) ψ(cid:17)
u = 1+2|D|cos (16a)
1
3 3
1(cid:16) ψ−2π(cid:17)
u = 1+2|D|cos (16b)
2
3 3
1(cid:16) ψ+2π(cid:17)
u = 1+2|D|cos . (16c)
3
3 3
Theserootscanbeassociatedtotheradiusofperiapsisr =2M/u (typesA,D)
p 2
and apoapsis r =2M/u (type D only). Since the argument of arctan in (15) is
a 3
positive, it follows that 0≤ψ ≤π, therefore u >u >u .
1 2 3
Using the substitution [24]
u(χ)=u −(u −u )cos2χ , (17)
2 2 3
equations (5)–(6b) are transformed into Legendre form of elliptic integrals and
integrated to obtain orbital variables λ, t, and τ as functions of χ:
λ(χ)=n(F(χ|m)−K(m)) (18)
2na(cid:34)(cid:18) n2−m (cid:19) u2
t(χ)= 1+u + 1 Π(n ;χ|m)+ 3 Π(n ;χ|m)
u2 3 2(m−n )(n −1) 1 1−u 2
3 1 1 3
n /2 (cid:18) (cid:16) m(cid:17)
+ 1 E(χ|m)− 1− F(χ|m)
(m−n )(n −1) n
1 1 1
(cid:112) (cid:33)(cid:35)
n sin2χ 1−msin2χ
− 1
2(1−n sin2χ)
1
(19)
(cid:18) (cid:19)
1 2n u
τ(χ)= t(χ)− Π(n ;χ|m)+ 3 Π(n ;χ|m) , (20)
E ˜lu3 1 1−u3 2
where:
u −u
m= 2 3 (21a)
u −u
1 3
2
n= √ (21b)
u −u
1 3
u
n =1− 2 (21c)
1
u
3
u −u
n = 2 3 . (21d)
2
1−u
3
8 UroˇsKosti´c
250
5
200
-30 -20 -10
-5 M@D 150
ΤM,D
@
t
100
5
50
-30 -20 -10
-5
0
10 15 20 25 30 35
r@MD
Fig. 7 Left: Orbit of type A with E = 1.01 and ˜l = 2.2. Black and red dots correspond to
pointsattimeintervals∆t=5M andpropertimeintervals∆τ =5M,respectively.Theblack
circle represents the Schwarzschild radius. Right: Time (black) and proper time (red) as a
functionofcoordinater,measuredfromtheinitialpointatr=34M.
Inverting (18) by χ(λ) = am(K(m)+λ/n |m) and using (17) one can also write
the solution to the orbit equation (5) as a function of true anomaly in the form:
u(λ)=u −(u −u )cn2(cid:0)K(m)+ λ|m(cid:1) . (22)
2 2 3
n
FortypeD orbits,bothλandχcangofrom−∞to+∞.FortypeA,thevalues
(cid:112)
of χ are in the interval χ ∈ (χ ,χ ), where χ = arccos( u /(u −u ))
min max min 2 2 3
(cid:112)
and χ = arccos(− u /(u −u )), while the values of λ are in the interval
max 2 2 3
λ/n∈(F(χ |m)−K(m),F(χ |m)−K(m)).Thevaluesofχatperiapsisand
min max
apoapsisareπ/2and0respectively,whileλ=0atperiapsis.Definitionsofelliptic
integrals and functions are from Wolfram [25].
In Figures 7 and 8, an example of solutions for a type A and type D orbits
are shown, with black and red dots marking equal time and proper time intervals
∆t=∆τ =5M. As expected, the dots are more widely spaced when closer to the
black hole, and the lengths of sections corresponding to proper time intervals ∆τ
are longer than those corresponding to time intervals ∆t (which is also clear in
the t=t(r) and τ =τ(r) plots of Fig.7 and 8, where t(r)>τ(r) for all r). Note
thatincaseoftypeDorbit,thedotsareplottedonlyforoneorbitalperiod,while
the orbit is plotted for 3 periods to show the periapsis precession.
In order to compare the efficiency and accuracy of the analytical expression
(19)forttoadirectnumericalintegration,equation(6a)hasbeenintegratedusing
fourth-orderRunge-Kuttamethodwithadaptivestep-sizecontrol[26].Theelliptic
integrals in equation (19) were calculated by Carlson’s algorithm [27], while the
Jacobi elliptic functions in (18) were from [26].
For type A orbit (E = 1.01, ˜l = 2.2) the integration limits were r =
min
6.15313M ≈r andr =50M,whilefortypeDorbits(E =0.9704,˜l=1.888),
p max
the limits were r =5.04581M ≈r and r =25.436M ∼r . In both cases,
min p max a
the numerical integration fails, if r gets too close to either r or r since these are
p a
the zeroes of the polynomial P(u) in (6a). Taking e.g. r = r (1+10−8) and
min p
Analyticaltime-likegeodesics 9
15
10 500
5
400
-5 5 10 15 20 25
-5
-10 MD 300
@
-15 ΤM,D
@
15 t
200
10
5
100
-5 5 10 15 20 25
-5
-10 0
5 10 15 20 25
-15 r@MD
Fig. 8 Left:OrbitoftypeDwithE =0.9704and˜l=1.888.Blackandreddotscorrespond
to points at time intervals ∆t = 5M and proper time intervals ∆τ = 5M, respectively. The
blackcirclerepresentstheSchwarzschildradius.Right:Time(black)andpropertime(red)as
afunctionofcoordinater,measuredfromtheinitialpointatr=rp∼5.045M.
r =r (1−10−8)andthusavoidingthedivergence,5itturnsoutthatnumerical
max a
integrationis∼50−80timesslowerthananalyticalsolution(19).Inaddition,the
relative error δt/t is ∼4 orders of magnitude and ∼2 orders of magnitude larger
for numerical integration than for analytical solution (19) in case of type A and
type D orbits, respectively.
Itshouldbealsonotedthatsomeadditionaleffortisrequiredwhennumerically
integrating(6a):iftheorbitpasseseitherr orr ,e.g.atypeDorbitspansmany
a p
periods (or even just one!), or a type A orbit passes the periapsis, some book-
keeping of periapsis and apoapsis passages has to be done in order to obtain
the correct solution, e.g. by adding the correct number of half-periods. If using
analytical solution, no such additional work is necessary, since (19) is essentially
expressed with an angle along the orbit.
2.2.2 Type B
The polynomial P(u) can be factorized as P(u) = (u−u )(u2+pu+q), where
1
u <0 is the only real root (see Fig.2). The coefficients p, q, and the root u are
1 1
5 Consequently,theperiapsisandapoapsispassagetimeshavetobecalculatedinadifferent
manner.
10 UroˇsKosti´c
expressed as [23,18]:
√
D=(α− D)1/3 (23a)
√
D(cid:101) =(α+ D)1/3 (23b)
1(cid:16) (cid:17)
u1 = 1+D+D(cid:101) (23c)
3
p=u −1 (23d)
1
q =−b+pu . (23e)
1
Using the substitution [24]
(cid:113) χ
u(χ)=u + u2+pu +qtan2 , (24)
1 1 1 2
equations (5)–(6b) are transformed into Legendre form of elliptic integrals and
integrated to obtain orbital variables λ, t, and τ as functions of χ:
λ(χ)=n(F(χ|m)−F(χ |m)) (25)
∞
(cid:40) 1 (cid:34) (cid:32) (cid:18) 1 (cid:19)2(cid:18) n2−m (cid:19) (cid:33) (cid:18) 1 (cid:19)(cid:35)
t(χ)=2a α n2 1+ 2+ 1 −1 −α n 1+ Π(n ;χ|m)
k2 2 1 k n (m−n ) 1 1 k 1
1 1 1 1 1
(cid:20) (cid:18) (cid:18) (cid:19)(cid:18) (cid:19) (cid:19) (cid:21)
(n −1)(1+k ) 1 1−m
+ 1 1 α n 1+ 1+ −2 −α ln|x |
2(cid:112)|n1−m| 2 1 k1 n1−m 1 1
(cid:20) (cid:18) (cid:19) (cid:21)
1 α α
+ 2α (n −1) 1+ − 1 − 3 F(χ|m)
2 1
k k k
1 1 2
α (cid:18)n (cid:18) 1 (cid:19)(cid:19)2(cid:34) n sinχ(cid:112)1−msin2χ (cid:35)
+ 2 1 1+ E(χ|m)− 1 (1+k cosχ)
m−n1 k1 k1 k1 1−n1sin2χ 1
(cid:18) (cid:19)(cid:20) (cid:21)(cid:41)
1 k
+α (1−n ) 1+ Π(n ;χ|m)+ √ 2 ln|x |
3 2 2 2
k2 n2−m
(26)
2α (cid:40)(cid:32)(cid:18) k (cid:19)2(cid:18) m(n −1) (cid:19) (cid:33)
τ(χ)= 2 1 1− 1 −1 Π(n ;χ|m)
˜lk12 k1−1 n1(n1−m) 1
(cid:18) (cid:19)
n (1+k ) k n +1−2m 2k
+ 1 1 1 1 −2 ln|x |+ 1 F(χ|m)
2(cid:112)|n1−m| k1−1 n1−m 1 k1−1
1 (cid:18) k (cid:19)2(cid:34) n sinχ(cid:112)1−msin2χ (cid:35)(cid:41)
+ 1 E(χ|m)− 1 (1+k cosχ) ,
m−n1 k1−1 k1 1−n1sin2χ 1
(27)
