Slicing conditions for axisymmetric gravitational collapse of Brill waves

TL;DR

This paper introduces a new 'quasi-maximal' slicing condition for numerical relativity that improves the evolution of near-critical Brill wave data, outperforming traditional gauges in stability and behavior.

Contribution

A novel slicing condition based on 1+log slicing with an added source term from maximal slicing, enhancing numerical stability for Brill wave simulations.

Findings

The new slicing condition demonstrates improved stability for near-critical Brill wave evolutions.

Constructed spacetimes show black holes settling into Schwarzschild solutions.

Gauge-independent quantities confirm the robustness of the new method.

Abstract

In numerical relativity, spacetimes involving compact strongly gravitating objects are constructed as numerical solutions of Einstein's equations. Success of such a process strongly depends on the availability of appropriate coordinates, which are typically constructed dynamically. A very robust coordinate choice is a so-called moving puncture gauge, commonly used for numerical simulations of black hole spacetimes. Nevertheless it is known to fail for evolving near-critical Brill wave data. We construct a new `quasi-maximal' slicing condition and demonstrate that it exhibits better behavior for such data. This condition is based on the 1+log slicing with an additional source term derived from maximal slicing. It is relatively simple to implement in existing moving puncture codes and computationally inexpensive. We also illustrate the properties of constructed spacetimes based on gauge-independent quantities in compactified spacetime diagrams. These invariants are also used to show how created black holes settle down to a Schwarzschild black hole.

Equations

$${\gamma }_{ij}={\psi }^4\left[e^{2q}\left(\mathrm{d}{\rho }^2+\mathrm{d}{z}^2\right)+{\rho }^2\mathrm{d}{\phi }^2\right],$$

$$q(\rho ,z)=A{\left(\frac{\rho }{\sigma }\right)}^2{e}^{-({\rho }^2+z^2)/{\sigma }^2}.$$

$$\Delta \psi +\frac{1}{4}\left({\partial }_{\rho \rho }q+{\partial }_{zz}q\right)\psi =0,$$

$$\mathrm{d}{s}^2=\left(-{\alpha }^2+{\beta }_i{\beta }^i\right)\mathrm{d}{t}^2+2{\beta }_i\mathrm{d}t\mathrm{d}{x}^i+{\gamma }_{ij}\mathrm{d}{x}^i\mathrm{d}{x}^j.$$

$$\left({\partial }_t-{\mathcal{L}}_{\beta }\right){\gamma }_{ij}$$

$$\left({\partial }_t-{\mathcal{L}}_{\beta }\right)K_{ij}$$

$${}^{(3)}R+K^2-K_{ij}{K}^{ij}$$

$${\mathrm{D}}_j{K}^{ij}-{\mathrm{D}}^{i}K$$

$$\phi$$

$$K$$

$${\bar{\gamma }}_{ij}$$

$${\bar{A}}_{ij}$$

$${\bar{\Gamma }}^i$$

$$\left({\partial }_t-{\mathcal{L}}_{\beta }\right)\phi$$

$$\left({\partial }_t-{\mathcal{L}}_{\beta }\right){\bar{\gamma }}_{ij}$$

$$\left({\partial }_t-{\mathcal{L}}_{\beta }\right)K$$

$$\left({\partial }_t-{\mathcal{L}}_{\beta }\right){\bar{A}}_{ij}$$

$$\left({\partial }_t-{\mathcal{L}}_{\beta }\right){\bar{\Gamma }}^i$$

$${\partial }_{t}K=-{\mathrm{D}}^i{\mathrm{D}}_i\alpha +K_{ij}{K}^{ij}\alpha =0,$$

$$\left({\partial }_t-{\mathcal{L}}_{\beta }\right)\alpha =-2\alpha K.$$

$$\left({\partial }_t-{\mathcal{L}}_{\beta }\right){\beta }^i=\xi {\bar{\Gamma }}^i-\eta {\beta }^i,$$

$${\mathcal{K}}_{\mathrm{Schwarzschild}}=\frac{48M^2}{R^6},$$

$$M_{\mathcal{K}}(R)=\sqrt{\frac{\mathcal{K}{R}^6}{48}},$$

$${\bar{\rho }}^2={\eta }_{\mu }{\eta }^{\mu }={\gamma }_{yy}{x}^2,$$

$$4{\bar{\rho }}^{\prime 2}={\left|\nabla {\bar{\rho }}^2\right|}^2=({\mathrm{D}}_i{\bar{\rho }}^2)({\mathrm{D}}^i{\bar{\rho }}^2)-4(K_{ij}{\eta }^i{\eta }^j{)}^2.$$

$$\chi =\frac{{\bar{\rho }}^{\prime 2}}{{\bar{\rho }}^2}.$$

$$M_{\bar{\rho }}=\frac{{\bar{\rho }}^2-{\bar{\rho }}^{\prime 2}}{2\bar{\rho }}.$$

$$0=\left({\partial }_t-{\mathcal{L}}_{\beta }\right)\left[{\gamma }^{ij}{\mathrm{D}}_i{\mathrm{D}}_j\alpha -K_{ij}{K}^{ij}\alpha \right],$$

$${\gamma }^{ij}{\mathrm{D}}_i{\mathrm{D}}_j\left[({\partial }_t-{\mathcal{L}}_{\beta })\alpha \right]-K_{ij}{K}^{ij}\left[({\partial }_t-{\mathcal{L}}_{\beta })\alpha \right]=$$

$$-\left[\left({\partial }_t-{\mathcal{L}}_{\beta }\right){\gamma }^{ij}\right]{\mathrm{D}}_i{\mathrm{D}}_j\alpha +{\gamma }^{ij}\left({\partial }_t{\Gamma }_{ij}^k\right){\partial }_k\alpha$$

$$-\left({\gamma }^{ij}{\mathrm{D}}_i{\mathrm{D}}_j{\beta }^k\right){\mathrm{D}}_k\alpha -{\beta }^j{R}_j^i{\mathrm{D}}_i\alpha +\alpha \left({\partial }_t-{\mathcal{L}}_{\beta }\right)\left(K_{ij}{K}^{ij}\right).$$

$$W=({\partial }_t-{\mathcal{L}}_{\beta })\alpha$$

$${\mathrm{D}}^i{\mathrm{D}}_{i}W-K_{ij}{K}^{ij}W=$$

$$-2\alpha K^{ij}{\mathrm{D}}_i{\mathrm{D}}_j\alpha +{\gamma }^{ij}\left({\partial }_t{\Gamma }_{ij}^k\right){\partial }_k\alpha -\left({\gamma }^{ij}{\mathrm{D}}_i{\mathrm{D}}_j{\beta }^k\right){\mathrm{D}}_k\alpha$$

$$-{\beta }^j{R}_j^i{\mathrm{D}}_i\alpha +\alpha \left(2\dot{K_{ij}}{K}^{ij}+4\alpha K_j^i{K}_i^k{K}_k^j\right),$$

$$({\partial }_t-{\mathcal{L}}_{\beta })\alpha =-2\alpha K+\kappa W_a,$$

$$\psi \left(x^0,x^1\right)=1+\sum _{k=0}^{N-1}\sum _{l=0}^{M-1}{C}_{kl}{B}_k^0\left(x^0\right)B_l^1\left(x^1\right).$$

Full text

Keywords: Numerical relativity, gravitational collapse, axial symmetry

1 Introduction

During the last decade numerical relativity became tremendously successful at simulating many astrophysically important processes in which general relativity plays a crucial role. A significant portion of this progress can be attributed to the discovery of coordinates suitable for such wildly dynamic spacetime geometries. Present codes can handle mergers both of purely vacuum black hole binaries and neutron star binaries (followed by a collapse). It may thus seem surprising that a certain theoretically important situation — gravitational wave packet collapse — still defies these codes.

Numerical treatment of gravitational waves collapsing to form a black hole dates back to the early age of numerical relativity. Eppley in his pioneering work [1] successfully constructed what is known as the Brill waves [2] — a family of vacuum initial data describing axisymmetric wave packets at the moment of time symmetry — and located apparent horizons in them. Another breakthrough was the paper of Abrahams and Evans [3] which used another family of axisymmetric data called the Teukolsky waves [4]. They found evidence of critical behavior — discrete self-similarity and power-law scaling of various quantities close to the threshold of black hole formation — analogous to the results of Choptuik for a scalar field [5, 6].

Since then, as computer performance and numerical methods improved, several attempts were made to investigate the nonlinear regime of the Brill data. Alcubierre et al. [7] used a combination of the BSSN evolution system and maximal slicing to put rough bounds on the critical value of the amplitude parameter. That result was confirmed in [8] using a mixed elliptic-hyperbolic reduction of the evolution equations, again with maximal slicing. Ultimately, insufficient performance of contemporary computers constrained what could be achieved with these attempts.

As the moving puncture gauge [9, 10] showed its strength for black hole mergers, Hilditch et al. [11] attempted to use it also for gravitational wave collapse. One of their discoveries was the pathological behavior of this gauge for near-critical Brill waves, producing what they conjectured to be coordinate singularities. Recently [12, 13] they managed to avoid these problems with a new pseudospectral code using a generalized harmonic formulation.

In this paper, we follow up on the work of [11], using the same initial data and similar numerical techniques (although a completely different code). We analyze the way in which the moving puncture gauge — a combination of the ‘1+log’ slicing and the ‘$\Gamma$-driver’ shift condition — breaks, and modify the slicing condition. Our modification takes the form of an extra source term pushing the constant-time hypersurfaces closer towards maximal slicing, which is known to guarantee smooth solutions and has shown promising results in earlier work [7, 8, 14]. The result is a new slicing condition — which we call ‘quasi-maximal slicing’ — that no longer exhibits the problems reported in [11]. While quasi-maximal slicing involves solving an elliptic equation, solver accuracy only affects the distance from the preferred gauge and does not give rise to ADM constraint violations — in contrast to the usual implementations of maximal slicing.

Using this slicing, we illustrate features of gravitational wave collapse that are present even farther away from the critical amplitude, such as formation of trapped surfaces and the event horizon, and construct compactified diagrams of considered spacetimes. These were previously not available, since earlier techniques are not able to handle long-term simulations easily.

This paper is laid out as follows. In Section 2 we briefly summarize the continuum equations on which our numerics is based. In Section 3 we introduce our modified slicing condition. Section 4 describes the key points of the numerical codes used and why we trust their output. Finally in Section 5 we present the results for long-term evolution of the considered class of initial data. We use geometrized units $G=c=1$ in this paper.

2 Field equations

2.1 Brill wave initial data

The Brill waves [2] are a family of vacuum axially symmetric initial data defined by the spatial metric in cylindrical coordinates $\left\{\rho,z,\varphi\right\}$

[TABLE]

and the condition of time symmetry, which implies initially vanishing extrinsic curvature. Here $\psi$ is the conformal factor and $q=q(\rho,z)$ is the so-called ‘seed function’, which needs to have certain regularity and decay properties [2], but can otherwise be chosen arbitrarily. For ease of comparison, we use the same form of $q$ as was used in [11], specifically

[TABLE]

The parameter $A$ determines the amplitude of the waves, with $A=0$ being flat space. We only consider non-negative values of $A$ in this paper, though it is worth noting that $A<0$ also produces valid initial data. The critical point $A^{*}$ is the smallest value of $A$ for which a gravitational singularity is formed. For the above choice of $q$ it is bracketed by our simulations to $A^{*}\in\left[4.69,4.7\right]$, which is compatible with tighter bounds given in [13].

The scale $\sigma$ (with dimension of length) fixes the units of all physical quantities. In our code we set $\sigma=1$, so all numerical values of time, length or mass are assumed to be in units of $\sigma$. This allows us to compare our results directly with [11], where the same choice is made. For quantities that do not clearly have dimension of length we print the appropriate power of $\sigma$ explicitly.

Due to time symmetry the momentum constraints are trivial, so to construct the data we only need to solve the Hamiltonian constraint. It reduces to

[TABLE]

a linear elliptic equation for the conformal factor $\psi$, where $\Delta$ is the flat space Laplacian. Our code for solving this equation is described in Section 4.1.

2.2 Evolution system

We follow the standard 3+1 splitting procedure, with the 4-dimensional metric $g_{\mu\nu}$ decomposed into the spatial metric $\gamma_{ij}$, the lapse $\alpha$ and the shift vector $\beta^{i}$

[TABLE]

The vacuum Einstein equations can then be written as a set of evolution equations

[TABLE]

and a set of constraints

[TABLE]

Here $K_{ij}$, defined by (5a), is the extrinsic curvature of the spatial slices, $\mathcal{L}_{\beta}$ is the Lie derivative along $\beta^{i}$, $\mathrm{D}_{i}$, $R_{ij}$ and ${}^{(3)}R$ are respectively the covariant derivative, the Ricci tensor and the scalar curvature associated with $\gamma_{ij}$.

Introducing new evolved variables

[TABLE]

where $\bar{\Gamma}_{jk}^{i}$ are the Christoffel symbols associated with $\bar{\gamma}_{ij}$, and using the constraints, we can write the evolution equations as

[TABLE]

where TF denotes the trace-free part of the bracketed expression. The conformally related metric $\bar{\gamma}_{ij}$ and its inverse $\bar{\gamma}^{ij}$ are used to raise and lower indices on all the quantities with a bar. Terms involving $\mathrm{D}_{i}$, $R_{ij}$ and ${}^{(3)}R$ are now assumed to be computed from $\bar{\gamma}_{ij}$ and $\varphi$ and their derivatives in a straightforward manner.

This is known as the BSSN [15, 16] formulation, which (given appropriate gauge choice) is known to be strongly hyperbolic, and has a solid track record in black hole simulations.

One of the slicing conditions we consider is the maximal slicing, defined by requiring the volume elements associated with normal observers to remain constant. This implies that the trace of the extrinsic curvature $K$ is identically zero at all times. Provided that the initial slice satisfies $K=0$, equation (5fghc) immediately gives us

[TABLE]

a linear elliptic equation for the lapse $\alpha$ that needs to be solved at each step of the evolution. For its many desirable properties [17], this slicing has been used since early times of numerical relativity, and specifically for Brill waves e.g. in [7, 8]. Its main disadvantage is the fact that solving an elliptic equation at each time step is usually impractical in 3D and tends to be very resource-intensive even in 2D.

According to [18], maximal slicing is only well-posed when the constraint $K=0$ is enforced, which is easily achieved with BSSN by just not evolving $K$.

The other common lapse choice we make use of is the 1+log slicing [9], given by a hyperbolic evolution equation for the lapse

[TABLE]

This condition is intended to approximate the main advantages of the maximal slicing at a lower computational cost.

For the shift we mainly use the $\Gamma$-driver condition (original version introduced in [10], the flavor we use is from [19])

[TABLE]

a hyperbolic equation related to elliptic minimal-distortion conditions. Here the dimensionless parameter $\xi$ determines the velocities of the longitudinal and transversal gauge waves, while $\eta$ is a damping parameter (with a dimension of inverse length) that seems to be necessary to avoid shock formation in the shift. In some cases we also set the shift to be identically zero.

Together 1+log slicing and $\Gamma$-driver shift are known as the moving puncture gauge and are widely used in numerical relativity, among other thing for black-hole spacetimes.

2.3 Invariants

To extract coordinate-independent information from our simulations we look at spacetime invariants. Since the end-state of the super-critical runs has to be a Schwarzschild black hole, the invariants’ behavior in Schwarzschild is of special interest to us.

The most prominent scalars are those constructed from curvature, e.g. the Kretschmann scalar $\mathcal{K}=R_{\mu\nu\alpha\beta}R^{\mu\nu\alpha\beta}$ (where $R_{\mu\nu\alpha\beta}$ is the four-dimensional Riemann tensor). For a Schwarzschild black hole with mass $M$ we have

[TABLE]

where $R$ is the usual areal radius. This allows us to define the ‘Kretschmann mass’ as

[TABLE]

where $\mathcal{K}$ and $R$ are evaluated at some point in the equatorial plane $z=0$.

One practical issue with $\mathcal{K}$ is that it contains second derivatives of the evolved quantities and its stationary limit (on which the mass estimate (5fghm) is based) has a rather fast $1/R^{6}$ falloff. When the evolved metric contains noise — e.g. caused by reflections on the grid refinement boundaries — the spurious terms decay less steeply and can become of comparable size to the stationary value.

Other invariants, which are simpler and easier to calculate, can be obtained from the fact that our spacetimes are axially symmetric. That implies the existence of an angular Killing vector $\eta^{\mu}=\frac{\partial}{\partial\varphi}^{\mu}$. When appropriately scaled, its norm is the circumferential radius

[TABLE]

where the second equality holds in the $y=0$ plane. Clearly in spherical symmetry we have $\bar{\rho}^{2}=R^{2}$ in the equatorial plane. Another related scalar we can construct is the norm of the gradient of $\bar{\rho}^{2}$

[TABLE]

Note that despite the notation, $\bar{\rho}^{\prime 2}$ can be negative — e.g. in Schwarzschild the $R=\mathrm{const.}$ surfaces are spacelike above the horizon and timelike below it.

From the above two scalars we can construct a dimensionless quantity

[TABLE]

In our simulations we can use $\chi$ as a measure of spacetime deformation in the equatorial plane. For the Schwarzschild solution $\chi=1-\frac{2M}{R}\sin^{2}\theta$, so in super-critical spacetimes the black hole horizon will settle down to the $\chi=0$ hypersurface. Additionally, in stationary regions we can form another mass estimate in the equatorial plane

[TABLE]

In future studies of the critical collapse of axisymmetric gravitational waves, the quantity $\chi$ may be of interest, since for a spacetime with discrete self-similarity such a dimensionless quantity should repeat with the same range of values on each echoing period [6].

3 Modified 1+log slicing

Though maximal slicing appears to be well-behaved for near-critical Brill waves, implementing it in a numerical code entails practical difficulties — solving elliptic equations at each time step is computationally demanding and errors involved immediately lead to constraint violations. However, since the lapse is a gauge function, we do not actually insist on it being exactly the maximal one. A ‘good-enough’ approximation that still avoids the problems of the 1+log slicing works just as well for our purposes.

In obtaining such an approximation, we start with the maximal slicing condition (5fghi) and take its time derivative

[TABLE]

which with some straightforward manipulations transforms into

[TABLE]

Now we replace the time derivative of the lapse with a new function $W$

[TABLE]

and use the evolution equations to replace the time derivatives of $\gamma_{ij}$ and $K_{ij}$ with purely spatial quantities known on one slice. We obtain

[TABLE]

where $\dot{K}_{ij}$ is a shorthand for the right-hand side of (5b). This is an elliptic equation for the new function $W$, with the same structure as the maximal slicing condition (5fghi), but a more complicated right-hand side.

Our starting point (5fghr) is nothing else than the demand that the acceleration of $K$ vanishes. So if we now solved (5fghu) exactly at each time step and used its solution $W_{e}$ to evolve the lapse according to (5fght), then — assuming that initially $K=0$ and the lapse is maximal — we would get precisely the maximal slicing (leaving aside the question of stability of such a system). Of course in a numerical simulation we will not have an exact solution and good approximations to it do not come cheaply. One can hope, however, that even a rough approximation to the solution of (5fghu) can be exploited to get closer to maximal slicing.

The basic idea is then as follows — at each time step we compute $W_{a}$, a (very approximate) solution to (5fghu). Then we use this quantity as an extra source term in the 1+log slicing

[TABLE]

where $\kappa$ is a function of coordinate time used to switch between slicing conditions — we set it equal to one when the extra term is to be active and smoothly send it to zero when we want to recover the original 1+log slicing.

The motivation for this extra term is as follows. In the weak-field regime, the Bona-Massó class of slicing conditions yields a hyperbolic $(K,\alpha)$ subsystem — a scalar sector corresponding to so-called lapse gauge waves. These exhibit the usual wave behavior and naturally drive weak-field configurations to $K=0$ as waves (assuming their finite ‘energy’) disperse outwards. The aim of our modification (5fghv) is to push the slicing towards $K=0$ even in a non-linear regime.

It is known that for PDE systems such as BSSN we are not free to modify even the gauge conditions arbitrarily — adding combinations of evolved functions and their derivatives may lead to loss of well-posedness of the system. The source term added in (5fghv) does not modify the principal part of the PDEs linearized around flat space, so it does not affect the hyperbolicity of the PDE system. This follows from the fact that (5fghu) has a trivial principal part — $\Delta W=0$. Also — as described later — in our approach we do not assume that the elliptic solver should provide a solution close to the continuum limit.

Further in this paper we show numerical evidence that slicing condition (5fghv) is not only stable for small-amplitude Brill waves, but also that for larger amplitudes the extra term acts as a driver that pushes the slicing closer to the maximal one, curing the problems of 1+log slicing. We call this condition the ‘quasi-maximal’ slicing.

4 Numerics

4.1 Implementation

Initial data

Due to axial symmetry, the equation (3) only needs to be solved in two dimensions. We use a pseudo-spectral method, writing the conformal factor as a series

[TABLE]

Using (5fghw) in (3) and demanding that the equation be satisfied exactly at $NM$ collocation points gives us $NM$ linear equations for the coefficients $C_{kl}$. We solve those using LU decomposition, which then allows us to reconstruct $\psi$ at arbitrary points through (5fghw).

Specifically, we choose polar coordinates $\left\{x^{0},x^{1}\right\}=\left\{r,\theta\right\}$ and basis functions [20]

[TABLE]

Here $L$ is a constant that determines the compactification scale and needs to be tuned to the problem being solved. In our simulations we always use empirically determined value $L=3$ (as usual, in the units of $\sigma$) for the initial data.

Functions $\mathrm{SB}_{k}(x)$ decay as $\frac{1}{x}$ towards infinity, which is the behavior we expect from $\psi$. Taken together this combination of basis functions automatically satisfies the symmetry and decay properties of the solution, so we need not impose any explicit boundary conditions. One property this basis set does not guarantee is regularity of the solution at origin — that would further require that the sub-series corresponding to each $B_{l}^{1}$ has a $2l$-th order root at $r=0$ (this is sometimes called the parity theorem [21]). However it turns out in practice that it is not necessary to impose these conditions explicitly and the solution is regular anyway.

We implemented the above procedure as a stand-alone library written in C/x86 assembler, which is then called from the evolution code to construct the initial data. The linear system is solved using LAPACK [22].

Evolution

Our evolution code is based on the Einstein Toolkit [23, 24], which bundles together the Cactus framework [25, 26] and other packages relevant to numerical relativity. The Carpet [27] code provides Berger-Oliger-style [28] fixed mesh refinement. The BSSN evolution equations are implemented by the McLachlan project [29, 30, 31].

To exploit the axial symmetry of our systems, we use Cartesian coordinates $\left\{x,y,z\right\}$ with the analytic Cartoon method, as described in detail in [12]. Its main idea is replacing the $y$ derivatives of arbitrary tensors with analytically derived combinations of $x$ and $z$ derivatives. Compared to the ‘classical’ Cartoon method [32] implemented in the Einstein Toolkit, which uses a thin layer of points in the $y$-direction filled by interpolation, the analytic method has significantly reduced memory requirements and is much faster. It should also be more accurate, though this was not our main motivation for using it. As our initial data is also symmetric with respect to reflection through the $z=0$ plane, we evolve just the $z\geq 0$ region.

The simulation domain is thus composed of a set of nested squares. Each one encloses an equidistant grid, with the step size $\Delta x$ halving per each nesting (refinement) level. We use the method of lines for time evolution, with 4th order Runge-Kutta as the time integrator. Spatial derivatives in the BSSN equations are approximated with 8th order finite differences — upwind in the advection terms, centered elsewhere. Kreiss-Oliger dissipation of the 9th order is applied to damp high-frequency noise. When (quasi-)maximal slicing is used, there is also a separate pseudo-spectral grid present, as described later.

Maximal slicing

As the maximal slicing condition (5fghi) is a linear elliptic PDE, same as the equation (3) we solve for initial data construction, we reuse the basic methods (and some of the code) from the initial data solver, with a number of changes.

The most obvious of those stem from the fact that the equation is not stand-alone, but a part of a system — it needs to be solved at each intermediate step of the time integrator and requires the evolved metric variables as input. Since the evolution happens on the equidistant finite difference grid, we need to interpolate the metric variables onto the pseudo-spectral grid. Additionally, due to mesh refinement, not all data will always be available at the required time level, which means we have to interpolate in time as well. We use fourth order Lagrange interpolation in space and linear in time.

Another change concerns the coordinates. When using polar coordinates we observed that the regularity conditions at origin stop being satisfied automatically during evolution and the lapse develops a discontinuity at $r=0$. We believe this happens because the data injected from the finite difference grid is not sufficiently smooth for the spectral solver. To avoid this issue we use Cartesian coordinates $\left\{x,z\right\}$ and $B_{k}^{0}=B_{k}^{1}=\mathrm{SB}_{2k}$. The scaling constant $L$ is chosen in such a way that the outermost collocation point remains causally disconnected from the outer boundary during evolution.

Since the equation is now solved many times, performance of the solver becomes important. For this reason we augment LU decomposition with the BiCGSTAB iterative method [33]. Since the pseudo-spectral matrix changes slowly with time, we run the LU decomposition once per $S$ solves and use it to compute the exact inverse. This inverse is then used as the preconditioner for the next $S$ BiCGSTAB iterative solves. Typically we take $S=1024$, then on average around 5 iterations are needed for the solution to converge to within $10^{-15}$ absolute error, which takes about $20\times$ less time to compute than the LU decomposition.

A significant weakness of this approach is that it involves solving a dense linear system of $N^{2}$ equations, where $N$ is the number of the basis functions in one dimension. Our memory requirements thus grow as $N^{4}$, so the maximum practically achievable number of collocation points is much lower than the number of points on the finite differencing grid.

With our chosen basis, the collocation points are closely clustered near the origin and become very sparse further out. So the way lack of resolution manifests in practice is that as the initial wave pulse travels away from the origin, it enters the area where the distance between collocation points is too large for the waves to be resolved. The higher frequencies then become invisible to the pseudospectral solver, strongly reducing the accuracy of the solution. Since we enforce $K=0$ for maximal slicing, the errors manifest as large violations of the ADM constraints. Similar issues have also been reported in [8].

We might be able to mitigate this issue somewhat with a better choice of the basis, or even resolve it fully through some sort of a multi-domain method. However that would require a significant amount of additional programming and, since maximal slicing is only an intermediate step of our work, we do not pursue such efforts any further.

Quasi-maximal slicing

The obvious approach to computing $W_{a}$ is reusing the method we used for maximal slicing — solve on each time step, using time interpolation to fill in missing data points. That, however, consumes a large amount of computational time and, as previously mentioned, we have considerable freedom to sacrifice accuracy for efficiency in solving (5fghu).

So we use another method, which requires far less resources and turns out to work well in practice. At its core is the pseudo-spectral elliptic solver we used for maximal slicing, but it is not run at each time step. It works as follows:

• •

We pick the finest refinement level $L$ that encloses all the collocation points used by the pseudo-spectral solver.

• •

At simulation time $t_{n}^{L}$, before we start the recursive Berger-Oliger evolution to time $t_{n+1}^{L}$, we have the final BSSN quantities for level $L$ and all the finer levels. After interpolating this data onto the pseudo-spectral grid, we execute the pseudo-spectral solver to solve (5fghu), obtaining a set of spectral coefficients $C_{n}$.

• •

During the following recursive time stepping we need the values of $W_{a}$ at times between $t_{n}^{L}$ and $t_{n+1}^{L}$, i.e. in the future from the solutions we already have. So we perform linear extrapolation from the last two solutions — $C_{n-1}$ and $C_{n}$ — to predict the coefficients at time $t_{n+1}^{L}$. We denote this prediction $C_{n+1}^{\prime}$.

• •

To evaluate $W_{a}$ at times in the interval $\left(t_{n}^{L},t_{n+1}^{L}\right)$, we linearly interpolate between the two most recent sets of predicted coefficients — $C_{n}^{\prime}$ and $C_{n+1}^{\prime}$. This ensures that $W_{a}$ is piecewise linear in time.

This process is sketched in Figure 1.

Finally, we discovered empirically that simply evaluating $W_{a}$ as a sum of the spectral series tends to introduce strong high-frequency components to the finite difference grid. It seems that those consequently propagate to $\alpha$ and from there to the other evolved functions and due to non-linearities get aliased as spurious low-frequency noise, manifesting as ADM constraint violations.

We deal with this issue by applying a low-pass filter to the basis functions. In our current multi-grid setup, the Nyquist frequency at a given grid point depends on the distance from origin $r$. For each basis function $B_{k}(x)$ we estimate its ‘base’ frequency from the location of its first zero $x_{k}^{0}$, which implies it is no longer accurately resolved at the radius $r_{k}^{0}$ where $\Delta x\approx x_{k}^{0}$. Therefore when evaluating $W_{a}$, we damp each basis function such that [34]

[TABLE]

$J=36$ ensures that for $r=r_{0}$ the exponential evaluates to the double-precision $\varepsilon$, while the empirically determined choice of $d=6$ makes sure the transition area is never smaller than the roughest step size (the results do not seem to depend on the precise value of $d$ very strongly). Altogether the filter ensures that $W_{a}$ remains non-zero only close to the origin and does not affect the grid further out, where it is not needed.

At late simulation times we usually switch from quasi-maximal back to 1+log slicing, since the added term no longer seems to be necessary for stability and for super-critical spacetimes the spectral solver can be ill-behaved in the presence of punctures. The switch is implemented by using a damping function similar to (5fghxy) for $\kappa$ in (5fghv), specifically

[TABLE]

where $t_{0}$ is the time at which the switch starts, determined empirically for a given simulation. The exact choice of parameters did not make a qualitative difference in our runs, $t_{0}$ is rougly the time when the apparent horizon forms for super-critical runs or the time when the largest peak in $\mathcal{K}$ forms and dissolves.

4.2 Code validation

Super-critical Brill waves in moving puncture gauge

We validate our moving puncture gauge setup by reproducing the slicing failure for $A=5$ Brill waves from [11]. We use 6 levels of mesh refinement, with $2048\times 2048$ grid points on each, the coarsest level using a spatial step size of $\Delta x=\frac{1}{4}$. The $\Gamma$-driver parameters are set to $\xi=1;\eta=11.4/\sigma$. Our results are fully consistent with the cited paper.

The slicing issues are illustrated in Fig. 2. The circumferential radius $\bar{\rho}$ should be a smooth spacetime invariant, so for a regular hypersurface $t=\mathrm{const.}$ the derivative of $\bar{\rho}$ along the unit normal $n^{\mu}$ should also be smooth. In the left panel we show that at times $t\lesssim 4.9$ — just before the simulations fail (independently of the spatial resolution) — this function develops a feature resembling a discontinuity. As can be seen in the inset, the gradient becomes higher with increasing number of grid points. When testing with a number of different resolutions from 384 to 2048 grid points, we found that the maximum of the gradient of the plotted function exhibits non-convergent behavior consistent with computing finite differences across a step function. We thus conjecture that the spatial slices themselves become non-smooth, which leads to simulation failure. The right panel then shows the ‘break’ in the same quantity in the $x-z$ plane.

Quasi-maximal slicing

First, we verify our solver for the equation (5fghu) does indeed converge to $\partial_{tt}K=0$. To that end, we run a series of simulations with weak waves ($A=0.1$) and no mesh refinement, where $W$ is evaluated at each full step of the time integrator (i.e. not using the extrapolation/interpolation procedure from Section 4.1). The simulations differ only in the time step $\Delta t$, which halves for each finer run — the coarsest run makes just one full time step, while the finest one makes 128. As the left panel of Figure 3 shows, the maximum norm of $K$ goes linearly to zero as the time step decreases.

Next, we run a series of simulations for stronger waves ($A=1$), with the grid setup as for the previously described 1+log runs. The reference run uses the 1+log slicing, and we compare it against several quasi-maximally sliced runs (now as described in Section 4.1) with varying CFL factor. As can be seen from the top panel of Figure 3, even when $W_{a}$ is a very crude approximation to the solution of $(\ref{eq:qms2})$, $K$ is significantly lower for quasi-maximal slicing than for 1+log, and further decreases with smaller time step.

The vanishing of $K$ as the time step goes down is slow and hard to quantify, so we find it more useful to regard quasi-maximal slicing merely as a procedure that gives us some slicing which is ‘close’ to the maximal one, but distinct from it. What we want to verify then is that when we hold the parameters of the quasi-maximal solver fixed and increase the resolution of the finite difference grid, the simulations converge to the same underlying spacetime, though possibly in different slicing. For that purpose, we look at spacetime invariants.

We run three simulations, again with $A=1$, and $\left\{256,385,512\right\}$ grid points per refinement level. For this data, the circumferential radius $\bar{\rho}$ in the equatorial plane $z=0$ remains a monotonous function of $x$ throughout the simulation, so it can be used as a radial coordinate (this is not the case for stronger Brill waves). For each simulation we calculate the value of $\bar{\rho}^{\prime 2}$ and $\mathcal{K}$ seen by the observer at $\bar{\rho}=1$. To work with purely invariant quantities, we also compute that observer’s proper time $\tau$ by integrating $\mathrm{d}\tau=-(g_{\mu\nu}\mathrm{d}x^{\mu}\mathrm{d}x^{\nu})^{\frac{1}{2}}$ and plot the aforementioned invariants as functions of $\tau$. As we can see from Figure 4, those values converge to each other with the fourth order, which is the order of the time integrator. We can thus conclude that quasi-maximally sliced simulations (with fixed pseudo-spectral solver order) do indeed converge to the same physical spacetime. As previously mentioned, it would be hard to state precisely the measure of convergence towards the maximal slicing.

5 Evolution of Brill waves with quasi-maximal slicing

Near-critical case: $A=5$

Our main result is demonstrating that quasi-maximal slicing cures the singularities arising with 1+log slicing for $A=5$ Brill waves. This near-critical initial data, with ADM mass $0.70$, was used as the failing case in [11]. Though quasi-maximal slicing resolves the singularities from the $\alpha$/$K$ sector, we have encountered additional issues with the shift. When using the $\Gamma$-driver, the shift seems to develop a shock profile around $t\approx 10$. For that reason, we use zero shift in the simulations described here. While that allows us to run the simulation further, it is well known that a non-zero shift vector is required in black-hole spacetimes to avoid ‘slice stretching’ — we are thus still unable to evolve this data indefinitely.

Around simulation time $t\approx 10.25$ we first find an apparent horizon with mass $M_{AH}=0.56$, which grows to $M_{AH}=0.58$ by $t=20$. The masses agree with values reported in [13]. By integrating the apparent horizon to the past as a null surface we construct the event horizon, seen in Figure 5. Its shape at early times turns out to be that of a disk with a non-smooth rim. This can be understood from the curious fact that if we shoot photons from the origin at $t=0$, those along the $z$-axis will escape to infinity, while those along the $x$-axis will end up below the horizon. The event horizon thus does not intersect the initial slice and appears around $t=1.7$, when the first photon in the $z$ direction fails to escape. Its worldline along the $x$-axis then has a spatial character until it connects with its generating null geodesic, as is sketched in the right panel of Figure 5. This means the event horizon can briefly exist as a torus, which we observed for a run with quasi-maximal slicing and $A=4.8$.

There are various sources of uncertainty in determination of the event horizon. The most important one is the error of the simulation, which can be judged by relative variation of late-time MOTS areas between various simulations and appears to be of order $10^{-3}$. Of similar order is the error with which we can trace null geodesics forming the event horizon (this is more accurate along the $x$ and $z$ axes, where the data available to the tracing code is much better resolved in time). Then there is the difference between the event horizon and MOTS and the uncertainty of future evolution of the spacetime, but due to unstable nature of future-oriented null geodesics near the event horizon, we find that at $t=10$ a spatial variation of $\Delta x\approx 10^{-5}$ will — during the runtime of the simulation — either direct the null ray inside the black hole or put it on an escape trajectory.

To demonstrate that the quasi-maximal run actually does cover a larger chunk of spacetime than the crashing 1+log run (rather than merely collapsing the lapse in the critical area), we compare the two slicings in double-null coordinates $\left\{U,V\right\}$. Those are constructed by integrating null geodesics along the $x$-axis in, respectively, the positive and negative $x$ directions. The values of $U$ and $V$ in the initial slice $t=0$ are both set equal to half the proper distance from the origin. Figure 6 clearly shows that quasi-maximal slicing advances farther towards the physical singularity in the central area.

Super-critical case: $A=5.5$

The ADM mass of this data is $0.84$. Unlike the previous case, we are able to use the $\Gamma$-driver without any pathologies, which allows the simulations to run much longer. At $t\approx 7.5$ we find an apparent horizon with mass $M_{AH}=0.63$, which quickly increases to its final value of $0.73$.

At later times, the spacetime should evolve into a central region that rapidly relaxes into a Schwarzschild black hole, plus an outward-travelling wave packet. We illustrate this fact in Figure 7 by plotting mass estimates (5fghm) and (5fghq). These, as mentioned, have the meaning of mass only in the static end-state, so we can see them oscillate wildly as the waves radiate away. At late simulation times they closely approach the value of $M_{AH}$, which shows that we really obtain a Schwarzschild black hole with expected properties.

The $\Gamma$-driver behaves in a manner similar to that known for the evolutions of Schwarzschild starting with ‘wormhole’ slices, effectively excising the central region from the numerical grid. At late times the numerical slices terminate at $\bar{\rho}\approx 1.31$.

Other near-critical amplitudes

We run a number of additional simulations with both sub- and super-critical initial data until, respectively, the waves disperse in the central region or an apparent horizon forms. We observe no slicing singularities in any of them.

For analyzing the spacetime dynamics we rely on the invariants introduced in Section 2.3. Figure 8 shows values of the scalar $\chi$, defined by (5fghp), for long-term evolutions of sub- and super-critical waves. As previously mentioned, $\chi$ is negative when the gradient of the circumferential radius is timelike. In Schwarzschild this is true inside the event horizon, with $\nabla^{\mu}R$ pointing to the future in the ‘white hole’ region and to the past in the ‘black hole’ region. It is interesting that for Brill waves, the orientation is to the future in the first negative-$\chi$ region (i.e. ‘the wrong way’) and then switches direction, possibly multiple times, which shows how wildly dynamic those spacetimes are.

A drawback of $\chi$ is that it only has a clear mapping to Schwarzschild in the equatorial plane, since physically meaningful values of $\theta$ are not available. One way around that is the scalar $\frac{1-\chi}{\bar{\rho}^{2}}$, which is equal to $\frac{2M}{R^{3}}$ in Schwarzschild. It is, however, no longer dimensionless. This quantity is shown alongside $\mathcal{K}$ in Figure 9 for a sub-critical run at the time when those scalars attain its highest values. Interestingly this happens on the $z$ axis some distance away from the origin (cf. the results of [13]).

To compare the spacetime geometries with differing amplitude parameter $A$ side-by-side, we mimic the usual compactification of the Schwarzschild manifold by a transformation that maps the positive part of the $x$-axis through $X=\mathrm{arctan}\left[\frac{3}{2}\mathrm{sinh}\left(\frac{x}{6}\right)\right]$ into the interval $(0,\frac{\pi}{2})$ and then maps the null geodesics in the $x-t$ plane into straight $45$-degree lines. This produces the conformal diagrams in Figure 10, which illustrate propagation of waves towards future null infinity.

6 Summary

We run numerical simulations of sub- and super-critical Brill wave initial data. We manage to verify previously described slicing problems for strong data with a different code and illustrate them with figures showing that the $t=\mathrm{const.}$ slices develop a cusp-like shape. These problems turn out to be removable by amending the 1+log slicing with a source term derived from the maximal slicing condition.

This new slicing condition, which we call the ‘quasi-maximal’ slicing, improves the regularity of evolved spatial slices by pushing them closer towards the maximal ones. We describe long-term evolution of the spacetimes arising from near-critical waves and construct their conformal diagrams. Of special note is the non-smooth shape of the event horizon for weaker data. This data is not obtainable with the 1+log slicing because of the aforementioned pathologies.

As mentioned in the introduction, Hilditch et al. are currently working on the same problem with a new pseudospectral code using a generalized harmonic formulation [12, 13], so a brief comparison is in order. Since our code is based on finite differences and the moving puncture approach, it inherits their advantages, such as simplicity, robustness and the ability to run long-term simulations containing black holes without the need for excision. The price to pay is mainly lower efficiency, i.e. higher resolution required for the same accuracy. Of course, the fact that both these approaches produce the same results (as confirmed in private communication), is greatly encouraging.

The code used for running the simulations described in this paper is available for download through the git protocol at git://git.khirnov.net/qms_source_2017.

We would like to thank David Hilditch for very productive discussions. This work is supported by the Charles University in Prague, project GA UK No 2000314, GA UK No 1176217 and SVV-260211. T.L. acknowledges the support from the Czech Science Foundation grant No 14-37086G (A. Einstein Center). Computational resources were provided by the CESNET LM2015042 and the CERIT Scientific Cloud LM2015085, provided under the programme ‘Projects of Large Research, Development, and Innovations Infrastructures’. We would also like to thank the developers of the Einstein Toolkit for making open-source numerical relativity possible.