# Event and trapping horizons in Vaidya spacetime¶

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

These computations are based on SageManifolds (v0.9)

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

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

In [1]:
%display latex

## Spacetime¶

We declare the spacetime manifold $M$:

In [2]:
M = Manifold(4, 'M')
print(M)
4-dimensional differentiable manifold M

We use coordinates $(t,r,\theta,\phi)$ analogous to the 3+1 Eddington-Finkelstein coordinates in Schwarzschild spacetime, i.e. coordinates such that the advanced time $v = t+r$ is constant on the ingoing radial null geodesics:

In [3]:
X.<t,r,th,ph> = M.chart(r't r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\phi')
X
Out[3]:
In [4]:
X.coord_range()
Out[4]:

## Metric tensor¶

The metric tensor corresponding to the Vaidya solution is:

In [5]:
m = function('m')
g = M.lorentzian_metric('g')
g[0,0] = -(1-2*m(t+r)/r)
g[0,1] = 2*m(t+r)/r
g[1,1] = 1+2*m(t+r)/r
g[2,2] = r^2
g[3,3] = (r*sin(th))^2
g.display()
Out[5]:
In [6]:
nice_derivatives(False)

## Einstein equation¶

Let us compute the Ricci tensor associated with the metric $g$:

In [7]:
Ric = g.ricci()
Ric.display()
Out[7]:
In [8]:
Ric.display()
Out[8]:

The Ricci scalar is vanishing:

In [9]:
g.ricci_scalar().display()
Out[9]:

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

In [10]:
T = 1/(8*pi)*Ric
T.display()
Out[10]:

Since $v=t+r$, we have $\mathrm{d}v = \mathrm{d}t + \mathrm{d}r$:

In [11]:
dv = M.scalar_field({X: t+r}, name='v').differential()
dv.display()
Out[11]:

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

In [12]:
v = var('v')
mp(v) = diff(m(v),v)
mp(v)
Out[12]:
In [13]:
T == 1/(4*pi)*mp(t+r)/r^2 * dv*dv
Out[13]:

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

In [14]:
ks = -dv.up(g)
print(ks)
ks.display()
Vector field on the 4-dimensional differentiable manifold M
Out[14]:
In [15]:
g(ks,ks).display()
Out[15]:

Let us consider the vector field:

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

It is a null vector:

In [17]:
g(l,l).display()
Out[17]:

Moreover it is a pregeodesic vector field:

In [18]:
nab = g.connection()
acc = nab(l).contract(l)
acc.display()
Out[18]:
In [19]:
acc[0]/l[0]
Out[19]:
In [20]:
acc[1]/l[1]
Out[20]:
In [21]:
acc[0]/l[0] == acc[1]/l[1]
Out[21]:

## Integration of the outgoing radial null geodesics¶

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

In [22]:
drdt = (l[1] / l[0]).expr()
drdt
Out[22]:

Let us choose a simple function $m(v)$:

In [23]:
v0 = var('v0')
m0 = var('m0')
h(v) = (1+sgn(v))/2  # the Heaviside function
m1(v) = m0* (h(v)*h(v0-v)*v/v0 + h(v-v0))
m1(v)
Out[23]:
In [24]:
plot(m1(v).subs({m0: 1, v0: 3}), (v,-3, 6), thickness=3,
axes_labels=[r'$v/m_0$', r'$m(v)/m_0$'])
Out[24]:

We plug this function into the expression of $\frac{\mathrm{d}r}{\mathrm{d}t}$ found above:

In [25]:
drdt1 = drdt.substitute_function(m, m1)
drdt1
Out[25]:

and we perform a numerical integration for $v_0 = 3 m_0$

In [26]:
drdt0 = drdt1.subs({m0: 1, v0: 3})
drdt0
Out[26]:
In [27]:
outgeods = []
tmax = 4
for t0 in [-6, -5, -4, -3, -2.8, -2.5, -2, -1.5, -1]:
sol = desolve_rk4(drdt0, r, ivar=t, ics=[t0,0.01], end_points=[t0, tmax], step=0.02)
outgeods.append(line([[s[1],s[0]] for s in sol if s[1]>0], color='green'))
In [28]:
graph = Graphics()
rmax = 6
nl = 20
for i in range(nl+1):
t0 = float(i)/float(nl)*3
graph += line([(0,t0), (rmax, t0-rmax)], color='yellow')
In [29]:
for geod in outgeods:
graph += geod
show(graph, aspect_ratio=1, ymin=-4, axes_labels=[r'$r/m_0$', r'$t/m_0$'])

The event horizon:

In [30]:
t0 = -2.6
sol = desolve_rk4(drdt0, r, ivar=t, ics=[t0,0.01], end_points=[t0, tmax], step=0.02)
hor = line([[s[1],s[0]] for s in sol if s[1]>0], color='black', thickness=3)
graph += hor
In [31]:
ingeods = []
rmax = 6
for t0 in range(-6, 10):
ingeods.append(line([(0,t0), (rmax, t0-rmax)], color='green'))
for geod in ingeods:
graph += geod

The trapping horizon:

In [32]:
trap = line([(0,0), (2,1), (2,tmax)], color='red', thickness=2)
graph += trap
In [33]:
for t0 in [0.5, 1.5, 2.5]:
graph += line([(0,t0), (rmax, t0-rmax)], color='green')
show(graph, aspect_ratio=1, ymin=-4, ymax=4, axes_labels=[r'$r/m_0$', r'$t/m_0$'])
In [34]:
graph.save("vaidya.pdf", aspect_ratio=1, ymin=-4, ymax=4, axes_labels=[r'$r/m_0$', r'$t/m_0$'])

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.

In [35]:
show(graph, aspect_ratio=1, ymin=-0.8, ymax=1.2, xmax=2.5, axes_labels=[r'$r/m_0$', r'$t/m_0$'])
In [ ]: