Black hole formation from the collision of plane-fronted gravitational waves

This paper introduces a new effort to study the collision of plane-fronted gravitational waves in four-dimensional, asymptotically flat spacetime, using numerical solutions of the Einstein equations. The pure vacuum problem requires singular, Aichelburg-Sexl-type sources to achieve finite energy solutions, which are problematic to treat both mathematically and numerically. Instead then, we use null (massless) particles to source nontrivial geometry within the initial wave fronts. The main purposes of this paper are to (a) motivate the problem, (b) introduce methods for numerically solving the Einstein equations coupled to distributions of collisionless massless or massive particles, and (c) present a first result on the formation of black holes in the head-on collision of axisymmetric distributions of null particles. Regarding the last-named, initial conditions are chosen so that a black hole forms promptly, with essentially no matter escaping the collision. This can be interpreted as approaching the ultrarelativistic collision problem from within an infinite boost limit, but where the matter distribution is spread out, and thus nonsingular. We find results that are consistent with earlier perturbative calculations of the collision of Aichelburg-Sexl singularities, as well as numerical studies of the high-speed collision of boson stars, black holes, and fluid stars: a black hole is formed containing most of the energy of the spacetime, with the remaining $15\pm 1\%$ of the initial energy radiated away as gravitational waves. The methods developed here could be relevant for other problems in strong-field gravity and cosmology that involve particle distributions of matter.

Figure 1

A spacetime diagram depicting a null-fronted plane wave (shaded region) traveling in the $+x$ direction. $u=t-x$ and $v=t+x$ are null coordinates; the two coordinates $y$, $z$ transverse to the wave are not shown. On either side of the pulse $u<0$ and $u>\mathrm{\Delta }u$, the geometry is that of Minkowski spacetime.

Figure 2

A spacetime diagram depicting the collision of left and right moving null-fronted plane waves. The nonlinear interaction occurs in the region $(u>0,v>0)$, depicted by the darker shaded region. To the past of this, the spacetime geometry is either that of one of the null waves (lighter shaded regions), or Minkowski spacetime.

Figure 3

Log of the L2 norm of the constraints (27) over the inner $x\in [-125M,125M]$, $y\in [0,125M]$ portion of the computational domain (units/normalization arbitrary up to a constant shift). Shown are three characteristic resolutions, with the number of particles $n_0$ scaled as explained in the text when changing the base resolution $h_0$. For the $h_0$ case, data from two additional runs are shown with different numbers of particles, demonstrating that here we are essentially in the domain where the grid-based truncation error is dominating the solution error. The trends with resolution are broadly consistent with second order convergence; some of the spikes in the higher resolution curves, particularly noticeable for the highest resolution (solid black) curve at early times, are due to the mesh-refinement algorithm temporarily dropping a highest resolution grid, and so during that time in the region about the black hole (which dominates the constraint error) the grid resolutions of the $h_0$ and $h_0/2$ simulations are actually the same. Of course, we could have required the hierarchies to be identical amongst the runs to give cleaner looking convergence plots, but as that is not typically how we do runs, the above is a more representative example of the convergence behavior of the code.

Figure 4

Log of the L2 norm over all geodesics of the normalization constraint ${\ell }_a{\ell }^a=0$ when performing a free evolution of the geodesic equation (units/normalization arbitrary up to a constant shift). The curves end once an apparent horizon is first detected, after which we remove any geodesics (the vast majority of them for this initial data) in the then excised region of the domain. Constrained evolution gives lower norms for the corresponding error calculated when adjusting the ${\ell }^t$ component of each geodesic to enforce the constraint, though for these short geodesic evolutions the difference in the affect on the spacetime truncation error is negligible.

Figure 5

Four snapshots of ${}^{(r)}{\mathrm{\Psi }}_4\times rM$ from the collision described in the text. Times from top left to bottom right are $t=0M$, $5M$, $20M$, and $50M$. The symmetry ($x$) axis is horizontal and passes through the center of each panel; we only simulate the top ($y\ge 0$) half of the plane, but for visualization purposes also show $y\le 0$. The size of each panel is $150M\times 150M$. The central black dots in the panels after $t=0$ indicate the black hole that formed, and are specifically the excised regions, set to 60% of the size of the apparent horizon at each time. By $t=5M$, $>99.9\%$ of the null particles sourcing the initial wave fronts have fallen into the black hole, and to excellent approximation then and in subsequent panels the spacetime is vacuum exterior to the excised regions. (Note that the implied amplitude of the initial wave fronts in the top left panel is misleading. At $t=0$, the two pulses are moving in the $x$ direction, but are offset from $y=0$, so that ${}^{(4)}{\mathrm{\Psi }}_4$ measures some radiation is simply because the $r$ tetrad direction has some overlap with the propagation direction. However, that the magnitude appears not to decay with distance from origin at $t=0$ is entirely due to resolution effects [${}^{(x)}{\mathrm{\Psi }}_4$ given by (21) drops off like $1/{\rho }^2$]: the mesh-refinement algorithm successively lowers the resolution moving outward, and because of how thin in $x$ the initial waves are, far from the origin these features are significantly under-resolved. This is exacerbated in the calculation of ${}^{(r)}{\mathrm{\Psi }}_4$, which requires second gradients of the metric. As evolution proceeds, the Kreiss-Oliger numerical dissipation we use smooths out these under-resolved features.)

Figure 6

${}^{(r)}{\mathrm{\Psi }}_4(t)\times rM$ measured at three radii and two angles relative to the origin (the top panel is at a right angle relative to the collision axis, the bottom panel is at $\varphi =\pi /32$ from the collision axis); see Fig. 7 for the same data on a logarithmic vertical scale. The $r=75M$ and $r=100M$ date were shifted in time by the amounts shown in the legend to account for the different propagation times (and note that these shifts are different for the two panels). Though the simulation was run to $t=175M$, the data for the $r=100M$ point was truncated at $t=150M$ to avoid possible artifacts coming from the computational boundaries at $|x|=y=250M$, assuming a unit coordinate light speed. Note that, particularly near the axis, we are not yet sufficiently far from the collision that we see the expected $1/r$ decay of the wave. (In the wave zone, given the scaling of ${\mathrm{\Psi }}_4$ by $rM$, the shifted waves should have the same amplitudes.)

Figure 7

The natural logarithm of the magnitude of the same data shown in Fig. 6. For reference, overlayed (orange double-dot dash) is a straight line segment with slope $-0.089$, which is the expected QNM decay rate of the least damped $\ell =2$ quadrupolar mode (see, e.g., [80]); note that this line is not a fit to the data.

Figure 8

The mass of the black hole estimated from the area of the apparent horizon, for three different resolutions. The apparent horizon is first detected around $t=1.1M$ in all cases, though between $t\sim 4–7M$ the apparent horizon finder fails to find it to within a reasonable tolerance (we use a flow method, assuming a “star shaped” surface, and though the horizon is quite deformed here, it is not close to violating this assumption). We are not sure why this happens, however, the temporary loss of this surface that guides excision (the excision surface is set to the same shape as the apparent horizon, but 60% its size, and is frozen in shape when the finder fails) does not affect stability at the excision surface, implying characteristics remain ingoing there.

Figure 9

The ratio of proper equatorial to proper polar circumference $C_{eq}/C_p$ of the apparent horizon from the $h_0$ resolution case, to illustrate the initial dynamics of the horizon and its ringdown to a Schwarzschild black hole. For the default resolution we use to resolve the horizon, namely uniformly discretized in $\varphi$ with $N_{ah}=33$ points, $C_{eq}/C_p$ is not calculated with a similar accuracy as the area; see Fig. 8 (for which the estimated fractional error in $M_{ah}$ is 0.8%, compared to 3.3% for $C_{eq}/C_p$). However, this is not a reflection of the underlying accuracy of the spacetime solution, as demonstrated by the higher resolution apparent horizon finder curves, with the same $h_0$ spacetime resolution. $C_{eq}/C_p$ is converging to the expected value of one at first order [dominated by integration truncation error on the axis ($\varphi =0$, $\pi$)]. We do not typically run with higher apparent horizon resolutions, as the flow method’s computation time scales poorly with $N_{ah}$, and is also more challenging to robustly find a very dynamical horizon to within low tolerance, as indicated by additional times between $t\sim 18–22M$ when the higher resolution cases fail to find the horizon.

Figure 10

The volume average of the L2 norm of the constraints (27) over the domain at several different spatial resolutions and particle numbers. At these resolutions, the constraint violation is dominated by the spatial discretization truncation error, as opposed to that coming from the finite number of particles (the $N^{1/3}=512$ case is not continued as long as the others for this reason). The decrease in the constraint violation with increasing spatial resolution is faster than second order convergence, suggesting that the error from reconstructing the stress-energy tensor is subdominant to that coming from evolving the metric and particles.

Figure 11

Comparison of the density contrast measured by an observer comoving with the matter at the points of maximum (labeled OD) and minimum (labeled UD) density vs a local measure of the scale factor (see [22] for exact definition). The results obtained with particle evolution closely track those obtained in [22] by evolving the fluid equations.