#!/usr/bin/env python # coding: utf-8 # # Naked singularity in Vaidya collapse # # This Jupyter/SageMath notebook is relative to the lectures # [Geometry and physics of black holes](https://relativite.obspm.fr/blackholes/). # # The computations make use of tools developed through the [SageManifolds project](https://sagemanifolds.obspm.fr). # In[1]: version() # In[2]: get_ipython().run_line_magic('display', 'latex') # In[3]: M = Manifold(4, 'M', structure='Lorentzian') print(M) # In[4]: XN. = M.chart(r'v:(0,+oo) r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\varphi:periodic') XN # ### The two roots $(x_1, x_2)$ of the polynomial $\alpha x^2 - x +2$ for $0<\alpha<\frac{1}{8}$ # In[5]: graph = plot((1 - sqrt(1 - 8*x))/(2*x), (x, 0.001, 1/8), color='blue', legend_label=r'$x_1$', axes_labels=[r'$\alpha$', '']) graph += plot((1 + sqrt(1 - 8*x))/(2*x), (x, 0.08, 1/8), color='red', legend_label=r'$x_2$') graph.show(frame=True, gridlines=True) # In[6]: x1 = var('x1', latex_name=r'x_1', domain='real') assume(x1 > 2, x1 < 4) x2 = var('x2', latex_name=r'x_2', domain='real') assume(x2 > 4) # In[7]: alpha = 2/(x1*x2) alpha # ### Vaidya metric in terms of the $(v, r, \theta,\varphi)$ coordinates # In[8]: g = M.metric() g[0,0] = -(1 - alpha*v/r) g[0,1] = 1 g[2,2] = r^2 g[3,3] = (r*sin(th))^2 g.display() # ## Double null-coordinate system in the region $r > v / x_2$ # In[9]: N1 = M.open_subset('N1', latex_name=r'N_1', coord_def={XN: r > v/x2}) # In[10]: X1. = N1.chart(r'u v:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\varphi:periodic') X1 # In[11]: assume(x2*r - v > 0) # In[12]: XN_to_X1 = XN.transition_map(X1, ((v/x1 - r)/((r - v/x2))^(x1/x2), v, th, ph)) XN_to_X1.display() # In[13]: g.display_comp(X1.frame(), chart=XN.restrict(N1)) # In[14]: g[X1.frame(),0,1].expr() # To simplify the components of $g$, let us substitute $x_2$ by its expression # in terms of $x_1$, i.e. $x_2 = \frac{2 x_1}{x_1 - 2}$: # In[15]: xx2 = 2*x1/(x1 - 2) g.apply_map(lambda x: x.subs({x2: xx2}).simplify_full(), frame=X1.frame(), chart=XN.restrict(N1), keep_other_components=True) g.display_comp(X1.frame(), chart=XN.restrict(N1)) # Alternative form of $g_{uv}$: # In[16]: guv = g[X1.frame(),0,1].expr() guv # In[17]: guv_alt = - x2/(x2 - x1)/r*(r - v/x2)^(x1/2) guv_alt # Test: # In[18]: s = guv - guv_alt.subs({x2: xx2}) s.simplify_full().canonicalize_radical() # Inverse metric in terms of the coordinates $(u,v,\theta,\varphi)$, taking into account the identity $x_2 = \frac{2 x_1}{x_1 - 2}$: # In[19]: g.inverse().apply_map(lambda x: x.subs({x2: xx2}).simplify_full(), frame=X1.frame(), chart=XN.restrict(N1), keep_other_components=True) g.inverse()[X1.frame(),:,XN.restrict(N1)] # We note that $g^{uu} = 0$ and $g^{vv} = 0$, which proves that the hypersurfaces of constant $u$ are null, as well as the hypersurfaces of constant $v$, i.e. that $(u,v,\theta,\varphi)$ is a **double-null coordinate system** on $N_1$. # ### Special case $\alpha = 1/9$ # In[20]: u_vr = XN_to_X1(v, r, th, ph)[0] u_vr # In[21]: u_vr1 = u_vr.subs({x1: 3, x2: 6}) u_vr1 # In[22]: u_vr1^2 # Solving for $r$ in terms of $(u,v)$: # In[23]: eq = u^2 == u_vr1^2 eq # In[24]: solve(eq, r) # In[25]: ruv = solve(eq, r)[0].rhs() ruv # Recovering Fig. 3a of B. Waugh and K. Lake, [Phys. Rev. D **34**, 2978 (1986)](https://doi.org/10.1103/PhysRevD.34.2978): # In[26]: contour_plot(ruv, (v, 0, 9), (u, -3, 9.5), cmap=['black'], contours=(0.1, 0.2, 0.5, 1., 5), fill=False, axes_labels=(r'$v$', r'$u$'), axes=True) # In[27]: plot3d(ruv, (u, -2, 9.5), (v, 0, 9), axes_labels=('u', 'v', 'r')) # In[28]: ruv.series(v, 3) # ## Double null-coordinate system in the region $r < v / x_1$ # In[29]: N2 = M.open_subset('N2', latex_name=r'N_2', coord_def={XN: r < v/x1}) # In[30]: X2. = N2.chart(r'u v:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\varphi:periodic') X2 # In[31]: assume(v - x1*r > 0) # In[32]: XN_to_X2 = XN.transition_map(X2, ((v/x2 - r)/((v/x1 - r))^(x2/x1), v, th, ph)) XN_to_X2.display() # In[33]: g.display_comp(X2.frame(), chart=XN.restrict(N2)) # In[34]: g[X2.frame(),0,1].expr() # To simplify the components of $g$, let us substitute $x_1$ by its expression # in terms of $x_2$, i.e. $x_1 = \frac{2 x_2}{x_2 - 2}$: # In[35]: xx1 = 2*x2/(x2 - 2) g.apply_map(lambda x: x.subs({x1: xx1}).simplify_full(), frame=X2.frame(), chart=XN.restrict(N2), keep_other_components=True) g.display_comp(X2.frame(), chart=XN.restrict(N2)) # Alternative form of $g_{uv}$: # In[36]: guv = g[X2.frame(),0,1].expr() guv # In[37]: guv_alt = - x1/(x2 - x1)/r*(v/x1 - r)^(x2/2) guv_alt # Test: # In[38]: s = guv - guv_alt.subs({x1: xx1}) s.simplify_full().canonicalize_radical() # Inverse metric in terms of the coordinates $(u,v,\theta,\varphi)$, taking into account the identity $x_1 = \frac{2 x_2}{x_2 - 2}$: # In[39]: g.inverse().apply_map(lambda x: x.subs({x1: xx1}).simplify_full(), frame=X2.frame(), chart=XN.restrict(N2), keep_other_components=True) g.inverse()[X2.frame(),:,XN.restrict(N2)] # We note that $g^{uu} = 0$ and $g^{vv} = 0$, which proves that the hypersurfaces of constant $u$ are null, as well as the hypersurfaces of constant $v$, i.e. that $(u,v,\theta,\varphi)$ is a **double-null coordinate system** on $N_2$. # ## Special case $\alpha = 1/9$ # In[40]: u_vr = XN_to_X2(v, r, th, ph)[0] u_vr # In[41]: u_vr1 = u_vr.subs({x1: 3, x2: 6}) u_vr1 # Solving for $r$ in terms of $(u,v)$: # In[42]: eq = u == u_vr1 solve(eq, r) # In[43]: ruv = solve(eq, r)[1].rhs() ruv # In[44]: ruv_alt = v/3 + (sqrt(1 - 2*u*v/3) - 1)/(2*u) ruv_alt # In[45]: (ruv - ruv_alt).simplify_full().canonicalize_radical() # In[46]: ruv_alt.series(u, 2) # In[47]: ruv_alt.series(v, 3) # Recovering Fig. 3b of B. Waugh and K. Lake, [Phys. Rev. D **34**, 2978 (1986)](https://doi.org/10.1103/PhysRevD.34.2978): # In[48]: contour_plot(ruv, (v, 0, 14), (u, -10, 5), cmap=['black'], contours=(0, 0.05, 0.1, 0.2, 0.5, 1., 2.), fill=False, axes_labels=(r'$v$', r'$u$'), axes=True) # In[49]: plot3d(ruv, (u, -10, 5), (v, 0, 14), axes_labels=('u', 'v', 'r'))