Kerr_extremal_throat_emb.ipynb

Kernel: SageMath 9.4.beta3

Isometric embedding of the extremal Kerr throat

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

version()

'SageMath version 9.4.beta3, Release Date: 2021-06-21'

%display latex

r = var('r')
mu = var('mu', latex_name=r'\mu')
rho(r, mu) = sqrt(1 - mu^2)*sqrt(r^2 + 1 + 2*r*(1 - mu^2)/(r^2 + mu^2))
rho

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

diff(rho(r, mu), r).simplify_full()

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\frac{\sqrt{-{\mu}^{2} + 1} {\mu}^{4} r + 2 \, \sqrt{-{\mu}^{2} + 1} {\mu}^{2} r^{3} + \sqrt{-{\mu}^{2} + 1} r^{5} + {\left({\mu}^{2} - 1\right)} \sqrt{-{\mu}^{2} + 1} r^{2} - {\left({\mu}^{4} - {\mu}^{2}\right)} \sqrt{-{\mu}^{2} + 1}}{{\left({\mu}^{4} + 2 \, {\mu}^{2} r^{2} + r^{4}\right)} \sqrt{\frac{r^{4} + {\left({\mu}^{2} + 1\right)} r^{2} + {\mu}^{2} - 2 \, {\left({\mu}^{2} - 1\right)} r}{{\mu}^{2} + r^{2}}}}$$

zp(r, mu) = sqrt( (r^2 + mu^2)/(r - 1)^2 - (diff(rho(r), r))^2 ).simplify_full()
zp(r, mu)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\sqrt{\frac{{\mu}^{2} r^{12} - 2 \, {\left({\mu}^{2} - 1\right)} r^{11} + 2 \, {\left(2 \, {\mu}^{4} + {\mu}^{2}\right)} r^{10} + 2 \, {\mu}^{10} - 2 \, {\left(3 \, {\mu}^{4} - {\mu}^{2} - 2\right)} r^{9} + {\left(6 \, {\mu}^{6} + 4 \, {\mu}^{4} + 9 \, {\mu}^{2} - 4\right)} r^{8} - 3 \, {\mu}^{8} - 2 \, {\left(5 \, {\mu}^{6} - {\mu}^{4} - 3 \, {\mu}^{2} - 1\right)} r^{7} + {\left(4 \, {\mu}^{8} + 9 \, {\mu}^{6} + 9 \, {\mu}^{4} - {\mu}^{2} - 1\right)} r^{6} + 3 \, {\mu}^{6} - 2 \, {\left(5 \, {\mu}^{8} - 6 \, {\mu}^{4} + 2 \, {\mu}^{2} - 1\right)} r^{5} + {\left({\mu}^{10} + 10 \, {\mu}^{8} + 5 \, {\mu}^{6} - 5 \, {\mu}^{4} + 5 \, {\mu}^{2} - 1\right)} r^{4} - {\mu}^{4} - 2 \, {\left(2 \, {\mu}^{10} + {\mu}^{6} - 5 \, {\mu}^{4} + 2 \, {\mu}^{2}\right)} r^{3} + {\left(7 \, {\mu}^{10} - 9 \, {\mu}^{8} + 13 \, {\mu}^{6} - 7 \, {\mu}^{4} + 2 \, {\mu}^{2}\right)} r^{2} - 2 \, {\left(3 \, {\mu}^{10} - 6 \, {\mu}^{8} + 4 \, {\mu}^{6} - {\mu}^{4}\right)} r}{r^{12} - 10 \, {\mu}^{2} r^{9} + 2 \, {\left(2 \, {\mu}^{2} + 1\right)} r^{10} - 2 \, r^{11} + 6 \, {\mu}^{8} r^{2} + 3 \, {\left(2 \, {\mu}^{4} + 4 \, {\mu}^{2} - 1\right)} r^{8} + {\mu}^{8} - 2 \, {\left(9 \, {\mu}^{4} + 2 \, {\mu}^{2} - 1\right)} r^{7} + 4 \, {\left({\mu}^{6} + 6 \, {\mu}^{4} - 2 \, {\mu}^{2}\right)} r^{6} - 2 \, {\left(7 \, {\mu}^{6} + 6 \, {\mu}^{4} - 3 \, {\mu}^{2}\right)} r^{5} + {\left({\mu}^{8} + 20 \, {\mu}^{6} - 6 \, {\mu}^{4}\right)} r^{4} - 2 \, {\left(2 \, {\mu}^{8} + 6 \, {\mu}^{6} - 3 \, {\mu}^{4}\right)} r^{3} - 2 \, {\left(2 \, {\mu}^{8} - {\mu}^{6}\right)} r}}$$

zp2 = (2*r^7 + 4*r^5 - 4*r^4 + 2*r^3 - r^2 + 2*r - 1)/(r^3*(r^3 + r + 2)*(r - 1)^2)
zp2

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

bool(zp(r, 0)^2 == zp2)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\mathrm{True}$$

s = (zp(r, mu)^2).numerator()
s.collect(r)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}{\mu}^{2} r^{12} - 2 \, {\left({\mu}^{2} - 1\right)} r^{11} + 2 \, {\left(2 \, {\mu}^{4} + {\mu}^{2}\right)} r^{10} + 2 \, {\mu}^{10} - 2 \, {\left(3 \, {\mu}^{4} - {\mu}^{2} - 2\right)} r^{9} + {\left(6 \, {\mu}^{6} + 4 \, {\mu}^{4} + 9 \, {\mu}^{2} - 4\right)} r^{8} - 3 \, {\mu}^{8} - 2 \, {\left(5 \, {\mu}^{6} - {\mu}^{4} - 3 \, {\mu}^{2} - 1\right)} r^{7} + {\left(4 \, {\mu}^{8} + 9 \, {\mu}^{6} + 9 \, {\mu}^{4} - {\mu}^{2} - 1\right)} r^{6} + 3 \, {\mu}^{6} - 2 \, {\left(5 \, {\mu}^{8} - 6 \, {\mu}^{4} + 2 \, {\mu}^{2} - 1\right)} r^{5} + {\left({\mu}^{10} + 10 \, {\mu}^{8} + 5 \, {\mu}^{6} - 5 \, {\mu}^{4} + 5 \, {\mu}^{2} - 1\right)} r^{4} - {\mu}^{4} - 2 \, {\left(2 \, {\mu}^{10} + {\mu}^{6} - 5 \, {\mu}^{4} + 2 \, {\mu}^{2}\right)} r^{3} + {\left(7 \, {\mu}^{10} - 9 \, {\mu}^{8} + 13 \, {\mu}^{6} - 7 \, {\mu}^{4} + 2 \, {\mu}^{2}\right)} r^{2} - 2 \, {\left(3 \, {\mu}^{10} - 6 \, {\mu}^{8} + 4 \, {\mu}^{6} - {\mu}^{4}\right)} r$$

s = (zp(r, mu)^2).denominator()
s.factor()

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

Specific choice of $\mu := \cos\theta$

#mu0 = sqrt(3)/2  # theta = pi/6
mu0 = 0           # theta = pi/2

zpf = fast_callable(zp(r, mu0), vars=[r])

def zz(r0, r1, verbose=False):
    numint = numerical_integral(zpf, r0, r1, algorithm='qags')
    error = numint[1]
    if verbose:
        print("Error = {}".format(error))
    if error > 1e-3: 
        print("Warning: error = {}".format(error))
    return numint[0]   

zz(2, 3, verbose=True)

Error = 1.6403995233359683e-14

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

zz(1.5, 3, verbose=True)

Error = 5.314212959789082e-11

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

zz(1.001, 3, verbose=True)

Error = 2.628708897735641e-09

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

zz(1.00001, 3, verbose=True)

Error = 8.104079286898108e-07

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

r_0 = 2
rmin = 1.0001
plot(lambda r: zz(r_0, r), (r, rmin, 4), 
     axes_labels = [r'$r/m$', r'$Z(r)/m$'])

rhof = fast_callable(rho(r, mu0), vars=[r])

plot(rhof, (1, 4),
     axes_labels = [r'$r/m$', r'$P(r)/m$'])

rmin = 1.00001
r_0 = 2

def xe(r, ph):
    return rhof(r)*cos(ph)

def ye(r, ph):
    return rhof(r)*sin(ph)

def ze(r, ph):
    return zz(r_0, r)

rmax = 10
graph = parametric_plot3d([xe, ye, ze], (rmin, rmax), (0, 2*pi), 
                          color='lightsteelblue')
graph

Constant $r$ curve:

def const_r(r0, color='blue', thickness=5):
    return parametric_plot3d([lambda ph: xe(r0, ph), lambda ph: ye(r0, ph), 
                              lambda ph: ze(r0, ph)], (0, 2*pi), 
                              color=color, thickness=thickness)

Adding the ergosurface, i.e. the curve $r=m(1 + \sin\theta)$:

graph += const_r(1 + sqrt(1 - mu0^2), thickness=7)

Adding the curves $r/m = 1 + 10^{-5},\ldots, 1 + 10^{-1}$:

for i in range(1, 6):
    graph += const_r(1 + 10^(-i), color='black')

show(graph, axes_labels_style=dict(fontsize=12, fontfamily='Liberation Sans'))