Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News AboutSign UpSign In
| Download
Views: 133
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')
θ\renewcommand{\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)
−r3(u(r)−1)r3−2 b2+r3u(r)r3+2 b2\renewcommand{\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
1\renewcommand{\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))
Image in a Jupyter notebook

Derivatives of the mass function

md(b, r) = diff(m(b, r), r) md(b, r)
3 r5(u(r)−1)(r3−2 b2)2−3 r5u(r)(r3+2 b2)2+r3δ(r)r3+2 b2−r3δ(r)r3−2 b2−3 r2(u(r)−1)r3−2 b2+3 r2u(r)r3+2 b2\renewcommand{\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)
3 r5(u(r)−1)(r3−2 b2)2−3 r5u(r)(r3+2 b2)2−3 r2(u(r)−1)r3−2 b2+3 r2u(r)r3+2 b2\renewcommand{\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))
Image in a Jupyter notebook
mdd(b, r) = remove_dirac(diff(md(b, r), r)) mdd(b, r)
−18 r7(u(r)−1)(r3−2 b2)3+18 r7u(r)(r3+2 b2)3+24 r4(u(r)−1)(r3−2 b2)2−24 r4u(r)(r3+2 b2)2−6 r(u(r)−1)r3−2 b2+6 ru(r)r3+2 b2\renewcommand{\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))
Image in a Jupyter notebook

Horizons

D=r^2+a^2-2*m(b,r)*r plot(D.subs({a:a0},{b:b0}),(r,-5,5))
Image in a Jupyter notebook
#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)\renewcommand{\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)
−3 r5(r3+2 b2)2+3 r2r3+2 b2\renewcommand{\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}))
ar∂∂rM(r)+−r∂∂rM(r)+M(r)r32−aM(r)a2r∂∂rM(r)+r3−a2M(r)\renewcommand{\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()
a2r∂∂rM(r)+r3+2 −r∂∂rM(r)+M(r)arM(r)−(a2+2 r2)M(r)−a2r3∂∂rM(r)2+r5−3 a2rM(r)2+2 ((a3+3 ar2)M(r)−(a3r+ar3)∂∂rM(r))−r∂∂rM(r)+M(r)r−3 (a2r2+r4)M(r)+(3 a2r3+r5+4 a2r2M(r))∂∂rM(r)r\renewcommand{\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()
−(a3+3 ar2)M(r)−(a3r+ar3)∂∂rM(r)−(a2r2+r4+2 a2rM(r))−r∂∂rM(r)−M(r)r(a2r∂∂rM(r)+r3−a2M(r))−a2r4∂∂rM(r)2−r6+3 a2r2M(r)2+3 (a2r3+r5)M(r)−(3 a2r4+r6+4 a2r3M(r))∂∂rM(r)−2 ((a3r2+3 ar4)M(r)−(a3r3+ar5)∂∂rM(r))−r∂∂rM(r)−M(r)ra4r2∂∂rM(r)2+r6−2 a2r3M(r)+a4M(r)2+2 (a2r4−a4rM(r))∂∂rM(r)\renewcommand{\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()
−(a3+3 ar2)M(r)−(a3r+ar3)∂∂rM(r)+(a2r2+r4+2 a2rM(r))−r∂∂rM(r)−M(r)r(a2r∂∂rM(r)+r3−a2M(r))−a2r4∂∂rM(r)2−r6+3 a2r2M(r)2+3 (a2r3+r5)M(r)−(3 a2r4+r6+4 a2r3M(r))∂∂rM(r)+2 ((a3r2+3 ar4)M(r)−(a3r3+ar5)∂∂rM(r))−r∂∂rM(r)−M(r)ra4r2∂∂rM(r)2+r6−2 a2r3M(r)+a4M(r)2+2 (a2r4−a4rM(r))∂∂rM(r)\renewcommand{\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
−3 r5(u(r)−1)(r3−2 b2)2+3 r5u(r)(r3+2 b2)2+3 r2(u(r)−1)r3−2 b2−3 r2u(r)r3+2 b2−r3(u(r)−1)r3−2 b2−r3u(r)r3+2 b2r\renewcommand{\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
1r\renewcommand{\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
r4+2 a3−3 a2r+6 ar2−3 r3−3 a2\renewcommand{\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
r4−3 r3\renewcommand{\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: ''
Image in a Jupyter notebook
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.52502568341426−0.0409208183571186i\renewcommand{\Bold}[1]{\mathbf{#1}}8.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)\renewcommand{\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
((a2+r2)E−La)2−(m02r2+(Ea−L)2+Q)(a2+2 (r3(u(r)−1)r3−2 b2−r3u(r)r3+2 b2)r+r2)\renewcommand{\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
2 m02r3+(E2−m02)r4−Qa2+(E2a2−a2m02−L2−Q)r2+2 (E2a2−2 ELa+L2+Q)r\renewcommand{\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
2 m02r3+(E2−m02)r4−(L2+Q)r2+2 (L2+Q)r\renewcommand{\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
−2 (a2+2 (r3(u(r)−1)r3−2 b2−r3u(r)r3+2 b2)r+r2)m02r+4 ((a2+r2)E−La)Er−2 (m02r2+(Ea−L)2+Q)(r3(u(r)−1)r3−2 b2−r3u(r)r3+2 b2−3 (r5(u(r)−1)(r3−2 b2)2−r5u(r)(r3+2 b2)2−r2(u(r)−1)r3−2 b2+r2u(r)r3+2 b2)r+r)\renewcommand{\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
2 E2a2+6 m02r2+4 (E2−m02)r3−4 ELa+2 L2+2 (E2a2−a2m02−L2−Q)r+2 Q\renewcommand{\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
6 m02r2+4 (E2−m02)r3+2 L2−2 (L2+Q)r+2 Q\renewcommand{\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
−8 (r3(u(r)−1)r3−2 b2−r3u(r)r3+2 b2−3 (r5(u(r)−1)(r3−2 b2)2−r5u(r)(r3+2 b2)2−r2(u(r)−1)r3−2 b2+r2u(r)r3+2 b2)r+r)m02r+8 E2r2−2 (a2+2 (r3(u(r)−1)r3−2 b2−r3u(r)r3+2 b2)r+r2)m02+2 (m02r2+(Ea−L)2+Q)(6 r5(u(r)−1)(r3−2 b2)2−6 r5u(r)(r3+2 b2)2−6 (3 r7(u(r)−1)(r3−2 b2)3−3 r7u(r)(r3+2 b2)3−4 r4(u(r)−1)(r3−2 b2)2+4 r4u(r)(r3+2 b2)2+r(u(r)−1)r3−2 b2−ru(r)r3+2 b2)r−6 r2(u(r)−1)r3−2 b2+6 r2u(r)r3+2 b2−1)+4 ((a2+r2)E−La)E\renewcommand{\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
2 E2a2−2 a2m02+12 m02r+12 (E2−m02)r2−2 L2−2 Q\renewcommand{\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
12 m02r+12 (E2−m02)r2−2 L2−2 Q\renewcommand{\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
2 a2e2+12 (e2−1)r2−2 a2−2 l2−2 q+12 r\renewcommand{\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
12 (e2−1)r2−2 l2−2 q+12 r\renewcommand{\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})
NaN\renewcommand{\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'
Image in a Jupyter notebook
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: ''
Image in a Jupyter notebook
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: ''
Image in a Jupyter notebookImage in a Jupyter notebook
find_root(Condition2_pos.subs({a:a0},{b:b0},{th:pi/2}),1.001,2)
0.0\renewcommand{\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'
Image in a Jupyter notebook
Condition3_pos.subs({a:a0},{b:b0},{th:pi/2},{r:1.04})
1.62683109949525+1.36936138033619i\renewcommand{\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')
θ\renewcommand{\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
(a,b,r) ↦ 0.00107187558536188\renewcommand{\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
(a,b,r) ↦ −0.00197809424213708\renewcommand{\Bold}[1]{\mathbf{#1}}\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')
m0\renewcommand{\Bold}[1]{\mathbf{#1}}m_{0}
R.subs({b:0},{m0:1}).simplify_full() # Kerr limit
(E2−1)r4−Qa2+((E2−1)a2−L2−Q)r2+2 r3+2 (E2a2−2 ELa+L2+Q)r\renewcommand{\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
E2r4−(L2+Q)r2+2 (L2+Q)r\renewcommand{\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)
Image in a Jupyter notebook
R.subs({m0:0},{a:a2},{b:b2},{E:e2},{L:l2},{Q:q2})
(r2−0.990000000000000)2−0.420000000000000 (r3(u(r)−1)r3−2−r3u(r)r3+2)r−0.210000000000000 r2−0.170100000000000\renewcommand{\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})
−0.240000000000000\renewcommand{\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)
Image in a Jupyter notebook
#plot2.save('R_photons_dark_zone_b=0.pdf')