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')
θ{\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)
r3(u(r)1)r32b2+r3u(r)r3+2b2-\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
11
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)
3r5(u(r)1)(r32b2)23r5u(r)(r3+2b2)2+r3δ(r)r3+2b2r3δ(r)r32b23r2(u(r)1)r32b2+3r2u(r)r3+2b2\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)
3r5(u(r)1)(r32b2)23r5u(r)(r3+2b2)23r2(u(r)1)r32b2+3r2u(r)r3+2b2\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)
18r7(u(r)1)(r32b2)3+18r7u(r)(r3+2b2)3+24r4(u(r)1)(r32b2)224r4u(r)(r3+2b2)26r(u(r)1)r32b2+6ru(r)r3+2b2-\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)

EE, LL 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')
(m0,b)\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)
3r5(r3+2b2)2+3r2r3+2b2-\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}))
arrM(r)+rrM(r)+M(r)r32aM(r)a2rrM(r)+r3a2M(r)\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()
a2rrM(r)+r3+2rrM(r)+M(r)arM(r)(a2+2r2)M(r)a2r3rM(r)2+r53a2rM(r)2+2((a3+3ar2)M(r)(a3r+ar3)rM(r))rrM(r)+M(r)r3(a2r2+r4)M(r)+(3a2r3+r5+4a2r2M(r))rM(r)r\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()
(a3+3ar2)M(r)(a3r+ar3)rM(r)(a2r2+r4+2a2rM(r))rrM(r)M(r)r(a2rrM(r)+r3a2M(r))a2r4rM(r)2r6+3a2r2M(r)2+3(a2r3+r5)M(r)(3a2r4+r6+4a2r3M(r))rM(r)2((a3r2+3ar4)M(r)(a3r3+ar5)rM(r))rrM(r)M(r)ra4r2rM(r)2+r62a2r3M(r)+a4M(r)2+2(a2r4a4rM(r))rM(r)-\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()
(a3+3ar2)M(r)(a3r+ar3)rM(r)+(a2r2+r4+2a2rM(r))rrM(r)M(r)r(a2rrM(r)+r3a2M(r))a2r4rM(r)2r6+3a2r2M(r)2+3(a2r3+r5)M(r)(3a2r4+r6+4a2r3M(r))rM(r)+2((a3r2+3ar4)M(r)(a3r3+ar5)rM(r))rrM(r)M(r)ra4r2rM(r)2+r62a2r3M(r)+a4M(r)2+2(a2r4a4rM(r))rM(r)-\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
3r5(u(r)1)(r32b2)2+3r5u(r)(r3+2b2)2+3r2(u(r)1)r32b23r2u(r)r3+2b2r3(u(r)1)r32b2r3u(r)r3+2b2r-\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
1r\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
r4+2a33a2r+6ar23r33a2r^{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
r43r3r^{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})
8.525025683414260.0409208183571186i8.52502568341426 - 0.0409208183571186i

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

Defining RR

var('E L Q m0')
(E,L,Q,m0)\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
((a2+r2)ELa)2(m02r2+(EaL)2+Q)(a2+2(r3(u(r)1)r32b2r3u(r)r3+2b2)r+r2){\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
2m02r3+(E2m02)r4Qa2+(E2a2a2m02L2Q)r2+2(E2a22ELa+L2+Q)r2 \, 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
2m02r3+(E2m02)r4(L2+Q)r2+2(L2+Q)r2 \, 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
2(a2+2(r3(u(r)1)r32b2r3u(r)r3+2b2)r+r2)m02r+4((a2+r2)ELa)Er2(m02r2+(EaL)2+Q)(r3(u(r)1)r32b2r3u(r)r3+2b23(r5(u(r)1)(r32b2)2r5u(r)(r3+2b2)2r2(u(r)1)r32b2+r2u(r)r3+2b2)r+r)-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
2E2a2+6m02r2+4(E2m02)r34ELa+2L2+2(E2a2a2m02L2Q)r+2Q2 \, 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
6m02r2+4(E2m02)r3+2L22(L2+Q)r+2Q6 \, 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
8(r3(u(r)1)r32b2r3u(r)r3+2b23(r5(u(r)1)(r32b2)2r5u(r)(r3+2b2)2r2(u(r)1)r32b2+r2u(r)r3+2b2)r+r)m02r+8E2r22(a2+2(r3(u(r)1)r32b2r3u(r)r3+2b2)r+r2)m02+2(m02r2+(EaL)2+Q)(6r5(u(r)1)(r32b2)26r5u(r)(r3+2b2)26(3r7(u(r)1)(r32b2)33r7u(r)(r3+2b2)34r4(u(r)1)(r32b2)2+4r4u(r)(r3+2b2)2+r(u(r)1)r32b2ru(r)r3+2b2)r6r2(u(r)1)r32b2+6r2u(r)r3+2b21)+4((a2+r2)ELa)E-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
2E2a22a2m02+12m02r+12(E2m02)r22L22Q2 \, 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
12m02r+12(E2m02)r22L22Q12 \, 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
2a2e2+12(e21)r22a22l22q+12r2 \, 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
12(e21)r22l22q+12r12 \, {\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})
NaNNaN
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)
0.00.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})
1.62683109949525+1.36936138033619i1.62683109949525 + 1.36936138033619i

Frequency of the orbits at the ISCO

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

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

reset('m0')
var('m0')
m0m_{0}
R.subs({b:0},{m0:1}).simplify_full() # Kerr limit
(E21)r4Qa2+((E21)a2L2Q)r2+2r3+2(E2a22ELa+L2+Q)r{\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
E2r4(L2+Q)r2+2(L2+Q)rE^{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})
(r20.990000000000000)20.420000000000000(r3(u(r)1)r32r3u(r)r3+2)r0.210000000000000r20.170100000000000{\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})
0.240000000000000-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')