Black Hole Formation from the Collision of Plane-Fronted Gravitational Waves

# Black Hole Formation from the Collision of Plane-Fronted Gravitational Waves

Frans Pretorius Department of Physics, Princeton University, Princeton, New Jersey 08544, USA. CIFAR, Cosmology & Gravity Program, Toronto, Ontario M5G 1Z8, Canada    William E. East Perimeter Institute for Theoretical Physics, Waterloo, Ontario N2L 2Y5, Canada
###### Abstract

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 non-trivial 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 ultra-relativistic collision problem from within an infinite boost limit, but where the matter distribution is spread out, and thus non-singular. 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 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.

## I Introduction

As a theory of gravity, one aspect of general relativity that bares stark contrast to Newton’s theory is the non-linear nature of the Einstein field equations. This imbues general relativity with a strong-field regime, loosely defined as the class of solutions where the non-linearities of the theory are manifest, exhibiting qualitatively different properties and phenomena than weak-field solutions. These include black holes, cosmological solutions describing the global structure and evolution of the universe, and gravitational collapse. The last-named refers to situations where, beginning from non-singular initial data (even weak-field data where linearity is satisfied to good approximation early on), evolution eventually leads to the formation of some kind of singularity in the structure of spacetime.

Barring questions about the very early universe (the pre-inflationary epoch, or bounce in cyclic models), and whether dark energy or dark matter might be due to a failure of general relativity describing gravity on these scales, knowledge of strong-field gravity relevant to the observable universe has become quite thorough over the past several decades. This includes the geometries of neutron stars and black hole exteriors, the spacetime dynamics in mergers of such compact objects, and knowledge of how the exterior geometry settles to a stationary Kerr solution when gravitational collapse occurs (see, e.g., Shibata and Taniguchi (2011); Fryer and New (2011); Faber and Rasio (2012); Blanchet (2014); Lehner and Pretorius (2014); Paschalidis and Stergioulas (2017) for some review articles). Of course, many details of these processes remain to be understood, and fundamental open questions remain where matter plays an important or dominant role—core-collapse supernova, self-consistent models connecting binary neutron star mergers to the host of observed/hypothesized electromagnetic counter parts, how accreting black holes power observed jets, etc.

However, many other aspects of strong-field general relativity (of arguably less astrophysical importance, but still of significant physical and mathematical interest) remain poorly understood, despite the relevant open questions having been identified decades ago. These include the interior structure of black holes and the generic nature of singularities formed in gravitational collapse, critical phenomena at the threshold of black hole formation (in particular in situations lacking symmetries), the classification and properties of solutions with event horizons in higher dimensional spacetimes, gravity in asymptotically Anti-de Sitter (AdS) spacetimes, and the ultra-relativistic scattering problem (see e.g. Gundlach (2003); Berger (2002); Emparan and Reall (2008); Dafermos and Rodnianski (2013); Chesler and Yaffe (2013); Choptuik et al. (2015); Garfinkle (2017); Dafermos and Luk (2017); Barack et al. (2018); Coley (2018)). The goal of this paper is to introduce a new research effort to tackle aspects of this last-named problem.

As part of this effort, we develop methods for evolving distributions of collisionless particles coupled to Einstein gravity. Though the main focus of this work is the null particle case, these same methods can be used for massive particles, as we demonstrate in passing in this work. This could be relevant for tackling a number of other problems in strong field gravity, such as studies of cosmic censorship Shapiro and Teukolsky (1991) or inhomogeneous cosmologies. Recently, there has been interest in performing fully general-relativistic calculations of structure formation scenarios Bentivegna and Bruni (2016); Giblin et al. (2016a, b); Macpherson et al. (2017); East et al. (2018); Macpherson et al. (2018) to study effects not captured by standard Newtonian N-body calculations. However, current approaches rely on fluid descriptions that breakdown in the presence of multistream regions that arise, e.g., during halo formation, or use weak-field approximations to the Einstein equations Adamek et al. (2016).

Before outlining the main content of the remainder of the paper, we describe the ultra-relativistic scattering problem in more detail, giving some background, and a list of open questions we would eventually like to address.

### i.1 The Ultra-relativistic scattering problem

In the context of strong-field gravity, the ultra-relativistic scattering problem refers to understanding the dynamics of spacetime following the interaction of two distributions of energy, initially approaching each other from opposite directions near or at the speed of light as measured in the center-of-momentum frame. The source of energy could be pure spacetime itself, i.e. plane-fronted gravitational waves or black holes (including Aichelburg-Sexl singularities in the infinite boost limit), or some classical model of a particle, such as boson stars, fluid stars, or even black holes. Two seminal works in this area where initiated by Khan and Penrose Khan and Penrose (1971), and Penrose Penrose (1974) almost 50 years ago, essentially addressing two opposite extremes in the landscape of the scattering problem : collision of waves of infinite transverse extent, and the collision of point-particles.

#### i.1.1 Collision of planar gravitational waves

Khan and Penrose Khan and Penrose (1971) studied the collision of two infinite, plane-fronted gravitational waves, and found the interaction always resulted in a naked curvature singularity, regardless of how weak (in terms of curvature) each individual wave is. The formation of the singularity can be understood as the result of the focusing of the geometry of one wave front by the other, and vice versa. The weaker the initial curvature is, the longer it takes to focus to a singularity, though it always does. One can argue that formation of a singularity in this scenario is an artifact of both the perfect focusing caused by the planar geometry of the wavefront, as well as its infinite transverse extent (hence the spacetime has infinite energy and is not asymptotically flat). Moreover, even in the absence of any dynamics resulting from the collision of two such wave fronts, a single wave, that is locally exactly Minkowski spacetime on either side of the front, is not globally hyperbolic: any spacelike hypersurface of the single wave geometry evolved forward in time using the Einstein equations will encounter a Cauchy horizon Penrose (1965). In a sense, uniform plane-fronted gravitational waves are infinitely powerful lenses, capable of focusing all of Minkowski spacetime down to a point. We illustrate this more explicitly in Sec. II.2.

The open questions pertaining to this limit of the ultra-relativistic scattering problem then relate to how these “pathological” outcomes might change if the apparent sources of the pathology —a uniform distribution and infinite total energy— are removed. If the energy density in each wave has finite transverse extent, one would expect the focusing to end within a transverse light-crossing time. If this is longer than the time for the singularity to form, will a black hole eventually form to censor the singularity? How do inhomogeneities in the energy density affect the evolution of the wave-front post collision?111These questions could be of relevance to certain bubble formation/collision scenarios in the early universe Hawking et al. (1982). A bubble of true vacuum nucleated in a false vacuum will expand, with the bubble wall gaining all the energy of the false vacuum swept up. The wall quickly becomes relativistic, and eventually strongly self-gravitating. The wall will therefore begin acting as a strong lens focusing matter it encounters, and conversely inhomogeneity in the matter will back-react to create inhomogeneity in the bubble wall. If two such bubbles collide, in a region about the initial point of contact much smaller than the radius of curvature of each bubble, to good approximation the collision could be treated as the collision of two plane-fronted waves. The curvature in the walls around the collision is not expected to be strong enough to focus to singularities before the non-planar nature of the walls influence the dynamics (if that were the case, the bubbles would have individually collapsed to black holes prior to the collision(see e.g. Deng et al. (2018))). Though self-gravitating bubble collisions have been studied extensively before (see e.g. Wainwright et al. (2014a, b); Johnson et al. (2016) and the references therein), these studies have not included the effect of inhomogeneities, which would be interesting to explore. See Braden et al. (2015a, b); Bond et al. (2015) for related work suggesting domain walls (neglecting self-gravity) are unstable to perturbations.

#### i.1.2 Collision of point-particles

At the other end of the spectrum, Penrose Penrose (1974) initiated the study of the collision of plane-fronted waves sourced by singular, point-like distributions of energy. Such a plane wave geometry can be obtained by taking the Schwarzschild solution with rest mass , boosting it with Lorentz factor , and considering the limit that and while the energy remains fixed. This yields the Aichelburg-Sexl (AS) solution Aichelburg and Sexl (1971). The AS geometry is Minkowski on either side of the shock front, with all the curvature confined to the transverse plane. The magnitude of the gravitational wave as measured by a Newman-Penrose scalar drops like with transverse distance from the origin, and though not asymptotically flat in the strictest sense of the definition D’eath and Payne (1992a), the spacetime still approaches Minkowski at large , and does not suffer the infinite energy/focusing pathologies of the previously discussed uniform plane wave solutions. This however comes at the expense of the origin now being a naked, curvature singularity. Nevertheless, Penrose was able to show that an apparent horizon exists in the spacetime formed by the superposition of two colliding AS shocks; the exact solution to the causal future of the collision is unknown, but this implies a black hole forms. Though an extreme example of relaxing the infinite extent problems of uniform plane wave collisions, this is suggestive that the pathologies with the latter scenario are indeed due to infinite extent, and not properties of the non-linear interaction in gravitational wave scattering.

The study of the point-particle limit was reinvigorated a couple of decades ago following the realization that if extra dimensions exist, the true Planck scale could be orders of magnitude lower than the scale naively inferred from 4D dimensional analysis Arkani-Hamed et al. (1998); Antoniadis et al. (1998); Randall and Sundrum (1999). It was argued one “natural” magnitude for the Planck energy is . This renewed interest in the gravitational scattering problem because particle collisions above the Planck energy are generically expected to form black holes ’t Hooft (1987); Amati et al. (1987); Banks and Fischler (1999), so if the low-TeV Planck-scale scenario is true, the Large Hadron Collider, as well as high energy cosmic ray collisions with the Earth’s atmosphere, could then form black holes Dimopoulos and Landsberg (2001); Giddings and Thomas (2002); Feng and Shapere (2002). No evidence for this has been found to date; see e.g. CMS Collaboration (2018).

At sufficiently high energies in particle collisions, the black holes that form would be large enough to censor any details of the collisions, and hence it is conjectured that in the ultra-relativistic limit “matter does not matter”; i.e., gravity dominates the interaction, and moreover, the geometry of each boosted particle is dominated by its kinetic energy, hence the ultrarelativistic limit is uniquely captured by the collision of two AS shock waves. This argument certainly makes intuitive sense, though from a geometric perspective, given the rather stark differences between the geometries of large-but-finite boost time-like compact objects and the singular, plane-fronted AS shock-wave solution222For example, all polynomial invariant scalars of the Riemann tensor, such as the Kretschmann scalar, vanish for null, plane-fronted gravitational wave spacetimes Kramer et al. (1980), including the AS spacetime., it would be rather remarkable if this conjecture were true. Nevertheless, evidence has been gathered in its favor from simulations of ultra-relativistic boson star Choptuik and Pretorius (2010), fluid star East and Pretorius (2013), and black hole collisions Sperhake et al. (2008): for collisions between material objects, black holes do form above a threshold roughly in line with hoop-conjecture arguments Thorne (1972), and in all cases the gravitational wave energy emission in head-on collisions, extrapolated to infinite boosts, agrees with perturbative calculations of that produced in the collision of two AS waves D’eath and Payne (1992b) (the last-named finding of the initial mass of the spacetime being radiated).

Many open questions remain here. Regarding the apparent limiting behavior of large boosts, one challenge to establish the connection with AS is that detailed comparisons will require explicit solutions, and it is unclear how to deal with the AS naked singularity, in particular in a numerical solution scheme. The approach we propose here is to begin at the infinite boost limit, i.e. with null plane-fronted waves, and replace the singularity with a smooth null matter source (such solutions are sometimes called gyratons, the first examples of which were discovered by Bonnor Bonnor (1969) and Peres Peres (1959, 1960)). This would also allow one to address how “stable” the AS singularity is in the first place, by studying the stability of any family of regular spacetimes that approach the AS spacetime in some continuous limit. One conceivable outcome of instability is that perturbations generically lead to black hole formation.

Not much is known about the dynamics of off-axis collisions. Perturbative calculations of large impact parameter scattering suggests the ultra-relativistic two body problem is significantly simpler than the more astrophysical, rest-mass dominated regime Galley and Porto (2013); Gruzinov and Veneziano (2016). This is consistent with studies that indicate some of the leading order physics in these interactions, even black hole formation, can be captured by appealing to geodesic motion on a single, relevant background geometry Kaloper and Terning (2008); East and Pretorius (2013). Another result from perturbative calculations of large impact parameter deflections is the radiation produced is highly beamed Galley and Porto (2013). This, together with geodesic focusing of the emitted radiation, could explain the rather striking growth in black hole mass noted during a moderate-boost () grazing encounter simulation of two black holes Sperhake et al. (2013). It would be interesting to explore these grazing encounters at much higher boosts.

A further set of questions relate to the threshold of black hole formation. In particular, if critical phenomena Choptuik (1993) is present, and if so, when tuning to threshold, which critical solution is revealed: that of the underlying matter source, or of vacuum gravity.

Of course, not all questions pertaining to the ultra-relativistic scattering problem need to try to make a connection with one of the two extreme limits : point-particle scattering or infinite uniform plane-wave scattering. There is potentially a vast landscape of interesting phenomenology in between, worthy of study in its own right.

### i.2 Outline of the remainder of the paper

In Sec. II we describe the formalism in more detail: the Einstein equations coupled to null matter, certain properties of plane-fronted wave solutions, and the similarities/differences between the pure vacuum versus matter sourced cases. In Sec. III we describe a new code designed to solve the Einstein-collisionless particle system, in particular highlighting the issues required to self-consistently and efficiently compute the stress energy tensor summed over the contributions from the discrete collection of particles. In Sec. IV we present results from simulations of the head-on collision of two axisymmetric null-particle distributions, including convergence tests. The main result is that we find a spacetime consistent with prior approaches to studying the AS collision limit: a black hole forms containing the initial energy of the spacetime, with the rest escaping as gravitational radiation (modulo in energy of the tail end of the null particle distribution that did not fall into the black hole). Regarding the structure of the waves, we find it is highly beamed about the collision axis. Also, the plane-fronted shock-like features of the initial geometries are trapped by the black hole, continually propagate about it, albeit with an amplitude that exponentially decreases with time. We conclude in Sec. V by mentioning some improvements to the code that would be needed before all the topics discussed above can be addressed. In the appendix we further validate the particle code, showing how it can be applied to the massive particle case, by applying it to a simple inhomogeneous cosmology setup.

We use geometric units where Newton’s constant and the speed of light are set to unity.

## Ii Formalism

In this section we discussion the Einstein equations coupled to a null fluid (Sec. II.1); the problems with homogeneous, plane-fronted gravitational wave spacetimes (Sec. II.2); how these problems can be resolved when using a null fluid as the source of the plane-wave geometry (Sec. II.3); and the issues of using a null fluid to study colliding, regular plane-wave geometries and how null particle distributions can alleviate them (Sec. II.4).

### ii.1 The Einstein field equations with a null fluid source

We first consider the Einstein equations with a pressureless fluid as a source:

 Rab−12Rgab = 8πTab ≡ 8πρeℓaℓb,

where is the Ricci tensor, the Ricci scalar, the metric tensor, the stress energy tensor, and is the energy density flowing along the direction . For a timelike pressureless fluid is the four velocity of the fluid (with ) and has an unambiguous interpretation as the rest-frame density of the fluid. However, for a null fluid (), as we consider below, is the energy density in some chosen frame. This frame then determines the normalization of , which is equivalent to fixing the affine parameter in the parametric representation of the null vector, (or vice versa: a frame where the normalization of has been chosen can be thought of as the frame defining the interpretation of ).

Later we will generalize this to a distribution of non-interacting null particles, though for now the null fluid is more convenient to illustrate plane-fronted wave solutions sourced by matter. Moreover, prior to the collision of two such fronts, there is a simple one-to-one correspondence between the fluid and particle solutions, and the former is more convenient to use to provide initial data for the collision.

### ii.2 Plane-fronted waves

Consider a plane-fronted wave traveling in the direction—see Fig. 1. The metric, in Brinkmann-like form, is

 ds2=−dudv+B2(u)[dy2+dz2]+H(u,y,z)du2, (2)

where is a null coordinate, is null when , and are the Cartesian-like coordinates transverse to the wave. There are several different coordinate systems often used to represent plane-fronted metrics (see Kramer et al. (1980)); the above is convenient for our discussion related to null-sourced waves, and to motivate the eventual form of the metric we use to set initial conditions for numerical evolution (in Sec. IV.1 below). The only non-trivial Einstein equation (II.1) for this metric ansatz is

 ∇2HB2+4¨BB+8πρe=0, (3)

where is the 2-dimensional flatspace Laplacian, and the overdot () denotes . With respect to orthonormal basis vectors aligned with the coordinate directions, the only non-zero Newman-Penrose scalar for this metric is , defined to measure gravitational waves propagating in the direction333Hence the -superscript, to differentiate it from shown in the results section, which is calculated with a tetrad to measure radiation propagating in the radial direction.:

 (x)Ψ4=H,zz−H,yy2B2−iH,yzB2, (4)

where here and below we use a comma to denote the ordinary (partial) derivative operator.

To study pure gravitational wave spacetimes () it is most convenient to set and ; then the field equations allow an arbitrary amplitude profile , with the transverse structure constrained by , and it is manifest that the spacetime is Minkowski away from the wavefront (when ). This form of the metric is adequate to capture the entire class of vacuum, plane-fronted gravitational wave solutions Kramer et al. (1980), and also highlights why these spacetimes are problematic to address some of the questions identified in the introduction: uniqueness properties of dictate the only way to obtain finite transverse extent, or inhomogeneity in the transverse plane, is via AS-like singularities, or boundary conditions at infinity. The latter is essentially how inhomogeneities are introduced in asymptotically AdS spacetimes, were gravitational wave scattering is used to model heavy ion collisions using gauge-gravity dualities (see e.g. Chesler et al. (2015)), though this option is not available in asymptotically flat spacetime. Before describing how null dust can circumvent these problems, it is useful to consider the homogeneous dust-sourced spacetimes, with .

The simplest homogeneous dust solutions are those without any “pure” gravitational waves (i.e. all non-trivial curvature is in the trace-full part of the Riemann tensor constrained by the Einstein equations, with the Weyl tensor identically zero) : , with satisfying

 ¨B+2πρeB=0. (5)

Suppose we have a matter pulse with support in (Fig. 1), and we choose initial conditions to (5) such that ; i.e., the metric [(2) with ] is in standard Minkowski form prior to passage of the wave. If is positive, will begin to decrease as increases, and always reach zero, either within the pulse for sufficiently high amplitude and/or wide pulses, or sometime after the passage of the pulse. For example, for a constant density pulse with within , otherwise, the solution to (5) is

 B(u) = 1,   u<0 (6) = cos(ω0u),   0≤u≤Δu = cos(ω0Δu) −(u−Δu)(ω0sin(ω0Δu)),   u>Δu,

with . For the coordinate singularity to occur behind the pulse, we must have ). To see that this is a coordinate singularity, consider the following coordinate transformation ( is unchanged):

 y = ¯y/B(u), z = ¯z/B(u), v = ¯v+L(u)[¯y2+¯z2] (7)

If we choose , then away from the matter where , the above transforms (2) with to .

On the other hand, is more than a coordinate singularity. Consider the focusing caused by the passage of the wave on a set of timelike geodesics, initially at rest. In particular, looking at the transverse coordinates of any such geodesic, if , , then the geodesic equation says , ; i.e. in these coordinates such geodesics remain at fixed transverse coordinate locations. However, from (2) the proper transverse (geodesic) distance between any pair of geodesics with coordinate separation () at some is . In other words, the proper distance between all these timelike geodesics goes to zero at , when . In that sense the plane wavefront is an infinitely powerful lens. Moreover, given that curves are geodesics, and and are Killing vectors of the spacetime, we can toroidally compactify the space by identifying with , and with for some constants . The geodesic focusing then tells us that the entire compactified transverse space is focused to a point at . This illustrates that it is more appropriate to think of the plane wavefront as focusing all of spacetime, and not merely a class of geodesics. Similarly, considering the regular coordinate chart, except for where (), the entire range is mapped to when . This demonstrates that these regular coordinates cannot be used to specify complete, Cauchy data for the region .

### ii.3 Regular, plane-fronted waves

As mentioned above, to obtain vacuum plane-fronted gravitational waves that have finite total energy and are asymptotically flat transverse to the wavefront requires singular sources. Using a null fluid instead can remedy the problem. The form of the metric ansatz (2) used above is not ideal for numerical evolution, in that for a localized source asymptotically , with . Moreover, the focusing can still lead to coordinate singularities near the wake of the fluid. A better suited coordinate system can be found appealing to the coordinate transformation (II.2), which essentially “unfocuses” the metric following the wave (however unlike (II.2) for a homogeneous wave, the finite wave has no Cauchy horizon problems). For the axisymmetric collisions explored later, it will also be more convenient to use cylindrical coordinates. We therefore use the following ansatz for the metric and fluid source propagating along characteristics, with the specific form of the functions chosen to simplify the resulting field equations:

 ds2 = −dudv−8πf(u)h(ρ,θ)du2 (8) +2β(u)q(ρ)dudρ+dρ2+ρ2dθ2, Tab = ρe(u,ρ,θ)ℓaℓb, (9)

with

 β(u) = 4π∫u−∞f(~u)d~u, (10) ρe(u,ρ,θ) = f(u)g(ρ,θ), (11)

and null vector normalized to . With this ansatz, the one non-trivial Einstein equation is:

 ∇2h≡h,ρρ+h,ρρ+h,θθρ2=g−q,ρ−qρ (12)

If we further decompose the transverse dependence of the energy density into cylindrical harmonics

 g(ρ,θ)≡g0(ρ)+∞∑m=1gm(ρ)cos(mθ), (13)

then we can use to solve for the monopole contribution to (12):

 q′+qρ=g0, (14)

where denotes differentiation with respect to . then captures the metric response to a non-axisymmetric source, which can readily be solved via a similar decomposition

 h(ρ,θ)≡∞∑m=1hm(ρ)cos(mθ), (15)

where for each the remaining portion of the field equation (12) reduces to the following ordinary differential equation (ODE)

 h′′m+h′mρ−m2hmρ2=gm. (16)

When solving the above equations, we are free to choose the energy density via the functions and , and then solve for the remaining metric functions using (10),(14),(15) and (16). We will restrict the class of free initial data to that which is regular in the limit , as well as an asymptotically flat spacetime in the limit . In the limit , regular solutions require

 g0(ρ) = α0+O(ρ2) (17) gm(ρ) = γ0ρm+O(ρm+2), (18)

where and are constants. In the limit we require sufficiently rapidly to give the spacetime finite total mass. For the energy density profiles, we either use compactly supported pulses ( for ), or a Gaussian (), the latter which we use in the numerical evolution. Then, the metric (8) asymptotes to Minkowski as , (though we have for the axisymmetric numerical solutions discussed later).

We will use the ADM (Arnowitt-Deser-Misner) mass Arnowitt et al. (1962) as a measure of the spacetime energy, integrated on a cylinder at (arbitrarily) , with , centered about the pulse, and taking the size of the cylinder to :

 MADM = β(∞)8[∫∞0(dq(~ρ)d~ρ~ρ+q(~ρ)) d~ρ (19) +limρ→∞ρq(ρ)].

In reaching the above form for , we have assumed the pulse has finite extent in , and that at least as fast as ; thus the asymmetric piece of the metric does not contribute to the ADM mass. We can simplify the expression using (14), giving

We are not aware of studies that have explored the validity of the ADM mass measure for this class of spacetime in these coordinates. In Sec. IV.1 we show that with the family of initial data we use, taking the AS limit with fixed as defined above does give the AS solution with mass parameter equal to . Interestingly though, if we directly evaluate (20) with the exact AS solution, we obtain . For the regular, limiting sequence of solutions there is equal contribution of from the integral along the end-cap of the cylinder (first term in (20)) and the integral over the barrel of the cylinder (second term in (20)). However, the former piece identically vanishes in vacuum. The AS solution is singular on the axis (when ) in these coordinates, which may ultimately be the source of the discrepancy.

The equivalent expression for (4) is more complicated in the coordinates (8). We will not reproduce the full expression here, though it is illuminating to consider the simpler case of an axisymmetric wave, namely (12) and (13):

 (x)Ψ4,m=0 (21) =2πf(u) (2q(ρ)ρ−g0(ρ))[cos(2θ)+isin(2θ)]

where we have substituted in (14). Note that the spatial tetrad with respect to which we have defined is still aligned with the coordinates, as with (4), hence the dependence in (21) is due to the rotation of this tetrad relative to the coordinate basis (i.e. the metric distortion induced by this gravitational wave is axisymmetric by construction, and would be manifestly so as measured by built out of a tetrad aligned with the coordinates). For an infinite, homogeneous matter source with a constant , and, as before, vanishes and the spacetime has no Weyl curvature. However, for a source with finite transverse extent, outside the region of matter where , and so . Thus, a localized, plane-fronted null fluid wave does act as a source of “pure” plane-fronted gravitational waves, it is just that the Weyl curvature happens to vanish within an inner core about the axis if that core has constant matter density, and the matter distribution remains axisymmetric throughout the spacetime.

### ii.4 Fluid to particle distributions

Though a null fluid is useful to understand plane-fronted gravitational wave spacetimes propagating in a single direction, this matter model is not ideal to explore collisions of such waves. The reason is such matter easily forms caustics, where the assumption that a single, unique velocity vector field describing the fluid flow exists, breaks down. Correspondingly, the Euler equations governing the flow break down at the caustic, and the solution cannot be uniquely extended beyond the caustic. This problem of lack of a unique velocity vector is actually more pronounced than merely being associated with caustics, as can be seen imagining the case where we collide two identical distributions head on, as follows. A pressureless fluid does not self-interact, and in the limit of weak gravity where the fluid is propagating in flat space, we thus expect the two opposing streams to simply pass through each other. This implies that, in the lab frame where both streams are observed to have identical energy profiles, there will be a moment of time symmetry as they pass through each other where instantaneously the superposed profiles have zero velocity. This is impossible to realize using a single null fluid with stress energy tensor of the form (II.1) (i.e., a null vector is incapable of describing a moment of time symmetry). The problem can be remedied for this particular scenario by considering two independent, non-interacting fluids, one describing the right propagating pulse, the other the left. The stress energy tensor is a sum of the two individual fluid stress energy tensors, and this sum will accurately reflect the moment of time symmetry. In particular there will be no momentum density, but there will be anisotropic pressure (pressure along the collision axis, none transverse to it)444The situation is somewhat different if a finite- collision is considered with time-like fluid stars. If the two stars are modeled with separate ideal fluids that only interact gravitationally, at the moment of time symmetry in a collision the net pressure tensor will also be anisotropic. However, here, a single fluid model can still be adequate to model the collision (as used in East and Pretorius (2013)), since at the moment of time symmetry the time-like vector describing the fluid flow will simply coincide with the lab-frame’s velocity 4-vector. Though in this latter case the pressure tensor must still be isotropic. In the high speed limit there will be a much steeper pressure gradient along the collision axis (due to the length contraction of each star in the lab frame) than transverse to it, which post-collision will cause evolution of the fluid to proceed in a manner at least qualitatively consistent with the case of two non-interacting fluids (though of course the two models will give very different outcomes in the low velocity limit).. However, when including gravity, a two-fluid model would only be a temporary fix, as the gravitational interaction between the streams will cause focusing, eventually leading to caustics in each flow.

We therefore decide to treat the matter as a distribution of collisionless particles instead, i.e. as governed by the collisionless Boltzmann (or Vlasov) equation. We approximate this using a particle treatment, where the energy distribution is given by the sum of a large number of particles, each locally traveling along a geodesic. This easily resolves the uniqueness issues associate with a single, global fluid vector field. However, this introduces a new problem of what exactly we mean by a “particle.” In Newtonian gravity, the treatment of a point-like distribution of energy is straightforward. In contrast, in general relativity, there is no simple analogue: putting all the matter at a point, however small the rest mass, produces a black hole. That suggests one viable particle model is a collection of black holes for timelike particles, and presumably this could be extended to null particles by locally taking AS limits for each black hole. However, a rigorous (without approximation) implementation of this idea shifts all “matter” to “geometry” on the left hand side of the Einstein equations, and certainly is not something practical to implement numerically, even for one AS particle as discussed above, let alone a large enough number to approximate a continuum energy distribution.

The second, more practical option, is to treat each particle as if it were some form of solitonic matter, and then define some kind of averaging operation that consistently adds the contribution of each particle to the discrete representation of the stress energy tensor used in the Einstein field equations. (Or, alternatively, we can think of our particles as a discrete sampling of the underlying continuum distribution described by the Boltzmann equation, and the problem is how to reconstruct the stress energy tensor from this sampling.) The easiest way to do this would be to assume the particle’s characteristic radius is much smaller than any mesh cell we will use in a numerical scheme, and the particle’s energy is sufficiently small that any self-force effects are negligible compared to numerical truncation error. In principle, one could consider larger particles, as in Newtonian smooth particle hydrodynamic codes. However, finite size effects might then need to be considered. Moreover, for particles moving relativistically length contraction must be taken into account, which complicates the averaging operation for particles that span multiple cells.

We finish this section by writing the relevant definitions and equations for the particle model, with the averaging for the stress tensor implied, but not explicitly stated; in Sec. III.5 we describe the particular averaging procedure we use in the code. The stress energy tensor for a collection of particles schematically takes the form

 Tab=N∑i=1ϵiℓaiℓbi, (22)

where labels the particle, traveling along a geodesic curve with the affine parameter along the curve, and the tangent to the curve:

 ℓai≡dxai(λ)dλ,   ℓaiℓia=−m2i, (23)

(no summation over particle label ). Here is a function that is related to the energy density of the particle, the choice of affine parameter, and also the averaging process. For null particles, the rest mass is zero, while for the massive case, the affine parameter is chosen such that the proper time . Each particle follows a geodesic of the spacetime

 ℓbi∇bℓai=0. (24)

In an appropriate limit of an infinite number of particles, this model should reduce to a collisionless Boltzmann model. Of course, directly discretizing the position-momentum phase space and evolving the density function would also be a viable approach to study the ultra-relativistic scattering problem. The computational difficulty to solving the Boltzmann equation is the additional dimensions required to represent the distribution in phase space. On the flip side, the computational shortcoming of a particle model is the slow convergence to the desired continuum limit.

For distributions that can be described by a single-valued velocity vector field, the particle model is equivalent to the fluid model discussed above for single wavefronts. In the code, as described in the following section, we use the fluid model to provide initial data, then evolve using the particle model. We always use a sufficient number of particles in the latter so that the large approximation to the continuum limit is smaller than numerical truncation error from the discretization of spacetime, as judged by convergence of the constraint equations.

## Iii Numerical Code

In this section we describe the numerical code, reviewing the generalized harmonic formalism we use for representing the Einstein equations (Sec. III.1); some aspects of initial data (Sec.III.2) and gauge choice (Sec. III.3), leaving details to Sec. IV; how we integrate geodesics (Sec. III.4); and how we compute the stress energy tensor from a distribution of geodesics (Sec. III.5). We focus on the case of colliding null waves, but besides the discussion of initial data and gauge conditions, the code is generally applicable to both null and time-like (massive) particles, and we demonstrate an application to an inhomogeneous matter-filled universe in the appendix.

### iii.1 Metric evolution

For the metric tensor, we solve the Einstein equations using the generalized harmonic formalism:

 gcdgab,cd + gcd,(agb)c,d+2H(a,b) − 2HdΓdab+2ΓcdbΓdca (25) = −8π(2Tab−gabT),

where is the metric connection, the trace of the stress energy tensor, and are the so-called source functions, defined by

 Ha≡□xa. (26)

Here round brackets denote symmetrization and . During evolution, the source functions are treated as independent functions, and can be thought of as encoding the coordinate freedom of the spacetime. Therefore, additional conditions/evolution equations need to be supplied for them, with the definition (26) then becoming a constraint

 Ca≡Ha−□xa=0 . (27)

In fact, the time derivative of this constraint essentially gives the usual Hamiltonian and momentum constraints (see, for example, Lindblom et al. (2006)). We numerically solve (III.1) with constraint damping terms Gundlach et al. (2005) (together with the gauge evolution equations discussed later) utilizing a order accurate finite difference code with adaptive mesh refinement (AMR), using the methods and techniques described in Pretorius (2005a, b); East et al. (2012a). (Though since we presently only employ a second order accurate calculation of the stress energy tensor from the particle distribution, as described below, the overall accuracy of the code is second order in the continuum limit.) We use Cartesian-like coordinates where , and for the axisymmetric evolutions presented here use the modified Cartoon method Alcubierre et al. (2001) introduced in Pretorius (2005a). The modified Cartoon approach still uses the Cartesian form of the metric, though only evolves a single slice of the spacetime.

One difference for these simulations compared to earlier studies performed with this code, is here we do not spatially compactify the coordinates. The reason is that with our initial data (see Sec. IV.1 for the explicit solutions) the asymptotic form of the metric, though asymptotically Minkowski, is not in the usual trivial Cartesian form of . Compactification transforms the metric into a representation that is singularity at infinity, but beginning from this form of the metric, the singular portion is easy to factor out analytically, with the code then only storing the regular part. Though this would in principle be possible to do with the null wave initial data as well, given the non-trivial, initial-data-dependent structure of the metric at infinity, it would have required significant updates to the code. The disadvantage to not compactifying is then specifying consistent, physically correct outer boundary conditions becomes challenging. We bypass this issue by placing the outer boundaries in the uncompactified code sufficiently far away that they are out of causal contact with the inner region of the domain where we will measure properties of the solution.

### iii.2 Initial data

For initial data, we superpose two null plane-fronted waves, one propagating to the right (the direction), the other to the left (the direction); see Fig. 2. The right (left) moving wave has compact support in () at (i.e. they do not overlap), and each is Minkowski spacetime on either side of the pulse. For the pulse moving to the right we use the form of the metric and stress tensor discussed in Sec. II.3, transforming to Cartesian coordinates using555As described in Sec. IV.1, for the runs presented here we further transform the metric far to the left (right) of the right (left) moving pulse to alleviate some of the resolution issues that might otherwise arise, but that is immaterial to the discussion here.

 u=t−x, v=t+x, y=ρsinθ,z=ρcosθ. (28)

The boundary conditions for the integral defining (10) are chosen so that the metric to the right of the pulse is in standard Minkowski form (as discussed in Sec. II.3, the form of the metric which is on both sides of the pulse has a divergence within the plane of the pulse). For the pulse moving to the left, an analogous solution is used, but with the non-trivial metric and matter functions depending on instead of , and boundary conditions for the analogous flipped so that the metric is to the left of the pulse. Then, consistent initial data for the Cauchy evolution performed by the code is trivial : at the solution for is exactly that of the right moving pulse, and for is exactly that of the left moving pulse.

### iii.3 Gauge conditions

For the source functions that define the gauge in the generalized harmonic formalism, we begin with the gauge of the initial, superposed exact plane-fronted wave solutions. This superposed gauge is not adequate to use after the interaction of the waves, and so, within an inner volume of the domain where the interaction takes place, we smoothly transition to a variant of the damped harmonic gauge used in earlier high speed soliton collision simulations Choptuik and Pretorius (2010); East and Pretorius (2013). The explicit form of the initial gauge source functions and source function evolution equations are given in Sec. IV.2.

### iii.4 Geodesic evolution

We solve the geodesic equation (24) for each particle as follows (dropping the particle label here for clarity). Instead of evolving the position as a function of affine parameter , we directly integrate as a function of coordinate time , using . Then the geodesic equation can be reduced to the following system of seven first order ordinary differential equations (ODEs), where we use the overdot () to denote :

 ˙xj=ℓj/ℓt ˙ℓa=−Γabcℓbℓc/ℓt (29)

where the metric connection is evaluated at the location of the geodesic . In the code, we integrate these equations using a order Runge-Kutta scheme. We have the option to enforce the normalization condition after each time step; there is no unique way to do this, and we have chosen to solve the normalization constraint for . If the constraint is not explicitly enforced during evolution, then we can use it as a diagnostic to check that it does converge to zero at the expected rate; similarly if we do enforce it, we can use the correction induced to after each time step as the diagnostic quantity that should converge to zero. In tests we have conducted, both schemes perform similarly over the relatively short time scale geodesic evolutions presented here (which for most geodesics is much shorter than the net run time as they are removed from the domain when they enter the excised region inside the apparent horizon), though constrained evolution seems to produce more accurate results at fixed resolution for longer time evolutions.

Initial conditions and are chosen so that the initial distribution of particles produce a stress energy tensor (as discussed in the following section) giving a consistent sampling of the desired fluid continuum limit (9).

### iii.5 Calculating the stress energy tensor

There are two issues that we need to address when calculating the effective stress energy tensor used in the Einstein equations coming from the distribution of particles. First, how to add the contribution of a particle to the stress energy tensor of the cell containing it. Second, how to efficiently incorporate this averaging process into the Berger and Oliger (BO) AMR algorithm we use, specifically as it relates to time sub-cycling, which naively would seem to require integrating the same geodesic multiple times on all resolution levels it overlaps (as the BO algorithm does for continuum evolution equations).

#### iii.5.1 Averaging

The averaging procedure is how we convert the stress energy tensor of a single particle into an equivalent cell-based representation such that we will obtain the same solution to the Einstein equations in the continuum limit as adding the contributions from a smooth distribution of particles. This can be quite complicated if we need high order accuracy, which nominally would entail distributing a finite-sized model of the particle smoothly over a set of cells, taking the variation in the spacetime geometry into account. For these initial studies, we are not concerned with high order convergence; second order will suffice, and so we can avoid all these complications by simply smoothing a particle to a single containing cell.

To implement this, we demand that the stress energy tensor at some moment of coordinate time , integrated over the proper volume of the cell containing the geodesic, gives the same energy/momentum that an observer in the reference frame of the simulation would measure the particle to have. Here, the relevant observer is that traveling along the unit timelike vector normal to surfaces, and the corresponding proper volume element is , where is the determinant of the spatial metric (not to be confused with the metric function used earlier in the discussion of plane wave solutions), and is the coordinate volume element. A straightforward calculation shows that the contribution of a single particle with 4-momentum to the effective continuum stress energy tensor of the cell containing it is

 (c)Tab=1√hΔVℓaℓb(−ℓdnd) . (30)

To check this, recall that an observer with 4-vector measures the energy density of the stress tensor to be , and measures the total energy of the particle to be . The net continuum stress energy tensor is then just obtained by summing the contribution from all particles in the cell. Referring back to the schematic form of the particle stress tensor written in Eq. (22), the equation above defines exactly what we mean by . Note that it is not a constant, and its particular value is affected by our choice of affine parameter, defined to let us interpret as the physical, 4-momentum of the particle. For a set of timelike particles, again, is the 4-momentum of a particle of rest mass and 4-velocity .

In our code, we discretize the Einstein equations using a vertex-centered mesh. To obtain a second order accurate representation of the vertex-centered stress tensor, we take the above computed for each particle, and distribute it to the surrounding vertices, linearly weighting the contribution to each vertex based on the distance of the geodesic from it.

#### iii.5.2 An efficient evolution scheme within the Berger and Oliger time sub-cycling algorithm

The BO AMR algorithm Berger and Oliger (1984) uses a grid hierarchy consisting of nested grids, where finer resolution (child) grids are entirely contained within coarser resolution (parent) grids. A collection of grids at the same resolution is called a level. Hyperbolic differential equations discretized on such a hierarchy are evolved in time using the following algorithm (for a more detailed description and pseudo-code see Pretorius and Choptuik (2006), for example). For simplicity, assume the refinement ratio between successive levels is two, and the same Courant-Fredrichs-Lewy (CFL) factor is used on all levels (as used in our code). Beginning from the same starting time, on any given level, one time step of size is taken on all the grids at that level before two times steps of size are taken on the grids of the child level. This rule is applied recursively.

The reason for this scheme is so that the solution from the parent level can be used to set boundary conditions at child level interior boundaries, via interpolation in time from the parent solution. [At computational domain boundaries the relevant partial differential equation (PDE) boundary conditions are always applied.] Also, this time sub-cycling is optimally efficient for hyperbolic PDEs in the sense that each level is evolved with the same CFL factor (as opposed to schemes were all levels, or a single level but with non-uniform grid cells, are evolved with the same time step; the effective CFL factor for coarse levels/cells in such an approach could be much smaller than necessary for stable evolution). Note also that on parent levels the unigrid evolution is applied everyone, even at regions where a higher resolution child grid is available; the subsequent solution from the child grid evolution is injected back into the parent level after they both are in sync again, so the discrete solution after each global time step is always single-valued. This might seem like a waste of computational resources, though it is a relatively minor expense (most of the computation happens on the finest levels covering any region), and significantly simplifies development of AMR capable codes in that the underlying PDE evolution scheme can be almost completely ignorant of the mesh hierarchy (it only needs to know which boundaries are physical vs interior, and for interior boundaries it simply leaves the corresponding points untouched).

We would now like to include our geodesic integration within the BO AMR algorithm, and following the spirit of the algorithm, we want to solve the coupled Einstein-particle equations on each grid in a manner which relies minimally on knowledge of what mesh hierarchy the grid is part of. However, the problem here is because the geodesics are lines through the spacetime, and are integrated via ODEs, they do not have the same natural multiscale representation on a BO mesh hierarchy that continuum functions, such as the metric or stress tensor, have. Yet, we still want the solution to the geodesics to be computed using metric values from the finest mesh containing the geodesic. How then do we evolve geodesics, and use them to define the stress energy tensor, within the recursive, time sub-cycling algorithm? One simple option is just restart each geodesic for every finer level containing it, with the last then giving the most accurate solution that is kept. The problem with this is it will be very computationally inefficient for a deep hierarchy, as we do not have a multiscale representation of the set of geodesics; i.e. each geodesic is integrated times if the depth of the hierarchy is at the location of the geodesic. (This problem is mitigated for mesh-based PDE evolution because of the multiscale representation—effectively the number of mesh points where the PDEs are multiply-integrated drops as per level if the refinement ratio is , and there are spatial dimensions.)

One workaround would be to introduce a multiscale sampling of the continuum matter distribution where, say, the geodesic number density per cell is kept fixed going from level to level, and a coarse-level geodesic is some average of fine-level geodesics. This would complicated the structure of the AMR driver significantly.

The solution that we take instead is to adapt the scheme proposed in Pretorius and Choptuik (2006) for evolving a system of elliptic-hyperbolic PDEs within the BO time-stepping framework. A similar problem to the above arises for the elliptic equations, and without going into details here, the solution is to not solve the elliptic equations when descending the recursive tree (going from course to fine), when the hyperbolic equations are solved. Rather, then values for the elliptic variables are extrapolated from prior time levels, and instead the elliptics are solved when ascending the recursive tree (when hyperbolic variables are injected from fine to coarse levels). The way this algorithm is adapted to particles is when descending the recursive tree and the Einstein equations are solved, the stress energy tensor used on the right hand side in each cell is either (a) extrapolated from past time levels in regions where finer levels exists, and in these cells no geodesics are integrated, or (b) the geodesics within the cell are integrated (this is thus the finest level containing them) and used to compute the stress energy tensor for the neighboring vertices as described above. When ascending the tree, the most accurately available stress energy tensor values are injected back up to coarser levels.

The reason this approach was simpler for us is the AMR driver code we use (PAMR/AMRD code PAMR and libraries (2017)) already has infrastructure to handle the extrapolation: we simply define the stress energy tensor as if it were an elliptic variable in the code, and all the unigrid geodesic integration code needs to know is whether a given cell is the finest cell: if so, the geodesics within it are integrated and the stress tensor computed there. In the present code we save two past time levels, which is adequate to allow second order accurate extrapolation, and maintain overall second order accuracy of the evolution (the metric and geodesics are still evolved with order accurate Runge-Kutta, but the stress energy energy calculation discussed above is only order accurate, independent of the way we evolve the geodesics within the BO framework).

In the remainder of the main text, we will focus on an axisymmetric, null particle case. However, we have also tested this code for 3D (i.e. non-axisymmetric) calculations, as well as with massive (timelike) particles, as we demonstrate in the appendix with results from a simple inhomogeneous cosmology setup.

## Iv Results: Black Hole Formation in axisymmetric, head-on collisions

Here we present results from simulations of the head-on collision of two axisymmetric null dust-sourced plane gravitational waves. We select initial data to closely link to the Aichelburg-Sexl limit. As discussed in the introduction, Penrose first looked at the AS case. Assuming a mass based on the initial apparent horizon area is a lower limit to the mass the black hole settles down to, and hence the difference between the initial spacetime and apparent horizon mass is an upper limit to the gravitational radiation emitted during the collision, Penrose’s computation gives . D’Eath D’Eath (1978), and D’Eath and Payne D’eath and Payne (1992b) perturbatively explored the far-field region of the AS limit, and were able to directly estimate the gravitational wave emission, timelike compact object collisions, Sperhake et al. (2008) first collided black holes up to , and extrapolating the results to , found . This scenario was extended to in  Healy et al. (2016), where they estimated 666This might seem in mild tension with our quoted number, though looking at Figs. 9 and 10 of  Healy et al. (2016) suggests there is an effective systematic uncertainty associated with the different classes of initial data and codes they use that might warrant a slightly more conservative error estimate, and so we think there is not yet any significant indication that the two different approaches will not reach the same AS collision spacetime in their respective limits.. In East and Pretorius (2013), compact fluid stars with up to were collided, and it was found that .

Similar to Penrose’s first investigation where he found an apparent horizon at the moment of impact, earlier studies of null-radiation (gyraton) interactions were able to show apparent horizon formation Yoshino et al. (2007). (These authors considered more general models of gyratons including rotation Frolov and Fursaev (2005); the null-particle case maps to non-rotating gyratons.) Here we are able to follow the evolution of spacetime through the formation of an apparent horizon which eventually settles down toward a Schwarzschild black hole, together with gravitational radiation streaming away from the collision. Regarding the net gravitational wave energy emitted, as explained below, we are not able to characterize the initial decay of the waves close to the collision (where we can measure them) in a manner that allows extrapolation to infinity in order to accurately calculate the energy they contain. Instead then, we assume is the difference between the initial spacetime mass and late-time apparent horizon mass, the latter which we can compute accurately, and using a conservative upper bound of of energy in particles that escaped entrapment by the black hole as an additional source of uncertainty, obtain . As discussed more below, we expect this scenario to give that is a little below the AS limit by .

The structure of the remainder of this section is as follows. In Sec. IV.1 we describe the specific initial data we use, in Sec. IV.2 we discuss our gauge source function evolution equations, and in Sec. IV.3 we present the results from the simulations, including convergence tests.

### iv.1 Initial data

Here we give the particular form of a single, right () propagating null fluid wavefront we use for initial data in the numerical evolution (the left propagating front has identical form, with , as discussed in Sec. III.2). For simplicity, here we show the pulse centered at , though it is trivial to shift it to any desired starting location at . We begin with the metric in the form (8), and, since we are restricting to axisymmetric spacetimes here, :

 ds2=−dudv+2β(u)q(ρ)dudρ+dρ2+ρ2dθ2 . (31)

For the matter profile, we use a piecewise polynomial function in , and a Gaussian in . Defining , , with and constant parameters that define the scale of the profile,

 ρe(u,ρ) = f(u)g0(ρ), (32) f(u) = ⎧⎪⎨⎪⎩0,¯u<−1(¯u2−1)2,−1≤¯u≤10,¯u>1 (33) g0(ρ) = Ae−¯ρ2, (34)

where is a parameter controlling the amplitude of the wave (note here that the -extent of the pulse is , a factor of 2 different from the discussion in Sec. II.2). Equations (10) and (14) then give

 β(u)4πΔu = (35) q(ρ) = (36)

The ADM mass for this spacetime evaluates to

 ¯M=A8π15ΔuΔρ2. (37)

If we take the point-particle limit , keeping fixed, (31) becomes, for ,

 ds2=−dudv+8¯MρΘ(u)dudρ+dρ2+ρ2dθ2, (38)

where is the Heaviside step function. For =0, the metric in (31) always has by regularity, but if instead we define (38) as the metric including , then this is the AS solution if we identify its total energy with . Interestingly, as discussed in Sec. II.3, if we were to take (38) (with ) as a vacuum solution to the Einstein equations and directly use it in the ADM formula (19), the latter would give an ADM mass of . 777Note that for a pure vacuum case one could still use the metric ansatz (8) with in axisymmetry, but would then not impose the conditions (10) and (11) that otherwise seem to link to ; instead can then be considered the arbitrary function we choose to specify the longitudinal extent of the pulse. For the matter sourced spacetimes with the above energy profile, each line in the formula (19) contributes , but the first line is identically zero for vacuum spacetimes.

A further curious property of this solution is, transforming to Cartesian coordinates via (28), the lapse function is . The lapse becomes singular when , implying the hypersurface fails to be globally spacelike, and hence cannot be used as initial data in a Cauchy evolution fixed to this time coordinate. For the AS solution, this occurs to the left of the shock for ; for matter solutions this implies we cannot use initial data in these coordinates if . Of course, this is just a coordinate singularity, though what is curious about it is it seems to “anticipate” black hole formation in a collision: colliding two identical pulses each with mass , the total spacetime mass will be , and the hoop conjecture argues a black hole should form following collision if all the matter is focused inward within a hoop of radius , the limiting estimate for for each pulse to have a well defined lapse.

An issue with the above coordinate system that relates to computational cost is that to the left of the pulse there is non-trivial structure in the metric within a region of size of the axis, all the way to . Numerical experiments show that for long-term stable evolution this region needs to be resolved with essentially the same resolution as the inner part of the domain where the collision occurs, even though the part of the spacetime is simply Minkowski, with no dynamics. This is somewhat wasteful, and suggests we should transform back to the representation of Minkowski here. Though, as discussed in Sec. II.2, we suspect one cannot transform to exactly this representation after the wave while also maintaining as the form of Minkowski ahead of the wave (which is desired for the simplicity of initial data construction) and not introducing a divergence as within the plane of the matter. However, we can partly transform back to post-shock, in particular in a region about the axis. In a sense, this can spread the non-trivial structure over a larger region in that then needs much less resolution to resolve. To do this, before transforming to Cartesian form from double null coordinates, we rescale following a generalization of (II.2):

 v=¯v+G(ρ)L(u). (39)

This transforms (31) to

 ds2=−dud¯v+(2βq−G′L)dudρ−GL′du2+dρ2+ρ2dθ2, (40)

where, for simplicity, we do not show the functional arguments, and prime () denotes a derivative with respect to the function’s argument. To not affect the form of the solution within the wave or ahead of it, we choose for some , and let smoothly increase to one over the region (in the numerical calculation we use a piece-wise order polynomial function for this). For we choose , where , , and is a transition function that is one for , and smoothly (again via a piece-wise order polynomial in the numerical calculation) goes to zero over the region . The axis resolution issues arose from the term in the metric; with this transformation we therefore eliminate this term in the region , and slowly reintroduce it over a longer scale in controlled by . The transformation creates a new piece of the metric within , though its size can likewise be controlled by .

### iv.2 Gauge evolution

At the initial time, and for a short time thereafter, the gauge source functions are simply set to the superposition of those of the exact solutions; i.e. (26) evaluated with the relevant version of (40), after each is shifted in or to have the pulses at the desired starting positions, and then transformed to Cartesian coordinates. This simple gauge prescription actually works even through collision for weak pulses that are not close to forming black holes, but for strong pulses leads to a coordinate singularity some time after the interaction. Therefore, at a specified time before the collision we smoothly transition the gauge, over a chosen time period, and within a chosen spatial volume of the origin, to the damped harmonic gauge described in Choptuik and Pretorius (2010) with amplitude parameter (which is equivalent to A15 in Lindblom and Szilágyi (2009) with and ).

The evolution is not particularly sensitive to exact parameter values, with ( is the total mass of the spacetime), and as long as the transition in time takes place within of the collision, and within a volume that comfortably encloses the apparent horizon. For the simulations presented below, we use , transition to damped harmonic within (with a order piece-wise polynomial in time), and within a coordinate sphere of radius of the origin. From to the gauge smoothly (again with a order polynomial) transitions to the superposed initial data gauge, and remains that from to our outer boundary at . The simulations are run till , and the two pulses are set apart so that the interaction beings within .

### iv.3 Case study

Here we present the results from simulations of the head-on collision of two plane-fronted gravitational waves, sourced by null particles. We choose initial conditions to connect to the shock-front AS collision limit; i.e. large amplitude waves that promptly form a single encompassing black hole, and whose characteristic width is smaller than its transverse length, the latter corresponding roughly to the size of the eventual black hole. Specifically, each initial pulse is given a profile of the form (32)-(34), with , and . (This is close to the maximum compaction in we can choose in these coordinates; is a less than the rough estimate of given in Sec. IV.1, which would be exact for a constant density with sharp -cutoff.) The total mass of the spacetime is . The facing edges of each pulse are set apart, so that they begin to interact within , and an apparent horizon is first detected at .

Choosing large amplitude data with results in essentially all the matter falling into the black hole (some particles in the Gaussian tail do escape), and so the post-collision spacetime closely connects to the vacuum AS collision problem. The true limit would take both and to zero, with the matter distribution approaching a delta function. As discussion above, we cannot make arbitrarily small with this class of initial data; however, for , a black hole forms promptly, and the gravitational interaction that happens outside the black hole transverse to the collision is then, by causality, going to be independent of how small is within the black hole (the asymptotic fall-off in outside the matter region is insensitive to for a sequence with fixed and ). So in that sense, in terms of , we believe we are quite close to the AS limit. Our initial data allows us to make arbitrarily small (for fixed and ). However, the smaller is, the more computationally expensive the simulations become, as this length scale needs to be resolved. In the limit , the geometry becomes shock-like, and this feature will persist (at least along the leading edge of the shock) post-collision. In terms of starting to resolve a shock-like feature, we are quite far from the AS limit. (In a relative sense, any finite is always infinitely far from ; a more physical measure might be the width of each initial wave divided by the final black hole diameter , which here equals .) Regarding net gravitational wave emission, we ran a preliminary survey decreasing to as low as , and extrapolating, this suggests the case underestimates the AS limit net energy emission by . However, these runs were performed with the same resolution (see below) as the case presented in detail here, which effectively means successively worse resolution for smaller . Hence, the number should be considered a rough estimate. Again, it would take significant computational resources to explore the AS limiting sequence in detail and accurately, and we leave that for future work.

Another practical reason for studying this point in parameter space for this first application of the null particle code, is the black hole that forms cleanly hides extreme focusing onto the axis with these axisymmetric profiles. The matter is pressureless (except for the effective anisotropic pressure that arises in multistream regions), and the initial data has no angular momentum about the symmetry axis, so even if no singularities form, nothing prevents the distribution from focusing to very small length scales, and to maintain convergence when this happens requires computationally expensive, deep mesh hierarchies.

Our base case resolution (which we label with ) is such that the coarsest level of the AMR hierarchy covers the entire domain , with cell-widths of size , and has up to seven levels of 2:1 refinement, giving a minimum possible cell width of . The hierarchy is dynamically generated via truncation error estimates. We evolve until , and when measuring some property of the solution restrict the measurement domain to be out of causal contact with the boundary (assuming unit light speed), to mitigate inconsistencies that may arise due to our outer boundary conditions (which are given by superposed exact null wave solutions). We sample each matter distribution with geodesics, which prior experimentation has shown is adequate to have the error be smaller than the discrete mesh truncation error (and we also demonstrate that in Fig. 3).

To check convergence—see Figs. 3 and 4 for norms of the Einstein and geodesic equation constraints respectively—and compute error estimates, we also ran simulations with 1.6 and finer ( and ) base level resolution, adjusting the truncation error threshold in the AMR algorithm to generate finer levels according to the expectation of overall order convergence. We also in tandem change the number of particles by when the spatial resolution is scaled by . One factor of keeps the particle density, hence error the same in each cell, the additional factor of is then to further increase/decrease the resolution of the sampling within the cell to match the scaling of the mesh-based truncation error. This highlights how computationally expensive a particle-based code is to achieve high accuracy, and it is not clear in this respect that it offers any advantage over directly discretizing the Boltzmann equation in phase space (our main reason for going the former route is simplicity of implementation).

See Fig.5 for a plot illustrating the collision in terms of the gravitational radiation produced, measured with , the Newman-Penrose scalar with tetrad adapted to measure outgoing radiation propagating in the radial coordinate direction . Evident in the figure is that the radiation, roughly speaking, appears to be composed of two components. One is the longer wavelength feature expected from the dominant quadrupolar quasi-normal mode (QNM) oscillation of the black hole created during the collision. The other can be associated with the initial wave front, the inner part of which is trapped by the black hole, while the remainder propagates about the black hole ad infinitum, forming the concentric, circular-like features of characteristic wavelength close to . In the AS limit where , the leading edge of this feature presumably remains shock-like.

The tetrad used to compute is completed with vectors tangent to the sphere ; we define the corresponding angles on the sphere so that measures the angle from the axis (so the simulation plane is ), () coincides with the plane (), and () coincides with the plane (). Thus, on the positive axis, this tetrad will be exactly that used to define (21), and similarly on the negative axis but flipped to measuring radiation propagating outward along the direction. This then gives us one way to easily understand the feature in Fig. 5 that vanishes on the collision axis. Even though the presence of the matter forces to vanish on the axis in the initial data [see the discussion after equation (21)], most of the matter region immediately falls into the black hole, and essentially all the radiation that reaches the axis can be traced back to vacuum regions along the initial wave fronts. Here, the polarization relative to an direction propagation vector [indicated by the and terms in (21)] required for the wave to be axisymmetric about the collision axis is such that when a ring of waves focuses to the axis, their sum must be zero.

Despite the fact that symmetry forces to be exactly zero on the axis, an interesting property of the radiation is how strongly beamed it is about the axis (this is consistent with perturbative calculations of the scattering problem, and also noted in high speed collisions of black holes Sperhake et al. (2008); Healy et al. (2016)). This is further illustrated in Fig. 6 showing at three radii (, , and ), and two angles relative to the collision axis (at a right-angle , and close to it with ). Figure 7 is the same data plotted on a logarithmic scale, to highlight the decay of the waves. Interestingly, the latter plot shows that, transverse to the collision axis, the decay rate is broadly consistent with the decay of the least damped QNMs, but less so near the axis. Also, the time-shifts required to align the waves on the plots indicate different effective propagation speeds along the axis versus transverse to it; how much of this is simply due to gauge (time slicing, how the coordinate relates to some geometric radius, etc.) as opposed to slightly different geometries the different parts of the wave are propagating through is unclear.

The above observations suggest at least part of the outgoing radiation might be not be captured as a sum of Schwarzschild QNMs (which do not form a complete basis), as in a sense the initial wave-fronts are part of the initial data, and not “produced” by a perturbed black hole. However, since the decay of such a wave-front propagating about the black hole will be controlled by the unstable photon orbit, similar to QNMs, at least at late-times there may be no practical distinction between what is a QNM versus what is a remnant of waves from the initial data.

In Fig. 8 we show an estimate of the black hole mass that forms as a function of time, calculated by computing the area of the apparent horizon. In Fig. 9 we show the ratio of proper equatorial to polar circumference of the horizon, illustrating the early time dynamics and subsequent decay to a Schwarzschild black hole. From the area, we estimate the mass of the remnant black hole to be . Given that the particles that escaped being trapped by the black hole collectively contain less than of energy, we can infer that the energy emitted in gravitational waves is of the initial spacetime mass .

Integration of measured on coordinate spheres gives order of magnitude consistent values, specifically , , and integrated on spheres of radius , , and , respectively, from the run. However, this does not seem to be sufficiently far into the wave-zone, in that we do not yet see the expected decay in the waveform, which would allow extrapolation to for an accurate estimate. This mostly seems to be due to near-axis beaming of the radiation, illustrated in Fig. 6. In addition, a couple of other effects, one gauge, the other numerical, hinder trying to fit the energy to a more complicated series expansion, and so for now we will take the apparent horizon based estimate as the more accurate measure of total energy radiated. The numerical issue is related to the short wavelength component of the radiation, proportional to the width of the initial data, and is responsible for the concentric rings evident in Fig. 5. This is not well resolved by the grid in the wave extraction zone, and numerical dissipation attenuates this feature more rapidly than the longer wavelength components of the wave. The gauge issue is related to the significant gauge dynamics we have as we transition from the initial data coordinates to the damped harmonic gauge post collision. This is more pronounced closer to the origin, and is exacerbated by our naive procedure of measuring radiation on coordinate spheres, and simply using as their geometric area (which is only correct asymptotically). Certainly a combination of improved resolution, farther extraction, and a more geometrically sound construction of extraction spheres could alleviate these issues, though that will take considerable effort and computational resources, and we leave it to future work.

## V Conclusions

We have described a formalism for studying the ultra-relativistic scattering problem using plane-fronted distributions of null particles as a matter source. We have developed a numerical code based on this, and as a first application presented a study of black hole formation in head-on axisymmetric collisions, with parameters of the particle sources chosen to give a post-collision spacetime close to that expected to be produced by the collision of two Aichelburg-Sexl shock waves. We find results broadly consistent with prior studies of this limit, whether perturbatively, or via full numerical solution but using finite boost compact object models of particles. Specifically, based on the area of the resultant black hole that forms, we infer of the initial mass of the spacetime is radiated as gravitational waves during the collision. We leave it to future work to explore in detail where this particular case fits into a limiting sequence reaching the AS limit, though we estimate in terms of net gravitational wave emission this is slightly below the AS limit by (hence the quoted value of in the abstract as the implied value for the AS limit).

In terms of other future directions, there are numerous avenues that can be pursued, many of them outlined in the introduction, so we will not repeat that discussion here. Rather, we will briefly mention a few outstanding issues in the code that would need to be addressed before the full breadth of applications could be tackled. First is simply the computational expense of the simulations, in part due to the scaling of particle number required to achieve a desired level of accuracy (which will be more severe for 3D applications than the axisymmetric 2D case presented here), and in part that the collisionless particle model generically allows focusing to caustic regions. In principle, the latter problem could be alleviated by including some form of self-interaction between the particles that produces an effective pressure; however, it is unclear to what extent that kind of matter could be used as a consistent source for plane-fronted gravitational wave spacetimes. A second issue is that earlier studies (in particular East and Pretorius (2013)), and preliminary investigations with this code, show the gauge conditions (i.e. source function evolution equations) that work well for prompt black hole formation are not adequate for slightly weaker interactions where the spacetime is highly dynamical but does not immediately (or will not ever) form horizons. What typically happens with the current gauges is a coordinate singularity forms. This is also a problem that has plagued attempts to study vacuum critical collapse (see e.g. Hilditch et al. (2017)), and if the critical solution is universal it would not be surprising if a single class of novel gauge condition could solve the coordinate issues for both these applications.

###### Acknowledgements.
F.P. acknowledges support from NSF grant PHY-1607449, the Simons Foundation, and the Canadian Institute For Advanced Research (CIFAR). This research was further supported in part by the Perimeter Institute for Theoretical Physics. Research at the Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development Canada and by the Province of Ontario through the Ministry of Research, Innovation and Science. Computational resources were provided by XSEDE under grant TG-PHY100053 and the Perseus cluster at Princeton University.

## Appendix: Application to Inhomogeneous Cold Matter Cosmology

In this appendix, we briefly present some results of applying the methods discussed in the main text for evolving the Einstein equations, coupled to particle distributions, to the case of inhomogeneous cold matter in an expanding universe. We do this both to demonstrate that these same methods can be applied to massive (timelike) particles, and in three spatial dimensions, as well as to further validate the code by comparing it to known results. In particular, we repeat a calculation from East et al. (2018) where we start with a (timelike) dust-filled, expanding universe (i.e. the Einstein-de Sitter model) and include some initially small inhomogeneities that are all at a wavelength that is four times the initial Hubble radius, and with velocity given by the Zel’dovich approximation. (In particular, this is the case labeled in East et al. (2018); see that reference for details.)

To construct initial data, we begin with a solution to the constraint equations obtained using the code of East et al. (2012b), as described in East et al. (2018), which uses a fluid description of the matter. We then create a particle distribution that is consistent with the density field in the following manner. We begin with a uniformly spaced lattice of particles that we then perturb from their positions according to the Zel’dovich approximation, as is typically done in N-body calculations. In particular, we apply the shift in position given by Eq. (31) in Chisari and Zaldarriaga (2011). We assign the velocity of the particle by interpolating the fluid velocity field to the particle position.

By default the particles would have uniform masses given by , where is the initial homogeneous density, is the volume of the domain, and is the number of particles. However, we then calculate a small, nonlinear correction to the mass of each particle by first finding the density given by , with calculated from (30) using uniform particle masses, and taking the ratio with the desired density field from the fluid representation at each particle position . We then rescale the mass of each particle by this factor .

The evolution is performed as described in the main text, except that our domain is three-dimensional and periodic (particles that exit the domain at one boundary are wrapped around to the opposite boundary), and the particles have non-zero rest mass. In Fig. 10, we demonstrate the convergence of the constraints (27) by performing the calculation at several resolutions ranging from to grid points covering the domain. As in the main text, the error from finite particle number is subdominant to the grid discretization error for the parameters considered here.

We can also compare the evolution of the density contrast, which increases and becomes nonlinear (as signaled by the divergence of the under and overdensities) as the inhomogeneities enter the horizon, to the results from East et al. (2018), which were obtained by evolving a pressureless fluid. As long as we restrict to times before multistream regions form, the two treatments should give the same answer. As shown in Fig. 11, the density contrast from the two cases is indeed a good match.

The advantage of the particle treatment is that it allows one to evolve through the formation of multistream regions that will arise during structure formation, and thus be more directly comparable to Newtonian N-body simulations. However, we leave a study of this to future work.

## References

You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters