| Download

SageMath notebooks associated to the Black Hole Lectures (https://luth.obspm.fr/~luthier/gourgoulhon/bh16)

Project: BHLectures

Path: BHLectures/sage / ges_null_geod.ipynb

Views: ²⁰¹¹⁰

Kernel: SageMath 9.1.beta0

Null geodesics in Schwarzschild spacetime

This Jupyter/SageMath notebook is relative to the lectures Geometry and physics of black holes.

Click here to download the notebook file (ipynb format). To run it, you must start SageMath with sage -n jupyter.

A version of Sage at least equal to 9.0 is required to run this notebook:

In [1]:

version()

'SageMath version 9.1.beta0, Release Date: 2020-01-10'

In [2]:

%display latex

Schwarzschild metric

In [3]:

M = Manifold(4, 'M', structure='Lorentzian')
X.<t, r, th, ph> = M.chart(r't r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):periodic:\phi')
X

In [4]:

g = M.metric()
m = 1
g[0,0] = -(1-2*m/r)
g[1,1] = 1/(1-2*m/r)
g[2,2] = r^2
g[3,3] = (r*sin(th))^2
g.display()

In [5]:

Xcart.<t, x, y, z> = M.chart()
X_to_Xcart = X.transition_map(Xcart, [t, r*sin(th)*cos(ph), r*sin(th)*sin(ph), r*cos(th)])
X_to_Xcart.display()

In [6]:

M.top_charts()

In [7]:

M._top_charts = [X]
M.top_charts()

In [8]:

M.identity_map().coord_functions(X, Xcart)

Some selected null geodesics

In [9]:

def initial_vector(r0, b, ph0=0, E=1, inward=False):
    t0, th0 = 0, pi/2
    L = b*E
    vt0 = 1/(1-2*m/r0)
    vr0 = sqrt(E^2 - L^2/r0^2*(1-2*m/r0))
    if inward:
        vr0 = - vr0
    vth0 = 0
    vph0 = L / r0^2
    p0 = M((t0, r0, th0, ph0), name='p_0')
    return M.tangent_space(p0)((vt0, vr0, vth0, vph0), name='v_0')

The photon orbit

In [10]:

r0 = 3
ph0 = 0
b = 3*sqrt(3)
E = 1
print("b = {} m".format(float(b)))

b = 5.196152422706632 m

In [11]:

v0 = initial_vector(r0, b)
v0.display()

In [12]:

p0 = v0.parent().base_point()
g.at(p0)(v0, v0)

In [13]:

s = var('s')
geod = M.integrated_geodesic(g, (s, 0, 100), v0)
geod

In [14]:

sol = geod.solve(step=0.01, method="ode_int")       # numerical integration
interp = geod.interpolate()              # interpolation of the solution for the plot

In [15]:

graph = circle((0, 0), 2*m, edgecolor='black', fill=True, facecolor='grey', alpha=0.5)
graph += geod.plot_integrated(chart=Xcart, ambient_coords=(x,y), plot_points=500, 
                              aspect_ratio=1, color='green')        
graph

Other null geodesics

In [16]:

U(r) = (1 - 2/r)/r^2

In [17]:

U1 = 0.02
b = 1/sqrt(U1)
print("b = {} m".format(b))
r0 = 10
v0 = initial_vector(r0, b, inward=True)
v0.display()

b = 7.07106781186548 m

In [18]:

geod = M.integrated_geodesic(g, (s, 0, 50), v0)
sol = geod.solve(step=0.01, method="ode_int") 
interp = geod.interpolate()   
graph += geod.plot_integrated(chart=Xcart, ambient_coords=(x,y), plot_points=500, 
                              aspect_ratio=1, color='red', thickness=1.5,
                              display_tangent=True, color_tangent='red', 
                              plot_points_tangent=4, scale=1) 
graph

In [19]:

U1 = 0.04
b = 1/sqrt(U1)
print("b = {} m".format(b))
r0 = 10
v0 = initial_vector(r0, b, inward=True)
v0.display()

b = 5.00000000000000 m

In [20]:

p0 = v0.parent().base_point()
g.at(p0)(v0, v0)

In [21]:

geod = M.integrated_geodesic(g, (s, 0, 13.25), v0)
sol = geod.solve(step=0.001, method="ode_int") 
interp = geod.interpolate()   
graph += geod.plot_integrated(chart=Xcart, ambient_coords=(x,y), plot_points=500, 
                              aspect_ratio=1, color='brown', thickness=1.5,
                              display_tangent=True, color_tangent='brown', 
                              plot_points_tangent=3, scale=0.2) 
graph

In [22]:

show(graph, xmin=-10)

In [23]:

U1 = 0.046
b = 1/sqrt(U1)
print("b = {} m".format(b))
r0 = 2.1
v0 = initial_vector(r0, b)
v0.display()

b = 4.66252404120157 m

In [24]:

p0 = v0.parent().base_point()
g.at(p0)(v0, v0)

In [25]:

geod = M.integrated_geodesic(g, (s, 0, 20), v0)
sol = geod.solve(step=0.001, method="ode_int") 
interp = geod.interpolate()   
graph += geod.plot_integrated(chart=Xcart, ambient_coords=(x,y), plot_points=500, 
                              aspect_ratio=1, color='orange', thickness=1.5,
                              display_tangent=True, color_tangent='orange', 
                              plot_points_tangent=3, scale=0.1)       
graph

In [26]:

show(graph, xmin=-10)

In [27]:

U1 = 0.035
b = 1/sqrt(U1)
print("b = {} m".format(b))
r0 = 2.4
v0 = initial_vector(r0, b, ph0=-pi/2)
v0.display()

b = 5.34522483824849 m

In [28]:

geod = M.integrated_geodesic(g, (s, 0, 3.5), v0)
sol = geod.solve(step=0.001, method="ode_int") 
interp = geod.interpolate()   
graph += geod.plot_integrated(chart=Xcart, ambient_coords=(x,y), plot_points=500, 
                              aspect_ratio=1, color='grey', thickness=1.5,
                              display_tangent=True, color_tangent='grey', 
                              plot_points_tangent=3, scale=0.1) 
graph

In [29]:

graph += text('1', (5, 4.5), color='red' ) \
         + text('2', (3, 3.8), color='brown' ) \
         + text('3', (-6, 3.5), color='orange' ) \
         + text('4', (0.7, -2.7), color='grey' )

In [30]:

show(graph, xmin=-10, ymin=-4, axes_labels=[r'$x/m$', r'$y/m$'])

In [31]:

graph.save("ges_null_geod.pdf", xmin=-10, ymin=-4, axes_labels=[r'$x/m$', r'$y/m$'])

Geodesics arising from a source at large distance

In [32]:

r0 = 100  # distance of the source

In [33]:

b_sel = [12., 8., 6., 5.355, 5.23, 5.2025, 5.19643, 5.196155]
#b_sel = [15., 10., 7., 6., 5.5, 5.23, 5.22, 5.2025, 5.198, 5.197, 5.1965, 5.19632]
b_pec = [5.355, 5.2025, 5.19632, 5.196155]

In [34]:

#b_sel = [5.19643]
for b in b_sel:
    print('b = {:.6f} m'.format(float(b)))
    graph = circle((0, 0), 2*m, edgecolor='black', fill=True, facecolor='grey', alpha=0.5)
    graph += circle((0, 0), 3*m, color='lightgreen', thickness=1., linestyle='--', zorder=1)
    graph += text(r'$b = {:.6f} \, m$'.format(float(b)), (-8, 12), fontsize=14, 
                  color='black')
    ph0 = asin(b/r0)
    v0 = initial_vector(r0, b, ph0=ph0, inward=True)
    geod = M.integrated_geodesic(g, (s, 0, 150), v0)
    sol = geod.solve(step=0.001, method="odeint", rtol=1.e-12, atol=1.e-12) 
    interp = geod.interpolate()
    first_point = geod(0)
    print("s=0:   {}".format(first_point.coord()))
    last_point = geod(150)
    print("s=150: {}".format(last_point.coord()))
    graph += geod.plot_integrated(chart=Xcart, ambient_coords=(x,y), plot_points=500, 
                                  aspect_ratio=1, color='green', thickness=1.5, zorder=2,
                                  display_tangent=True, color_tangent='green', 
                                  plot_points_tangent=12, scale=0.1) 
    show(graph, xmin=-12, xmax=12, ymin=-12, ymax=12, axes_labels=[r'$x/m$', r'$y/m$'])
    filename = 'ges_null_b_{:.6f}'.format(float(b))
    filename = filename.replace('.', '_') + '.pdf'
    graph.save(filename, xmin=-12, xmax=12, ymin=-12, ymax=12, axes_labels=[r'$x/m$', r'$y/m$'])

b = 12.000000 m
s=0:   (0.0, 100.0, 1.5707963267948966, 0.12028988239478806)
s=150: (161.7386305994338, 49.599207321700575, 1.5707963267948966, 3.3499140843834243)

b = 8.000000 m
s=0:   (0.0, 100.0, 1.5707963267948966, 0.08008558003365901)
s=150: (165.0262559276518, 47.87567237981162, 1.5707963267948966, 3.832461957010324)

b = 6.000000 m
s=0:   (0.0, 100.0, 1.5707963267948966, 0.06003605844527842)
s=150: (169.51340790400954, 46.21354032584246, 1.5707963267948966, 4.730793903988559)

b = 5.355000 m
s=0:   (0.0, 100.0, 1.5707963267948966, 0.05357562643499787)
s=150: (175.13467942245063, 43.837643428812676, 1.5707963267948966, 6.172293011040655)

b = 5.230000 m
s=0:   (0.0, 100.0, 1.5707963267948966, 0.052323872006442895)
s=150: (180.39625091484754, 41.32000814165003, 1.5707963267948966, 7.665165130714838)

b = 5.202500 m
s=0:   (0.0, 100.0, 1.5707963267948966, 0.052048497112969744)
s=150: (186.05222271416335, 38.48872884422706, 1.5707963267948966, 9.317135314376204)

b = 5.196430 m
s=0:   (0.0, 100.0, 1.5707963267948966, 0.05198771489670833)
s=150: (196.57550821738923, 33.13574656928606, 1.5707963267948966, 12.42140053051524)

b = 5.196155 m
s=0:   (0.0, 100.0, 1.5707963267948966, 0.05198496117647258)
s=150: (212.1623795944796, 25.165821003889214, 1.5707963267948966, 17.043633460770014)

Check of the law $b - b_{\rm crit} \sim 648(7\sqrt{3} - 12) m \, \mathrm{e}^{-\Delta\varphi}$

In [35]:

b_c = n(3*sqrt(3))
a = n(648*(7*sqrt(3) - 12))

In [36]:

b = 5.230000
Dphi = 7.665165130714838 - 0.052323872006442895
b - b_c, a*exp(-Dphi), abs((b - b_c)/(a*exp(-Dphi)) - 1)

In [37]:

b = 5.202500
Dphi = 9.317135314376204 - 0.052048497112969744
b - b_c, a*exp(-Dphi), abs((b - b_c)/(a*exp(-Dphi)) - 1)

In [38]:

b = 5.196430
Dphi = 12.42140053051524 - 0.05198771489670833
b - b_c, a*exp(-Dphi), abs((b - b_c)/(a*exp(-Dphi)) - 1)

In [39]:

b = 5.196155
Dphi = 17.043633460770014 - 0.05198496117647258
b - b_c, a*exp(-Dphi), abs((b - b_c)/(a*exp(-Dphi)) - 1)

In [ ]: