Aspherical deformations of the Choptuik spacetime
Thomas W. Baumgarte





Department of Physics and Astronomy, Bowdoin College, Brunswick, ME 04011, USA


We perform dynamical and nonlinear numerical simulations to study critical phenomena in the gravitational collapse of massless scalar fields in the absence of spherical symmetry.  We evolve axisymmetric sets of initial data and examine the effects of deviation from spherical symmetry.  For small deviations we find values for the critical exponent and echoing period of the discretely self-similar critical solution that agree well with established values; moreover we find that such small deformations behave like damped oscillations whose damping coefficient and oscillation frequencies are consistent with those predicted in the linear perturbation calculations of Martín-García and Gundlach.  However, we also find that the critical exponent and echoing period appear to decrease with increasing departure from sphericity, and that, for sufficiently large departures from spherical symmetry, the deviations become unstable and grow, confirming earlier results by Choptuik et.al..   We find some evidence that these growing modes lead to a bifurcation, similar to those reported by Choptuik et.al., with two centers of collapse forming on the symmetry axis above and below the origin.  These findings suggest that nonlinear perturbations of the critical solution lead to changes in the effective values of the critical exponent, echoing period and damping coefficient, and may even change the sign of the latter, so that perturbations that are stable in the linear regime can become unstable in the nonlinear regime.




Aspherical deformations of the Choptuik spacetime
    Thomas W. Baumgarte
    
=================================================






§ INTRODUCTION




Critical phenomena in gravitational collapse were first reported in Choptuik's seminal study of massless scalar fields <cit.>.  Specifically, Choptuik performed numerical simulations of the dynamical evolution of massless scalar fields coupled to Einstein's equations.  He considered different families of spherically symmetric initial data, parameterized by a parameter η, say.  Sufficiently strong initial data collapse to form a black hole, while sufficiently weak data disperse to infinity and leave behind flat space.   We refer to a value of the parameter η that separates supercritical from subcritical data as a critical parameter η_*.


Critical phenomena appear close to the black-hole threshold, in the vicinity of a critical parameter η_* (see, e.g., <cit.> for reviews).  For many different systems and matter models, and with η close to η_*, the dynamical evolution will, during an intermediate time, approximately follow a critical solution that contracts in a self-similar fashion.  This self-similar solution focusses on an accumulation event that occurs at a finite proper time τ_* as measured by an observer at the center.   At any proper time τ, the spatial scale R_ sc of the critical solution is therefore proportional to τ_* - τ,

    R_ sc≃τ_* - τ

(here and throughout we assume geometrized units, in which c = 1 and G = 1).   For massless scalar fields we refer to this self-similar critical solution as the Choptuik spacetime.  The better the fine-tuning to criticality, i.e. the closer η to η_*, the longer the evolution will follow the critical solution, and the later the time τ at which the evolution will depart from the critical solution.  The scale R_ sc≃τ_* - τ of the critical solution at this time τ imprints a length scale on the subsequent evolution, resulting, for example, in the famous scaling law

    M ≃ |η - η_*|^γ

for the black-hole mass M in supercritical evolutions.  Here γ is the critical exponent, which Choptuik found to be about γ≃ 0.37 in his original numerical experiments, independently of the parametrization of the initial data.

The critical exponent can also be found by considering perturbations of the critical solution.  In spherical symmetry, the departure from the critical solution is caused by an unstable, spherically symmetric mode that grows exponentially with exp(κ T), 
where

    T = - ln( τ_* - τ).

Note that τ→τ_* corresponds to T →∞.  The critical exponent can then be shown to be the inverse of the Lyapunov exponent κ, i.e. γ = 1/κ.  For a massless scalar field, Gundlach found γ = 0.374(1) from such a perturbative calculation <cit.>, in excellent agreement with Choptuik's numerical results.

Choptuik's discovery launched an entire new field of research, triggering a large number of studies, both numerical and perturbative, of critical phenomena for different matter models, spacetime dimensions and asymptotics (see <cit.> for reviews).  In particular, these studies showed that for some matter models, for example radiation fluids <cit.>, the critical solution displays a continuous self-similarity, while for others, including massless scalar fields, a discrete self-similarity.   A discretely self-similar solution performs an oscillation as it contracts towards the accumulation event.  In terms of the logarithmic time T, this oscillation has a period Δ, which can be determined by casting the problem as an eigenvalue problem <cit.>.  The value computed by <cit.> is Δ = 3.445452402(3).  At a time T + Δ, the self-similar solution will take the same shape as at time T, but on a radial scale that is smaller than that at T by a factor of 

    R_ sc(T + Δ)/R_ sc(T)    =    τ_* - τ(T + Δ)/τ_* - τ(T) = e^-(T + Δ)/e^-T
        =     e^- Δ∼1/31.3.

The reappearance of the solution at time intervals Δ is referred to as echoing, and Δ as the echoing period.  

Until recently most numerical studies of critical collapse assumed spherical symmetry, which is helpful for resolving the small spatial structures that emerge close to criticality (a notable early exception is the study by Abrahams and Evans of critical phenomena in the collapse of gravitational waves <cit.>).   A number of important questions, however, cannot even be addressed under the assumptions of spherical symmetry.  One such question concerns the stability of the critical solution to aspherical modes.  

Adopting a perturbative treatment, Martín-García and Gundlach <cit.> (hereafter MGG) found that all aspherical modes in the collapse of massless scalar fields are stable, leading them to conclude that “all nonspherical perturbations of the Choptuik spacetime decay".  We note, however, that the decay rate of the most slowly damped mode, an ℓ = 2 mode, was found to be quite small in magnitude, corresponding to a slow damping.

In <cit.>, Choptuik et.al. (hereafter CHLP) adopted a code in cylindrical coordinates (see also <cit.>) to study critical collapse of massless scalar fields in axisymmetry.  CHLP adopted a two-parameter family of initial data, with η parametrizing the overall strength of the data, and ϵ the departure from spherical symmetry (see Eq. <ref> below).   For a given value of ϵ, CHLP then fine-tuned η to the black-hole threshold (to about η/η_* - 1 ≃ 10^-15 for small values of ϵ) and studied the properties of the emerging critical solution.  

Briefly summarized, CHLP reported two key results.  One of their findings is that the critical exponent γ and the period Δ appear to depend on ϵ, with both of them decreasing with increasing departure from sphericity (see Table I in CHLP).  The other, perhaps more important finding is that, for large values of ϵ and for exquisite fine-tuning to criticality, the collapsing region in the center of the spacetime appears to bifurcate into two collapsing regions, located along the axis of symmetry.  This result indicates that there exists a non-spherical growing mode that dominates the evolution at late times T, in apparent conflict with the findings of MGG.  

Given the seemingly contradictory results of MGG and CHLP it is of interest to verify whether the results of CHLP can be reproduced with an independent numerical code.  However, only few codes have been able to simulate critical collapse in the absence of spherical symmetry, as most recent numerical relativity codes in three spatial dimensions have been designed for simulations of binary problems (but see <cit.> for attempts with such codes in the context of vacuum evolutions, as well as <cit.>  for examples of recent codes specifically designed for simulations of critical collapse in vacuum spacetimes).  In <cit.>, Healy and Laguna simulated critical collapse of scalar fields in three spatial dimensions (see also <cit.>).   Their calculations achieved more modest fine-tuning (to about η/η_* - 1 ≃ 10^-4), which was sufficient to measure the critical exponent γ and possibly the echoing period Δ, but not to make statements about the stability of the critical solution.  More recently, Deppe et.al. <cit.> performed similar simulations of critical collapse of scalar fields in three spatial dimensions, and achieved better fine-tuning (to about η/η_* - 1 ≃ 10^-6).  Interestingly, their reported values for the echoing period differ slightly from those found in previous numerical solutions as well as the semi-analytical value of <cit.>.  Deppe et.al. did not observe any growing modes, but it also is not clear whether their fine-tuning would be sufficient to find such modes.  

In this paper we adopt a code in spherical polar coordinates <cit.> to study critical collapse of massless scalar fields in the absence of spherical symmetry.  The code has been used previously for studies of critical collapse in ultrarelativistic fluids <cit.>.  In essence, we follow a suggestion made in the Conclusion section of CHLP, namely “... one could write a code adapted to the spherical critical solution (for instance using spherical polar coordinates with a logarithmic radial coordinate)."  We evolve the same axisymmetric initial data as those considered by CHLP (see Eq. <ref> below), and study the properties of near-critical solutions, fine-tuned to about η/η_* - 1 ≃ 10^-12 or better, for different deviations ϵ from sphericity.    

For small values of ϵ we find values of the critical exponent γ and the echoing period Δ that agree well with those found in most previous studies; furthermore we find that small deformations from sphericity behave as damped oscillations with a decay rate κ and oscillation frequencies similar to those computed by MGG.  We also confirm CHLP's finding that both γ and Δ appear to decrease with increasing ϵ.  For sufficiently large values of ϵ we find evidence for a growing unstable mode, in further agreement with CHLP's results.   We also present some evidence that this growing mode results in a bifurcation similar to that reported by CHLP, with two centers of oscillation forming on the symmetry axis above and below the origin.   These results suggest that, in the presence of nonlinear perturbations of the spherically symmetric critical solution, the evolution can be described in terms of effective values of γ, Δ and κ that depend on the departure from spherical symmetry.  The decay rate κ, which started out quite close to zero, may even change sign, thereby making a mode that was stable in the linear regime, as predicted by MGG, unstable for sufficiently large departures from spherical symmetry, as observed by CHLP.




This paper is organized as follows.  We present basic equations, our numerical method, initial data and diagnostics in Section <ref>. In Section <ref> we then present the results from our simulations, first under the assumption of spherical symmetry and then relaxing this assumption.  We briefly summarize in Section <ref>.  




§ SETUP OF THE PROBLEM







 §.§ Basic equations




We solve Einstein's equations 

    G_ab = 8 π T_ab,

where G_ab is the Einstein tensor, under the assumption that the stress-energy tensor T_ab is that of a massless scalar field (see Eq. (<ref>) below).  Specifically, we adopt the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation of Einstein's equations <cit.> in spherical polar coordinates (r,θ,φ) with the help of a reference-metric formulation <cit.>.  The BSSN formalism is based on a “3+1" decomposition of the spacetime, so that the line element takes the form

    ds^2 = g_ab dx^a dx^b = - α^2 dt^2 + γ_ij (dx^i + β^i dt)(dx^j + β^j dt).

Here g_ab is the spacetime metric, α the lapse function, γ_ij the spatial metric, and β^i the shift vector.   The normal on the spatial slices is given by

    n_a = (-α,0,0,0)  n^a = α^-1 (1, - β^i).

We further introduce a conformal decomposition of the spatial metric

    γ_ij = ψ^4 γ̅_ij,

where ψ is the conformal factor and γ̅_ij the conformally related metric.  For our applications in spherical coordinates we also identify the reference metric γ̂_ij with the flat metric η_ij expressed in spherical coordinates.  We refer the reader to <cit.> for details of this formalism, as well as its implementation in our code.

The massless scalar field ϕ satisfies the wave equation

    g^ab∇_a ∇_b ϕ = 0,

where ∇_a is the covariant derivative associated with the spacetime metric g_ab, and its
stress-energy tensor is given by

    T_ab = (∇_a ϕ)(∇_b ϕ) - 1/2 g_ab (∇_c ϕ)(∇^c ϕ).

Our convention for T_ab agrees with that used in many places in the literature, including <cit.> (as well as, for example, <cit.>), but differs from that of CHLP by a factor of 2 (see their Eq. 3).   In comparisons of our results with theirs the scalar field ϕ therefore needs to be rescaled with a factor of √(2).  




 §.§ Initial data




We impose initial data at an initial moment of time symmetry and assume that, at this time, the scalar field is given by the axisymmetric two-parameter family

    ϕ = ηexp(- (r / r_0)^2 (sin^2 θ + (1 - ϵ^2) cos^2 θ) )

(compare Eq. 9 in CHLP).  The parameter η determines the overall strength of the initial data, and ϵ^2 the departure from sphericity.   The constant r_0 has units of length; in practice we set r_0 = 1, meaning that we report all dimensional results in units of r_0.  



We also assume the spatial metric to be conformally flat initially, γ̅_ij = γ̂_ij, and, consistent with time symmetry, choose the shift vector β^i as well as time derivatives of the spatial metric and the scalar field to vanish.  The initial conformal factor ψ can then be found by solving the Hamiltonian constraint

    D̂^2 ψ = - 2 πψ^5 ρ

subject to Robin boundary conditions at the outer boundary.  Here D̂^2 is the flat Laplace operator associated with γ̂_ij, and ρ is the energy density

    ρ = n_a n_b T^ab.

We use an iterative process, alternately solving Eq. (<ref>) and updating (<ref>), to construct simultaneous solutions to the two equations.



We complete the specification of the initial data by choosing a “pre-collapsed lapse" α = ψ^-2 as the initial data for the lapse function α, as well as zero shift, β^i = 0.  




 §.§ Numerical Evolution




Our numerical code is very similar to that used in <cit.>, except that here we couple Einstein's equations to a massless scalar field instead of a fluid.  We evolve the gravitational fields using the code described in <cit.>.  The code uses finite-difference methods to solve the BSSN equations <cit.>, recast with the help of a reference-metric formulation <cit.>, in spherical coordinates.  We also rescale all tensorial variables, so that the coordinate singularities at the origin and on the axis can be handled analytically.  

We define a new auxiliary variable 

    Π≡ - n^a ∇_a ϕ = - 1/α (∂_t ϕ - β^i ∇_i ϕ )

in order to cast the wave equation (<ref>) as a pair of two equations that are first order in time.  Unlike in many other applications, however, we do not introduce variables that absorb first spatial derivatives of the scalar field ϕ, and instead leave the spatial derivatives in (<ref>) in terms of second derivatives – similar to the treatment of the gravitational fields in the BSSN equations.

The code does not make any symmetry assumptions, but for the axisymmetric simulations presented in this paper we set to zero all derivatives with respect to φ and use only a single grid point in the azimuthal direction.   We also impose an equatorial symmetry, and evolve only one of the two hemispheres.  For spherically symmetric simulations we use the minimum number of grid points possible in the θ-direction, N_θ = 2, while in the absence of spherical symmetry we use N_θ = 12 uniformly allocated angular grid points for ϵ^2 = 0.01 and 0.5, and N_θ = 14 for ϵ^2 = 0.75 (unless noted otherwise).

The radial grid points are allocated logarithmically (see Appendix A in <cit.>), so that each each grid cell is larger than its inner neighbor by a factor of c = 1.025.  For most simulations we use N_r = 312 grid points, and initially place the outer boundary at r_ out^ init = 64 (in our code units).   We allow the code to regrid whenever the length-scale l = (ϕ_,rr)^-1/2 becomes smaller than 25 Δ r, where Δ r is the innermost grid size.  In each regridding, the code variables are interpolated to a new grid that extends to a smaller outer boundary r_ out.  We allow up to 20 regrids, shrinking the outer boundary by equal factors down to r_ out^ final = 0.32.  We terminate all runs before the center of the simulation comes into causal contact with the outer boundary.   For larger ϵ^2 the accumulation event occurs at a later proper time τ_* (see Table <ref> below), and we therefore placed the outer boundary at r_ out^ init = 128, used N_r = 340 radial grid points (in order to achieve the same resolution at the center), and allowed 25 regrids to the same final outer boundary location of r_ out^ final = 0.32.  

For this latter setup, the ratio between the radial grid size at the outer boundary and that at the origin is c^N_r≃ 4430.  Further taking into account the regridding, we find that the ratio between initial grid size at the outer boundary and the final grid size at the origin is (r_ out^ init / r_ out^ final )   c^N_r≃ 1.77 × 10^6.   For comparison, a code with adaptive-mesh refinement would require about 20 levels of refinement to cover a similar range of scales, assuming refinements by factors of two between each level.  

We impose “moving puncture" coordinate condition during our evolution calculations.  Specifically, we adopt the 1+log slicing condition <cit.>

    (∂_t - β^i ∂_i) α = - 2 α K

for the lapse, where K is the mean curvature, and choose a Gamma-driver condition <cit.> for the shift, namely the version presented in <cit.>. 




 §.§ Diagnostics




We use several diagnostics to analyze the properties of our numerical solutions.  

As discussed in Section <ref>, we allowed regridding to small values of the outer boundary in order to be able to follow the self-similar solution to late times.  The disadvantage of this approach is that we had to terminate the simulations before the horizons of newly formed black holes had enough time to settle down to approximate equilibrium.  In this paper we therefore do not analyze the scaling of black hole masses for supercritical data.  As pointed out by <cit.>, however, scaling can also be observed for subcritical data.  The maximum encountered value of the central density ρ_c, for example, has units of inverse length, and therefore must satisfy a scaling law

    ρ_c^ max≃ |η_* - η|^-2 γ.

We note that the density ρ, defined in (<ref>), is in general not an invariant quantity, since it depends on the spacetime slicing through the normal vector n^a.  At the origin of the coordinate system, however, the symmetries in the problem single out a preferred normal vector, so that the values of the central density ρ_c reported here do take on an invariant meaning.

The scaling laws (<ref>) and (<ref>) describe the overall behavior of quantities close to criticality, but both Gundlach <cit.> and Hod and Piran <cit.> realized that, for scalar fields, the discretely self-similar nature of the critical solution imprints a periodic fine structure on the scaling law.  The scaling for the maximum encountered central density can then be written as

    lnρ_c^ max = C - 2 γln|η_* - η| + f(ln|η_* - η|),

where C is the logarithm of the constant of proportionality between ρ_c and the term |η_* - η|^-2γ in (<ref>), and where f(x) is a function that is periodic in x with angular frequency

    ω = Δ/2 γ.

In the fits presented in Section <ref> we assume that, to leading order, we can approximate f(x) as

    f(x) = A sin(ω x + ϕ_ ph),

where A is the amplitude of the fine structure and ϕ_ ph a phase, so that (<ref>) becomes

    lnρ_c^ max = C - 2 γln|η_* - η| + A sin(ωln|η_* - η| + ϕ_ ph).

Fits to numerical data then provide values for η_*, the critical exponent γ, as well as ω.  Using (<ref>) we can then compute Δ from ω and γ – the first of three different approaches to computing the period Δ.

We also adopt an approach similar to that of CHLP to monitor the departure from sphericity.  At logarithmic times T_0 (see Eq. <ref>) we launch pairs of photons from the center, one along the axis (θ = 0) and one in the equatorial plane (θ = π/2).   For each pair of photons we initialize the affine parameter λ to zero, and the derivative of λ to that of T, i.e. d λ/dt  = dT/dt (other normalizations are possible, but with this normalization λ inherits the self-similar nature of T).  We follow the photons' trajectories during the subsequent evolution, and record the value of the scalar field ϕ at each photon's current location as a function of its affine parameter λ.  We then compute the difference

    δϕ(T_0, λ) = 
    ϕ_ ax(λ) - ϕ_ eq(λ)

for a pair of photons launched at T_0, as well as the maxima of these differences

    Δϕ (T_0) ≡max_λ | δϕ(T_0, λ) |,

for each pair of photons emitted at logarithmic time T_0.  

For ϵ^2 > 0, a projection of our initial data (<ref>) into Legendre polynomials would result in contributions to modes of all even order ℓ.  However, we expect higher-order modes to decay more rapidly than lower-order modes (see Table I in MGG) so that at sufficiently late times our diagnostics δϕ and Δϕ become a measure of ℓ = 2 modes.




§ RESULTS







 §.§ Spherical symmetry




In order to both test and calibrate our code, and for better comparison with our aspherical results in Section <ref>, we start with a discussion of results in spherical symmetry, i.e. for ϵ = 0 in the initial data (<ref>).  Essentially the results presented in this Section reproduce those of <cit.> and numerous follow-up publications, including <cit.> who, similar to our treatment here, adopted the BSSN formalism and moving-puncture coordinates. 



We bracket the critical parameter η_* using bisection, except that in each step we refine the bracketing interval into ten new intervals instead of two by adding one more digit to η.   For subcritical data the scalar field disperses, leaving behind flat space, and the lapse function α approaches unity, while for supercritical data the minimum of the lapse function drops to zero, indicating the formation of a black hole.  We show an example in Fig. <ref>, which already suggests the presence of a pattern that keeps repeating on increasingly small scales.

While we can bracket the critical parameter η_* to about 13 significant digits, we would like to emphasize that these numbers do depend on the numerical grid.  The true uncertainty in η_* is significantly larger than the 13 digits suggest; see our discussion below.



In Fig. <ref> we show numerical results for the maximum encountered central density ρ_c as a function of η for subcritical data.   We fit the numerical data to (<ref>) and obtain values for η_*, γ and ω (as well as C, A and ϕ_ ph).  From (<ref>) we can then obtain our first estimate of the echoing period Δ. 

Several different sources of error affect the uncertainty in the parameters determined in these fits.   The critical amplitude η_* is mostly affected by numerical finite-difference error.  Based on comparisons with different grid setups we  expect this value to be accurate to within about 0.1% for the spherically symmetric simulations considered here.  The coefficients γ and Δ also depend more sensitively on which data are included in the fits – evolutions too close to the critical point develop features that can no longer be resolved on our grids, and therefore lead to large errors, while for evolutions far way from the critical point the expected power-law scalings no longer apply.  We crudely estimate the resulting errors to be about 1%; slightly less for γ, and slightly more for Δ.





We next consider the central value ϕ_c of the scalar field for a subcritical evolution close to the critical point.  When plotted as a function of proper time of an observer at the origin, as in Fig. <ref>, ϕ_c displays oscillations with decreasing period, as one would expect for a self-similarly contracting solution whose time and length scales continuously decrease.  The oscillations accumulate at the accumulation time τ_*.  

We can determine the accumulation time τ_* by looking for periodic behavior in the logarithmic time T (see Eq. (<ref>)).  In practice we identify the proper times τ_n of zero-crossings of ϕ_c.  We then consider a pair of subsequent zero-crossings, say τ_n and τ_n+1, and a second pair τ_m and τ_m+1.  Assuming that the advance of logarithmic time T for each pair is equal – namely Δ/2 – results in the estimate

    τ_* = τ_n τ_m+1 - τ_n+1τ_m/τ_n - τ_n+1 - τ_m + τ_m+1.

Equating the advance of T for each pair to Δ/2 then yields 

    Δ = 2 ln( τ_* - τ_n/τ_* - τ_n+1),

our second approach to determining Δ.  In practice we compute Δ for all pairs during the self-similar part of the evolution, and report average values in Table <ref>.  For different pairs, the values of Δ differ by up to approximately 1%, which we therefore take to be the approximate error in this value.   In Fig. <ref> we show the central value of the scalar field as a function of T, which clearly reveals the periodic behavior during an intermediate regime.  



Fig. <ref> helps distinguish three different phases during the evolution.  In Phase 1, the initial data approach the self-similar, critical solution.  This phase depends on the specifics of the initial data, but with sufficient fine-tuning it will result in the critical solution, the Choptuik spacetime, plus a perturbation whose amplitude depends on the degree of fine-tuning.     During Phase 2, the evolution can then be described as the critical solution plus the slowly growing perturbation.  Once the perturbation has become sufficiently large, the evolution departs from critical solution nonlinearly, marking the transition from Phase 2 to Phase 3, and the scalar field either disperses to infinity or collapses into a black hole.  The intermediate regime during which the periodic behavior can be observed in Fig. <ref> corresponds to Phase 2.



In Fig. <ref> we also label every other zero-crossing of ϕ_c, i.e. times T that, during Phase 2, are separated by a whole period Δ, as opposed to the half-periods discussed above.  We list both coordinate times and proper times of these zero-crossings in Table <ref>.    Note that N = 0 occurs just before the evolution enters Phase 2, while for N=3 the resolution at the center becomes quite poor.  According to the properties of the discretely self-similar solution, radial profiles of ϕ at these times should be echoes of each other, meaning that they should be identical to each other (in those parts of the solution that have reached self-similarity already), except rescaled by a factor e^-Δ between each echo (see Eq. <ref>).  In Fig. <ref> we show these profiles as a function of proper distance R from the origin, which we find by integrating from the origin to each grid point along lines of constant angles θ and φ.  For those parts of the profiles in the self-similar regime (and with sufficient numerical resolution) we find excellent agreement if we choose Δ = 3.45.  This provides us with a third measure of the echoing period Δ.   The agreement becomes visually worse if we increase or decrease Δ by about 1%.  

As an aside we note that it was not clear a priori that the discretely self-similar symmetry of the Choptuik spacetime would be revealed in our time slicing, but apparently the 1+log slicing (<ref>) does reflect the symmetry of the spacetime (see also the discussion in <cit.> in the context of ultra-relativistic fluids; see also <cit.> for proposals for symmetry-seeking spacetime coordinates).  





We report all our results in Table <ref>.  In particular we find that our different approaches to determining Δ yield values that are within our estimated margins of error.  Our values are also in good agreement with most previous numerical studies (e.g. <cit.>) as well as the values provided by <cit.>, who found Δ by casting the problem as an eigenvalue problem.  Similarly, our values of the critical exponent agree well with previous numerical results, as well as the perturbative values found by <cit.>.  In the following, including in Table <ref>, we will refer to the results of <cit.> as “quasi-analytical" results, even though solving the eigenvalue problem or the perturbative equations involves numerical work.




 §.§ Aspherical deformations







  §.§.§ Critical parameters



We next consider aspherical evolutions with ϵ > 0 in the initial data (<ref>).   For a given value of ϵ we bracket the critical value η_* as described in Section <ref>.  






For subcritical evolutions we can again fit the values of the maximum encountered density to the scaling relation (<ref>).  The result of these fits is shown in Fig. <ref>, which generalizes Fig. <ref> to include aspherical evolutions.  The overall trend - namely that the critical exponent γ appears to decrease with increasing ϵ – can already be observed in Fig. <ref>, but it is easier to see this after subtracting the overall scaling - 2 γ_0 ln(η_* - η)  (see Eq. <ref>) from each curve.   Here γ_0 = 0.374 is the value that we found in our spherical evolutions (see Table <ref>).  The resulting curves are shown in Fig. <ref> – note the similarity with the corresponding Fig. 3 in CHLP.   For ϵ^2 = 0.75 we included data only for relatively large values of ln(η_* - η).   Closer to the black-hole threshold we found larger deviations from these fits, which we believe are caused by a growing aspherical mode, as we will discuss in more detail below.

From the fits to (<ref>) we determine both γ and, via (<ref>), the echoing period Δ (see Table <ref>).   For small deformations, we estimate the errors to be similar to those that we found for spherically symmetric data.  For larger values of ϵ, two effects lead to increasing errors.  A stronger angular dependence leads to larger numerical finite-difference errors for a given angular resolution.  Moreover, for ϵ^2 = 0.75 we included fewer data points in the fits, as discussed above, which also increases the error.  We estimate the latter by varying the number of points included in the fits and find that, for large ϵ, the errors in γ and Δ are closer to about 5%.  We nevertheless note that both values decrease with increasing ϵ, in good agreement with the findings of CHLP (see their Table I).   

The critical exponent γ and the echoing period Δ describe properties of the spherically symmetric critical solution and its linear perturbations; strictly speaking, therefore, they take unique values that are independent of ϵ.   Apparently, however, the dynamics of nonlinear perturbations of the critical solution can still be described heuristically in terms of similar parameters γ and Δ, except that they now take on effective values that do depend on ϵ.  When we discuss the dependence of γ and Δ on ϵ in the following, we mean these effective values, rather than those defined in the context of linear perturbation theory.



As for our spherically symmetric evolutions in Section <ref>, we next consider the time dependence of the central value of ϕ for a near-critical evolution.   Identifying the proper times of zero-crossings during Phase 2 we can determine the accumulation time τ_* from (<ref>), and the echoing period Δ from (<ref>).  Our results for different values of ϵ are listed in Table <ref>.  For ϵ^2 = 0.01 and 0.5 we estimate the error in Δ to be again approximately 1%.  For ϵ^2 = 0.75, however, we observe fewer zero-crossings during Phase 2, and therefore believe that our error is close to about 5%.   As before we find that Δ decreases with increasing ϵ, in good agreement with our previous determination of Δ from the scaling of ρ_c^ max.  In Fig. <ref> we show ϕ_c as a function of logarithmic time T for different values of ϵ; this graph also shows that the period Δ becomes smaller for increasing ϵ.  





Next we again consider echoing in radial profiles of the scalar field ϕ at moments of alternating zero-crossings of its central value ϕ_c.   We again list the times of these zero-crossings in Table <ref>.  The echoing is now harder to identify than for the spherical evolutions, because the new angular-dependency leads to additional oscillations with periods different from Δ, meaning that the echoing is no longer exact (see also Eq. <ref> below).  For small values of ϵ the overall time dependence of ϕ is still dominated by that of the spherical Choptuik spacetime, so that the echoing can still be seen very cleanly, and Δ can still be determined quite accurately.  For ϵ^2 = 0.5 the deviations are larger, as we show in Fig. <ref>.  We now show each echo in its own panel, and include profiles of ϕ both in the axial direction (θ = 0, solid lines) and the equatorial direction (θ = π/2, dotted lines) in each panel.   While the profiles in the two directions do not agree, the differences do not appear to decay or grow with time (i.e. from echo to echo).  Moreover we can still identify the overall features of the Choptuik spacetime, which match quite well if we rescale each echo with Δ = 3.35.

The behavior becomes more complicated for ϵ^2 = 0.75, as shown in Fig. <ref>.  It now appears that the differences between the axial and equatorial profiles do grow – in fact, these deviations appear to cause an early departure from self-similar solution, so that the third echo N=3 is no longer in what we had identified as Phase 2 (see also Fig. <ref>).   We still observe some resemblance with the Choptuik solution for larger values of R, and were able to match these features reasonably well with Δ = 2.75 – clearly, however, this value should be considered a crude estimate only.   At smaller values of R, however, and on the axis, a new feature appears to have emerged by the time of the third echo.  

In Fig. <ref> we included results for two angular resolutions N_θ = 14 and N_θ = 12.  The two resolutions cannot be distinguished in the figure at early times, as well as at sufficiently large radius.  The new features that emerge late in the evolution, however, appear on the axis, and therefore have a strong angular dependence that is difficult to resolve in spherical coordinates.  While this leads to increasing quantitative differences between the different angular resolutions, the qualitative behavior remains very similar, as we discuss in more detail in the following Section.  (The “double peak" peak at e^3 Δ R ≃ 10 for N=3, θ = 0 and N_θ = 12 leads to a single peak very similar to that for N_θ = 14, just at a slightly later time; compare Fig. <ref>.)




  §.§.§ Properties of the deformations








In order to analyze the growth or decay of deviations from spherical symmetry more carefully we emit pairs of photons from the origin in the axial and equatorial directions, as described in Section <ref>.  We then record the differences δϕ in the scalar field ϕ that these pairs observe as a function of the photons' affine parameter λ (see Eq. <ref>), and also compute the maximum values of the differences Δϕ (see Eq. <ref>).   In Fig. <ref> we show results for δϕ for a near-critical evolution with ϵ^2 = 0.01, and in Fig. <ref> maximum values Δϕ for different values of ϵ (rescaled with ϵ^2).

For spherical evolutions, with ϵ = 0, the differences Δϕ always remain close to truncation error (about 5 × 10^-14 or less).  This behavior is different from that reported by CHLP, who find a growth in Δϕ even for ϵ = 0, although these values of Δϕ decrease with improving numerical accuracy (see their Fig. 6).  This effect might be related to a small drift of the center of symmetry with respect to the origin in the simulations of CHLP, which also converges to zero as numerical accuracy is improved.   We do not observe such a drift in our simulations with spherical polar coordinates and equatorial symmetry, nor a growth in Δϕ for ϵ = 0.  

According to MGG, linear perturbations of the Choptuik spacetime can be described by functions of the form

    u(T) ≃ e^κ Tcos(ω_p T + ϕ_ ph) g(T),

where κ determines the rate of perturbation's growth or decay.  The function g(T) is a function with the period Δ of the discretely self-similar background solution.   Values of κ and ω_p for different modes ℓ are tabulated in Table I of MGG.  In particular, MGG found that all κ are negative, leading them to conclude that “all nonspherical perturbations of the Choptuik spacetime decay."  

In general, the period 2 π / ω_p of the perturbation is not commensurate with the period Δ of the background, making it difficult to analyze theses functions, especially if data are available for only a small number of periods Δ (see also the discussion in MGG).   This situation is quite different from the critical collapse of ultrarelativistic fluids, for which the critical solution is continuously self-similar, so that that perturbations behave like (<ref>) but without the function g(T) (see <cit.>).  In this case the perturbations behave like damped or growing oscillations, and the coefficients κ and ω_p can be determined from fits to the numerical data (see <cit.>). 

While it is significantly harder to determine κ and ω_p for scalar fields, the curve for ϵ^2 = 0.01 in Fig. <ref> does appear like a slowly damped oscillation (the same damping is present in Fig. <ref>, but more difficult to see in a surface plot).  For comparison, we included the exponential exp(κ_ pert T) with κ_ pert  = - 0.07 / Δ, the value for an even ℓ = 2 mode according to MGG (see their Table I).   The slope of this curve appears to agree quite well with the overall decay of our numerical curve for Δϕ.  



Following MGG, we now multiply δϕ for ϵ^2 = 0.01 with exp(- κ_ pert T_0) in order to remove the overall exponential decay, and then take the Fourier transform of exp(- κ_ pert T_0)  δϕ for fixed values of λ > 0 (λ = 0 corresponds to the origin, where δϕ = 0 identically).  The low-frequency end of these Fourier transforms is shown in Fig. <ref>; for larger frequencies the spectra quickly decay to zero.  Given that we have data for T_0 ≲ 3 Δ only, the resolution of these spectra is quite crude.   We can nevertheless observe well-defined peaks, and moreover these peaks agree well with those predicted in the perturbative calculations of MGG, marked by the vertical dotted lines in Fig. <ref>.   We note that the value of ω_0 Δ / (2 π) = 0.3 provided in their Table I is by definition in the range between 0 and 2 (i.e. it is defined modulo two).  Peaks in the spectrum of the scalar field perturbations can then be found by adding odd multiples of 2 π / Δ to ω_0, yielding values of ω_p Δ / (2 π) = 1.3, 3.3, 5.3 …(compare their Fig. 2), while peaks in the spectrum of the metric perturbation can be found by adding even multiples of 2 π / Δ  (see also the discussion around their Eq. 173).  We conclude that, for small deformations ϵ^2, our results are consistent with both the slow exponential decay and the oscillation rates of MGG.

For larger deformations ϵ^2, however, we find deviations.  For ϵ^2 = 0.5, the curve in Fig. <ref> does not seem to decay at all, which is qualitatively consistent with our observations from Fig. <ref>.  For ϵ^2 = 0.75 the deviations grow, consistent with Fig. <ref>, which suggests that κ has changed sign and is now positive.

It may seem surprising that κ depends on the size ϵ of deviations from sphericity, since it describes linear perturbations of the unperturbed, spherically symmetric background solution, the Choptuik spacetime, and hence reflects properties of the latter.  However, we have already confirmed the observation of CHLP that other coefficients describing the properties of the background solution, namely γ and Δ, also appear to depend on ϵ.  In this sense the ϵ-dependence of κ is consistent with that of the others.  Unlike the others, the decay constant κ changes sign, but that may be a result of it starting out very close to zero.   We also note that a similar effect was observed for rotating ultrarelativistic fluids, where the decay rate of ℓ = 1 modes also appeared to show some dependence on the rotation rate (see <cit.>).   Strictly speaking, the coefficients κ, γ and Δ should be defined for small deviations from spherical symmetry only, so that the evolution can be approximated as the Choptuik spacetime plus a linear perturbation.  The results of CHLP, <cit.>, and ours here suggest that for larger values of the deviation from sphericity the evolution can, heuristically, still be described in terms of similar parameters, except that nonlinear terms in the deviation lead to changes in the effective values of these coefficients.  



We next present some evidence that the growing modes lead to a bifurcation in the solution, similar to that described by CHLP.  A first hint is offered by the lapse function α.  In Fig. <ref> we show graphs of the lapse function for subcritical evolutions close to criticality, both for ϵ^2 = 0.5 (top panels) and ϵ^2 = 0.75 (bottom panels).  In all panels we show the value of the lapse at the center as a function of coordinate time (α_c, blue solid lines), as well as the minimum values of the lapse along the axis (α_ min^ ax, red dashed lines) and in the equatorial plane (α_ min^ eq, green dotted lines).  For ϵ^2 = 0.5, the graph is qualitatively very similar to that in Fig. <ref> for spherical evolutions.  In particular we see that all three curves coincide at times when the lapse takes a minimum, indicating that the lapse takes its smallest values at the origin.   For ϵ^2 = 0.75, on the other hand, the behavior is qualitatively different.   The evolution no longer shows the same pattern repeating at increasingly small scales; instead the pattern changes at late times.  Moreover, we now see periods of time when the minimum value of the lapse along the axis no longer coincides with its value at the center.  This indicates that the lapse takes a smaller value somewhere on the symmetry axis, suggesting that a new collapsing region emerges away from the center.





This is consistent with the behavior of the scalar field ϕ, which, at late times, starts to develop large oscillations at a point on the symmetry axis.  The example shown in Fig. <ref> shows such an oscillation at times shortly after the N= 3 panel in Fig. <ref> where, as we discussed before, a similar feature can be seen.  In Fig. <ref> we show ϕ at the same times as surface plots.

Recall that we imposed equatorial symmetry in our simulations.  This means that if a feature develops at a certain distance from the origin on the symmetry axis (θ = 0 in Fig. <ref>, or the z̅-axis in Fig. <ref>), then the same feature also develops at the same distance on the other side of the equatorial plane.  Our simulations therefore suggest the formation of two centers of collapse on the symmetry axis, one on each side of the origin.  This behavior was first observed by CHLP, who referred to this process as a bifurcation.

The plots in Fig. <ref> also demonstrate, however, that our numerical grid is not sufficient to accurately resolve a small feature developing away from the origin – exactly as CHLP had anticipated: “Of course, spherical polar coordinates would not be well suited to following a solution beyond a bifurcation, but it should be adequate to study the growth or decay of perturbations."  We showed most results for ϵ^2 = 0.75 with N_θ = 14, but also compared with results for N_θ = 12, the resolution used in our other aspherical simulations.  In Fig. <ref>, for example, the two resolutions can hardly be distinguished except in the presence of the new feature that grows on the axis, and which leads to an increasing angular dependence.  We similarly include results for both angular resolutions in Fig. <ref>.  Note that these simulations differ not only in angular resolution, but also in precise distance from the critical parameter η_*, which makes a direct comparison difficult.  Despite the evident quantitative differences, we find very similar qualitative behavior, namely the formation of new centers of oscillation on the symmetry axis away from the origin.  While our numerical resolution is not sufficient to analyze these features in greater detail, our simulations suggest that, for large initial deviations from sphericity, perturbations of the Choptuik spacetime can grow, and that these growing perturbations may lead to a bifurcation in the solution as previously observed by CHLP.






































§ SUMMARY




We perform numerical simulations of the critical collapse of massless scalar fields in the absence of spherical symmetry.   We use a numerical relativity code that adopts spherical coordinates with a logarithmic radial coordinate and, following CHLP, evolve axisymmetric initial data whose deviation from sphericity is parameterized by ϵ.  

For small ϵ we find values for the critical exponent γ and the echoing period Δ that agree well with established values.  We also find that small deviations from spherical symmetry behave like slowly damped oscillations, with a decay rate κ and frequencies ω_p similar to those found by MGG in perturbative calculations.   

For larger deviations from spherical symmetry we find effective values for the critical exponent γ and the echoing period Δ that decrease with increasing ϵ, in good agreement with the findings of CHLP.    Moreover, we find that, for sufficiently large departures from spherical symmetry, the deviations grow, rather than decay, also confirming earlier results of CHLP.  This suggests that the effective value for the decay rate κ, like those for γ and Δ, depends on ϵ also, and, in fact, changes sign.  We find some evidence that, at late times, these growing modes lead to a bifurcation of the critical solution, with large oscillations of the scalar field developing around two points on the symmetry axis, one above and one below the origin.   While our numerical grid does not provide sufficient resolution away from the center to follow the evolution of these oscillations, and possibly their collapse to black holes (unlike the code of CHLP, our code does not employ AMR), our simulations nevertheless display the growth of these features, and hence support the original discovery of this bifurcation by CHLP.

A similar phenomenon was recently discussed by <cit.> in the context of collapse of axisymmetric gravitational waves (see also <cit.>).  Specifically, <cit.> evolved Brill wave initial data and  found that, close to criticality, the Kretschmann scalar takes maxima on the symmetry axis away from the origin, reminiscent of the results of CHLP as well as ours here.   Moreover, <cit.> find that, in their equatorially symmetric simulations, two black holes with disjoint apparent horizons form on the axis, suggesting that near-critical Brill wave initial data lead to the formation of two black holes which then collide head-on.   The authors further suggest that this behavior may dominate whenever axisymmetric modes dominate over spherically symmetric modes.

Unlike data for gravitational waves, the initial data (<ref>) can be produced as an expansion in ϵ about spherically symmetric data.  To linear order, one would therefore expect agreement with the perturbative results of MGG, and indeed, for small ϵ we observe a decay rate κ and oscillation frequencies ω_p that are similar to the values computed by MGG.  For larger ϵ, however, we observe shifts in the effective values of the critical exponent γ, the echoing period Δ as well as the decay rate κ.  Similar behavior was also observed by <cit.> for rotating ultrarelativistic fluids, even though there the effect of increasing rotation was to stabilize the perturbations, whereas here the effect of increasing deformations is to destabilize the perturbation [There exists a well-known equivalence between irrotational stiff fluids and scalar fields, which might mislead one to conclude that critical phenomena in one system ought to behave exactly as those in the other.  This is not the case, however; in fact, as clarified by <cit.>, the critical solution in one system is not observed in the other, and vice versa.].  The coefficients γ and κ have a well-defined meaning only in the context of linear perturbations.  Heuristically, however, it appears that a similar description of the dynamics is still possible in the presence of nonlinear effects, except that the effective coefficients γ and κ then depend on the size of the perturbation.  Exploring this possibility quantitatively would require (at least) a second-order perturbation analysis, which, unfortunately, is quite challenging. 





It is a great pleasure to thank Carsten Gundlach for numerous elucidating conversations and for suggesting the spectral analysis in Fig. <ref>.  I would also like to thank him, David Hilditch, Steve Liebling and Frans Pretorius for detailed and very helpful comments on earlier versions of this manuscript.   This work was supported in part by NSF grants PHY-1402780 and 1707526 to Bowdoin College, and through sabbatical support from the Simons Foundation (Grant No. 561147 to TWB).  Numerical simulations were performed on the Bowdoin Computational Grid.
44
natexlab#1#1bibnamefont#1#1bibfnamefont#1#1citenamefont#1#1url<#>1urlprefixURL 


[Choptuik(1993)]Cho93
authorM. W. Choptuik,
  journal volume70, pages9
  (year1993).

[Gundlach(2003)]Gun03
authorC. Gundlach,
  journalPhys. Rept. volume376,
  pages339 (year2003).

[Gundlach and Martín-García(2007)]GunM07
authorC. Gundlach and
  authorJ. M. Martín-García,
  journalLiving Rev. Rel. volume10,
  pages5 (year2007).

[Gundlach(1997)]Gun97
authorC. Gundlach,
  journal volume55, pages695
  (year1997).

[Evans and Coleman(1994)]EvaC94
authorC. R. Evans and
  authorJ. S. Coleman,
  journal volume72, pages1782
  (year1994).

[Gundlach(1995)]Gun95
authorC. Gundlach,
  journal volume75, pages3214
  (year1995).

[Martín-García and
  Gundlach(2003)]MarG03
authorJ. M.
  Martín-García and
  authorC. Gundlach,
  journal volume68, eid024011
  (year2003).

[Abrahams and Evans(1993)]AbrE93
authorA. M. Abrahams
  and authorC. R.
  Evans, journal
  volume70, pages2980 (year1993).

[Martín-García and Gundlach(1999)]MarG99
authorJ. M. Martín-García
  and
  authorC. Gundlach,
  journal volume59,
  pages064031/1 (year1999).

[Choptuik
  et al.(2003a)Choptuik, Hirschmann,
  Liebling, and Pretorius]ChoHLP03b
authorM. W. Choptuik,
  authorE. W. Hirschmann,
  authorS. L. Liebling,
  and
  authorF. Pretorius,
  journal volume68,
  pages044007/1 (year2003a).

[Choptuik
  et al.(2003b)Choptuik, Hirschmann,
  Liebling, and Pretorius]ChoHLP03a
authorM. W. Choptuik,
  authorE. W. Hirschmann,
  authorS. L. Liebling,
  and
  authorF. Pretorius,
  journal volume20, pages1857
  (year2003b).

[Alcubierre et al.(2000)Alcubierre,
  Allen, Brügmann, Lanfermann, Seidel, Suen, and
  Tobias]AlcABLSST00
authorM. Alcubierre,
  authorG. Allen,
  authorB. Brügmann,
  authorG. Lanfermann,
  authorE. Seidel,
  authorW.-M. Suen,
  and authorM. Tobias,
  journal volume61,
  pages041501(R) (year2000).

[Hilditch et al.(2013)Hilditch,
  Baumgarte, Weyhausen, Dietrich, Brügmann, Montero, and
  Müller]HilBWDBMM13
authorD. Hilditch,
  authorT. W. Baumgarte,
  authorA. Weyhausen,
  authorT. Dietrich,
  authorB. Brügmann,
  authorP. J. Montero,
  and
  authorE. Müller,
  journal volume88,
  pages103009/1 (year2013).

[Sorkin(2010)]Sor10
authorE. Sorkin,
  journal volume81,
  pages084062 (year2010).

[Sorkin(2011)]Sor11
authorE. Sorkin,
  journal volume28,
  pages025011 (year2011).

[Hilditch et al.(2016)Hilditch,
  Weyhausen, and Brügmann]HilWB15
authorD. Hilditch,
  authorA. Weyhausen,
  and
  authorB. Brügmann,
  journal volume93, eid063006
  (year2016).

[Hilditch et al.(2017)Hilditch,
  Weyhausen, and Brügmann]HilWB17
authorD. Hilditch,
  authorA. Weyhausen,
  and
  authorB. Brügmann,
  journal volume96, eid104051
  (year2017).

[Healy and Laguna(2014)]HeaL14
authorJ. Healy and
  authorP. Laguna,
  journalGen. Relativ. Gravit. volume46,
  pages1722/1 (year2014).

[Clough and Lim(2016)]CloL16
authorK. Clough and
  authorE. A. Lim,
  journalArXiv e-prints  (year2016),
  1602.02568.

[Deppe et al.(2018)Deppe, Kidder,
  Scheel, and Teukolsky]DepKST18
authorN. Deppe,
  authorL. E. Kidder,
  authorM. A. Scheel,
  and authorS. A.
  Teukolsky, journalArXiv e-prints
  (year2018), 1802.08682.

[Baumgarte et al.(2013)Baumgarte,
  Montero, Cordero-Carrión, and Müller]BauMCM13
authorT. W. Baumgarte,
  authorP. J. Montero,
  authorI. Cordero-Carrión,
  and
  authorE. Müller,
  journal volume87,
  pages044026/1 (year2013).

[Baumgarte et al.(2015)Baumgarte,
  Montero, and Müller]BauMM15
authorT. W. Baumgarte,
  authorP. J. Montero,
  and
  authorE. Müller,
  journal volume91,
  pages064035/1 (year2015).

[Baumgarte and Montero(2015)]BauM15
authorT. W. Baumgarte
  and authorP. J.
  Montero, journal
  volume92, pages124065/1
  (year2015).

[Baumgarte and Gundlach(2016)]BauG16
authorT. W. Baumgarte
  and
  authorC. Gundlach,
  journal volume116,
  pages221103 (year2016).

[Gundlach and Baumgarte(2016)]GunB16
authorC. Gundlach and
  authorT. W. Baumgarte,
  journal volume94,
  pages084012 (year2016).

[Gundlach and Baumgarte(2018)]GunB18
authorC. Gundlach and
  authorT. W. Baumgarte,
  journal volume97,
  pages064006 (year2018).

[Nakamura et al.(1987)Nakamura,
  Oohara, and Kojima]NakOK87
authorT. Nakamura,
  authorK. Oohara, and
  authorY. Kojima,
  journalProg. Theor. Phys. Suppl. volume90,
  pages1 (year1987).

[Shibata and Nakamura(1995)]ShiN95
authorM. Shibata and
  authorT. Nakamura,
  journal volume52, pages5428
  (year1995).

[Baumgarte and Shapiro(1998)]BauS98
authorT. W. Baumgarte
  and authorS. L.
  Shapiro, journal
  volume59, pages024007/1
  (year1998).

[Bonazzola et al.(2004)Bonazzola,
  Gourgoulhon, Grandclément, and Novak]BonGGN04
authorS. Bonazzola,
  authorE. Gourgoulhon,
  authorP. Grandclément,
  and authorJ. Novak,
  journal volume70,
  pages104007/1 (year2004).

[Shibata et al.(2004)Shibata,
  Uryū, and Friedman]ShiUF04
authorM. Shibata,
  authorK. Uryū,
  and authorJ. L.
  Friedman, journal
  volume70, pages044044/1
  (year2004).

[Brown(2009)]Bro09
authorJ. D. Brown,
  journal volume79,
  pages104029/1 (year2009).

[Gourgoulhon(2012)]Gou12
authorE. Gourgoulhon,
  title3+1 formalism in general relativity
  (publisherSpringer, New York, year2012).

[Hawking and Ellis(1973)]HawE73
authorS. W. Hawking
  and authorG. F. R.
  Ellis, titleThe large scale structure of
  space-time (publisherCambridge University Press, Cambridge,
  year1973).

[Akbarian and Choptuik(2015)]AkbC15
authorA. Akbarian and
  authorM. W. Choptuik,
  journal volume92,
  pages084037/1 (year2015).

[Celestino and Baumgarte(2018)]CelB18
authorJ. Celestino and
  authorT. W. Baumgarte,
  journalArXiv e-prints  (year2018),
  1805.10442.

[Bona et al.(1995)Bona, Massó,
  Seidel, and Stela]BonMSS95
authorC. Bona,
  authorJ. Massó,
  authorE. Seidel, and
  authorJ. Stela,
  journal volume75, pages600
  (year1995).

[Alcubierre et al.(2003)Alcubierre,
  Brügmann, Diener, Koppitz, Pollney, Seidel, and
  Takahashi]AlcBDKPST03
authorM. Alcubierre,
  authorB. Brügmann,
  authorP. Diener,
  authorM. Koppitz,
  authorD. Pollney,
  authorE. Seidel, and
  authorR. Takahashi,
  journal volume67,
  pages084023 (year2003).

[Thierfelder et al.(2011)Thierfelder,
  Bernuzzi, and Brügmann]ThiBB11
authorM. Thierfelder,
  authorS. Bernuzzi,
  and
  authorB. Brügmann,
  journal volume84,
  pages044012 (year2011).

[Garfinkle and Duncan(1998)]GarD98
authorD. Garfinkle and
  authorG. C. Duncan,
  journal volume58,
  pages064024/1 (year1998).

[Hod and Piran(1997)]HodP97b
authorS. Hod and
  authorT. Piran,
  journal volume55, pagesR440
  (year1997).

[Garfinkle and Gundlach(1999)]GarG99
authorD. Garfinkle and
  authorC. Gundlach,
  journal volume16, pages4111
  (year1999).

[Gundlach(2002)]Gun02
authorC. Gundlach,
  journal volume65,
  pages084021/1 (year2002).

[Brady et al.(2002)Brady, Choptuik,
  Gundlach, and Neilsen]BraCGN02
authorP. R. Brady,
  authorM. W. Choptuik,
  authorC. Gundlach,
  and authorD. W.
  Neilsen, journal
  volume19, pages6359 (year2002).




