Stable_circular_orbits.ipynb

Kernel: SageMath (stable)

Study of stable circular orbits in the equatorial plane

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

Spacetime

a0 = 0.7
b0 = 1

var('r a b')
var('th', latex_name=r'\theta')

m(b,r) = (1-unit_step(r))* r^3/(r^3 - 2*b^2) + \
         unit_step(r)*     r^3/(r^3 + 2*b^2)
m(b,r)

m(0,r)  # Kerr limit

plot(m(1, r), (r, -6, 6), axes_labels=[r"$r$", r"$m(r)$"], 
     title=r"$b=1$", title_pos=(0.9, 1.1))

Derivatives of the mass function

md(b, r) = diff(m(b, r), r)
md(b, r)

Because of the step function, the derivative contains Dirac delta functions. We get rid of them by means of the following Python function:

def remove_dirac(f):
    return f.subs({dirac_delta(r): 0}).simplify()

md(b, r) = remove_dirac(md(b, r))
md(b, r)

plot(md(1, r), (r, -6, 6), axes_labels=[r"$r$", r"$m'(r)$"],
     title=r"$b=1$", title_pos=(0.9, 1))

mdd(b, r) = remove_dirac(diff(md(b, r), r))
mdd(b, r)

plot(mdd(1, r), (r, -6, 6), axes_labels=[r"$r$", r"$m''\,(r)$"],
     title=r"$b=1$", title_pos=(0.9, 1))

Horizons

D=r^2+a^2-2*m(b,r)*r
plot(D.subs({a:a0},{b:b0}),(r,-5,5))

#find_root(D.subs({a:a0},{b:b0}),0,1)

Circular orbits in the equatorial plane (Condition1>0 & Condition2>0)

$E$, $L$ and $\Omega$ of a test particle in the equatorial plane

#t_dot(a,b,E,L,r,th) = ((r^2+a^2)*(E*(r^2+a^2)-a*L)/D - a*(a*E*sin(th)^2-L))/Sigma
#phi_dot(a,b,E,L,r,th) = (a*(E*(r^2+a^2)-a*L)/D - (a*E-L/sin(th)^2))/Sigma
#mu0(a,b,E,L,r,th) = -(g00*t_dot^2 + 2*g03*t_dot*phi_dot + g33*phi_dot^2) #normalization (m0=0 for a photon)

var('m0','b')

M = function('M')(r)
#M=r^3/(r^3+2*b^2)
Sigma = r^2+a^2*cos(th)^2
Delta = r^2+a^2-2*M*r
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

diff(r^3/(r^3+2*b^2),r)

var('e l q');

dg00=diff(g00,r)
dg03=diff(g03,r)
dg33=diff(g33,r)
Mdot=diff(M,r)
omega_pos=(-dg03+sqrt(dg03^2-dg00*dg33))/dg33 
omega_neg=(-dg03-sqrt(dg03^2-dg00*dg33))/dg33

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

E_pos=-(g00+g03*omega_pos)/sqrt(-(g00+2*g03*omega_pos+g33*omega_pos^2)) 
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_neg=(g03+g33*omega_neg)/sqrt(-(g00+2*g03*omega_neg+g33*omega_neg^2))

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

L_pos.subs({th:pi/2,b:0}).simplify_full()

L_neg.subs({th:pi/2,b:0}).simplify_full()

Conditions 1 & 2

reset('M')

M=m(b,r)
Md = md(b ,r)
Mdd = mdd(b, r)

Condition1 = M/r - Md
Condition1

Condition1.subs({b:0})  # Kerr limit

Condition2_pos=-(a^2*r^2*Md^2 - r^4 + 3*a^2*M^2 + 3*(a^2*r + r^3)*M - (3*a^2*r^2 + r^4 + 4*a^2*r*M)*Md - 2*((a^3 + 3*a*r^2)*M - (a^3*r + a*r^3)*Md*sqrt(-(Md - M/r))))

Condition2_neg=-(a^2*r^2*Md^2 - r^4 + 3*a^2*M^2 + 3*(a^2*r + r^3)*M - (3*a^2*r^2 + r^4 + 4*a^2*r*M)*Md + 2*((a^3 + 3*a*r^2)*M - (a^3*r + a*r^3)*Md*sqrt(-(Md - M/r))))

Condition2_pos.subs({b:0})  # Kerr limit

Condition2_pos.subs({a:0},{b:0})  # Schwarzschild limit

p1 = plot(20, -1.587,0, fill=True, fillalpha=0.3)
p2 = plot(20, 1.587, 5, fill=True, fillalpha=0.3)
p3 = plot(20, 1.31, 1.34, fill=True, fillcolor='black', fillalpha=1)
p4 = plot(20, 0.89, 0.92, fill=True, fillcolor='black', fillalpha=1)

p = plot((Condition1.subs({b:b0}),
          Condition2_pos.subs({a:a0},{b:b0}),
          Condition2_neg.subs({a:a0},{b:b0})), (r,-3,5), ymax=12, ymin=-1.5,
         legend_label=['$A(r)$','$B_+(r)$','$B_-(r)$'], color=('blue','red','green'), 
         axes_labels=[r'$r/m$',''])
graph = p+p1+p2#+p3+p4
show(graph)

verbose 0 (3749: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 76 points.
verbose 0 (3749: plot.py, generate_plot_points) Last error message: ''
verbose 0 (3749: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 75 points.
verbose 0 (3749: plot.py, generate_plot_points) Last error message: ''

graph.save('Circular_orbits_a=07_b=1.pdf')

#find_root(Condition2_pos.subs({a:a0},{b:b0}),-1,-0.3)

Condition2_pos.subs({a:a0},{b:b0},{r:1.587})

Stable circular orbits in the equatorial plane (Condition1>0, Condition2>0 & Condition 3<0), ISCO

Defining $R$

var('E L Q m0')

reset('M')
M=m(b,r)
#Md = md(b ,r)
#Mdd = mdd(b, r)
Sigma = r^2+a^2*cos(th)^2
Delta = r^2+a^2-2*M*r
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)
omega_pos=(-dg03+sqrt(dg03^2-dg00*dg33))/dg33 
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_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_neg=(g03+g33*omega_neg)/sqrt(-(g00+2*g03*omega_neg+g33*omega_neg^2))

R = ((r^2+a^2)*E - a*L)^2 - Delta*((a*E-L)^2 + m0^2*r^2 + Q)
R

R.subs({b:0}).simplify_full()   # Kerr limit

R.subs({a:0},{b:0}).simplify_full()  # Schwarzschild limit

Rd = remove_dirac(diff(R, r))
Rd

Rd.subs({b:0}).simplify_full() # Kerr limit

Rd.subs({a:0},{b:0}).simplify_full()  # Schwarzschild limit

Rdd = remove_dirac(diff(Rd, r))
Rdd

Rdd.subs({b:0}).simplify_full() # Kerr limit

Rdd.subs({a:0},{b:0}).simplify_full()  # Schwarzschild limit

Rd2=Rdd.subs({E:e},{L:l},{Q:q})
Rd2_normalized=Rd2.subs({m0:1}).simplify() #normalization by m0

Rd2_normalized.subs({b:0}).simplify_full()  # Kerr limit

Rd2_normalized.subs({a:0},{b:0}).simplify_full()  # Schwarzschild limit

Conditions 1, 2 & 3

q1=0 #value of q=Q/m0^2 corresponding to a motion in the equatorial plane 
e1_pos=E_pos.subs({a:a0},{b:b0}) #value of e=E/m0 corresponding to a circular orbit
e1_neg=E_neg.subs({a:a0},{b:b0})
l1_pos=L_pos.subs({a:a0},{b:b0}) #value of l=L/m0 corresponding to a circular orbit
l1_neg=L_neg.subs({a:a0},{b:b0})

Condition3_pos = Rd2_normalized.subs({e:e1_pos},{l:l1_pos},{q:q1}).simplify()
Condition3_neg = Rd2_normalized.subs({e:e1_neg},{l:l1_neg},{q:q1}).simplify() 

Condition3_neg.subs({a:a0},{b:b0},{th:pi/2},{r:1})

plot(Condition3_neg.subs({a:a0},{b:b0},{th:pi/2}),(r,0,10))
#pc3.save('ISCO_a0_b0.pdf')

verbose 0 (3749: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 67 points.
verbose 0 (3749: plot.py, generate_plot_points) Last error message: 'unable to convert -488.7212108566773 + 1.1209327348532101e-13*I to float; use abs() or real_part() as desired'

plot((Condition1.subs({b:b0}), 
      Condition2_neg.subs({a:a0},{b:b0}),
      Condition3_neg.subs({a:a0},{b:b0},{th:pi/2})), (r,-3,10), ymax=10, ymin=-10,
     legend_label=['A(r)','B(r)','C(r)'], color=('blue','red','green'), 
     title='Counter-rotating', axes_labels=[r'$r$',''])

verbose 0 (3749: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 47 points.
verbose 0 (3749: plot.py, generate_plot_points) Last error message: ''
verbose 0 (3749: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 72 points.
verbose 0 (3749: plot.py, generate_plot_points) Last error message: ''

p_pos = plot((Condition1.subs({b:b0}), 
              Condition2_pos.subs({a:a0},{b:b0}),
              Condition3_pos.subs({a:a0},{b:b0},{th:pi/2})), (r,-3,10), ymax=10, ymin=-10,
             legend_label=['A(r)','B(r)','C(r)'], color=('blue','red','green'), 
             title='Co-rotating',axes_labels=[r'$r$',''])
p_neg = plot((Condition1.subs({b:b0}), 
              Condition2_neg.subs({a:a0},{b:b0}),
              Condition3_neg.subs({a:a0},{b:b0},{th:pi/2})), (r,-3,10), ymax=10, ymin=-10,
             legend_label=['A(r)','B(r)','C(r)'], color=('blue','red','green'), 
             title='Counter-rotating', axes_labels=[r'$r$',''])
show(p_pos)
show(p_neg)

verbose 0 (3749: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 46 points.
verbose 0 (3749: plot.py, generate_plot_points) Last error message: ''
verbose 0 (3749: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 57 points.
verbose 0 (3749: plot.py, generate_plot_points) Last error message: ''
verbose 0 (3749: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 46 points.
verbose 0 (3749: plot.py, generate_plot_points) Last error message: ''
verbose 0 (3749: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 73 points.
verbose 0 (3749: plot.py, generate_plot_points) Last error message: ''

find_root(Condition2_pos.subs({a:a0},{b:b0},{th:pi/2}),1.001,2)

plot(Condition3_pos.subs({a:a0},{b:b0},{th:pi/2}),(r,0,2))

verbose 0 (3749: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 157 points.
verbose 0 (3749: plot.py, generate_plot_points) Last error message: 'unable to convert -8.32736020802459 + 2.683660004211231*I to float; use abs() or real_part() as desired'

Condition3_pos.subs({a:a0},{b:b0},{th:pi/2},{r:1.04})

Frequency of the orbits at the ISCO

%display latex

var('r a b M0 m')
var('th', latex_name=r'\theta')

%display latex
var('r a b M0 m')
var('th', latex_name=r'\theta')

m0=4.31*10^6*2*10^30
M(b,r) = m0*r^3/(r^3 + 2*m0*b^2)
#M=function('M',r)
Sigma = r^2 + a^2*cos(th)^2
Delta = r^2 - 2*M*r + a^2
g00(a,b,r) = -(1 - 2*r*M/Sigma)
g03(a,b,r) = -2*a*r*sin(th)^2*M/Sigma
g33(a,b,r) = (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)
omega_pos=(-dg03+sqrt(dg03^2-dg00*dg33))/dg33 
omega_neg=(-dg03-sqrt(dg03^2-dg00*dg33))/dg33

a1=m0*0.9
b1=m0*1
r1_pos=m0*1.59
r1_neg=m0*8.43

om_pos=omega_pos.subs({a:a1},{b:b1},{r:r1_pos},{th:pi/2})
(3*10^8)^3/(6.67*10^(-11))*om_pos

om_neg=omega_neg.subs({a:a1},{b:b1},{r:r1_neg},{th:pi/2})
(3*10^8)^3/(6.67*10^(-11))*om_neg

General case: only $R(r,E,L,Q) \geq 0$ allowed

reset('m0')

var('m0')

R.subs({b:0},{m0:1}).simplify_full() # Kerr limit

R.subs({a:0},{b:0},{m0:0}).simplify_full()  # Schwarzschild limit

a2 = 0.9
b2 = 1
e2 = 1
l2 = 2
q2 = -1
plot0=plot(R.subs({m0:0},{a:a2},{b:b2},{E:e2},{L:l2},{Q:q2}), (r,-2,2), axes_labels = [r'$r/m$', r'$R(r,E_1,L_1,Q_1)$'])
show(plot0)

R.subs({m0:0},{a:a2},{b:b2},{E:e2},{L:l2},{Q:q2})

#find_root(R_null(a2,b2,e2,l2,q2,r),-2,2)

R.subs({m0:0},{a:a2},{b:b2},{E:e2},{L:l2},{Q:q2},{r:1})

#p1 = plot(11, -0.665666,0.713731,ymax=8, fill=True, fillalpha=0.3)
p2 = plot(25, -0.054222, 2, ymax=22, ymin=-10, fill=True, fillalpha=0.3)
p3 = plot(11,-2, -1.446841, fill=True, fillalpha=0.3)
p4 = plot(25, 1.425, 1.44, fill=True, fillcolor='black', fillalpha=1)
p5 = plot(25, 0.545, 0.56, fill=True, fillcolor='black', fillalpha=1)
plot2=plot0+p2+p4+p5
show(plot2)

#plot2.save('R_photons_dark_zone_b=0.pdf')