Geodesic Motion in the Schwarzschild Metric

The Geodesic Equations for Massive Particles

The geodesic motion of freely falling particles can be obtained from the solution of the geodesic equations. We choose the coordinates so that the motion is confined in the θ=(π/2) plane. Therefore (OverDot[θ])=0, sin(θ)=1. Then we have to solve the system
(Overscript[r, ¨])=((((OverDot[r])^2)-(e^2))/((r^2)(1-(2/r))))+ ((l^2)/(r^3))(1-(2/r))
(OverDot[t])=(e/(1-(2/r)))
(OverDot[ϕ])=(l/(r^2))
where (OverDot[t])=(d/(dτ))t(τ) etc, and the conserved quantities
e=(1-(2/r))(OverDot[t])
l=(r^2) (OverDot[ϕ])
are given. The radial motion “conservation of energy”  is given by
ℰ= (((e^2)-1)/2)= (1/2) ((OverDot[r])^2)+(Subscript[V, eff])(r)
(Subscript[V, eff])(r) =  (1/2)[(1-(2/r))(1+((l^2)/(r^2)))-1]= -(( 1)/r)+ ((l^2)/(2(r^2)))- ((l^2)/(r^3))
The turning points of the radial motion are calculated from the equation          (type ℰ by [esc]scE[esc])
(Subscript[V, eff])(r)=ℰ

In the above equations, dimensionless quantities are used. To go back to geometrized units, use the dictionary:
 τ-> (τ/M)       t-> (t/M)        r      -> (r/M)             ϕ-> ϕ  
 l-> (l/M)       ℰ-> ℰ        (Subscript[V, eff])-> (Subscript[V, eff])          e-> e

Effective potential for radial motion of massive particle

Plot (Subscript[V, eff])(r) for various angular momenta

l=.;
Veff[r_]:= -(( 1)/r)+((l^2)/(2 (r^2)))- ((l^2)/(r^3));
Plot[Evaluate[{
Veff[r]/.l-> 3.2,
Veff[r]/.l-> 4.0,
Veff[r]/.l-> 4.6,
Veff[r]/.l-> 5.0
}],{r,2,80},PlotRange->{-0.075,0.1},PlotLegends->{"3.2","4.0","4.6","5.0"},ImageSize->Large]

Compute turning points: 
l is l/M

l=4.3; rmax=40.;
Veff[r_]:= -(( 1)/r)+((l^2)/(2 (r^2)))- ((l^2)/(r^3));
Plot[Evaluate[{
Veff[r],
Veff[rmax]
}],{r,2,80},
PlotRange      ->{-0.05,0.05},
PlotLegends ->{"l="<>ToString[l],"E= " <>ToString[Veff[rmax]]},
PlotLabel     -> "Turning Points: "<>ToString[r/.NSolve[Veff[r] ==Veff[rmax],r] ],
ImageSize     ->Large]

Geodesic Equations

Veff[r_]:= -(( 1)/r)+((l^2)/(2 (r^2)))- ((l^2)/(r^3));

solveGeodesicEqs[e_,l_,ϕ0_,r0_,vr0_, τfinal_] := NDSolve[{
t'[τ]    == (e/(1-(2/(r[τ])))),
ϕ'[τ]    == (l/((r[τ])^2))   ,
r''[τ] == -(( (e^2))/(((r[τ])^2)(1-(2/(r[τ])))))+((((r'[τ]))^2)/(((r[τ])^2)(1-(2/(r[τ])))))+(1-(2/(r[τ]))) ((l^2)/((r[τ])^3)),
t[0]==0,ϕ[0]==ϕ0,r[0]==r0,r'[0]==vr0     (* initial conditions               *)
},{t,ϕ,r},{τ,0,τfinal}
];

Plotting function

Skip this section if you are not interested in the plotting commands.
This is the function that plots the results.

eRangeDef={{0,40},{-0.04,0.05}};               (* use this variable for eRange below for most plots *)
plotResults[xyRange_,eRange_] :=                     (* xyRange is the range of the x-y   plot            *)
Module[      {xyrange =xyRange},                          (* eRange  is the range in the ℰ/(Subscript[V, eff])-plot            *)
GraphicsGrid[{                                                           (* Arrange the plots in a grid                       *)
{  (* 1st row *)
Show[{
ParametricPlot[Evaluate[
{r[τ] Cos[ϕ[τ]],r[τ] Sin[ϕ[τ]]}   (* x=r cos(ϕ) y=r sin(ϕ)                             *)
/.sol],{τ,0.1,τmax},AxesLabel->{"x","y"}, PlotRange->xyrange],
Graphics[{Red,Disk[{0,0},2]}]      (* a red disk in the black hole area                 *)
}],
Plot[{
Veff[r],Energy                                           (* The effective potential with the energy level E   *)
},{r,0,40},
PlotRange-> eRange,
PlotLabel-> "Turning Points: "           (* compute the turning points                        *)
<>ToString[Select[r/.NSolve[Veff[r] ==Energy ,r],#∈ Reals && #>0&]] 
<> "     ℰ = "<>ToString[Energy], (*  display the value of E on the plot             *)

AxesLabel->{"r","(Subscript[V, eff])[r]"}]
},(* 1st row *)
{    (* 2nd row *)
Plot[Evaluate[{
r[τ]                                                               (* radial  coordinate as a function of λ            *)
}/.sol],{τ,0,τmax},AxesLabel->{"τ","r[τ]"}, PlotRange->All],
Plot[Evaluate[{
{Mod[ϕ[τ],2π],2π}                            (* angular coordinate as a function of λ            *)
}/.sol],{τ,0,τmax},AxesLabel->{"τ","ϕ[τ]"}, PlotRange->All]
},(* 2nd row *)
{    (* 3rd row *)
Plot[Evaluate[{
t[τ]/τ                                                          (* time  coordinate as a function of λ             *)
}/.sol],{τ,0,τmax},AxesLabel->{"τ","γ[τ]≡((t[τ])/τ)"}, PlotRange->All],
Plot[Evaluate[                                                  (* check numerical errors in ℰ                     *)
( 0.5 ((r'[τ])^2)+Veff[r[τ]]-Energy)/.sol
],{τ,0,τmax},AxesLabel->{"τ","Δℰ"},PlotLabel->"Error in ℰ"]
}   (* 3rd row *)

},ImageSize->Full]
]

Solutions to the geodesic equations for massive particles

Set:
r1:                            one of the turning points
l  :                             the angular momentum
r0:                            initial radial position
τmax:                     maximum time for numerical integration
radialdirection: the sign of the initial radial velocity (OverDot[r])(0)

r1 will determine the energy ℰ  and the conserved quantity e, which in turn will determine |(OverDot[r])(0)|
v0 = (radialdirection)×|(OverDot[r])(0)|

Bound trajectories

r1     = 30.0   ;l=4.3;  τmax=10300 ;
r0     = 29.5   ; radialdirection = -1; (* Set it to +/- 1. The sign of (OverDot[r])(0)  *)

Energy=Veff[r1];e=(Sqrt[1+2 Energy]);
v0= radialdirection(Sqrt[2(Energy-Veff[r0])]);
ϕ0=0;

sol = solveGeodesicEqs[e,l,ϕ0,r0,v0, τmax];gg= plotResults[All,eRangeDef]

r1= 500.0;l=4.3;
τmax=132000;
r0=490.5; radialdirection = -1; (* Set it to +/- 1. The sign of (OverDot[r])(0)  *)

Energy=Veff[r1];e=(Sqrt[1+2 Energy]);
v0= radialdirection(Sqrt[2(Energy-Veff[r0])]);
ϕ0=0;

sol = solveGeodesicEqs[e,l,ϕ0,r0,v0, τmax];gg= plotResults[All,eRangeDef]

Scattering

Set:
b: the impact parameter
x0: initial x coordinate
V0: (Subscript[v, 0]), the x-component  of the velocity


Then:
(Subscript[r, 0])        =(Sqrt[((Subscript[x, 0])^2)+(b^2)])
(Subscript[(OverDot[r]), 0])        =(Subscript[v, 0]) cos((Subscript[ϕ, 0]))              We assume (r≫2): (v^r)=(dr/dt)=(dr/(dτ))((dτ)/dt)≈(dr/(dτ))=(OverDot[r])
(Subscript[r, 0])  (Subscript[(OverDot[ϕ]), 0])=(Subscript[v, 0]) sin((Subscript[ϕ, 0]))                                                                 ((dϕ)/dt)≈((dϕ)/(dτ))
l           = ((Subscript[r, 0])^2) (Subscript[(OverDot[ϕ]), 0])                    = (Subscript[r, 0]) (Subscript[v, 0]) sin((Subscript[ϕ, 0]))
ℰ        = (1/2)((Subscript[(OverDot[r]), 0])^2)+(Subscript[V, eff])((Subscript[r, 0]))

b=15.;x0=500.;V0=0.30;
τmax=1800;
r0 =(Sqrt[(x0^2)+(b^2)]);
ϕ0 = ArcTan[(b/x0)];
v0 =     -V0  Cos[ϕ0];
l   =  r0 V0  Sin[ϕ0];

Energy= 0.5  (v0^2) + Veff[r0];e=(Sqrt[1+2 Energy]);

sol = solveGeodesicEqs[e,l,ϕ0,r0,v0, τmax];gg= plotResults[{{-10,100},{-50,50}},Automatic];

Print[
"Scattering Angle: (degrees)\n"
"θ= (tan^(-1))\!\(\*FractionBox[\(\(\\ \)\(\*SuperscriptBox[\(v\), \(y\)] \((∞)\)\)\), \(\*SuperscriptBox[\(v\), \(x\)] \((∞)\)\)]\) = ",
(180/π)(ArcTan[(* ArcTan[x,y] gives angle in the range (-π,π]*)
r'[τmax] Cos[ϕ[τmax]]-ϕ'[τmax]r[τmax] Sin[ϕ[τmax]]/.sol,r'[τmax] Sin[ϕ[τmax]]+ϕ'[τmax]r[τmax] Cos[ϕ[τmax]]/.sol 
]) 
];
Show[gg]

b=15.;x0=500.;V0=0.30;
τmax=1800;
r0 =(Sqrt[(x0^2)+(b^2)]);
ϕ0 = ArcTan[(b/x0)];
v0 =     -V0  Cos[ϕ0];
l   =  r0 V0  Sin[ϕ0];

Energy= 0.5  (v0^2) + Veff[r0];e=(Sqrt[1+2 Energy]);

sol = solveGeodesicEqs[e,l,ϕ0,r0,v0, τmax];gg= plotResults[{{-10,500},{-10,30}},Automatic];

Print[
"Scattering Angle: (degrees)\n"
"θ= (tan^(-1))\!\(\*FractionBox[\(\(\\ \)\(\*SuperscriptBox[\(v\), \(y\)] \((∞)\)\)\), \(\*SuperscriptBox[\(v\), \(x\)] \((∞)\)\)]\) = ",
(180/π)(ArcTan[(* ArcTan[x,y] gives angle in the range (-π,π]*)
r'[τmax] Cos[ϕ[τmax]]-ϕ'[τmax]r[τmax] Sin[ϕ[τmax]]/.sol,r'[τmax] Sin[ϕ[τmax]]+ϕ'[τmax]r[τmax] Cos[ϕ[τmax]]/.sol 
]) 
];
Show[gg]

Motion starting from r<(Subscript[r, max])

Here a trajectory that cannot escape: ℰ<(Subscript[V, eff])((Subscript[r, max]))
We try to move outwards: (OverDot[r])(0)>0

r1= 5;l=4.3;
τmax=6.12;
r0=3; radialdirection = +1; (* Set it to +/- 1. The sign of (OverDot[r])(0)  *)

Energy=Veff[r1];e=(Sqrt[1+2 Energy]);
v0= radialdirection(Sqrt[2(Energy-Veff[r0])]);
ϕ0=0;

sol = solveGeodesicEqs[e,l,ϕ0,r0,v0, τmax];gg= plotResults[All,eRangeDef]

Now we try to escape: Set ℰ>(Subscript[V, eff])((Subscript[r, max]))

Energy=0.10;l=4.3;
τmax    =100.4;
r0        =3; radialdirection = +1; (* Set it to +/- 1. The sign of (OverDot[r])(0)  *)

e  =(Sqrt[1+2 Energy]);
v0= radialdirection(Sqrt[2(Energy-Veff[r0])]);
ϕ0=0;

sol = solveGeodesicEqs[e,l,ϕ0,r0,v0, τmax];gg= plotResults[All,eRangeDef]

But if we make the mistake to start with (OverDot[r])(0)<0, then we fall into the BH

Energy=0.10;l=4.3;
τmax=1.58;
r0=3; radialdirection = -1; (* Set it to +/- 1. The sign of (OverDot[r])(0)  *)

e=(Sqrt[1+2 Energy]);
v0= radialdirection(Sqrt[2(Energy-Veff[r0])]);
ϕ0=0;

sol = solveGeodesicEqs[e,l,ϕ0,r0,v0, τmax];gg= plotResults[All,eRangeDef]

Precession of the perihelion

r1= 500.0;l=4.3;
τmax=132000;
r0=490.5; radialdirection = -1; (* Set it to +/- 1. The sign of (OverDot[r])(0)  *)

Energy=Veff[r1];e=(Sqrt[1+2 Energy]);
v0= radialdirection(Sqrt[2(Energy-Veff[r0])]);
ϕ0=0;

sol = solveGeodesicEqs[e,l,ϕ0,r0,v0, τmax];gg= plotResults[All,eRangeDef]

Compute the precession angle of the motion:

(*Define ordinary functions for r and ϕ,with capital letters: *)
Φ[τ_]:= ϕ[τ]/.sol[[1,2]];
R[τ_]:= r[τ]/.sol[[1,3]];

(* from the plots above determine approximate τ near maxima and use them below: *)
τ1= τ/.Last[FindMaximum[R[τ],{τ,25000}]];
τ2= τ/.Last[FindMaximum[R[τ],{τ,50000}]];
τ3= τ/.Last[FindMaximum[R[τ],{τ,70000}]];
τ4= τ/.Last[FindMaximum[R[τ],{τ,90000}]];

(* compute the corresponding angles: *)
ϕ1=Φ[τ1];
ϕ2=Φ[τ2];
ϕ3=Φ[τ3];
ϕ4=Φ[τ4];

Print[
"(Subscript[τ, 1])= ",τ1,"  (Subscript[τ, 2])= ",τ2, " (Subscript[τ, 3])= ",τ3, " (Subscript[τ, 4])= ",τ4,"\n",
"(Subscript[ϕ, 1])= ", ϕ1, " = ", Mod[ϕ1,2π] (180/π),"°     ",
"(Subscript[ϕ, 2])= ", ϕ2, " = ", Mod[ϕ2,2π] (180/π), "°     ",
"(Subscript[ϕ, 3])= ",ϕ3, " = ", Mod[ϕ3,2π] (180/π), "°     ",
"(Subscript[ϕ, 4])= ",ϕ4, " = ", Mod[ϕ4,2π] (180/π), "°     \n",
"Define: (Subscript[δϕ, i])= (Subscript[ϕ, i + 1])-(Subscript[ϕ, i])-2π\n",
"(Subscript[δϕ, 1])= ",ϕ2-ϕ1-2π," = ",(ϕ2-ϕ1-2π)(180/π), "°    ",
  "(Subscript[δϕ, 2])= ",ϕ3-ϕ2-2π," = ",(ϕ3-ϕ2-2π)(180/π), "°    ",
  "(Subscript[δϕ, 3])= ",ϕ4-ϕ3-2π," = ",(ϕ4-ϕ3-2π)(180/π), "°    "
]

Acknowledgements

This notebook has been programmed by Konstantinos Anagnostopoulos, Physics Department, National Technical University of Athens, Greece, while he was an instructor of the 4th year undergraduate course "General Relativity and Cosmology". It was created for fun, but it may turn out to be useful to everyone studying the General Theory of Relativity for the first time.
Email: konstant@mail.ntua.gr
Web:    http://physics.ntua.gr/konstant 
It is offered under a GPL/CC BY 4.0 license (in that order, depending on whether they apply on the programming part or the text part of the notebook).