Kerr_shadow.ipynb

Kernel: SageMath 9.3.beta7

Shadow and critical curve of a Kerr black hole

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

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$.

a, r = var('a r') 

lsph(a, r) = (r^2*(3 - r) - a^2*(r + 1))/(a*(r -1))
lsph

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left( a, r \right) \ {\mapsto} \ -\frac{a^{2} {\left(r + 1\right)} + {\left(r - 3\right)} r^{2}}{a {\left(r - 1\right)}}$$

qsph(a, r) = r^3 / (a^2*(r - 1)^2) * (4*a^2 - r*(r - 3)^2)
qsph

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left( a, r \right) \ {\mapsto} \ -\frac{{\left({\left(r - 3\right)}^{2} r - 4 \, a^{2}\right)} r^{3}}{a^{2} {\left(r - 1\right)}^{2}}$$

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

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

rph_ss(a) = 1/2 + cos(2/3*asin(a) + 2*pi/3)
rph_ss

$$\renewcommand{\Bold}[1]{\mathbf{#1}}a \ {\mapsto}\ \cos\left(\frac{2}{3} \, \pi + \frac{2}{3} \, \arcsin\left(a\right)\right) + \frac{1}{2}$$

rph_s(a) = 4*cos(acos(-a)/3 + 4*pi/3)^2
rph_s

$$\renewcommand{\Bold}[1]{\mathbf{#1}}a \ {\mapsto}\ 4 \, \cos\left(\frac{4}{3} \, \pi + \frac{1}{3} \, \arccos\left(-a\right)\right)^{2}$$

rph_p(a) = 4*cos(acos(-a)/3)^2
rph_p

$$\renewcommand{\Bold}[1]{\mathbf{#1}}a \ {\mapsto}\ 4 \, \cos\left(\frac{1}{3} \, \arccos\left(-a\right)\right)^{2}$$

rph_m(a) = 4*cos(acos(a)/3)^2
rph_m

$$\renewcommand{\Bold}[1]{\mathbf{#1}}a \ {\mapsto}\ 4 \, \cos\left(\frac{1}{3} \, \arccos\left(a\right)\right)^{2}$$

We add the radius of the marginally stable orbit:

rph_ms(a) = 1 - (1 - a^2)^(1/3)
rph_ms

$$\renewcommand{\Bold}[1]{\mathbf{#1}}a \ {\mapsto}\ -{\left(-a^{2} + 1\right)}^{\frac{1}{3}} + 1$$

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

$$\renewcommand{\Bold}[1]{\mathbf{#1}}a \ {\mapsto}\ 2 \, \sqrt{-\frac{1}{3} \, a^{2} + 1} \cos\left(\frac{1}{3} \, \arccos\left(-\frac{a^{2} - 1}{{\left(-\frac{1}{3} \, a^{2} + 1\right)}^{\frac{3}{2}}}\right)\right) + 1$$

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

$$\renewcommand{\Bold}[1]{\mathbf{#1}}a \ {\mapsto}\ 2 \, \sqrt{-\frac{1}{3} \, a^{2} + 1} \cos\left(\frac{2}{3} \, \pi + \frac{1}{3} \, \arccos\left(-\frac{a^{2}}{{\left(-\frac{1}{3} \, a^{2} + 1\right)}^{\frac{3}{2}}} + \frac{1}{{\left(-\frac{1}{3} \, a^{2} + 1\right)}^{\frac{3}{2}}}\right)\right) + 1$$

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

$$\renewcommand{\Bold}[1]{\mathbf{#1}}r_{\rm ph}^{**} = -0.477673658836338$$

$$\renewcommand{\Bold}[1]{\mathbf{#1}}r_{\rm ph}^{\rm ms} = 0.539741795874205$$

$$\renewcommand{\Bold}[1]{\mathbf{#1}}r_{\rm ph}^{*} = 0.658372153864346$$

$$\renewcommand{\Bold}[1]{\mathbf{#1}}r_- = 0.687750100080080$$

$$\renewcommand{\Bold}[1]{\mathbf{#1}}r_+ = 1.31224989991992$$

$$\renewcommand{\Bold}[1]{\mathbf{#1}}r_{\rm ph}^+ = 1.38628052846298$$

$$\renewcommand{\Bold}[1]{\mathbf{#1}}r_{\rm ph}^- = 3.95534731767268$$

$$\renewcommand{\Bold}[1]{\mathbf{#1}}r_{\rm ph}^{\rm pol} = 2.49269429554008$$

$$\renewcommand{\Bold}[1]{\mathbf{#1}}r_{\rm ph}^{\rm pol,in} = -0.399338575773941$$

r1 = rph_p(a0)
r2 = rph_m(a0)

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

r1, r2

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(1.38628052846298, 3.95534731767268\right)$$

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

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\theta_{\mathscr{O}} = 1.57079632679490$$

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

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(0.000000000000000, 0.000000000000000\right)$$

fa(rmin), fa(rmax)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(-2.58221244493685, 6.91641650063537\right)$$

fb(rmin), fb(rmax)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(5.73939650993649 \times 10^{-24} + 9.37314581466670 \times 10^{-8}i, 5.90445633372987 \times 10^{-8}\right)$$

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

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

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)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}4.87509181305858$$

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

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(5.19615242270663, 4.82842712474619\right)$$

Relative amplitude:

1 - R_a_1 / R_a_0 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}0.0707687665884322$$

Shadow for $a=1$

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

Screen coordinates for $a=m$

th = var('th', latex_name=r'\theta_{\mathscr{O}}')
r0 = var('r0', latex_name=r'r_0')

alpha(1, th, r0)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\frac{{r_0}^{2} - 2 \, {r_0} - 1}{\sin\left({\theta_{\mathscr{O}}}\right)}$$

beta(1, th, r0).simplify_full()

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\sqrt{\frac{{r_0}^{4} + \cos\left({\theta_{\mathscr{O}}}\right)^{4} - 4 \, {r_0}^{3} + 2 \, {\left({r_0}^{2} + 2 \, {r_0}\right)} \cos\left({\theta_{\mathscr{O}}}\right)^{2}}{\cos\left({\theta_{\mathscr{O}}}\right)^{2} - 1}}$$

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

alpha(1, pi/2, r0)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}{r_0}^{2} - 2 \, {r_0} - 1$$

beta(1, pi/2, r0).simplify_full()

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\sqrt{-{r_0}^{4} + 4 \, {r_0}^{3}}$$

Special values for $r_0 = m$:

alpha(1, th, 1)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}-\frac{2}{\sin\left({\theta_{\mathscr{O}}}\right)}$$

beta(1, th, 1).simplify_full()

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\sqrt{\frac{\cos\left({\theta_{\mathscr{O}}}\right)^{4} + 6 \, \cos\left({\theta_{\mathscr{O}}}\right)^{2} - 3}{\cos\left({\theta_{\mathscr{O}}}\right)^{2} - 1}}$$

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)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(\arcsin\left(\sqrt{3} - 1\right), 0.821327461377416, 47.0585971351201\right)$$

(cos(th_crit)^4 + 6*cos(th_crit)^2 - 3).simplify_full()

$$\renewcommand{\Bold}[1]{\mathbf{#1}}0$$

beta(1, th_crit, 1).simplify_full()

$$\renewcommand{\Bold}[1]{\mathbf{#1}}0$$

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

$$\renewcommand{\Bold}[1]{\mathbf{#1}}1.73205080756888$$

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

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

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

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

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left( a, {\ell}, q, r \right) \ {\mapsto} \ r^{4} - a^{2} q + {\left(a^{2} - {\ell}^{2} - q\right)} r^{2} + 2 \, {\left({\left(a - {\ell}\right)}^{2} + q\right)} r$$

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

R1(q, r) = R(1, 2, q, r)
R1

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left( q, r \right) \ {\mapsto} \ r^{4} - {\left(q + 3\right)} r^{2} + 2 \, {\left(q + 1\right)} r - q$$

R1(q, r).factor()

$$\renewcommand{\Bold}[1]{\mathbf{#1}}{\left(r^{2} - q + 2 \, r\right)} {\left(r - 1\right)}^{2}$$

R1(3, r).factor()

$$\renewcommand{\Bold}[1]{\mathbf{#1}}{\left(r + 3\right)} {\left(r - 1\right)}^{3}$$

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

Comparing the critical curves at the same inclination angle

$\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

Comparing the critical curves at the same value of $a$

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