SharedStable_circular_orbits.ipynbOpen in CoCalc

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')
\newcommand{\Bold}[1]{\mathbf{#1}}{\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)
\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{r^{3} {\left(\mathrm{u}\left(r\right) - 1\right)}}{r^{3} - 2 \, b^{2}} + \frac{r^{3} \mathrm{u}\left(r\right)}{r^{3} + 2 \, b^{2}}
m(0,r)  # Kerr limit
\newcommand{\Bold}[1]{\mathbf{#1}}1
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)
\newcommand{\Bold}[1]{\mathbf{#1}}\frac{3 \, r^{5} {\left(\mathrm{u}\left(r\right) - 1\right)}}{{\left(r^{3} - 2 \, b^{2}\right)}^{2}} - \frac{3 \, r^{5} \mathrm{u}\left(r\right)}{{\left(r^{3} + 2 \, b^{2}\right)}^{2}} + \frac{r^{3} \delta\left(r\right)}{r^{3} + 2 \, b^{2}} - \frac{r^{3} \delta\left(r\right)}{r^{3} - 2 \, b^{2}} - \frac{3 \, r^{2} {\left(\mathrm{u}\left(r\right) - 1\right)}}{r^{3} - 2 \, b^{2}} + \frac{3 \, r^{2} \mathrm{u}\left(r\right)}{r^{3} + 2 \, b^{2}}

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)
\newcommand{\Bold}[1]{\mathbf{#1}}\frac{3 \, r^{5} {\left(\mathrm{u}\left(r\right) - 1\right)}}{{\left(r^{3} - 2 \, b^{2}\right)}^{2}} - \frac{3 \, r^{5} \mathrm{u}\left(r\right)}{{\left(r^{3} + 2 \, b^{2}\right)}^{2}} - \frac{3 \, r^{2} {\left(\mathrm{u}\left(r\right) - 1\right)}}{r^{3} - 2 \, b^{2}} + \frac{3 \, r^{2} \mathrm{u}\left(r\right)}{r^{3} + 2 \, b^{2}}
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)
\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{18 \, r^{7} {\left(\mathrm{u}\left(r\right) - 1\right)}}{{\left(r^{3} - 2 \, b^{2}\right)}^{3}} + \frac{18 \, r^{7} \mathrm{u}\left(r\right)}{{\left(r^{3} + 2 \, b^{2}\right)}^{3}} + \frac{24 \, r^{4} {\left(\mathrm{u}\left(r\right) - 1\right)}}{{\left(r^{3} - 2 \, b^{2}\right)}^{2}} - \frac{24 \, r^{4} \mathrm{u}\left(r\right)}{{\left(r^{3} + 2 \, b^{2}\right)}^{2}} - \frac{6 \, r {\left(\mathrm{u}\left(r\right) - 1\right)}}{r^{3} - 2 \, b^{2}} + \frac{6 \, r \mathrm{u}\left(r\right)}{r^{3} + 2 \, b^{2}}
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')
\newcommand{\Bold}[1]{\mathbf{#1}}\left(m_{0}, b\right)
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)
\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{3 \, r^{5}}{{\left(r^{3} + 2 \, b^{2}\right)}^{2}} + \frac{3 \, r^{2}}{r^{3} + 2 \, b^{2}}
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}))
\newcommand{\Bold}[1]{\mathbf{#1}}\frac{a r \frac{\partial}{\partial r}M\left(r\right) + \sqrt{-r \frac{\partial}{\partial r}M\left(r\right) + M\left(r\right)} r^{\frac{3}{2}} - a M\left(r\right)}{a^{2} r \frac{\partial}{\partial r}M\left(r\right) + r^{3} - a^{2} M\left(r\right)}
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()
\newcommand{\Bold}[1]{\mathbf{#1}}\frac{a^{2} r \frac{\partial}{\partial r}M\left(r\right) + r^{3} + 2 \, \sqrt{-r \frac{\partial}{\partial r}M\left(r\right) + M\left(r\right)} a \sqrt{r} M\left(r\right) - {\left(a^{2} + 2 \, r^{2}\right)} M\left(r\right)}{\sqrt{-a^{2} r^{3} \frac{\partial}{\partial r}M\left(r\right)^{2} + r^{5} - 3 \, a^{2} r M\left(r\right)^{2} + 2 \, {\left({\left(a^{3} + 3 \, a r^{2}\right)} M\left(r\right) - {\left(a^{3} r + a r^{3}\right)} \frac{\partial}{\partial r}M\left(r\right)\right)} \sqrt{-r \frac{\partial}{\partial r}M\left(r\right) + M\left(r\right)} \sqrt{r} - 3 \, {\left(a^{2} r^{2} + r^{4}\right)} M\left(r\right) + {\left(3 \, a^{2} r^{3} + r^{5} + 4 \, a^{2} r^{2} M\left(r\right)\right)} \frac{\partial}{\partial r}M\left(r\right)} \sqrt{r}}
L_pos.subs({th:pi/2,b:0}).simplify_full()
\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{{\left(a^{3} + 3 \, a r^{2}\right)} M\left(r\right) - {\left(a^{3} r + a r^{3}\right)} \frac{\partial}{\partial r}M\left(r\right) - {\left(a^{2} r^{2} + r^{4} + 2 \, a^{2} r M\left(r\right)\right)} \sqrt{-\frac{r \frac{\partial}{\partial r}M\left(r\right) - M\left(r\right)}{r}}}{{\left(a^{2} r \frac{\partial}{\partial r}M\left(r\right) + r^{3} - a^{2} M\left(r\right)\right)} \sqrt{-\frac{a^{2} r^{4} \frac{\partial}{\partial r}M\left(r\right)^{2} - r^{6} + 3 \, a^{2} r^{2} M\left(r\right)^{2} + 3 \, {\left(a^{2} r^{3} + r^{5}\right)} M\left(r\right) - {\left(3 \, a^{2} r^{4} + r^{6} + 4 \, a^{2} r^{3} M\left(r\right)\right)} \frac{\partial}{\partial r}M\left(r\right) - 2 \, {\left({\left(a^{3} r^{2} + 3 \, a r^{4}\right)} M\left(r\right) - {\left(a^{3} r^{3} + a r^{5}\right)} \frac{\partial}{\partial r}M\left(r\right)\right)} \sqrt{-\frac{r \frac{\partial}{\partial r}M\left(r\right) - M\left(r\right)}{r}}}{a^{4} r^{2} \frac{\partial}{\partial r}M\left(r\right)^{2} + r^{6} - 2 \, a^{2} r^{3} M\left(r\right) + a^{4} M\left(r\right)^{2} + 2 \, {\left(a^{2} r^{4} - a^{4} r M\left(r\right)\right)} \frac{\partial}{\partial r}M\left(r\right)}}}
L_neg.subs({th:pi/2,b:0}).simplify_full()
\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{{\left(a^{3} + 3 \, a r^{2}\right)} M\left(r\right) - {\left(a^{3} r + a r^{3}\right)} \frac{\partial}{\partial r}M\left(r\right) + {\left(a^{2} r^{2} + r^{4} + 2 \, a^{2} r M\left(r\right)\right)} \sqrt{-\frac{r \frac{\partial}{\partial r}M\left(r\right) - M\left(r\right)}{r}}}{{\left(a^{2} r \frac{\partial}{\partial r}M\left(r\right) + r^{3} - a^{2} M\left(r\right)\right)} \sqrt{-\frac{a^{2} r^{4} \frac{\partial}{\partial r}M\left(r\right)^{2} - r^{6} + 3 \, a^{2} r^{2} M\left(r\right)^{2} + 3 \, {\left(a^{2} r^{3} + r^{5}\right)} M\left(r\right) - {\left(3 \, a^{2} r^{4} + r^{6} + 4 \, a^{2} r^{3} M\left(r\right)\right)} \frac{\partial}{\partial r}M\left(r\right) + 2 \, {\left({\left(a^{3} r^{2} + 3 \, a r^{4}\right)} M\left(r\right) - {\left(a^{3} r^{3} + a r^{5}\right)} \frac{\partial}{\partial r}M\left(r\right)\right)} \sqrt{-\frac{r \frac{\partial}{\partial r}M\left(r\right) - M\left(r\right)}{r}}}{a^{4} r^{2} \frac{\partial}{\partial r}M\left(r\right)^{2} + r^{6} - 2 \, a^{2} r^{3} M\left(r\right) + a^{4} M\left(r\right)^{2} + 2 \, {\left(a^{2} r^{4} - a^{4} r M\left(r\right)\right)} \frac{\partial}{\partial r}M\left(r\right)}}}

Conditions 1 & 2

reset('M')
M=m(b,r)
Md = md(b ,r)
Mdd = mdd(b, r)
Condition1 = M/r - Md
Condition1
\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{3 \, r^{5} {\left(\mathrm{u}\left(r\right) - 1\right)}}{{\left(r^{3} - 2 \, b^{2}\right)}^{2}} + \frac{3 \, r^{5} \mathrm{u}\left(r\right)}{{\left(r^{3} + 2 \, b^{2}\right)}^{2}} + \frac{3 \, r^{2} {\left(\mathrm{u}\left(r\right) - 1\right)}}{r^{3} - 2 \, b^{2}} - \frac{3 \, r^{2} \mathrm{u}\left(r\right)}{r^{3} + 2 \, b^{2}} - \frac{\frac{r^{3} {\left(\mathrm{u}\left(r\right) - 1\right)}}{r^{3} - 2 \, b^{2}} - \frac{r^{3} \mathrm{u}\left(r\right)}{r^{3} + 2 \, b^{2}}}{r}
Condition1.subs({b:0})  # Kerr limit
\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{r}
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
\newcommand{\Bold}[1]{\mathbf{#1}}r^{4} + 2 \, a^{3} - 3 \, a^{2} r + 6 \, a r^{2} - 3 \, r^{3} - 3 \, a^{2}
Condition2_pos.subs({a:0},{b:0})  # Schwarzschild limit
\newcommand{\Bold}[1]{\mathbf{#1}}r^{4} - 3 \, r^{3}
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})
\newcommand{\Bold}[1]{\mathbf{#1}}8.52502568341426 - 0.0409208183571186i

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

Defining R

var('E L Q m0')
\newcommand{\Bold}[1]{\mathbf{#1}}\left(E, L, Q, m_{0}\right)
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
\newcommand{\Bold}[1]{\mathbf{#1}}{\left({\left(a^{2} + r^{2}\right)} E - L a\right)}^{2} - {\left(m_{0}^{2} r^{2} + {\left(E a - L\right)}^{2} + Q\right)} {\left(a^{2} + 2 \, {\left(\frac{r^{3} {\left(\mathrm{u}\left(r\right) - 1\right)}}{r^{3} - 2 \, b^{2}} - \frac{r^{3} \mathrm{u}\left(r\right)}{r^{3} + 2 \, b^{2}}\right)} r + r^{2}\right)}
R.subs({b:0}).simplify_full()   # Kerr limit
\newcommand{\Bold}[1]{\mathbf{#1}}2 \, m_{0}^{2} r^{3} + {\left(E^{2} - m_{0}^{2}\right)} r^{4} - Q a^{2} + {\left(E^{2} a^{2} - a^{2} m_{0}^{2} - L^{2} - Q\right)} r^{2} + 2 \, {\left(E^{2} a^{2} - 2 \, E L a + L^{2} + Q\right)} r
R.subs({a:0},{b:0}).simplify_full()  # Schwarzschild limit
\newcommand{\Bold}[1]{\mathbf{#1}}2 \, m_{0}^{2} r^{3} + {\left(E^{2} - m_{0}^{2}\right)} r^{4} - {\left(L^{2} + Q\right)} r^{2} + 2 \, {\left(L^{2} + Q\right)} r
Rd = remove_dirac(diff(R, r))
Rd
\newcommand{\Bold}[1]{\mathbf{#1}}-2 \, {\left(a^{2} + 2 \, {\left(\frac{r^{3} {\left(\mathrm{u}\left(r\right) - 1\right)}}{r^{3} - 2 \, b^{2}} - \frac{r^{3} \mathrm{u}\left(r\right)}{r^{3} + 2 \, b^{2}}\right)} r + r^{2}\right)} m_{0}^{2} r + 4 \, {\left({\left(a^{2} + r^{2}\right)} E - L a\right)} E r - 2 \, {\left(m_{0}^{2} r^{2} + {\left(E a - L\right)}^{2} + Q\right)} {\left(\frac{r^{3} {\left(\mathrm{u}\left(r\right) - 1\right)}}{r^{3} - 2 \, b^{2}} - \frac{r^{3} \mathrm{u}\left(r\right)}{r^{3} + 2 \, b^{2}} - 3 \, {\left(\frac{r^{5} {\left(\mathrm{u}\left(r\right) - 1\right)}}{{\left(r^{3} - 2 \, b^{2}\right)}^{2}} - \frac{r^{5} \mathrm{u}\left(r\right)}{{\left(r^{3} + 2 \, b^{2}\right)}^{2}} - \frac{r^{2} {\left(\mathrm{u}\left(r\right) - 1\right)}}{r^{3} - 2 \, b^{2}} + \frac{r^{2} \mathrm{u}\left(r\right)}{r^{3} + 2 \, b^{2}}\right)} r + r\right)}
Rd.subs({b:0}).simplify_full() # Kerr limit
\newcommand{\Bold}[1]{\mathbf{#1}}2 \, E^{2} a^{2} + 6 \, m_{0}^{2} r^{2} + 4 \, {\left(E^{2} - m_{0}^{2}\right)} r^{3} - 4 \, E L a + 2 \, L^{2} + 2 \, {\left(E^{2} a^{2} - a^{2} m_{0}^{2} - L^{2} - Q\right)} r + 2 \, Q
Rd.subs({a:0},{b:0}).simplify_full()  # Schwarzschild limit
\newcommand{\Bold}[1]{\mathbf{#1}}6 \, m_{0}^{2} r^{2} + 4 \, {\left(E^{2} - m_{0}^{2}\right)} r^{3} + 2 \, L^{2} - 2 \, {\left(L^{2} + Q\right)} r + 2 \, Q
Rdd = remove_dirac(diff(Rd, r))
Rdd
\newcommand{\Bold}[1]{\mathbf{#1}}-8 \, {\left(\frac{r^{3} {\left(\mathrm{u}\left(r\right) - 1\right)}}{r^{3} - 2 \, b^{2}} - \frac{r^{3} \mathrm{u}\left(r\right)}{r^{3} + 2 \, b^{2}} - 3 \, {\left(\frac{r^{5} {\left(\mathrm{u}\left(r\right) - 1\right)}}{{\left(r^{3} - 2 \, b^{2}\right)}^{2}} - \frac{r^{5} \mathrm{u}\left(r\right)}{{\left(r^{3} + 2 \, b^{2}\right)}^{2}} - \frac{r^{2} {\left(\mathrm{u}\left(r\right) - 1\right)}}{r^{3} - 2 \, b^{2}} + \frac{r^{2} \mathrm{u}\left(r\right)}{r^{3} + 2 \, b^{2}}\right)} r + r\right)} m_{0}^{2} r + 8 \, E^{2} r^{2} - 2 \, {\left(a^{2} + 2 \, {\left(\frac{r^{3} {\left(\mathrm{u}\left(r\right) - 1\right)}}{r^{3} - 2 \, b^{2}} - \frac{r^{3} \mathrm{u}\left(r\right)}{r^{3} + 2 \, b^{2}}\right)} r + r^{2}\right)} m_{0}^{2} + 2 \, {\left(m_{0}^{2} r^{2} + {\left(E a - L\right)}^{2} + Q\right)} {\left(\frac{6 \, r^{5} {\left(\mathrm{u}\left(r\right) - 1\right)}}{{\left(r^{3} - 2 \, b^{2}\right)}^{2}} - \frac{6 \, r^{5} \mathrm{u}\left(r\right)}{{\left(r^{3} + 2 \, b^{2}\right)}^{2}} - 6 \, {\left(\frac{3 \, r^{7} {\left(\mathrm{u}\left(r\right) - 1\right)}}{{\left(r^{3} - 2 \, b^{2}\right)}^{3}} - \frac{3 \, r^{7} \mathrm{u}\left(r\right)}{{\left(r^{3} + 2 \, b^{2}\right)}^{3}} - \frac{4 \, r^{4} {\left(\mathrm{u}\left(r\right) - 1\right)}}{{\left(r^{3} - 2 \, b^{2}\right)}^{2}} + \frac{4 \, r^{4} \mathrm{u}\left(r\right)}{{\left(r^{3} + 2 \, b^{2}\right)}^{2}} + \frac{r {\left(\mathrm{u}\left(r\right) - 1\right)}}{r^{3} - 2 \, b^{2}} - \frac{r \mathrm{u}\left(r\right)}{r^{3} + 2 \, b^{2}}\right)} r - \frac{6 \, r^{2} {\left(\mathrm{u}\left(r\right) - 1\right)}}{r^{3} - 2 \, b^{2}} + \frac{6 \, r^{2} \mathrm{u}\left(r\right)}{r^{3} + 2 \, b^{2}} - 1\right)} + 4 \, {\left({\left(a^{2} + r^{2}\right)} E - L a\right)} E
Rdd.subs({b:0}).simplify_full() # Kerr limit
\newcommand{\Bold}[1]{\mathbf{#1}}2 \, E^{2} a^{2} - 2 \, a^{2} m_{0}^{2} + 12 \, m_{0}^{2} r + 12 \, {\left(E^{2} - m_{0}^{2}\right)} r^{2} - 2 \, L^{2} - 2 \, Q
Rdd.subs({a:0},{b:0}).simplify_full()  # Schwarzschild limit
\newcommand{\Bold}[1]{\mathbf{#1}}12 \, m_{0}^{2} r + 12 \, {\left(E^{2} - m_{0}^{2}\right)} r^{2} - 2 \, L^{2} - 2 \, Q
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
\newcommand{\Bold}[1]{\mathbf{#1}}2 \, a^{2} e^{2} + 12 \, {\left(e^{2} - 1\right)} r^{2} - 2 \, a^{2} - 2 \, l^{2} - 2 \, q + 12 \, r
Rd2_normalized.subs({a:0},{b:0}).simplify_full()  # Schwarzschild limit
\newcommand{\Bold}[1]{\mathbf{#1}}12 \, {\left(e^{2} - 1\right)} r^{2} - 2 \, l^{2} - 2 \, q + 12 \, r

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})
\newcommand{\Bold}[1]{\mathbf{#1}}NaN
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)
\newcommand{\Bold}[1]{\mathbf{#1}}0.0
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})
\newcommand{\Bold}[1]{\mathbf{#1}}1.62683109949525 + 1.36936138033619i

Frequency of the orbits at the ISCO

%display latex
var('r a b M0 m')
var('th', latex_name=r'\theta')
\newcommand{\Bold}[1]{\mathbf{#1}}{\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
\newcommand{\Bold}[1]{\mathbf{#1}}\left( a, b, r \right) \ {\mapsto} \ 0.00107187558536188
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
\newcommand{\Bold}[1]{\mathbf{#1}}\left( a, b, r \right) \ {\mapsto} \ -0.00197809424213708

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

reset('m0')
var('m0')
\newcommand{\Bold}[1]{\mathbf{#1}}m_{0}
R.subs({b:0},{m0:1}).simplify_full() # Kerr limit
\newcommand{\Bold}[1]{\mathbf{#1}}{\left(E^{2} - 1\right)} r^{4} - Q a^{2} + {\left({\left(E^{2} - 1\right)} a^{2} - L^{2} - Q\right)} r^{2} + 2 \, r^{3} + 2 \, {\left(E^{2} a^{2} - 2 \, E L a + L^{2} + Q\right)} r
R.subs({a:0},{b:0},{m0:0}).simplify_full()  # Schwarzschild limit
\newcommand{\Bold}[1]{\mathbf{#1}}E^{2} r^{4} - {\left(L^{2} + Q\right)} r^{2} + 2 \, {\left(L^{2} + Q\right)} r
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})
\newcommand{\Bold}[1]{\mathbf{#1}}{\left(r^{2} - 0.990000000000000\right)}^{2} - 0.420000000000000 \, {\left(\frac{r^{3} {\left(\mathrm{u}\left(r\right) - 1\right)}}{r^{3} - 2} - \frac{r^{3} \mathrm{u}\left(r\right)}{r^{3} + 2}\right)} r - 0.210000000000000 \, r^{2} - 0.170100000000000
#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})
\newcommand{\Bold}[1]{\mathbf{#1}}-0.240000000000000
#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')