Code C4: Schwarzschild inspirals
Last edited by Jens Boos (jens.boos@kit.edu), Oct 7, 2024.
(For definitions and conventions refer to the lecture notes.)
Website: http://www.spintwo.net/Courses/EMRI-101/

(* Overall units: GM = 1 *)

(* mass ratio (make it larger here to increase the effect) *)
μ=10^-2;

(* particle energy and angular momentum *)
EE=0.97;
LLmin=Sqrt[12];
LL=LLmin+0.4;

(* bound orbit: this quantity has to be negative *)
check=EE^2-2(4-LL(LL-Sqrt[LL^2-12]))^2/(LL(LL-Sqrt[LL^2-12])^3)

(* simulation duration *)
τmax = 150000;

(* initial radial position (in units of GM=1) *)
r0=20;

(* convert these into initial 4-momentum *)
tdot0=EE/(1-2/r0);
rdot0=Sqrt[EE^2-(1+LL^2/r0^2)(1-2/r0)];
ϕdot0=LL/r0^2;

(* SANITY CHECK: verify that initial conditions describe timelike particle *)
test=(-1)(1-2/r0)tdot0^2+rdot0^2/(1-2/r0)+r0^2ϕdot0^2//N

(* basic metric function f = f(r) and its derivative fp = f'(r) *)
f:=1-2/r[τ];
fp:=2/r[τ]^2;

(* Christoffel symbols, index position: Γttr = (Subscript[((OverTilde[Γ])^t), tr]) etc *)
Γttr:=fp/(2f);
Γrtt:=f fp/2;
Γrrr:=-fp/(2f);
Γrθθ:=-f r[τ];
Γrϕϕ:=-f r[τ]Sin[θ[τ]]^2;
Γθrθ:=1/r[τ];
Γϕrϕ:=1/r[τ];
Γθϕϕ:=-Sin[θ[τ]]Cos[θ[τ]];
Γϕθϕ:=Cos[θ[τ]]/Sin[θ[τ]];
θ[τ_]:=π/2;

(* quadrupole formula *)
II = μ r[τ]^2{{Cos[ϕ[τ]]^2-1/3,Cos[ϕ[τ]]Sin[ϕ[τ]],0},{Cos[ϕ[τ]]Sin[ϕ[τ]],Sin[ϕ[τ]]^2-1/3,0},{0,0,-1/3}};
IIdd=FullSimplify[1/t'[τ] D[1/t'[τ] D[II,τ],τ]];
IIddd=FullSimplify[1/t'[τ] D[1/t'[τ] D[1/t'[τ] D[II,τ],τ],τ]];
EEdot=(1/5)FullSimplify[Tr[IIddd.IIddd]];
LLdot=(2/5)FullSimplify[(IIdd.IIddd)[[1,2]]-(IIdd.IIddd)[[2,1]]];

(* integrate EOMs *)
sol=NDSolve[
{t''[τ]+2Γttr EEE[τ]/f r'[τ]==0,
r''[τ]+Γrtt EEE[τ]^2/f^2+Γrrr r'[τ]^2+Γrϕϕ LLL[τ]^2/r[τ]^4==0,
ϕ''[τ]+2Γϕrϕ r'[τ] LLL[τ]/r[τ]^2==0,
EEE'[τ]==-t'[τ]EEdot,
LLL'[τ]==-t'[τ]LLdot,
t[0]==0,r[0]==r0,ϕ[0]==0,
EEE[0]==EE,LLL[0]==LL,
t'[0]==tdot0,r'[0]==rdot0,ϕ'[0]==ϕdot0
},{t,r,ϕ,EEE,LLL},{τ,0,τmax},PrecisionGoal->25
][[1]];
rmax=1.5r0;
curve[ττ_]:={r[τ]Cos[ϕ[τ]],r[τ]Sin[ϕ[τ]]}/.sol/.τ->ττ;
rpos[ττ_]:=r[τ]/.sol/.τ->ττ;
Veff[r_,τ_]:=(1-2/r)(1+LLL[τ]^2/r^2)/.sol;
p[ττmax_]:=ParametricPlot[curve[τ],{τ,0,ττmax},AspectRatio->1,PlotStyle->{Opacity[0.5]},PlotRange->{{-rmax,rmax},{-rmax,rmax}},Frame->True,GridLines->Automatic,FrameLabel->{"x/GM","y/GM"},PlotPoints->1000];
v[τ_]:=Show[
Plot[{Veff[r,τ],EEE[τ]^2/.sol},{r,0,rmax},AxesLabel->{"r/GM"},AspectRatio->1,PlotRange->{{0,rmax},{0.9,1}},PlotLegends->Placed[LineLegend[{"(Subscript[V, eff])","(E^2)"}],{Right,Top}]],
Graphics[{Red,Circle[{rpos[τ],Veff[rpos[τ],τ]},{rmax/60,rmax/60/270}]}]
];
movie[τ0_]:=
GraphicsRow[{
Show[p[τ0],Graphics[{Disk[{0,0},1]}],Graphics[{Red,Disk[curve[τ0],0.75 (rmax/60)]}]],
v[τ0]},AspectRatio->1,ImageSize->800];
Animate[movie[τ],{τ,0.01,τmax},AnimationRate->50]

(* total number of orbits *)
ϕ[τmax]/(2π)/.sol

(* it is useful to convert proper time into coordinate time *)
tt=t/.sol[[1]];
ττ=InverseFunction[tt];
Show[
Plot[tt[τ],{τ,0,0.95τmax},PlotRange->All,PlotStyle->ColorData[97,"ColorList"][[1]]],
Plot[ττ[t],{t,tt[0],tt[0.95τmax]},PlotRange->All,PlotStyle->ColorData[97,"ColorList"][[2]]],
AxesOrigin->{0,0},AspectRatio->0.625]