ges_null_images.ipynb

Kernel: SageMath 9.1.beta0

Multiple images 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.

version()

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

%display latex

Schwarzschild metric

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

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

For plots:

E.<x, y, z> = EuclideanSpace()
phi = M.diff_map(E, [r*sin(th)*cos(ph), r*sin(th)*sin(ph), r*cos(th)])
phi.display()

Bunch of geodesics

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

r0 = 100  # distance of the source
s = var('s')  # affine parameter along a geodesic

b_max = 12.
b_min = - b_max
nb = 31
db = (b_max - b_min) / (nb - 1)
b_sel = [(b_min + i*db, 'green') for i in range(nb)]
#
b_min = 5.0
b_max = 5.4
nb = 3
db = (b_max - b_min) / (nb - 1)
b_sel += [(b_min + i*db, 'olive') for i in range(nb)]
#
b_min = -5.4
b_max = -5.0
nb = 3
db = (b_max - b_min) / (nb - 1)
b_sel += [(b_min + i*db, 'olive') for i in range(nb)]
b_sel

Computation:

geodesics = {}
for b, color in b_sel:
    print('b = {:.6f} m'.format(float(b)))
    ph0 = asin(b/r0)
    v0 = initial_vector(r0, b, ph0=ph0, inward=True)
    geod = M.integrated_geodesic(g, (s, 0, 200), v0, across_charts=True)
    sol = geod.solve_across_charts(step=0.05) 
    interp = geod.interpolate()   
    geodesics[b] = geod

b = -12.000000 m
b = -11.200000 m
b = -10.400000 m
b = -9.600000 m
b = -8.800000 m
b = -8.000000 m
b = -7.200000 m
b = -6.400000 m
b = -5.600000 m
b = -4.800000 m
b = -4.000000 m
b = -3.200000 m
b = -2.400000 m
b = -1.600000 m
b = -0.800000 m
b = 0.000000 m
b = 0.800000 m
b = 1.600000 m
b = 2.400000 m
b = 3.200000 m
b = 4.000000 m
b = 4.800000 m
b = 5.600000 m
b = 6.400000 m
b = 7.200000 m
b = 8.000000 m
b = 8.800000 m
b = 9.600000 m
b = 10.400000 m
b = 11.200000 m
b = 12.000000 m
b = 5.000000 m
b = 5.200000 m
b = 5.400000 m
b = -5.400000 m
b = -5.200000 m
b = -5.000000 m

Images of a generic source $S$

graph = circle((0, 0), 2*m, edgecolor='black', thickness=2, fill=True, facecolor='grey', 
               alpha=0.5)
graph += circle((0, 0), 3*m, color='lightgreen', thickness=1.5, linestyle='--', zorder=1)
b_spec = [7.2, 5.2, -6.4]  # geodesics through S
for b, color in b_sel:
    thickness = 2 if any(abs(b - bs)<1e-14 for bs in b_spec) else 1 
    graph += geodesics[b].plot_integrated(mapping=phi, ambient_coords=(x,y), plot_points=500,
                                         across_charts=True, aspect_ratio=1, color=[color], 
                                         thickness=thickness, zorder=2) 
 
show(graph, xmin=-12, xmax=20, ymin=-12, ymax=12, axes_labels=[r'$x/m$', r'$y/m$'])

graphp = graph + circle((-6.45, 2.05), 0.25, fill=True, color='red') \
         + text(r"$S$", (-6.2, 3.), color='red', fontsize=20) \
         + circle((20.32, 5.2), 0.25, thickness=2, color='red') \
         + text("3", (19.5, 5.7), color='red', fontsize=16) \
         + circle((20.32, 7.2), 0.25, thickness=2,color='red') \
         + text("1", (19.5, 7.7), color='red', fontsize=16) \
         + circle((20.32, -6.4), 0.25, thickness=2, color='red') \
         + text("2", (19.5, -5.9), color='red', fontsize=16)
show(graphp, xmin=-12, xmax=20, ymin=-12, ymax=12, axes_labels=[r'$x/m$', r'$y/m$'],
     frame=True, axes=False, figsize=10)

graphp.save('ges_mult_images.pdf', xmin=-12, xmax=20, ymin=-12, ymax=12, 
            axes_labels=[r'$x/m$', r'$y/m$'], frame=True, axes=False, figsize=10)

Image of an aligned source (source on the $x$-axis)

graph = circle((0, 0), 2*m, edgecolor='black', thickness=2, fill=True, facecolor='grey', 
               alpha=0.5)
graph += circle((0, 0), 3*m, color='lightgreen', thickness=1.5, linestyle='--', zorder=1)
b_spec = [8, -8, 5.2, -5.2]  # geodesics through A
for b, color in b_sel:
    thickness = 2 if any(abs(b - bs)<1e-14 for bs in b_spec) else 1 
    graph += geodesics[b].plot_integrated(mapping=phi, ambient_coords=(x,y), plot_points=500,
                                         across_charts=True, aspect_ratio=1, color=[color], 
                                         thickness=thickness, zorder=2) 
 
show(graph, xmin=-12, xmax=20, ymin=-12, ymax=12, axes_labels=[r'$x/m$', r'$y/m$'])

graphp = graph + circle((-10.35, 0.), 0.25, fill=True, color='red') \
         + text(r"$A$", (-10.3, 1.1), color='red', fontsize=20) \
         + circle((20.32, 5.2), 0.25, thickness=2, color='red') \
         + text("3", (19.5, 5.7), color='red', fontsize=16) \
         + circle((20.32, -5.2), 0.25, thickness=2, color='red') \
         + text("4", (19.5, -4.7), color='red', fontsize=16) \
         + circle((20.32, 8.0), 0.25, thickness=2,color='red') \
         + text("1", (19.5, 8.5), color='red', fontsize=16) \
         + circle((20.32, -8.0), 0.25, thickness=2, color='red') \
         + text("2", (19.5, -7.5), color='red', fontsize=16)
show(graphp, xmin=-12, xmax=20, ymin=-12, ymax=12, axes_labels=[r'$x/m$', r'$y/m$'],
     frame=True, axes=False, figsize=10)

graphp.save('ges_images_aligned.pdf', xmin=-12, xmax=20, ymin=-12, ymax=12, 
            axes_labels=[r'$x/m$', r'$y/m$'], frame=True, axes=False, figsize=10)