Simulations of gravitational collapse in null coordinates: II. Critical collapse of an axisymmetric scalar field
Carsten Gundlach Thomas W. Baumgarte David Hilditch









Mathematical Sciences, University of Southampton,
  Southampton SO17 1BJ, United Kingdom 

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

CENTRA, Departamento de Física, Instituto Superior
  Técnico IST, Universidade de Lisboa UL, Avenida Rovisco Pais 1,
  1049 Lisboa, Portugal 





  
We present the first numerical simulations in null coordinates of  the collapse of nonspherical regular initial data to a black  hole. We restrict to twist-free axisymmetry, and re-investigate the  critical collapse of a non-spherical massless scalar field. We find  that the Choptuik solution governing scalar field critical collapse  in spherical symmetry persists when fine-tuning moderately  non-spherical initial data to the threshold of black hole  formation. The non-sphericity evolves as an almost-linear  perturbation until the end of the self-similar phase, and becomes  dominant only in the final collapse to a black hole. We compare with  numerical results of Choptuik et al, Baumgarte, and Marouda et al,  and conclude that they have been able to evolve somewhat more  non-spherical solutions. Future work with larger deviations from  spherical symmetry, and in particular vacuum collapse, will require  a different choice of radial coordinate that allows the null  generators to reconverge locally.


Simulations of gravitational collapse in null coordinates:
  II. Critical collapse of an axisymmetric scalar field
    David Hilditch
    24 April 2024
==================================================================================================================







§ INTRODUCTION





In a companion paper <cit.> (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
centre, intending to simulate gravitational collapse with them. 

As a first physics application, here we re-investigate 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
<cit.>.

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 <cit.> it was claimed, based on the numerical
solution of a mode ansatz, that all non-spherical 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 <cit.> appeared to indicate a
slowly growing l=2 perturbation of slightly non-spherical critical
collapse that eventually leads to the formation of two centres of
collapse.

This tension between the results of <cit.> and
<cit.> was resolved by fresh simulations, in
axisymmetry in spherical coordinates by Baumgarte
<cit.>, which indicated that sufficiently small
non-spherical perturbations decay (in agreement with the mode
stability claimed in <cit.>) but that larger
perturbations do grow (in agreement with <cit.>).

Interestingly, both <cit.> and <cit.>
observe that in a regime of small finite deviations from spherical
symmetry, the spherical Choptuik solution
<cit.> is still observed as an
approximate critical solution but with a period that decreases from
Δ≃ 3.44 with increasing non-sphericity, and that the
critical exponent also decreases from γ≃ 0.374.

The authors of <cit.> and <cit.>, in
axisymmetry, achieved fine-tuning to |p-p_*|∼ 10^-15 and
|p-p_*|∼ 10^-13, respectively (except that <cit.>
achieved only 10^-8 at the largest deviation from spherical
symmetry). Other investigations of non-spherical scalar field collapse
had achieved much less fine-tuning: <cit.> achieved
∼ 10^-3.5, and <cit.> ∼ 10^-5.5, both for
non-axisymmetric initial datap, and <cit.> achieved
∼ 10^-1 in axisymmetry. Incidentally,
<cit.> also evolve initial
data that are non-smooth at R=0.

Reid and Choptuik <cit.> have fine-tuned axisymmetric,
time-symmetric, antisymmetric in z initial data to ∼ 10^-14
and found that two (mirror-symmetric) centres of collapse form (as one
would expect). These start out in a highly non-spherical 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 non-spherical
mode (in the evolution of each separate centre of collapse).

Marouda et al <cit.> have examined twist-free
axisymmetric collapse of a complex massless scalar field, fine-tuning
to ∼ 10^-10. As they treat a slightly different matter model,
none of their families of initial data coincides exactly with those of
<cit.> and <cit.>, but their results
are consistent with those of <cit.> and
<cit.> in the sense that δΔ and
δγ observed in their two non-spherical families fit into a
monotonic ordering of families by (inferred) non-sphericity. See
already Table <ref> below.

The achievable level of fine-tuning matters, as the l=2 unstable
mode claimed in <cit.> would be very slowly-growing:
its authors see no bifurcation for non-sphericity parameter
ϵ^2≤ 2/3 [see Eq. (<ref>) below for a definition],
but they see one after 2 and 1.5 echos, respectively, for
ϵ^2=3/4 and 5/6, where 3 echos is the best that can be
achieved in double-precision arithmetic, with |p-p_*|∼ 10^-15.
<cit.> see about 4 echos 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 <cit.>, whereas our choice
of null gauge should be able to resolve the critical solution well,
even without the adaptive mesh refinement of
<cit.>. It turns out that this is true.

The plan of the paper is as follows: In Sec. <ref> we
briefly summarise our metric ansatz, gauge choice and diagnostics,
give our two-parameter family (strength and non-sphericity) of initial
data, and describe the bisection to the threshold of collapse. In
Sec. <ref> 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. 

Sec. <ref> gives our results for non-spherical
(twist-free axisymmetric) critical collapse. Our initial data, being
set on an outgoing null cones, cannot be compared directly with those
of <cit.>. 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 non-spherical family, ϵ_2=0.75, compares to
ϵ^2=0.5 of <cit.>. With the
gauge choice presented here we cannot fine-tune families of initial
data that are more non-spherical.

We conclude in Sec. <ref> with a discussion of the
physical results and an outlook for our code. 





§ SETUP







 §.§ Metric and field equations




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

    R_ab=8π∇_aψ∇_bψ,

and the massless, minimally coupled wave equation

    ∇^a∇_aψ=0.

We use units where c=G=1.

We write the general twist-free axisymmetric metric in the form 

    ds^2   =   -2G du dx-H du^2 
          +R^2[e^2S fS^-1(dy+S b du)^2
    +e^-2S fS dφ^2].

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
(<ref>) we have used the shorthand

    S:=1-y^2=sin^2θ.


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 parameterised generators of N^+_u.  The ingoing
future-directed null vector
normal to S_u,x is

    Ξ:=∂_u -B∂_x-Sb ∂_y,

where we have defined the shorthand

    B:=H 2G.

As we see from (<ref>) 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
normalised relative to each other as Ξ^aU_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).






 §.§ Gauge choice




Our numerical domain is 0≤ x≤ x_max, 0≤ u≤
u_max. 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 x=x_max. We set
H=G=1 at x=R=0, 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 lightcone of the accumulation point of scale echos. 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
<cit.>, using regridding, and in
<cit.> 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” flavour, defined by the shift B defined in (<ref>)
taking the value

    B_lsB4=-Ξ R-min_y(Ξ R) R_,x
    +(1-x x_0)(-Ξ R R_,x)_x=0.

(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 continous but has discontinuous
transverse derivative across the hypersurfaces x=x_*(u) where the
y-location of min_y(Ξ R) jumps.





 §.§ Diagnostics




For any closed spacelike 2-surface S, we define its “Hawking
compactness”

    C( S):=1+1 2π∫_ Sρ_+ρ_- dS,

where ρ_+ and ρ_- are the outgoing and ingoing null
divergence. It is related to the well-known Hawking mass by

    M( S):=1 2√(A( S) 4π) C( S),

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_xC(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_xC(u,x)≤ 0.01 as an indicator of
dispersion. In spherical symmetry, 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
discretised in different ways, see Paper I for details.






 §.§ Initial data




The authors of <cit.> investigate non-spherical
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 <cit.>, ours is optimised for a single
centre 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 <cit.>.

The time-symmetric initial data of <cit.> and
<cit.> at t=0, written in our notation, are

    ϕ(0,r,y)=p e^-r^2(1-ϵ_2y^2) d^2,     ϕ_,t(0,r,y)=0.

With ϵ_2>0, the initial data are elongated along the symmetry
axis y=±1, and with ϵ_2<0 they are squashed. In
<cit.>, ϵ_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,
<cit.> 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. As an approximate equivalent, we consider the function

    ϕ(t,r,y)   :=   1 2r[
    (t+r)e^-(t+r)^2 d^2(1-ϵ_2y^2)
          -(t-r)e^-(t-r)^2 d^2(1-ϵ_2y^2)].

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 (<ref>) at t=0. In our null
code, we now set the free initial data at u=0 to

    ψ(0,R,y)=p ϕ(t_0+R,R,y)

for some constant t_0, so that u=0 corresponds to
t=t_0+R. We also set

    f(0,x,y) = 0.

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 2,

which is essentially a gauge choice.

Formally, as p→ 0 and ϵ_2→ 0, the solution
from these data corresponds to that arising from (<ref>), with
the identification u=t-R. We also initialise 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, (<ref>)
is no longer an exact solution of the curved-spacetime wave equation,
and the solution evolved from the null initial data
(<ref>) no longer has a moment of exact time symmetry, but it
should approximate the solutions arising from (<ref>), 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
(<ref>) represents a Gaussian that is well separated from the
centre and still essentially ingoing. In the flat-space, spherical
limit ϵ_2=0, p→ 0 the wave reaches the centre 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 realised only near the completion of this paper, from the
convergence tests in Appendix <ref>, that
for ϵ_2 0 and t 0 (<ref>) is not actually
single-valued at the origin r=0, but depends on y. Our null
initial data (<ref>) are therefore also not single-valued at the
origin. In Appendix <ref> we present an
analytic solution of the flat space wave equation valid for all values
of ϵ_2 that also reduces to (<ref>) for all values of
ϵ_2, and null initial data (<ref>) derived
from this solution. These are the initial data we intended to
create. 

However, convergence tests in the strong-gravity regime of the
non-analytic null data (<ref>) with (<ref>) 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 corrrected null data
(<ref>) also converge, but preliminary tests show
that at ϵ_2=0.75 they require very large values of
x_max in near-critical evolutions. For both these reasons —
the error in the non-analytic data does not seem to matter, and the
analytic data are much harder to evolve — we present here only
evolutions with the non-analytic initial data, and we believe that the
following results are still reliable.





 §.§ Bisection to the threshold of black hole formation




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(p-p_*)≃-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 echos 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|p-p_*| 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. We take this combination as an indicator of collapse, even if
C=0.99 has not yet been reached. At the largest non-sphericity 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 mass-like 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|p-p_*|, with 30 values per power of 10, to
give us better-resolved plots of the mass and curvature scaling laws.





§ RESULTS IN SPHERICAL SYMMETRY





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 lightcone of the accumulation point of scaling echos 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 N̅_y=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 non-spherical 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 the Ricci
scaling (see Fig. <ref>) and at about super12 for
the mass scaling (see Fig. <ref>).

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.





 §.§ Fine structure of the mass and curvature scaling laws




From Appendix <ref>, which clarifies the
  derivation given in <cit.>, we expect the scaling laws

    -12ln T   ≃   A+γln(p_*-p)
           +f_T(ln(p_*-p)+B), 
    ln M   ≃   A+γln(p-p_*)
          
    +f_M(ln(p-p_*)+B),

for

    T:=max_whole spacetime|∇_aψ∇^aψ|,

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 a given DSS-compatible time slicing (such as our null
cone slicing), but does depend on the slicing, as well as our collapse
criterion.

Note that the Ricci scalar is |Ric|=8π T, so
ln|Ric| and ln T differ only by a constant.

To look for these scaling laws, we plot ln M-(A+γln|p-p_*|)
against lnln|p-p_*|, 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|p_*-p|). 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 (<ref>) for our family
of spherically symmetric initial data.









Our numerical measurement of f_T with these fitting parameters is
shown in Fig. <ref>. 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. <ref>. 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 centre, 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 <cit.> 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
<cit.> 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 <cit.> we estimate the
amplitude in ln M as ≃ 0.4, similar again to
<cit.>, 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 discontinous, 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 <cit.> plot the
fine-structure in lnmax|Ric|, and from this plot we
estimate the amplitude of f_T as 0.3 in lnmax|Ric|, or
0.15 in (1/2)lnmax|T|. 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
<cit.> has a fine-structure with amplitude of about 0.35 in
lnmaxρ, where ρ=n^an^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 <cit.> for
lnmax|Ric|.





 §.§ Self-similarity of the critical solution




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|∇_aψ∇^aψ|

are periodic in the DSS-adapted coordinates

    ξ   :=   R  u_*-u, 
    τ   :=   -ln(u_*-u).

The null slicing is already DSS-compatible, in the sense of
<cit.>, but we call this specific coordinate system
DSS-adapted 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 ψ_*(ξ,τ+Δ/2)=-ψ_*(ξ,τ)
and R_ab=8πψ_,aψ_,b is invariant under ψ→-ψ.
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 echos. 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
echos, 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. <ref>,
<ref> and <ref> for our closest
subcritical evolution.

















 §.§ Comparison with gsB gauge





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 

    (B)=ṡ(u)x-Ξ R s(u),

and we have settled on the particular flavor given by

    ṡ(u)=1 x_0min_yΞ R(u,x_0,y).

 
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.

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 non-monotonic 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.





§ RESULTS FOR NONSPHERICAL INITIAL DATA





As already noted, the strong-field convergence tests with
ϵ_2=0.75 in Appendix <ref>
showed an initial burst of non-convergent 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 non-analytic 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 parameterised by (p,ϵ_2) is equivalent to the
family of Cauchy data parameterised by (p,ϵ^2) of
<cit.> 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 <ref> 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 N̅_y=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 (<ref>), we expect that for sufficiently small
ϵ_2, the non-sphericity evolves as a linear perturbation of
spherical symmetry, dominated by l=2, and that this linear
perturbation decays. For somewhat larger ϵ^2,
<cit.> also found changes to Δ
and γ. For even larger ϵ^2, in particular
ϵ^2=0.75, Baumgarte <cit.> also found evidence
of the (nonlinear) instability of the Choptuik solution reported by
Choptuik et al <cit.> to lead to two centres 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 non-sphericity is more similar to the
evolutions of Baumgarte <cit.> with his
ϵ^2=0.5.










 §.§ The final collapse phase





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 non-spherical, 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 Ξ 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

    R̅_,u=Ξ R+BR_,x

tries to counteract the non-spherical 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
non-spherical 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 non-spherical 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 non-decreasing by construction. Of course
either M̃ or M should be considered reliable only to the
extent they agree with each other.





 §.§ Scaling and echoing




We look for three types of evidence that the non-spherical
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 <cit.>. We examine evidence for
all this in the next two subsections, starting here with the scaling
laws.

Fig. <ref> and
Fig. <ref> 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 non-spherical
initial data, Figs. <ref> and
Fig. <ref> 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. <ref> and
Fig. <ref> 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. <ref> is real but, we do not
have enough numerical resolution for quantitatively reliable results.

Our plots of f_T in Fig. <ref> at different
values of ϵ_2 can be compared directly with Fig. 3 of
<cit.> and Fig. 7 of <cit.>, with the
differences that we plot -1/2ln T+A rather than ln T on the
vertical axis, and ln(p_*-p)+B rather ln(p_*-p) on the
horizontal axis. In other words, <cit.> and
<cit.> 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
<cit.>. In our most
non-spherical evolutions, ϵ_2=0.75, we only see a change of
δγ≃ -0.0016. By comparison, <cit.> 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 <ref>).
 
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. <ref> and
<ref>, 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. <ref>
for a fit of f_T at ϵ_2=0 to a sine wave of period 3.44, and
Fig. <ref> 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, <cit.> find δΔ=-0.05...-0.12 at their
ϵ^2=0.5 (see also Table <ref>).

Unfortunately, <cit.> and <cit.> do not
present plots of f_M. Our own results in
Fig. <ref> 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 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.






























 §.§ Self-similarity and evolution of the nonsphericity




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. 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 linearised equations by
angular dependence l, and for each l and dimensionless field
ϕ make a mode ansatz

    δϕ_l(ξ,τ)=e^λτ δϕ̃_lλ(ξ,τ),

with δϕ̃_lλ defined to be periodic in τ
with period Δ, and regular at the centre and at the past
lightcone.  Here ϕ stands for any scale-invariant quantity, such
as ψ or T. From this ansatz, the complex mode function
δϕ̃_lλ and corresponding complex Lyapunov
exponent

    λ:=κ+iω

are then determined as the eigenfunctions and eigenvalues of a
(singular) linear boundary value problem <cit.>. 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λ(ξ,τ)] 
       =   Ae^κτ[cos(ωτ+α)  Re δϕ̃_lλ(ξ,τ) 
          -sin(ωτ+α)  Im δϕ̃_lλ(ξ,τ)],

for arbitrary positive real amplitude A and phase 0≤α<2π.  (Note
this is slightly incorrect in Eq. (25) of Ref. <cit.>).
Unless ω is 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 <cit.>, using different
similarity coordinates from the ones defined here. We can therefore
compare the λ with <cit.>, but not directly
the δϕ̃_lλ(τ,x). As expected for a critical
solution, <cit.> finds a single growing l=0 mode,
with λ and the corresponding mode function real. All other
spherical modes, and all non-spherical modes, are complex and decay,
that is, κ<0. The least damped (most slowly decaying)
non-spherical 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 single-mode formula
(<ref>) 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 non-spherical 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 ϵ_2e^κτ, all
l=4 components to be proportional to ϵ_2^2e^2κτ,
and so on for higher l.

To demonstrate the expected scaling with both τ and and
ϵ_2, in Fig. <ref> we plot the maximum and
minimum of ϵ_2^-1e^-κτψ_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. <ref> for f_2. We find that
ϵ_2^-1e^-κτψ_2 is essentially the same at
ϵ_2=10^-4, 10^-2 and 0.1, and so is
ϵ_2^-1e^-κτ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 <cit.> 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 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. <ref> and <ref>.

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^-2e^-2κτ. This is demonstrated in
Figs. <ref>-<ref>.


































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
(<ref>) is

    ω_N:=ω+N2πΔ,     N∈ Z,

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
(<ref>), but for a random α,
<cit.> found that the peak of the spectrum was at
N=5. The corresponding period in τ is

    P_N:=2π|ω_N|=Δ|N+ωΔ 2π|≃3.44 |N+0.3|.

If we count the peaks in Fig. <ref>, 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. <ref>, 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.)





 §.§ Comparison with the evolutions of Choptuik et al,
  Baumgarte, and Marouda et al




In order to compare the non-sphericity of our evolutions with those of
<cit.> (and by implication also of
<cit.> for the same initial data), we note that
Baumgarte has measured the difference in ψ on outgoing null
geodesics at y=± 1 and y=0 (poles and equator). Fig. 11 of
<cit.> 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, normalised to λ=0 and dλ/dR=1 at the
origin, and u_* is fitted as described above. This measure is
gauge-invariant.

We note that for infinitesimal deviations from spherical symmmetry,
and assuming the deviation only has an l=2 component,

    δψ(u,x)   :=   ψ(u,x,±
    1)-ψ(u,x,0) 
       =   [P_2(±1)-P_2(0)]ψ_2(u,x) 
       =   3 2ψ_2(u,x),

where ψ_2 denotes the l=2 component of ψ. λ along
outgoing null geodesics is given by

    λ(u,x,y)=∫_0^x G(u,x',y) dx'


The equivalent of Fig. 11 of <cit.> created from our
data for ϵ_2=10^-2 is
Fig. <ref>. 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
<cit.>, but recall that we set initial data on different
slices.

Fig. <ref> shows δψ, approximated as
(3/2)ψ_2 for our ϵ_2=0.75. For comparison,
Figs. <ref> and
<ref> show plots of δψ
created from the simulations of <cit.> 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 non-spherical as ϵ^2=0.5 of <cit.>. In
both families δψ does not grow noticeably. By contrast,
Fig. <ref> shows clear growth of
δψ with τ for ϵ^2=0.75 and an end to the
perturbative regime around τ=6.

Table <ref> compares the black hole threshold in
selected families taken from <cit.>,
<cit.>, <cit.> 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|δψ| from the null code
and that of <cit.> are read off from plots shown in the
present paper. The values for ϵ^2=0 and ϵ^2=2/3 of
<cit.> are read off from Figs. 6 and 7 of that paper,
respectively. In the spherically symmetric case we include
max|δψ| 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
<cit.> are read off from plots communicated by the
authors.

A range given for max|δψ| 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 <cit.>, 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|δψ| (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 change the values of Δ and γ, as
well as the decay rate κ of the l=2 mode, such that for
non-sphericities of max|δψ|≳ 0.1 it actually
grows. (Compare this with max|ψ_*|≃ 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 non-sphericity. For lack 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.  



























 §.§ Numerical problems at larger non-sphericity




We have already mentioned that to complete the bisection in p for
ϵ_2=0.75 at our highest resolution N̅_y=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 N̅_y=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.





§ CONCLUSIONS





This work was born out of an attempt to generalise the method of
Garfinkle <cit.> 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) <cit.> seemed to have reconciled the linear
perturbation results of Martín-García and Gundlach
<cit.> with the collapse simulations of Choptuik et al <cit.> (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 centres of collapse. This conclusion is also
supported by the more recent evolutions of <cit.> and
(less quantitatively) of <cit.>.

Specifically, <cit.> and <cit.> find
that non-sphericities decay for values of their non-sphericity
parameter ϵ^2≤ 0.5, and <cit.> still find
this at ϵ^2=2/3. Both <cit.> and
<cit.> find a bifurcation for ϵ^2=0.75, and
<cit.> also at ϵ^2=5/6. (We note that
<cit.> also observe bifurcation in their most
non-spherical data for the complex scalar field).
<cit.>, <cit.> and
<cit.> also found that Δ and γ (measured
before any bifurcation occurs) decrease with increasing ϵ_2, see
Table 1 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 non-sphericity. In particular, the echoing period and critical
exponent remain at Δ≃ 3.44 and γ≃ 0.374, their
values in spherical symmetry, and the l=2 non-spherical components
of all fields decay at roughly the rate predicted for linear
perturbations in <cit.>. 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
non-sphericity 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
<cit.> predicts l=4 to decay much more rapidly than
l=2.

However, we cannot evolve the same family of initial data as
<cit.> and <cit.>, because we set data
on an outgoing null cone. We have compared a gauge-invariant measure
of the non-sphericity of the scalar field, in the phase where a
near-critical evolution is approximately self-similar, with the same
measure for the evolutions of <cit.>, and find that this
comparable between our ϵ_2=0.5 and ϵ^2=0.75 of
<cit.> (and by implication of
<cit.>). If we compare our ϵ_2=0.5 with
ϵ^2=0.75 of <cit.>, 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 compactness 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
<cit.> 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 non-sphericity 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 non-sphericity to fully duplicate the results of
<cit.>: at large
non-sphericity, 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 lightcones is trying to become negative
at some points on our last null cone, but the generalised Bondi
coordinates we use break down if this happens. Therefore, the next
step will be to change to a generalised affine radial coordinate. This
will then allow us to investigate if our null cones themselves remain
regular in critical collapse, or form caustics. 





We would like to thank Bernd Brügmann, Tomáš Ledvinka, Anton
Khirnov, Daniela Cors and Krinio Marouda for helpful discussions, and
the Mathematical Research Institute Oberwolfach for supporting this
work through its “Oberwolfach Research Fellows” scheme.  DH was
supported in part by FCT (Portugal) Project No. UIDB/00099/2020. TWB
was supported in part by National Science Foundation (NSF) grants
PHY-2010394 and PHY-2341984 to Bowdoin College.













§ FAMILY-DEPENDENT PARAMETERS IN THE FINE STRUCTURE OF THE SCALING LAWS





In this Appendix, we show that the parameters A and B in
Eqs. (<ref>) and (<ref>) are family-dependent, and
independent of each other.

Focus first on the late-time evolution of the exactly-critical member
of a given 1-parameter family of initial data. If we fix a
(family-independent) small length scale L, then we can, for example,
measure the phase of the scalar field ϕ at the centre when the
Ricci scalar at the centre takes exactly the value L^-2. This is
equivalent to our A.

On the other hand, to derive the scaling laws we use the fact that in
near-critical evolutions the amplitude of the growing perturbation is,
to leading order, proportional to (p-p_*)e^λ_0τ. The constant of
proportionality must again be family-dependent, already for the
trivial reason that p can have any dimension. Its logarithm is our B.

Note that for perfectly-critical initial data the growing perturbation
of the critical solution is by construction absent, while the scaling
laws rely on the amplitude of the growing mode. Hence the
family-dependent constants A and B are independent of each
other. 





§ CONVERGENCE TESTS IN THE STRONG SUBCRITICAL REGIME





For convergence testing we choose our ϵ_2=0.75 initial data
with the two amplitudes that we have also used as an initial bracket
of the black hole threshold, namely p=0.2 (subcritical) and p=0.3
(supercritical). We use the grid parameters x_max=50 and
x_0=33 throughout, which we also used in our near-critical evolutions.

For p=0.2 gravity is strong but not close to collapse, with a
maximum Hawking compactness of C̅≃ 0.08 in the initial
data, and reaching C̅≃ 0.2 during the evolution (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 has 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 timesteps. 

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
|p_*-p|≃ 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 N̅_y 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={0.1, 0.05,0.025,0.0125} with fixed angular
resolution N̅_y=17, l_max=16; and at
N̅_y={9,13, 17,25,33,49,65}, l_max=N̅_y-1,
with fixed radial resolution Δ x=0.05. Their intersection is
our baseline solution, chosen because the discretisation 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 non-converging error at small u and x, in all
ψ_l with l>2. We believe this is due to the fact that our
initial data (<ref>-<ref>) is not actually single-valued
at the origin. The non-convergent error disappears quickly, and we
believe this is because the boundary conditions at the centre 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
discretisation error in y decreases with resolution only as
l_max^-2.  Second-order convergence in the L^2 norm with
Δ x and 1/l_max is demonstrated for ψ_0 in
Fig. <ref>, for ψ_4 in
Fig. <ref>, and for C̅ in
Fig. <ref>. The convergence is
actually pointwise. The global maximal errors max_u ||
  E_x||_L^2 and max_u|| E_y||_L^2 for the dispersing
are shown in Table <ref>. Given second-order convergence,
we can accurately estimate the error at all 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 N̅_y≥ 25, see the
lower plot in Fig. <ref>.











































§ CONVERGENCE TESTS WITH INITIAL DATA THAT ARE REGULAR AT THE ORIGIN





To check the effect of our non-analytic 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, 
    ϕ_(2k),t(0,r,y)   =   0,

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 generalised 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∑_k=0^5 (ϵ_2r^2y^2)^k k!
       =    
    p e^-r^2(1-ϵ_2y^2) +O(ϵ^6 r^12 y^12), 
    ϕ_,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 expression numerically in the code, so we settle for the
above truncation at k=5.

Table <ref> 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 non-analytic 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 exponentional 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 <ref>, and L^2 norms over x of the errors in
ψ_0, ψ_4 and C̅ are shown in
Figs. <ref>,
<ref> and
<ref>.  For completeness,
we show the unscaled errors in y, for the variable C̅, and for
both the analytic and non-analytic ϵ_2=0.75, p=0.2 data in
Fig. <ref>.




99




paper1 C. Gundlach, D. Hilditch and T. W. Baumgarte,
  Simulations of gravitational collapse in null coordinates I:
  Formulation and weak-field tests in generalised Bondi gauges,
  https://arxiv.org/abs/2404.15105arXiv:2404.15105.

GundlachLRR C. Gundlach and J. M. Martín-García,
  Critical phenomena in gravitational collapse,
https://link.springer.com/article/10.12942/lrr-2007-5Living
  Rev. Relativ. 7, 05 (2007).

JMMGundlach1999 J. M. Martín-García and C. Gundlach,
  All nonspherical perturbations of the Choptuik spacetime decay, 
  https://journals.aps.org/prd/abstract/10.1103/PhysRevD.59.064031Phys. Rev. D 59, 064031
  (1999). 

Choptuiketal2003 M. W. Choptuik, E. W. Hirschmann,
  S. L. Liebling and F. W. Pretorius, Critical collapse of the
  massless scalar field in axisymmetry, https://journals.aps.org/prd/abstract/10.1103/PhysRevD.68.044007Phys. Rev. D 68, 044077
  (2003). 

Baumgarte2018 T. W. Baumgarte, Aspherical deformations of
  the Choptuik spacetime, https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.084012Phys. Rev. D 98, 084012 (2018). 

Choptuik1993 M. W. Choptuik, Universality and scaling in
  gravitational collapse of a massless scalar field,
  https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.70.9Phys. Rev. Lett. 10, 1103 (1993).

MartinGundlach2003 J. M. Martín-García and
  C. Gundlach, Global structure of Choptuik’s critical solution in
  scalar field collapse,
  https://journals.aps.org/prd/abstract/10.1103/PhysRevD.68.024011Phys. Rev. D
    68, 024011 (2003).

Healeyetal2014 J. Healy and P. Laguna, Critical collapse of
  scalar fields beyond axisymmetry, https://link.springer.com/article/10.1007/s10714-014-1722-2Gen. Relativ. Gravit. 46, 1722
    (2014).

Deppetelal2019 N. Deppe, L. E. Kidder, M.
  A. Scheel and S. A. Teukolsky, Critical behavior in 3D
  gravitational collapse of massless scalar fields, https://journals.aps.org/prd/abstract/10.1103/PhysRevD.99.024018Phys. Rev. D
    99, 024018 (2019).

CloughLim2016 K. Clough and E. A. Lim, Critical Phenomena in
  Non-spherically Symmetric Scalar Bubble Collapse, unpublished, https://arxiv.org/abs/1602.02568https://arxiv.org/abs/1602.02568.

Reid2023 G. D. Reid and M. W. Choptuik, Universality in the
  critical collapse of the Einstein-Maxwell system,
  https://journals.aps.org/prd/abstract/10.1103/PhysRevD.108.104021Phys. Rev. D
    108, 104021, (2023).

Maroudaetal2024 K. Marouda, D. Cors, H. R. Rüter, F. Atteneder and D. Hilditch, Twist-free axisymmetric critical collapse of a complex scalar field,
https://arxiv.org/abs/2402.06724arXiv:2402.06724 

Garfinkle1995 D. Garfinkle, Choptuik scaling in null
  coordinates, https://journals.aps.org/prd/abstract/10.1103/PhysRevD.51.5558Phys. Rev. D 51, 5558
  (1995).

ymscalar C. Gundlach, T. W. Baumgarte and D. Hilditch, Critical
  phenomena in gravitational collapse with two competing massless
  matter fields, https://journals.aps.org/prd/abstract/10.1103/PhysRevD.100.104010Phys. Rev. D 100, 104010 (2019).

Rinne2020 O. Rinne, Type II critical collapse on a single
  fixed grid: a gauge-driven ingoing boundary method, https://link.springer.com/article/10.1007/s10714-020-02768-xGen. Relativ. Gravit. 52, 117
  (2020). 

PortoGundlach2022 B. Porto-Veronese and C. Gundlach,
  Critical phenomena in gravitational collapse with competing scalar
  field and gravitational waves in 4+1 dimensions,
  https://journals.aps.org/prd/abstract/10.1103/PhysRevD.106.104044Phys. Rev. D
    106, 104044 (2022).

Gundlach1997 C. Gundlach, Understanding critical collapse of
  a scalar field, https://journals.aps.org/prd/abstract/10.1103/PhysRevD.55.695Phys. Rev. D 55, 695
  (1996). 

PuerrerHusaAichelburg M. Pürrer, S. Husa and
  P. C. Aichelburg, News from critical collapse: Bondi mass, tails,
  and quasinormal modes, https://journals.aps.org/prd/abstract/10.1103/PhysRevD.71.104005Phys. Rev. D 71, 104005
  (2005).
 
CrespoOliveiraWinicour2019 J. A. Crespo, H. P. de Oliveira and
  J. Winicour, The affine-null formulation of the gravitational
  equations: spherical case, https://journals.aps.org/prd/abstract/10.1103/PhysRevD.100.104017Phys. Rev. D 100, 104017 (2019). 

GarfinkleDuncan1998 D. Garfinkle and G. C. Duncan, Scaling
  of curvature in subcritical gravitational collapse, https://journals.aps.org/prd/abstract/10.1103/PhysRevD.58.064024Phys. Rev. D
  58, 064024 (1998).

Corsetal2023 D. Cors, S. Renkhoff, H. R. Rüter,
  D. Hilditch and B. Brügmann, Formulation improvements for critical
  collapse simulations,
  https://journals.aps.org/prd/abstract/10.1103/PhysRevD.108.124021
Phys. Rev. D 108, 124021, (2023).

Siebeletal2003 F. Siebel, J. A. Font, E. Müller and
  P. Papadopoulos, Axisymmetric core collapse simulations using
  characteristic numerical relativity,
  https://journals.aps.org/prd/abstract/10.1103/PhysRevD.67.124018Phys. 
Rev. D 67, 124018 (2003).







