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.
version()
%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)
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
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
forget(r<0) assumptions()
g = M.lorentzian_metric('g')
a, b = var('a b') Sigma = r^2 + a^2*cos(th)^2
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())
g.display_comp(XBL1.frame())
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())
g.display_comp(XBL2.frame())
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
gm1 = g.inverse() gm1.display(XBL1.frame())
gm1.display(XBL2.frame())
#g.christoffel_symbols_display(XBL1)
#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())
#g.christoffel_symbols_display(XBL2)
#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())
gtt component
g.restrict(M1)[0,0]
g.restrict(M2)[0,0]
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
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 gtt=−1/N2:
M.top_charts()
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()
graph = plot_profile(NN, {a: 0.9, b: 1}, -8, 8, r'$N$') graph.save('lapse.pdf') graph
gtϕ component
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
graph = plot_profile(g, {a: 2, b: 0.5}, -4, 4, r'$g_{\phi\phi}$', comp=(3,3)) graph.save('g_phiphi.pdf') graph
g.restrict(M1)[3,3].expr().factor()
g.restrict(M2)[3,3].expr().factor()
graph = plot_profile(g, {a: 0.9, b: 1}, -4, 4, r'$g_{rr}$', comp=(1,1)) graph.save('g_rr.pdf') graph
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
Ric = g.ricci() ; print(Ric)
Ric.display_comp(XBL1.frame())
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() Rscal.display()
graph = plot_profile(Rscal, {a: 0.9, b: 1}, -3, 3, r'$R$') graph.save('ricci_scalar_HT.pdf') graph
Riemann tensor
R = g.riemann() ; print(R)
R[0,1,2,3]
Kretschmann scalar
The tensor R♭, of components Rabcd=gamR bcdm:
dR = R.down(g); print(dR)
The tensor R♯, of components Rabcd=gbpgcqgdrR pqra:
uR = R.up(g); print(uR)
The Kretschmann scalar K:=RabcdRabcd:
Kr_scalar = uR['^{abcd}']*dR['_{abcd}'] Kr_scalar.display()
K = Kr_scalar.expr() K
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→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 g3=2b2).
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)
graph = plot_profile(Kr_scalar, {a: 0.9, b: 1}, -3, 3, r'$K$') graph.save('Kretschmann_scalar_HT.pdf') graph
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()
Chern-Pontryagin scalar
We start by getting the Levi-Civita 4-vector ϵabcd:
epsv = g.volume_form(contra=4) print(epsv)
epsv.display()
The dual Riemann tensor is computed as ∗Rabcd=21ϵabmnRmnrsgrcgsd with 21ϵabmnRmnrs as a first step:
tmp = 1/2*epsv['^{abmn}']*dR['_{mnrs}'] print(tmp)
dual_Riem = tmp.up(g) print(dual_Riem)
dual_Riem[0,1,2,3]
The Chern-Pontryagin scalar is ∗RabcdRabcd:
CP_scalar = dual_Riem['^{abcd}']*dR['_{abcd}'] CP_scalar.display()
CP = CP_scalar.expr().factor() CP
The Kerr value is obtained for b=0:
CP.subs(b=0).factor()
num = _.numerator()/(96*a*cos(th)*r) num
num.expand()
The above value of the Chern-Pontryagin scalar coincides with −K2 as given by Eq. (32) with Q=0 of Cherubini et al., IJMPD 11, 827 (2002).
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)
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 RabRab:
Ric2_scalar = Ric['_ab'] * Ric.up(g)['^ab'] Ric2_scalar.display()
Ric2 = Ric2_scalar.expr().factor() Ric2
The Kerr value is obtained for b=0; we check that it is zero:
Ric2.subs({b: 0})
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)
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
E_scalar=-Kr_scalar+4*Ric2_scalar-Rscal*Rscal
graph = plot_profile(E_scalar, {a: 0.9, b: 1}, -3, 3, r'$E$') graph.save('Euler_scalar.pdf') graph