Comparison of E, L and Omega.ipynb

Kernel: SageMath (stable)

Comparisons of $E$, $L$ and $\Omega$ between Bardeen et al, Toshmatov et al and Lamy et al

This Jupyter/SageMath notebook is related to the article Lamy et al, arXiv:1802.01635.

version()

'SageMath version 8.1, Release Date: 2017-12-07'

%display latex

from sage.manifolds.utilities import simplify_chain_real as simpl

$E$, $L$ and $\Omega$ in Bardeen '72 (Kerr case)

var('r a b m0', domain='real')
var('th', latex_name=r'\theta', domain='real')
assume(r>0)
assume(m0>0)
assume(a>=0)

omega_pos_B=sqrt(m0)/(r^(3/2)+a*sqrt(m0)) # Omega computed in Bardeen et al '72 (Kerr) for co-rotating orbits
omega_neg_B=-sqrt(m0)/(r^(3/2)-a*sqrt(m0))# Omega computed in Bardeen et al '72 (Kerr) for counter-rotating orbits

E_pos_B=(r^(3/2)-2*m0*sqrt(r)+a*sqrt(m0))/(r^(3/4)*sqrt(r^(3/2)-3*m0*sqrt(r)+2*a*sqrt(m0))) # E computed in Bardeen et al '72 (Kerr) for co-rotating orbits
E_neg_B=(r^(3/2)-2*m0*sqrt(r)-a*sqrt(m0))/(r^(3/4)*sqrt(r^(3/2)-3*m0*sqrt(r)-2*a*sqrt(m0))) # E computed in Bardeen et al '72 (Kerr) for counter-rotating orbits

L_pos_B=(sqrt(m0)*(r^2-2*a*sqrt(m0)*sqrt(r)+a^2))/(r^(3/4)*sqrt(r^(3/2)-3*m0*sqrt(r)+2*a*sqrt(m0))) # L computed in Bardeen et al '72 (Kerr) for co-rotating orbits
L_neg_B=(-sqrt(m0)*(r^2+2*a*sqrt(m0)*sqrt(r)+a^2))/(r^(3/4)*sqrt(r^(3/2)-3*m0*sqrt(r)-2*a*sqrt(m0))) # L computed in Bardeen et al '72 (Kerr) for counter-rotating orbits

ratio_pos_B=L_pos_B/E_pos_B

ratio_pos_B.simplify_full()

$E$, $L$ and $\Omega$ in Toshmatov et al '17 (Kerr with $M(r)$)

Definitions

M=m0*r^3/(r^3+2*m0*b^2) #Let's assume the Hayward form for M(r), which boils down to Kerr for b=0
Sigma = r^2 + a^2*cos(th)^2
Delta = r^2 - 2*M*r + a^2
g00=-(1 - 2*r*M/Sigma)
g03=-2*a*r*sin(th)^2*M/Sigma
g33=(r^2 + a^2 + 2*r*(a*sin(th))^2*M/Sigma)*sin(th)^2
dg00=diff(g00,r)
dg03=diff(g03,r)
dg33=diff(g33,r)
Mdot=diff(M,r)

omega_pos_T=(sqrt(M/r-Mdot)/(a*sqrt(M/r-Mdot)+r)) # Omega computed in Toshmatov et al '17
omega_neg_T=sqrt(M/r-Mdot)/(a*sqrt(M/r-Mdot)-r)

E_pos_T=(r-2*M+a*sqrt(M/r-Mdot))/(r*(r-3*M+2*a*sqrt(M/r-Mdot)+r*Mdot)) # E computed in Toshmatov et al '17
E_neg_T=-(r-2*M+a*sqrt(M/r-Mdot))/(r*(r-3*M+2*a*sqrt(M/r-Mdot)+r*Mdot))

L_pos_T=(2*a*M-(r^2+a^2)*sqrt(M/r-Mdot))/(r*(r-3*M+2*a*sqrt(M/r-Mdot)+r*Mdot)) # L computed in Toshmatov et al '17
L_neg_T=-(2*a*M-(r^2+a^2)*sqrt(M/r-Mdot))/(r*(r-3*M+2*a*sqrt(M/r-Mdot)+r*Mdot))

Comparison with Bardeen '72

$\Omega$: agreement

omega_pos_B

simpl(omega_pos_T.subs({b:0, th:pi/2}))

bool(omega_pos_B == simpl(omega_pos_T.subs({b:0, th:pi/2})))

$E$: disagreement

simpl(E_pos_B)

simpl(E_pos_T.subs({b:0, th:pi/2}))

E_pos_B.subs({a:1},{m0:1},{r:3}).n()

E_pos_T.subs({a:1},{b:0},{m0:1},{r:3},{th:pi/2}).n() #b=0 boils down to Kerr 

E_neg_B.subs({a:1},{m0:1},{r:3}).n()

E_neg_T.subs({a:1},{b:0},{m0:1},{r:3},{th:pi/2}).n() #b=0 boils down to Kerr 

$L$: disagreement

simpl(L_pos_B)

simpl(L_pos_T.subs({b:0, th:pi/2}))

L_pos_B.subs({a:1},{m0:1},{r:3}).n()

L_pos_T.subs({a:1},{b:0},{m0:1},{r:3},{th:pi/2}).n() #b=0 boils down to Kerr 

$E$, $L$ and $\Omega$ in Lamy et al '18 (Kerr with $M(r)$)

Definitions

omega_pos=(-dg03+sqrt(dg03^2-dg00*dg33))/dg33 # Omega computed in Lamy et al '18
omega_neg=(-dg03-sqrt(dg03^2-dg00*dg33))/dg33

E_pos=-(g00+g03*omega_pos)/sqrt(-(g00+2*g03*omega_pos+g33*omega_pos^2)) #E computed in Lamy et al '18 
E_neg=-(g00+g03*omega_neg)/sqrt(-(g00+2*g03*omega_neg+g33*omega_neg^2))

L_pos=(g03+g33*omega_pos)/sqrt(-(g00+2*g03*omega_pos+g33*omega_pos^2)) # L computed in Lamy et al '18
L_neg=(g03+g33*omega_neg)/sqrt(-(g00+2*g03*omega_neg+g33*omega_neg^2))

Comparison with Bardeen '72

$\Omega$: agreement

simpl(omega_pos_B)

simpl(omega_pos.subs({b:0, th:pi/2}))

bool(simpl(omega_pos_B) == simpl(omega_pos.subs({b:0},{th:pi/2})))

$E$: agreement

E_pos_B1 = simpl(E_pos_B)
E_pos_B1

E_pos1 = simpl(E_pos.subs({b:0, th:pi/2})) #b=0 boils down to Kerr 
E_pos1

bool(E_pos1^2 == E_pos_B1^2)

$L$: agreement

L_pos_B1 = simpl(L_pos_B)
L_pos_B1

L_pos1 = simpl(L_pos.subs({b:0, th:pi/2})) #b=0 boils down to Kerr 
L_pos1

bool(L_pos1^2 == L_pos_B1^2)