¹⁹⁹ views

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.

In [1]:

version()

Out[1]:

'SageMath version 8.9, Release Date: 2019-09-29'

In [2]:

%display latex

In [3]:

from sage.manifolds.utilities import simplify_chain_real as simpl

Spacetime

In [4]:

a0 = 0.7
b0 = 1

In [5]:

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

Out[5]:

In [6]:

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)

Out[6]:

In [7]:

m(0,r)  # Kerr limit

Out[7]:

In [8]:

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

Out[8]:

Derivatives of the mass function

In [9]:

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

Out[9]:

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

In [10]:

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

In [11]:

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

Out[11]:

In [12]:

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

Out[12]:

In [13]:

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

Out[13]:

In [14]:

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

Out[14]:

Horizons

In [15]:

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

Out[15]:

In [16]:

#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

In [17]:

#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)

In [18]:

var('m0','b')

Out[18]:

In [19]:

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

In [20]:

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

Out[20]:

In [21]:

var('e l q');

In [22]:

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

In [23]:

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

Out[23]:

In [24]:

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))

In [25]:

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))

In [26]:

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

Out[26]:

In [27]:

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

Out[27]:

In [28]:

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

Out[28]:

Conditions 1 & 2

In [29]:

reset('M')

In [30]:

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

In [31]:

Condition1 = M/r - Md
Condition1

Out[31]:

In [32]:

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

Out[32]:

In [33]:

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))))

In [34]:

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))))

In [35]:

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

Out[35]:

In [36]:

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

Out[36]:

In [37]:

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)

In [38]:

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)

Out[38]:

verbose 0 (3698: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 74 points.
verbose 0 (3698: plot.py, generate_plot_points) Last error message: 'Unable to compute f(1.53702130825)'
verbose 0 (3698: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 75 points.
verbose 0 (3698: plot.py, generate_plot_points) Last error message: 'Unable to compute f(1.55952294993)'

In [0]:

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

In [0]:

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

In [0]:

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$

In [0]:

var('E L Q m0')

In [0]:

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))

In [0]:

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

In [0]:

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

In [0]:

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

In [0]:

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

In [0]:

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

In [0]:

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

In [0]:

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

In [0]:

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

In [0]:

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

In [0]:

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

In [0]:

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

In [0]:

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

Conditions 1, 2 & 3

In [0]:

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})

In [0]:

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()

In [0]:

In [0]:

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

In [0]:

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

In [0]:

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$',''])

In [0]:

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)

In [0]:

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

In [0]:

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

In [0]:

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

Frequency of the orbits at the ISCO

In [0]:

%display latex

In [0]:

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

In [0]:

%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

In [0]:

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

In [0]:

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

In [0]:

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

In [0]:

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

In [0]:

reset('m0')

In [0]:

var('m0')

In [0]:

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

In [0]:

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

In [0]:

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)

In [0]:

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

In [0]:

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

In [0]:

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

In [0]:

#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)

In [0]:

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

Study of stable circular orbits in the equatorial plane

Spacetime

Derivatives of the mass function

Horizons

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

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

Conditions 1 & 2

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

Defining $R$

Conditions 1, 2 & 3

Frequency of the orbits at the ISCO

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

Product

Resources

Company

Study of stable circular orbits in the equatorial plane

Spacetime

Derivatives of the mass function

Horizons

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

EEE, LLL and Ω\OmegaΩ of a test particle in the equatorial plane

Conditions 1 & 2

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

Defining RRR

Conditions 1, 2 & 3

Frequency of the orbits at the ISCO

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

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

Defining $R$

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