This This Jupyter/SageMath notebook is relative to the lectures Geometry and physics of black holes. It explores the global Near-Horizon Extremal Kerr (NHEK) spacetime $(\mathcal{N}, h)$.
This notebook requires a version of SageMath at least equal to 9.0:
version()
'SageMath version 9.5.beta2, Release Date: 2021-09-26'
First we set up the notebook to display mathematical objects using LaTeX rendering, and, to speed up computations, we ask for running them in parallel on 8 threads:
%display latex
Parallelism().set(nproc=8)
N = Manifold(4, 'N', latex_name=r'\mathscr{N}', structure='Lorentzian',
metric_name='h')
print(N)
4-dimensional Lorentzian manifold N
Global coordinate chart $(\tau, y, \theta, \psi)$:
X.<ta, y, th, ps> = N.chart(r"ta:\tau y th:(0,pi):\theta ps:(0,2*pi):periodic:\psi")
print(X)
X
Chart (N, (ta, y, th, ps))
X.coord_range()
The coordinate 1-forms:
X.coframe()[:]
dta, dy, dth, dps = X.coframe()[:]
The mass parameter $m$ appears in the the NHEK metric $h$ as an overall scale parameter $m^2$. In what follows, we set $m=1$ for simplicity. We then declare $h$ as
h = N.metric()
dps_ydta = dps + y*dta
h.set((1 + cos(th)^2)*(- (1 + y^2)*(dta*dta)+ (dy*dy)/(1 + y^2) + dth*dth)
+ 4*sin(th)^2/(1 + cos(th)^2)*(dps_ydta*dps_ydta))
h.display()
The NHEK metric is a solution of the vacuum Einstein equation:
h.ricci().display()
h.inverse()[:]
Some simplifications are in order:
h.inverse().apply_map(factor)
h.inverse()[:]
h00 = h.inverse()[0,0].expr().subs({sin(th)^2: 1 - cos(th)^2})
h00
h33 = h.inverse()[3,3].expr().subs({sin(th)^2: 1 - cos(th)^2})
h33
Hence we set:
h.inverse()[0,0] = h00
h.inverse()[3,3] = h33
and we get
h.inverse()[:]
Two obvious Killing vectors are $\partial/\partial\psi$ and $\partial/\partial\tau$, since the components of $h$ do not depend on $\psi$ or $\tau$:
eta = N.vector_field(0, 0, 0, 1, name='eta', latex_name=r'\eta')
eta.display()
h.lie_derivative(eta).display()
J0 = N.vector_field(1, 0, 0, 0, name='J_0')
J0.display()
h.lie_derivative(J0).display()
A third Killing vector arises from the isometry expressed by $(T,R)\mapsto(\alpha T, R/\alpha)$ in Poincaré-type coordinates:
J1 = N.vector_field(y*sin(ta)/sqrt(1 + y^2), -cos(ta)*sqrt(1 + y^2),
0, sin(ta)/sqrt(1 + y^2), name='J_1')
J1.display()
h.lie_derivative(J1).display()
Finally a fourth Killing vector is $\partial/\partial T$ of the Poincaré-type coordinates. We actually consider the linear combination $\partial/\partial T - \partial/\partial\tau$ since it has slightly simpler components in terms of the global coordinates $(\tau,y,\theta,\psi)$:
J2 = N.vector_field(y*cos(ta)/sqrt(1 + y^2), sin(ta)*sqrt(1 + y^2),
0, cos(ta)/sqrt(1 + y^2), name='J_2')
J2.display()
h.lie_derivative(J2).display()
First of all we notice that $\eta$ commutes with all the three other Killing vectors:
[eta.bracket(j).display() for j in (J0, J1, J2)]
To write easily the other commutation relations, we supplement $(J_0, J_1, J_2)$ by $\partial/\partial\theta$ to make it a vector frame on $\mathscr{N}$:
J = N.vector_frame('J', (J0, J1, J2, X.frame()[2]))
J0, J1, J2 = J[:-1]
J
J2.display(J)
for v in J:
show(v.display())
Then we ask for the display of the commutators in that frame:
J0.bracket(J1).display(J)
J0.bracket(J2).display(J)
J1.bracket(J2).display(J)
These commutation relations are not the standard ones for $\mathfrak{sl}(2,\mathbb{R})$; in order to get these, let us introduce the following linear combinations:
K0 = 2*J1
K0.set_name('K_0')
K0.display(J)
K1 = J2 - J0
K1.set_name('K_1')
K1.display(J)
K2 = J2 + J0
K2.set_name('K_2')
K2.display(J)
K = N.vector_frame('K', (K0, K1, K2, X.frame()[2]))
K0, K1, K2 = K[:-1]
K
for v in K:
show(v.display())
We have then
K0.bracket(K1).display(K)
K0.bracket(K2).display(K)
K1.bracket(K2).display(K)
The above commutation relations are exactly those of $\mathfrak{sl}(2,\mathbb{R})$ in the standard matrix representation:
sl2 = lie_algebras.sl(QQ, 2, representation='matrix') # QQ instead of RR to deal with an exact field
E,F,H = sl2.gens()
E,F,H
all([H.bracket(E) == 2*E,
H.bracket(F) == -2*F,
E.bracket(F) == H])
The spacetime total angular momentum $J$ is computed from the axisymmetric Killing vector $\eta$ via the Komar integral: $$ J = \frac{1}{16\pi} \oint_{\mathscr{S}} \star (\mathrm{d}\underline{\eta}) $$ where
In vacuum, the Komar integral $J$ is independent of the choice of $\mathscr{S}$. Let us choose $\mathscr{S}$ to be a sphere $(\tau, y) = \mathrm{const}$. We have then $$ J = \frac{1}{16\pi} \int_0^{2\pi} \int_0^\pi \left(\star (\mathrm{d}\underline{\eta}) \right)_{\theta\psi} \mathrm{d}\theta\; \mathrm{d}\psi$$
The Killing vector $\eta$:
eta.display()
The 1-form $\underline{\eta}$:
ueta = eta.down(h)
ueta.set_name('ueta', latex_name=r'\underline{\eta}')
ueta.display()
The 2-form $\mathrm{d}\underline{\eta}$:
deta = diff(ueta)
deta.display()
The 2-form $\star (\mathrm{d}\underline{\eta})$:
sdeta = deta.hodge_dual()
sdeta.display()
The $\theta\psi$ component:
sdeta23 = sdeta[2,3].expr()
sdeta23
This can be simplified to $\frac{8\sin^3\theta}{(1 + \cos^2\theta)^2}$:
A = 8*sin(th)^3/(1 + cos(th)^2)^2
A
(sdeta23^2 - A^2).simplify_full()
Hence we have
J_Komar = 1/(16*pi)*integrate(integrate(A, (th, 0, pi)), (ps, 0, 2*pi))
J_Komar
If we restore the metric scale factor $m^2$, we have $h \to m^2 h$, $\underline{\eta} \to m^2 \underline{\eta}$, $\mathrm{d}\underline{\eta} \to m^2 \mathrm{d}\underline{\eta}$ and $\star(\mathrm{d}\underline{\eta}) \to m^2 \star(\mathrm{d}\underline{\eta})$. The last property holds because in dimension 4, the Hodge dual is conformally invariant on 2-forms (as it can be easily checked from its definition).
Hence we conclude that
$$ J = m^2 $$
In other words, the angular momentum of the NHEK spacetime is identical to the angular momentum of the Kerr spacetime from which it derives.
For the record:
integrate(A, (th, 0, pi))
integrate(A, th)
bool(diff(_, th) == A) # a check
Let us introduce the "conformal" global coordinates $(\tau,\chi,\theta,\psi)$ such that $y = \tan\chi$:
XC.<ta, ch, th, ps> = N.chart(r"ta:\tau ch:(-pi/2,pi/2):\chi th:(0,pi):\theta ps:(0,2*pi):periodic:\psi")
print(XC)
XC
Chart (N, (ta, ch, th, ps))
XC.coord_range()
X_to_XC = X.transition_map(XC, [ta, atan(y), th, ps])
X_to_XC.display()
X_to_XC.inverse().display()
h.display(XC)
for v in J[:-1]:
show(v.display(XC))
for v in K[:-1]:
show(v.display(XC))
NP = N.open_subset('NP', latex_name=r'\mathscr{N}_{\rm P}',
coord_def={XC: [-ch - pi/2 < ta, ta < - ch - 3*pi/2]})
print(NP)
Open subset NP of the 4-dimensional Lorentzian manifold N
XP.<T,R,th,Ph> = NP.chart(r"T R th:(0,pi):\theta Ph:(0,2*pi):periodic:\Phi")
print(XP)
XP
Chart (NP, (T, R, th, Ph))
Conf_to_Poinc = XC.transition_map(XP, [sin(ta)/(cos(ta) + sin(ch)),
(cos(ta) + sin(ch))/cos(ch),
th,
ps + log((cos(ta)*cos(ch) + sin(ta)*sin(ch))/(sin(ta) + cos(ch)))])
Conf_to_Poinc.display()
Conf_to_Poinc.set_inverse(atan2(2*R^2*T, R^2*(1 - T^2) +1) + pi*unit_step(-R),
atan((R^2*(1 + T^2) - 1)/(2*R)),
th,
Ph - ln(((1 - T*R)^2 + R^2)/sqrt(((1 + T^2)*R^2 - 1)^2 + 4*R^2)))
Check of the inverse coordinate transformation: ta == pi*unit_step(-(cos(ta) + sin(ch))/cos(ch)) + arctan2(2*(cos(ta) + sin(ch))*sin(ta)/cos(ch)^2, 2*(cos(ta)^2 + cos(ta)*sin(ch))/cos(ch)^2) **failed** ch == arctan(sin(ch)/cos(ch)) **failed** th == th *passed* ps == ps *passed* T == (R^3*T^2*sin(pi*unit_step(-R)) - 2*R^3*T*cos(pi*unit_step(-R)) - R^3*sin(pi*unit_step(-R)) - R*sin(pi*unit_step(-R)))/(2*R^3*T*sin(pi*unit_step(-R)) - R^3*cos(pi*unit_step(-R)) + (R^3*cos(pi*unit_step(-R)) - R^2*abs(R))*T^2 - (R^2 - 1)*abs(R) - R*cos(pi*unit_step(-R))) **failed** R == -1/2*(2*R^3*T*sin(pi*unit_step(-R)) - R^3*cos(pi*unit_step(-R)) + (R^3*cos(pi*unit_step(-R)) - R^2*abs(R))*T^2 - (R^2 - 1)*abs(R) - R*cos(pi*unit_step(-R)))/(R*abs(R)) **failed** th == th *passed* Ph == Ph + log((R^2*T^2*abs(R)*sin(pi*unit_step(-R)) - 2*(R^2*cos(pi*unit_step(-R)) - R*sin(pi*unit_step(-R)))*T*abs(R) - (R^2*sin(pi*unit_step(-R)) + 2*R*cos(pi*unit_step(-R)) - sin(pi*unit_step(-R)))*abs(R))/(R^3*T^2*sin(pi*unit_step(-R)) - 2*R^3*T*cos(pi*unit_step(-R)) - R^3*sin(pi*unit_step(-R)) - 2*R*abs(R) - R*sin(pi*unit_step(-R)))) **failed** NB: a failed report can reflect a mere lack of simplification.
Conf_to_Poinc.inverse().display()
graph_XP = XP.plot(chart=XC, ambient_coords=(ch, ta), fixed_coords={th: pi/2, Ph: 0},
ranges={T: (-9, 9), R: (-9, 8)}, color={T: 'red', R: 'grey'},
number_values=17, plot_points=200)
show(graph_XP, aspect_ratio=1)
NP1 = NP.open_subset('NP1', latex_name=r'\mathscr{N}_{\rm P}^+',
coord_def={XP: R>0})
NP1
with assuming(R>0):
for v in J[:-1]:
show(v.display(XP))
with assuming(R>0):
for v in K[:-1]:
show(v.display(XP))
graph_J1 = J1.plot(chart=XC, ambient_coords=(ch, ta), chart_domain=XC,
fixed_coords={th: pi/2, ps: 0}, ranges={ta: (-pi, 2*pi)},
number_values={ta: 16, ch: 7},
color='purple', scale=0.4, arrowsize=2)
graph = graph_XP + graph_J1
show(graph, figsize=8)
graph_K1 = K1.plot(chart=XC, ambient_coords=(ch, ta), chart_domain=XC,
fixed_coords={th: pi/2, ps: 0}, ranges={ta: (-pi, 2*pi)},
number_values={ta: 16, ch: 7},
color='blue', scale=0.4, arrowsize=2)
graph = graph_XP + graph_K1
show(graph, figsize=8)
taf = (Conf_to_Poinc.inverse()(T, R, 0, 0)[0]).function(T, R)
taf
chf = (Conf_to_Poinc.inverse()(T, R, 0, 0)[1]).function(T, R)
chf
def plot_p(n, T_values=None, R_min=0.0001, R_max=100,
color_T='dimgray', linestyle_T='-',
R_values=None, T_min=-100, T_max=100,
color_R='red', linestyle_R='-'):
n2 = 2*n
res = polygon([(pi/2, (n2 - 1)*pi), (-pi/2, n2*pi), (pi/2, (n2 + 1)*pi)],
color='white', edgecolor='black')
res += line([(pi/2, (n2 - 1)*pi), (-pi/2, n2*pi), (pi/2, (n2 + 1)*pi)],
color='black', thickness=3)
if T_values is not None:
for T0 in T_values:
res += parametric_plot((chf(T0, R), taf(T0, R) + n2*pi), (R, R_min, R_max),
color=color_T, linestyle=linestyle_T)
if R_values is not None:
for R0 in R_values:
res += parametric_plot((chf(T, R0), taf(T, R0) + n2*pi), (T, T_min, T_max),
color=color_R, linestyle=linestyle_R)
return res
def plot_m(n, T_values=None, R_min=-100, R_max=-0.0001,
color_T='dimgray', linestyle_T='-',
R_values=None, T_min=-100, T_max=100,
color_R='red', linestyle_R='-'):
n2 = 2*n
res = polygon([(-pi/2, n2*pi), (pi/2, (n2 + 1)*pi), (-pi/2, (n2 + 2)*pi)],
color='cornsilk', edgecolor='black')
res += line([(-pi/2, n2*pi), (pi/2, (n2 + 1)*pi), (-pi/2, (n2 + 2)*pi)],
color='black', thickness=3)
if T_values is not None:
for T0 in T_values:
res += parametric_plot((chf(T0, R), taf(T0, R) + n2*pi), (R, R_min, R_max),
color=color_T, linestyle=linestyle_T)
if R_values is not None:
for R0 in R_values:
res += parametric_plot((chf(T, R0), taf(T, R0) + n2*pi), (T, T_min, T_max),
color=color_R, linestyle=linestyle_R)
return res
R_val_p = [0.2, 0.4, 0.6, 0.8, 1]
R_val_m = [-1, -0.8, -0.6, -0.4, -0.2]
T_val = [-2, -1.5, -1, -0.5, 0, 0.5, 1, 1.5, 2]
graphics_array([plot_p(0, T_values=T_val, R_max=R_val_p[-1], R_values=R_val_p),
plot_m(0, T_values=T_val, R_min=R_val_m[0], R_values=R_val_m)])
graph = Graphics()
for n in [-1, 0, 1]:
graph += plot_p(n, T_values=T_val, R_max=R_val_p[-1], R_values=R_val_p)
graph += plot_m(n, T_values=T_val, R_min=R_val_m[0], R_values=R_val_m)
graph += line([(-pi/2, -4), (-pi/2, 8)], color='darkmagenta', thickness=2) \
+ line([(pi/2, -4), (pi/2, 8)], color='darkmagenta', thickness=2)
graph.axes_labels([r'$\chi$', r'$\tau$'])
show(graph, ymin=-pi, ymax=2*pi, frame=True, axes=False, figsize=10)
Adding some labels:
graph += text(r"$R=\!+\infty$", (pi/2+0.25, 0), rotation=90, color='darkmagenta',
fontsize=16) \
+ text(r"$R=\!-\infty$", (-pi/2-0.2, pi), rotation=90, color='darkmagenta',
fontsize=16) \
+ text(r"$R=1$", (0.75, 1), rotation=55, color='red', fontsize=16) \
+ text(r"$R=\!-1$", (-0.8, 2), rotation=55, color='red', fontsize=16) \
+ text(r"$R=0$", (-0.15, pi/2+0.15), rotation=45, color='black', fontsize=16) \
+ text(r"$\mathscr{N}_{\rm P}^+$", (0.9, 0), color='black', fontsize=18) \
+ text(r"$\mathscr{N}_{\rm P}^-$", (-0.9, pi), color='black', fontsize=18)
graph.axes_labels([r'$\chi$', r'$\tau$'])
graph.save("exk_NHEK_spacetime.pdf", xmin=-2.2, xmax=2.2, ymin=-pi, ymax=2*pi,
frame=True, axes=False, figsize=9)
show(graph, xmin=-2.2, xmax=2.2, ymin=-pi, ymax=2*pi, frame=True, axes=False, figsize=9)