# Evolutions of centered Brill waves with a pseudospectral method

###### Abstract

The pseudospectral code bamps is used to evolve axisymmetric gravitational waves. We consider a one-parameter family of Brill wave initial data, taking the seed function and strength parameter of Holz et. al. A careful comparison is made to earlier work. Our results are mostly in agreement with the literature, but we do find that some amplitudes reported elsewhere as subcritical evolve to form apparent horizons. Related to this point we find that by altering the slicing condition, the position of the peak of the Kretschmann scalar in these supercritical data can be controlled so that it appears away from the symmetry axis before the method fails, demonstrating that such behavior is at least partially a coordinate effect. We are able to tune the strength parameter to an interval of range around the onset of apparent horizon formation. In this regime we find that large spikes appear in the Kretschmann scalar on the symmetry axis but away from the origin. From the supercritical side disjoint apparent horizons form around these spikes. On the subcritical side, down to this range, evidence of power-law scaling of the Kretschmann scalar is not conclusive, but the data can be fitted as a power-law with periodic wiggle.

###### pacs:

95.30.Sf, 04.25.D-## I Introduction

Motivated partially by the findings of HilBauWey13 (), in particular by our difficulties in evolving gravitational wave data close to the critical threshold of black hole formation with the moving puncture gauge, we turned to an alternative formulation and a more accurate numerical method. We implemented the generalized harmonic formulation in a pseudospectral code, bamps, which was recently described in detail in HilWeyBru15 (); BugDieBer15 (); HilHarBug16 (). Presently we use this tool to evolve Brill wave initial data Bri59 () in the form most often treated numerically. Primarily we choose such data for ease of comparison with the literature, but additionally since it is axisymmetric it allows us to run the code most efficiently. Ultimately we hope to obtain a proper understanding of, and a robust numerical method for gravitational waves close to the threshold of black hole formation. This study, like HilHarBug16 (), is another step in that direction.

The key results in the literature on critical collapse of gravitational waves are those of Abrahams and Evans AbrEva93 (); AbrEva94 (), who considered one-parameter families of Teukolsky wave initial data in axisymmetry, and found that the resulting black hole mass scales as a power law in a neighborhood of the critical threshold. They also found more tentative evidence for ‘echoing’, or periodicity in spacetime scale, of the solution. Primarily because of its simplicity, most subsequent studies have focused on Brill wave initial data. In evolutions of Brill waves Sorkin Sor10 () found evidence for another critical solution in which the peak of the curvature appears on concentric rings around the symmetry axis. This study employed the generalized harmonic formulation in axisymmetry, and so is the natural starting point for us. For a detailed discussion of critical phenomena in gravitational collapse see GunGar07 ().

## Ii Setup

The bamps code uses a pseudospectral method to evolve a first order formulation LinSchKid05 () of the generalized harmonic gauge formulation, with

(1) |

in the standard notation. We parametrize the free scalars by and , with and constant. We employ radiation controlling constraint preserving boundary conditions like those described in Rin06a (); RinLinSch07 (), imposed via the Bjørhus method Bjo95 (), but modified to minimize reflections caused by the use of constraint damping, which can otherwise cause the code to crash with our gauge conditions. The numerical method is similar to that of SpEC spec (), employing many subpatches across which data is communicated using a penalty approach. For our grids we presently take either cubed-sphere or cubed shells RonIacPao96 (); Tho04 (). We discretize in space using Chebyschev polynomials, filtering the highest order coefficients SziLinSch09 (). Although bamps is fully 3d, to evolve efficiently in axisymmetry we have implemented the analytic-cartoon method Pre05 (); AlcBraBru99a (), in which all angular spatial derivatives are evaluated using the Killing vector. The resulting regularity conditions are imposed on axis and at the origin. The code is parallelized with MPI at the subpatch level, with one or more subpatches per core, resulting in excellent strong-scaling for up to several thousand cores. Presently we evolve Brill wave initial data Bri59 (); Epp77 (), in which the spatial metric takes the form,

(2) |

and the extrinsic curvature vanishes. We choose always the seed function of HolMilWak93 (),

(3) |

with , and . We call this data a centered geometrically prolate Brill wave. To build initial data we use the solver of DieBru13 (). Our main tool for classifying spacetimes as supercritical is the axisymmetric apparent horizon finder AHloc, which we run in postprocessing. A detailed description of the setup can be found in HilWeyBru15 (); Wey14 ().

## Iii Centered Geometrically prolate Brill waves

In this section we present our evolutions of centered geometrically prolate Brill waves. The first numerical evolutions of Brill waves that we are aware of were presented in Epp79 (). Since then Brill wave initial data have been considered multiple times in the numerical relativity literature AlcAllBru99 (); GarDun00 (); Rin06 (); San06 (); HilBauWey13 (). Therefore to ensure that bamps is performing properly we start with a detailed comparison of the evolutions performed with the seed function (3) away from criticality. We then compare with the results of Sor10 () before going towards criticality.

### iii.1 Review and comparison with earlier work

#### Alcubierre et. al.:

In AlcAllBru99 () Brill waves were evolved numerically for the first time in 3d. The BSSNOK formulation was used in combination with maximal slicing and vanishing shift. Using the given data (3) with , , which we consider throughout, it was found that the critical point lies between and , and furthermore that this finding could be refined to , although the data for this latter claim were never presented. Supercriticality was diagnosed by finding an apparent horizon, which occurred for the data at . In GunGar06 () it was shown that BSSNOK combined with this gauge choice results in an ill-posed PDE system, meaning that this approach should be either abandoned or modified if we are to achieve accurate results that converge to the continuum solution as resolution is increased.

#### Garfinkle and Duncan:

In GarDun00 () it was found, evolving Brill wave initial data with as in (3), again taking , , that the critical amplitude lies between and . The data was classified either by evolving until the spacetime was close to flat and subsequent collapse seemed implausible, or by explicitly finding an apparent horizon. The formulation employed was explicitly axisymmetric, and consisted of a mixed elliptic-hyperbolic system with maximal slicing , well-posedness of which, to the best of our knowledge, has not been studied. We agree with the findings of both AlcAllBru99 () and GarDun00 (). Because our method employs a different gauge however, it is difficult to make a side-by-side comparison beyond classifying the spacetimes as sub or supercritical. The effect of changing the shape of the initial data parameters and was also studied in GarDun00 (), but we have not yet followed up on this.

#### Rinne:

In the PhD thesis Rin06 () evolutions of centered geometrically prolate Brill waves were presented with a free-evolution and a partially constrained scheme, both in explicit axisymmetry. We focus on the free-evolution scheme, the Z formulation, since that is where we are able to make the clearest comparisons. In that case harmonic slicing was taken with vanishing shift. This choice is convenient for our comparison because although we can not choose the same shift condition, harmonic lapse is a pure slicing condition, which means that we should obtain the same foliation of the same spacetime (starting from the same initial lapse) albeit with different spatial coordinates. Since there is a preferred observer, namely that at the origin, we can compare quantities explicitly there. Fortunately the work Rin06 () contains several plots along this worldline. In the upper left panel of Fig. 1 we plot the Kretschmann scalar

(4) |

at the origin as a function of time for centered Brill wave data, which should be compared with Fig. 9.2 of Rin06 (). In this test we evolved with harmonic slicing and the damped harmonic shift and otherwise our standard setup. The agreement, at least by eye, is extremely good. Taking it was found that with sufficient resolution a sharp feature in the gradient of the lapse could be resolved. It was found that the data was, in agreement with GarDun00 (), subcritical. We see the same result. In the upper right panel of Fig. 1 we show at . This is the time at which Rinne finds the largest peak in this quantity. A similar plot is Fig. in Rin06 (), which can not be directly compared because of the differing spatial coordinates, although the qualitative agreement is very clear. We find that the largest peak appears at around , but the magnitudes in differ by less than across these times. In the lower two panels of Fig. 1, to be compared with Fig. of Rin06 (), the left panel shows the logarithm of the lapse at the origin as a function of coordinate time for amplitudes . In the right hand panel we plot the value of the Kretschmann scalar at the origin. At lower resolutions we find that the Kretschmann scalar exhibits high-frequency oscillatory behavior, but that these wiggles converge away rapidly. The most challenging data evolved in Rinne’s experiments was the wave, for which sub and supercriticality was not discerned using the free-evolution algorithm, partially because the data was still oscillatory at the time the method failed at around . All of the different resolutions we tried with in this suite of tests crashed at coordinate time . We ran our apparent horizon finder on this data, the result of which is plotted in Fig. 2. We find the apparent horizon first at and thereafter until the evolution fails. Rinne correctly concluded that the critical amplitude lies below , although no apparent horizon could be found in his data. Instead his classification was made by observing that the Kretschmann scalar was blowing-up. This diagnostic is flawed because as one approaches the critical point we expect to generate arbitrarily large curvature scalars even in subcritical data. In the absence of an apparent horizon or event horizon however, other diagnostics may be similarly flawed, and the Kretschmann scalar is at least a spacetime scalar, so one might prefer it as a diagnostic to ‘collapse of the lapse’ Alc08 () which is clearly gauge dependent. Running our apparent horizon finder on the data we find an apparent horizon at and later.

#### Previous studies with BAM:

In San06 () the BAM finite differencing code HilBerThi12 (); ThiBerBru11 (); BruGonHan06 () was used to evolve centered geometrically prolate Brill waves with the BSSNOK formulation combined with several different gauge conditions, maximal slicing, harmonic slicing, and the moving puncture gauge condition. The main complications were reported to be constraint violation, which was likely caused by lack of resolution, and which did vary significantly from one gauge to another. To understand how much further, if at all, standard modern numerical relativity methods can go beyond those previously discussed, recently in HilBauWey13 () the BAM code was once again used, this time alongside the code of BauMonCor12 () with the focus purely on the moving puncture gauge. Starting with (that is, weak) data, it was found that the lapse initially decreased, but rapidly returned back to unity; we find qualitatively the same behavior despite the different gauge used in bamps, although the lapse function decreases by a smaller amount in the new data. The Kretschmann scalar decreases from a maximum of about 216 at to zero and reaches a second maximum at . The maximum value of the Kretschmann scalar in the domain immediately decreases and never grows beyond the initial value as the wave propagates away. In the upper left panel of Fig. 1 we also plot the central value of the Kretschmann scalar obtained in this BAM experiment. Since the time coordinates used differ, the horizontal axes would be different, but the values of the maxima should be the same. The BAM value is about which agrees extremely well with the bamps experiment. The Hamiltonian constraint violation in this particular BAM run is of the order at the time of the second local in time maximum, whereas the roughly analogous constraint inside bamps is less than . This is not a fair comparison, because we are not considering computational cost whatsoever, but does indicate that the bamps data is superior in this case. Therefore one expects that as more resolution is added to the BAM grid the result would converge to the bamps result. Taking data that is stronger, for example , the method of HilBauWey13 () failed as an incoming pulse in the lapse became evermore sharp, resulting in what seemed to be a coordinate singularity. This would be acceptable if an apparent horizon could be found before the code crashed, but this was not the case. Going to higher amplitudes still, similar failures occurred, and the conclusion was drawn that moving puncture coordinates were not suitable for managing this initial data. We have already seen in our comparison with Rin06 () that using bamps the data can be classified supercritical, and we did not see any sign of a coordinate singularity before an apparent horizon was discovered despite the two lapse functions appearing qualitatively similar at the beginning.

### iii.2 Comparison with Sorkin

Having collected a bank of evidence that our numerical results are correct we now compare with Sor10 (). We demonstrate firstly that we can obtain qualitatively the same type of behavior described therein, namely that the peak of the Kretschmann scalar appears away from the symmetry axis. This we achieve however just by evolving supercritical Brill wave data and changing the gauge source function we use to evolve it. We find that the position of the peak of the Kretschmann scalar can be controlled by the choice of gauge source function. This can be understood geometrically. Secondly, by locating an apparent horizon in the time development, we demonstrate explicitly that the amplitude , evolved and classified in Sor10 () as subcritical, is in fact supercritical. We looked at evolutions of this data at several resolutions, both of the numerical spacetime and the apparent horizon computed on top of the data, and find that the outcome is robust.

#### Simulation setup:

The simulations of this subsection have been carried out on a cubed-ball grid, as described in HilWeyBru15 (); Wey14 (), with the following setup. The inner cube extends from to and is divided into subpatches with gridpoints in each dimension. For the transition shell from to we use only one shell with points in the radial direction. From here we go to the outer boundary at using outer shells with also radial collocation points. These tests were carried out in 3d using the octant symmetry mode of bamps.

#### Position of the peak curvature:

Consider a collapse spacetime with the standard causal structure. Different foliations of the spacetime, in which the time coordinate tends to tick more or less slowly in a region of high curvature depending on some singularity avoidance parameter will have different profiles in a spacetime diagram. If we are given a patch of this spacetime up to a finite time coordinate, as in a numerical relativity simulation, the specific observer that encounters the largest curvature before the code crashes depends, among other things, on the singularity avoidance parameter. In our case such a parameter is given by . We performed evolutions of a centered Brill wave with . This amplitude is shown to be supercritical in the next paragraph. We evolved with fixed and one of or . In the left plot of Fig. 3 we show the lapse in the evolutions at coordinate times , which demonstrate the effect of the singularity avoidance parameter . In the pure harmonic slicing case we have the strongest singularity avoidance, and find that the peak of the Kretschmann scalar appears at , where it simply grows until the numerics fail. Increasing the parameter to we find that the peak of the Kretschmann appears at a coordinate radius of and respectively, before the numerics fail.

#### Apparent horizon formation:

We evolved the same centered data with different resolutions, fixing the gauge parameters and , so the largest choice of from above. We find rapid convergence of the constraints. In particular we used our standard cubed-sphere setup with , and points per subpatch, and find that, for example, the maximum of the component of the Harmonic constraints along the x-axis are approximately and respectively at . We then searched for apparent horizons using the method described in HilWeyBru15 (); Wey14 (). We first find an apparent horizon at around . On a fixed numerical spacetime data set we find perfect fourth order convergence in the apparent horizon data consistent with the Runge-Kutta method employed. Comparing the apparent horizons discovered on the different data we find perfect qualitative agreement. Furthermore we see behavior consistent with rapid convergence when evaluating the differences between the discovered horizons. The apparent horizons of different time slices are plotted in Fig. 3, where the coordinate expansion of the horizon can clearly be seen. For comparison we again evolved the same initial data inside the BAM finite differencing code, but were unable with our current setup to find apparent horizons from that data. One issue is that the constraint violation is many orders of magnitude greater in the finite differencing code. Qualitatively however we find good agreement in the evolution, at least initially.

#### Summary:

We are unable to reproduce the results of Sor10 (), despite using, to the best of our knowledge, identical initial data and gauge. In fact our results appear to contradict the earlier study. The reason for the disagreement is presently not clear. It is possible that an apparent horizon search on Sorkin’s older finite differencing data was too challenging because of numerical error, or perhaps even that the numerical dissipation was sufficient to let these strong data spuriously settle down to flat-space, although the latter does not seem likely. The power-law scaling obtained in Sor10 () in the rapid oscillations of the curvature is nevertheless interesting, and would be good to properly understand in the future. Given our earlier code validation HilWeyBru15 () and the literature comparison in section III.1 however, we have no reason to doubt the bamps results, so for now we move on.

### iii.3 Discussion

#### Critical collapse of the scalar-field:

Before moving to harder experiments, consider the case of the minimally coupled massless scalar field. Working in spherical symmetry, evolving families of initial data with strength parameter , Choptuik found Cho93 () compelling numerical evidence for the existence of a critical solution at serving as the boundary between dispersion and collapse. The appearance of critical phenomena with in a neighborhood of was neatly explained by the conjecture that the critical solution is an attractor of co-dimension one in phase space. In other words in a neighborhood of the critical solution there should to be just one growing mode. Working in perturbation theory around the critical solution, but crucially allowing for aspherical mode perturbations, Gundlach and Martín-Gárcía found strong numerical evidence for this conjecture GarGun98 (). They found that the most slowly decaying mode came with an eigenvalue of associated with a spherical harmonic. For this matter model the single unstable mode has eigenvalue . Because of the exponential decay of the former it was conjectured that the qualitative picture obtained in spherical symmetry using the full Einstein equations would not be altered for generic initial data close to criticality. This picture, roughly speaking, says that in a neighborhood of critical collapse the fields should strongly interact in a confined region for a finite, but ever longer time as the critical solution is approached, and ultimately either collapse or disperse. Interestingly for the current study, Choptuik and collaborators ChoHirLie03 () then studied axisymmetric configurations and found evidence of a second growing mode associated with a spherical harmonic, causing the strong field solution to bifurcate into two strongly interacting regions, in apparent contradiction with the earlier perturbative result.

#### Issues with the apparent horizon as a diagnostic:

The apparent horizon is used for classification of the spacetime as sub or supercritical. In earlier work the apparent horizon mass has also been used to indicate power-law scaling in the critical regime. Therefore it is worth noting explicitly the weaknesses of this approach. Firstly, the appearance of an apparent horizon guarantees the existence of an event horizon only in strongly asymptotically predictable spacetimes, which arise from generic asymptotically flat initial data only if the weak-cosmic censorship conjecture holds HawEll73 (); Wal84 (). This fact makes our approach blind to violations of weak-cosmic censorship. It is also possible that no apparent horizon appears in our particular foliation, even if there is a black hole region WalIye91 (). Secondly the foliation dependence of the apparent horizon makes it difficult to trust the black hole masses obtained from the horizon area, particularly away from spherical symmetry. We may try to diagnose power-law scaling by looking at horizon mass, but the foliation could be such that the apparent horizon always forms with large area, and we have to choose when to evaluate the mass for each member of the one-parameter family. In fact in preliminary experiments we found that when approaching the critical regime, behavior resembling power-law scaling in the initial horizon masses could appear with a particular choice of gauge source function, but then completely disappear once we altered the choice to avoid coordinate singularities. These difficulties make us strongly prefer to study the subcritical approach to black hole formation, because there we can, as in Sor10 (), unambiguously consider the spacetime maximum of the Kretschmann scalar, provided that we are able to evolve the spacetime long enough to be confident that it is really subcritical.

### iii.4 Towards the critical regime

Sweep | Changes | |||||

/ | ||||||

X | X | X | ||||

/ | ||||||

X | X | X | ||||

/ | ||||||

X | X | X | ||||

Bisect | ||||||

1 | X | X | X | |||

2 | / | |||||

3 | X | X | X | |||

4 | X | X | X | |||

5 | / | |||||

6 | / |

Grid | |||||||
---|---|---|---|---|---|---|---|

G0 | |||||||

G1 | |||||||

G2 | |||||||

G3 |

#### Search strategy:

Our previous tests demonstrate that apparent horizon formation first occurs between and . We therefore searched for a critical amplitude in this range. Running exclusively in cartoon mode, we started with that bracketing and within it sampled at equally spaced values, thus obtaining a new bracket times smaller. We did this in three times, with for . At each level, once the critical amplitude is bracketed we increase resolution on the bracketing amplitudes to check convergence and be sure that we truly have bracketed . After the third such sweep we went to a straightforward bisection search, increasing resolution, adjusting grids and gauge sources as seemed appropriate. Our current best bound for the critical point is that , an interval of width . Afterwards we performed additional runs away from the critical point to help understand the behavior of the Kretschmann scalar and the initial apparent horizon masses as a function of . We started with the gauge source parameters

(5) |

and take our defaults for the GHG formulation settled on in HilWeyBru15 (); Wey14 (). The parameters of our base grid (G0 in Tab. 1) result in a total of subpatches in cartoon mode. Spreading this over the maximum number of cores possible on SuperMUC (i.e. one subpatch per core), the code computes at around time units per hour at the lowest resolution. Modifications to the gauge and grid are summarized in Tab. 1, which gives the results both of the bracketing for each sweep and in the bisection search. Grids G1, G2 and G3 consist of and subpatches, respectively. The ADM masses of each initial data set are also given in the table, although one should be careful to remember that we impose boundary conditions at a finite coordinate radius, which could affect the dynamics of the evolution, and furthermore makes the interpretation of non-trivial.

#### Termination of search:

It may be possible to push to a better bound by brute force with the current method, but we stopped our bisection search at the range as we prefer to conserve resources to attack alternative initial data. A main issue preventing us from going further economically is the lack of mesh-refinement inside bamps. The implementation of a true pseudospectral adaptive mesh-refinement algorithm is a major undertaking due to its complexity, and because efficient parallelization then becomes challenging. Another issue is that our initial data solver can only solve the Hamiltonian constraint down to around the level, which is presumably caused by the use of irregular coordinates, the Chebyschev-Fourier-Fourier discretization and simply machine precision Wey14 (), but this level of error is certainly not the leading order in our present simulations. It is curious that in studies with various matter models it has been possible to fine-tune to much higher accuracy, even in cases where the basic accuracy of the numerical method is much lower than in ours. This can often be achieved by keeping the numerical resolution and everything else, except the amplitude fixed. However, the fine-tuning error in for a given resolution can be much smaller than the drift towards convergence when increasing the resolution. As we explain below, we believe that in the present simulations the main issue preventing us from going further in the fine-tuning of is related to gauge choice.

#### Description of dynamics:

The basic dynamics from each of the initial data are initially rather similar. At first a pulse in the Kretschmann scalar propagates out from the origin predominantly in the direction. The pulse then propagates more slowly, eventually turning around and traveling towards the origin. As it propagates in, the pulse is smeared out parallel to the -axis. As the pulse hits the axis, there is a rapid growth resulting in a maximum at some . In the first sweep, for data an apparent horizon is found around or just after the time of this growth. In the run, this peak in the Kretschmann occurs at , with a value of . The feature then starts to propagate away, predominantly in the -direction, and no apparent horizon is found until later. Instead the evolution continues until a second large peak appears in the Kretschmann scalar around on the symmetry axis, as the wave content leftover from the first big peak again crashes onto the axis. An apparent horizon is found shortly afterwards and consistently until the evolution fails at , as the Kretschmann scalar starts to grow ever more rapidly around these peaks. In the second sweep we switch the slicing condition to to avoid spikes appearing in the lapse. The largest amplitude of this sweep, , is subcritical. The peak of the Kretschmann in this run is around and appears at ; a movie of the dynamics can be found online HilWebsite (). As we go closer to the critical point we see the appearance of yet-more spikes, which then propagate up and down the axis. This behavior is discussed in more detail in the following paragraphs.

#### Disjoint apparent horizons:

Starting with the third sweep, , we find in supercritical data two disjoint apparent horizons, centered roughly around the position of the second large feature growing on the axis at . An example is shown in Fig. 4. In other words, close to the critical point these initial data produce what seem to be axisymmetric binary black hole spacetimes! Obviously in the case of gravitational waves spherical decay is impossible, so this bifurcation is perhaps the generic near-critical behavior. Evidence for this could be sought by evolving different families of data. We expect that the reflection symmetry about the x-axis plays a role here in the outcome however. For generic axisymmetric supercritical data, with one loosely defined strong-field region, one might expect that the bifurcation happens but that apparent horizon formation appears on only one side. This result means that with our current setup we are not able to evolve to a final end-state, as bamps does not have a moving-excision setup or the control mechanism of SpEC SchGieHem14 (); HemSchKid13 (). Within the lifetime of our simulations we do not find a common horizon surrounding the disjoint MOTS. But the lifetime of the simulations after horizon formation is short, so this to be anticipated. To be sure that these spacetimes really do contain two black holes it will be necessary to search for an event horizon, but the short lifetime prohibits this also. Nevertheless the fact that as we get closer to the critical solution, the initial apparent horizon size gets smaller, whilst the coordinate distance between the horizons remains roughly constant hints that the spacetimes do contain two black holes. In the context of critical collapse we are predominantly interested in the strong-field region near to the critical threshold, so we continue with the search, focusing primarily on the subcritical regime. The detailed study of supercritical data is left for future work.

#### Scaling of the Kretschmann scalar:

According to GarDun98 (), if critical phenomena are present during gravitational collapse, then one should see power-law scaling of curvature invariants in the subcritical regime . Since we are working in vacuum any scalar built from the Ricci curvature is unavailable so, as in Sor10 (), we focus on the Kretschmann scalar. In the right hand panel of Fig. 4 we plot the maximum value of the Kretschmann scalar in the spacetimes as a function of in a log-log plot. There is a plateau in the maximum before it starts to increase rapidly in . The rapid increase occurs as a later implosion of the wave onto the axis starts to dominate over the previous implosion. The resulting curve can be fit as,

(6) |

with and of period around . Over one period the maximum in the Kretschmann scalar increases by a factor of , with . The first of these numbers agrees with that obtained in AbrEva93 (), whilst our value of does not. It is not yet clear how seriously these numbers should be taken because so far we observe only one full period in the Kretschmann scalar. Therefore we postpone the assignment of error bars and for now simply advise caution against over-interpretation of the finding.

#### Comparison with the scalar-field:

Our results are reminiscent of the bifurcation of the scalar field ChoHirLie03 () discussed earlier in section III.3. It is possible that there is a direct relationship; our results indicate that in vacuum axisymmetry, near the critical solution, decay proceeds by the aforementioned bifurcation. On the other hand, we know empirically that in spherical symmetry dispersion of the scalar field is determined by a single unstable spherical mode. Take a spacetime with a single strong-field region with scalar field and gravitational wave content. Imagine fixing, in some sense, the ratio of gravitational wave and scalar field content, and heading towards the threshold of black hole formation. By continuity, we expect the critical solution to interpolate between the two scenarios of bifurcation (driven by the gravitational wave content) and spherical decay (driven by the scalar-field content) as the ratio of the two is adjusted. This idea is compatible with ChoHirLie03 (), in which gravitational wave content is added to the initial data by placing axisymmetric perturbations on the metric. It is also compatible with the perturbative results of GarGun98 (). By construction the Choptuik solution is absent of gravitational waves, so a linear analysis could not spot a complicated nonlinear admixture of the decay mechanisms. Alternatively, it may not matter whether the strong-field is formed by a gravitational wave or other source. Comparing spherically symmetric with axisymmetric evolutions, the main difference is that it is possible to form multiple centers of collapse in axisymmetry. With sufficiently large asphericity this could be generic.

#### Wishlist for future work:

Evidence for the above suggestion could be sought in several obvious ways. First, evolution of different families of axisymmetric vacuum data must be performed to see whether or not the bifurcation behavior really is generic. Next, it would be good to compute accurately the value of in vacuum (assuming that the tentative behavior persists) and to compare with the value obtained by ChoHirLie03 () as the critical solution is deformed by gravitational wave content. At first glance our value appears consistent with that of ChoHirLie03 (), but at this stage nothing is certain. Another possibility is to work in second order perturbation theory about the Choptuik critical solution and to look for evidence of the bifurcation behavior. Finally, more results for the critical axisymmetric scalar field are also highly desirable.

#### Spikes in the Kretschmann Scalar and code failure:

For our final test, which we can not yet classify, close to the critical point we evolved initial data with amplitude . If this experiment were successful it would correspond to the next bisection step. We find that there are a sequence of large spikes on the symmetry axis as the gravitational wave implodes, then propagates up and down the symmetry axis before imploding once more. Each of the large spikes is finer and therefore requires more resolution for accuracy. It is tempting to label the sequence of strong oscillations ‘echoes’, but again, perhaps because of a suboptimal gauge, we can not quantify this claim and therefore resist. Fig. 5 shows the run-up to and the evolution of the final spike before the code crashes. In practice the numerics fail not as such a spike forms but rather as it dissipates away. As this happens we see along the -axis that a sharp feature suddenly forms in the metric component and causes the code to crash. The difficulty is in classifying the spacetime rather than the code crash per se, as do not find an apparent horizon in the data before the crash. The cause of the feature is unclear, but possible candidates are simple numerical error, the formation of an apparent horizon that the present method is unable to unveil, the formation of a coordinate singularity, or even the seemingly unlikely formation of a naked singularity. Increasing resolution very substantially hardly affects the appearance of the feature or crash-time of the code, so in this case it is doubtful that even mesh-refinement could address the problem directly. The main suspect is therefore the formation of a coordinate singularity. This view is further enforced by the fact that coordinate problems, although of a different specific form, also occurred in simulations with the qualitatively similar 1+log slicing using the BSSNOK formulation San06 (); HilBauWey13 (). To investigate this we have evolved with different gauge source parameters. Informed by earlier experience we increased . This however has the unfortunate side-effect of allowing the strong-field region to bleed out from the central cube into the transition shell where we have lower resolution, and so results in other problems. Going to slightly lower wave amplitudes, the sudden spike in persists even with a large value of , albeit at a later coordinate time. Ultimately a radical change of coordinates may be needed. Addressing the problem with an improved continuum formulation and numerical method is a priority.

## Iv Conclusions

We have continued our study of gravitational waves in the regime separating dispersion from black hole formation. To maximize overlap with earlier results we focused exclusively on Brill waves with the seed function (3), and evolved only prolate (), geometrically prolate () centered () data. Our main findings are first that, while our results are in agreement with several other publications, we are unable to reproduce those of Sor10 (), despite performing evolutions of the same initial data with the same gauge conditions. In particular we unambiguously find apparent horizons in data classified there as subcritical. The reason for this difference is not clear. Moving closer to the threshold of black hole formation, surprisingly, we find that two disjoint apparent horizons are found centered around some non-zero on the symmetry axis, indicating the likelihood that the Brill wave collapses to form a head-on collision of two black holes. The fact of, and time-scale for the merger of these horizons is to be determined. Finally we have bounded the critical amplitude within a range of about . This is an improvement over the previous bound by some orders of magnitude. We see evidence of power-law scaling, since the maximum of the Kretschmann scalar is well described by the form (6), as expected if critical phenomena are present GarDun98 (). The power-law exponent, at least, appears consistent with the value of Abrahams and Evans AbrEva93 (), but our tentative value of is very different from theirs . On the other hand our value of is compatible with that of the mixed axisymmetric gravitational wave and scalar-field data ChoHirLie03 (). Since we only see one period of the wiggle however, we must warn against premature jubilation. The wiggle could disappear, or the period of subsequent wiggles may differ substantially with more tuning of , so further work is definitely needed. Closer to the critical point we find more and more extreme behavior in the Kretschmann scalar. Particularly interesting are the ever-finer spikes that rapidly form on the symmetry axis. Superficially this even seems evocative of BKL type behavior.

Close to the critical point our current method suffers from larger errors, particularly in the form of constraint violation around the spikes and, we suspect, because coordinate singularities form. Evidently there is still much to understand. The next steps will include looking at different initial data, including Brill waves with , prolate and off-centered seed functions, along with the Teukolsky waves of AbrEva93 (), which suit better the original spirit of Cho93 () since they consist of incoming colliding waves. In the future we furthermore hope that the combination of mesh-refinement and the use of the dual-foliation Hil15 (); HilHarBug16 () approach will help to allay our current difficulties.

###### Acknowledgements.

We are grateful to David Garfinkle, Sascha Husa, Anton Khirnov, Tomáš Ledvinka and Hannes Rüter for interesting discussions, and especially to Carsten Gundlach for helpful comments on the manuscript. This work was supported in part by the FCT (Portugal) IF Program IF/00577/2015, the Deutsche Forschungsgemeinschaft (DFG) through its Transregional Center SFB/TR7 “Gravitational Wave Astronomy”, by the DFG Research Training Group 1523/2 “Quantum and Gravitational Fields”, and the Graduierten-Akademie Jena. Computations were performed primarily at the LRZ (Munich), but also on the Jena ARA cluster.## References

- [1] David Hilditch, Thomas W. Baumgarte, Andreas Weyhausen, Tim Dietrich, Bernd Brügmann, Pedro J. Montero, and Ewald Müller. Collapse of nonlinear gravitational waves in moving-puncture coordinates. Phys.Rev., D88(10):103009, 2013.
- [2] David Hilditch, Andreas Weyhausen, and Bernd Brügmann. Pseudospectral method for gravitational wave collapse. Phys. Rev., D93(6):063006, 2016.
- [3] Marcus Bugner, Tim Dietrich, Sebastiano Bernuzzi, Andreas Weyhausen, and Bernd Brügmann. Solving 3D relativistic hydrodynamical problems with WENO discontinuous Galerkin methods. Phys. Rev., D94(8):084004, 2016.
- [4] David Hilditch, Enno Harms, Marcus Bugner, Hannes Rüter, and Bernd Brügmann. The evolution of hyperboloidal data with the dual foliation formalism: Mathematical analysis and wave equation tests. 2016.
- [5] D. S. Brill. On the positive definite mass of the Bondi-Weber-Wheeler time-symmetric gravitational waves. Ann. Phys. (N. Y.), 7:466–483, 1959.
- [6] Andrew M. Abrahams and Charles R. Evans. Critical behavior and scaling in vacuum axisymmetric gravitational collapse. Phys. Rev. Lett., 70:2980, 1993.
- [7] A. M. Abrahams and C. R. Evans. Universality in axisymmetric vacuum collapse. Phys. Rev., D49:3998–4003, 1994.
- [8] Evgeny Sorkin. On critical collapse of gravitational waves. Class. Quant. Grav., 28:025011, 2011.
- [9] Carsten Gundlach and José M. Martín-García. Critical phenoena in gravitational collapse. Living Reviews in Relativity, 10(5), 2007.
- [10] Lee Lindblom, Mark A. Scheel, Lawrence E. Kidder, Robert Owen, and Oliver Rinne. A new generalized harmonic evolution system. Class. Quant. Grav., 23:S447–S462, 2006.
- [11] Oliver Rinne. Stable radiation-controlling boundary conditions for the generalized harmonic Einstein equations. Class. Quant. Grav., 23:6275–6300, 2006.
- [12] Oliver Rinne, Lee Lindblom, and Mark A. Scheel. Testing outer boundary treatments for the Einstein equations. Class. Quant. Grav., 24:4053–4078, 2007.
- [13] M. Bjørhus. The ODE Formulation of Hyperbolic PDEs Discretized by the Spectral Collocation Method. SIAM J. Sci. Comput., 16(3):542–557, 1995.
- [14] SpEC - Spectral Einstein Code, http://www.black-holes.org/SpEC.html.
- [15] C. Ronchi, R. Iacono, and P.S. Paolucci. The âcubed sphereâ: A new method for the solution of partial differential equations in spherical geometry. J. Comput. Phys., 124(1):93 – 114, 1996.
- [16] Jonathan Thornburg. Black hole excision with multiple grid patches. Classical Quantum Gravity, 21(15):3665–3691, 7 August 2004.
- [17] Bela Szilagyi, Lee Lindblom, and Mark A. Scheel. Simulations of binary black hole mergers using spectral methods. Phys. Rev., D80:124010, 2009.
- [18] Frans Pretorius. Evolution of binary black hole spacetimes. Phys. Rev. Lett., 95:121101, 2005.
- [19] M. Alcubierre, S. R. Brandt, B. Brügmann, D. Holz, E. Seidel, R. Takahashi, and J. Thornburg. Symmetry without symmetry: Numerical simulation of axisymmetric systems using Cartesian grids. Int. J. Mod. Phys. D, 10(3):273–289, 2001.
- [20] Kenneth R. Eppley. Evolution of time-symmetric gravitational waves: Initial data and apparent horizons. Phys. Rev. D, 16(6):1609–1614, 1977.
- [21] D. Holz, W. Miller, M. Wakano, and J. Wheeler. In B. Hu and T. Jacobson, editors, Directions in General Relativity: Proceedings of the 1993 International Symposium, Maryland; Papers in honor of Dieter Brill, Cambridge, England, 1993. Cambridge University Press.
- [22] Tim Dietrich and Bernd Brügmann. Solving the Hamiltonian constraint for 1+log trumpets. Phys.Rev., D89:024014, 2014.
- [23] Andreas Weyhausen. Numerical Methods for Collapsing Gravitational Waves. PhD thesis, Friedrich-Schiller-Universität Jena, December 2014.
- [24] K. Eppley. Pure gravitational waves. In L. L. Smarr, editor, Sources of Gravitational Radiation, pages 275–291. Cambridge University Press, Cambridge, 1979.
- [25] M. Alcubierre, G. Allen, B. Brügmann, E. Seidel, and W.-M. Suen. Towards an understanding of the stability properties of the 3+1 evolution equations in general relativity. Phys. Rev. D, 62:124011, 2000.
- [26] David Garfinkle and G. Comer Duncan. Numerical evolution of Brill waves. Phys. Rev. D, 63:044011, 2001.
- [27] Oliver Rinne. Axisymmetric Numerical Relativity. PhD thesis, University of Cambridge, Cambridge, England, 13 September 2005. gr-qc/0601064.
- [28] Lucia Santamaria. Nonlinear 3d evolutions of Brillwave sapcetimes and critical phenomena. Master’s thesis, Friedrich-Schiller-Universität Jena, October 2006.
- [29] Carsten Gundlach and Jose M. Martin-Garcia. Well-posedness of formulations of the Einstein equations with dynamical lapse and shift conditions. Phys. Rev. D, 74:024016, 2006.
- [30] Miguel Alcubierre. Introduction to 3+1 Numerical Relativity. Oxford University Press, Oxford, 2008.
- [31] David Hilditch, Sebastiano Bernuzzi, Marcus Thierfelder, Zhoujian Cao, Wolfgang Tichy, and Bernd Brügmann. Compact binary evolutions with the Z4c formulation. Phys. Rev. D, 88:084057, 2013.
- [32] Marcus Thierfelder, Sebastiano Bernuzzi, and Bernd Brügmann. Numerical relativity simulations of binary neutron stars. Phys. Rev. D, 84:044012, 2011.
- [33] Bernd Brügmann, José A. González, Mark Hannam, Sascha Husa, Ulrich Sperhake, and Wolfgang Tichy. Calibration of Moving Puncture Simulations. Phys. Rev. D, 77:024027, 2008.
- [34] T. W. Baumgarte, P. J. Montero, I. Cordero-Carrión, and E. Müller. Numerical relativity in spherical polar coordinates: Evolution calculations with the BSSN formulation. Phys. Rev. D, 87(4):044026, 2013.
- [35] Matthew W. Choptuik. Universality and scaling in gravitational collapse of massless scalar field. Phys. Rev. Lett., 70:9, 1993.
- [36] J. M. Martin-Garcia and C. Gundlach. All nonspherical perturbations of the choptuik spacetime decay. Phys. Rev. D, 59:064031, 1999.
- [37] Matthew W. Choptuik, Eric W. Hirschmann, Steven L. Liebling, and Frans Pretorius. Critical collapse of the massless scalar field in axisymmetry. Phys. Rev. D, 68:044007, 2003.
- [38] Stephen W. Hawking and George F. R. Ellis. The large scale structure of spacetime. Cambridge University Press, Cambridge, England, 1973.
- [39] Robert M. Wald. General relativity. The University of Chicago Press, Chicago, 1984.
- [40] Robert M. Wald and Vivek Iyer. Trapped surfaces in the Schwarzschild geometry and cosmic censorship. Phys. Rev. D, 44:R3719–R3722, 1991.
- [41] David Garfinkle and G. Comer Duncan. Scaling of curvature in subcritical gravitational collapse. Phys.Rev., D58:064024, 1998.
- [42] https://www.tpi.uni-jena.de/tiki-view_tracker_item.php?itemId=254.
- [43] Mark A. Scheel, Matthew Giesler, Daniel A. Hemberger, Geoffrey Lovelace, Kevin Kuper, Michael Boyle, Béla Szilágyi, and Lawrence E. Kidder. Improved methods for simulating nearly extremal binary black holes. Classical Quantum Gravity, 32(10):105009, 2015.
- [44] Daniel A. Hemberger, Mark A. Scheel, Lawrence E. Kidder, Bela Szilagyi, Geoffrey Lovelace, et al. Dynamical Excision Boundaries in Spectral Evolutions of Binary Black Hole Spacetimes. Class.Quant.Grav., 30:115001, 2013.
- [45] David Hilditch. Dual Foliation Formulations of General Relativity. 2015.