In a companion paper [1] (from now on Paper I), we have reviewed various formulations of the Einstein equations in twist-free axisymmetry on outgoing null cones emanating from a regular center, intending to simulate gravitational collapse with them.
As a first physics application, here we reinvestigate the problem of twist-free axisymmetric scalar field critical collapse. By critical collapse we understand the investigation of the threshold, along a 1-parameter family of initial data, between regular data that do and do not form a black hole in their time evolution. In "type-II" critical collapse, increased fine-tuning to the threshold p ¼ p à gives rise to arbitrarily small black hole masses, arbitrarily large curvatures, and in the limit, it is conjectured, to a naked singularity [2].
We choose scalar field critical collapse because it provides a continuous bridge from spherical symmetry, where simulations on null cones work well, to the twist-free axisymmetric vacuum case that is our long-term goal. But it is also of physical interest in its own right.
In [3] it was claimed, based on the numerical solution of a mode ansatz, that all nonspherical linear perturbations of the spherical scalar field critical solution decay, with the least damped l ¼ 2 perturbation only decaying quite slowly. Subsequent numerical time evolutions in axisymmetry in cylindrical coordinates by Choptuik et al. [4] appeared to indicate a slowly growing l ¼ 2 perturbation of slightly nonspherical critical collapse that eventually leads to the formation of two centers of collapse.
This tension between the results of [3 and [4] was resolved by fresh simulations, in axisymmetry in spherical coordinates by Baumgarte [5], which indicated that sufficiently small nonspherical perturbations decay (in agreement with the mode stability claimed in [3]) but that larger perturbations do grow (in agreement with [4]).
Interestingly, both [4] and [5] observe that in a regime of small finite deviations from spherical symmetry, the spherical Choptuik solution [6,7] is still observed as an approximate critical solution but with a period that decreases from Δ ≃ 3.44 with increasing nonsphericity, and that the critical exponent also decreases from γ ≃ 0.374.
The authors of [4] and [5], in axisymmetry, achieved fine-tuning to jpp à j ∼ 10 -15 and jpp à j ∼ 10 -13 , respectively (except that [5] achieved only 10 -8 at the largest deviation from spherical symmetry). Other investigations of nonspherical scalar field collapse had achieved much less fine-tuning: [8] achieved ∼10 -3.5 , and [9] ∼10 -5. 5 , both for nonaxisymmetric initial datap, and [10] achieved ∼10 -1 in axisymmetry. Incidentally, [8][9][10] also evolve initial data that are nonsmooth at R ¼ 0.
Reid and Choptuik [11] have fine-tuned axisymmetric, time-symmetric, antisymmetric in z initial data to ∼10 -14 and found that two (mirror-symmetric) centers of collapse form (as one would expect). These start out in a highly nonspherical state (although this was not quantified). However, critical phenomena were observed with Δ and γ as in spherical symmetry, and it was concluded that there is no evidence of a growing nonspherical mode (in the evolution of each separate center of collapse).
Marouda et al. [12] have examined twist-free axisymmetric collapse of a complex massless scalar field, finetuning to ∼10 -10 . As they treat a slightly different matter model, none of their families of initial data coincides exactly with those of [4] and [5], but their results are consistent with those of [4] and [5] in the sense that δΔ and δγ observed in their two nonspherical families fit into a monotonic ordering of families by (inferred) nonsphericity. See already Table II below.
The achievable level of fine-tuning matters, as the l ¼ 2 unstable mode claimed in [4] would be very slowly growing: its authors see no bifurcation for nonsphericity parameter ϵ 2 ≤ 2=3 [see Eq. ( 10) below for a definition], but they see one after 2 and 1.5 echoes, respectively, for ϵ 2 ¼ 3=4 and 5=6, where 3 echoes is the best that can be achieved in double-precision arithmetic, with jpp à j ∼ 10 -15 . [12] see about 4 echoes in their family IV before bifurcation.
In principle, being in polar coordinates, our code, as presented in detail in Paper I, should resolve small deviations from spherical symmetry as well as that of [5], whereas our choice of null gauge should be able to resolve the critical solution well, even without the adaptive mesh refinement of [4]. It turns out that this is true.
The plan of the paper is as follows: In Sec. II we briefly summarize our metric ansatz, gauge choice and diagnostics, give our two-parameter family (strength and nonsphericity) of initial data, and describe the bisection to the threshold of collapse. In Sec. III we describe in detail our results for spherical critical collapse of a massless scalar field. In particular, we plot the fine structure of the curvature and black hole mass scaling laws with high accuracy and compare with the available literature.
Section IV gives our results for nonspherical (twist-free axisymmetric) critical collapse. Our initial data, being set on an outgoing null cones, cannot be compared directly with those of [4,5]. However, by comparing the amplitude of the l ¼ 2 deviation of the scalar field ψ from spherical symmetry during the phase of the evolution where the solution is well approximated by the Choptuik solution, we infer that our most nonspherical family, ϵ 2 ¼ 0.75, compares to ϵ 2 ¼ 0.5 of [4,5]. With the gauge choice presented here we cannot fine-tune families of initial data that are more nonspherical.
We conclude in Sec. V with a discussion of the physical results and an outlook for our code.
We state here the equations we solve numerically, with full details of the formulation and its discretization given in Paper I. The field equations we want to solve are the Einstein equations
and the massless, minimally coupled wave equation
We use units where c ¼ G ¼ 1.
We write the general twist-free axisymmetric metric in the form
We assume that the central worldline R ¼ 0 is at x ¼ 0 and that spacetime is regular there. We have defined y ≔cos θ, so that the range 0 ≤ θ ≤ π corresponds to -1 ≤ y ≤ 1. The azimuthal angle φ has range 0 ≤ φ < 2π.
The Killing vector generating the axisymmetry is ∂ φ . (We use the convention of equating vector fields with derivative operators.) In (3) we have used the shorthand
Each surface N þ u of constant u is an outgoing null cone, assumed to have a regular vertex. Each surface S u;x of constant u and x is assumed to be spacelike, and has topology S 2 . The outgoing future-directed null vector field normal to S u;x is U ≔ G -1 ∂ x , and is also tangent to the affinely parametrized generators of N þ u . The ingoing future-directed null vector normal to S u;x is
where we have defined the shorthand
As we see from (5) B plays the role of a "shift" in the x-direction, with B and b relating our time direction ∂ u to the ingoing null direction Ξ. U and Ξ are normalized relative to each other as Ξ a U a ¼ -1. See Fig. 1 of Paper I for a sketch of our coordinate null cones and these vector fields.
As described in Paper I, we can solve a subset of the Einstein equations and the wave equation (which we call the hierarchy equations) for G, b, ΞR, Ξf and Ξψ on one time slice, given R, f and ψ there, and this can be done by explicit integration in x in the right order. H and ∂ u do not appear in the hierarchy equations except in the combination Ξ, and H remains undetermined.
This means that we have a time evolution scheme where R, f and ψ are specified at u ¼ 0. We then find G, b, ΞR, Ξf and Ξψ, choose B freely and thus obtain R ;u , f ;u and ψ ;u , and then evolve R, f and ψ forward in u.
We discretize in x and u using finite differencing, and in y using a pseudospectral expansion in Legendre polynomials (that is, spherical harmonics restricted to axisymmetry).
We choose H within a class of gauge choices with the property that R ¼ 0 remains at x ¼ 0 (and is timelike and regular), and that x ¼ x max is future spacelike. This means that our numerical domain is a subset of the domain of dependence of the initial data on 0 ≤ x ≤ x max , u ¼ 0, and no boundary condition is required on the outer boundary
so that the proper time along this "central wordline" is given by u.
Furthermore, we impose that x ¼ x 0 is marginally future spacelike, in the sense that Hðu; x 0 ; yÞ ≤ 0, with equality for some y, for some fixed constant x 0 in the range 0 < x 0 ≤ x max . Our expectation is that, if a part of the spacetime we are constructing numerically is approximated by a self-similar critical spacetime (with accumulation point on the central worldline x ¼ 0), then we can (by trial and error) adjust x 0 so that x ¼ x 0 remains close to the past light cone of the accumulation point of scale echoes. The resulting coordinate system will then zoom in on the accumulation point and maintain resolution of the ever smaller spacetime features as it is approached, without the need for adaptive mesh refinement. This approach to critical collapse in null coordinates was used in [13,14], using regridding, and in [15,16] using a shift term, all in spherical symmetry.
More specifically, for this paper we have settled on a class of gauges we call "local shifted Bondi gauge" (from now on, lsB gauge), defined by Rðu; x; yÞ ¼ Rðu; xÞ. Within that class, we choose the "lsB4" flavor, defined by the shift B defined in (6) taking the value
(Note that in Paper I we extensively tested the slightly different lsB2 gauge.) In this gauge, every surface of constant x > x 0 is future spacelike. We refer the reader to Paper I for more details. Note that Bðu; x; yÞ is continuous but has discontinuous transverse derivative across the hypersurfaces x ¼ x à ðuÞ where the y-location of min y ðΞRÞ jumps.
For any closed spacelike 2-surface S, we define its "Hawking compactness"
where ρ þ and ρ -are the outgoing and ingoing null divergence. It is related to the well-known Hawking mass by
where AðSÞ is the area of S.
The standard indicator of black hole formation used in numerical relativity on a spacelike time slicing is the appearance of a marginally outer-trapped surface (from now on, MOTS) embedded in a time slice. An outer trapped surface is defined by ρ þ ≤ 0 at every point (and hence C > 1), and a MOTS by ρ þ ¼ 0 (and hence C ¼ 1). As we have discussed in Paper I, we cannot expect the ingoing past light cone of a MOTS to converge to a point, and so we cannot expect to find any MOTS embedded in a coordinate null cone.
Instead we use the Hawking compactness Cðu; xÞ of our coordinate 2-surfaces S u;x as an indicator of black-hole formation. When max x Cðu; xÞ ≥ 0.99 is first reached during an evolution, we use the Hawking mass M of the S u;x where that happens as an indication of the "initial mass" of the black hole. Similarly, we use max x Cðu; xÞ ≤ 0.01 as an indicator of dispersion. In spherical symmetry, FIG. 1. ϵ 2 ¼ 0: Numerical determination of the curvature scaling function f T . Here and in all following plots of f T , the horizontal axis shows lnðp à -pÞ þ B and the vertical axis shows -1=2 ln T -A. To make the plots larger, we omit axis labels in all line plots throughout. Circles indicate our data points, 30 per decade in pp à . We show only the range where the function is approximately periodic. The gray line represents a fit as a sine wave of amplitude 0.147. For comparison, in all following plots of f T and f M the range of lnðp à -pÞ þ B, will be the same as here, namely ½-30.8; -5.8.
this works perfectly well: C ¼ 1 is actually equivalent to a MOTS, and this MOTS is approached uniformly at all angles ðθ; φÞ. Beyond spherical symmetry, for current lack of a better alternative, we still use Cðu; xÞ to distinguish between dispersion and collapse.
However, CðSÞ ≥ 1 is clearly only necessary, not sufficient, for S to be outer-trapped, and we are not aware of any rigorous results linking 2-surfaces with C ≃ 1 to black holes, nor of their previous use in the numerical relativity literature beyond spherical symmetry.
In lsB gauge, where R ¼ Rðu; xÞ, the Hawking mass Mðu; xÞ of the coordinate 2-surfaces S u;x obeys M ;x ðu; xÞ ≥ 0. There are separate integrals for Cðu; xÞ and Mðu; xÞ, and this gives to two separate ways of computing C, whose numerical results we denote by C and C, and similarly M and M. Their agreement is strong test of numerical accuracy, as they are discretized in different ways, see Paper I for details.
The authors of [4] investigate nonspherical scalar field collapse in axisymmetry, using two families of initial data on a Cauchy surface. One is time-symmetric, with the initial value of the scalar field a Gaussian elongated along the rotation axis, with equatorial symmetry. The other is approximately ingoing, with ψ antisymmetric under reflections through the equator. Like the code of [5], ours is optimized for a single center of collapse, while the antisymmetric data have two, so we evolve only a family of data designed to be similar to the first family of [4].
The time-symmetric initial data of [4] and [5] at t ¼ 0, written in our notation, are ϕð0; r; yÞ ¼ pe
Here r is the standard radial coordinate in a 3-metric assumed to be conformally flat at t ¼ 0. With ϵ 2 > 0, the initial data are elongated along the symmetry axis y ¼ AE1, and with ϵ 2 < 0 they are squashed. In [4,5], ϵ 2 is written as ϵ 2 and only positive values are considered, but there is no reason not to consider ϵ 2 < 0, hence our change of notation. Given the initial data for the scalar field, [4,5] obtain full initial data for the Einstein-scalar system by making the initial 3-metric conformally flat and solving the Hamiltonian constraint for the conformal factor.
As we set initial data on an outgoing null cone u ¼ 0 rather than a time slice t ¼ 0, we cannot construct initial data giving exactly the same solutions. Instead we identify the two datasets as different cross sections of an (approximate) solution of the flat spacetime wave equation. For this purpose, we consider the function ϕðt; r; yÞ ≔ 1 2r
For ϵ 2 ¼ 0 only, this a spherically symmetric analytic solution of the scalar wave equation on flat spacetime. For any ϵ 2 , it also reduces to (10) at t ¼ 0. We now identify the null cone u ¼ 0 where we impose initial data with the conical cross section r ¼ tt 0r ¼ 0 of (11) for some constant t 0 , and we identify r with R on this cross section. Finally, we set an overall amplitude p. Hence we set ψð0; R; yÞ
where ϕ is defined by (11). We also set
This is simply for lack of an alternative likely to be closer to the corresponding data at t ¼ 0 being conformally flat, and is at least consistent with spacetime being flat in the limit p → 0. Finally, we set Rð0; x; yÞ ¼ x
which is essentially a gauge choice. Formally, as p → 0 and ϵ 2 → 0, the solution from these data corresponds to that arising from (10), with the identification u ¼ t -R. We also initialize Rð0; x; yÞ ¼ x=2, which can be considered our initial choice of an lsB gauge.
Obviously, with gravity and/or with ϵ 2 ≠ 0, (11) is no longer an exact solution of the curved-spacetime wave equation, and the solution evolved from the null initial data (12) no longer has a moment of exact time symmetry, but it should approximate the solutions arising from (10), with ϵ 2 playing a quantitatively similar role.
Without loss of generality we set d ¼ 1 to fix an overall scale. We also choose t 0 ¼ -5, which means that the null initial data (12) represents a Gaussian that is well separated from the center and still essentially ingoing. In the flatspace, spherical limit ϵ 2 ¼ 0, p → 0 the wave reaches the center at u ¼ -t 0 , and the moment of time symmetry t ¼ 0 corresponds to u þ R ¼ -t 0 . Hence (always in this limit) our null initial data u ¼ 0 intersect t ¼ 0 at R ¼ -t 0 . For ϵ 2 ¼ 0, we shall choose x max ¼ 11, so that R ¼ -t 0 intersects our initial data surface.
We realized only near the completion of this paper, from the convergence tests in Appendix B, that for ϵ 2 ≠ 0 and t ≠ 0 (11) is not actually single-valued at the origin r ¼ 0, but depends on y. Our null initial data (12) are therefore also not single-valued at the origin. In Appendix C we present an analytic solution of the flat space wave equation valid for all values of ϵ 2 that also reduces to (10) for all values of ϵ 2 , and null initial data (C6) derived from this solution. These are the initial data we intended to create.
However, convergence tests in the strong-gravity regime of the nonanalytic null data (12) with (11) show that the breakdown of convergence is limited to a small region near the origin, at early times. The code seems to smooth out the data to something that then converges. The corrected null data (C6) also converge, but preliminary tests show that at ϵ 2 ¼ 0.75 they require very large values of x max in nearcritical evolutions. For both these reasons-the error in the nonanalytic data does not seem to matter, and the analytic data are much harder to evolve-we present here only evolutions with the nonanalytic initial data, and we believe that the following results are still reliable.
Beginning from a value of p that disperses and one that collapses, we find the threshold value p à by bisection. We use supern as shorthand for lnðpp Ã Þ ≃ -n and subn for lnðp à -pÞ ≃ -n. (This makes sense because p is dimensionless and p à ∼ 1.) We can find p à to machine precision (15 digits), and the accumulation point of echoes u à to about 7 digits, but our numerical values of p à , u à and x à will agree with their continuum value to fewer digits. This is a well-known feature of numerical simulations of critical collapse. Below we give their values at different resolutions as a rough indication of numerical error.
Beyond spherical symmetry, during bisection in p à , or later when we sample ln jpp à j more finely, we occasionally obtain inconclusive evolutions. For small ϵ 2 at least there is a heuristic workaround: when the evolution stops as the time step goes below an acceptable threshold, we always find that max x;y ρ -> 0 as well. (The expansion of ingoing null cones, which is negative in flat spacetime, has become positive at least one point, presumably because the spacetime and our coordinate null cones have become highly nonspherical.) We take this combination as a (completely heuristic) indicator of collapse, even if C ¼ 0.99 has not yet been reached. At the largest nonsphericity and highest resolution, this collapse criterion also fails, and we replace it by max x Cðu; xÞ ≥ 0.8.
We cannot be sure that what we thus classify as collapse or dispersion really is, but if we then see subcritical curvature scaling on the "dispersion" side, and approximate self-similarity on both sides, we take this to be a strong indication that we really are bisecting to the collapse threshold. Furthermore, any masslike quantity in near-critical evolutions will also scale, in particular the Hawking mass on the first surface S u;x where C ¼ 0.99 (or C ¼ 0.8), and seeing this scaling at least provides further confirmation that we are bisecting to the black hole threshold, even if the mass we measure is related to the black hole mass only by a factor of order one.
Once we have an estimate for p à (which depends on all the numerical parameters), we evolve super and subcritical values of p that are equally spaced in ln jpp à j, with 30 values per power of 10, to give us better-resolved plots of the mass and curvature scaling laws.
We begin with the spherical case ϵ 2 ¼ 0, and run with N y ¼ 1 in lsB4 gauge. By experimentation we find that x 0 ¼ 8.24 allows bisection down to machine precision, while keeping the first appearance of C ≥ 0.99 somewhere in the middle of the grid as the bisection proceeds (in order to resolve the critical solution). In other words, the past light cone of the accumulation point of scaling echoes is at x à ≃ 8.24. We also find that x max ¼ 11 is required in order to capture the location where C first reaches the threshold value of 0.99 that we take as an indication of collapse. We find p à ≃ 0.2402.
We have checked that our code shows clean pointwise second-order convergence with Δx in the two evolutions forming our initial bracket: the clearly subcritical p ¼ 0.2, and the clearly supercritical p ¼ 0.3 (up to a little before our collapse criterion stops the supercritical evolution).
As a basic check, we have evolved spherically symmetric initial data also with Ny ¼ 3. This makes no visible difference to the scaling plots or critical solution (except where they dissolve into round-off noise). However, round-off error in the spectral operations triggers unphysical random nonspherical perturbations. These are then smoothed out by the shrinking of the numerical grid, and then, once their x-dependence is resolved, continue to evolve under the continuum perturbation equations. The l ¼ 2 spherical harmonic components of ψ and f, which from now on we shall denote ψ 2 and f 2 , reach an amplitude ∼10 -12 and then decay in sub15, while they blow up in super15 during the final collapse. This blowup may be explained by the physical instability of gravitational collapse to nonspherical perturbations, even if these perturbations were triggered here by numerical error.
We stop the bisection to p à after 50 iterations, where we are essentially at machine precision. However, the scaling laws for the initial black hole mass and maximum Ricci curvature become noisy well before we reach machine precision, at about sub13 for the Ricci scaling (see Fig. 1) and at about super12 for the mass scaling (see Fig. 2).
In spherical symmetry and evolving with N y ¼ 1, the source of the randomness is presumably round-off error, becoming important in the evolutions already three or so orders of magnitude before we reach machine precision in p itself.
From Appendix A, which clarifies the derivation given in [17], we expect the scaling laws
and the mass M of the first MOTS detected on our time slicing. We consider T -1=2 because it has dimension length, like M, and we use the natural logarithm for all our plots.
Here Δ ≃ 3.44 and γ ≃ 0.374 are universal constants, and A and B are family-dependent constants. f T is a universal periodic function with period Δ=ð2γÞ ≃ 4.60 in lnðp à -pÞ. f M is also periodic with the same period, and is universal with respect to initial data in any time slicing (such as our null cone slicing) that is compatible with discrete self-similarity (from now on, DSS). However, f M does depend on the slicing, as well as on our collapse criterion.
Note that the Ricci scalar is jRicj ¼ 8πT, so ln jRicj and ln T differ only by a constant.
To look for these scaling laws, we plot ln M -ðA þ γ ln jpp à jÞ against ln ln jpp à j, and similarly for T. We then expect to see only the fine-structures f M and f T . To uniquely fix A for a given family, we define f T to have zero mean. To fix B modulo periodicity, we let the minima of -ð1=2Þ ln T coincide with the minima of sinðγ ln jp à -pjÞ. We find that for our family of spherical initial data A ≃ 1.113, B ¼ 0.592, γ ¼ 0.374 and Δ ¼ 3.44 provide a good fit to (15) for our family of spherically symmetric initial data.
Our numerical measurement of f T with these fitting parameters is shown in Fig. 1. We see accurate periodicity from sub2.5 to sub13. The function f T is well approximated by a fundamental sine wave of amplitude 0.147, except for a widened top and a clear discontinuity in the derivative at the left edge of this flattened top. This kink is expected: it corresponds to the appearance of a new local maximum of T in each scale period.
Our numerical measurement of f M is shown in Fig. 2. Recall that M is the Hawking mass of the first surface S u;x with Hawking compactness C ¼ 0.99. f M is periodic with an amplitude of about 0.026 and mean -1.4835. There is one large discontinuity in each period, corresponding to the formation of the first MOTS one scale period later. Again this is expected: when a new local maximum of C takes over, by definition it has the same C, but not the same M. We see accurate periodicity with the same γ and Δ as for the curvature scaling law, over a slightly smaller range from sub4 to sub12. The range may be smaller because the MOTS occurs near the past light cone of the singularity of the underlying critical solution, whereas the maximum of the Ricci curvature occurs at the center, which may be less affected by details of the initial data.
We are aware of the following previous plots of scaling law fine-structures in the critical collapse of a spherical scalar field in the literature.
Pürrer, Husa and Aichelburg [18] extend compactified Bondi coordinates to future null infinity, and so can read off the asymptotic black hole mass. They show a periodic fine structure of this mass which appears to be continuous, with an amplitude of about 0.4 in ln M. Crespo, Oliveira and Winicour [19] use an affine radial coordinate also compactified to null infinity and show a similar fine structure, but with an amplitude that seems to be closer to 0.3 in ln M. These results are in the same time slicing, so they should agree, and they roughly do. From a plot in Rinne [15] we estimate the amplitude in ln M as ≃0.4, similar again to [18,19], although this is not the asymptotic mass. By contrast, our f M has an amplitude of only 0.026 in ln M, and has a jump of similar amplitude. Our best guess is that the Hawking mass of our collapse diagnostic is still very far from the final black hole mass. We have explained why our f M is discontinuous, but by contrast we have no theoretical understanding of why the MðpÞ measured by the other researchers is continuous.
Switching now to the subcritical scaling of the maximum of the Ricci curvature, Garfinkle and Duncan [20] plot the fine-structure in ln max jRicj, and from this plot we estimate the amplitude of f T as 0.3 in ln max jRicj, or 0.15 in ð1=2Þ ln max jTj. This is consistent with our value of 0.147. The bottom of their curve, corresponding to the top of ours, is flattened, again consistent with our observation. Baumgarte [5] has a fine-structure with amplitude of about 0.35 in ln max ρ, where ρ ¼ n a n b T ab is the energy density measured by an observer normal to the time slicing. This depends on the slicing, but the amplitude of 0.35 is again similar to the 0.3 observed by us and [20] for ln max jRicj.
We demonstrate self-similarity of the critical solution by showing that the scalar field ψ, the compactness C and the scaled curvature scalar T ≔ ðu à -uÞ 2 j∇ a ψ∇ a ψj ð 18Þ are periodic in the DSS-adapted coordinates
The null slicing is already DSS-compatible, in the sense of [21], but we call this specific coordinate system DSSadapted as the parameter u à needs to be adjusted to the accumulation point of the specific DSS. Unlike the scalar field, the metric actually has period Δ=2 because the Choptuik solution has the property that
The parameter u à is a family-dependent constant. A rough initial guess of u à is given by the time that near-critical supercritical evolutions stop. We then refine it by making ψ, C and T as periodic as possible. For our family, we find u à ≃ 5.609. Note that adjusting u à slightly will affect τðuÞ the more the closer u gets to u à . In practice, we can therefore adjust the τ-location of the last echo of, say, ψðx; τÞ by adjusting u à , almost without moving the other echoes. By contrast, this arbitrariness is absent when we use u à to make the amplitude of the last echo of T ðx; τÞ equal to that of the preceding echoes, and is therefore what we use to determine u à .
We use our closest super-and subcritical evolutions as a proxy for the critical solution. They agree with each other until one disperses and the other one forms a black hole. We see clear echoing in ψ, C and T for 2 ≲ τ ≲ 10. The results are shown in Figs. 3, 4, and 5 for our closest subcritical evolution.
We have also bisected ϵ 2 ¼ 0, N y ¼ 1 in "global shifted Bondi" (gsB) gauge. Recall from Paper I that this is given by R ¼ sðuÞx, giving
and we have settled on the particular flavor given by
We needed x max ¼ 20 in gsB gauge (instead of 11 in lsB gauge), and correspondingly N x ¼ 800 to maintain Δx ¼ 0.025. Comparing the zoomed-in scaling laws with lsB and gsB, they agree as expected down to about super/sub12, but then begin to disagree as both become noisy. FIG. 3. ϵ 2 ¼ 0: Surface plot of the scalar field ψðx; τÞ against the similarity coordinates ξ and τ, in our closest subcritical evolution (sub15). The plot has been cropped to 0 ≤ ξ ≤ 1, but the entire range of τ is shown, starting from the initial data at u ¼ 0 at the right edge of the plot, and ending at the left edge when C ≤ 0.01 indicates dispersion. The center is at ξ ¼ 0 (front edge), and the past light cone of the singularity at ξ ≃ 1 (back edge). Approximate DSS is seen for the range 2 ≲ τ ≲ 10. Our numerical data form a mesh made up of lines of constant u (constant τ) and constant x (not constant ξ). Both have been down-sampled for visual clarity. FIG. 4. ϵ 2 ¼ 0: Surface plot of the compactness C, otherwise as described in Fig. 3. FIG. 5. ϵ 2 ¼ 0: Surface plot of the scaled curvature diagnostic T , otherwise as in Fig. 3, except that the plot has been rotated slightly to the right for a clearer view.
In spite of this apparent success, we see a fundamental problem with gsB gauge already here in spherical scalar field collapse. While B is a monotonically decreasing function in weakly curved spacetimes, going through B ¼ 1 at x ¼ 0 and B ¼ 0 at x ¼ x 0 , in near-critical spacetimes it becomes nonmonotonic and may cross zero at x ¼ x 0 with positive slope. This means that some surfaces of constant x > x 0 are no longer future spacelike but timelike. This could include the outer boundary, and we would then not evolve on the domain of dependence. It is possible that B becomes negative again for sufficiently large x, and indeed this happens for x max ¼ 20 in this family of spherical initial data, but it is not something we can rely on.
Given that the only free parameter of gsB gauge is x 0 , and that this needs to be set to x à , which in turn is set by the family of initial data, in order to resolve the critical solution this seems to be a fundamental problem with gsB gauge. Indeed, gsB gauge fails for this reason in nonspherical evolutions, and we present no results in this gauge.
As already noted, the strong-field convergence tests with ϵ 2 ¼ 0.75 in Appendix B showed an initial burst of nonconvergent error caused by our initial data not being analytic at the origin, but we see second-order convergence with Δx and l max for all u ≳ 0.2. Here we present results for these nonanalytic initial data, as we believe the initial error does not affect our results in the critical regime.
We have bisected in the families with ϵ 2 ¼ 0, 10 -4 , 10 -2 , 0.1, 0.5 and 0.75, using lsB4 gauge. Note our family of null data parametrized by ðp; ϵ 2 Þ is equivalent to the family of Cauchy data parametrized by ðp; ϵ 2 Þ of [4,5] only in the limit of small ϵ 2 and small p, as we rely on a spherically symmetric solution of the wave equation on flat spacetime to relate the scalar ψ field at t ¼ 0 to u ¼ 0. Table I lists the numerical parameters and outcomes of our successful bisections. To estimate the error from discretization in x and y, we have run ϵ 2 ¼ 0.75 at four resolutions, with Δx ¼ 0.05 and 0.025, and Ny ¼ 9 and 17 for both. However, because of the large number of long evolutions involved, we have not plotted the fine structures of the scaling laws for the combination of the two higher resolutions. Moreover, at higher fine-tuning during the bisection, we had to reduce our criterion for collapse from C ≥ 0.99 to C ≥ 0.8. This would of course affect the mass scaling law (if we computed it), but appears to be sufficient for the bisection, and hence finding the approximate critical solution.
From (11), we expect that for sufficiently small ϵ 2 , the nonsphericity evolves as a linear perturbation of spherical symmetry, dominated by l ¼ 2, and that this linear perturbation decays. For somewhat larger ϵ 2 , [4,5] also found changes to Δ and γ. For even larger ϵ 2 , in particular ϵ 2 ¼ 0.75, Baumgarte [5] also found evidence of the (nonlinear) instability of the Choptuik solution reported by Choptuik et al [4] to lead to two centers of collapse. We will see, by contrast, that for our ϵ 2 ¼ 0.75 evolutions the nonspherical perturbations appear to still be in the linear regime, and indeed their nonsphericity is more similar to the evolutions of Baumgarte [5] with his ϵ 2 ¼ 0.5.
Our best supercritical evolution remains, by definition, close to our best subcritical one until the end when one decides to disperse and the other to collapse. However, in its final collapse phase the best supercritical solution becomes much more nonspherical, and this is challenging numerically.
In flat spacetime ΞR ¼ -1=2, Rρ þ ¼ 1 and 2Rρ -¼ -1. By contrast, in critical collapse, the divergence ρ þ of the outgoing null geodesics remains approximately spherical and positive in the self-similar phase, but in the final collapse becomes very small at large x. Similarly, ρ -and TABLE I. Table showing parameters of our families of initial data. ϵ 2 , in the first column, is the only physical input parameter. The second group of columns shows numerical (input) parameters, and the third group physical (observed) parameters (fitted to the numerical data, and rounded to four significant digits for this table). These families, or the best subcritical evolutions in them, have been used in the figures listed in the last column. For visual clarity, a * means the entry is unchanged from the previous row. Half-frequency filtering was applied except in the first row, and the collapse criterion was C ≥ 0.99 except in the last row, where it was C ≥ 0.8.
.24 11. 0.2402 5.609 1.113 0.592 No filtering, Figs. 1-7, 12-13 10 -4 * 9 8 * * * * * * Figs. 12-13 10 -2 * 3=5=9 2=4=8 8.275 15. * 5.610 1.115 * Figs. 6-9, 12-15, 18-19, 22 0.1 * 3=5=9 2=4=8 8.7 * 0.2406 5.620 1.131 * Figs. 6-7, 12-15, 18-19 0.5 * 5=9 4=8 12.5=13.3 * 0.2452 5.657 1.230 0.737 * 0.75 0.05 9 8 24. 50. 0.2555 5.642 1.250 1.059 Figs. 10-11 * * 17 16 30. * 0.2556 * * * * * 0.025 9 8 24.5 * 0.2535 5.651 * * * * * 17 16 33. * * 5.652 * * C ≥ 0.8, Figs. 6-7, 10-21, 23
ΞR remain negative and approximately spherical in the self-similar phase, but become very large, and positive for a range of angles intermediate between the poles and equator, in the final collapse phase. We also see a fundamental numerical challenge of highly nonspherical evolutions in lsB gauge: The x-shift B in
tries to counteract the nonspherical part of ΞR to keep R independent of y but has little to act on as R ;x → 0, while ΞR also becomes larger. Therefore B becomes large, and by the Courant condition applied to the upwinded shift terms, the time step becomes small. In spherical symmetry this is less of a problem because R ;x → 0 also indicates an approaching MOTS, whereas in a nonspherical spacetime it does not. Another consequence of the asphericity of the final collapse phase is that we do not then have sufficient angular resolution to obtain reliable values for M. In the continuum C ¼ C and M ¼ M, but in the final collapse phase their numerical values differ strongly. There is little hope that we can overcome this problem with more angular resolution, as physically we expect the S u;x to become arbitrarily nonspherical inside the black hole. As a heuristic indicator of collapse for the purpose of bisecting to the black hole threshold, we therefore use C rather than C, as M is at least nondecreasing by construction. Of course either M or M should be considered reliable only to the extent they agree with each other.
We look for three types of evidence that the nonspherical perturbations are in the linear regime: for sufficiently small ϵ 2 , the spherical part of the critical solution should be the same critical solution as for ϵ 2 ¼ 0, and its perturbations should decay and oscillate as predicted by linear perturbation theory about the spherical critical solution. The growth rate λ 0 > 0 of the unique growing mode should be unchanged, and hence the scaling laws, including γ, f T and f M should be as for ϵ 2 ¼ 0. Finally, the decay rates λ of the leading perturbation modes should be independent of ϵ 2 for ϵ 2 sufficiently small, with the values predicted in perturbation theory in [17]. We examine evidence for all this in the next two subsections, starting here with the scaling laws.
Figures 6 and 7 compare the fine structures of T -1=2 and M scaling, respectively, for different ϵ 2 , using our highest resolution data at each ϵ 2 . The same γ ¼ 0.374 power law has been taken out for all values of ϵ 2 . We are treating each value of ϵ 2 as a different 1-parameter family of initial data for the purpose of fitting A and B.
To estimate the accuracy of our scaling results for nonspherical initial data, Figs. 8 and 10 show f T at different resolutions for ϵ 2 ¼ 10 -2 and ϵ 2 ¼ 0.75, respectively. At ϵ 2 ¼ 10 -2 we have very good accuracy for f T . At ϵ 2 ¼ 0.75 we have some discretization error in the amplitude of f T , but essentially none in its period.
On the other hand, Figs. 9 and 11 show significant resolution dependence in the numerically measured f M already at ϵ 2 ¼ 10 -2 , and similar, but not much larger, at ϵ 2 ¼ 0.75. We believe the variation of f M with ϵ 2 shown in Fig. 7 is real but, we do not have enough numerical resolution for quantitatively reliable results.
Our plots of f T in Fig. 6 at different values of ϵ 2 can be compared directly with Fig. 3 of [4] and Fig. 7 of [5], with the differences that we plot -1=2 ln T þ A rather than ln T on the vertical axis, and lnðp à -pÞ þ B rather lnðp à -pÞ FIG. 6. Comparison of the curvature scaling functions f T for the ϵ 2 ¼ 0 (black line), ϵ 2 ¼ 10 -2 (blue line), ϵ 2 ¼ 0.1 (cyan line), ϵ 2 ¼ 0.5 (light green line) and ϵ 2 ¼ 0.75 (dark green circles) families. Each f T is computed at the highest numerical resolution shown in Table I. Here and in the following plots, the horizontal range is exactly the same in Figs. 1 and 2, and the distribution of data points is the same as indicated by circles there, but for clarity we no longer show those individual data points, except for ϵ 2 ¼ 0.75, where our computation at the highest resolution has gaps due to lack of computing time. FIG. 7. Comparison of the mass scaling function f M for the same values of ϵ 2 as in Fig. 6: ϵ 2 ¼ 0 (black), ϵ 2 ¼ 10 -2 (blue), ϵ 2 ¼ 0.1 (cyan), ϵ 2 ¼ 0.5 (light green) and ϵ 2 ¼ 0.75 (dark green), at the same resolutions as there.
on the horizontal axis. In other words, [4] and [5] do not fit our family-dependent parameters A and B, but in effect set them to zero.
Setting that aside, our results show only a tiny effect of ϵ 2 on γ, compared to [4,5,12]. In our most nonspherical evolutions, ϵ 2 ¼ 0.75, we only see a change of δγ ≃ -0.0016. By comparison, [4,5] find δγ ¼ -0.007 and -0.005, respectively, at their ϵ 2 ¼ 0.5, which appears to be the nearest comparator to our ϵ 2 ¼ 0.75 (see also Table II).
We also see only a tiny deviation of the period of f T from its theoretical value Δ=ð2γÞ (with Δ ¼ 3.44 and γ ¼ 0.374) for all values ϵ 2 ¼ 0…0.75. This is demonstrated most clearly in Figs. 12 and 13, which show ψð0; τÞ (technically, ψ 0 extrapolated to x ¼ 0, which is not on the grid). Our best fits are δ ln Δ ≃ 5 × 10 -4 , -8 × 10 -4 and -64 × 10 -4 at ϵ 2 ¼ 0.1, 0.5 and 0.75. These are not even monotonic in ϵ 2 and therefore more likely due to numerical errors than a physical change.
For more indirect evidence of Δ from the periodicity of fine structure of the scaling laws, see also Fig. 1 for a fit of f T at ϵ 2 ¼ 0 to a sine wave of period 3.44, and Fig. 6 for a comparison of f T with values from ϵ 2 ¼ 0 up to ϵ 2 ¼ 0.75.
We conclude that Δ does not change by more than 10 -3 up to our ϵ 2 ¼ 0.75. By comparison, [4,5] find δΔ ¼ -0.05… -0.12 at their ϵ 2 ¼ 0.5 (see also Table II).
Unfortunately, [4] and [5] do not present plots of f M . Our own results in Fig. 6 show that f M differs significantly from ϵ 2 ¼ 0 already at ϵ 2 ¼ 10 -2 , and has become very much larger at ϵ 2 ¼ 0.75. Recall however that our M is only a very rough indication of the final black hole mass, particularly so in our best resolution for ϵ 2 ¼ 0.75 where we have measured the area of a two-surface with Hawking FIG. 9. The mass scaling function f M for ϵ 2 ¼ 10 -2 , at the same four resolutions as in Fig. 8. The differences in f M between resolutions are very large, compared to f T in Fig. 8. FIG. 10. The curvature scaling function f T for ϵ 2 ¼ 0.75, at four different resolutions: Ny ¼ 17, Δx ¼ 0.005 (our best resolution, circles), Ny ¼ 9, Δx ¼ 0.0025, x 0 ¼ 24.5 (dashed lines), and Ny ¼ 17, Δx ¼ 0.005, x 0 ¼ 30 (dotted), and Ñy ¼ 9, Δx ¼ 0.005, x 0 ¼ 24 (dot-dashed). The red straight lines represent a fit by eye to the local minima and maxima of f T . They have slope δγ ¼ -0.0016. FIG. 8. The curvature scaling function f T for ϵ 2 ¼ 10 -2 (blue), at three different angular resolutions: Ny ¼ 3 (solid) Ny ¼ 5 (dashed) and Ny ¼ 9 (dotted), and l max ¼ 2, 4, 8, respectively, all with Δx ¼ 0.025, x 0 ¼ 8.275, x max ¼ 15. The curves are indistinguishable, indicating that the errors from discretizing in y are small. FIG. 11. The mass scaling function f M for ϵ 2 ¼ 0.75, at the same four resolutions as in Fig. 10. The red straight lines have slope δγ ¼ -0.0016, taken from the fit in Fig. 10. f M is both not periodic enough and too resolution-dependent for a meaningful fit to its average slope, but δγ ¼ -0.0016 is consistent with the data. compactness of only 0.8. We therefore do not believe our results for f M are accurate, in contrast to our results for f T .
We now test the hypothesis that there is a regime of small ϵ 2 where near-critical time evolutions can be approximated by the spherical critical solution plus small perturbations, and attempt to find the limit of its validity.
We begin with the spherical part of the scalar field and metric. Under the assumption that T is still dominated by the derivatives of ψ 0 , we adjust u à to make T as periodic as possible, and in particular with as many maxima and minima taking the same values, just as in the spherical case.
TABLE II. Comparison of families of initial data near the black hole threshold, combining data from our evolutions with those of [4,5] and [12]. Families corresponding to the same initial data have been grouped together. A dash means no data are available. A range given for max jδψj means that its value increases noticeably during the phase where the spherical part of the solution is approximated by the Choptuik solution. By contrast, the ranges given for δγ and δΔ express uncertainty, see the main text for details. Footnotes: (e) relative to exact values Δ ≃ 0.374 and Δ ≃ 3.44 in spherical symmetry; (n) relative to numerical value obtained in spherical symmetry; (i) initial and (f) final value in the approximately Choptuik phase; (1) first, (2) second and (3) third method of determining Δ; (3,n1), for example means, the numerical value of Δ in a nonspherical evolution obtained by the third method, relative to the numerical value determined by the first method in spherical symmetry.
Family max jδψj δγ δΔ Bif. Choptuikþ ϵ 2 ¼ 0 0 ðiÞ …0.05 ðfÞ 0.008 ðeÞ 0 ð1;eÞ …0.05 ð2;eÞ No Baumgarte ϵ 2 ¼ 0 0 0 ðeÞ 0.02 ð3;eÞ …0.03 ð1;eÞ No Maroudaþ I --0.004 ðeÞ -0.02 ð3;eÞ …0.01 ð1;eÞ No This paper ϵ 2 ¼ 0 0 0 ðeÞ 0 ðeÞ No Baumgarte ϵ 2 ¼ 10 -2 0.005 0 ðeÞ -0.02 ð3;n12Þ …0.03 ð1;eÞ No This paper ϵ 2 ¼ 10 -2 0.008 < 10 -3 < 10 -3 No Maroudaþ II 0.02 -0.002 ðeÞ …0.002 ðnÞ -0.07 ð3;eÞ …0.03 ð2;eÞ No Maroudaþ III 0.06 -0.010 ðeÞ … -0.006 ðnÞ -0.12 ð3;eÞ … -0.05 ð2;eÞ No This paper ϵ 2 ¼ 0.5 0.07 < 10 -3 < 10 -3 No Choptuikþ ϵ 2 ¼ 0.5 --0.007 ðnÞ …0.001 ðeÞ -0.12 ð1;n2Þ … -0.05 ð2;eÞ No Baumgarte ϵ 2 ¼ 0.5 0.08 -0.005 ðeÞ -0.12 ð3;n1Þ … -0.05 ð1;eÞ No Choptuikþ ϵ 2 ¼ 2=3 0.06 ðiÞ …0.12 ðfÞ -0.036 ðnÞ … -0.028 ðeÞ -0.41 ð2;n2Þ -0.31 ð1;eÞ No This paper ϵ 2 ¼ 0.75 0.09 -0.0016 ðeÞ -0.023 No Maroudaþ IV 0.1 ðiÞ …0.2 ðfÞ -0.045 ðeÞ … -0.041 ðnÞ -0.49 ð1;eÞ … -0.35 ð2;eÞ Yes Choptuikþ ϵ 2 ¼ 0.75 --0.069 ðeÞ … -0.061 ðnÞ -0.62 ð1;n2Þ … -0.41 ð2;eÞ Yes Baumgarte ϵ 2 ¼ 0.75 0.11 ðiÞ …0.34 ðfÞ -0.068 ðeÞ -0.72 ð3;n1Þ … -0.57 ð1;eÞ Yes Choptuikþ ϵ 2 ¼ 5=6 --0.102 ðnÞ … -0.094 ðeÞ -2.49 ð2;n2Þ … -1.44 ð1;eÞ Yes FIG.
12. The scalar field at the center, ψ 0 ð0; τÞ, at our best numerical resolution and best subcritical value of p, for each of ϵ 2 ¼ 0 (black), 10 -4 (purple), 10 -2 (blue), 0.1 (cyan), 0.5 (lightgreen) and 0.75 (dark-green). The horizontal axis shows τ and the vertical axis ψ at the center.
FIG. 13. The same plot as in Fig. 12, but a linear transformation applied to τ applied to ϵ 2 ¼ 0.1, 0.5 and 0.75 (but not the smaller values), in order to align the first and third full maxima. For this, τ is stretched by factors of 0.9995, 1.0008 and 1.0068, respectively, as well as shifted.
We then find that ψ 0 ðx; τÞ, and Cðx; τÞ, with the same fitted value of u à , are essentially identical from ϵ 2 ¼ 0 through to ϵ 2 ¼ 0.75. See again the limit on a possible change of Δ of the critical solution documented above. We next address the hypothesis that the deviations from spherical symmetry are small and essentially linear.
To find the linear perturbations of a given spherical DSS critical solution ϕ Ã , one can separate the linearized equations by angular dependence l, and for each l and dimensionless field ϕ make a mode ansatz δϕ l ðξ; τÞ ¼ e λτ δ φlλ ðξ; τÞ; ð24Þ with δ φlλ defined to be periodic in τ with period Δ, and regular at the center and at the past light cone. Here ϕ stands for any scale-invariant quantity, such as ψ or T . From this ansatz, the complex mode function δ φlλ and corresponding complex Lyapunov exponent
are then determined as the eigenfunctions and eigenvalues of a (singular) linear boundary value problem [3].
As the background solution and its perturbations are real (the complex mode ansatz is only for convenience), λ and the corresponding δ φlλ are either real or form complex conjugate pairs. From these complex modes, one can construct the corresponding real perturbations as δϕ l ðξ; τÞ ¼ Re  Ae iα e ðκþiωÞτ δ φlλ ðξ; τÞ Ã
for arbitrary positive real amplitude A and phase 0 ≤ α < 2π. (Note this is slightly incorrect in Eq. (25) of Ref. [5]). Unless ω a rational multiple of 2π=Δ, the product e -κτ δϕ l ðξ; τÞ, while not growing or decaying, is nevertheless not periodic in τ but only quasiperiodic. The perturbation modes of the spherical scalar field critical solution were found numerically in [3], using different similarity coordinates from the ones defined here. We can therefore compare the λ with [3], but not directly the δ φlλ ðτ; xÞ. As expected for a critical solution, [3] finds a single growing l ¼ 0 mode, with λ and the corresponding mode function real. All other spherical modes, and all nonspherical modes, are complex and decay, that is, κ < 0. The least damped (most slowly decaying) nonspherical mode was found to have l ¼ 2 angular dependence, with κ ≃ -0.07=Δ and ω=ð2πÞ ≃ 0.3=Δ.
We expect that initial data which are almost spherical and fine-tuned to the threshold of collapse, but otherwise generic, evolve into something that can be approximated by the spherical critical solution plus a linear perturbation, and that the perturbation can be represented as sum of modes each with its own complex amplitude Ae iα determined by the initial data.
During the intermediate range of τ where the solution is approximated by the critical solution plus small perturbations, we expect the least damped l ¼ 2 perturbation to dominate the l ¼ 2 component of the solution. The singlemode formula (26) should then approximately describe the l ¼ 2 component of near-critical evolutions in this regime. Moreover, its amplitude should be proportional to ϵ 2 , as the l ¼ 2 component of the initial data is proportional to ϵ 2 to leading order.
In our initial data, the nonspherical part of ψ is, to leading order in ϵ 2 , purely l ¼ 2, and so, going to quadratic order, we expect the initial data for the l ¼ 4 spherical harmonic component ψ 4 to be proportional to ϵ 2 2 . During the time evolution ψ 4 is sourced by terms that are linear in f 4 or quadratic in f 2 and ψ 2 . f 4 is initially zero and then sourced by terms linear in ψ 4 and quadratic in f 2 and ψ 2 . Perturbatively in ϵ 2 , we therefore expect all l ¼ 2 components to be proportional to ϵ 2 e κτ , all l ¼ 4 components to be proportional to ϵ 2 2 e 2κτ , and so on for higher l. To demonstrate the expected scaling with both τ and ϵ 2 , in Fig. 14 we plot the maximum and minimum of ϵ -1 2 e -κτ ψ 2 ðx; τÞ over x in the range 0 ≤ x ≤ 1 against τ, for different ϵ 2 , with our best subcritical p at our best numerical resolution. The same is done in Fig. 15 for f 2 . We find that ϵ -1 2 e -κτ ψ 2 is essentially the same at ϵ 2 ¼ 10 -4 , 10 -2 and 0.1, and so is ϵ -1 2 e -κτ f 2 . At ϵ 2 ¼ 0.5 and 0.75, they are a bit larger (by a factor of ∼1.5 at ϵ 2 ¼ 0.75), but still agree qualitatively feature by feature. In other words, ψ 2 and f 2 are almost linear up to ϵ 2 ¼ 0.1, and still approximately linear up to 0.75.
In these plots, we have fitted κ ¼ -0.03 by eye. Compare this with the value of κ ≃ -0.07=Δ ≃ -0.02 determined in [3] from a numerical construction of the eigenvalues of the operator evolving a linear perturbation from τ to τ þ Δ. Neither is very accurate, but we believe that our fitted value FIG. min x and max x of ϵ -1 2 e -κτ ψ 2 ðx; τÞ, plotted against τ, for ϵ 2 ¼ 10 -4 (purple), ϵ 2 ¼ 10 -2 (blue), ϵ 2 ¼ 0.1 (cyan), ϵ 2 ¼ 0.5 (light green) and ϵ 2 ¼ 0.75 (dark green).
is compatible with the theoretically expected one, but not with zero. In other words, we find that ψ 2 and f 2 decay as expected from the linear perturbation mode analysis, up to ϵ 2 ¼ 0.75.
Having established these scalings, we give the full surface plots for ψ 2 ðx; τÞ and also f 2 ðx; τÞ, for ϵ 2 ¼ 0.75 in Figs. 16 and 17.
The l ¼ 4 modes of ψ and f in the best subcritical evolutions also seem to be quasiperiodic in τ and roughly independent of ϵ 2 when rescaled with ϵ -2 2 e -2κτ . This is demonstrated in Figs. 18-21. FIG. 17. Surface plot of the scaled metric component ϵ -1 2 e -κτ f 2 ðx; τÞ, otherwise as in Fig. 16.
x and max x of ϵ -2 2 e -2κτ ψ 4 ðx; τÞ, plotted against τ, otherwise as in Fig. 14. FIG. 15. min x and max x of ϵ -1 2 e -κτ f 2 ðx; τÞ, plotted against τ, otherwise as in Fig. 14. FIG. 19. min x and max x of ϵ -2 2 e -2κτ f 4 ðx; τÞ, plotted against τ, otherwise as in Fig. 14. Where the minimum is exactly zero, f 4 ≥ 0, with equality only at the origin. FIG. 16. Surface plot of the scaled nonspherical scalar field component ϵ -1 2 e -κτ ψ 2 ðx; τÞ, in the sub15 evolution of the ϵ 2 ¼ 0.75 family. Otherwise as described in Fig. 3. Note that the dark green curves in Fig. 14 FIG. 20. Surface plot of the scaled nonspherical scalar field component ϵ -2 2 e -2κτ ψ 4 ðx; τÞ against the similarity coordinates ξ and τ, in the sub15 evolution of the ϵ 2 ¼ 0.75 family.
The mode frequency ω could in principle be obtained by a (discrete) Fourier transform of our data in τ. We have not attempted this as the region where the background spherical solution is clearly DSS and the different perturbation amplitudes agree (after rescaling) is only 1 ≲ τ ≲ 9, giving us only about two periods of the background solution. The discrete set of angular frequencies present in the Fourier transform with respect to τ of (26) is
where the cited ω is the value for N ¼ 0, that is, the smallest positive frequency in the spectrum. Because of the symmetry ψðx; τ þ Δ=2Þ ¼ -ψðx; τÞ, N takes odd values for the scalar field and even values for the metric. The spectrum obviously depends on the value of the parameter α in (26), but for a random α, [3] found that the peak of the spectrum was at N ¼ 5. The corresponding period in τ is
If we count the peaks in Fig. 14, we find about 12 in the range 0 ≤ τ ≤ 12, giving us a "pseudoperiod" of about 1, similar to P 3 ≃ 1.0. If we count the peaks in Fig. 15, with find about 8 in the same range, giving us a pseudoperiod of about 2=3, similar to P 5 ≃ 0.65. The numerical measurement of both κ and ω would be improved in proportion to increasing the range of τ over which our fine-tuned numerical solutions approximate the critical solution. Going any further in this would require quadruple precision, not only so that p can be represented to more significant figures, but more importantly so that round-off error during the numerical solution is suppressed better.
We note in passing that ψ 2 =x 2 , f 2 =x 2 and f 4 =x 4 , which should be regular at x ¼ 0 in the continuum, are also numerically regular at x ¼ 0, but ψ 4 =x 4 is not. (It is generally hard to enforce f l ∼ r l in any freely evolved variable, without explicitly taking out the factor r l .)
Baumgarte, and Marouda et al.
In order to compare the nonsphericity of our evolutions with those of [5] (and by implication also of [4] for the same initial data), we note that Baumgarte has measured the difference in ψ on outgoing null geodesics at y ¼ AE1 and y ¼ 0 (poles and equator). Fig. 11 of [5] shows this δψ for ϵ 2 ¼ 10 -2 during the approximately self-similar phase, plotted against the similarity variables τ ¼lnðu à -uÞ and ξ λ ≔ λ=ðu à -uÞ, where λ is the affine parameter, normalized to λ ¼ 0 and dλ=dR ¼ 1 at the origin, and u à is fitted as described above. This measure is gaugeinvariant.
We note that for infinitesimal deviations from spherical symmmetry, and assuming the deviation only has an l ¼ 2 component, δψðu; xÞ ≔ ψðu; x; AE1Þψðu; x; 0Þ ð 29Þ
where ψ 2 denotes the l ¼ 2 component of ψ. λ along outgoing geodesics is given by
The equivalent of Fig. 11 of [5] created from our data for ϵ 2 ¼ 10 -2 is Fig. 22. Note that in this plot we have not applied the factor of e -κτ . (This makes little difference, as κ is so small.) At small ϵ 2 our solutions are a little more nonspherical than those of [5], but recall that we set initial data on different slices.
Figure 23 shows δψ, approximated as ð3=2Þψ 2 for our ϵ 2 ¼ 0.75. For comparison, Figs. 24 and 25 show plots of δψ created from the simulations of [5] for ϵ 2 ¼ 0.5 and 0.75, respectively. This suggests that, during the approximately self-similar phase, our ϵ 2 ¼ 0.75 family is about as nonspherical as ϵ 2 ¼ 0.5 of [5]. In both families δψ does not grow noticeably. By contrast, Fig. 25 shows clear growth of δψ with τ for ϵ 2 ¼ 0.75 and an end to the perturbative regime around τ ¼ 6.
Table II compares the black hole threshold in selected families taken from [4,5,12] and the present paper. Families that should be directly comparable are grouped together. For comparison have translated δψ=ϵ 2 and ψ 2 into δψ (by multiplying by ϵ 2 and 3=2, respectively). The values for max jδψj from the null code and that of [5] are read off from plots shown in the present paper. The values for ϵ 2 ¼ 0 and ϵ 2 ¼ 2=3 of [4] are read off from Figs. 6 and 7 of that paper, respectively. In the spherically symmetric case we include max jδψj as an estimate of the numerical error generated by evolving a spherically symmetric solution in cylindrical coordinates. Finally, the values for families II, III and IV of [12] are read off from plots communicated by the authors.
A range given for max jδψj means that its value increases noticeably from the beginning to the end of the phase where the spherical part of the solution is approximated by the Choptuik solution. A single value implies that the amplitude does not grow noticeably. (This is estimated by eye, as the perturbations are only quasiperiodic.)
By contrast, the ranges given for δγ and δΔ express uncertainty: in the four papers, up to three different methods have been used to estimate Δ (but only one for γ). Independently, one can use either the numerical values of γ and Δ obtained in spherical symmetry in each paper, or their known exact values as reference values. We cite here the largest and smallest of these (up to six) possible differences, in order to indicate a plausible interval for δγ and δΔ for each family. We do not include the fitting errors given in [4,12], as these appear to be smaller than the systematic errors suggested by the intervals just mentioned.
We have attempted to order the families across papers, taking as our first ordering criterion the bifurcation (yes or no) of the critical solution, as our second criterion max jδψj (where we have data), and otherwise δγ and δΔ.
The fact that we are able to order the families of initial data consistently (within the estimated plausible ranges) is one of the main physics results of this paper. It supports the hypothesis that all near-critical evolutions go through a phase where the solution can be approximated by the Choptuik solution plus the growing spherical mode and the least damped l ¼ 2 mode. Large nonsphericity seems to FIG. 23. As in Fig. 22, but for ϵ 2 ¼ 0.75. FIG. 24. Contour plot of δψ=ϵ 2 against τ and x λ from the best subcritical ϵ 2 ¼ 0.5 evolution of Baumgarte [5]. This is the ϵ 2 ¼ 0.5 equivalent of Fig. 11 of [5]. ðτ; x λ Þ are called ðT 0 ; λÞ in this plot. We have chosen a contour plot rather than a surface plot as this is clearer, given the coarse resolution of the available data. FIG. 22. Surface plot of ð3=2Þ=ϵ 2 ψ 2 ≃ δψ=ϵ 2 for ϵ 2 ¼ 10 -2 against similarity coordinates τ and ξ λ . This should be compared to Fig. 11 of [5]. The color map has been chosen purely for visual agreement with that figure. The range of ξ λ is the same as there, but our range of τ is larger, [0, 14] rather than [0, 8], as we can fine-tune more closely to the black hole threshold.
change the values of Δ and γ, as well as the decay rate κ of the l ¼ 2 mode, such that for nonsphericities of max jδψj ≳ 0.1 it actually grows. (Compare this with max jψ Ã j ≃ 0.6 for the critical solution.)
The exception from this consistent picture are the values of δγ and δΔ obtained in this paper. Up to our ϵ 2 ¼ 0.5 family, these are zero within our estimated numerical errors. For our ϵ 2 ¼ 0.75 family, we have a signifant but small change of δΔ ∼ 10 -3 (but not of δγ), but this is still an order of magnitude smaller than δΔ observed in the other three papers at what we take to be comparable nonsphericity. For of a better explanation we suspect that this is due to some systematic numerical error, more likely in the null code than the other three codes.
We have already mentioned that to complete the bisection in p for ϵ 2 ¼ 0.75 at our highest resolution Ny ¼ 17, Δx ¼ 0.025, we had to lower our heuristic diagnostic of collapse to C ≥ 0.8.
We have tried to bisect at ϵ 2 ¼ 0.79, with Ny ¼ 17, Δx ¼ 0.05, and again using C ≥ 0.8 as the collapse diagnostic, starting again with the initial bracket 0.2 < p à < 0.3. Already at the third bisection step, the code stops because B becomes very large at the boundary, and so the time step becomes very small. On closer inspection, B becomes very irregular at late times at the outermost few points, and very much larger at the outermost grid point than anywhere else. This cannot be a boundary instability in the usual sense because the outer boundary is treated exactly like an ordinary grid point. Hence the instability must be one that grows much faster at larger x, and so most rapidly at the outer boundary.
We believe what is at fault here is a combination of two problems we have already mentioned. One is that in our gauge, we evaluate min y ΞRðu; x; yÞ, and this means that the B ;x becomes discontinuous. This lack of smoothness then propagates to the other variables. A separate problem is that R ;x becomes small as outgoing light cones are trying to recollapse, while lsB gauge (or any gauge which solves the Raychaudhuri equation along null generators for G rather than for R) requires R ;x > 0. We believe this would still give rise to very large values of B even if we found a better gauge within the lsB family that avoided the unsmoothness.
This work was born out of an attempt to generalize the method of Garfinkle [13] for simulating critical collapse without the need for adaptive mesh refinement beyond spherical symmetry. We have demonstrated that this is possible, in the example of axisymmetric scalar field critical collapse.
We have attempted to duplicate previous physics results as well. Recall that the collapse simulations of Baumgarte (in spherical polar coordinates) [5] seemed to have reconciled the linear perturbation results of Martín-García and Gundlach [3] with the collapse simulations of Choptuik et al. [4] (in cylindrical coordinates): small deviations from spherical symmetry decay when one fine-tunes to the threshold of collapse, while large perturbations grow and lead to the formation of two centers of collapse. This conclusion is also supported by the more recent evolutions of [12] and (less quantitatively) of [11].
Specifically, [5] and [4] find that nonsphericities decay for values of their nonsphericity parameter ϵ 2 ≤ 0.5, and [4] still find this at ϵ 2 ¼ 2=3. Both [5] and [4] find a bifurcation for ϵ 2 ¼ 0.75, and [4] also at ϵ 2 ¼ 5=6. (We note that [12] also observe bifurcation in their most nonspherical data for the complex scalar field). [4,5] and [12] also found that Δ and γ (measured before any bifurcation occurs) decrease with increasing ϵ 2 , see Table I for representative numerical values.
By contrast, we have found that even in the evolution of our ϵ 2 ¼ 0.75 family of data perturbation theory is a good model for the nonsphericity. In particular, the echoing period and critical exponent remain at Δ ≃ 3.44 and γ ≃ 0.374, their values in spherical symmetry, and the l ¼ 2 nonspherical components of all fields decay at roughly the rate predicted for linear perturbations in [3]. We also find that the amplitude of the l ¼ 2 field components is almost linear in ϵ 2 for values between 10 -4 and 0.75, enhanced only by a factor ∼1.5 at ϵ 2 ¼ 0.75. Similarly, the l ¼ 4 nonsphericity scales with ðϵ 2 Þ 2 , as one would expect from perturbation theory to quadratic order, given that l ¼ 4 is absent in the initial data to linear order in ϵ 2 , and that [3] predicts l ¼ 4 to decay much more rapidly than l ¼ 2.
However, we cannot evolve the same family of initial data as [5] and [4], because we set data on an outgoing null cone. We have compared a gauge-invariant measure of the nonsphericity of the scalar field, in the phase where a nearcritical evolution is approximately self-similar, with the same measure for the evolutions of [5], and find that this comparable between our ϵ 2 ¼ 0.5 and ϵ 2 ¼ 0.75 of [5] (and by implication of [4]). If we compare our ϵ 2 ¼ 0.5 with ϵ 2 ¼ 0.75 of [4,5], a mild tension remains, as these authors claim that Δ ≃ 3.44 is reduced by ∼0.08 and γ ≃ 0.374 by ∼0.007, whereas we see much smaller reductions. This could potentially still be explained by numerical error.
In order to accurately measure the critical exponent γ, we have fitted not only the power laws with exponent γ, but also the small periodic fine-structure superimposed on them. This has only been done in spherical symmetry before. Our fine-structure of the curvature scaling laws (on the subcritical side) agrees well with published results in spherical symmetry, but our fine structure of the mass scaling law does not. The reason may be that our collapse criterion is not the first appearance of a trapped surface (in a given time slicing) but the first appearance of a coordinate 2-surface S u;x with Hawking above a threshold value, and we evaluate the Hawking mass of that surface.
We believe this is the first simulation of gravitational collapse beyond spherical symmetry in null coordinates. (The paper [22] investigates axisymmetric supernova core collapse, but not black hole formation.) We have applied our methods to the challenging problem of axisymmetric scalar field critical collapse, and at moderate nonsphericity have been able to fine-tune our initial data to the black-hole threshold to machine precision, without the need for mesh refinement.
We have not yet been able to simulate critical collapse for large enough nonsphericity to fully duplicate the results of [4,5,12]: at large nonsphericity, our evolutions stop before we can classify them as either forming a black hole or dispersing. The proximate cause for this is that the divergence of the null generators of our coordinate light cones is trying to become negative at some points on our last null cone, but the generalized Bondi coordinates we use break down if this happens. Therefore, the next step will be to change to a generalized affine radial coordinate. This will then allow us to investigate if our null cones themselves remain regular in critical collapse, or form caustics.
(compared to C ≃ 0.6 in the critical solution and C ¼ 1 on a black-hole horizon). We run to u ¼ 6, when the solution largely dispersed, with C < 0.01, but the numerical domain has contracted from R max ¼ 25 at the initial time only to R max ≃ 14, thus avoiding very small time steps.
The evolutions with p ¼ 0.3 start from C ≃ 0.18 and reach our heuristic collapse criterion C ¼ 0.8 at u ≃ 3.9, at which point we stop the evolution.
Both p ¼ 0.2 and p ¼ 0.3 are far enough from p à ≃ 0.25 that it remains meaningful to compare numerical solutions with different resolutions at the same coordinate time u and amplitude p, as jp Ãpj ≃ 0.05 is then much larger than the variation of p à with different numerical parameters. By contrast, in the critical regime we would need to compare different resolutions at the same (small) p=p à -1, with p à resolution-dependent, and plot them against the similarity-adapted coordinates ðτ; ξÞ, which also depend sensitively on numerical parameters through the fitted value of u à .
We output all fields as l-components (rather than against y) to make comparisons at different Ny and l max simpler. We have tested for self-convergence to second order in Δx and in l max , see Paper I for details of how the errors presented here are defined.
We have evolved at two intersecting families of resolutions: at Δx ¼ f0.1; 0.05; 0.025; 0.0125g with fixed angular resolution Ny ¼ 17, l max ¼ 16; and at Ny ¼ f9; 13; 17; 25; 33; 49; 65g, l max ¼ Ny -1, with fixed radial resolution Δx ¼ 0.05. Their intersection is our baseline solution, chosen because the discretization errors in x and y are roughly similar. Note that our adaptive timestep Δu is approximately proportional to Δx and approximately independent of angular resolution.
We find a large and nonconverging error at small u and x, in all ψ l with l > 2. We believe this is due to the fact that our initial data ( 11)-( 12) is not actually single-valued at the origin. The nonconvergent error disappears quickly, and we believe this is because the boundary conditions at the center force the solution to become single-valued. At all other ðu; xÞ we find the expected second-order convergence in Δx. For a spectral method acting on a smooth solution, we would expect approximately exponential convergence, but in fact the discretization error in y decreases with resolution only as l -2 max . Second-order convergence in the L 2 norm with Δx and 1=l max is demonstrated for ψ 0 in Fig. 26, for ψ 4 in Fig. 27, and for C in Fig. 28. The convergence is actually pointwise. The global maximal errors max u jjE x jj L 2 and max u jjE y jj L 2 for the dispersing are shown in Table IV. Given second-order convergence, we can accurately estimate the error at all FIG. 27. As in Fig. 26, but now for ψ 4 . The large spike at u ¼ 0 (not fully shown, and not convergent) is a signal of the irregularity of our initial conditions. All ψ l with l > 2 have this problem. FIG. 28. As in Fig. 26, but now for C.
higher resolutions by scaling these numbers by factors ðΔx=0.05Þ 2 and ðl max =16Þ -2 , respectively. The computation of the diagnostic C in the collapsing solution seems to be particularly challenging. Here we see clear second-order convergence, even at early times, only for Ny ≥ 25, see the lower plot in Fig. 28.
To check the effect of our nonanalytic original choice of initial data (where ψ is not single-valued at the origin) on convergence, we have created a second family of initial data that are analytic.
We computed the exact solution ϕ ð2kÞ ðt; r; yÞ of the flat spacetime scalar wave equation with Cauchy data ϕ ð2kÞ ð0; r; yÞ ¼ e -r 2 ð-ryÞ 2k ; ðC1Þ ϕ ð2kÞ;t ð0; r; yÞ ¼ 0; ðC2Þ for k ¼ 0; 1; …5. (Recall that y ≔cos θ. We have set the width of the Gaussian to d ¼ 1 to fix an overall scale). These can be found as linear combinations of the generalized d'Alembert solutions ϕ l ðt; rÞP l ðyÞ constructed in Paper I, for even l up to 2k, each with an ansatz for the free function χ l ðrÞ of that is expð-r 2 Þ times an odd polynomial in r of order l þ 1.
We then use the ϕ ð2kÞ ðt; r; yÞ as building blocks to construct the solution with initial data φð0; r; yÞ ¼ pe -r 2 X 5
¼ pe -r 2 ð1-ϵ 2 y 2 Þ þ Oðϵ 6 r 12 y 12 Þ; ðC4Þ ϕ ;t ð0; r; yÞ ¼ 0:
We then read off the desired null data at u ¼ 0 from the full solution as ψð0; r; yÞ ¼ φðr; r; yÞ:
We were able to compute up to k ¼ 6 in Mathematica, and the error in the Taylor series is then below machine precision. However, in contrast to the k ≤ 5 terms, the k ¼ 6 term would need to be approximated near r ¼ 0 to avoid large roundoff error when evaluating the exact TABLE III. Successful initial bracketings of the collapse threshold with the analytic initial data (C6). In each case, p ¼ 0.2 disperses and p ¼ 0.3 collapses. x max needs to be this large only for the collapsing solutions. x max is twice as large as for the nonanalytic data (10) for ϵ 2 ¼ 0.5, and three times as large for ϵ 2 ¼ 0.75. x 0 has not yet been fine-tuned for bisection to the black hole threshold. We have used C ≥ 0.8 as the collapse threshold and C ≤ 0.05 as the dispersion threshold. 2 Δx Ny l max x 0 x max 0.5 0.025 9=17 8=16 29 30 0.6 0.025 9=17 8=16 39 40 0.7 0.05 17 16 99 100 0.75 0.05 17 16 149 150 TABLE IV. Table of the maximum in u of the L 2 norm in x and y of the discretization errors in x and y, in the three variables ψ 0 , ψ 0 and C, in the p ¼ 0.2 (dispersing) solution, at the baseline resolution Δx ¼ 0.0.5, Ny ¼ 17, l max ¼ 16. Top: nonanalytic data. Bottom: analytic data. Nonanalytic ψ 0 ψ 4 C max u jjE x jj L 2 1.3 × 10 -3 8 × 10 -5 3.5 × 10 -3 max u jjE y jj L 2 1.7 × 10 -4 3 × 10 -5 3 × 10 -4 Analytic ψ 0 ψ 4 C max u jjE x jj L 2 2 × 10 -3 5 × 10 -4 7 × 10 -3 max u jjE y jj L 2 1.2 × 10 -2 4 × 10 -3 2 × 10 -2 FIG. 29. Error in ψ 0 as in Fig. 26, but now for the analytic initial data (C6) with ϵ 2 ¼ 0.75, and with errors in y now scaled by exp½-κðl max -16Þ with κ ≃ 0.280 for l max ¼ 8, 12, 16 and κ ≃ 0.085 for l max ¼ 24, 32, 48.
expression numerically in the code, so we for the above truncation at k ¼ 5.
Table III shows successful initial brackets of the black hole threshold with the analytic data. We have not bisected the analytic data to the black hole threshold.
For convergence testing, we have evolved ϵ 2 ¼ 0.75, p ¼ 0.2 and p ¼ 0.3, at the same two intersecting families of resolution as for our nonanalytic data. We again see second-order convergence in Δx. Power-law convergence in l max is no longer a good fit. A better fit is exponential convergence in l max with a break in the exponent, namely exp½-κðl max -16Þ with κ ≃ 0.280 for l max ¼ 8, 12, 16 and κ ≃ 0.085 for l max ¼ 24, 32, 48.
Numerical values for the errors at the baseline resolution are shown in Table IV, and L 2 norms over x of the in ψ 0 , ψ 4 and C are shown in Figs. 29, 30, and 31. For completeness, we show the unscaled errors in y, for the variable C, and for both the analytic and nonanalytic ϵ 2 ¼ 0.75, p ¼ 0.2 data in Fig. 32. FIG. 31. As in Fig. 29, but now for C. FIG. 32. Direct comparison of the unscaled discretization errors in y for the nonanalytic (dashed lines) and analytic (solid lines) ϵ 2 ¼ 0.75, p ¼ 0.2 data. As before Ny ¼ f9; 13; 17; 25; 33; 49g are shown in purple, green, blue, orange, yellow, red, and we estimate the error by subtracting the Ny ¼ 65 evolutions.