Shadow and critical curve of a Kerr black hole

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

In [1]:
version()
Out[1]:
'SageMath version 9.3.beta7, Release Date: 2021-02-07'
In [2]:
%display latex

Functions $\ell_{\rm c}(r_0)$ and $q_{\rm c}(r_0)$ for critical null geodesics

We use $m=1$ and denote $r_0$ simply by $r$.

In [3]:
a, r = var('a r') 
In [4]:
lsph(a, r) = (r^2*(3 - r) - a^2*(r + 1))/(a*(r -1))
lsph
Out[4]:
In [5]:
qsph(a, r) = r^3 / (a^2*(r - 1)^2) * (4*a^2 - r*(r - 3)^2)
qsph
Out[5]:

The radii $r_+$ and $r_-$ of the two horizons:

In [6]:
rp(a) = 1 + sqrt(1 - a^2)
rm(a) = 1 - sqrt(1 - a^2)

Critical radii $r_{\rm ph}^{**}$, $r_{\rm ph}^*$, $r_{\rm ph}^+$, $r_{\rm ph}^-$, $r_{\rm ph}^{\rm ms}$ and $r_{\rm ph}^{\rm pol}$

In [7]:
rph_ss(a) = 1/2 + cos(2/3*asin(a) + 2*pi/3)
rph_ss
Out[7]:
In [8]:
rph_s(a) = 4*cos(acos(-a)/3 + 4*pi/3)^2
rph_s
Out[8]:
In [9]:
rph_p(a) = 4*cos(acos(-a)/3)^2
rph_p
Out[9]:
In [10]:
rph_m(a) = 4*cos(acos(a)/3)^2
rph_m
Out[10]:

We add the radius of the marginally stable orbit:

In [11]:
rph_ms(a) = 1 - (1 - a^2)^(1/3)
rph_ms
Out[11]:

as well as the radius of outer and inner polar orbits:

In [12]:
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
Out[12]:
In [13]:
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
Out[13]:
In [14]:
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)))
In [15]:
r1 = rph_p(a0)
r2 = rph_m(a0)

#r1 = n(rph_ss(a0))
#r2 = 0

r1, r2
Out[15]:
In [16]:
r = var('r')
plot(qsph(a0, r), (r, r1, r2)) + plot(r^3*(4 - r), (r, r1, r2), color='red')
Out[16]:
In [17]:
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,))
In [18]:
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)
In [19]:
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))
In [20]:
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)
In [21]:
plot(Theta(a0, th_obs, r), (r, r1, r2))
Out[21]:
In [22]:
rmin, rmax = r0_bounds(a0, th_obs)
rmin - r1, rmax - r2
Out[22]:
In [23]:
fa(rmin), fa(rmax)
Out[23]:
In [24]:
fb(rmin), fb(rmax)
Out[24]:
In [25]:
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
In [26]:
shadow_plot(0, pi/2)
Out[26]:
In [27]:
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
Out[27]:
In [28]:
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
Out[28]:
In [29]:
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
Out[29]:

Case $\theta_{\mathscr{O}} = 0$

In [30]:
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
Out[30]:

Shadow radius for the on-axis observer:

In [31]:
def R_shad(a):
    r0 = rph_pol(a)
    return 2*sqrt(r0*(r0^2 - a^2))/(r0 - 1)
In [32]:
R_shad(a0)
Out[32]:
In [33]:
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
Out[33]:

Extremal values:

In [34]:
R_a_0 = n(3*sqrt(3))
R_a_1 = n(2*(sqrt(2)+1))
R_a_0, R_a_1
Out[34]:

Relative amplitude:

In [35]:
1 - R_a_1 / R_a_0 
Out[35]:

Shadow for $a=1$

The partial shadow boundary:

In [36]:
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
Out[36]:

With the NHEK line:

In [37]:
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
Out[37]:

As soon as $a\neq m$, the shadow boundary is a smooth closed curve:

In [38]:
shadow_plot(0.9999, pi/2, fill=False)
rmin :  1.01637428079421   rmax :  3.99991107028894
rmin_col :  1.01637427063047   rmax_col :  3.9999111102880467
Out[38]:

Screen coordinates for $a=m$

In [39]:
th = var('th', latex_name=r'\theta_{\mathscr{O}}')
r0 = var('r0', latex_name=r'r_0')
In [40]:
alpha(1, th, r0)
Out[40]:
In [41]:
beta(1, th, r0).simplify_full()
Out[41]:

Special values for $\theta_{\mathscr{O}} = \pi/2$:

In [42]:
alpha(1, pi/2, r0)
Out[42]:
In [43]:
beta(1, pi/2, r0).simplify_full()
Out[43]:

Special values for $r_0 = m$:

In [44]:
alpha(1, th, 1)
Out[44]:
In [45]:
beta(1, th, 1).simplify_full()
Out[45]:

Critical inclination for the existence of the NHEK line:

In [46]:
th_crit = asin(sqrt(3) - 1)
th_crit, n(th_crit), n(th_crit/pi*180)
Out[46]:
In [47]:
(cos(th_crit)^4 + 6*cos(th_crit)^2 - 3).simplify_full()
Out[47]:
In [48]:
beta(1, th_crit, 1).simplify_full()
Out[48]:
In [49]:
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
Out[49]:
In [50]:
shadow_plot(1, 1.1*th_crit)
rmin :  1.00000001000000   rmax :  3.67516782408563
rmin_col :  1   rmax_col :  4.0
Out[50]:
In [51]:
shadow_plot(1, 0.9*th_crit)
rmin :  1.13415225868399   rmax :  3.50327928110302
rmin_col :  1   rmax_col :  4.0
Out[51]:
In [52]:
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
Out[52]:
In [53]:
n(sqrt(3))
Out[53]:

Cardiod shape for $\theta_{\mathscr{O}} = \pi/2$

In [54]:
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
Out[54]:
In [55]:
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
Out[55]:

Radial polynomial $\mathcal{R}(r)$ for $a=m$ and $\ell = 2m$

Generic form of $\mathcal{R}(r)$:

In [56]:
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
Out[56]:

Specialization to $a=m$ and $\ell = 2m$:

In [57]:
R1(q, r) = R(1, 2, q, r)
R1
Out[57]:
In [58]:
R1(q, r).factor()
Out[58]:
In [59]:
R1(3, r).factor()
Out[59]:
In [60]:
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
Out[60]:

Comparing the critical curves at the same inclination angle

$\theta_{\mathcal{O}} = \frac{\pi}{2}$

In [61]:
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}$

In [62]:
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}$

In [63]:
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$

In [64]:
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