Spin - orbital coupling and motion demonstration

In this notebook we will demonstrate some properties of the equations of motion of spinning binaries. First, the 1.5PN equations of motion are examined and solved. Second, we plot some Kerr geodesics.

This notebook was prepared for the Selected chapters of GR lecture in October 2023 by Vojtěch Witzany.

1.5 PN Hamiltonian

Vectors in Hamiltonian

xv = {r Cos[ϕ] Sin[θ], r Sin[ϕ] Sin[θ],r Cos[θ]};
pv = {pr,((pϕ)/((r^2) ((Sin[θ])^2))) ,((pθ)/(r^2)) }.{D[xv,r],D[xv,ϕ],D[xv,θ]}//Simplify
Lv = Cross[xv,pv]
S1v ={(Sqrt[(S1^2)-((ξ1)^2)])Sin[ψ1],(Sqrt[(S1^2)-((ξ1)^2)])Cos[ψ1],ξ1}
S2v ={(Sqrt[(S2^2)-((ξ2)^2)])Sin[ψ2],(Sqrt[(S2^2)-((ξ2)^2)])Cos[ψ2],ξ2}
Seffv = (1 + (3/4)(m2/m1))S1v + (1 + (3/4)(m1/m2))S2v;

Hamiltonian

H0 = μ(((pv.pv)/2) - (M/r));
H1 = μ(((3ν-1)/8)(((pv.pv))^2) + ((M^2)/(2(r^2)))-(M/(2r))((3+ν)pv.pv+ν (pr^2)));
H15 = (2/(r^3))Lv .Seffv;
H = Simplify[(H0+H1+H15)/.{μ->1,M->1}]

Derivation of equations of motion

momenta = {pr,pϕ,pθ,ξ1,ξ2};
coords = {r,ϕ,θ,ψ1,ψ2};
momEOM = Table[-D[H,coords[[i]]],{i,1,5}]//Simplify;
coordEOM = Table[D[H,momenta[[i]]],{i,1,5}]//Simplify;
funSub = Join[Table[coords[[i]]->coords[[i]][t],{i,1,5}] ,Table[momenta[[i]]->momenta[[i]][t],{i,1,5}] ];
RHS = Join[momEOM,coordEOM]/.funSub;
LHS = Join[ Table[momenta[[i]]'[t],{i,1,5}],Table[coords[[i]]'[t],{i,1,5}]];
EOM = Table[LHS[[i]]==RHS[[i]],{i,1,10}]

Parameters, initiation

(*Units so that total m is 1*)
q = 1/100;
χ1 = 1;
χ2 =1; 
(*Mass ratio q, χs equal to 1 correspond to maximally spinning BHs. q smaller than 1 means m1 is heavier*)
(*Spins are angular momenta, which scale with the respective mass squared (one for size of object, another one for the momentum they carry)*)
parSub = {m1 -> (1/(1+q)),m2->(q/(1+q)),ν-> (q/(((1+q))^2)),S1->((( (1/(1+q))))^2)χ1,S2->((((q/(1+q))))^2)χ2}
(*Note that xi1 and xi2 have to be bounded by S1 S2*)
r0 = 30;
Init = {
r[0]==r0,
ϕ[0]==0,
θ[0]==π/2,
pr[0] == 0,
pϕ[0] == 0.9(Sqrt[r0]),
pθ[0]==0.2(Sqrt[r0]),
ξ1[0] ==0.99 S1,
ξ2[0] ==0.2 S2,
ψ1[0] == 0,
ψ2[0] == π
}

Integration

traj = NDSolve[Join[EOM,Init]/.parSub,Join[Table[coords[[i]][t],{i,1,5}],Table[momenta[[i]][t],{i,1,5}]],{t,0,100 (r0^(3/2))}]

ParametricPlot3D[xv/.funSub/.traj, {t,0,100 (r0^(3/2))}]
Show[Graphics3D[Sphere[{0,0,0}]],ParametricPlot3D[(S2v/S2)/.parSub/.funSub/.traj, {t,0,100 (r0^(3/2))}]]
Show[Graphics3D[Sphere[{0,0,0}]],ParametricPlot3D[(S1v/S1)/.parSub/.funSub/.traj, {t,0,100 (r0^(3/2))}]]
Show[Graphics3D[Sphere[{0,0,0}]],ParametricPlot3D[(Lv/(Sqrt[Lv.Lv]))/.parSub/.funSub/.traj, {t,0,100 (r0^(3/2))}]]

(Lv/(Sqrt[Lv.Lv]))/.parSub/.funSub/.traj/.{t->1}

Kerr geodesics

Test particle Hamiltonian

Σ = (r^2) + (a^2) ((Cos [θ])^2);
Δ = (r^2) + (a^2) - 2M r;
gutt = ((-Σ ((r^2)+(a^2)) - 2M r (a^2) ((Sin [θ])^2))/(Δ Σ));
gutϕ = ((-2M r a)/(Δ Σ));
guϕϕ= ((Δ- (a^2) ((Sin[θ])^2))/(Δ Σ ((Sin[θ])^2)));
gurr = (Δ/Σ) ;
guθθ = (1/Σ) ;
HK = ((gutϕ)/gutt) pϕ + (Sqrt[(((((gutϕ)/gutt) pϕ))^2) - (1/gutt)(gurr (pr^2)+guϕϕ ((pϕ)^2)+guθθ ((pθ)^2)+1)])
(*We could integrate the EOM, but instead we will use the Kerrgeodesics implementation here.*)

KerrGeodesics implementation

<< KerrGeodesics`

(*The KerrGeodesics package implements analytical solutions for Kerr geodesics, for more see bhptoolkit.org and https://github.com/BlackHolePerturbationToolkit/KerrGeodesics.
Install as follows:
 PacletSiteRegister["https://pacletserver.bhptoolkit.org","Black Hole Perturbation Toolkit Paclet Server"]
PacletSiteUpdate["https://pacletserver.bhptoolkit.org"]
PacletFindRemote["*",<|"Location"->"https://pacletserver.bhptoolkit.org"|>]
PacletInstall["KerrGeodesics"];
*)

(*
KerrGeoOrbit generates a bound geodesic parametrized by its turning points. It takes the arguments:
 - spin a, 
  - semi-latus rectum p, 
  - eccentricity e, 
  - and the parameter x ∈ (-1,1) such that i = arcsin((Sqrt[1-(x^2)])) is the inclination angle with respect to the equatorial plane. (negative values of x correspond to counter-rotating orbits).
*)
orbit=KerrGeoOrbit[0.99,15,0.4,0.7]
{t,r,θ,φ}=orbit["Trajectory"];
Show[ParametricPlot3D[{r[λ] Sin[θ[λ]] Cos[φ[λ]],r[λ] Sin[θ[λ]] Sin[φ[λ]],r[λ] Cos[θ[λ]]},{λ,0,20},Boxed->False,Axes->False,PlotStyle->Red,PlotRange->All],Graphics3D[{Black,Sphere[{0,0,0},1+Sqrt[1-0.99^2]]}]]
Show[ParametricPlot[{r[λ] Sin[θ[λ]], r[λ] Cos[θ[λ]]},{λ,0,20},PlotStyle->Red,PlotRange->All],Graphics[{Black,Disk[{0,0},1+Sqrt[1-0.99^2]]}]]
Show[ParametricPlot[{r[λ] Sin[φ[λ]], r[λ] Cos[φ[λ]]},{λ,0,20},PlotStyle->Red,PlotRange->All],Graphics[{Black,Disk[{0,0},1+Sqrt[1-0.99^2]]}]]
Show[ParametricPlot[{r[λ],  Cos[θ[λ]]},{λ,0,20},PlotStyle->Red,AxesLabel->{"r","cos(θ)"},AspectRatio->1]]
Show[ParametricPlot3D[{Sin[φ[λ]]Cos[θ[λ]],Cos[φ[λ]]Cos[θ[λ]],Sin[θ[λ]]},{λ,0,20},PlotStyle->Red,PlotRange->All,AxesLabel->{"x","y","sin(θ)"}],Graphics3D[Sphere[{0,0,0}]]]
ClearAll[t,r,θ,φ]