Rotating_Hayward_metric_curvature.ipynb

Kernel: SageMath (stable)

Rotating Hayward metric

This Jupyter/SageMath notebook is related to the article Lamy et al, arXiv:1802.01635.

The metric is that obtained by Bambi & Modesto, Phys. Lett. B 721, 329 (2013) by applying the Newman-Janis transformation to the (non-rotating) Hayward metric for regular black holes (Hayward, PRL 96, 031103 (2006))

version()

'SageMath version 8.1, Release Date: 2017-12-07'

%display latex

To speed up the computation of the Riemann tensor, we ask for parallel computations on 8 threads:

Parallelism().set(nproc=1)  # use nproc=1 on CoCalc

M = Manifold(4, 'M')
print(M)

4-dimensional differentiable manifold M

XBL.<t,r,th,ph> = M.chart(r't r th:(0,pi):\theta ph:(0,2*pi):\phi')
XBL

g = M.lorentzian_metric('g')
b, a = var('b a')
m = r^3/(r^3+2*b^2)
Sigma = r^2+a^2*cos(th)^2
Delta = r^2-2*m*r+a^2
#m=function('m')(r)
#Sigma=function('Sigma')(r,th)
#Delta=function('Delta')(r,th)
g[0,0] = -(1 - 2*r*m/Sigma)
g[0,3] = -2*a*r*sin(th)^2*m/Sigma
g[1,1] = Sigma/Delta
g[2,2] = Sigma
g[3,3] = (r^2 + a^2 + 2*r*(a*sin(th))^2*m/Sigma)*sin(th)^2
g.display()

g.display_comp()

gm1 = g.inverse()
gm1.display()

g.christoffel_symbols_display()

$g_{tt}$ component

g[0,0]

graph = plot(g[0,0].expr().subs(a=0.9, b=1, th=0), (r, -2, 8), color='lightblue', legend_label=r'$\theta=0$')
graph += plot(g[0,0].expr().subs(a=0.9, b=1, th=pi/4), (r, -2, 8), color='magenta', legend_label=r'$\theta=\pi/4$')
graph += plot(g[0,0].expr().subs(a=0.9, b=1, th=pi/2), (r, -2, 8), color='blue', legend_label=r'$\theta=\pi/2$')
show(graph, axes_labels=[r'$r/m$', r'$g_{tt}$'], ymin=-10, ymax=10)

graph.save('g_tt.pdf', axes_labels=[r'$r/m$', r'$g_{tt}$'], ymin=-10, ymax=10)

p1 = plot3d(g[0,0].expr().subs(a=0.9, b=1), (r,0,1), (th,0,2*pi), color='red')
p2 = plot3d(0, (r,0,1), (th,0,2*pi), color='blue')
show(p1+p2, viewer='threejs', online=True, axes_labels=['r','theta','g_tt'])

MIME type unknown not supported

Lapse function

The lapse function $N$ is deduced from the standard formula $g^{tt} = - 1/N^2$:

NN = sqrt(- 1 / g.inverse()[0,0])
NN

graph = plot(NN.expr().subs(a=0.9, b=1, th=0), (r, -2, 8), color='lightblue', legend_label=r'$\theta=0$')
graph += plot(NN.expr().subs(a=0.9, b=1, th=pi/4), (r, -2, 8), color='magenta', legend_label=r'$\theta=\pi/4$')
graph += plot(NN.expr().subs(a=0.9, b=1, th=pi/2), (r, -2, 8), color='blue', legend_label=r'$\theta=\pi/2$')
show(graph, axes_labels=[r'$r/m$', r'$N$'])

$g_{t\phi}$ component

graph = plot(g[0,3].expr().subs(a=0.9, b=1, th=0), (r, -2, 8), color='lightblue', legend_label=r'$\theta=0$')
graph += plot(g[0,3].expr().subs(a=0.9, b=1, th=pi/4), (r, -2, 8), color='magenta', legend_label=r'$\theta=\pi/4$')
graph += plot(g[0,3].expr().subs(a=0.9, b=1, th=pi/2), (r, -2, 8), color='blue', legend_label=r'$\theta=\pi/2$')
show(graph, axes_labels=[r'$r/m$', r'$g_{t\phi}$'], ymin=-4, ymax=4)

p1 = plot3d(g[0,3].expr().subs(a=0.9, b=1), (r,0,1), (th,0,2*pi), color='red')
p2 = plot3d(0, (r,0,1), (th,0,2*pi), color='blue')
show(p1+p2, viewer='threejs', online=True, axes_labels=['r','theta','g_03'])

MIME type unknown not supported

Other metric components

graph = plot(g[3,3].expr().subs(a=0.9, b=1, th=0), (r, -2, 4), color='lightblue', legend_label=r'$\theta=0$')
graph += plot(g[3,3].expr().subs(a=0.9, b=1, th=pi/4), (r, -2, 4), color='magenta', legend_label=r'$\theta=\pi/4$')
graph += plot(g[3,3].expr().subs(a=0.9, b=1, th=pi/2), (r, -2, 4), color='blue', legend_label=r'$\theta=\pi/2$')
show(graph, axes_labels=[r'$r/m$', r'$g_{\phi\phi}$'], ymin=-20, ymax=20)

g[1,1].expr().factor()

graph = plot(g[1,1].expr().subs(a=0.9, b=1, th=0), (r, -2, 4), color='lightblue', legend_label=r'$\theta=0$')
graph += plot(g[1,1].expr().subs(a=0.9, b=1, th=pi/4), (r, -2, 4), color='magenta', legend_label=r'$\theta=\pi/4$')
graph += plot(g[1,1].expr().subs(a=0.9, b=1, th=pi/2), (r, -2, 4), color='blue', legend_label=r'$\theta=\pi/2$')
show(graph, axes_labels=[r'$r/m$', r'$g_{rr}$'], ymin=-10, ymax=10)

graph.save("g_rr.pdf", axes_labels=[r'$r/m$', r'$g_{rr}$'], ymin=-10, ymax=10)

graph = plot(g[2,2].expr().subs(a=0.9, b=1, th=0), (r, -2, 4), color='lightblue', legend_label=r'$\theta=0$')
graph += plot(g[2,2].expr().subs(a=0.9, b=1, th=pi/4), (r, -2, 4), color='magenta', legend_label=r'$\theta=\pi/4$')
graph += plot(g[2,2].expr().subs(a=0.9, b=1, th=pi/2), (r, -2, 4), color='blue', legend_label=r'$\theta=\pi/2$')
show(graph, axes_labels=[r'$r/m$', r'$g_{\theta\theta}$'])

Ricci tensor

Ric = g.ricci() ; print(Ric)

Field of symmetric bilinear forms Ric(g) on the 4-dimensional differentiable manifold M

Ric.display_comp()

We check that for $b=0$, we are dealing with a solution of the vacuum Einstein equation:

all([all([Ric[i,j].expr().subs(b=0) == 0 for i in range(4)]) for j in range(4)])

The Ricci scalar:

Rscal = g.ricci_scalar().expr()
Rscal

latex(Rscal.factor())

graph = plot(Rscal.subs(a=0.9, b=1, th=0), (r, -2, 2), color='lightblue', legend_label=r'$\theta=0$')
graph += plot(Rscal.subs(a=0.9, b=1, th=pi/4), (r, -2, 2), color='magenta', legend_label=r'$\theta=\pi/4$')
graph += plot(Rscal.subs(a=0.9, b=1, th=pi/2), (r, -2, 2), color='blue', legend_label=r'$\theta=\pi/2$')
show(graph, axes_labels=[r'$r/m$', r'$R$'], ymin=-100, ymax=100)

graph.save("Ricci_scalar_Hayward.pdf", axes_labels=[r'$r/m$', r'$R$'], ymin=-100, ymax=100)

plot3d(Rscal.subs({a: 0.9, b: 1}), (r, -1/2, 1/2), (th, 0, pi/2), plot_points=200,
       aspect_ratio=[2, 1, 0.1], axes_labels=['r', 'theta', 'R'],
       viewer='threejs', online=True)

MIME type unknown not supported

Riemann tensor

R = g.riemann() ; print(R)

Tensor field Riem(g) of type (1,3) on the 4-dimensional differentiable manifold M

R[0,1,2,3]

Kretschmann scalar

The tensor $R^\flat$, of components $R_{abcd} = g_{am} R^m_{\ \, bcd}$:

dR = R.down(g); print(dR)

Tensor field of type (0,4) on the 4-dimensional differentiable manifold M

The tensor $R^\sharp$, of components $R^{abcd} = g^{bp} g^{cq} g^{dr} R^a_{\ \, pqr}$:

uR = R.up(g); print(uR)

Tensor field of type (4,0) on the 4-dimensional differentiable manifold M

The Kretschmann scalar $K := R^{abcd} R_{abcd}$:

Kr_scalar = uR['^{abcd}']*dR['_{abcd}']
Kr_scalar.display()

K = Kr_scalar.expr()
K

graph = plot(Kr_scalar.expr().subs(a=0.9, b=1, th=0), (r, -2, 2), color='lightblue', legend_label=r'$\theta=0$')
graph += plot(Kr_scalar.expr().subs(a=0.9, b=1, th=pi/4), (r, -2, 2), color='magenta', legend_label=r'$\theta=\pi/4$')
graph += plot(Kr_scalar.expr().subs(a=0.9, b=1, th=pi/2), (r, -2, 2), color='blue', legend_label=r'$\theta=\pi/2$')
show(graph, axes_labels=[r'$r/m$', r'$K$'], ymin=-100, ymax=100)

graph.save("Kretschmann_scalar_H.pdf", axes_labels=[r'$r/m$', r'$K$'], ymin=-100, ymax=100)

K_0 = K.subs(r=0)
K_0

The equatorial value of the Kretschmann scalar is

K_eq = K.subs(th=pi/2).simplify_full()
K_eq

taylor(K_eq, r, 0, 6)

The limit $r\rightarrow 0$:

K_eq.subs(r=0)

We recover the same value as that given by Eq. (24) of Bambi & Modesto, Phys. Lett. B 721, 329 (2013) (note that the quantity $g$ used by Bambi & Modesto is related to our $b$ by $g^3 = 2 b^2$).

K1 = K.subs({a: 0.9, b: 1})
plot3d(K1, (r, -1/2, 1/2), (th, 0, pi/2), plot_points=200,
       aspect_ratio=[2, 1, 0.1], axes_labels=['r', 'theta', 'Kr'],
       viewer='threejs', online=True)

MIME type unknown not supported

Non-rotating limit

Ka0 = K.subs(a=0).simplify_full()
Ka0

taylor(Ka0, r, 0, 6)

Check: we recover Schwarzschild value when $b=0$:

Ka0.subs(b=0).simplify_full()

latex(K.simplify_full())