Introduction to (bound) geodesics in Kerr spacetime and the KerrGeodesic` package

Maarten van de Meent

BHPToolkit Spring 2020 workshop

Geodesics in Kerr spacetime

The Kerr metric

Kerr metric (G = M = c = 1):

(Subscript[g, μν])=(({{-1+((2 r)/Σ), 0, 0, ((2 a r (-1+(z^2)))/Σ)}, {0, (Σ/Δ), 0, 0}, {0, 0, (Σ/(1-(z^2))), 0}, {((2 a r (-1+(z^2)))/Σ), 0, 0, (((-1+(z^2)) (2 (a^2) r (-1+(z^2))-Σ((a^2)+(r^2)) ))/Σ)}}))

Σ == (r^2)+(a^2)(z^2)
Δ== (r-2)r+(a^2)

Geodesics

Trajectory :

(x^μ)==(t,r,z,ϕ)
(u^μ)== (d/(dτ))(x^μ)

We want (x^μ) to solve the geodesic equation:

(u^α)(Subscript[∇, α])(u^μ)=  (((d^2)(x^μ))/((dτ)^2))-(Subscript[Γ, αβ]^μ)((dx^α)/(dτ))((dx^β)/(dτ))==0

Constants of Motion

Norm

Normalization of 4 - velocity

-1 == (u^μ)(Subscript[u, μ])

Energy and Angular Momentum

Kerr spacetime has two Killing vector fields:

(kt^μ):=-(∂/(∂t))     and     ((kϕ)^μ):=(∂/(∂ϕ))

Two more conserved quantities, the (specific) Energy and Angular Momentum :

ℰ ==  (Subscript[kt, μ])(u^μ)

ℒ ==(Subscript[kϕ, μ])(u^μ)

Carter Constant

Killing tensor :

(Subscript[K, μν]):= 2Σ (Subscript[l, (μ])(Subscript[n, ν)])+(r^2)(Subscript[g, μν])

A fourth conserved quantity, the Carter Constant :

𝒬 +(((ℒ-a ℰ))^2) == (u^μ)(Subscript[K, μν])(u^ν)

Equations of Motion

Solve (1)-(4):

(((Σ((d r)/(dτ))))^2)== -(1-(ℰ^2))(r^4)+2 (r^3) -((a^2) (1-(ℰ^2))+(ℒ^2)+𝒬) (r^2)+2((a^2) (ℰ^2)-2 a ℰ ℒ+(ℒ^2)+𝒬) r -(a^2)𝒬
(((Σ((d z)/(dτ))))^2)==(a^2) (1-(ℰ^2)) (z^4)-((a^2) (1-(ℰ^2))+(ℒ^2)+𝒬) (z^2)+𝒬
Σ((d t)/(dτ)) == ((r (2 (a^2)  ℰ-2 a  ℒ+(a^2) ℰ r+ℰ (r^3)))/Δ)+(a^2) ℰ (z^2)
Σ((d ϕ)/(dτ)) ==((a (2  ℰ r-a ℒ))/Δ)+(ℒ/(1-(z^2)))

Introduce "Mino time"

dλ ==  ((dτ)/Σ)

(((((d r)/(dλ))))^2)== (-1+(ℰ^2))(r^4)+2 (r^3) +((a^2) (-1+(ℰ^2))-(ℒ^2)-𝒬) (r^2)+2((a^2) (ℰ^2)-2 a ℰ ℒ+(ℒ^2)+𝒬) r -(a^2)𝒬 ==R(r)
(((((d z)/(dλ))))^2)==-(a^2) (-1+(ℰ^2)) (z^4)+((a^2) (-1+(ℰ^2))-(ℒ^2)-𝒬) (z^2)+𝒬== Z(z)
((d t)/(dλ)) == ((r (2 (a^2)  ℰ-2 a  ℒ+(a^2) ℰ r+ℰ (r^3)))/Δ)+(a^2) ℰ (z^2)
((d ϕ)/(dλ)) ==((a (2  ℰ r-a ℒ))/Δ)+(ℒ/(1-(z^2)))

Turning points :

Radial equation

R(r)== (-1+(ℰ^2))(r-(Subscript[r, 1]))(r-(Subscript[r, 2]))(r-(Subscript[r, 3]))(r-(Subscript[r, 4]))

(for bound orbits)
(Subscript[r, 1])≥ (Subscript[r, 2])≥ (Subscript[r, 3])>(SubPlus[r])>(SubMinus[r])>(Subscript[r, 4])≥ 0

With[
{apex=Sequence[.9,3.5,.4,.95]},
Module[
{r1,r2,r3,r4,ℰ,a},
a=First@List[apex];
{r1,r2,r3,r4}=KerrGeodesics`OrbitalFrequencies`Private`KerrGeoRadialRoots[apex];
ℰ=KerrGeoEnergy[apex];
Plot[(-1+(ℰ^2))(r-r1)(r-r2)(r-r3)(r-r4),{r,0,6},PlotRange->{-3,7.5},PlotTheme->"Detailed",Axes->True,PlotLegends->None,Epilog->{Red,Point[{{r1,0},{r2,0},{r3,0},{r4,0}}]},
FrameLabel->{"r","R[r]"},
ImageSize->500,
FrameTicks->{{Automatic,Automatic},{{0,1,2,3,4,5,6,7,8},{{r1,"(Subscript[r, 1])"},{r2,"(Subscript[r, 2])"},{r3,"(Subscript[r, 3])"},{r4,"(Subscript[r, 4])"},{1+(Sqrt[1-(a^2)]),"(SubPlus[r])"},{1-(Sqrt[1-(a^2)]),"(SubMinus[r])"}}}}
]
]
]

Polar equation

Z(z)==((z^2)-(Subscript[z, 1]^2))((a^2)(-1+(ℰ^2))(z^2)-(Subscript[z, 2]^2))

(for bound orbits)
0≤ (Subscript[z, 1])≤ 1≤ (Subscript[z, 2])

With[
{apex=Sequence[.999,3.5,.9,.25]},
Module[
{z1,z2,ℰ,a},
a=First@List[apex];
{z2,z1}=KerrGeodesics`OrbitalFrequencies`Private`KerrGeoPolarRoots[apex];
ℰ=KerrGeoEnergy[apex];
Plot[((z^2)-(z1^2))((a^2)(1-(ℰ^2))(z^2)-(z2^2)),{z,-(z2/(Sqrt[(a^2)(1-(ℰ^2))])),(z2/(Sqrt[(a^2)(1-(ℰ^2))]))},PlotRange->All,PlotTheme->"Detailed",Axes->True,PlotLegends->None,Epilog->{Red,Point[{{z1,0},{(z2/(Sqrt[(a^2)(1-(ℰ^2))])),0},{-z1,0},{-(z2/(Sqrt[(a^2)(1-(ℰ^2))])),0}}]},
FrameLabel->{"z","Z[z]"},
ImageSize->500,
FrameTicks->{{Automatic,Automatic},{Range[-15,15],{{z1,"(Subscript[z, 1])"},{(z2/(Sqrt[(a^2)(1-(ℰ^2))])),"(Subscript[z, 2])"},{-z1,"-(Subscript[z, 1])"},{-(z2/(Sqrt[(a^2)(1-(ℰ^2))])),"-(Subscript[z, 2])"}}}}
]
]
]

"apex" parameters

The turning points (Subscript[r, 1]), (Subscript[r, 2]), and (Subscript[z, 1])uniquely identify a bound orbit. In the KerrGeodesics` package they expressed in terms of the semi-latus rectum p, eccentricity e, and inclination parameter x.

(Subscript[r, 1])== (p/(1-e))
(Subscript[r, 2])== (p/(1+e))
(Subscript[z, 1])==(Sqrt[1-(x^2)])

or

p==2(((Subscript[r, 1])(Subscript[r, 2]))/((Subscript[r, 1])+(Subscript[r, 2])))
e== (((Subscript[r, 1])-(Subscript[r, 2]))/((Subscript[r, 1])+(Subscript[r, 2])))
x ==  Sign[ℒ](Sqrt[1-(Subscript[z, 1]^2)])

(circular) 0<e<1 (parabolic)
(retrograde) -1 < x < 1 (prograde)
p> (Subscript[p, LSO])(e,x)

KerrGeodesics` functions

<<KerrGeodesics`

Constants of Motion

Analytic

KerrGeoEnergy[a,p,e,x]

Numerical Evaluation:

With[
{
a=.9,p=10.,e=.5,x=.5
},
{
KerrGeoEnergy[a,p,e,x],
KerrGeoAngularMomentum[a,p,e,x],
KerrGeoCarterConstant[a,p,e,x]
}
]

With[
{
a=.9,p=10.,e=.5,x=.5
},
KerrGeoConstantsOfMotion[a,p,e,x]
]

Numerical Evaluation follows the precision of the input.

With[
{
a=.9`100,p=10.`100,e=.5`100,x=.5`100
},
KerrGeoConstantsOfMotion[a,p,e,x]
]

Frequencies

apex=Sequence[.9,10.,.5,.5];

KerrGeoFrequencies[apex,"Time"-> "Mino"]
KerrGeoFrequencies[apex,"Time"-> "Proper"]
KerrGeoFrequencies[apex,"Time"-> "BoyerLindquist"]

Orbits

Analytic solutions

The equations of motion can be solve analytically.

r== (((Subscript[r, 3])((Subscript[r, 1]) - (Subscript[r, 2]))((JacobiSN[EllipticK[(Subscript[k, r])]/π (Subscript[q, r]),(Subscript[k, r])])^2)-(Subscript[r, 2])((Subscript[r, 1])-(Subscript[r, 3])))/(((Subscript[r, 1])-(Subscript[r, 2]))((JacobiSN[EllipticK[(Subscript[k, r])]/π (Subscript[q, r]),(Subscript[k, r])])^2)-((Subscript[r, 1])-(Subscript[r, 3]))))
z== (Subscript[z, 1]) JacobiSN[EllipticK[(Subscript[k, z])] 2/π (Subscript[q, z]),(Subscript[k, z])]

with

(Subscript[k, r])== ((((Subscript[r, 1])-(Subscript[r, 2]))((Subscript[r, 3])-(Subscript[r, 4])))/(((Subscript[r, 1])-(Subscript[r, 3]))((Subscript[r, 2])-(Subscript[r, 4]))))
(Subscript[k, z])==(a^2)(1-(ℰ^2)) ((Subscript[z, 1]^2)/(Subscript[z, 2]^2))

Phases

(Subscript[q, r])== (Subscript[Υ, r])λ+(Subscript[q, r,0])
(Subscript[q, z])== (Subscript[Υ, z])λ+(Subscript[q, z,0])

(Subscript[q, r]) ==  0 -> r==(Subscript[r, 2])
(Subscript[q, z])== 0 -> z==0, z'>0

Evaluation

With[
{
a=.9,p=10.,e=.5,x=.5
},
orbit=KerrGeoOrbit[a,p,e,x,Method->"Analytic"]
]

Evaluate orbit

orbit[λ]

orbit[2.]

Extract constants of motion etc.

orbit["ConstantsOfMotion"]

orbit["Frequencies"]

All Keys

Keys@orbit[[-1]]

Fast spectral implementation

With[
{
a=.9,p=10.,e=.5,x=.5
},
fastorbit=KerrGeoOrbit[a,p,e,x,Method->"FastSpec"]
];//Timing

orbit[10.]//RepeatedTiming
fastorbit[10.]//RepeatedTiming

((orbit[10.])/(fastorbit[10.]))-1

Plot

Module[
{to,ro,θo,ϕo,λ,Λr},
{to,ro,θo,ϕo}=fastorbit[λ];
Λr=2π/orbit["RadialFrequency"];
Show[
ParametricPlot3D[
ro{Cos[ϕo]Sin[θo],Sin[ϕo]Sin[θo],Cos[θo]},{λ,0,1Λr}
,Lighting->"Neutral",PlotPoints->20,PlotRange->20],
Graphics3D[{Black,Specularity[.5],Sphere[{0,0,0},1+(Sqrt[1-((.9)^2)])]}]
]
]

Special Orbits

Circular Orbits

e ==0 
(Subscript[r, 1])==(Subscript[r, 2])

Innermost bound spherical orbit (IBSO)

(Subscript[r, 1])==(Subscript[r, 2]) && ℰ==1

Analytic for polar orbits

KerrGeoIBSO[a,0]

Generic

KerrGeoIBSO[.9,.5]

Photon sphere

ℰ==∞
ℒ == ∞
e==0

Schwarzschild

KerrGeoPhotonSphereRadius[0,x]

Kerr equatorial

KerrGeoPhotonSphereRadius[a,1]
KerrGeoPhotonSphereRadius[a,-1]

Generic

ContourPlot[
KerrGeoPhotonSphereRadius[a,x],{x,-1,1},{a,-1,1},PlotTheme->"Detailed",FrameLabel->{"x","a"}]

Parabolic or Marginally bound orbits

ℰ == 1
e == 1
(Subscript[r, 1])== ∞

Currently not implemented :

KerrGeoEnergy[a,p,1,x]

Equatorial orbits

x==±1
(Subscript[z, 1])==0
𝒬 ==0

Polar orbits

x==0
(Subscript[z, 1])==1
ℒ==0

Constants of motion work :

KerrGeoConstantsOfMotion[.9,10,.1,0]

Frequencies are broken

KerrGeoFrequencies[.9,10,.1,0,"Time"->"Mino"]

Last Stable orbit or Separatrix

(Subscript[r, 2])== (Subscript[r, 3])

KerrGeoSeparatrix[.9,.5,.5]

Innermost stable circular (equatorial) orbit (ISCO)

(Subscript[r, 1])==(Subscript[r, 2])== (Subscript[r, 3]) && (Subscript[z, 1])==0

Analytic :

KerrGeoISCO[a,1]
KerrGeoISCO[a,-1]

Plot[{KerrGeoISCO[a,1],KerrGeoISCO[a,-1]},{a,-1,1},PlotTheme->"Detailed",ImageSize->700]

Innermost stable spherical  orbit (ISSO)

(Subscript[r, 1])==(Subscript[r, 2])== (Subscript[r, 3])

KerrGeoISSO[.9,.5]

ContourPlot[
KerrGeoISSO[a,x],{x,-1,1},{a,(10^(-3)),1-(10^(-3))},PlotTheme->"Detailed",FrameLabel->{"x","a"}]

Generic

KerrGeoSeparatrix[(*a*).9,(*e*).5,(*x*).5]

Parallel Transport

Parallel transport equation :

(u^α)(Subscript[∇, α])(Subscript[((Subscript[e, a])), μ])==0

The 4 - velocity (Subscript[u, α]) itself is parallel transported, so can use it as a first leg of a tetrad.

(Subscript[((Subscript[e, 1])), μ])==(Subscript[u, μ])

We can use the Killing - Yano tensor (Subscript[ℱ, μν]) to construct another parallel propogated tetrad leg.

(Subscript[((Subscript[e, 4])), μ])==(1/𝒦)(Subscript[ℱ, μα])(u^α)

With two legs fixed the remain degree of freedom of the tetrad is a rotation of the plane orthogonal to the first two legs. The corresponding angle satisfies a simple first order equation of motion.

((dψ)/(dλ))== (Sqrt[𝒦])(((((r^2)+(a^2))ℰ-a ℒ)/(𝒦+(r^2)))+a ((ℒ-a(1-(z^2))ℰ)/(𝒦-(a^2)(z^2))))

This separable equation can be solved analytically.

Implementation

The solutions have been implemented in the KerrGeodesics`ParallelTransport` subpackage.

<<"KerrGeodesics`ParallelTransport`"

M=1;
a0=0.9;

frame=KerrParallelTransportFrame[apex,Method->"Analytic"]

The KerrParallelTransportFrameFunction stores everything related to the parallel transported frame and the orbit. We extract the orbital trajectory using the same Key as KerrGeoOrbitFunction:

{t1,r1,θ1,ϕ1}=frame["Trajectory"];

Or the constants of motion :

frame["ConstantsOfMotion"]

The precession frequency (Subscript[Υ, ψ]) is calculated as part of the object:

frame["PrecessionFrequency"]

The frame itself is obtained by evaluating the KerrParallelTransportFrameFunction at an argument λ, which can be symbolic:

frame[λ]//MatrixForm

or numerical

frame[42.]//MatrixForm

Keys@frame[[-1]]

Nice plot

To generate a nice plot, I will use an example from a notebook created by Niels. This uses the GeneralRelativityTensors package to manipulate the generated frame.

Needs["GeneralRelativityTensors`"]

We will use the second leg of the tetrad for this example.

g=ToMetric["Kerr"]
e2=ToTensor["(e^((2)))",g,frame[λ][[2]],{-α}]

We will transform to Cartesian coordinates to make a nice plot

x=r Sin[θ] Cos[ϕ];
y=r Sin[θ] Sin[ϕ];
z=r Cos[θ];

transform={{1,0,0,0},{0,D[x,r],D[y,r],D[z,r]},{0,D[x,θ],D[y,θ],D[z,θ]},{0,D[x,ϕ],D[y,ϕ],D[z,ϕ]}};

e2Cartesian=TensorValues[e2[α]].transform;

norm3=Sqrt[e2Cartesian[[2;;4]].e2Cartesian[[2;;4]]];

Create the vectors from the spatial components of the tetrad leg

CreateVector[λ0_]:=Block[{a=a0,λ=λ0,r=r1[λ],θ=θ1[λ],ϕ=ϕ1[λ],vec},
vec=e2Cartesian;
Line[{{x-2vec[[2]]/norm3,y-2vec[[3]]/norm3,z-2vec[[4]]/norm3},{x,y,z}}//Re]
]

Do[line[λ0]=CreateVector[λ0];,{λ0,0,3.,0.1}]

Create a list of the position of the secondary

Table[pos[λ0]={x,y,z}/.{r->r1[λ0],θ->θ1[λ0],ϕ->ϕ1[λ0]}//Re,{λ0,0,3,0.1}];

Make a nice figure

Show[
ParametricPlot3D[{x,y,z}/.{r->r1[λ],θ->θ1[λ],ϕ->ϕ1[λ]}//Re,{λ,-3,3},ImageSize->700,Boxed->False,Axes->False,PlotStyle->Red,PlotRange->All,Lighting->"Neutral"],
Graphics3D[{Black,Specularity[.5],Sphere[{0,0,0},1+Sqrt[1-(a0^2)]]},Lighting->"Neutral"],
Graphics3D[Table[{Black,Specularity[.5],Sphere[pos[λ],0.2]},{λ,0,2,0.1}]],
Graphics3D[Table[{Blue,Thick,line[λ]},{λ,0.,2,0.1}]]
]

Exercises

Exercise 1: parameter evolution

Typically a code for producing Teukolsky fluxes will output (OverDot[ℰ]), (OverDot[ℒ]), and (OverDot[𝒬]) as a function of a, p, e, and x. To evolve an inspiral we need (OverDot[p]), (OverDot[e]), and (OverDot[x]).

Easy

For circular orbits in Schwarzschild, write a function that takes as an input a value for the radius r and energy flux (OverDot[ℰ]), and outputs r and (OverDot[r]).

Medium

For eccentric orbits in Schwarzschild, write a function that takes as an input a value for the semilatus rectum p and eccentricity e, together with the corresponding fluxes (OverDot[ℰ]) and (OverDot[ℒ]), and outputs (OverDot[p]) and (OverDot[e]).

Hard

For generic bound orbits in Kerr, write a function that takes as an input a, p, e, and x, together with the corresponding fluxes (OverDot[ℰ]), (OverDot[ℒ]), and (OverDot[𝒬]), and outputs (OverDot[p]), (OverDot[e]), and (OverDot[x]).

Exercise 2 : Orbital resonances

Resonant orbits are orbits for which two (or more) frequencies form a rational ratio, e.g. ((Subscript[Υ, r])/(Subscript[Υ, ϕ]))=(1/2). Besides being pretty, these orbits are of special interest.

Easy

In Schwarzschild, find a orbit with (Subscript[Υ, r]) : (Subscript[Υ, ϕ]) == 2 : 3. Plot this orbit.

Medium

In Kerr, find a orbit with  (Subscript[Υ, r]) : (Subscript[Υ, θ]) : (Subscript[Υ, ϕ]) ==1: 2 : 3. Plot this orbit.

Hard

The precession of a secondary spin introduces a fourth frequency. The following function allows easy access to this frequency from the parallelTransport subpackage.

In Kerr, find a orbit with  (Subscript[Υ, r]) : (Subscript[Υ, ψ]) : (Subscript[Υ, θ]) : (Subscript[Υ, ϕ]) ==2:3 : 4 : 6. Plot this orbit.

Mini - Project 1: edge - cases

Find edge cases (extreme spin, polar, marginally bound, etc.) where some of the KerrGeoXXX functions are broken. Implement exceptions for these cases.

Project 2: scattering orbits

Implement everything (or something) in the KerrGeodesics package for scattering orbits.

Project 3: plunging orbits

Implement everything (or something) in the KerrGeodesics package for plunging orbits.

Project 4: lightlike orbits

Implement everything (or something) in the KerrGeodesics package for lightlike orbits.