Gravity and Geodesics

We work out some consequences of the geodesic equations for gravity around a black hole using our understanding of the geodesic equations on surfaces of revolution. 

It turns out that a slice of the Schwartzchild metric around a black hole (at constant time) has the same geometry as Flamm’s paraboloid, which is the graph of the function 



over R^2. Here r_s is the “Schwartzchild radius” or “event horizon” of the black hole. We can plot this to take a look:

Manipulate[RevolutionPlot3D[2 Sqrt[Rs(R - Rs)],{R,Rs,3Rs}],{Rs,1,5}]

We can plot a geodesic on this surface which shares some features with the trajectory of real massless particles orbiting the black hole (we have to be careful-- since this is a slice of the metric at a single instant in time, only a particle with infinite velocity would follow such a track).

First, we need to express the paraboloid (let’s say rS = 1) as a surface of revolution. We have the curves

ϕ[v_] := v; ψ[v_] := 2 Sqrt[1(v - 1)];

From this, we can write down the geodesic equations from Clairaut’s relation. They are 



We can go ahead and compute the coefficients:

2 D[ϕ[v],v]/ϕ[v]

Simplify[D[ϕ[v],v]ϕ[v]/((D[ϕ[v],v]^2) + (D[ψ[v],v]^2))]

Simplify[(D[ϕ[v],v]D[ϕ[v],v,v] + D[ψ[v],v] D[ψ[v],v,v])/((D[ϕ[v],v]^2) + (D[ψ[v],v]^2))]

This gives some a system of nice ordinary differential equations. We’ll introduce some auxiliary variables uP[s] and vP[s] to represent the derivatives of u[s] and v[s]. Then we have the system

GeodesicSystem = {u'[s] == uP[s], v'[s] == vP[s], uP'[s] == -(2/v[s]) uP[s] vP[s], vP'[s] == (v[s] -1) uP[s]^2 - (1/(2v[s] - 2v[s]^2)) vP[s]^2,u[0] == 0, v[0] == 2, uP[0] == A, vP[0] == B};

GeodesicSystem /. {A->1,B->0}

sol = NDSolve[GeodesicSystem /. {A->0.1,B->0.1}, {u,v,uP,vP},{s,0,5}]

Geodesic2D[S_,sol_] := Module[{s},(Evaluate[{ϕ[v[s]] Cos[u[s]],ϕ[v[s]]Sin[u[s]]} /. sol] /. s-> S)[[1]]];

Geodesic3D[S_,sol_] := Module[{s},(Evaluate[{ϕ[v[s]] Cos[u[s]],ϕ[v[s]]Sin[u[s]],ψ[v[s]]} /. sol] /. s-> S)[[1]]];

DynamicModule[{sol},Manipulate[sol = NDSolve[GeodesicSystem /. {A->Cos[θ],B->Sin[θ]}, {u,v,uP,vP},{s,0,S}]; Show[{Graphics[Circle[{0,0},1]],ParametricPlot[Geodesic2D[s,sol],{s,0.1,S},AxesOrigin->{0,0},AspectRatio->Automatic]},PlotRange->{{-5,5},{-5,5}},GridLines->{Range[-5,5],Range[-5,5]},GridLinesStyle->Directive[Blue,Opacity[0.2]]],
{S,0.2,10},{θ,0,2Pi}]]

This looks even better in 3D:

DynamicModule[{sol},Manipulate[sol = NDSolve[GeodesicSystem /. {A->Cos[θ],B->Sin[θ]}, {u,v,uP,vP},{s,0,S}]; Show[{RevolutionPlot3D[2 Sqrt[1(R - 1)],{R,1,5},PlotStyle->{Opacity[0.2]}],ParametricPlot3D[Geodesic3D[s,sol],{s,0.1,S},AxesOrigin->{0,0},AspectRatio->Automatic,PlotStyle->{Directive[Thick]}]},PlotRange->{{-5,5},{-5,5},{0,4}}],
{S,0.2,10},{θ,0,2Pi}]]

Here’s an interesting example: suppose that geodesics (say, light rays), come from a point source (say, a star) in all directions. We see

Sols = Table[NDSolve[GeodesicSystem /. {A->Cos[θ],B->Sin[θ]}, {u,v,uP,vP},{s,0,10}],{θ,0,2Pi,0.1*Pi}];
Plots = ParametricPlot[Geodesic2D[s,#],{s,0,5},AxesOrigin->{0,0},AspectRatio->Automatic] & /@ Sols;
 ThreeDPlots = ParametricPlot3D[Geodesic3D[s,#],{s,0,5},AxesOrigin->{0,0},AspectRatio->Automatic,PlotStyle->{Thickness[0.0025]}] & /@ Sols;

Show[Graphics[Circle[{0,0},1]],Plots,PlotRange->{{-5,5},{-5,5}},GridLines->{Range[-5,5],Range[-5,5]},GridLinesStyle->Directive[Blue,Opacity[0.2]]]

Show[{RevolutionPlot3D[2 Sqrt[1(R - 1)],{R,1,5},PlotStyle->{Opacity[0.2]}],ThreeDPlots}]

This is something that we actually see in the sky! This is a galaxy whose image is bent by a foreground galaxy, captured by the Hubble telescope.

Thank you for taking this course!