Schwarzschild Orbits

Introduction

The geodesic equations and constants of the motion are derived in the "Schwarzschild Metric" notebook using MathTensor. The results are given in this notebook without requiring MathTensor.   The first order derivatives  are derived from the metric and the constants of the motion. The behavior of the first derivities are explored.  The geodesic equations are solved using NDSolve for various initial conditions.  Solutions and plots of several orbits are explored. For the most part this notebook  follows  "Exploring Black Holes" (EBH)  by Taylor and Wheeler and "Gravity" by Hartle. 
William C. Collins, Space Systems,USNRL(retired), bccollins1@comcast.net

\(\(\(\(Off[General::"spell"]\) \)\(;\)\(\(Off[
        General::"spell1"]\) \)\(;\)\(\(<< 
        Graphics`Graphics`\) \)\(;\)\(\(Off[General::"spell", 
        Replace::"rep"]\) \)\(;\)\(
\)\)\)

Geodesic equations for the equatorial plane

From the Schwarzschild Metric Notebook:θ= (π/2),θ'[s] = 0

SetAttributes[M, Constant]

geor[M_]=((M ((r'[τ])^2))/(2 M r[τ] - r[τ]^2))+((M (-2 M + r[τ]) ((t'[τ])^2))/(r[τ]^3))+(2 M - r[τ]) ((ϕ'[τ])^2)+r''[τ];

geoϕ=((2 r'[τ] ϕ'[τ])/(r[τ]))+ϕ''[τ];

geot[M_]=((2 M r'[τ] t'[τ])/(-2 M r[τ] + r[τ]^2))+t''[τ];

First derivatives with Em and Lm as arguments

Constants of the motion: Em & Lm

Defining β(r,M)

β[r_, M_] = 1 - ((2 M)/r);

Timelike metric:d(τ^2) =β(r, M)(dt^2)-((β(r, M))^(-1)) (dr^2)- ((r )^2)d(ϕ^2)

Energy per unit mass :  Em  =  β (r, 
          M) (dt/d)τ

Angular momentum per unit mass :   Lm  =  
      (r^2) d(ϕ/d)τ

These constants of the motion allow us to determine dt/dτ  and dϕ/dτ  while dr/dτ is obtained  from the metric in terms of these constants.

dt/dτ-> dtdτ[r, Em, M]

dtdτ[r_, Em_, M_] := Em /β[r, M]

dtdτ[r, Em, M]

We can get rid of the explicit M dependancy by expressing r in terms of M

dtdτ[r*M, Em, M]

dϕ/dt ->   dϕdτ[r, Lm]

dϕdτ[r_, Lm_] := (Lm/(r^2))

dϕdτ[r, Lm]

dr/dτ ->   drdτ[r, Em, Lm, M]

drdτ[r_, Em_, Lm_, 
        M_] := (Sqrt[(Em^2) - β[r, M] (1 + ((Lm^2)/(r^2)))])

drdτ[r, Em, Lm, M]

drdτ[r, Em, Lm, M]^2 // Expand

We can get rid of the explicit M dependancy by expressing r and Lm in terms of M

drdτ[r* M, Em, Lm* M, M]^2 // Expand

Effective potential: Vm

Relativistic energy  for flat space

\(\(The\ relativistic\ \ 4\ momentum\ has\ E\ as\ its\ time\ component\ and\ 3\ momentum\ as\ its\ space\ components\ ;\)\)
\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ p\^α\  = \ \(\((E, p\&_)\)\  = \ \((mγ, \ 
          mγ 
             v\&_)\)\), 
The\ magnitude\ is\ the\ mass, \ 
    an\ invariant . 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ p\^α\ p\_α\  = \ \(\(-m\^2\)\  = \ \(\(-E\^2\)\  + \ 
            p\&_\^2\  = \(-\ 
                m\^2\) γ\^2 + \ \ \(m\^2\) \(v\&_\^2\) γ\^2\)\)\)
\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ E\^2\  = \ 
      m\^2\  + \ \ p\&_\^\(\(2\)\(\ \)\)\)
Energy per unit mass
\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Em\^2\ \  = \ \ \ \(E\^2/m\^2\  = \ 
        1\  + \ \ v\&_\^2\)\ \)
\(To\ relate\ this\ to\ Newtonain\ mechanics\ we\ follow\ Hartle\ and\ define\ ξ\)
\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ξ = \ \ \ \(\(Em\^2\  - 1\)\/2 = \ \(1\/2\) v\&_\^2\)\)
\(With\ \ \((dr/dτ)\)\^2/2\  = \ \(\((\ Em\^2 - 1\ )\)/2\  - 
          Vm\  = ξ\  - \ Vm\)\ , \ 
    where\ \ ξ\  = \(1\/2\) \((\ 
          Em\^2 - 1\ )\)\ \ \ we\ get\ the\ Newtonian\ form; \ \ ξ\  = \ \(\(1\/2\) \((dr\/dτ)\)\^2 + 
          Vm . \ \ \ \ \ \ dr/dτ\  = \ \(0\ for\ ξ = 
          Vm\)\)\)

Define Effective Potential

\(Expand[\(Em\^2 - 1\)\/2 - drdτ[r, Em, Lm, M]\^2\/2]\)

Vm[r_, Lm_, M_] := -(((Lm^2) M)/(r^3)) + ((Lm^2)/(2 (r^2))) - 
        (M/r)

We can get rid of the explicit M dependancy by expressing r and Lm in terms of M

Expand[Vm[r*M, Lm*M, M]]

With this definition of the effective potential we make a direct comparison to the Newtonian limit. The last two terms are the Newtonian potential while the negative inverse cube term is due to the metric. The Newtonian potential,

Vn[r_, Lm_, M_] := ((Lm^2)/(2 (r^2))) - (M/r)

Compare the Newtonian potential Vn to the Schwarzschild effective potential Vm

Show[
        Plot[{Vm[r M, 4.2 M, M], Vn[r M, 4.2 M, M]}, {r, 2, 40}, 
          Frame -> True, 
          FrameLabel -> {"r/M", None, None, None}, 
          PlotRange -> {{0, 40}, {-0.04, 0.05}}, 
          ImageSize -> {400, Automatic}, 
          PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0]}, 
          DefaultFont -> {"Courier", 14}, 
          DisplayFunction -> Identity], 
        Graphics[{Text["Vn", {10, 0.01}], Text["Vm", {8, - ..03}], 
            Text["Lm=4.2M", {20, 0.04}]}], 
        DisplayFunction -> $DisplayFunction];

Extrema of the potential

Solve for the extrema of the potential

Solve[(∂_r) Vm[r, Lm, M] == 0, r]

rmax[Lm_, 
        M_] := ((Lm^2)/(2  M)) ( 
          1 + (Sqrt[1 - 12  ((M^2)/(Lm^2)) ]))

rmin[Lm_, 
        M_] := ((Lm^2)/(2  M)) ( 
          1 - (Sqrt[1 - 12  ((M^2)/(Lm^2)) ]))

The minimum for both rmin and rmax occur for (Lm^2)/(M^2) = 12

rmin[(Sqrt[12]) M, M]

rmax[(Sqrt[12]) M, M]

For a given angular momentum Lm, rmin corresponds to inner unstable circluar orbits while rmax corresponds to outer stable circular orbits.  The transition occurs at  Lm = (Sqrt[12])M with rmin = rmax = 6M.

\(\(Show[{Plot[rmax[Lm*M, M]/M, {Lm, \@12, 6}, PlotStyle -> Red, 
            DisplayFunction -> Identity], 
          Plot[rmin[Lm, 1], {Lm, \@12, 6}, PlotStyle -> Blue, 
            DisplayFunction -> Identity], 
          Graphics[{Text["rmin", {5.5, 6}], 
              Text["rmax", {5.5, 23}]}]}, 
        ImageSize -> {400, Automatic}, Frame -> True, 
        DefaultFont -> {"Courier", 14}, 
        FrameLabel -> {"Lm/M", \*""\!\(r\/M\)"", None, None}, 
        RotateLabel -> False, 
        DisplayFunction -> $DisplayFunction];\)\)

The points indicate the extrema of Vm

Show[Plot[{Vm[r M, 6 M, M], Vm[r M, 5 M, M], 
              Vm[r M, 4  M, M], Vm[r M, (Sqrt[12]) M, M]}, {r, 2, 40}, 
            Frame -> True, Axes -> None, 
            FrameLabel -> {"r/M", "Vm", None, None}, 
            RotateLabel -> False, 
            PlotRange -> {{0, 34}, {-0.2, 0.4}}, 
            ImageSize -> {500, Automatic}, 
            PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0], 
                RGBColor[0, 0, 1]}, DefaultFont -> {"Courier", 14}, 
            DisplayFunction -> Identity],Graphics[{Text["Extrema of the effective potential",{20, 0.3}],Text["Lm",{5, 0.35}],Text["6M",{5, 0.3}],Text["5M",{5, 0.12}],Text["4M",{5, 0.01}],Text["(Sqrt[12])M",{5, -0.1}],PointSize[0.015],
Red,Point[{rmin[(Sqrt[12]), 1], Vm[rmin[(Sqrt[12]), 1], (Sqrt[12]), 1]}],Point[{rmax[(Sqrt[12]), 1], Vm[rmax[(Sqrt[12]), 1], (Sqrt[12]), 1]}],Blue,Point[{rmin[4, 1], Vm[rmin[4, 1], 4, 1]}],Point[{rmax[4, 1], Vm[rmax[4, 1], 4, 1]}],Green,Point[{rmin[5, 1], Vm[rmin[5, 1], 5, 1]}],Point[{rmax[5, 1], Vm[rmax[5, 1], 5, 1]}],Red,Point[{rmin[6, 1], Vm[rmin[6, 1], 6, 1]}],Point[{rmax[6, 1], Vm[rmax[6, 1], 6, 1]}]}],DisplayFunction -> $DisplayFunction];

Differential relationships: dtdts, drdrs

Local  relationships follows EBH:   dtdts[r, M] ,       drdrs[r, M].

From the metric:  d(τ^2) = (dts^2) - (drs^2) - (r^2)d(ϕ^2)  =  β (dt^2) - (dr^2)/β - (r^2)d(ϕ^2)
 dts and drs corresponds to the "shell" coordinates from EBH

dtdts[r_, 
          M_]:=β[r, M]^(-1/
              2)

dtdts[r, M]

drdrs[r_, M_] := β[r, M]^(1/2)

drdrs[r, M]

Relationship between proper time and local time for locally measured radial velocity v from special relativity;  d(τ^2) = (dts^2) - (drs^2) where drs/dts = v.  dtsdτ = γ[v]

γ[
          v_] = (((1 - (v^2)))^(-1/2));

dtsdτ[v_]:=γ[v]

Tangential component of velocity with respect to local time for an object with velocity v at angle a from outward radial:   dϕdts[r,v,a]

dϕdts[r_, v_, 
          a_] := ((v Sin[a Degree ])/r)

Tangential component of velocity with respect to coordinate time: dϕdt[r,Em,Lm,M]

dϕdt[r_, Em_, Lm_, M_] := 
      dϕdτ[r, Lm]/dtdτ[r, Em, M]

dϕdt[r, Em, Lm, M]

EM[r, v, M]

Radial velocity with respect to coordinate time:  drdt[r, Em, Lm, M, sign]

drdt[r_,  Em_,  Lm_,  M_]  := 
        drdτ[r,  Em,  Lm,  M]/dtdτ[r,  Em,  M]

drdt[r, Em, Lm, M]

Simplify[
      drdt[r, Em, Lm, M, 1] == β[r, M]/(Em ) 
          drdτ[r, Em, Lm, M, 1]]

We can get rid of the explicit M dependancy by expressing r and Lm in terms of M

drdt[r*M, Em, Lm*M, M]

Radial velocity with respect to  local time: drdts[r, Em, Lm, M, sign]

drdts[r_, Em_, Lm_, M_] := drdt[r, Em, Lm, M]*dtdts[r, M]

drdts[r, Em, Lm, M]

We can get rid of the explicit M dependancy by expressing r and Lm in terms of M

drdts[r*M, Em, Lm*M, M]

Radial local velocity as observed by local time:  drsdts[r.Em, Lm,M, sign]

drsdts[r_, Em_, Lm_, M_] := 
        dtdts[r, M]/drdrs[r, M]*drdt[r, Em, Lm, M]

drsdts[r, Em, Lm, M]

note that:

dtdts[r, M]/drdrs[r, M]

Combined velocities

Proper time frame

vτ[r_, v_, a_, 
        M_] := √(drdτ[r, EM[r, v, M], LM[r, v, a], M]^2 + 
            (r^2) dϕdτ[r, LM[r, v, a]]^2)

Distant frame

vt[r_, v_, a_, 
        M_] := √(drdt[r, EM[r, v, M], LM[r, v, a], M]^2 + 
            (r^2) dϕdt[r, EM[r, v, M], LM[r, v, a], M]^2)

Local frame

vts[r_, v_, a_, 
        M_] := (Sqrt[drdts[r, EM[r, v, M], LM[r, v, a], M]^2 + (r^2) dϕdts[r, v, a]^2])

purely tangential and radial velocities

Simplify[vτ[r, v, 90, M]^2]

Simplify[vτ[r, v, 0, M]^2]

Simplify[vτ[r, v, 45, M]^2]

Simplify[vt[r, v, 90, M]^2]

Simplify[vt[r, v, 0, M]^2]

Simplify[vts[r, v, 90, M]^2]

Simplify[vts[r, v, 0, M]^2]

Plot Routines

3D  coordinates (rCos(ϕ), rSin(ϕ), t):   rtphiplot[s0,s1,data]

Options: p-> track points, k -> displayfunction,ft->default font , pp->PlotPoints, vp-> viewpoint, br->BoxRatios

Options[rtphiplot] = {p -> 20, k -> $DisplayFunction, 
        pp -> Automatic, ft -> {"Courier", 14}, 
        vp -> {1, 0, 0}, br -> {1, 1, 1}}; 
    rtphiplot[s0_, s1_, data_, opts___] := 
      rtphiplotopt[s0, s1, 
        data, p /. {opts} /. Options[rtphiplot], k /. {opts} /. 
          Options[rtphiplot], pp /. {opts} /. 
          Options[rtphiplot], ft /. {opts} /. 
          Options[rtphiplot], vp /. {opts} /. 
          Options[rtphiplot], br /. {opts} /. Options[rtphiplot]]; 
    rtphiplotopt[s0_, s1_, data_, p_, k_, pp_, ft_, vp_, br_] := 
      Show[Table[
          Graphics3D[{RGBColor[1, 0, 0], PointSize[0.02], 
              Point[Flatten[{t[j*s1/20] /. data, 
                    r[s1*j/20] Cos[ϕ[s1*j/20]] /. data, 
                    r[s1*j/20] Sin[ϕ[s1*j/20]] /. data}]]}], {j, p, 
            20}], ParametricPlot3D[
          Evaluate[{t[s], r[s] Cos[ϕ[s]], r[s] Sin[ϕ[s]], 
                RGBColor[1, 0, 0]} /. data, {s, s0, s1}, 
            DisplayFunction -> Identity, PlotPoints -> pp]], 
        Axes -> True, AxesLabel -> {t, x, y}, DisplayFunction -> k, 
        DefaultFont -> ft, ViewPoint -> vp, 
        BoxRatios -> br];

2D  coordinates (rCos(ϕ), rSin(ϕ)):      rphiplot[s0,s1,data].

Options: p-> track points, k -> displayfunction, mb--> MaxBend, pd->PlotDivision, pp->PlotPoints

Remove[rphiplot]

Options[rphiplot] = {p -> 20, df -> $DisplayFunction, 
          pr -> Automatic, gl -> None, 
          ft -> {"Courier", 14}, mb -> 10, pd -> 30, 
          pp -> 100, h -> 300};
rphiplot[s0_, s1_, data_, opts___] := 
rphiplotopt[
          s0, s1, data, p /. {opts} /. 
            Options[rphiplot], df /. {opts} /. 
            Options[rphiplot], pr /. {opts} /. 
            Options[rphiplot], gl /. {opts} /. 
            Options[rphiplot], ft /. {opts} /. 
            Options[rphiplot], mb /. {opts} /. 
            Options[rphiplot], pd /. {opts} /. 
            Options[rphiplot], pp /. {opts} /. 
            Options[rphiplot], h /. {opts} /. 
            Options[rphiplot]];
\(\(\(rphiplotopt[s0_, s1_, data_, p_, df_, pr_, gl_, ft_, mb_, pd_, pp_, 
          h_] := Show[
          Graphics[{RGBColor[0, 1, 0], 
              Circle[{0, 0}, 2]}], 
		Graphics[{RGBColor[0, 1, 0], 
              PointSize[0.05], 
              Point[{0, 0}]}], 
		Graphics[{RGBColor[0, 0, 1], 
              PointSize[0.02], 
              Point[Flatten[{Chop[r[0] Cos[ϕ[0]]] /. data, 
                    Chop[r[0] Sin[ϕ[0]]] /. data}]]}], 
	
          Table[Graphics[{RGBColor[1, 0, 0], PointSize[0.02], 
                Point[Flatten[{Chop[r[j*s1/20] Cos[ϕ[j*s1/20]]] /. data, 
                      Chop[r[j*s1/20] Sin[ϕ[j*s1/20]]] /. data}]]}], {j, 
              p, 20}], 
          ParametricPlot[
            Evaluate[
              Chop[{r[s] Cos[ϕ[s]], r[s] Sin[ϕ[s]]}] /. data], {s, 
              s0, s1}, PlotStyle -> {RGBColor[1, 0, 0]}, MaxBend -> mb, 
            PlotDivision -> pd, PlotPoints -> pp, 
            DisplayFunction -> Identity], Axes -> True, Frame -> True, 
          DefaultFont -> ft, AxesLabel -> {x, y}, 
          AspectRatio -> Automatic, DisplayFunction -> df, 
          GridLines -> gl, PlotRange -> pr, 
          ImageSize -> {h, Automatic}];\)\(
\)
    \)\)

2D r and t : rtplot[s0,s1,data];

Options: track points->p,  df->$DisplayFunction, MaxBend->mb, PlotDivision->pd, PlotPoints->pp

Options[rtplot] = {p -> 20, df -> $DisplayFunction, 
          ft -> {"Courier", 14}, mb -> 10, pd -> 30, 
          pp -> 25};
rtplot[s0_, s1_, data_, opts___] := 
        rtplotopt[s0, s1, 
          data, p /. {opts} /. Options[rtplot], df /. {opts} /. 
            Options[rtplot], ft /. {opts} /. 
            Options[rtplot], mb /. {opts} /. 
            Options[rtplot], pd /. {opts} /. 
            Options[rtplot], pp /. {opts} /. 
            Options[rtplot]];
\(rtplotopt[s0_, s1_, data_, p_, df_, ft_, mb_, pd_, pp_] := 
      Show[
		Graphics[{RGBColor[0, 1, 0], PointSize[0.05], 
            Point[{0, 0}]}], 
		Graphics[{RGBColor[0, 0, 1], 
            PointSize[0.02], Point[Flatten[{t[0] /. data, r[0] /. data}]]}], 
        Graphics[{RGBColor[0, 1, 0], 
            Line[{{0, 2}, Flatten[{t[s1] /. data, 2}]}]}], 
	
        Table[	Graphics[{RGBColor[1, 0, 0], PointSize[0.02], 
              Point[Flatten[{t[j*s1/20] /. data, r[j*s1/20] /. data}]]}], {j, 
            p, 20}], 
        ParametricPlot[
          Evaluate[{t[s], r[s]} /. data, {s, s0, s1}, 
            PlotStyle -> {RGBColor[1, 0, 0]}, MaxBend -> mb, 
            PlotDivision -> pd, DefaultFont -> ft, 
            PlotPoints -> pp], DisplayFunction -> Identity], 
        Axes -> True, Frame -> True, PlotRange -> All, 
        AxesLabel -> {t, r}, DisplayFunction -> df]\)

Constants of the motion: EM, LM in terms of initial conditions

Angular momentum:   LM[r0, v0, a]

Angular momentum per unit mass(constant of the motion) in terms of initial  r0,v0 and initial angle a:  LM[r0, v0, a]  = (r0^2) dϕ0/dτ.

LM[r0_,  v0_,  a_] := ((r0)^2)dϕdts[r0, v0, a]*dtsdτ[v0];

LM[r0, v0, a]

Energy: EM[r0, v0, M]

Energy  per unit mass for an object initially at r0 moving with initial velocity v0.

EM[r0_, v0_, M_] := β[r0, M] dtdts[r0, M]*
        dtsdτ[v0]

EM[r0, v0, M]

Differentials in terms of initial conditions:

dt/dτ

dtdτ[r, EM[r0, v0, M], M]

dr/dτ

drdτ[r, EM[r0, v0, M], LM[r0, v0, a0], M]

dϕ/dτ

dϕdτ[r, LM[r0, v0, a]]

dr/dt

drdt[r, EM[r0, v0, M], LM[r0, v0, a0], M]

drs/dts

drsdts[r, EM[r0, v0, M], LM[r0, v0, a0], M]

Solutions using NDSolve

Setting the option M=1, scales r in units of M.

Explicit initial conditions: solve[r0, vr0, vt0, ϕ0, vϕ0].

Initial Conditions: r[0]=r0,r'[0]=vr0,ϕ[0]=,ϕ0,ϕ'[0]=vϕ0,t'[0]=0.

Options[solve] = {M -> 1, ms -> Infinity, 
          ag -> Automatic, wp -> MachinePrecision, 
          pg -> Automatic};
solve[s_, r0_, vr0_, ϕ0_, vϕ0_, vt0_, opts___] := 
      solveopt[s, r0, vr0, ϕ0, vϕ0, 
        vt0, M /. {opts} /. Options[solve], ms /. {opts} /. 
          Options[solve], ag /. {opts} /. 
          Options[solve], wp /. {opts} /. 
          Options[solve], pg /. {opts} /. Options[solve]]; 
    solveopt[s_, r0_, vr0_, ϕ0_, vϕ0_, vt0_, M_, ms_, ag_, wp_, 
        pg_] := NDSolve[{geor[M] == 0, geoϕ == 0, geot[M]  == 0, 
          r[0] == r0, r'[0] == 
            vr0, ϕ[0] == ϕ0 Degree, ϕ'[0] == 
            vϕ0, t[0] == 0, t'[0] == vt0}, {r, ϕ, t}, {τ, 
          0, s}, MaxSteps -> ms, AccuracyGoal -> ag, 
        WorkingPrecision -> wp, PrecisionGoal -> pg];

Parametic solution: gensolve[s, r, ϕ, v, a]

Here we define the initial conditions with r[0],ϕ[0], v0,a0:

Options[gensolve] = {M -> 1, rs -> 1, ms -> Infinity, 
          ag -> Automatic, wp -> MachinePrecision, 
          pg -> Automatic};
gensolve[s_, r0_, ϕ0_, v0_, a0_, 
          opts___] := 
gensolveopt[s, r0, ϕ0, v0, 
          a0, rs /. {opts} /. Options[gensolve], M /. {opts} /. 
            Options[gensolve], ms /. {opts} /. 
            Options[gensolve], ag /. {opts} /. 
            Options[gensolve], wp /. {opts} /. 
            Options[gensolve], pg /. {opts} /. 
            Options[gensolve]];
gensolveopt[s_, r0_, ϕ0_ , v0_, a0_, rs_, M_, ms_, ag_, wp_, 
          pg_] := 
NDSolve[{geor[M]   == 0, 
            geoϕ == 0, geot[M]  == 0, 
            r[0] == r0, r'[0] == 
              rs*drdτ[r0, EM[r0, v0, M], LM[r0, v0, a0 ], M], ϕ[
                0] == ϕ0 Degree, ϕ'[0] == 
              dϕdτ[r0, LM[r0, v0, a0 ]], 
            t[0] == 0, t'[0] == 
              dtdτ[r0, EM[r0, v0, M], M]}, {r, ϕ, t}, {τ, 0, 
            s}, MaxSteps -> ms, AccuracyGoal -> ag, 
          WorkingPrecision -> wp, 
          PrecisionGoal -> pg];

Constants of the motion:   cmsolve[s, r, ϕ, Em,Lm]

Here we define the initial conditions with r[0],ϕ[0], Em, Lm.

Options[cmsolve] = {rs -> 1, M -> 1, ms -> Infinity, 
          ag -> Automatic, wp -> MachinePrecision, 
          pg -> Automatic};
cmsolve[s_, r0_, ϕ0_, Em_, Lm_, 
          opts___] := 
cmsolveopt[s, r0, ϕ0, Em, 
          Lm, rs /. {opts} /. Options[cmsolve], M /. {opts} /. 
            Options[cmsolve], ms /. {opts} /. 
            Options[cmsolve], ag /. {opts} /. 
            Options[cmsolve], wp /. {opts} /. 
            Options[cmsolve], pg /. {opts} /. 
            Options[cmsolve]];

    
cmsolveopt[s_, r0_, ϕ0_ , Em_, Lm_, rs_, M_, ms_, ag_, wp_, 
          pg_] := 
NDSolve[{geor[M]   == 0, 
            geoϕ == 0, geot[M]  == 0, 
            r[0] == r0, r'[0] == 
              rs*drdτ[r0, Em, Lm, M], ϕ[
                0] == ϕ0 Degree, ϕ'[0] == 
              dϕdτ[r0, Lm], 
            t[0] == 0, t'[0] == dtdτ[r0, Em, M]}, {r, ϕ, 
            t}, {τ, 0, s}, MaxSteps -> ms, AccuracyGoal -> ag, 
          WorkingPrecision -> wp, 
          PrecisionGoal -> pg];

Radial motion

Falling from rest at infinity: Em = 1, Lm = 0.

EM/Em is unity for an object initially at rest at infinity and  remains unity. Note that our Newtonian ξ = ((Em^2)-1)/2 goes to zero at infinity as Em goes to 1

EM[r0, v0, M]

EM[∞, 0, M]

As the object falls from infinity to r its energy remains unity so v varies as,

Solve[EM[r, v, M] == 1, v]

If (v^2) = 2M/r then EM =1

EM[r, (Sqrt[((2  M)/r)]), M]

Radial derivitives

-drdτ[r, 1, 0, M]

-drsdts[r, 1, 0, M]

-drdt[r, 1, 0, M]

-drdts[r, 1, 0, M]

Thus for free fall from ∞, drs/dts == dr/dτ -> c while  drdt and drdts go to zero at the horizon.

\(\(Show[{Plot[{drdτ[r\ M, 1, 0, M], drdt[r\ M, 1, 0, M], 
              drdts[r\ M, 1, 0, M], drsdts[r\ M, 1, 0, M]}, {r, 2, 10}, 
            Frame -> True, 
            FrameLabel -> {"r/M", None, None, None}, 
            RotateLabel -> False, ImageSize -> {500, Automatic}, 
            PlotStyle -> {{RGBColor[1, 0, 0]}, {RGBColor[0, 1, 0]}, {RGBColor[
                    0, 0, 1]}, {RGBColor[1, 0, 1]}}, 
            DisplayFunction -> Identity]}, 
        Graphics[{Text[\*""(dr/d)τ=(drs/dts)"", {5,  ..7}], Text[\*""(dr/dt)"", {3.4, drdt[3, 1, 0, 1]}], 
            Text[\*""(dr/dts)"", {3.5,  .54}], 
            Text["Radial derivitives for fall from ∞: Em=1 and Lm=0", {6, 0.9}]}], DisplayFunction -> $DisplayFunction, 
        PlotRange -> {{2, 10}, {0, 1}}, 
        DefaultFont -> {"Courier", 12}, 
        ImageSize -> {500, Automatic}];\)\)

Falling from rest at r0:

Radial derivitives

drdτ[r, EM[r0, 0, M], 0, M]

drdt[r, EM[r0, 0, M], 0, M]

drsdts[r, EM[r0, 0, M], 0, M]

Plots of radial derivitives

\(\(Show[{Plot[{drdτ[r, EM[10, 0, 1], 0, 1], 
              drdt[r, EM[10, 0, 1], 0, 1], drdts[r, EM[10, 0, 1], 0, 1], 
              drsdts[r, EM[10, 0, 1], 0, 1]}, {r, 2, 10}, 
            PlotRange -> {{2, 10}, {0, 1}}, 
            PlotStyle -> {{RGBColor[1, 0, 0]}, {RGBColor[0, 1, 0]}, {RGBColor[
                    0, 0, 1]}, {RGBColor[1, 0, 1]}}, Frame -> True, 
            FrameLabel -> {"r", None, None, None}, 
            RotateLabel -> False, DisplayFunction -> Identity], 
          Graphics[Line[{{0, 0}, {1, 0}}]], 
          Graphics[{Text[\*""(dr/d)τ"", {2.7, 0.65}], 
              Text[\*""(dr/dt)"", {4.5, 0.25}], 
              Text[\*""(dr/dts)"", {2.5,  .46}], 
              Text[\*""(drs/dts)"", {5.5,  .52}], 
              Text["Radial derivitives for a fall from r=10: Em=1,Lm=0", {6,  .9}]}]}, DisplayFunction -> $DisplayFunction, 
        ImageSize -> {500, Automatic}, 
        DefaultFont -> {"Courier", 12}];\)\)

Acceleration due to gravity at r0

Earth surface using Newtonian physics

Quit[]

<< Miscellaneous`PhysicalConstants`

<< Miscellaneous`Units`

((EarthMass*GravitationalConstant)/(EarthRadius^2))

For the Schwarzschild metric

For the locally observed acceleration we need the second derivitive of rs with respect to ts

gr[r0_, M_:  1] := 
      ReplaceAll[
        ReplaceAll[(∂_ts) 
            drsdts[r[ts], EM[r0, 0, M], 0, M], r'[ts] -> 
            drdrs[r[ts], M] drsdts[r[ts], EM[r0, 0, M], 0, M]], 
        r[ts] -> r0]

gr[r0, M]

gr has units (Meter^(-1)) .   At the earth' 
      s surface in geometrized units,

gr[EarthRadius /Meter, 0.004435]/Meter

In conventional units we multiply by (c^2),  
    which gives Meter/(Second^2)

c = SpeedOfLight

grc[r_, M_] := gr[r, M]*(c^2)

grc[EarthRadius/Meter, 0.004435 ]/Meter

Tidal forces

We look at the gradient of gravitational forces

grad[r_, M_] = (∂_r) gr[r, M]

grad has units 1/(Meter^2)

grad[EarthRadius/Meter, 0.004435]/(Meter^2) // N

In conventional units

grad[EarthRadius/Meter, 0.004435] (c^2)/(Meter^2) // N

\(This\ is\ Acceleration/Meter\ or\ \(Meter/Second\^2\)/Meter . \ Thus, 
    at\ the\ earth' s\ surface, 
    gravitational\ acceleration\ is\ changing\ by\ approximately\ 3  
        ppm\ per\ meter . \ Now\ consider\ a\ black\ hole\ with\ \ M\  = \ 
      10000\ meters . \ 
          At\ what\ distance\ will\ the\ change\ in\ acceleration\ per\ meter\ equal\ that\ at\ the\ earth' s\ \(\(surface\)\(.\)\)\)

grad[5.685*10^6, 10000 ]*(c^2) (Second^2)/(Meter^2) // N

This means that at this distance the acceleration is increasing(tidal forces) by one g every meter.

Falling with initial  inward  radial velocity v0 at r0 ;

Energy

Energy constant of the motion EM(r0,v0,M=1) depends on initial velocity v0 and position r0 for a given mass and goes to zero at r0 = 2M

EM[r, v, M]

\(Show[
      Plot[Evaluate[Table[EM[r0\ M, v0, M], {v0, 0, 0.9,  .3}]], {r0, 2, 
          10}, 
PlotStyle -> {{RGBColor[1, 0, 0], 
              AbsoluteDashing[{1, 10}], 
              Thickness[0.007]}, {RGBColor[0, 1, 0], 
              AbsoluteDashing[{1, 10}], 
              Thickness[0.007]}, 
{RGBColor[0, 0, 1], 
              AbsoluteDashing[{1, 10}], Thickness[0.007]}}, 
        Frame -> 
          True, 
FrameLabel -> {"r0/M", "EM", 
            None, None}, 
RotateLabel -> False, 
        DefaultFont -> {"Courier", 12}, 
        ImageSize -> {500, 
            Automatic}, 
DisplayFunction -> Identity], 
      Graphics[{Text["v0=0", {6, 0.75}], Text["0.3", {6, 0.9}], 
          Text["0.6", {6, 1.1}], Text["v0=0.9", {6, 1.8}], 
          Text["EM for Initial Conditions:r0/M,v0", {6, 
              0.4}]}], 
DisplayFunction -> $DisplayFunction]\)

Radial derivitives

dr/dτ

drdτ[r, Em, 0, M, ]

drdτ[r, Em, 0, M]

drdτ[r, Em, 0, M, -1]

drdτ[r, EM[r0, v0, M], 0, M, -1]

drdτ[r, Limit[EM[r0, v0, M], r0 -> ∞], 0, 
      M, -1]

dr/dt

drdt[r, EM[r0, v0, M], 0, M, -1]

drdt[r, Limit[EM[r0, v0, M], r0 -> ∞], 0, 
      M, -1]

drs/dts

drsdts[r, EM[r0, v0, M], 0, M, -1]

drsdts[r, Limit[EM[r0, v0, M], r0 -> ∞], 0, 
      M, -1]

Radial derivitives for initial v0=0.4 from r0=10M

Show[{Plot[{drdτ[r M, EM[10 M, 0.4, M], 0, M], 
              drdt[r M, EM[10  M, 0.4, M], 0, M], 
              drdts[r M, EM[10  M, 0.4, M], 0, M], 
              drsdts[r M, EM[10  M, 0.4, M], 0, M]}, {r, 2, 10}, 
            PlotRange -> {{2, 10}, {0, 1}}, 
            PlotStyle -> {{RGBColor[1, 0, 0]}, {RGBColor[0, 1, 0]}, {RGBColor[
                    0, 0, 1]}, {RGBColor[1, 0, 1]}}, Frame -> True, 
            GridLines -> {None, {0.4}}, 
            FrameLabel -> {"r/M", None, None, None}, 
            RotateLabel -> False, DisplayFunction -> Identity], 
          Graphics[Line[{{0, 0}, {1, 0}}]], 
          Graphics[{Text[(dr/d)τ, {4.5, 0.7}], Text[(dr/dt), {4.5,  ..3}], 
              Text[(dr/dts), {2.7,  .5}], Text[(drs/dts), {3.5, 0.65}], 
              Text["Initial Conditions: v0=0.4 at r0=10M", {6,  ..9}]}]}, 
        DisplayFunction -> $DisplayFunction, 
        ImageSize -> {500, Automatic}, 
        DefaultFont -> {"Courier", 12}];

Radial derivitives for initial v0=0.6 from r0=10M

Show[{Plot[{drdτ[r*M, EM[10*M, 0.6, M], 0, M], 
              drdt[r*M, EM[10*M, 0.6, M], 0, M], 
              drdts[r*M, EM[10*M, 0.6, M], 0, M], 
              drsdts[r*M, EM[10*M, 0.6, M], 0, M]}, {r, 2, 10}, 
            PlotRange -> {{2, 10}, {0, 1}}, 
            PlotStyle -> {{RGBColor[1, 0, 0]}, {RGBColor[0, 1, 0]}, {RGBColor[
                    0, 0, 1]}, {RGBColor[1, 0, 1]}}, Frame -> True, 
            GridLines -> {None, {0.6}}, 
            FrameLabel -> {"r/M", None, None, None}, 
            RotateLabel -> False, DisplayFunction -> Identity], 
          Graphics[Line[{{0, 0}, {1, 0}}]], 
          Graphics[{Text[(dr/d)τ, {4.5, 0.9}], Text[(dr/dt), {4.5,  ..35}], 
              Text[(dr/dts), {2.7,  .5}], Text[(drs/dts), {3.5, 0.75}], 
              Text["Initial conditions: v0=0.6 at r0=10", {7, 0.2}]}]}, 
        DisplayFunction -> $DisplayFunction, 
        ImageSize -> {500, Automatic}, 
        DefaultFont -> {"Courier", 12}];

rs'[ts] plots for initial conditions v0 from r0=10M

\(Show[
      Plot[Evaluate[
          Table[drsdts[r*M, EM[10  M, v0, M], 0, M], {v0, 0, 4/5, 1/5}]], {r, 
          2, 10}, ImageSize -> {500, Automatic}, 
        PlotStyle -> {{RGBColor[1, 0, 0]}, {RGBColor[0, 1, 0]}, {RGBColor[0, 
                0, 1]}}, DefaultFont -> {"Courier", 12}, 
        Frame -> True, 
        FrameLabel -> {"r/M", \*""\!\(drs\/dts\)""}, 
        RotateLabel -> False, DisplayFunction -> Identity], 
      Graphics[{Text["v0=0", {9, 0.13}], Text["v0=0.2", {9, 0.3}], 
          Text["v0=0.4", {9, 0.45}], Text["v0=0.6", {9, 0.65}], 
          Text["v0=0.8", {9, 0.85}], 
          Text["Initial conditions: r0=10,v0", {5, 0.15}]}], 
      DisplayFunction -> $DisplayFunction]\)

r'[t] plots for initial conditions v0 from r=10M

\(\(Show[
        Plot[Evaluate[
            Table[drdt[r*M, EM[10  M, v0, M], 0, M], {v0, 0, 4/5, 1/5}]], {r, 
            2, 10}, PlotRange -> {0, 0.7}, 
          ImageSize -> {500, Automatic}, 
          PlotStyle -> {{RGBColor[1, 0, 0]}, {RGBColor[0, 1, 0]}, {RGBColor[
                  0, 0, 1]}}, DefaultFont -> {"Courier", 12}, 
          Frame -> True, 
          FrameLabel -> {"r/M", \*""\!\(dr\/dt\)""}, 
          RotateLabel -> False, DisplayFunction -> Identity], 
        Graphics[{Text["v0=0", {9, 0.15}], Text["v0=1/5", {9, 0.22}], 
            Text["v0=2/5", {9, 0.35}], Text["v0=3/5", {9, 0.51}], 
            Text["v0=0.8", {9, 0.65}], 
            Text["Initial conditions: r0=10,v0", {5, 0.15}]}], 
        DisplayFunction -> $DisplayFunction];\)\)

r'[τ] plots  initial conditions v0 from r0=10M

\(\(Show[
        Plot[Evaluate[
            Table[drdτ[r*M, EM[10  M, v0, M], 0, M], {v0, 0, 4/5, 
                1/5}]], {r, 2, 10}, PlotRange -> {0, 1.5}, 
          ImageSize -> {500, Automatic}, 
          PlotStyle -> {{RGBColor[1, 0, 0]}, {RGBColor[0, 1, 0]}, {RGBColor[
                  0, 0, 1]}}, DefaultFont -> {"Courier", 12}, 
          Frame -> True, 
          FrameLabel -> {"r/M", \*""\!\(dr\/dτ\)""}, 
          RotateLabel -> False, DisplayFunction -> Identity], 
        Graphics[{Text["v0=0", {9, 0.1}], Text["v0=1/5", {9, 0.3}], 
            Text["v0=2/5", {9, 0.45}], Text["v0=3/5", {9, 0.74}], 
            Text["v0=4/5", {9, 1.25}], 
            Text["Initial conditions: r0=10,v0", {5, 0.15}]}], 
        DisplayFunction -> $DisplayFunction];\)\)

dr/dτ exceeding unity

dr/dτ will exceed unity anytime v exceeds 1/(Sqrt[2]) ,regardless of the presence of a black hole.  Far from a black hole, dr/dt is just a particle's velocity which never exceeds unity.

Simplify[Limit[drdt[r*M, EM[r, v, M], 0, M], r -> ∞], 
      1 ≤ v ≥ 0]

The proper velocity is obtained through multiplying by the time dilation factor, i.e. dt/dτ.

Limit[dtdτ[r*M, EM[r, v, M], M], r -> ∞]

Thus

Solve[
      Limit[drdτ[r, EM[r, v, M], 0, M], r -> ∞] == 
        1, v]

The "proper" velocity of a photon is infinite since its proper time is zero.

Outward radial trajectory

For example with drs/dts= v0 = 1/5,  r0 = 10, ϕ = Pi/4, a = 0, the energy is less than unity.

Solve[drsdts[10  M, Em, 0, M] == 1/5, Em]

drdτ[10  M, (Sqrt[(5/6)]), 0, M]

dtdτ[10  M, (Sqrt[(5/6)]), M]

Numerical solution

radial1 = gensolve[69, 10, 45, 1/5, 0]

radial1a = gensolve[69, 10, 45, 1/5, 0, M -> 2]

rphiplot[0, 68, radial1, p -> 1]

Find the the turnaround point.

FindMaximum[r[s] /. radial1[[1]], {s, 0}]

The energy remains constant: v=0 at the turnaround r=12.

EM[12, 0]

EM[10, 1/5]

Compare r'[s] and r'[t]

Plot[{r'[s] /. radial1, r'[s]/t'[s] /. radial1},{s, 0, 68.08},Frame -> True,PlotRange -> {{0, 68}, {-1 ,  .2}},PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0]},FrameLabel->{"s","velocity","r'[s] and r'[t]",""},GridLines -> Automatic]

Escape from a black hole:  Em=1  -> v0  = (Sqrt[((2  M)/r0)])

Solve for escape velocity

Solve[EM[r, v, M] == 1, v]

vx[r0_, M_:  1] := (Sqrt[((2  M)/r0)])

drsdts[r0*M, 1, 0, M]

Numerical solution for escape from r0 =8M

vx[8]

escape1 = gensolve[10^12, 8, 45, 1/2, 0]

LogPlot[Evaluate[r'[τ] /. escape1], {τ, 0, 10^12}, 
      Frame -> True, GridLines -> Automatic, 
      PlotRange -> All]

Circular orbits

minimum stable orbit : rmin = rmax:  L = (Sqrt[12])M

Determine initial conditions

The maximum and minimun of the effective potential occur at rmin=rmax

Solve[rmin[L, M] == rmax[L, M], L]

Plot[Vm[r*M, (Sqrt[12]) M, M], {r, 3, 12}]

This minimum circular orbit has r = 6:

rmax[(Sqrt[12]) M, M]

Velocity in the minimum orbit

Solve[LM[6 M, v, 90] == (Sqrt[12]) M, v]

Angular momentum in the minimum orbit

LM[6 M, 1/2, 90]

Energy in the minimum orbit

EM[6  M, 1/2, M]

The locally observed orbital velocity in the minimum orbit

6  M dϕdts[6 M, 1/2, 90]

The distant observed orbital velocity in the minimum orbit:

?dϕdt

6  M dϕdt[6  M, ((2 (Sqrt[2]))/3), (Sqrt[12]) M, M]

The observed orbital velocity in the minimum orbit:

6  M dϕdτ[6  M, (Sqrt[12]) M]

Numerical solution

mincircle = gensolve[100, 6, 0, 1/2, 90]

rphiplot[0, 50, mincircle, p -> 1, n -> 20];

Circular orbit for L=4

Determine initial conditions

rmin[8 , 2]

Solve[LM[12, v, 90] == 4, v]

EM[12, 1/(Sqrt[10]), 1]

Numerical solution

circle2 = gensolve[200, 12, 45, 1/(Sqrt[10]), 90]

rphiplot[0, 200, circle2];

Precessing orbits

precessing orbits near circluar:

r0=12,v0=1/3

LM[12, 1/3, 90]

EM[12, 1/3, 1]

ps1 = cmsolve[1000, 12, 0, ((Sqrt[15])/4), 3 (Sqrt[2])];

rphiplot[0, 1000, ps1, pp -> 100];

highly elliptical orbits

LM[12, 2/5, 90] // N

EM[12, 2/5, 1]

drdτ[237, EM[12, 2/5, 1], LM[12, 2/5, 90], 1]

12  dϕdτ[12, LM[12, 2/5, 90]]*(Sqrt[1 - 4/25])

ps3 = 
        cmsolve[100000, 12, 0, EM[12, 2/5, 1], LM[12, 2/5, 90]];

rphiplot[0, 100000, ps3];

Unstable orbits:   Em = Maximum Vm

Determine r and v  for given L which gives an unstable orbit

r is determined by maximizing Vm for a given L while v is determined from Em = MaxVm

unstablevelocity[L_] := 
      Solve[EM[r /. Maximize[Vm[r, L], 2 < r, r][[2]], 
            v] == Maximize[Vm[r, L], 2 < r, r][[1]], 
        v]

L=4

unstablevelocity[4]

r /. Maximize[Vm[r, 4], 2 < r, r][[2]]

a /. Maximize[Vm[4, LM[4, 1/(Sqrt[2]), a]], a][[2]] // N

Initial conditions for exact  position  r = 4

unstable1 = gensolve[10000, 4, 0, (1/(Sqrt[2])), 90]

rphiplot[9900, 10000, unstable1];

Initial conditions for machine precision  position  rmin = 4.0  for L = 4

unstable2 = gensolve[100, 4.0, 0, (1/(Sqrt[2])), 90]

rphiplot[0, 100, unstable2];

Photon trajectories

Constants of the motion for photons and radial motion

Constants of the motion

Photons have the same energy and angular momentum constants of the motion as massive particles due to the same symmetry conditions.

dtdλ[r_, Em_, M_] := ((Em r)/(2 M - r))

dϕdλ[r_, Lm_] := (Lm/(r^2))

From the Schwarzschild notebook:

ReplaceAll[MetricS . V . V, θ -> (π/2)]

Radial velocity

Photons differ from massive particles by having null geodesics:u . u=(Subscript[g, α])β((dx^α)/d)λ((dx^β)/d)λ =0. λ is the affine parameter.

Solve[
        ReplaceAll[MetricS . V . V, θ -> (π/2)] == 
          0, r'[τ]];
r'[τ] /. %;
drdλ[r_, Em_, Lm_, M_:  1] = 
      FullSimplify[
        ReplaceAll[%[[2]], {t'[τ] -> 
              dtdλ[r, Em, M], ϕ'[τ] -> 
              dϕdλ[r, Lm]}], 
        Assumptions -> (M ≥ 0 && 
              r ≥ 2  M)]

Due to the null condition,λ can multiplied by any constant and the equation is unchanged. Looking at dϕ/dλ and dtdλ and changing λ to kλ divides Em and Lm by the same constant and so only their ratio has physical significance. We form an expression containing only their ratio.

\(Simplify[drdλ[r, Em, Lm, M]\^2\/Lm\^2]\)

With b = Lm/Em this can be written as (1/(b^2))=(1/(Lm^2))((((dr/d)λ))^2)+ ((r - 2  M)/(r^3))  which is Hartle's(9.63). By analogy with ((((dr/d)τ))^2) = (Em^2)-(V^2) for massive particles,  (1/(b^2)) acts as the energy while the last is an effective potential for photons.

Veff[r_, M_:  1] := (Sqrt[((r - 2 M)/(r^3))])

Note: Hartle uses the square as the potential to more easily compare to Newtonian gravity.

Impact parameter b[r0, a0, M]

From EBH , b = L / p = L/(Sqrt[(E^2) -  (m^2)])  which goes to LM/EM as m -> 0. 
 LM/EM is obtained by letting v -> 1.

EM[r0, v0, M]

LM[r0, v0, a0]

Limit[Evaluate[LM[r0, v0, a0]/EM[r0, v0, M]], v0 -> 1]

b[r0_, a0_, 
        M_:  1] := ((r0 Sin[a0 Degree ])/(Sqrt[1 - ((2 M)/r0)]))

Distant observers: drldt[r, bl, M] , dϕdt[r,Em.Lm,M]

For a massive particle

drdt[r, Em, Lm, M]

This can be converted to massless particles by substituting Lm=bl*Em and letting the energy per unit mass go to infinity.

Limit[ReplaceAll[drdt[r, Em, Lm, M], Lm -> bl Em], 
      Em -> ∞]

dr/dt for light

drldt[r_, bl_, 
        M_:  1] := (1 - ((2 M)/r)) (Sqrt[1 - (1 - ((2 M)/r)) ((bl^2)/(r^2))])

The equivalent equation for dϕ/dt in the case of photons is to collect Lm/Em as the impact parameter b.

dϕdt[r, Em, Lm, M]

ReplaceAll[dϕdt[r, Em, Lm, M], Lm -> bl Em]

dϕldt[r_, bl_, 
        M_:  1] := (bl/(r^2)) (1 - ((2 M)/r))

Substituting for bl.

dϕldt[r, b[r0, a0, M], M]

Resulting light speed as measured by distant observers:  (Sqrt[(((dr/dt))^2) + (((rdϕ/dt))^2)])

lightspeed[r_, bl_, 
        M_:  1] := √(drldt[r, bl, M]^2   + 
            r^2*dϕldt[r, bl, M]^2)

lightspeed[r, bl, M]

Observed light speed goes to zero at r = 2M

Plot[
        Evaluate[Table[lightspeed[r, b[20, a0]], {a0, 0, 90, 10}]], {r, 2, 
          20}, Frame -> True, GridLines -> Automatic, 
        PlotLabel -> "Launch from r0 = 20 with varying impact \
parameters", 
        PlotStyle -> 
          Evaluate[Table[Hue[i/10], {i, 1, 9}]]];

Local observers: drlsdts[r, bl, M], dϕldts[r, b1, M]

drlsdts[r_, b1_, M_:  1] := (drldt[r, b1, M]*dtdts[r, M])/
        drdrs[r, M]

drlsdts[r, b1, M]

dϕldts[r_, b1_, M_:  1] := 
      dϕldt[r, b1, M]*dtdts[r, M]

dϕldts[r, b1, M]

Light speed is always unity as measured locally

√(drlsdts[r, b1, M, 1]^2  +  (r^2) 
            dϕldts[r, b1, M]^2)

Another impact parameter derivation: Consider a photon launched at angle a0 from the outward radial at some r0.  The tangential component is then Sin(a0) which must  equal rdϕ/dts as observed locally.

Solve[r dϕldts[r, bl, M] == Sin[a0], bl]

Unstable circular orbits at the critical point

For what value of b is the local radial velocity zero? This must correspond to a circular orbit .

Flatten[Solve[drlsdts[r, bl, 1, M]^2 == 0, bl]][[2, 
        2]]

blc[r_, M_:  1] := (r/(Sqrt[1 - ((2 M)/r)]))

Plot[b1c[r, 1], {r, 2, 10}]

Solve[(∂_r) b1c[r, 1] == 0, r]

b1c[3, 1]

Initial conditions  for NDSolve

The initial conditions are determined by speed of light =1 and  dτ =0.                 
               d(τ^2)   =  (1-((2  M)/r))(dt^2) - (((1-((2  M)/r)))^(-1))(dr^2) - (r^2)d(ϕ^2) = 0
               t'[λ]^2 = (((1-((2  M)/r)))^(-2))r'[λ]^2  + (r^2)(((1-((2  M)/r)))^(-1))ϕ'[λ]^2
 speed of light unity      r'[λ]^2 + r[λ]^2 ϕ'[λ]^2 = 1
 The two components of the light velocity can be expressed in terms of the launch angle.

drds[a0_] := Cos[a0 Degree ]

dϕds[r0_, a0_] := Sin[a0 Degree]/r0

dtds[r0_, a0_, 
        M_:  1] := √((((1 - ((2 M)/r0)))^(-2)) 
              drds[a0]^2 + (r0^2) (((1 - ((2 M)/r0)))^(-1)) 
              dϕds[r0, a0 ]^2)

Tangential launch

dtds[r0, 90]

Radial launch

dtds[r0, 0]

NDSolve for photons: lightsolve[s, r0, ϕ0, a0]

lightsolve[s1_, r0_, ϕ0_, a0_] := 
      solve[s1, r0, drds[a0], ϕ0 , dϕds[r0, a0], 
        dtds[r0, a0]]

Unstable circular orbit at the critical point

b[3,90]=(Sqrt[27]);  tangental velocity at critical point, r0=3.

critical1 = lightsolve[200, 3, 0, 90];

rphiplot[0, 130, critical1]

critical2 = lightsolve[200, 3 + 10^-6, 0, 90];

rphiplot[0, 40, critical2]

\(\(\(Consider\ a\ photon\ launched\ \ \ from\ \ some\ large\ distance\ r0\ at\ ϕ0\  = \ \(0\ \ with\ an\ impact\ parameter\ bl = \ \(3 \@ 3 . \ 
              Since\ r0 = \@\(y\^2 + 3\^2\)\ the\ photon\ launch\ angle\ must\ be\)\)\)\(\ \)\(,\)\)\)

ϕ0c[r0_] := 
        ArcSin[((3 (Sqrt[ 3]))/(Sqrt[(r0^2) - 9]))]*180/π;

ϕ0c[40] // N

180 - ϕ0c[40] // N

critical3 = lightsolve[200, 40, 0, 180 - ϕ0c[40]];

rphiplot[0, 200, critical3, pr -> {{-4, 40}, {-4, 20}}, 
      gl -> Automatic]

Changing the launch angle slightly

critical[p_] := 
        lightsolve[200, 40, 0, 180 - ϕ0c[40] + p];

critical[24397004/1000000000];

rphiplot[0, 200, %, pr -> {{-4, 4}, {-4, 4}}, 
      gl -> Automatic]

Array of photons

We want a uniform distribution of initial photons along y = constant with x  = -max to plus max step p.  To obtain the required   r0 , ϕ0 and a0 initial conditions we must solve for r0 and ϕ0.

Remove[ic, rs, rv, as, av]

Off[Solve::"ifun"]

ic[x_, y_] := 
      Solve[{r0 Sin[ϕ0] == y, 
          r0 Cos[ϕ0] == x}, {r0, ϕ0}, r0]

We retain only the positive r0 solutions and the corresponding ϕ0.

rs = 
      Table[Flatten[Table[ic[x, 40], {x, -40, 40, 5}]][[ 2  i - 1, 
            2]], {i, 20}]

To control how close we get to the event horizon we include a limiting sum parameter n.

rv[n_] := 
      Join[Table[rs[[i]], {i, 2, n, 2}], 
        Reverse[Table[rs[[i]], {i, 2, n, 2}]]]

rv[16]

as = 
      Table[Flatten[Table[ic[x, 40], {x, -40, 40, 5}]][[2 i, 
                2]], {i, 20}]*180/Pi // N

av[n_] := 
      Join[Table[as[[i]], {i, 2, n, 2}], 
        Reverse[Table[180 - as[[i]], {i, 2, n, 2}]]]

av[16]

Impact parameters

Table[
      b[rv[16][[i]], 270 - av[16][[i]]], {i, 1, 
        16}]

(Sqrt[27]) // N

Off[NDSolve::"ndsz"]

array1 = 
        Table[lightsolve[100, rv[16][[i]], av[16][[i]], 
            270 - av[16][[i]]], {i, 1, 16}];

Off[InterpolatingFunction::"dmval"]

Show[
      Table[rphiplot[0, 100, array1[[i]], df -> Identity], {i, 1, 16}], 
      DisplayFunction -> $DisplayFunction]