Project: Regular rotating black holes

Path: rotating_Hayward_metric_ext.ipynb

Views: ¹⁹⁰

Kernel: SageMath (stable)

Extended 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)), extended to cover the region $r<0$ .

In [1]:

version()

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

In [2]:

%display latex

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

In [3]:

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

In [4]:

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

4-dimensional differentiable manifold M

In [5]:

M1 = M.open_subset('M_1')
XBL1.<t,r,th,ph> = M1.chart(r't r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\phi')
XBL1

In [6]:

M2 = M.open_subset('M_2')
forget(r>0)
XBL2.<t,r,th,ph> = M2.chart(r't r:(-oo,0) th:(0,pi):\theta ph:(0,2*pi):\phi')
M._top_charts.append(XBL2)
XBL2

In [7]:

forget(r<0)
assumptions()

In [8]:

g = M.lorentzian_metric('g')

In [9]:

a, b = var('a b')
Sigma = r^2 + a^2*cos(th)^2

In [10]:

m = r^3/(r^3 + 2*b^2)
m1 = m
Delta = r^2 - 2*m*r + a^2
g1 = g.restrict(M1)
g1[0,0] = -(1 - 2*r*m/Sigma)
g1[0,3] = -2*a*r*sin(th)^2*m/Sigma
g1[1,1] = Sigma/Delta
g1[2,2] = Sigma
g1[3,3] = (r^2 + a^2 + 2*r*(a*sin(th))^2*m/Sigma)*sin(th)^2
g.display(XBL1.frame())

In [11]:

g.display_comp(XBL1.frame())

In [12]:

m = -r^3/(-r^3 + 2*b^2)
m2 = m
Delta = r^2 - 2*m*r + a^2
g2 = g.restrict(M2)
g2[0,0] = -(1 - 2*r*m/Sigma)
g2[0,3] = -2*a*r*sin(th)^2*m/Sigma
g2[1,1] = Sigma/Delta
g2[2,2] = Sigma
g2[3,3] = (r^2 + a^2 + 2*r*(a*sin(th))^2*m/Sigma)*sin(th)^2
g.display(XBL2.frame())

In [13]:

g.display_comp(XBL2.frame())

In [14]:

graph = plot(m1.subs(b=1), (r, 0, 8), axes_labels=[r'$r/m$', r'$M(r)/m$'], gridlines=True) 
graph += plot(m2.subs(b=1), (r, -8, 0))
graph

In [15]:

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

In [16]:

gm1.display(XBL2.frame())

In [17]:

#g.christoffel_symbols_display(XBL1)

In [18]:

#gam = g.christoffel_symbols(XBL1)
#for i in range(4):
#    for j in range(4):
#        for k in range(j,4):
#            print("---------------------")
#            print("Gamma^{}_{}{} for r >=0:".format(i,j,k))
#            if gam[i,j,k] == 0:
#                print(0)
#            else:
#                show(gam[i,j,k].expr().factor())
#                print(gam[i,j,k].expr().factor())

In [19]:

#g.christoffel_symbols_display(XBL2)

In [20]:

#gam = g.christoffel_symbols(XBL2)
#for i in range(4):
#    for j in range(4):
#        for k in range(j,4):
#            print("--------------------")
#            print("Gamma^{}_{}{} for r < 0:".format(i,j,k))
#            if gam[i,j,k] == 0:
#                print(0)
#            else:
#                show(gam[i,j,k].expr().factor())
#                print(gam[i,j,k].expr().factor())

$g_{tt}$ component

In [21]:

g.restrict(M1)[0,0]

In [22]:

g.restrict(M2)[0,0]

In [23]:

def plot_profile(field, param_values, rmin, rmax, y_label, comp=None, 
                 ymin=None, ymax=None, gridlines=True):
    graph = Graphics()
    rmin0 = max(rmin, 0)
    rmax0 = min(rmax, 0)
    if rmax > 0:
        f1 = field.restrict(M1)
        if comp is not None:
            f1 = f1[comp].expr()
        else:
            f1 = f1.expr()
        f1 = f1.subs(param_values)
        graph += plot(f1.subs(th=0), (r, rmin0, rmax), color='lightblue', 
                      legend_label=r'$\theta=0$', axes_labels = [r'$r/m$', y_label],
                      ymin=ymin, ymax=ymax, gridlines=gridlines)
        graph += plot(f1.subs(th=pi/4), (r, rmin0, rmax), color='magenta', 
                      legend_label=r'$\theta=\pi/4$')
        graph += plot(f1.subs(th=pi/2), (r, rmin0, rmax), color='blue', 
                      legend_label=r'$\theta=\pi/2$')
    if rmin < 0:
        f1 = field.restrict(M2)
        if comp is not None:
            f1 = f1[comp].expr()
        else:
            f1 = f1.expr()
        f1 = f1.subs(param_values)
        graph += plot(f1.subs(th=0), (r, rmin, rmax0), color='lightblue', 
                      axes_labels = [r'$r/m$', y_label],
                      ymin=ymin, ymax=ymax)
        graph += plot(f1.subs(th=pi/4), (r, rmin, rmax0), color='magenta')
        graph += plot(f1.subs(th=pi/2), (r, rmin, rmax0), color='blue')
    return graph

In [24]:

graph = plot_profile(g, {a: 0.9, b: 1}, -8, 8, r'$g_{tt}$', comp=(0,0))
graph.save('g_tt.pdf')
graph

Lapse function

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

In [25]:

M.top_charts()

In [26]:

NN = M.scalar_field(coord_expression={XBL1: sqrt(- 1 / g.restrict(M1).inverse()[0,0]),
                                      XBL2: sqrt(- 1 / g.restrict(M2).inverse()[0,0])},
                    name='N')
NN.display()

In [27]:

graph = plot_profile(NN, {a: 0.9, b: 1}, -8, 8, r'$N$')
graph.save('lapse.pdf')
graph

$g_{t\phi}$ component

In [28]:

graph = plot_profile(g, {a: 0.9, b: 1}, -8, 8, r'$g_{t\phi}$', comp=(0,3))
graph.save('g_tphi.pdf')
graph

Other metric components

In [29]:

graph = plot_profile(g, {a: 2, b: 0.5}, -4, 4, r'$g_{\phi\phi}$', comp=(3,3))
graph.save('g_phiphi.pdf')
graph

In [30]:

g.restrict(M1)[3,3].expr().factor()

In [31]:

g.restrict(M2)[3,3].expr().factor()

In [32]:

graph = plot_profile(g, {a: 0.9, b: 1}, -4, 4, r'$g_{rr}$', comp=(1,1))
graph.save('g_rr.pdf')
graph

In [33]:

graph = plot_profile(g, {a: 0.9, b: 1}, -2, 2, r'$g_{\theta\theta}$', comp=(2,2))
graph.save('g_thth.pdf')
graph

Ricci tensor

In [34]:

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

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

In [35]:

Ric.display_comp(XBL1.frame())

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

In [36]:

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

The Ricci scalar:

In [37]:

Rscal = g.ricci_scalar()
Rscal.display()

In [38]:

graph = plot_profile(Rscal, {a: 0.9, b: 1}, -3, 3, r'$R$')
graph.save('ricci_scalar_HT.pdf')
graph

Riemann tensor

In [39]:

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

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

In [40]:

R[0,1,2,3]

Kretschmann scalar

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

In [41]:

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

In [42]:

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

In [43]:

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

In [44]:

K = Kr_scalar.expr()
K

In [45]:

K_0 = K.subs(r=0)
K_0

The equatorial value of the Kretschmann scalar is

In [46]:

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

In [47]:

taylor(K_eq, r, 0, 6)

The limit $r\rightarrow 0$ :

In [48]:

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

In [49]:

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', 'K'],
       viewer='threejs', online=True)

MIME type unknown not supported

In [50]:

graph = plot_profile(Kr_scalar, {a: 0.9, b: 1}, -3, 3, r'$K$')
graph.save('Kretschmann_scalar_HT.pdf')
graph

Non-rotating limit

In [51]:

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

In [52]:

taylor(Ka0, r, 0, 6)

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

In [53]:

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

Chern-Pontryagin scalar

We start by getting the Levi-Civita 4-vector $\epsilon^{abcd}$ :

In [54]:

epsv = g.volume_form(contra=4)
print(epsv)

4-vector field on the 4-dimensional differentiable manifold M

In [55]:

epsv.display()

The dual Riemann tensor is computed as ${}^*R^{abcd} = \frac{1}{2} \epsilon^{abmn} R_{mnrs} g^{rc} g^{sd}$ with $\frac{1}{2} \epsilon^{abmn} R_{mnrs}$ as a first step:

In [56]:

tmp = 1/2*epsv['^{abmn}']*dR['_{mnrs}']
print(tmp)

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

In [57]:

dual_Riem = tmp.up(g)
print(dual_Riem)

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

In [58]:

dual_Riem[0,1,2,3]

The Chern-Pontryagin scalar is ${}^*R^{abcd} R_{abcd}$ :

In [59]:

CP_scalar = dual_Riem['^{abcd}']*dR['_{abcd}']
CP_scalar.display()

In [60]:

CP = CP_scalar.expr().factor()
CP

The Kerr value is obtained for $b=0$ :

In [61]:

CP.subs(b=0).factor()

In [62]:

num = _.numerator()/(96*a*cos(th)*r)
num

In [63]:

num.expand()

The above value of the Chern-Pontryagin scalar coincides with $-K_2$ as given by Eq. (32) with $Q=0$ of Cherubini et al., IJMPD 11, 827 (2002).

In [64]:

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

MIME type unknown not supported

In [65]:

graph = plot_profile(CP_scalar, {a: 0.9, b: 1}, -3, 3, r'$CP$')
graph.save('ChernPontryagin_scalar_HT.pdf')
graph

Ricci squared

The Ricci squared is $R_{ab} R^{ab}$ :

In [66]:

Ric2_scalar = Ric['_ab'] * Ric.up(g)['^ab']
Ric2_scalar.display()

In [67]:

Ric2 = Ric2_scalar.expr().factor()
Ric2

The Kerr value is obtained for $b=0$ ; we check that it is zero:

In [68]:

Ric2.subs({b: 0})

In [69]:

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

MIME type unknown not supported

In [70]:

graph = plot_profile(Ric2_scalar, {a: 0.9, b: 1}, -3, 3, r'$R_{ab} R^{ab}$')
graph.save('Ric2_scalar_HT.pdf')
graph

Euler invariant

In [71]:

E_scalar=-Kr_scalar+4*Ric2_scalar-Rscal*Rscal

In [73]:

graph = plot_profile(E_scalar, {a: 0.9, b: 1}, -3, 3, r'$E$')
graph.save('Euler_scalar.pdf')
graph

In [0]: