Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News Sign UpSign In
| Download
Project: BHLectures
Views: 51
Kernel: SageMath 9.6.rc2

Vaidya spacetime

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

The computations make use of tools developed through the SageManifolds project.

version()
'SageMath version 9.6.rc2, Release Date: 2022-04-25'

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

%display latex

Spacetime

We declare the spacetime manifold MM:

M = Manifold(4, 'M', structure='Lorentzian') print(M)
4-dimensional Lorentzian manifold M

We introduce coordinates (v,r,θ,φ)(v,r,\theta,\varphi) analogous to the ingoing null Eddington-Finkelstein coordinates in Schwarzschild spacetime, i.e. such that vv is constant along ingoing radial null geodesics:

XN.<v,r,th,ph> = M.chart(r'v r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\varphi:periodic') XN

(M,(v,r,θ,φ))\displaystyle \left(M,(v, r, {\theta}, {\varphi})\right)

XN.coord_range()

v: (,+);r: (0,+);θ: (0,π);φ: [0,2π](periodic)\displaystyle v :\ \left( -\infty, +\infty \right) ;\quad r :\ \left( 0 , +\infty \right) ;\quad {\theta} :\ \left( 0 , \pi \right) ;\quad {\varphi} :\ \left[ 0 , 2 \, \pi \right] \mbox{(periodic)}

Metric tensor

The metric tensor corresponding to the Vaidya solution is:

m = function('m') g = M.metric() g[0,0] = -(1 - 2*m(v)/r) g[0,1] = 1 g[2,2] = r^2 g[3,3] = (r*sin(th))^2 g.display()

g=(2m(v)r1)dvdv+dvdr+drdv+r2dθdθ+r2sin(θ)2dφdφ\displaystyle g = \left( \frac{2 \, m\left(v\right)}{r} - 1 \right) \mathrm{d} v\otimes \mathrm{d} v +\mathrm{d} v\otimes \mathrm{d} r +\mathrm{d} r\otimes \mathrm{d} v + r^{2} \mathrm{d} {\theta}\otimes \mathrm{d} {\theta} + r^{2} \sin\left({\theta}\right)^{2} \mathrm{d} {\varphi}\otimes \mathrm{d} {\varphi}

g.inverse().display()

g1=vr+rv+(r2m(v)r)rr+1r2θθ+1r2sin(θ)2φφ\displaystyle g^{-1} = \frac{\partial}{\partial v }\otimes \frac{\partial}{\partial r }+\frac{\partial}{\partial r }\otimes \frac{\partial}{\partial v } + \left( \frac{r - 2 \, m\left(v\right)}{r} \right) \frac{\partial}{\partial r }\otimes \frac{\partial}{\partial r } + \frac{1}{r^{2}} \frac{\partial}{\partial {\theta} }\otimes \frac{\partial}{\partial {\theta} } + \frac{1}{r^{2} \sin\left({\theta}\right)^{2}} \frac{\partial}{\partial {\varphi} }\otimes \frac{\partial}{\partial {\varphi} }

g.inverse()[:]

(01001r2m(v)r00001r200001r2sin(θ)2)\displaystyle \left(\begin{array}{rrrr} 0 & 1 & 0 & 0 \\ 1 & \frac{r - 2 \, m\left(v\right)}{r} & 0 & 0 \\ 0 & 0 & \frac{1}{r^{2}} & 0 \\ 0 & 0 & 0 & \frac{1}{r^{2} \sin\left({\theta}\right)^{2}} \end{array}\right)

Curvature

The Ricci tensor is

Ric = g.ricci() Ric.display()

Ric(g)=2mvr2dvdv\displaystyle \mathrm{Ric}\left(g\right) = \frac{2 \, \frac{\partial\,m}{\partial v}}{r^{2}} \mathrm{d} v\otimes \mathrm{d} v

It has zero trace, i.e. the Ricci scalar vanishes:

g.ricci_scalar().expr()

0\displaystyle 0

The Riemann tensor:

Riem = g.riemann() Riem.display_comp(only_nonredundant=True)

Riem(g)vvvrvvvr=2m(v)r3Riem(g)vθvθvθvθ=m(v)rRiem(g)vφvφvφvφ=m(v)sin(θ)2rRiem(g)rvvrrvvr=2(rm(v)2m(v)2)r4Riem(g)rrvrrrvr=2m(v)r3Riem(g)rθvθrθvθ=mvRiem(g)rθrθrθrθ=m(v)rRiem(g)rφvφrφvφ=sin(θ)2mvRiem(g)rφrφrφrφ=m(v)sin(θ)2rRiem(g)θvvθθvvθ=r2mv+rm(v)2m(v)2r4Riem(g)θvrθθvrθ=m(v)r3Riem(g)θrvθθrvθ=m(v)r3Riem(g)θφθφθφθφ=2m(v)sin(θ)2rRiem(g)φvvφφvvφ=r2mv+rm(v)2m(v)2r4Riem(g)φvrφφvrφ=m(v)r3Riem(g)φrvφφrvφ=m(v)r3Riem(g)φθθφφθθφ=2m(v)r\displaystyle \begin{array}{lcl} \mathrm{Riem}\left(g\right)_{ \phantom{\, v} \, v \, v \, r }^{ \, v \phantom{\, v} \phantom{\, v} \phantom{\, r} } & = & \frac{2 \, m\left(v\right)}{r^{3}} \\ \mathrm{Riem}\left(g\right)_{ \phantom{\, v} \, {\theta} \, v \, {\theta} }^{ \, v \phantom{\, {\theta}} \phantom{\, v} \phantom{\, {\theta}} } & = & -\frac{m\left(v\right)}{r} \\ \mathrm{Riem}\left(g\right)_{ \phantom{\, v} \, {\varphi} \, v \, {\varphi} }^{ \, v \phantom{\, {\varphi}} \phantom{\, v} \phantom{\, {\varphi}} } & = & -\frac{m\left(v\right) \sin\left({\theta}\right)^{2}}{r} \\ \mathrm{Riem}\left(g\right)_{ \phantom{\, r} \, v \, v \, r }^{ \, r \phantom{\, v} \phantom{\, v} \phantom{\, r} } & = & \frac{2 \, {\left(r m\left(v\right) - 2 \, m\left(v\right)^{2}\right)}}{r^{4}} \\ \mathrm{Riem}\left(g\right)_{ \phantom{\, r} \, r \, v \, r }^{ \, r \phantom{\, r} \phantom{\, v} \phantom{\, r} } & = & -\frac{2 \, m\left(v\right)}{r^{3}} \\ \mathrm{Riem}\left(g\right)_{ \phantom{\, r} \, {\theta} \, v \, {\theta} }^{ \, r \phantom{\, {\theta}} \phantom{\, v} \phantom{\, {\theta}} } & = & \frac{\partial\,m}{\partial v} \\ \mathrm{Riem}\left(g\right)_{ \phantom{\, r} \, {\theta} \, r \, {\theta} }^{ \, r \phantom{\, {\theta}} \phantom{\, r} \phantom{\, {\theta}} } & = & -\frac{m\left(v\right)}{r} \\ \mathrm{Riem}\left(g\right)_{ \phantom{\, r} \, {\varphi} \, v \, {\varphi} }^{ \, r \phantom{\, {\varphi}} \phantom{\, v} \phantom{\, {\varphi}} } & = & \sin\left({\theta}\right)^{2} \frac{\partial\,m}{\partial v} \\ \mathrm{Riem}\left(g\right)_{ \phantom{\, r} \, {\varphi} \, r \, {\varphi} }^{ \, r \phantom{\, {\varphi}} \phantom{\, r} \phantom{\, {\varphi}} } & = & -\frac{m\left(v\right) \sin\left({\theta}\right)^{2}}{r} \\ \mathrm{Riem}\left(g\right)_{ \phantom{\, {\theta}} \, v \, v \, {\theta} }^{ \, {\theta} \phantom{\, v} \phantom{\, v} \phantom{\, {\theta}} } & = & -\frac{r^{2} \frac{\partial\,m}{\partial v} + r m\left(v\right) - 2 \, m\left(v\right)^{2}}{r^{4}} \\ \mathrm{Riem}\left(g\right)_{ \phantom{\, {\theta}} \, v \, r \, {\theta} }^{ \, {\theta} \phantom{\, v} \phantom{\, r} \phantom{\, {\theta}} } & = & \frac{m\left(v\right)}{r^{3}} \\ \mathrm{Riem}\left(g\right)_{ \phantom{\, {\theta}} \, r \, v \, {\theta} }^{ \, {\theta} \phantom{\, r} \phantom{\, v} \phantom{\, {\theta}} } & = & \frac{m\left(v\right)}{r^{3}} \\ \mathrm{Riem}\left(g\right)_{ \phantom{\, {\theta}} \, {\varphi} \, {\theta} \, {\varphi} }^{ \, {\theta} \phantom{\, {\varphi}} \phantom{\, {\theta}} \phantom{\, {\varphi}} } & = & \frac{2 \, m\left(v\right) \sin\left({\theta}\right)^{2}}{r} \\ \mathrm{Riem}\left(g\right)_{ \phantom{\, {\varphi}} \, v \, v \, {\varphi} }^{ \, {\varphi} \phantom{\, v} \phantom{\, v} \phantom{\, {\varphi}} } & = & -\frac{r^{2} \frac{\partial\,m}{\partial v} + r m\left(v\right) - 2 \, m\left(v\right)^{2}}{r^{4}} \\ \mathrm{Riem}\left(g\right)_{ \phantom{\, {\varphi}} \, v \, r \, {\varphi} }^{ \, {\varphi} \phantom{\, v} \phantom{\, r} \phantom{\, {\varphi}} } & = & \frac{m\left(v\right)}{r^{3}} \\ \mathrm{Riem}\left(g\right)_{ \phantom{\, {\varphi}} \, r \, v \, {\varphi} }^{ \, {\varphi} \phantom{\, r} \phantom{\, v} \phantom{\, {\varphi}} } & = & \frac{m\left(v\right)}{r^{3}} \\ \mathrm{Riem}\left(g\right)_{ \phantom{\, {\varphi}} \, {\theta} \, {\theta} \, {\varphi} }^{ \, {\varphi} \phantom{\, {\theta}} \phantom{\, {\theta}} \phantom{\, {\varphi}} } & = & -\frac{2 \, m\left(v\right)}{r} \end{array}

The Kretschmann scalar K=RabcdRabcdK = R_{abcd} R^{abcd}:

K = Riem.down(g)['_{abcd}'] * Riem.up(g)['^{abcd}'] K.expr()

48m(v)2r6\displaystyle \frac{48 \, m\left(v\right)^{2}}{r^{6}}

Wave vector kk

XN.coframe()

(M,(dv,dr,dθ,dφ))\displaystyle \left(M, \left(\mathrm{d} v,\mathrm{d} r,\mathrm{d} {\theta},\mathrm{d} {\varphi}\right)\right)

dv = XN.coframe()[0] dv.display()

dv=dv\displaystyle \mathrm{d} v = \mathrm{d} v

k = - dv.up(g) k.set_name('k') k.display()

k=r\displaystyle k = -\frac{\partial}{\partial r }

Check that 𝑘 is a null vector:

g(k, k).expr()

0\displaystyle 0

Check that kk is a geodesic vector field, i.e. fulfils kk=0\nabla_k k = 0:

nabla = g.connection() acc = nabla(k).contract(k) acc.display()

0\displaystyle 0

Homethetic Killing vector for m(v)=(m0/v0)vm(v) = (m_0/v_0) v:

xi = M.vector_field(v, r, 0, 0, name='xi', latex_name=r'\xi') xi.display()

ξ=vv+rr\displaystyle \xi = v \frac{\partial}{\partial v } + r \frac{\partial}{\partial r }

Lg = g.lie_derivative(xi) - 2*g Lg.display()

2(vmvm(v))rdvdv\displaystyle \frac{2 \, {\left(v \frac{\partial\,m}{\partial v} - m\left(v\right)\right)}}{r} \mathrm{d} v\otimes \mathrm{d} v

If m(v)m(v) is a linear function, the above result is identically zero, showing that Lξg=2g\mathcal{L}_{\xi} g = 2 g in that case, i.e. that ξ\xi is a homethetic Killing vector.

Ingoing Eddington-Finkelstein coordinates (t,r,θ,φ)(t,r,\theta,\varphi)

Let us introduce a new chart (t,r,θ,φ)(t,r,\theta,\varphi) such that the advanced time t+rt+r is vv: v=t+rv = t + r; this is the analog of ingoing Eddington-Finkelstein (IEF) coordinates in Schwarzschild spacetime.

X.<t,r,th,ph> = M.chart(r't r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\varphi:periodic') X

(M,(t,r,θ,φ))\displaystyle \left(M,(t, r, {\theta}, {\varphi})\right)

X.coord_range()

t: (,+);r: (0,+);θ: (0,π);φ: [0,2π](periodic)\displaystyle t :\ \left( -\infty, +\infty \right) ;\quad r :\ \left( 0 , +\infty \right) ;\quad {\theta} :\ \left( 0 , \pi \right) ;\quad {\varphi} :\ \left[ 0 , 2 \, \pi \right] \mbox{(periodic)}

We declare the transition map between the (t,r,θ,φ)(t,r,\theta,\varphi) and (v,r,θ,φ)(v,r,\theta,\varphi) coordinates:

X_to_XN = X.transition_map(XN, (t + r, r, th, ph)) X_to_XN.display()

{v=r+tr=rθ=θφ=φ\displaystyle \left\{\begin{array}{lcl} v & = & r + t \\ r & = & r \\ {\theta} & = & {\theta} \\ {\varphi} & = & {\varphi} \end{array}\right.

X_to_XN.inverse().display()

{t=r+vr=rθ=θφ=φ\displaystyle \left\{\begin{array}{lcl} t & = & -r + v \\ r & = & r \\ {\theta} & = & {\theta} \\ {\varphi} & = & {\varphi} \end{array}\right.

Expression of the metric tensor in the IEF coordinates:

g.display(X)

g=(r2m(r+t)r)dtdt+2m(r+t)rdtdr+2m(r+t)rdrdt+(r+2m(r+t)r)drdr+r2dθdθ+r2sin(θ)2dφdφ\displaystyle g = \left( -\frac{r - 2 \, m\left(r + t\right)}{r} \right) \mathrm{d} t\otimes \mathrm{d} t + \frac{2 \, m\left(r + t\right)}{r} \mathrm{d} t\otimes \mathrm{d} r + \frac{2 \, m\left(r + t\right)}{r} \mathrm{d} r\otimes \mathrm{d} t + \left( \frac{r + 2 \, m\left(r + t\right)}{r} \right) \mathrm{d} r\otimes \mathrm{d} r + r^{2} \mathrm{d} {\theta}\otimes \mathrm{d} {\theta} + r^{2} \sin\left({\theta}\right)^{2} \mathrm{d} {\varphi}\otimes \mathrm{d} {\varphi}

From now on, we set the IEF chart X to be the default one on MM:

M.set_default_chart(X) M.set_default_frame(X.frame())

Then g.display(X) can be substituted by g.display():

g.display()

g=(r2m(r+t)r)dtdt+2m(r+t)rdtdr+2m(r+t)rdrdt+(r+2m(r+t)r)drdr+r2dθdθ+r2sin(θ)2dφdφ\displaystyle g = \left( -\frac{r - 2 \, m\left(r + t\right)}{r} \right) \mathrm{d} t\otimes \mathrm{d} t + \frac{2 \, m\left(r + t\right)}{r} \mathrm{d} t\otimes \mathrm{d} r + \frac{2 \, m\left(r + t\right)}{r} \mathrm{d} r\otimes \mathrm{d} t + \left( \frac{r + 2 \, m\left(r + t\right)}{r} \right) \mathrm{d} r\otimes \mathrm{d} r + r^{2} \mathrm{d} {\theta}\otimes \mathrm{d} {\theta} + r^{2} \sin\left({\theta}\right)^{2} \mathrm{d} {\varphi}\otimes \mathrm{d} {\varphi}

g.inverse()[:]

(r+2m(r+t)r2m(r+t)r002m(r+t)rr2m(r+t)r00001r200001r2sin(θ)2)\displaystyle \left(\begin{array}{rrrr} -\frac{r + 2 \, m\left(r + t\right)}{r} & \frac{2 \, m\left(r + t\right)}{r} & 0 & 0 \\ \frac{2 \, m\left(r + t\right)}{r} & \frac{r - 2 \, m\left(r + t\right)}{r} & 0 & 0 \\ 0 & 0 & \frac{1}{r^{2}} & 0 \\ 0 & 0 & 0 & \frac{1}{r^{2} \sin\left({\theta}\right)^{2}} \end{array}\right)

xi.display()

ξ=tt+rr\displaystyle \xi = t \frac{\partial}{\partial t } + r \frac{\partial}{\partial r }

Einstein equation

The Ricci tensor in terms of the IEF coordinates:

Ric.display()

Ric(g)=2m(r+t)r2dtdt+2m(r+t)r2dtdr+2m(r+t)r2drdt+2m(r+t)r2drdr\displaystyle \mathrm{Ric}\left(g\right) = \frac{2 \, \frac{\partial\,m}{\partial \left( r + t \right)}}{r^{2}} \mathrm{d} t\otimes \mathrm{d} t + \frac{2 \, \frac{\partial\,m}{\partial \left( r + t \right)}}{r^{2}} \mathrm{d} t\otimes \mathrm{d} r + \frac{2 \, \frac{\partial\,m}{\partial \left( r + t \right)}}{r^{2}} \mathrm{d} r\otimes \mathrm{d} t + \frac{2 \, \frac{\partial\,m}{\partial \left( r + t \right)}}{r^{2}} \mathrm{d} r\otimes \mathrm{d} r

The notation m(r+t)\frac{\partial m}{\partial(r+t)} to denote dmdv\frac{\mathrm{d}m}{\mathrm{d}v} is quite unfortunate (this shall be improved in a future version). The display of the corresponding symbolic expression is slightly better, D0(m)\mathrm{D}_0(m) standing for the derivative of function mm with respect to its first (index 00) and unique argument, i.e. D0(m)=dmdv\mathrm{D}_0(m) = \frac{\mathrm{d}m}{\mathrm{d}v}:

Ric[0,0].expr()

2D0(m)(r+t)r2\displaystyle \frac{2 \, \mathrm{D}_{0}\left(m\right)\left(r + t\right)}{r^{2}}

The Ricci scalar is vanishing:

g.ricci_scalar().display()

r(g):MR(v,r,θ,φ)0(t,r,θ,φ)0\displaystyle \begin{array}{llcl} \mathrm{r}\left(g\right):& M & \longrightarrow & \mathbb{R} \\ & \left(v, r, {\theta}, {\varphi}\right) & \longmapsto & 0 \\ & \left(t, r, {\theta}, {\varphi}\right) & \longmapsto & 0 \end{array}

The energy-momentum vector ensuring that the Einstein equation is fulfilled is then:

T = 1/(8*pi)*Ric T.set_name('T') T.display()

T=m(r+t)4πr2dtdt+m(r+t)4πr2dtdr+m(r+t)4πr2drdt+m(r+t)4πr2drdr\displaystyle T = \frac{\frac{\partial\,m}{\partial \left( r + t \right)}}{4 \, \pi r^{2}} \mathrm{d} t\otimes \mathrm{d} t + \frac{\frac{\partial\,m}{\partial \left( r + t \right)}}{4 \, \pi r^{2}} \mathrm{d} t\otimes \mathrm{d} r + \frac{\frac{\partial\,m}{\partial \left( r + t \right)}}{4 \, \pi r^{2}} \mathrm{d} r\otimes \mathrm{d} t + \frac{\frac{\partial\,m}{\partial \left( r + t \right)}}{4 \, \pi r^{2}} \mathrm{d} r\otimes \mathrm{d} r

Since v=t+rv=t+r, we have dv=dt+dr\mathrm{d}v = \mathrm{d}t + \mathrm{d}r:

dv.display()

dv=dt+dr\displaystyle \mathrm{d} v = \mathrm{d} t+\mathrm{d} r

The derivative of the function m(v)m(v):

mp(v) = diff(m(v), v) mp(v)

vm(v)\displaystyle \frac{\partial}{\partial v}m\left(v\right)

T == 1/(4*pi)*mp(t+r)/r^2 * dv*dv

True\displaystyle \mathrm{True}

The future-directed null vector along the ingoing null geodesics:

k.display()

k=tr\displaystyle k = \frac{\partial}{\partial t }-\frac{\partial}{\partial r }

Outgoing radial null geodesics

Let us consider the vector field:

l = M.vector_field(1, (r - 2*m(t+r))/(r+2*m(t+r)), 0, 0, name='l', latex_name=r'\ell') l.display()

=t+(r2m(r+t)r+2m(r+t))r\displaystyle \ell = \frac{\partial}{\partial t } + \left( \frac{r - 2 \, m\left(r + t\right)}{r + 2 \, m\left(r + t\right)} \right) \frac{\partial}{\partial r }

It is a null vector:

g(l,l).display()

g(,):MR(v,r,θ,φ)0(t,r,θ,φ)0\displaystyle \begin{array}{llcl} g\left(\ell,\ell\right):& M & \longrightarrow & \mathbb{R} \\ & \left(v, r, {\theta}, {\varphi}\right) & \longmapsto & 0 \\ & \left(t, r, {\theta}, {\varphi}\right) & \longmapsto & 0 \end{array}

Moreover \ell is a pregeodesic vector field, i.e. it obeys =κ\nabla_\ell \ell = \kappa \ell:

acc = nabla(l).contract(l) acc.display()

(4(rm(r+t)m(r+t))r2+4rm(r+t)+4m(r+t)2)t+(4(r(2m(r+t)+1)m(r+t)r2m(r+t)2m(r+t)2)r3+6r2m(r+t)+12rm(r+t)2+8m(r+t)3)r\displaystyle \left( -\frac{4 \, {\left(r \frac{\partial\,m}{\partial \left( r + t \right)} - m\left(r + t\right)\right)}}{r^{2} + 4 \, r m\left(r + t\right) + 4 \, m\left(r + t\right)^{2}} \right) \frac{\partial}{\partial t } + \left( \frac{4 \, {\left(r {\left(2 \, \frac{\partial\,m}{\partial \left( r + t \right)} + 1\right)} m\left(r + t\right) - r^{2} \frac{\partial\,m}{\partial \left( r + t \right)} - 2 \, m\left(r + t\right)^{2}\right)}}{r^{3} + 6 \, r^{2} m\left(r + t\right) + 12 \, r m\left(r + t\right)^{2} + 8 \, m\left(r + t\right)^{3}} \right) \frac{\partial}{\partial r }

kappa = acc[0]/l[0] kappa

4(rm(r+t)m(r+t))r2+4rm(r+t)+4m(r+t)2\displaystyle -\frac{4 \, {\left(r \frac{\partial\,m}{\partial \left( r + t \right)} - m\left(r + t\right)\right)}}{r^{2} + 4 \, r m\left(r + t\right) + 4 \, m\left(r + t\right)^{2}}

acc == kappa*l

True\displaystyle \mathrm{True}

local/var/lib/sage/venv-python3.8/lib/python3.8/site-packages/scipy/integrate/## Integration of the outgoing radial null geodesics

The outgoing radial null geodesics are the field lines of \ell; they thus obey to drdt=rt \frac{\mathrm{d}r}{\mathrm{d}t} = \frac{\ell^r}{\ell^t}. Hence the value of drdt\frac{\mathrm{d}r}{\mathrm{d}t}:

drdt = (l[1] / l[0]).expr() drdt

r2m(r+t)r+2m(r+t)\displaystyle \frac{r - 2 \, m\left(r + t\right)}{r + 2 \, m\left(r + t\right)}

Choice of function m(v)m(v)

Let us choose a simple function m(v)m(v), based on one of the following smoothstep functions:

S0(x) = x S1(x) = -2*x^3 + 3*x^2 S2(x) = 6*x^5 - 15*x^4 + 10*x^3 S(x) = S0(x) #S(x) = S2(x)
m_0 = var('m_0') v_0 = var('v_0') h(v) = (1+sgn(v))/2 # the Heaviside function mS(v) = m_0*(h(v)*h(v_0 - v)*S(v/v_0) + h(v - v_0)) mS(v)

14(v(sgn(v)+1)(sgn(v+v0)+1)v0+2sgn(vv0)+2)m0\displaystyle \frac{1}{4} \, {\left(\frac{v {\left(\mathrm{sgn}\left(v\right) + 1\right)} {\left(\mathrm{sgn}\left(-v + v_{0}\right) + 1\right)}}{v_{0}} + 2 \, \mathrm{sgn}\left(v - v_{0}\right) + 2\right)} m_{0}

NB: we don't use Sage's predefined heaviside function, since it is incompatible with SciPy numerical integrators.

Case of singularity entirely hidden under the event horizon

The singularity is entirely hidden under the event horizon for large energy densities of the radiation shell, i.e. for small values of the shell width v0v_0. The precise criterion is v0<16mv_0 < 16 m for S(x)=S0(x)S(x) = S_0(x). We select here v0=3mv_0 = 3m:

v0 = 3
m_num(v) = mS(v).subs({m_0: 1, v_0: v0}) m_num(v)

112v(sgn(v)+1)(sgn(v+3)+1)+12sgn(v3)+12\displaystyle \frac{1}{12} \, v {\left(\mathrm{sgn}\left(v\right) + 1\right)} {\left(\mathrm{sgn}\left(-v + 3\right) + 1\right)} + \frac{1}{2} \, \mathrm{sgn}\left(v - 3\right) + \frac{1}{2}

graph = plot(m_num(v), (v, -v0, 2*v0), color='blue', legend_label=r'$S(x) = x$', thickness=3, axes_labels=[r'$v/m$', r'$M(v)/m$'], frame=True, axes=False, gridlines=True) graph += plot(S2(v/v0) , (v, 0, v0), color='red', legend_label=r'$S(x) = 6x^5 - 15x^4 + 10x^3$', thickness=3) graph += line([(v0, 0), (v0, 1)], linestyle='--', color='grey') graph += text(r'$v_0$', (v0, -0.08), color='black', fontsize=16) graph.set_axes_range(ymin=0) graph.save("vai_mass_function.pdf") graph
Image in a Jupyter notebook

Numerical integration of the outgoing radial null geodesics

We plug the function m(v)m(v) into the expression of drdt\frac{\mathrm{d}r}{\mathrm{d}t} along the outgoing radial null geodesics found above:

drdt1 = drdt.substitute_function(m, m_num) drdt1

(r+t)(sgn(r+t)+1)(sgn(rt+3)+1)6r+6sgn(r+t3)+6(r+t)(sgn(r+t)+1)(sgn(rt+3)+1)+6r+6sgn(r+t3)+6\displaystyle -\frac{{\left(r + t\right)} {\left(\mathrm{sgn}\left(r + t\right) + 1\right)} {\left(\mathrm{sgn}\left(-r - t + 3\right) + 1\right)} - 6 \, r + 6 \, \mathrm{sgn}\left(r + t - 3\right) + 6}{{\left(r + t\right)} {\left(\mathrm{sgn}\left(r + t\right) + 1\right)} {\left(\mathrm{sgn}\left(-r - t + 3\right) + 1\right)} + 6 \, r + 6 \, \mathrm{sgn}\left(r + t - 3\right) + 6}

and we perform a numerical integration:

fdrdt = fast_callable(drdt1, vars=[t, r])
from scipy.integrate import ode def solve_ode(t0, t1, r0=0.001, dt=0.02): t0 = RDF(t0) t1 = RDF(t1) r0 = RDF(r0) dt = RDF(dt) forward = 1 if t1 > t0 else -1 rd = ode(solve_ode.rhs) rd.set_initial_value(r0, t0) rd.set_integrator('dopri5') sol = [] ts = t0 rs = r0 while rd.successful() and forward*ts < forward*t1 and rs > 0: sol.append((ts, rs)) ts = rd.t + forward*dt rs = RDF(rd.integrate(ts)[0]) return sol solve_ode.rhs = fdrdt

Plot parameters:

tmax = v0 + 1.5 ymin = -4 ymax = tmax - 0.5 rmax = 6

Outgoing radial null geodesics from the Minkowski region:

def geod_line(sol): return line([(s[1], s[0]) for s in sol], color='green') def outgeods_from_Mink(t0, tmax, dt0, t0_max=0): geods = [] print('t0 : ', end='') while t0 < t0_max: sol = solve_ode(t0, tmax, r0=1e-8, dt=0.05) print(t0, end=', ') geods.append(geod_line(sol)) t0 += dt0 return geods
outgeods = outgeods_from_Mink(-10., tmax, 1., t0_max=-3.)
t0 : -10.0000000000000, -9.00000000000000, -8.00000000000000, -7.00000000000000, -6.00000000000000, -5.00000000000000, -4.00000000000000,
dt0 = 0.5 if S == S0 else 0.49 outgeods += outgeods_from_Mink(-3., tmax, dt0)
t0 : -3.00000000000000, -2.50000000000000, -2.00000000000000, -1.50000000000000, -1.00000000000000, -0.500000000000000,

Drawing of the radiation region (yellow rays):

def draw_radiation_region(v0, rmax, nl): graph = Graphics() for i in range(nl+1): t0 = float(i)/float(nl)*v0 graph += line([(0, t0), (rmax, t0 - rmax)], color='yellow') return graph
graph = draw_radiation_region(v0, rmax, 20)

Adding the outgoing radial null geodesics:

for geod in outgeods: graph += geod show(graph, aspect_ratio=1, xmax=rmax, ymin=ymin, ymax=ymax, axes_labels=[r'$r/m$', r'$t/m$'], figsize=10)
Image in a Jupyter notebook

The ingoing null geodesics:

def ingoing_geods(tmin, tmax, rmax, v0, dt_out, dt_in): geods = [] t0 = tmin while t0 < 0: geods.append(line([(0, t0), (rmax, t0 - rmax)], color='green', linestyle='--')) t0 += dt_out t0 = 0 while t0 <= v0: geods.append(line([(0, t0), (rmax, t0 - rmax)], color='green', linestyle='--')) t0 += dt_in t0 = v0 + dt_out while t0 < tmax + rmax: geods.append(line([(0, t0), (rmax, t0 - rmax)], color='green', linestyle='--')) t0 += dt_out return geods
ingeods = ingoing_geods(-6., tmax, rmax, v0, 1., 0.5) for geod in ingeods: graph += geod

The curvature singularity (in orange):

def curvature_sing(tmax, nb): h = tmax/(nb - 1) return line([(0, h*i) for i in range(nb)], thickness=3, color='darkorange', marker='D', markersize=8)
graph += curvature_sing(tmax, 21)

The event horizon (in black):

def event_hor(tmax, tmin=-4.): sol = solve_ode(tmax, tmin, r0=2.) hor = line([(s[1], s[0]) for s in sol], color='black', thickness=4) # Refinement near the horizon birth point: t0, r0 = sol[-2][0], sol[-2][1] sol = solve_ode(t0, tmin, r0=r0, dt=0.001) hor += line([(s[1], s[0]) for s in sol], color='black', thickness=4) thb = sol[-1][0] dthb = sol[-2][0] - thb print("Coordinate t at the event horizon birth [unit: m]: ") print(thb, " +/- ", dthb) return hor
graph += event_hor(tmax) pA = (2., v0 - 2.) # point A
Coordinate t at the event horizon birth [unit: m]: -2.600999999999992 +/- 0.0009999999999998899

Check of the determination of thbt_{\rm hb} by comparison with the analytic formula for S(x)=S0(x):=xS(x) = S_0(x) := x:

if S == S0: aa = 1/sqrt(16/v0 - 1) thb = numerical_approx(-4*exp(-2*aa*atan (aa))) print("Analytic value of coordinate t at the event horizon birth [unit: m]:") print(thb) pB = (-thb/2, thb/2) # point B pC = (0, thb) # point C
Analytic value of coordinate t at the event horizon birth [unit: m]: -2.60135096372454

The trapping horizon (in red):

def trapping_hor(mv, tmax): return parametric_plot((2*mv, v - 2*mv), (v, 0, tmax + 2), color='red', thickness=2)
graph += trapping_hor(m_num(v), tmax)
graph_wo_vectors = copy(graph)

The vectors kk and \ell at some point pp:

p = M((1, 4, pi/2, 0), chart=X) l.at(p).display()

=t+(m(5)2m(5)+2)r\displaystyle \ell = \frac{\partial}{\partial t } + \left( -\frac{m\left(5\right) - 2}{m\left(5\right) + 2} \right) \frac{\partial}{\partial r }

graph += k.at(p).plot(ambient_coords=(r, t), color='green', fontsize=16) graph += l.at(p).plot(ambient_coords=(r, t), color='green', fontsize=16, parameters={m(5): 1}, label_offset=0.15)
if S == S0: graph += circle(pA, 0.08, fill=True, color='black') \ + text(r"$A$", (pA[0]+0.2, pA[1]+0.1), fontsize=16, color='black') graph += circle(pB, 0.08, fill=True, color='black') \ + text(r"$B$", (pB[0]+0.25, pB[1]), fontsize=16, color='black') graph += circle(pC, 0.08, fill=True, color='black') \ + text(r"$C$", (pC[0]+0.25, pC[1]), fontsize=16, color='black') show(graph, aspect_ratio=1, xmax=rmax, ymin=ymin, ymax=ymax, axes_labels=[r'$r/m$', r'$t/m$'], figsize=10)
Image in a Jupyter notebook
filename = "vai_diag_S0.pdf" if S == S0 else "vai_diag_S2.pdf" graph.save(filename, aspect_ratio=1, xmax=rmax, ymin=ymin, ymax=ymax, axes_labels=[r'$r/m$', r'$t/m$'], figsize=10)

A zoom on the trapping horizon in its dynamical part: notice that the "outgoing" null geodesics cross it with a vertical tangent, in agreement with the cross-sections of the trapping horizon being marginally trapped surfaces.

show(graph_wo_vectors, aspect_ratio=1, ymin=-0.8, ymax=1.2, xmax=2.5, axes_labels=[r'$r/m$', r'$t/m$'])
Image in a Jupyter notebook

Case of a naked singularity

The naked singularity is obtained for small energy densities of the radiation shell, i.e. for large values of the shell width v0v_0. The precise criterion is v0>16mv_0 > 16 m for S(x)=S0(x)S(x) = S_0(x). We select here v0=18mv_0 = 18m:

v0 = 18
m_num(v) = mS(v).subs({m_0: 1, v_0: v0}) m_num(v)

172v(sgn(v)+1)(sgn(v+18)+1)+12sgn(v18)+12\displaystyle \frac{1}{72} \, v {\left(\mathrm{sgn}\left(v\right) + 1\right)} {\left(\mathrm{sgn}\left(-v + 18\right) + 1\right)} + \frac{1}{2} \, \mathrm{sgn}\left(v - 18\right) + \frac{1}{2}

drdt1 = drdt.substitute_function(m, m_num) fdrdt = fast_callable(drdt1, vars=[t, r]) solve_ode.rhs = fdrdt
#tmax = v0 + 4.5 tmax = v0 + 3.5 ymax = tmax - 0.5 rmax = 16

Outgoing radial null geodesics from the Minkowski region:

outgeods = outgeods_from_Mink(-rmax - 3, tmax, 2.)
t0 : -19, -17.0000000000000, -15.0000000000000, -13.0000000000000, -11.0000000000000, -9.00000000000000, -7.00000000000000, -5.00000000000000, -3.00000000000000, -1.00000000000000,

Outgoing radial null geodesics in the black hole region:

def outgeods_in_BH(tmax): geods = [] sol = solve_ode(tmax, -1, r0=1) geods.append(geod_line(sol)) sol = solve_ode(2/3*tmax, -1, r0=1) geods.append(geod_line(sol)) sol = solve_ode(2/3*tmax, tmax, r0=1) geods.append(geod_line(sol)) sol = solve_ode(1/3*tmax, -1, r0=0.5) geods.append(geod_line(sol)) sol = solve_ode(1/3*tmax, tmax, r0=0.5) geods.append(geod_line(sol)) return geods
outgeods += outgeods_in_BH(tmax)

The Cauchy horizon (in blue):

if S == S0: a0 = 2 / v0 rC = n(4/(1 - sqrt(1 - 8*a0))) tC = v0 - rC pC = (rC, tC) # point C sol = solve_ode(tC, -1, r0=rC) cauchy_hor = line([(s[1], s[0]) for s in sol], color='blue', thickness=2) sol = solve_ode(tC, tmax, r0=rC) cauchy_hor += line([(s[1], s[0]) for s in sol], color='blue', thickness=2) else: print("approximate CH") sol = solve_ode(-1.e-3, tmax, r0=1e-8) cauchy_hor = line([(s[1], s[0]) for s in sol], color='blue', thickness=2)

Outgoing radial null geodesics emerging from the initial singularity:

if S == S0: rB = n(4/(1 + sqrt(1 - 8*a0))) pB = (rB, v0 - rB) # point B else: rB = 2.5 rC = 6.5 nl = 5 dr = (rC - rB)/(nl - 1) print("r0 : ", end='') for i in range(nl): rr = rB + i*dr print(rr, end=', ') tt = v0 - rr sol = solve_ode(tt, -1, r0=rr) outgeods.append(geod_line(sol)) sol = solve_ode(tt, tmax, r0=rr) outgeods.append(geod_line(sol))
r0 : 3.00000000000000, 3.75000000000000, 4.50000000000000, 5.25000000000000, 6.00000000000000,

Drawing of the radiation region (yellow rays):

graph = draw_radiation_region(v0, rmax, 40)

Adding the outgoing radial null geodesics:

for geod in outgeods: graph += geod graph += cauchy_hor show(graph, aspect_ratio=1, xmax=rmax, ymin=-4, ymax=ymax, axes_labels=[r'$r/m$', r'$t/m$'], figsize=10)
Image in a Jupyter notebook
show(graph, aspect_ratio=1, ymin=0, ymax=1.4, xmax=1, axes_labels=[r'$r/m$', r'$t/m$'])
Image in a Jupyter notebook

The ingoing null geodesics:

ingeods = ingoing_geods(-6., tmax, rmax, v0, 2., 2.) for geod in ingeods: graph += geod

The curvature singularity (in orange):

graph += curvature_sing(tmax, 21)

The event horizon (in black):

graph += event_hor(tmax) pA = (2, v0 - 2) # point A
Coordinate t at the event horizon birth [unit: m]: 0.001000000000369725 +/- 0.001

The trapping horizon (in red):

graph += trapping_hor(m_num(v), tmax)
graph_wo_points = copy(graph) if S == S0: graph += circle((0,0), 0.2, fill=True, color='black') \ + text(r"$O$", (0.5, -0.6), fontsize=16, color='black') graph += circle(pA, 0.2, fill=True, color='black') \ + text(r"$A$", (pA[0]+0.7, pA[1]+0.6), fontsize=16, color='black') graph += circle(pB, 0.2, fill=True, color='green') \ + text(r"$B$", (pB[0]+0.8, pB[1]+0.2), fontsize=16, color='green') graph += circle(pC, 0.2, fill=True, color='blue') \ + text(r"$C$", (pC[0]+0.8, pC[1]+0.2), fontsize=16) show(graph, aspect_ratio=1, xmax=rmax, ymin=ymin, ymax=ymax, axes_labels=[r'$r/m$', r'$t/m$'], figsize=10)
Image in a Jupyter notebook
filename = "vai_diag_naked_S0.pdf" if S == S0 else "vai_diag_naked_S2.pdf" graph.save(filename, aspect_ratio=1, xmax=rmax, ymin=ymin, ymax=ymax, axes_labels=[r'$r/m$', r'$t/m$'], figsize=10)

A zoom on the initial singularity:

graph = graph_wo_points graph += circle((0,0), 0.05, fill=True, color='black') \ + text(r"$O$", (0.12, -0.15), fontsize=20, color='black') show(graph, aspect_ratio=1, ymin=-0.8, ymax=2, xmax=2, #show(graph, aspect_ratio=0.05, ymin=0.5, ymax=3, xmax=0.2, axes_labels=[r'$r/m$', r'$t/m$'])
Image in a Jupyter notebook