Relativistic Orbits and Black HolesFrank Wang, fw@phys.columbia.eduupdated 5/20/2004Using Maple to derive equations of motion from the Lagrangian, to solve differential equations numerically, and to form graphs based on numerical solutions, we present trajectories of a particle under gravity based on Newtonian theory and general relativity with the Schwarzschild metric and the Kerr metric. The Lagrangian for the Kerr solution in general relativity, which describes a spherical object of mass NiMlIk1H and spin angular momentum NiMlIkpH , isNiMvJSJMRywqKigsJiIiIkYoKigiIiNGKCUiTUdGKCUickchIiJGLUYoKiQtJSVkaWZmRzYkLSUidEc2IyUkdGF1R0Y1RipGKEYqRi1GLSouRipGKCUiYUdGKEYrRihGLEYtRi9GKC1GMDYkLSUkcGhpR0Y0RjVGKEYtKihGLEYqKiZGKkYoJSZEZWx0YUdGKEYtLUYwNiQtRixGNEY1RipGKCooLCgqJEYsRipGKCokRjdGKkYoKipGKkYoRitGKEY3RipGLEYtRihGKCokRjhGKkYoRipGLUYo where NiMvJSJhRyomJSJKRyIiIiUiTUchIiI= and NiMvJSZEZWx0YUcsKCokJSJyRyIiIyIiIiooRihGKSUiTUdGKUYnRikhIiIqJCUiYUdGKEYp .If NiMvJSJhRyIiIQ== , namely the object does not rotate, the Kerr solution becomes the Schwarzschild solution. In the limit of NiMvJSJjRyUpaW5maW5pdHlH and NiMvJSJ0RyUkdGF1Rw== , the Schwarzschild solution further reduces to the Lagrangian for Newtonian mechanics, i.e., the difference between kinetic energy and potential energy. To find equations of motion, we simply employ the Euler-Lagrange equationNiMqJiUiZEciIiIlI2R0RyEiIg== NiMtJSVkaWZmRzYkJSJMRyYlI3EnRzYjJSJpRw== - NiMtJSVkaWZmRzYkJSJMRyYlInFHNiMlImlH = 0where NiMlInFH is a coordinate and NiMqJiUicUciIiIlIidHRiU= denotes its derivative with respect to time: NiMqJiUicUciIiIlIidHRiU= = NiMqJiUjZHFHIiIiJSNkdEchIiI= (conventional notation for this quantity as NiMlInFH with a dot above this letter).For additional information, see F. Wang, "Relativistic orbits with computer algebra," Am. J. Phys. 72, (2004 in press) and references therein; contact the author for preprint. By varying the initial conditions, one can discover many interesting orbits in relativity. I thank Professor Roman Koniuk (koniuk@yorku.ca) for suggesting the last three sections. restart: with(plots): with(plottools):<Text-field layout="_pstyle10" style="_cstyle11">Newtonian mechanics</Text-field>First define the kinetic energy, the potential energy, and the Lagrangian: T := 1/2*(diff(r(t),t)^2 + r(t)^2*diff(phi(t),t)^2);V := -G*M/r(t);L := T - V;L1 := subs({r(t)=var1, diff(r(t),t)=var2, phi(t)=var3, diff(phi(t),t)=var4}, L):Employ the Euler-Lagrange equation for the Newtonian Lagrangian (Epr11 through Eq26) Epr11 := diff(L1, var4):Epr12 := diff(L1, var3):Epr13 := subs({var1=r(t),var2=diff(r(t),t),var3=phi(t),var4=diff(phi(t),t)},Epr11):Eq14 := Epr13 = l;Epr21 := diff(L1, var2):Epr22 := diff(L1, var1):Epr23 := subs({var1=r(t), var2=diff(r(t),t), var3=phi(t), var4=diff(phi(t),t)}, Epr21):Epr24 := subs({var1=r(t), var2=diff(r(t),t), var3=phi(t), var4=diff(phi(t),t)}, Epr22):Epr25 := diff(Epr23,t):Eq26 := Epr25 - Epr24 = 0;Eq31 and Eq32 decouple the differential equations.Eq31 := isolate(Eq14, diff(phi(t),t));Eq32 := eval(Eq26, Eq31);Provide the initial values, and use Maple to solve these differential equations numerically and to form plots based on numerical solutions. G:=1; M:=1;Eq41 := r(0) = 26;Eq42 := D(r)(0) = 0;Eq43 := phi(0) = 0;Eq44 := D(phi)(0) = 0.0071;En := eval(T + V, {r(t)=rhs(Eq41), diff(r(t),t)=rhs(Eq42), diff(phi(t),t)=rhs(Eq44)});l := eval(lhs(Eq14), {r(t)=rhs(Eq41), diff(phi(t),t)=rhs(Eq44)});The eccentricity is not needed; it is listed as a reference.epsilon := sqrt(1 + 2*En*l^2/((G*M)^2)); ini1 := Eq41, Eq42, Eq43;Eq51:=dsolve({Eq31, Eq32, ini1}, {r(t), phi(t)}, numeric, output=listprocedure);polarplot([rhs(Eq51(t)[3]), rhs(Eq51(t)[2]), t=0..720], scaling=constrained, legend="Newtonian");pnewton := polarplot([rhs(Eq51(t)[3]), rhs(Eq51(t)[2]), t=0..720], scaling=constrained, color=black, legend="Newtonian", linestyle=3, thickness=2):<Text-field layout="_pstyle10" style="_cstyle11">Schwarzschild metric</Text-field>G := 'G': M := 'M':The procedure for the Schwarzschild metric is almost identical to the Newtonian mechanics; we merely change the Lagrangian.f := 1/2*(-(c^2 - 2*G*M/r(tau))*diff(t(tau),tau)^2
+ 1/(1 - 2*G*M/(c^2*r(tau)))*diff(r(tau),tau)^2 + r(tau)^2*diff(phi(tau),tau)^2);f1 := subs({t(tau)=var1, diff(t(tau),tau)=var2, r(tau)=var3, diff(r(tau),tau)=var4,
phi(tau)=var5, diff(phi(tau),tau)=var6}, f):Employ the Euler-Lagrange equations (Epr11 through Eq34)Epr11 := diff(f1, var2):Epr12 := diff(f1, var1):Epr13 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),
var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr11):Eq14 := Epr13 = -a;Epr21 := diff(f1, var4):Epr22 := diff(f1, var3):Epr23 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),
var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr21):Epr24 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),
var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr22):Epr25 := diff(Epr23, tau):Eq26 := Epr25 - Epr24 = 0;Epr31 := diff(f1, var6):Epr32 := diff(f1, var5):Epr33 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),
var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr31):Eq34 := Epr33 = h;Decouple the differential equations in the following three lines.Eq53 := isolate(Eq14, diff(t(tau),tau));Eq54 := isolate(Eq34, diff(phi(tau),tau));Eq55 := subs({Eq53, Eq54}, Eq26);Provide the initial conditions, and plot the numerical solutions.M := 1; G := 1; c := 1;Eq77 := r(0) = 26; Eq78 := D(r)(0) = 0;Eq79 := phi(0) = 0;Eq80 := D(phi)(0) = 0.0071; h := eval(lhs(Eq34), {r(tau)=rhs(Eq77), diff(phi(tau),tau)=rhs(Eq80)}); a := sqrt((1 - 2*M/rhs(Eq77))*(1 + h^2/rhs(Eq77)^2));The effective potential illuminates the turning points. Vsq := sqrt((1 - 2*M/r)*(1 + h^2/r^2)):plot([Vsq, a], r=2..40, 0.973..1.12, legend=["effective potential", "energy"]);fsolve(Vsq=a, r, 10..20);fsolve(Vsq=a, r, 20..30);ini := Eq77, Eq78, Eq79;Eq91 := dsolve({Eq54, Eq55, ini}, {r(tau), phi(tau)}, numeric,
output=listprocedure);psch1 := polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=0..2300],
scaling=constrained, axesfont=[TIMES, ROMAN, 12], legend="Schwarzschild"):psch2 := disk([15.47*cos(3.7797), 15.47*sin(3.7797)], 0.5, color=black):psch3 := disk([15.47*cos(11.3399), 15.47*sin(11.3399)], 0.5, color=black):psch4 := disk([15.47*cos(18.9008), 15.47*sin(18.9008)], 0.5, color=black):psch5 := disk([15.47*cos(26.4221), 15.47*sin(26.4221)], 0.5, color=black):pns := disk([0,0],5.8, color=orange):display([psch1, psch2, psch3, psch4, psch5, pns, pnewton]);pschwarz := polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=0..1000],
scaling=constrained, color=blue, legend=Schwarzschild):pschwarz3d := spacecurve([rhs(Eq91(s)[3]), rhs(Eq91(s)[2]),
sqrt(8*(rhs(Eq91(s)[3])-2)), s=0..2500],
coords=cylindrical, color=black, numpoints=400):<Text-field layout="_pstyle10" style="_cstyle11">Kerr metric</Text-field>M := 'M': l := 'l':The same procedure also applies to the Kerr metric; we only change the Lagrangian:a := J/M; Delta := r(tau)^2 - 2*M*r(tau) + a^2;f := 1/2*(-(1 - 2*M/r(tau))*diff(t(tau),tau)^2
+ r(tau)^2/(Delta)*diff(r(tau),tau)^2 + (r(tau)^2 + a^2 + 2*M*a^2/r(tau))*diff(phi(tau),tau)^2-4*a*M/r(tau)*diff(t(tau), tau)*diff(phi(tau), tau));f1 := subs({t(tau)=var1, diff(t(tau),tau)=var2, r(tau)=var3, diff(r(tau),tau)=var4,
phi(tau)=var5, diff(phi(tau),tau)=var6}, f):Employ the Euler-Lagrange equations (Epr11 through Eq34)Epr11 := diff(f1, var2):Epr12 := diff(f1, var1):Epr13 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),
var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr11):Eq14 := Epr13 = -E;Epr21 := diff(f1, var4):Epr22 := diff(f1, var3):Epr23 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),
var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr21):Epr24 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),
var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr22):Epr25 := diff(Epr23, tau):Eq26 := Epr25 - Epr24 = 0;Epr31 := diff(f1, var6):Epr32 := diff(f1, var5):Epr33 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),
var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr31):Eq34 := Epr33 = l;A similar decoupling manipulation (Eq53 through Eq55)Eq53 := isolate(Eq14, diff(t(tau),tau));Eq53p := subs(Eq53, Eq34);Eq54 := isolate(Eq53p, diff(phi(tau),tau));Eq54p := subs(Eq54, Eq53);Eq55 := subs({Eq53, Eq54}, Eq26);Solve the differential equations with initial conditions provided:M := 1; J := 0.37;l := 4.7996; E := 0.9772;Eq77 := r(0) = 26; Eq78 := D(r)(0) = 0;Eq79 := phi(0) = 0;Eq80 := D(phi)(0) = eval(rhs(Eq54), {r(tau)=rhs(Eq77)});ini := Eq77, Eq78, Eq79;Eq91 := dsolve({Eq54, Eq55, ini}, {r(tau), phi(tau)}, numeric,
output=listprocedure);polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=0..2500],
scaling=constrained, color=red, legend="Kerr direct");pkerrd := polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=0..1000],
scaling=constrained, color=red, legend="Kerr direct"):M := 1; J := -0.37;l := 4.7996; E := 0.976797;Eq77 := r(0) = 26; Eq78 := D(r)(0) = 0;Eq79 := phi(0) = 0;Eq80 := D(phi)(0) = eval(rhs(Eq54), {r(tau)=rhs(Eq77)}); ini := Eq77, Eq78, Eq79;Eq91 := dsolve({Eq54, Eq55, ini}, {r(tau), phi(tau)}, numeric,
output=listprocedure);polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=0..2500],
scaling=constrained, legend="Kerr retro");pkerro := polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=0..1000],
scaling=constrained, color=green, legend="Kerr retro", axesfont=[TIMES, ROMAN, 12]):pkerr1 := disk([16.31*cos(3.63467), 16.31*sin(3.63467)], 0.5, color=black):pkerr2 := disk([14.53*cos(3.923), 14.53*sin(3.923)], 0.5, color=black):pkerr3 := disk([16.31*cos(10.9765), 16.31*sin(10.9765)], 0.5, color=black):pkerr4 := disk([14.53*cos(11.7778), 14.53*sin(11.778)], 0.5, color=black):display([pnewton,pschwarz,pkerrd,pkerro,pkerr1,pkerr2,psch2,pkerr3,pkerr4,psch3,pns]);<Text-field layout="_pstyle10" style="_cstyle11">Embedding diagram</Text-field>We derive the embedding formula, see Misner, Thorne, and Wheeler, Gravitation, p. 613 for details. A star is assumed to have uniform density:M := 'M': m := rp^3/R^3*M;zin := int(sqrt(1/(1 - 2*m/rp) -1), rp=2*M..r) + C;zout := int(sqrt(1/(1 - 2*M/rp) -1), rp=2*M..r);Eq201 := eval(zin, r=R) = eval(zout, r=R):C := solve(Eq201, C):pin := plot3d([r*cos(phi), r*sin(phi), eval(zin, {M=1, R=5.8})], phi=0..2*Pi, r=0..5.8, grid=[25,6], orientation=[88,135], shading=zhue):pout := plot3d([r*cos(phi), r*sin(phi), eval(zout, {M=1, R=5.8})], phi=0..2*Pi, r=5.8..27, grid=[25,12]):We replot the trajectory in the section "Schwarzschild metric" in this embedding diagram:display([pin, pout, pschwarz3d], style=hidden);<Text-field layout="_pstyle22" style="_cstyle26">Example 1</Text-field>This example demonstrates an infalling particle with zero angular momentum encountering a Kerr black hole with maximal spin angular momentum NiMlIkpH (NiMvJSJhRyUiTUc=).restart:a := J/M; Delta := r(tau)^2 - 2*M*r(tau) + a^2;f := 1/2*(-(1 - 2*M/r(tau))*diff(t(tau),tau)^2
+ r(tau)^2/(Delta)*diff(r(tau),tau)^2 + (r(tau)^2 + a^2 + 2*M*a^2/r(tau))*diff(phi(tau),tau)^2-4*a*M/r(tau)*diff(t(tau), tau)*diff(phi(tau), tau));f1 := subs({t(tau)=var1, diff(t(tau),tau)=var2, r(tau)=var3, diff(r(tau),tau)=var4,
phi(tau)=var5, diff(phi(tau),tau)=var6}, f):Epr11 := diff(f1, var2):Epr12 := diff(f1, var1):Epr13 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),
var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr11):Eq14 := Epr13 = -E;Epr21 := diff(f1, var4):Epr22 := diff(f1, var3):Epr23 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),
var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr21):Epr24 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),
var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr22):Epr25 := diff(Epr23, tau):Eq26 := Epr25 - Epr24 = 0;Epr31 := diff(f1, var6):Epr32 := diff(f1, var5):Epr33 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),
var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr31):Eq34 := Epr33 = l;Eq53 := isolate(Eq14, diff(t(tau),tau));Eq53p := subs(Eq53, Eq34);Eq54 := isolate(Eq53p, diff(phi(tau),tau));Eq54p := subs(Eq54, Eq53);Eq55 := subs({Eq53, Eq54}, Eq26);M := 1; J := 1;l := 0; E := 1;Eq77 := r(0) = 5; Eq78 := D(r)(0) = -0.64;Eq79 := phi(0) = 0;Eq80 := D(phi)(0) = eval(rhs(Eq54), {r(tau)=rhs(Eq77)});ini := Eq77, Eq78, Eq79;Eq91 := dsolve({Eq54, Eq55, ini}, {r(tau), phi(tau)}, numeric,
output=listprocedure);with(plots):polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=0..4.49],
scaling=constrained, numpoints=2400);<Text-field layout="_pstyle22" style="_cstyle26">Example 2</Text-field>This example demonstrates the behavior of a particle near the unstable circular orbit. The particle orbits several times, then flies out. An animation of the motion is included. restart:f := 1/2*(-(c^2 - 2*G*M/r(s))*diff(t(s),s)^2 + 1/(1 - 2*G*M/(c^2*r(s)))*diff(r(s),s)^2 + r(s)^2*diff(phi(s),s)^2);f1 := subs({t(s)=var1, diff(t(s),s)=var2, r(s)=var3, diff(r(s),s)=var4, phi(s)=var5, diff(phi(s),s)=var6}, f):Epr11 := diff(f1, var2):Epr12 := diff(f1, var1):Epr13 := subs({var1=t(s), var2=diff(t(s),s), var3=r(s), var5=phi(s), var4=diff(r(s),s), var6=diff(phi(s),s)}, Epr11):Eq14 := Epr13 = -a/c;Epr21 := diff(f1, var4):Epr22 := diff(f1, var3):Epr23 := subs({var1=t(s), var2=diff(t(s),s), var3=r(s), var5=phi(s), var4=diff(r(s),s), var6=diff(phi(s),s)}, Epr21):Epr24 := subs({var1=t(s), var2=diff(t(s),s), var3=r(s), var5=phi(s), var4=diff(r(s),s), var6=diff(phi(s),s)}, Epr22):Epr25 := diff(Epr23, s):Eq26 := Epr25 - Epr24 = 0;Epr31 := diff(f1, var6):Epr32 := diff(f1, var5):Epr33 := subs({var1=t(s), var2=diff(t(s),s), var3=r(s), var5=phi(s), var4=diff(r(s),s), var6=diff(phi(s),s)}, Epr31):Eq34 := Epr33 = h/c;Eq53 := isolate(Eq14, diff(t(s),s));Eq54 := isolate(Eq34, diff(phi(s),s));Eq55 := eval(Eq26, {Eq53, Eq54});M := 1; G := 1; c := 1;Eq77 := r(0) = 3.76777; Eq78 := D(r)(0) = 0;Eq79 := phi(0) = 0;Eq80 := D(phi)(0) = 0.30290054; h := eval(lhs(Eq34), {r(s)=rhs(Eq77), diff(phi(s),s)=rhs(Eq80)}); a := sqrt((1 - 2*M/rhs(Eq77))*(1 + h^2/rhs(Eq77)^2));Vsq := sqrt((1 - 2*M/x)*(1 + h^2/x^2)):plot([Vsq, a], x=2..25, a-0.1*a..a+0.1*a, legend=["effective potential", "energy"]);ini := Eq77, Eq78, Eq79;Eq91 := dsolve({Eq54, Eq55, ini}, {r(s), phi(s)}, numeric, output=listprocedure);with(plots):polarplot([rhs(Eq91(s)[3]), rhs(Eq91(s)[2]), s=0..100], scaling=constrained);p1 := plot3d([x*cos(t), x*sin(t), sqrt(8*(x-2))], x=2..14, t=0..2*Pi, style=hidden): p2 := spacecurve([rhs(Eq91(s)[3]), rhs(Eq91(s)[2]), sqrt(8*(rhs(Eq91(s)[3])-2)), s=0..100], coords=cylindrical, color=black, numpoints=400):display([p1, p2]);N := 50; nstep := 2;for i from 0 by nstep to N*nstep do
pl[i] := polarplot([rhs(Eq91(s)[3]), rhs(Eq91(s)[2]), s=i..i+nstep], numpoints=5):
pl3d[i] := spacecurve([rhs(Eq91(s)[3]), rhs(Eq91(s)[2]), sqrt(8*(rhs(Eq91(s)[3])-2)), s=i..i+nstep], coords=cylindrical, numpoints=5, color=black):
plsum[i] := display(seq(pl[j*nstep], j=0..i/nstep)):
pl3dsum[i] := display({p1, seq(pl3d[j*nstep], j=0..i/nstep)}):
end do:display([seq(plsum[i*nstep], i=0..N)], insequence=true, scaling=constrained);display([seq(pl3dsum[i*nstep], i=0..N)], insequence=true);<Text-field layout="_pstyle22" style="_cstyle26">Example 3</Text-field>Literature on numerical integrations for relativistic orbits can be found in the early 1970's; see M. Johnston and R. Ruffini, "Generalized Wilkins effect and selected orbits in a Kerr-Newman geometry," Phys. Rev. D 10, 2324-2329 (1974). This example uses the same initial conditions for Fig. 6 of that paper (a particle initially orbits against a maximally rotating black hole). Our calculation shows interesting details of the last orbits that was omitted in that paper perhaps due to limited computational power of that era. restart:a := J/M; Delta := r(tau)^2 - 2*M*r(tau) + a^2;f := 1/2*(-(1 - 2*M/r(tau))*diff(t(tau),tau)^2
+ r(tau)^2/(Delta)*diff(r(tau),tau)^2 + (r(tau)^2 + a^2 + 2*M*a^2/r(tau))*diff(phi(tau),tau)^2-4*a*M/r(tau)*diff(t(tau), tau)*diff(phi(tau), tau));f1 := subs({t(tau)=var1, diff(t(tau),tau)=var2, r(tau)=var3, diff(r(tau),tau)=var4,
phi(tau)=var5, diff(phi(tau),tau)=var6}, f):Epr11 := diff(f1, var2):Epr12 := diff(f1, var1):Epr13 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),
var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr11):Eq14 := Epr13 = -E;Epr21 := diff(f1, var4):Epr22 := diff(f1, var3):Epr23 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),
var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr21):Epr24 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),
var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr22):Epr25 := diff(Epr23, tau):Eq26 := Epr25 - Epr24 = 0;Epr31 := diff(f1, var6):Epr32 := diff(f1, var5):Epr33 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau),
var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, Epr31):Eq34 := Epr33 = l;Eq53 := isolate(Eq14, diff(t(tau),tau));Eq53p := subs(Eq53, Eq34);Eq54 := isolate(Eq53p, diff(phi(tau),tau));Eq54p := subs(Eq54, Eq53);Eq55 := subs({Eq53, Eq54}, Eq26);M := 1; J := -1;l := 4; E := 0.97;A := (r^3 + a^2*r + 2*M*a^2);B := -4*a*M*l;C := -(r - 2*M)*l^2 - r*(r^2 - 2*M*r + a^2);Veff := (-B + sqrt(B^2 - 4*A*C))/(2*A);plot({Veff, E}, r=1..30, 0..1.2);fsolve(Veff=E, r=20..30);Eq77 := r(0) = %; Eq78 := D(r)(0) = 0;Eq79 := phi(0) = -5*Pi/4;Eq80 := D(phi)(0) = eval(rhs(Eq54), {r(tau)=rhs(Eq77)});ini := Eq77, Eq78, Eq79;Eq91 := dsolve({Eq54, Eq55, ini}, {r(tau), phi(tau)}, numeric,
output=listprocedure);with(plots): with(plottools):pbh := disk([0,0], 1, color=black):pergo := polarplot(2, phi=0..2*Pi, color=black, thickness=2):porb1 := polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=-180..207.4], scaling=constrained, axes=boxed):This plot resembles Fig. 6 of Johnston and Ruffini: the particle initially orbits counterclockwise, and the black hole spins clockwise. There is an indication of reversal of direction when the particle approaches the black hole, but we can see more details in the following. display([pbh, porb1]);pergo3d := sphereplot(M + sqrt(M^2 - a^2*cos(theta)^2), phi=0..2*Pi, theta=0..Pi, scaling=constrained, style=wireframe, grid=[16,16]):pkerrbh := sphere([0,0,0], 1, color=black):This plot shows the ergosphere of a Kerr black hole: the region between the static limit described by NiMvJSJyRywmJSJNRyIiIi0lJXNxcnRHNiMsJiokRiYiIiNGJyomJSJhR0YtLSUkY29zRzYjJSZ0aGV0YUdGLSEiIkYn and the event horizon described by NiMvJSJyRyUiTUc=. Within the ergosphere, a particle (even a photon) must corotate with the black hole. display([pergo3d, pkerrbh]);porb3d1 := spacecurve([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), 0, tau=200..207.4], coords=cylindrical, scaling=constrained, color=red, thickness=2):porb3d2 := spacecurve([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), 0, tau=207.4..207.648], coords=cylindrical, scaling=constrained, numpoints=800, color=red, thickness=2):porb2 := polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=200..207.4], scaling=constrained, axes=boxed):porb3 := polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=207.4..207.648], scaling=constrained, numpoints=400, axes=boxed):The numerical solution shows that the particle reverses its direction and corotates with the black hole as it spirals toward it. display([pkerrbh, pergo3d, porb3d1, porb3d2]);Viewing from the top:display([pbh, pergo, porb2, porb3]);