This Jupyter/SageMath notebook is relative to the lectures Geometry and physics of black holes.
version()
'SageMath version 9.3.beta7, Release Date: 2021-02-07'
%display latex
We use $m=1$ and denote $r_0$ simply by $r$.
a, r = var('a r')
lsph(a, r) = (r^2*(3 - r) - a^2*(r + 1))/(a*(r -1))
lsph
qsph(a, r) = r^3 / (a^2*(r - 1)^2) * (4*a^2 - r*(r - 3)^2)
qsph
The radii $r_+$ and $r_-$ of the two horizons:
rp(a) = 1 + sqrt(1 - a^2)
rm(a) = 1 - sqrt(1 - a^2)
rph_ss(a) = 1/2 + cos(2/3*asin(a) + 2*pi/3)
rph_ss
rph_s(a) = 4*cos(acos(-a)/3 + 4*pi/3)^2
rph_s
rph_p(a) = 4*cos(acos(-a)/3)^2
rph_p
rph_m(a) = 4*cos(acos(a)/3)^2
rph_m
We add the radius of the marginally stable orbit:
rph_ms(a) = 1 - (1 - a^2)^(1/3)
rph_ms
as well as the radius of outer and inner polar orbits:
rph_pol(a) = 1 + 2*sqrt(1 - a^2/3)*cos(1/3*arccos((1 - a^2)/(1 - a^2/3)^(3/2)))
rph_pol
rph_pol_in(a) = 1 + 2*sqrt(1 - a^2/3)*cos(1/3*arccos((1 - a^2)/(1 - a^2/3)^(3/2)) + 2*pi/3)
rph_pol_in
a0 = 0.95
# a0 = 1
show(LatexExpr(r'r_{\rm ph}^{**} = '), n(rph_ss(a0)))
show(LatexExpr(r'r_{\rm ph}^{\rm ms} = '), n(rph_ms(a0)))
show(LatexExpr(r'r_{\rm ph}^{*} = '), n(rph_s(a0)))
show(LatexExpr(r'r_- = '), n(rm(a0)))
show(LatexExpr(r'r_+ = '), n(rp(a0)))
show(LatexExpr(r'r_{\rm ph}^+ = '), n(rph_p(a0)))
show(LatexExpr(r'r_{\rm ph}^- = '), n(rph_m(a0)))
show(LatexExpr(r'r_{\rm ph}^{\rm pol} = '), n(rph_pol(a0)))
show(LatexExpr(r'r_{\rm ph}^{\rm pol,in} = '), n(rph_pol_in(a0)))
r1 = rph_p(a0)
r2 = rph_m(a0)
#r1 = n(rph_ss(a0))
#r2 = 0
r1, r2
r = var('r')
plot(qsph(a0, r), (r, r1, r2)) + plot(r^3*(4 - r), (r, r1, r2), color='red')
def alpha(a, th_obs, r0):
if a == 1:
ell = - r0^2 + 2*r0 + 1
else:
ell = lsph(a, r0)
return - ell / sin(th_obs)
def Theta(a, th_obs, r0):
if a == 1:
ell = - r0^2 + 2*r0 + 1
q = r0^3 * (4 - r0)
else:
ell = lsph(a, r0)
q = qsph(a, r0)
return q + cos(th_obs)^2 * (a^2 - ell^2/sin(th_obs)^2)
def beta(a, th_obs, r0, eps_theta=1):
return eps_theta * sqrt(Theta(a, th_obs, r0,))
def r0_bounds(a, th_obs, outer=True):
r"""
Return `(r0_min, r0_max)`
"""
if outer:
r1 = n(rph_p(a))
r2 = n(rph_m(a))
r3 = rph_pol(a)
else:
r1 = n(rph_ss(a))
r2 = 0
r3 = rph_pol_in(a)
#
# Computation of rmin:
try:
if a == 1:
th_crit = n(asin(sqrt(3) - 1))
if n(th_obs) < th_crit or n(th_obs) > n(pi) - th_crit:
rmin = find_root(lambda r: Theta(a, th_obs, r), r1, r3)
else:
rmin = 1
else:
rmin = find_root(lambda r: Theta(a, th_obs, r), r1, r3)
except TypeError: # special case th_obs = pi/2
rmin = r1
#
# Computation of rmax:
try:
rmax = find_root(lambda r: Theta(a, th_obs, r), r3, r2)
except TypeError: # special case th_obs = pi/2
rmax = r2
#
return (rmin, rmax)
th_obs = pi/2
#th_obs = n(160/180*pi) # M87 value
# th_obs = 1.e-3
show(LatexExpr(r'\theta_{\mathscr{O}} = '), n(th_obs))
fa = lambda r: alpha(a0, th_obs, r)
fb = lambda r: beta(a0, th_obs, r)
fbm = lambda r: beta(a0, th_obs, r, eps_theta=-1)
plot(Theta(a0, th_obs, r), (r, r1, r2))
rmin, rmax = r0_bounds(a0, th_obs)
rmin - r1, rmax - r2
fa(rmin), fa(rmax)
fb(rmin), fb(rmax)
def shadow_plot(a, th_obs, outer=True, color=None, number_colors=5,
range_col=None, thickness=2, linestyle='-', plot_points=200,
legend='automatic', legend_loc=(1.02, 0.36), fill=True,
fillcolor='grey', draw_NHEKline=True, frame=True, axes=True,
axes_labels='automatic', gridlines=True):
if a==0:
# Case a = 0:
rs = 3*sqrt(3)
if color is None:
color = 'black'
if legend == 'automatic':
legend = None
g = parametric_plot((rs*cos(x), rs*sin(x)), (x, 0, 2*pi),
color=color, thickness=thickness,
linestyle=linestyle, legend_label=legend,
fill=fill, fillcolor=fillcolor, frame=frame,
axes=axes, gridlines=gridlines)
else:
# Case a != 0
rmin, rmax = r0_bounds(a, th_obs, outer=outer)
if rmin > 0:
rmin = 1.00000001*rmin
rmax = 0.99999999*rmax
else:
rmin = 0.9999999*rmin
rmax = 1.0000001*rmax
print("rmin : ", rmin, " rmax : ", rmax)
fa = lambda r: alpha(a, th_obs, r)
fb = lambda r: beta(a, th_obs, r)
fbm = lambda r: beta(a, th_obs, r, eps_theta=-1)
if range_col is None:
range_col = r0_bounds(a, pi/2, outer=outer)
rmin_col, rmax_col = range_col
print("rmin_col : ", rmin_col, " rmax_col : ", rmax_col)
dr = (rmax_col - rmin_col) / number_colors
rm = rmin_col + int((rmin - rmin_col)/dr)*dr
r1s = rmin
r_ranges = []
while rm + dr < rmax:
col = hue((rm - rmin_col)/(rmax_col - rmin_col + 0.1))
r2s = rm + dr
r_ranges.append((r1s, r2s, col))
rm += dr
r1s = r2s
if color is None:
col = hue((rm - rmin_col)/(rmax_col - rmin_col + 0.1))
else:
col = color
r_ranges.append((r1s, rmax, col))
g = Graphics()
legend_label = None # a priori
if a == 1 and draw_NHEKline:
th_crit = asin(sqrt(3) - 1)
if th_obs > th_crit and th_obs < pi - th_crit:
# NHEK line
alpha0 = -2/sin(th_obs)
beta0 = sqrt(3 - cos(th_obs)**2 *(6 + cos(th_obs)**2))/sin(th_obs)
if legend == 'automatic':
legend_label = r"$r_0 = m$"
if color is None:
colNHEK = 'maroon'
else:
colNHEK = color
g += line([(alpha0, beta0), (alpha0, -beta0)], color=colNHEK,
thickness=thickness, linestyle=linestyle,
legend_label=legend_label)
if fill:
g += polygon2d([(fa(rmax), 0), (alpha0, beta0), (alpha0, -beta0)],
color=fillcolor, alpha=0.5)
for rg in r_ranges:
r1s, r2s = rg[0], rg[1]
col = rg[2]
if legend:
if legend == 'automatic':
if draw_NHEKline and abs(r1s - 1) < 1e-5:
legend_label = r"${:.2f}\, m < r_0 \leq {:.2f}\, m$".format(
float(r1s), float(r2s))
else:
legend_label = r"${:.2f}\, m \leq r_0 \leq {:.2f}\, m$".format(
float(r1s), float(r2s))
else:
legend_label = legend
g += parametric_plot((fa, fb), (r1s, r2s), plot_points=plot_points, color=col,
thickness=thickness, linestyle=linestyle,
legend_label=legend_label,
frame=frame, axes=axes, gridlines=gridlines)
g += parametric_plot((fa, fbm), (r1s, r2s), plot_points=plot_points, color=col,
thickness=thickness, linestyle=linestyle)
if fill:
g += parametric_plot((fa, fb), (rmin, rmax), fill=True, fillcolor=fillcolor,
thickness=0)
g += parametric_plot((fa, fbm), (rmin, rmax), fill=True, fillcolor=fillcolor,
thickness=0)
# end of case a != 0
if axes_labels:
if axes_labels == 'automatic':
g.axes_labels([r"$(r_{\mathscr{O}}/m)\; \alpha$",
r"$(r_{\mathscr{O}}/m)\; \beta$"])
else:
g.axes_labels(axes_labels)
g.set_aspect_ratio(1)
if legend:
g.set_legend_options(handlelength=2, loc=legend_loc)
return g
shadow_plot(0, pi/2)
g1 = shadow_plot(a0, pi/2)
g1.set_axes_range(-4, 8, -6, 6)
g1.save('gik_shadow_a95_th90.pdf')
g1
rmin : 1.38628054232578 rmax : 3.95534727811920 rmin_col : 1.3862805284629751 rmax_col : 3.95534731767268
g1 = shadow_plot(a0, pi/4, legend=None)
g1.set_axes_range(-4.5, 7.5, -6, 6)
g1.save('gik_shadow_a95_th45.pdf')
g1
rmin : 1.54632082701639 rmax : 3.53769021424612 rmin_col : 1.3862805284629751 rmax_col : 3.95534731767268
g1 = shadow_plot(a0, pi/6, legend=None)
g1.set_axes_range(-4.5, 7.5, -6, 6)
g1.save('gik_shadow_a95_th30.pdf')
g1
rmin : 1.76846188276999 rmax : 3.23697774665377 rmin_col : 1.3862805284629751 rmax_col : 3.95534731767268
g1 = shadow_plot(a0, 1e-3, legend=None)
g1.set_axes_range(-6, 6, -6, 6)
g1.save('gik_shadow_a95_th00.pdf')
g1
rmin : 2.49118715973461 rmax : 2.49420141399888 rmin_col : 1.3862805284629751 rmax_col : 3.95534731767268
Shadow radius for the on-axis observer:
def R_shad(a):
r0 = rph_pol(a)
return 2*sqrt(r0*(r0^2 - a^2))/(r0 - 1)
R_shad(a0)
g1 = plot(R_shad, (0, 1), thickness=2,
axes_labels=[r'$a/m$', r'$(r_{\mathscr{O}}/m) \; R_{\rm shad}$'],
frame=True, axes=False, gridlines=True)
g1.save('gik_shadow_radius_axis.pdf')
g1
Extremal values:
R_a_0 = n(3*sqrt(3))
R_a_1 = n(2*(sqrt(2)+1))
R_a_0, R_a_1
Relative amplitude:
1 - R_a_1 / R_a_0
The partial shadow boundary:
g1 = shadow_plot(1, pi/2, fill=False, draw_NHEKline=False)
g1.set_axes_range(-4, 8, -6, 6)
g1.save('gik_shadow_a1_th90_part.pdf')
g1
rmin : 1.00000001000000 rmax : 3.99999996000000 rmin_col : 1 rmax_col : 4.0
With the NHEK line:
g1 = shadow_plot(1, pi/2)
g1.set_axes_range(-4, 8, -6, 6)
g1.save('gik_shadow_a1_th90.pdf')
g1
rmin : 1.00000001000000 rmax : 3.99999996000000 rmin_col : 1 rmax_col : 4.0
As soon as $a\neq m$, the shadow boundary is a smooth closed curve:
shadow_plot(0.9999, pi/2, fill=False)
rmin : 1.01637428079421 rmax : 3.99991107028894 rmin_col : 1.01637427063047 rmax_col : 3.9999111102880467
th = var('th', latex_name=r'\theta_{\mathscr{O}}')
r0 = var('r0', latex_name=r'r_0')
alpha(1, th, r0)
beta(1, th, r0).simplify_full()
Special values for $\theta_{\mathscr{O}} = \pi/2$:
alpha(1, pi/2, r0)
beta(1, pi/2, r0).simplify_full()
Special values for $r_0 = m$:
alpha(1, th, 1)
beta(1, th, 1).simplify_full()
Critical inclination for the existence of the NHEK line:
th_crit = asin(sqrt(3) - 1)
th_crit, n(th_crit), n(th_crit/pi*180)
(cos(th_crit)^4 + 6*cos(th_crit)^2 - 3).simplify_full()
beta(1, th_crit, 1).simplify_full()
g1 = shadow_plot(1, th_crit, legend=False)
g1.set_axes_range(-4, 8, -6, 6)
g1.save('gik_shadow_a1_th_crit.pdf')
g1
rmin : 1.00000001000000 rmax : 3.59326048984047 rmin_col : 1 rmax_col : 4.0
shadow_plot(1, 1.1*th_crit)
rmin : 1.00000001000000 rmax : 3.67516782408563 rmin_col : 1 rmax_col : 4.0
shadow_plot(1, 0.9*th_crit)
rmin : 1.13415225868399 rmax : 3.50327928110302 rmin_col : 1 rmax_col : 4.0
g1 = shadow_plot(1, pi/3, legend=False)
g1.set_axes_range(-4, 8, -6, 6)
g1.save('gik_shadow_a1_th60.pdf')
g1
rmin : 1.00000001000000 rmax : 3.79787701838380 rmin_col : 1 rmax_col : 4.0
n(sqrt(3))
cardioid = parametric_plot([4*(1 + cos(x))*cos(x) - 1,
4*(1 + cos(x))*sin(x)], (x, 0, 2*pi),
linestyle=":", thickness=2.5)
cardioid
g1 = shadow_plot(1, pi/2, fill=False, number_colors=1, color='red',
legend=False) + cardioid
g1.axes_labels([r"$\bar{\alpha}$", r"$\bar{\beta}$"])
g1.set_axes_range(-3, 8, -6, 6)
g1.save('gik_shadow_cardioid.pdf')
g1
rmin : 1.00000001000000 rmax : 3.99999996000000 rmin_col : 1 rmax_col : 4.0
Generic form of $\mathcal{R}(r)$:
l = var('l', latex_name=r'\ell')
q = var('q')
R(a, l, q, r) = r^4 + (a^2 - l^2 - q)*r^2 + 2*(q + (l - a)^2)*r - a^2*q
R
Specialization to $a=m$ and $\ell = 2m$:
R1(q, r) = R(1, 2, q, r)
R1
R1(q, r).factor()
R1(3, r).factor()
g0 = plot(R1(1, r), (r, 0, 2), color='blue', thickness=2,
legend_label=r'$q = m^2$',
axes_labels=[r'$r/m$', r'$\mathcal{R}(r)/m^4$'],
frame=True, gridlines=True)
g0 += plot(R1(3, r), (r, 0, 2), color='red', thickness=2,
legend_label=r'$q = 3 m^2$')
g0 += plot(R1(5, r), (r, 0, 2), color='violet', thickness=2,
legend_label=r'$q = 5 m^2$')
g0 += line([(1, -2), (1, 2)], thickness=2, color='black')
g0 += point((1, 0), size=60, color='red', zorder=100)
g0 += point((sqrt(2) - 1, 0), size=60, color='blue', zorder=100)
g0 += point((sqrt(6) - 1, 0), size=60, color='violet', zorder=100)
g0 += line([(1, -0.9), (3, -0.9)], color='red', thickness=2,
linestyle=':')
g0 += line([(1, -1), (3, -1)], color='blue', thickness=2,
linestyle=':')
g0 += line([(sqrt(6) - 1, -0.8), (3, -0.8)], color='violet',
thickness=2, linestyle=':')
g0.set_legend_options(handlelength=2)
g0.set_axes_range(ymin=-1, ymax=1.5, xmax=2)
g0.save('gik_R_a_1_l_2m.pdf')
g0
$\theta_{\mathcal{O}} = \frac{\pi}{2}$
th_obs = pi/2
title = r'$\theta_{\mathscr{O}}={\pi}/{2}$'
g1 = shadow_plot(0, th_obs, number_colors=1, color='purple',
legend=r'$a=0$', fill=False)
g2 = shadow_plot(0.5, th_obs, number_colors=1, color='peru',
legend=r'$a=0.5\, m$', fill=False, plot_points=500)
g3 = shadow_plot(1, th_obs, number_colors=1, color='red',
legend=r'$a=m$', fill=False, plot_points=500)
g = g1 + g2 + g3
g.axes_labels(g1.axes_labels())
g.set_axes_range(-7, 7, -6, 6)
g.save('gik_shadow_comp_th90.pdf', title=title)
show(g, title=title)
rmin : 2.34729637880682 rmax : 3.53208885091707 rmin_col : 2.3472963553338606 rmax_col : 3.5320888862379562 rmin : 1.00000001000000 rmax : 3.99999996000000 rmin_col : 1 rmax_col : 4.0
$\theta_{\mathcal{O}} = \frac{\pi}{4}$
th_obs = pi/4
title = r'$\theta_{\mathscr{O}}={\pi}/{4}$'
g1 = shadow_plot(0, th_obs, number_colors=1, color='purple',
legend=None, fill=False)
g2 = shadow_plot(0.5, th_obs, number_colors=1, color='peru',
legend=None, fill=False, plot_points=500)
g3 = shadow_plot(1, th_obs, number_colors=1, color='red',
legend=None, fill=False, plot_points=500)
g = g1 + g2 + g3
g.axes_labels(g1.axes_labels())
g.set_axes_range(-7, 7, -6, 6)
#g.save('gik_shadow_comp_th45.pdf', title=title)
show(g, title=title)
rmin : 2.48552158277147 rmax : 3.33625142631788 rmin_col : 2.3472963553338606 rmax_col : 3.5320888862379562 rmin : 1.05826009412625 rmax : 3.55486581066046 rmin_col : 1 rmax_col : 4.0
$\theta_{\mathcal{O}} = \frac{\pi}{6}$
th_obs = pi/6
title = r'$\theta_{\mathscr{O}}={\pi}/{6}$'
g1 = shadow_plot(0, th_obs, number_colors=1, color='purple',
legend=None, fill=False)
g2 = shadow_plot(0.5, th_obs, number_colors=1, color='peru',
legend=None, fill=False, plot_points=500)
g3 = shadow_plot(1, th_obs, number_colors=1, color='red',
legend=None, fill=False, plot_points=500)
g = g1 + g2 + g3
g.axes_labels(g1.axes_labels())
g.set_axes_range(-7, 7, -6, 6)
g.save('gik_shadow_comp_th30.pdf', title=title)
show(g, title=title)
rmin : 2.59373553480987 rmax : 3.20003297615050 rmin_col : 2.3472963553338606 rmax_col : 3.5320888862379562 rmin : 1.50000001500015 rmax : 3.23205077524837 rmin_col : 1 rmax_col : 4.0
$\theta_{\mathcal{O}} = 0$
th_obs = 0.001
title = r'$\theta_{\mathscr{O}}=0$'
g1 = shadow_plot(0, th_obs, number_colors=1, color='purple',
legend=None, fill=False)
g2 = shadow_plot(0.5, th_obs, number_colors=1, color='peru',
legend=None, fill=False, plot_points=500)
g3 = shadow_plot(1, th_obs, number_colors=1, color='red',
legend=None, fill=False, plot_points=500)
g = g1 + g2 + g3
g.axes_labels(g1.axes_labels())
g.set_axes_range(-7, 7, -6, 6)
g.save('gik_shadow_comp_th00.pdf', title=title)
show(g, title=title)
rmin : 2.88260669342186 rmax : 2.88382889826327 rmin_col : 2.3472963553338606 rmax_col : 3.5320888862379562 rmin : 2.41250630313641 rmax : 2.41592046802216 rmin_col : 1 rmax_col : 4.0
a0 = 1
g1 = shadow_plot(a0, pi/2, number_colors=1, color='red', plot_points=500,
legend=r'$\theta_{\mathscr{O}}=\pi/2$', fill=False)
g2 = shadow_plot(a0, pi/6, number_colors=1, color='orange', plot_points=500,
legend=r'$\theta_{\mathscr{O}}=\pi/6$', fill=False)
g3 = shadow_plot(a0, 1e-4, number_colors=1, color='green', plot_points=500,
legend=r'$\theta_{\mathscr{O}}=0$', fill=False)
g = g1 + g2 + g3
g.axes_labels(g1.axes_labels())
g.set_axes_range(-7, 7, -6, 6)
g.set_legend_options(loc='upper left')
g
rmin : 1.00000001000000 rmax : 3.99999996000000 rmin_col : 1 rmax_col : 4.0 rmin : 1.50000001500015 rmax : 3.23205077524837 rmin_col : 1 rmax_col : 4.0 rmin : 2.41404287406783 rmax : 2.41438424713942 rmin_col : 1 rmax_col : 4.0