CoCalc Shared FilesBHLectures / sage / Einstein-Rosen_bridge.ipynbOpen in CoCalc with one click!
Author: Eric Gourgoulhon
Views : 33

Einstein-Rosen bridge

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

Click here to download the worksheet file (ipynb format). To run it, you must start SageMath with the Jupyter notebook, with the command sage -n jupyter

NB: a version of SageMath at least equal to 7.5 is required to run this worksheet:

In [1]:
version()
'SageMath version 8.0.beta6, Release Date: 2017-05-12'

First we set up the notebook to display mathematical objects using LaTeX formatting:

In [2]:
%display latex

We also define a viewer for 3D plots (use 'threejs' or 'jmol' for interactive 3D graphics):

In [3]:
viewer3D = 'threejs' # must be 'threejs', jmol', 'tachyon' or None (default)

The rescaled Lambert function:

In [4]:
W0(x) = lambert_w(x/RDF(e)) + 1

The minimal value of the Kruskal-Szekeres coordinate XX on a hypersurface T=constT=\mathrm{const}:

In [5]:
def Xmin(T): return sqrt(max(0, RDF(T^2-1)))

Plot of rr as a function of XX on a hypersurface T=T0T=T_0:

In [6]:
T0_list = [0., 0.5, 1., 1.5, 2.0001] g = Graphics() for T0 in T0_list: alp = 0.2 + 0.4*(2. - T0) g += plot(2*W0(x^2-T0^2), (x,-4,-Xmin(T0)), thickness=2, color='blue', alpha=alp, legend_label="$T_0={:02.1f}$".format(float(T0))) + \ plot(2*W0(x^2-T0^2), (x, Xmin(T0), 4), thickness=2, color='blue', alpha=alp) show(g, aspect_ratio=1, axes_labels=[r'$X$', r'$r/m$'])
In [7]:
g.save("max_SigmaT0_r_X.pdf", aspect_ratio=1, axes_labels=[r'$X$', r'$r/m$'])

Plot of hypersurfaces T=T0T=T_0 in the Kruskal diagram

We generate first plots for the curvature singularity and the bifurcate horizon:

In [8]:
sing = plot(sqrt(1+x^2), (x,-3,3), color='brown', thickness=4, linestyle='--') \ + text(r"$r=0$", (2.5, 3), rotation=45, fontsize=16, color='brown') \ + plot(-sqrt(1+x^2), (x,-3,3), color='brown', thickness=4, linestyle='--') \ + text(r"$r'=0$", (2.5, -3), rotation=-45, fontsize=16, color='brown') bifhor = line([(-3,-3), (3,3)], color='black', thickness=3) + \ line([(-3,3), (3,-3)], color='black', thickness=3) + \ text(r'$\mathscr{H}$', (3, 2.7), fontsize=20, color='black')

We add the inner limits of the isometring embeddings for T0>1|T_0|>1:

In [9]:
graph = sing var('T') Xemb(T) = abs(T)*sqrt(2*ln(abs(T))) graph += parametric_plot([Xemb(T), T], (T, 1.001, 2.5), color='grey', thickness=2, linestyle=':') \ + parametric_plot([-Xemb(T), T], (T, 1.001, 2.5), color='grey', thickness=2, linestyle=':') \ + parametric_plot([Xemb(T), T], (T, -2.5, - 1.001), color='grey', thickness=2, linestyle=':') \ + parametric_plot([-Xemb(T), T], (T, -2.5, - 1.001), color='grey', thickness=2, linestyle=':') show(graph, aspect_ratio=1, axes_labels=[r'$X$', r'$T$'])

The minimal value of the Kruskal-Szekeres coordinate XX on an embedded surface of constant TT:

In [10]:
def Xmin_emb(T): if abs(T) > 1: return Xemb(T) return 0

The plot of the hypersurfaces T=T0T=T_0 for selected values of T0T_0:

In [11]:
colorS = 'steelblue' colorS3 = 'lightsteelblue' #colorS = 'lightseagreen' #colorS = 'wheat' T0_list = [-2., -1.5, -1., -0.5, 0., 0.5, 1., 1.5, 2.0001] for T0 in T0_list: graph += line([(-3, T0), (-Xmin_emb(T0), T0)], color=colorS, thickness=2) \ + line([(-Xmin_emb(T0), T0), (-Xmin(T0), T0)], color=colorS, thickness=3, linestyle=':') \ + line([(Xmin(T0), T0), (Xmin_emb(T0), T0)], color=colorS, thickness=3, linestyle=':') \ + line([(Xmin_emb(T0), T0), (3, T0)], color=colorS, thickness=2) \ + text("$T_0={:02.1f}$".format(float(T0)), (3.5, T0), fontsize=12, color=colorS) graph += bifhor show(graph, aspect_ratio=1, xmin=-3, xmax=3, axes_labels=[r'$X$', r'$T$'])
In [12]:
graph.save('max_constant_T_slices.pdf', aspect_ratio=1, xmin=-3, xmax=3, axes_labels=[r'$X$', r'$T$'])

Isometric embeddings

The function Z(r)Z'(r):

In [13]:
var('r') zp(r,T) = sqrt((1-T^2*e^(-r/2))/(r/2-1+T^2*e^(-r/2))) zp
(r,T)  2T2e(12r)12T2e(12r)+r2\left( r, T \right) \ {\mapsto} \ \sqrt{2} \sqrt{-\frac{T^{2} e^{\left(-\frac{1}{2} \, r\right)} - 1}{2 \, T^{2} e^{\left(-\frac{1}{2} \, r\right)} + r - 2}}

The exact integral is not known in the general case:

In [14]:
z = integrate(zp(r,T), r) z
2T2e(12r)12T2e(12r)+r2dr\sqrt{2} \int \sqrt{-\frac{T^{2} e^{\left(-\frac{1}{2} \, r\right)} - 1}{2 \, T^{2} e^{\left(-\frac{1}{2} \, r\right)} + r - 2}}\,{d r}

The minimal value of rr on an embedded surface of constant TT:

In [15]:
def rmin(T): if T > 1: return RDF(4*ln(T)) return 2*W0(-T^2)
In [16]:
plot(rmin, (T,0,2), axes_labels=[r'$T_0$', r'$r_{\rm min}(T_0)/m$'], thickness=2)
In [17]:
rmin(1/2), rmin(3/2)
(1.79634313781,1.62186043243)\left(1.79634313781, 1.62186043243\right)

The integration of Z(r)Z'(r) to get Z(r)Z(r) is performed numerically, using an algorithm for an adaptive integration with (integrable) singularities (qags):

In [18]:
def zz(r1, T0): dzdr = zp(r, T0) numint = numerical_integral(dzdr, rmin(T0), r1, algorithm='qags') error = numint[1] if error > 1e-3: print("Warning: error = {}".format(error)) if T0 > 1: return numint[0] + 1 return numint[0]

Test that it works:

In [19]:
zz(4, 1/2), zz(4, 3/2)
(4.22684717142,2.4583755882)\left(4.22684717142, 2.4583755882\right)

Function Z(r)Z(r) for some selected values of T0T_0:

In [20]:
T0_list = [0., 0.5, 0.9, 1., 1.5, 2.] g = Graphics() for T0 in T0_list: g += plot(lambda r: zz(r, T0), rmin(T0), 20) show(g, aspect_ratio=1, xmin=0, xmax=20, axes_labels=[r'$r/m$', r'$z/m$'])

3D plots of the embeddings

In [21]:
var('ph', latex_name=r'\phi', domain='real')
ϕ{\phi}
In [22]:
from sage.manifolds.utilities import set_axes_labels
In [23]:
T0_list = [0., 0.5, 0.9, 1., 1.5, 2.] rmax = 8 # in units of m for T0 in T0_list: g1 = parametric_plot3d([lambda r,ph: r*cos(ph), lambda r,ph: r*sin(ph), lambda r,ph: zz(r, T0)], (rmin(T0), rmax), (0, 2*pi), color=colorS3) g2 = parametric_plot3d([lambda r,ph: r*cos(ph), lambda r,ph: r*sin(ph), lambda r,ph: -zz(r, T0)], (rmin(T0), rmax), (0, 2*pi), color=colorS3) g = set_axes_labels(g1 + g2, 'x', 'y', 'z') print("T_0 = {:02.1f}".format(float(T0))) show(g, aspect_ratio=1, viewer=viewer3D, online=True) # show(g, aspect_ratio=1, viewer='tachyon', frame=False, figsize=20) g.save_image("max_embedding_T0_{:02.1f}.png".format(float(T0)), aspect_ratio=1, figsize=10, viewer='jmol')
T_0 = 0.0
T_0 = 0.5
T_0 = 0.9