Geodesic Motion in the Schwarzschild Metric

The Geodesic Equations for Massless 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))))+ (1/(r^3))(1-(2/r))
(OverDot[t])=(e/(1-(2/r)))
(OverDot[ϕ])=(1/(r^2))
where (OverDot[t])=(d/(dλ))t(λ) etc, and the conserved quantities
e=(1-(2/r))(OverDot[t])
1=(r^2) (OverDot[ϕ])
are given. 

We have rescaled the affine parameter of the null geodesic, so that the angular momentum l=1.

The radial motion “conservation of energy”  is given by
(e^2)= (1/2) ((OverDot[r])^2)+(Subscript[W, eff])(r)
(Subscript[W, eff])(r) =  (1/(r^2))(1-(2/r))

The turning points of the radial motion are calculated from the equation    
(Subscript[W, eff])(r)=(e^2)

The maximum of (Subscript[W, eff])(r) is at (Subscript[r, max])=3 and (Subscript[W, eff])((Subscript[r, max]))= (1/27). Therefore, we have circular orbits when (e^2)=(1/27), b=(Sqrt[27])

The impact parameter is
(b^2)=(1/(e^2))

In the above equations, dimensionless quantities are used. To go back to geometrized units, use the dictionary:
First restore l:
λ-> l λ     (d/(dλ))-> (1/l) (d/(dλ))
(OverDot[t])-> (l^(-1)) (OverDot[t])  => e -> (l^(-1)) e    
Then restore M:
 τ-> (τ/M)       t-> (t/M)        r      -> (r/M)             ϕ-> ϕ  
 l-> (l/M)        (Subscript[W, eff])-> (M^2)(Subscript[W, eff])                       e-> e             b-> b

Effective potential for radial motion of massless particle

l=4.3; rmax=8.;
Weff[r_]:= (1/(r^2))(1-(2/r));
Plot[Evaluate[{
Weff[r],
Weff[rmax]
}],{r,2,40},
PlotRange      ->{-0.005,0.05},
PlotLegends ->{"l="<>ToString[l],"(e^2)= " <>ToString[Weff[rmax]]},
PlotLabel     -> "Turning Points: "<>
ToString[Select[r/.NSolve[Weff[r] ==Weff[rmax] ,r],#∈ Reals && #>0&]]<>
"  Impact parameter: b= (1/e) = "<>ToString[(1/(Sqrt[Weff[rmax]]))],
ImageSize     ->Large]

Geodesic Equations

Weff[r_]:= (1/(r^2))(1-(2/r));

solveGeodesicEqs[e_,ϕ0_,r0_,vr0_, λfinal_] := NDSolve[{
t'[τ]    == (e/(1-(2/(r[τ])))),
ϕ'[τ]    == (1/((r[τ])^2))   ,
r''[τ] == -(( (e^2))/(((r[τ])^2)(1-(2/(r[τ])))))+((((r'[τ]))^2)/(((r[τ])^2)(1-(2/(r[τ])))))+(1-(2/(r[τ]))) (1/((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.

plotResults[xyRange_] :=                                       (* xyRange is the range of the x-y plot *)
Module[{xyrange=xyRange},                                          
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[{
Weff[r],Energy                                                              (* The effective potential with the energy level E *)
},{r,0,40},
PlotRange->{-0.001,0.05},
PlotLabel-> "Turning Points: "                               (* compute the turning points                      *)
<>ToString[Select[r/.NSolve[Weff[r] ==Energy ,r],#∈ Reals && #>0&]] 
<> "     (e^2)= "<>ToString[Energy]                  (*  display the value of E on the plot              *)
<>"     b = " <>ToString[(1/e)],
AxesLabel->{"r","(Subscript[W, 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π]                                                              (* 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[λ]"}],
Plot[Evaluate[                                                                     (* check numerical errors in ℰ                     *)
( ((r'[τ])^2)+Weff[r[τ]]-Energy)/.sol
],{τ,0,λmax},AxesLabel->{"λ","Δℰ"},PlotLabel->"Error in ℰ ≡ (e^2)"]
}   (* 3rd row *)

},ImageSize->Full]
]

Solutions to the geodesic equations for massless particles

Set:
r1:                            one of the turning points
e  :                            (b^2)=1/(e^2)
r0:                            initial radial position
λmax:                     maximum time for numerical integration for the affine parameter λ
radialdirection: the sign of the initial radial velocity (OverDot[r])(0)

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

Circular trajectories

The unstable circular orbit for r(0)=(Subscript[r, max])=3, (e^2)=(Subscript[W, eff])((Subscript[r, max]))=(1/(Sqrt[27]))
The numerical integration has an instability for τ≈159

r1     = 3.  ; λmax=150 ;
r0     = 3.  ; radialdirection = -1; (* Set it to +/- 1. The sign of (OverDot[r])(0)  *)

Energy=Weff[r1];e=(Sqrt[ Energy]); b= 1/e;
v0= radialdirection(Sqrt[(Energy-Weff[r0])]);
ϕ0=0;

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

Instability of the circular orbits:
Slightly perturbed away from the (Subscript[r, max])=3, (Subscript[W, max])=(1/27) circular orbit

r1     = 3.0001   ; λmax=150 ;
r0     = 3.0001   ; radialdirection = -1; (* Set it to +/- 1. The sign of (OverDot[r])(0)  *)

Energy=Weff[r1];e=(Sqrt[ Energy]); b= 1/e;
v0= radialdirection(Sqrt[(Energy-Weff[r0])]);
ϕ0=0;

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

Scattering - Deflection of light

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[v, 0])        = 1 = (dx/dt)                  (far from r=0)
(v^r)        =-(Subscript[v, 0]) cos((Subscript[ϕ, 0])) = (dr/dt)= (dr/(dλ))/(dt/(dλ))=(OverDot[r])(0)/(((e(1-(2/r)))^(-1)))=> (OverDot[r])(0)=-((e(1-(2/r)))^(-1))  cos((Subscript[ϕ, 0])) => (OverDot[r])(0)=-(1/b)(((1-(2/r)))^(-1))  cos((Subscript[ϕ, 0])) 
(e^2)        = ((Subscript[(OverDot[r]), 0])^2)+(Subscript[W, eff])((Subscript[r, 0]))

To have finite deflection angle, we must have b>(Sqrt[27])≈5.196

b   = 8; x0 = 500.; λmax=12000 ;

r0 = (Sqrt[(x0^2)+(b^2)]);
ϕ0 = ArcTan[(b/x0)];
v0  = -(( 1)/b)(((1-(2/r0)))^(-1))Cos[ϕ0];

Energy = (v0^2) + Weff[r0];e=(Sqrt[ Energy]); b= 1/e; (* b must be redefined, it is only approximately equal to the impact parameter *)

sol=solveGeodesicEqs[e,ϕ0,r0,v0, λmax];gg=plotResults[{{-40,100},{-40,10}}];

Print[
"Scattering Angle: (degrees)\n",
"θ= (tan^(-1))\!\(\*FractionBox[\(\(\\ \)\(\*SuperscriptBox[\(v\), \(y\)] \((∞)\)\)\), \(\*SuperscriptBox[\(v\), \(x\)] \((∞)\)\)]\) = ",
Δθ=360+(180/π)(ArcTan[r'[λmax] Cos[ϕ[λmax]]-ϕ'[λmax]r[λmax] Sin[ϕ[λmax]]/.sol,r'[λmax] Sin[ϕ[λmax]]+ϕ'[λmax]r[λmax] Cos[ϕ[λmax]]/.sol ]) ,"\n",
"(Subscript[δθ, deflection])= ",Δθ -180
];
Show[gg]

b   = 5.3; x0 = 500.; λmax=12000 ;

r0 = (Sqrt[(x0^2)+(b^2)]);
ϕ0 = ArcTan[(b/x0)];
v0  = -(( 1)/b)(((1-(2/r0)))^(-1))Cos[ϕ0];
Energy = (v0^2) + Weff[r0];e=(Sqrt[ Energy]); b= 1/e;

sol=solveGeodesicEqs[e,ϕ0,r0,v0, λmax];gg  =plotResults[{{-10,100},{-10,70}}]

b   = 5.1; x0 = 500.; λmax=2567 ;

r0 = (Sqrt[(x0^2)+(b^2)]);
ϕ0 = ArcTan[(b/x0)];
v0  = -(( 1)/b)(((1-(2/r0)))^(-1))Cos[ϕ0];

Weff[r_]:= (1/(r^2))(1-(2/r));
Energy = (v0^2) + Weff[r0];e=(Sqrt[ Energy]); b= 1/e; (* b must be redefined, it is only approximately equal to the impact parameter *)

sol=solveGeodesicEqs[e,ϕ0,r0,v0, λmax];gg  =plotResults[{{-5,30},{-10,10}}]

b   = 5.2170155293; x0 = 500.; λmax=3000 ;

r0 = (Sqrt[(x0^2)+(b^2)]);
ϕ0 = ArcTan[(b/x0)];
v0  = -(( 1)/b)(((1-(2/r0)))^(-1))Cos[ϕ0];

Energy = (v0^2) + Weff[r0];e=(Sqrt[ Energy]); b= 1/e; (* b must be redefined, it is only approximately equal to the impact parameter *)

sol=solveGeodesicEqs[e,ϕ0,r0,v0, λmax];gg  =plotResults[{{-5,30},{-10,10}}];

Print[
"Scattering Angle: (degrees)\n",
"θ= (tan^(-1))\!\(\*FractionBox[\(\(\\ \)\(\*SuperscriptBox[\(v\), \(y\)] \((∞)\)\)\), \(\*SuperscriptBox[\(v\), \(x\)] \((∞)\)\)]\) = ",
Δθ=360+(180/π)(ArcTan[r'[λmax] Cos[ϕ[λmax]]-ϕ'[λmax]r[λmax] Sin[ϕ[λmax]]/.sol,r'[λmax] Sin[ϕ[λmax]]+ϕ'[λmax]r[λmax] Cos[ϕ[λmax]]/.sol ]) ,"\n",
"(Subscript[δθ, deflection])= ",Δθ -180
];
Show[gg]

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