To recap: to find a black hole apparent horizon we have to solve the boundary value problem
\begin{equation} \frac{d^2 h}{d \theta^2} = 2 h - \frac{\cot(\theta)}{C^2} \frac{d h}{d \theta} + \frac{4 h^2}{\psi C^2} \left( \frac{\partial \psi}{\partial r} - \frac{1}{h^2} \frac{\partial \psi}{\partial \theta} \frac{d h}{d \theta} \right) + \frac{3}{h} \left( \frac{d h}{d \theta} \right)^2, \end{equation}with boundary conditions
\begin{equation} \frac{d h}{d \theta} \left( \theta = 0 \right) = 0 = \frac{d h}{d \theta} \left( \theta = \pi \right), \end{equation}should a solution exist. If no solution exists, then there is no horizon. If multiple solutions exist then they do not cross and the AH is the solution with largest $h$ (which encodes the horizon radius $r$ at a particular angle $\theta$) for any angle.
There are many methods for solving ordinary differential equation boundary value problems such as this. The standard advice is to first use the shooting method.
In shooting we convert the boundary value problem to an initial value problem by providing a guess for the missing initial conditions. We then solve the initial value problem up to the boundary point. We can then use the value of the solution, compared to the boundary condition of the original problem, to modify the initial guess.
Explicitly, we write the problem above in the form of an initial value problem as
\begin{align} \frac{d}{d \theta} \begin{pmatrix} h \\ \frac{d h}{d \theta} \end{pmatrix} & = \frac{d}{d \theta} \begin{pmatrix} H_1 \\ H_2 \end{pmatrix} \\ & = \begin{pmatrix} \frac{d h}{d \theta} \\ 2 h - \frac{\cot(\theta)}{C^2} \frac{d h}{d \theta} + \frac{4 h^2}{\psi C^2} \left( \frac{\partial \psi}{\partial r} - \frac{1}{h^2} \frac{\partial \psi}{\partial \theta} \frac{d h}{d \theta} \right) + \frac{3}{h} \left( \frac{d h}{d \theta} \right)^2 \end{pmatrix} \\ & = \begin{pmatrix} H_2 \\ 2 H_1 - \frac{\cot(\theta)}{C^2} H_2 + \frac{4 H_1^2}{\psi C^2} \left( \frac{\partial \psi}{\partial r} - \frac{1}{H_1^2} \frac{\partial \psi}{\partial \theta} H_2 \right) + \frac{3}{H_1} H_2^2 \end{pmatrix} \end{align}with initial conditions
\begin{equation} h(\theta = 0) = H_1(\theta = 0) = h_0, \qquad \frac{d h}{d \theta}(\theta = 0) = H_2(\theta = 0) = 0. \end{equation}This defines an initial value problem for the state vector ${\bf H}$. The solution is a function of $\theta$ but depends on the guess for $h_0$, the horizon radius at the initial point $\theta = 0$. We therefore write the solution of the initial value problem as ${\bf H}(\theta ; h_0)$ to show that it is a function of $\theta$ depending on the parameter $h_0$.
Note that the additional terms in the initial value problem such as the auxilliary variable $C$ and the quantities encoding the spacetime ($\psi$ and its derivatives) are completely given in terms of $h$, $d h / d \theta$, and $\theta$. Therefore the initial value problem is in the form
\begin{equation} \frac{d}{d \theta} {\bf H} = {\bf F} \left( \theta, {\bf H} \right), \qquad {\bf H}(\theta = 0) = {\bf H}_0. \end{equation}There are many standard numerical methods for solving initial value problems. Within python the scipy library contains a number of routines. As a simple first test we will use scipy.integrate.odeint to solve the initial value problem.
Applying the initial value problem solvers we can construct ${\bf H}(\theta ; h_0)$. From the solution we can evaluate the solution at the endpoint: either at $\theta = \pi$ (if there is no symmetry in the problem) or at $\theta = \pi / 2$. This gives us, for example, ${\bf H}(\pi ; h_0)$, and most importantly gives us ${\bf H_2}(\pi ; h_0) = \frac{d h}{d \theta}(\pi ; h_0)$.
Now, the original boundary condition tells us that $\frac{d h}{d \theta}(\pi)$ should be identically zero. If our guess for $h_0$ is such that $\frac{d h}{d \theta}(\pi ; h_0) = 0$, then we have a solution of the boundary value problem. Of course, the chances of this happening by luck are minute. Instead we must repeatedly evaluate the solution for different values of $h_0$ in order to find a value that gives us a solution of the boundary value problem.
This step can be formalised in the following way. Define the function
\begin{equation} \phi(h_0) = \frac{d h}{d \theta}(\pi ; h_0). \end{equation}We are trying to find the value $\bar{h}_0$ such that $\phi(\bar{h}_0) = 0$. This is a standard nonlinear algebraic root-finding problem.
Again, there are many standard numerical methods for solving root-finding problems of this sort, and again python has some contained within scipy. In this first simple test we will use the recommended scipy.optimize.brentq method, which requires two initial guesses bracketing a single solution $\bar{h}_0$.
Once we have found a solution $\bar{h}_0$, then the apparent horizon is found by
converting to cartesian coordinates using
\begin{equation} x = h \sin(\theta), \qquad z = h \cos(\theta) \end{equation}
We can then plot the solution.
In principle it is necessary to find all solutions to the boundary value problem in order to locate the apparent horizon (which is the outermost solution). It can be very difficult to guarantee that you have truly found the apparent horizon just from the information given here.
In [1]:
%matplotlib inline
In [2]:
import numpy as np
import matplotlib
matplotlib.rcParams.update({'font.size': 14})
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.optimize import brentq
To implement the solver, we need the explicit formulas for the auxilliary variable $C$, which is defined as
\begin{equation} C = \left( 1 + \left( \frac{1}{h} \frac{d h}{d \theta} \right)^2 \right)^{-1/2}, \end{equation}and for the spacetime quantities
\begin{align} \psi & = 1 + \frac{1}{2} \sum_{i=1}^N \frac{m^{(i)}}{\left\| {\bf y} - {\bf y}^{(i)} \right\|_2}, \\ \frac{\partial \psi}{\partial r} & = \frac{1}{2} \sum_{i=1}^N \frac{m^{(i)} \left( z^{(i)} \cos(\theta) - r \right)}{\left\| {\bf y} - {\bf y}^{(i)} \right\|_2^3}, \\ \frac{\partial \psi}{\partial \theta} & = -\frac{1}{2} \sum_{i=1}^N \frac{m^{(i)} r z^{(i)} \sin(\theta)}{\left\| {\bf y} - {\bf y}^{(i)} \right\|_2^3}. \end{align}Note that we use coordinates $(r, \theta)$, where $(x, z) = (r \cos(\theta), r \sin(\theta))$, to match the symmetries of the spacetime. The $m^{(i)}$ are the bare masses of the singularities, and the ${\bf y}^{(i)} = (0, z^{(i)})$ are their locations (which, due to the axisymmetry, must lie on the $z$-axis).
We must note that the $\cot(\theta)$ term is problematic at $\theta = 0, \pi$. $\cot(\theta)$ blows up at those points. As $d h / d \theta$ should vanish there, it can be shown that the full term including $\cot(\theta)$ vanishes at those points. This must be imposed by hand, and can cause numerical problems.
First let us solve for the simplest case: one singularity of mass $1$. This situation corresponds to a Schwarzschild black hole - static and spherically symmetric. The singularity is located at the origin.
In this case we have
\begin{align} \psi & = 1 + \frac{1}{2 r}, \\ \frac{\partial \psi}{\partial r} & = -\frac{1}{2 r^2}, \\ \frac{\partial \psi}{\partial \theta} & = 0. \end{align}Thanks to the symmetry of the problem we can integrate from $\theta = 0$ to $\theta = \pi/2$, where the boundary condition $d h / d \theta = 0$ holds. This avoids the problem of trying to integrate up to $\theta = \pi$ where the coordinate singularity, indicated by the problem with the $\cot(\theta)$ term mentioned above, can cause trouble.
In [3]:
def horizon_function(H, theta):
"""Encoding the apparent horizon IVP equation."""
h = H[0]
dh = H[1]
psi = 1.0 + 1.0 / (2.0 * h)
dpsi_dr = -1.0 / (2.0 * h**2)
dpsi_dtheta = 0.0
C2 = 1.0 / (1.0 + (dh / h)**2)
# Impose that the term involving cot(theta) vanishes on axis.
if (abs(theta) < 1e-16) or (abs(theta - np.pi) < 1e-16):
cot_theta_dh_C2 = 0.0
else:
cot_theta_dh_C2 = dh / (np.tan(theta) * C2)
dHdtheta = np.zeros_like(H)
dHdtheta[0] = dh
dHdtheta[1] = 2.0*h - cot_theta_dh_C2 + 4.0*h**2/(psi*C2)*(dpsi_dr - dpsi_dtheta*dh/h**2) + 3.0*dh**2/h
return dHdtheta
In [4]:
def shooting_function(h0):
"""Encode the function phi in order to find the initial radius."""
theta = [0.0, np.pi / 2.0]
H0 = [h0, 0.0]
H = odeint(horizon_function, H0, theta, atol = 1.0e-6, rtol = 1.0e-4)
# Return dh/dtheta at pi / 2.
return H[-1, 1]
We now find the location of the horizon. The correct value is exactly $1/2$, but to test the algorithm we look for a solution in the range $h_0 \in [0.4, 0.7]$.
In [5]:
h0 = brentq(shooting_function, 0.4, 0.7)
print(r"Found horizon with radius at \theta=0 of {}".format(h0))
Now we can solve for the horizon and plot the solution.
In [6]:
theta = np.linspace(0.0, np.pi / 2.0)
H0 = [h0, 0.0]
H = odeint(horizon_function, H0, theta, atol = 1.0e-6, rtol = 1.0e-4)
h = H[:, 0]
x = h * np.sin(theta)
z = h * np.cos(theta)
fig = plt.figure(figsize = (8, 8))
ax = fig.add_subplot(111)
ax.plot(x, z)
ax.set_xlabel("$x$")
ax.set_ylabel("$z$")
ax.set_title("Horizon location for a Schwarzschild black hole, mass = 1")
ax.set_xbound(0.0, 0.55)
ax.set_ybound(0.0, 0.55)
plt.show()
Now look at the case where there are two singularities. For simplicity we will give these singularities unit mass $m^{(1,2)} = 1$ and locate them symmetrically across the axis, $z^{(1)} = -0.75$, $z^{(2)} = +0.75$. This means we still have reflection symmetry and only need to solve for $\theta \in [0, \pi/2]$.
We can use exactly the same techniques as in the Schwarzschild case.
In [7]:
def horizon_function_binary(H, theta):
"""Encoding the apparent horizon IVP equation for the binary case."""
h = H[0]
dh = H[1]
m_i = [1.0, 1.0]
z_i = [-0.75, 0.75]
distance_i = np.zeros_like(z_i)
for i in range(len(z_i)):
distance_i[i] = np.sqrt((h*np.sin(theta))**2 + (h*np.cos(theta) - z_i[i])**2)
psi = 1.0
dpsi_dr = 0.0
dpsi_dtheta = 0.0
for i in range(len(m_i)):
psi += 0.5*m_i[i]/distance_i[i]
dpsi_dr -= 0.5*m_i[i]*(h-z_i[i]*np.cos(theta))/distance_i[i] ** 3
dpsi_dtheta -= 0.5*m_i[i]*h*z_i[i]*np.sin(theta)/distance_i[i] ** 3
C2 = 1.0 / (1.0 + (dh / h)**2)
# Impose that the term involving cot(theta) vanishes on axis.
if (abs(theta) < 1e-16) or (abs(theta - np.pi) < 1e-16):
cot_theta_dh_C2 = 0.0
else:
cot_theta_dh_C2 = dh / (np.tan(theta) * C2)
dHdtheta = np.zeros_like(H)
dHdtheta[0] = dh
dHdtheta[1] = 2.0*h - cot_theta_dh_C2 + 4.0*h**2/(psi*C2)*(dpsi_dr - dpsi_dtheta*dh/h**2) + 3.0*dh**2/h
return dHdtheta
Note that the only difference for the shooting function is that the function it calls has changed from horizon_function to horizon_function_binary. This indicates that a general purpose routine should be constructed, but we ignore that for now.
In [8]:
def shooting_function_binary(h0):
"""Encode the function phi in order to find the initial radius."""
theta = [0.0, np.pi / 2.0]
H0 = [h0, 0.0]
H = odeint(horizon_function_binary, H0, theta, atol = 1.0e-8, rtol = 1.0e-6)
# Return dh/dtheta at pi / 2.
return H[-1, 1]
We can now attempt to find the solution. However, in order to do this we need a reasonable guess for the initial radius $h_0$. Given the solution from the Schwarzschild case we might expect a guess of $h_0 = 0.75 + 0.5 = 1.25$ to be reasonable if there was just one singularity. But there are two, of total mass $2$, suggesting the horizon might be as far out as $h_0 = 0.75 + 1.0 = 1.75$. But we can check by computing the solution for many different values of $h_0$ and just plotting the function $\phi$ whose root we are trying to find.
In [9]:
h0_range = np.linspace(1.2, 2.0, 100)
phi_range = np.zeros_like(h0_range)
for i in range(len(h0_range)):
phi_range[i] = shooting_function_binary(h0_range[i])
fig = plt.figure(figsize=(14,6))
ax = fig.add_subplot(121)
ax.plot(h0_range, phi_range)
ax.set_xlabel("$h_0$")
ax.set_ylabel("$\phi$")
ax.axhline(color = 'black', linestyle = '--')
ax = fig.add_subplot(122)
ax.plot(h0_range[:20], phi_range[:20])
ax.set_xlabel("$h_0$")
ax.set_ylabel("$\phi$")
ax.axhline(color = 'black', linestyle = '--')
plt.show()
We see that there are two solutions close together. The plot suggests that a reasonable range for the outer solution, corresponding to the apparent horizon, is $h_0 \in [1.25, 1.3]$.
In [10]:
h0_binary_outer = brentq(shooting_function_binary, 1.25, 1.3)
print(r"Found horizon with radius at \theta=0 of {}".format(h0_binary_outer))
We will also find the inner trapped surface for comparison.
In [11]:
h0_binary_inner = brentq(shooting_function_binary, 1.2, 1.25)
print(r"Found horizon with radius at \theta=0 of {}".format(h0_binary_inner))
Now we can solve for each trapped surface and plot.
In [12]:
theta = np.linspace(0.0, np.pi / 2.0)
H0 = [h0_binary_outer, 0.0]
H = odeint(horizon_function_binary, H0, theta, atol = 1.0e-6, rtol = 1.0e-4)
h = H[:, 0]
x_outer = h * np.sin(theta)
z_outer = h * np.cos(theta)
H0 = [h0_binary_inner, 0.0]
H = odeint(horizon_function_binary, H0, theta, atol = 1.0e-6, rtol = 1.0e-4)
h = H[:, 0]
x_inner = h * np.sin(theta)
z_inner = h * np.cos(theta)
fig = plt.figure(figsize = (8, 8))
ax = fig.add_subplot(111)
ax.plot(x_outer, z_outer, x_inner, z_inner)
ax.set_xlabel("$x$")
ax.set_ylabel("$z$")
ax.set_title("Horizon location for a symmetric binary, singularities at $z = \pm 0.75$.")
ax.set_xbound(0.0, 1.4)
ax.set_ybound(0.0, 1.4)
plt.show()
One question that could be asked is whether there is a common horizon at all, given symmetric singularities as above. By a common horizon we mean one that surrounds both singularities and so can be found using the simple algorithm above. It is intuitively clear that if the singularities are sufficiently far apart then each will be surrounded by an individual horizon, but not by the common horizon. This suggests the precise question: what is the value $z^*$ such that, for a spacetime with singularities of unit mass located at $\pm z$ there is a common horizon if $z < z^*$ but no common horizon if $z > z^*$?
This would appear to be another root finding problem. However, there is a difficulty: for sufficiently small separations the common horizon exists, but for large separations it does not. Thus there is no immediately obvious meaningful quantity associated with the horizon which can be used to find the root, as it will not be defined for large $z$.
However, the figure of $\phi$ against $h_0$ above suggests the following trick. The apparent horizon - the outermost trapped surface - is also associated with an inner trapped surface. As we move the singularities further apart, the surfaces come closer together. This moves the zeros of $\phi(h_0)$ closer together. When $\phi$ has no zeros, there is no solution to the boundary value problem, and hence no common horizon. At this point, the figure above suggests that the minimum value of $\phi$ will change from being negative to being positive. In other words, the critical value $z^*$ is such that the minimum value of $\phi$ is zero.
Formally, we define a new function
\begin{equation} \Phi(z) = \min_{h_0} \phi(h_0 ; z) \end{equation}where $\phi(h_0 ; z)$ encodes the shooting function given above, but explicitly notes that the precise value depends on the parameters, where the crucial parameter here is the separation between the singularities.
In order to test this, we need to generalize our functions a little.
In [13]:
def horizon_function_binary_general(H, theta, z):
"""Encoding the apparent horizon IVP equation for the binary case."""
h = H[0]
dh = H[1]
m_i = [1.0, 1.0]
z_i = [-z, z]
distance_i = np.zeros_like(z_i)
for i in range(len(z_i)):
distance_i[i] = np.sqrt((h*np.sin(theta))**2 + (h*np.cos(theta) - z_i[i])**2)
psi = 1.0
dpsi_dr = 0.0
dpsi_dtheta = 0.0
for i in range(len(m_i)):
psi += 0.5*m_i[i]/distance_i[i]
dpsi_dr -= 0.5*m_i[i]*(h-z_i[i]*np.cos(theta))/distance_i[i]**3
dpsi_dtheta -= 0.5*m_i[i]*h*z_i[i]*np.sin(theta)/distance_i[i]**3
C2 = 1.0 / (1.0 + (dh / h)**2)
# Impose that the term involving cot(theta) vanishes on axis.
if (abs(theta) < 1e-16) or (abs(theta - np.pi) < 1e-16):
cot_theta_dh_C2 = 0.0
else:
cot_theta_dh_C2 = dh / (np.tan(theta) * C2)
dHdtheta = np.zeros_like(H)
dHdtheta[0] = dh
dHdtheta[1] = 2.0*h - cot_theta_dh_C2 + 4.0*h**2/(psi*C2)*(dpsi_dr - dpsi_dtheta*dh/h**2) + 3.0*dh**2/h
return dHdtheta
Note that the horizon_function now takes the separation as an additional argument.
In [14]:
def shooting_function_binary_general(h0, z):
"""Encode the function phi in order to find the initial radius."""
theta = [0.0, np.pi / 2.0]
H0 = [h0, 0.0]
H = odeint(horizon_function_binary_general, H0, theta, atol = 1.0e-6, rtol = 1.0e-4, args = (z,))
# Return dh/dtheta at pi / 2.
return H[-1, 1]
Note again that we need an additional argument corresponding to the separation, and this is passed through to the horizon_function using the args parameter.
In [15]:
h0_range_1 = np.linspace(1.2, 1.4)
phi_range_1 = np.zeros_like(h0_range_1)
h0_range_2 = np.linspace(1.225, 1.4)
phi_range_2 = np.zeros_like(h0_range_2)
for i in range(len(h0_range_1)):
phi_range_1[i] = shooting_function_binary_general(h0_range_1[i], 0.75)
phi_range_2[i] = shooting_function_binary_general(h0_range_2[i], 0.775)
fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(111)
ax.plot(h0_range_1, phi_range_1, h0_range_2, phi_range_2)
ax.set_xlabel("$h_0$")
ax.set_ylabel("$\phi$")
ax.axhline(color = 'black', linestyle = '--')
plt.show()
We see that $0.75 < z^* < 0.775$. Now all we need to do is find the minimum of this function to define our new function $\Phi$, then find the root of that. Finding the minimum can be done by another standard scipy function.
In [16]:
from scipy.optimize import minimize_scalar
In [17]:
def critical_separation_function(z):
"""Function whose root gives the critical separation."""
result = minimize_scalar(shooting_function_binary_general, [1.23, 1.3], args=(z,))
return result.fun
Now we can find the root of this function, telling us the critical mass.
In [18]:
zstar = brentq(critical_separation_function, 0.75, 0.775)
In [19]:
print("The critical value of the separation is z* = {}".format(zstar))
Now we solve for the horizon at the critical separation. It is pretty much impossible to find the root using the brentq algorithm as the shooting function just touches zero (by design!). However, we know that the root occurs at the minimum location.
In [20]:
result_critical = minimize_scalar(shooting_function_binary_general, [1.23, 1.3], args=(zstar,))
h0_critical = result_critical.x
print(r"'Found' horizon with radius at \theta=0 of {}".format(h0_critical))
In [21]:
theta = np.linspace(0.0, np.pi / 2.0)
H0 = [h0_critical, 0.0]
H = odeint(horizon_function_binary_general, H0, theta, atol = 1.0e-6, rtol = 1.0e-4, args=(zstar,))
h = H[:, 0]
x_critical = h * np.sin(theta)
z_critical = h * np.cos(theta)
fig = plt.figure(figsize = (8, 8))
ax = fig.add_subplot(111)
ax.plot(x_critical, z_critical)
ax.set_xlabel("$x$")
ax.set_ylabel("$z$")
ax.set_title("Horizon location for the 'last' symmetric binary, singularities at $z = \pm {}$.".format(zstar))
ax.set_xbound(0.0, 1.4)
ax.set_ybound(0.0, 1.4)
plt.show()