Formulation Improvements for Critical Collapse Simulations

The precise tuning required to observe critical phenomena in gravitational collapse poses a challenge for most numerical codes. First, threshold estimation searches may be obstructed by the appearance of coordinate singularities, indicating the need for a better gauge choice. Second, the constraint violations to which simulations are susceptible may be too large and force searches to terminate prematurely. This is a particularly serious issue for first order formulations. We want our adaptive pseudospectral code bamps to be a robust tool for the study of critical phenomena so, having encountered both of these difficulties in work on the vacuum setting, we turn here to investigate these issues in the classic context of a spherically symmetric massless scalar field. We suggest two general improvements. We propose a necessary condition for a gauge choice to respect discrete self-similarity (DSS). The condition is not restricted to spherical symmetry and could be verified with any 3+1 formulation. After evaluating common gauge choices against this condition, we suggest a DSS-compatible gauge source function in generalized harmonic gauge (GHG). To control constraint violations, we modify the constraint damping parameters of GHG, adapting them to collapse spacetimes. This allows us to improve our tuning of the critical amplitude for several families of initial data, even going from 6 up to 11 digits. This is the most precise tuning achieved with the first order GHG formulation to date. Consequently, we are able to reproduce the well known critical phenomena as well as competing formulations and methods, clearly observing up to 3 echoes.


I. INTRODUCTION
The threshold of gravitational collapse separates spacetimes at the verge of black hole formation from those in which a black hole forms.Spacetimes near this region in solution space have large, dynamical curvature and, depending on the specific model under consideration, infinitesimal black holes.This makes them extremely interesting, albeit hard to treat.In 1993, Choptuik [1] tackled the problem by evolving massless scalar fields minimally coupled to the Einstein field equations in spherical symmetry.He observed three features bearing a strong resemblance to those observed near critical points in other fields of physics.He noted that, near the threshold, the spacetime becomes discretely self-similar (DSS), a periodic fractal like behavior also referred to as echoing.Second, as a result of scale invariance, scalars like the black hole mass M BH on one side of the threshold and the Kretschmann scalar I = R abcd R abcd on the other obey a power law of the form as a function to the parametric distance to the threshold, labeled here by the critical parameter p ⋆ .In the DSS case, this scaling behavior appears with a superimposed ∆-periodic wiggle [2,3].Third, by evolving different families of data in spherical symmetry, Choptuik concluded that these features are universal.In particular, the scaling exponents λ in the power laws (1), as well as the DSS echoing periods ∆, were independent of the initial data.Notably, all the families of initial data, when evolved close enough to the threshold, share a common configuration, now referred to as the Choptuik spacetime.A spacetime with this specific symmetry, has since been proven to exist [4].These three features (selfsimilarity, power law scaling and universality) constitute what is now known as critical phenomena in gravitational collapse.Critical phenomena in gravitational collapse have since then been observed in various matter models, mostly but not exclusively in spherical symmetry (see [5] for a detailed review).Families of initial data with a varying parameter p are evolved numerically such that for small values of p, the evolution eventually leads to flat spacetime, whereas large values lead to horizon formation.Following this end-state classification one can bisect towards a better estimation of the threshold parameter p * delimiting the threshold of collapse.In this way one can measure the distance to the threshold by the precision with which p * is tuned.The more digits are known, the closer to the threshold the spacetime is and the more likely we are to observe critical phenomena.
Choptuik achieved an exquisite 13 (or more) digit tuning of p * by using a maximally constrained formulation of general relativity (GR) with zero shift and areal radius in an adaptive mesh finite differencing code.Subsequently this setup was studied with a variety of formulations of GR tailored to spherical symmetry.For instance, Garfinkle and Duncan [6] used null coordinates, and Martín García and Gundlach [7] constructed coordinate systems adapted to self-similar spherically symmetric spacetimes.There are a number of studies of the scalar field model with aspherical perturbations.It was demonstrated numerically that at the linear level there is only one growing mode solution [8].In pioneering work [9], non-linear axisymmetric perturbations of the Choptuik spacetime were studied, and their findings were confirmed much more recently by Baumgarte in [10] with a 13 digit tuning.In the latter, the BSSN formulation of the Einstein field equations was used with 1+log slicing and Γ-driver shift condition.Fully 3d non-linear perturbations have been studied in [11,12], understandably with far less tuning to the threshold.Of these, Deppe et al.'s paper [12] is of particular interest to us, because they use a similar setup to ours.They studied the collapse of massless scalar fields using the generalized harmonic gauge (GHG) formulation of GR using a pseudospectral approximation.Coordinates that gradually zoom into the center of the computational domain and a new gauge source function were employed.They managed to properly reproduce Choptuik's results, for the first time with a pseudospectral code, tuning up to 6 digits.
Choptuik's results are relevant to our understanding of cosmic censorship, as his spacetime should contain a naked singularity, but generic physical data are not spherical, and so assessing the generality of his findings without spherical symmetry is fundamental.The 6 digit level of tuning of [12] gives echoing periods and scaling exponents consistent with those of spherical symmetry; whereas the finer tuned data of [9,10] find that these scalars are dependent on the aspherical deviation of the initial data.To assess whether the different results are due to the difference in tuning, we need to improve the best level of tuning with GHG and a pseudospectral setup.
Although non-spherical models are generally more difficult to treat, an understanding is emerging.In vacuum [13][14][15][16][17][18], in the presence of electromagnetic waves [19,20] and for scalar fields as mentioned above, numerical evidence now consistently suggests a deviation from the spherical phenomenology.Confidence in some of these results, for instance whether exact DSS occurs beyond spherical symmetry, relies on the degree of fine tuning to the threshold, as for instance a minimum number of periods needs to be observed to assess DSS.There is however room for improvement.For example, in the challenging case of gravitational waves, the best tuning [15,16] is around 6 digits.It should be noted that for one of the two common families evolved in these works, the oblate centred Brill wave, the tuning in [16] performed with first order GHG was a digit further from the threshold than the one performed with BSSN in [15].
Here we pave the way to improved studies of the threshold of collapse without symmetry, focusing on a subset of the obstructions to tuning we found in [13,16] for Brill waves, but retreating to the classic spherical scalar field setup as a testbed.In those searches we had to stop tuning either because our code failed to stably evolve the fields until they dispersed or formed an horizon, or because our post-processing apparent horizon finder AHloc3d [21] could not locate the presumably present, but strangely shaped horizons that formed.Work concerning apparent horizon finders is ongoing and will be reported on elsewhere.In this paper we focus instead on solving the problems that prevented our code from reliably evolving extreme data.Some families of data were prone to develop undesirable coordinate features.A shift with a large gradient would be the most common sign of this instance.Coordinate singularities have been reported to be the main obstacle of other studies of critical collapse whether performed with BSSN or GHG.For BSSN evolutions this motivated the use of "quasi-maximal" slicing [14] and of the shock-avoiding coordinate conditions [17,[22][23][24].In GHG, the remaining gauge freedom lies in the choice of gauge source functions.The authors of [12] report that to avoid divergences, they needed to employ substantially modified gauge source functions as compared with the standard choice [25] for binary spacetimes.Our first goal here is therefore to look for a suitable gauge.Similar in spirit to [26], we look for a coordinate choice compatible with DSS.In particular we find and check a necessary condition for compatibility.We then construct a gauge source function that satisfies this condition.
Large constraint violations were present in the failed evolutions of most unclassified data.Even using adaptive mesh refinement to help control error, the closer the spacetimes were to the threshold, the more constraint violations we observed.For the examples considered here this is a particularly serious problem, probably since the first order formulation we employ necessarily introduces a large number of constraints.Our second goal was therefore to reassess and adjust the constraint damping parameters of GHG, particularly focusing on the first order reduction constraints in the context of collapsing spacetimes.These adjustments make a very significant improvement to our best tuning and indeed enable our code to display the well known phenomena first discovered by Choptuik.In section II, we present the numerical and physical setup of our pseudospectral code bamps, which was used to carry out all the present simulations.In section III we derive a DSS-compatibility condition, evaluate most gauge choices against it and suggest a compatible gauge source function.In section IV we present a mode analysis and the resulting improved constraint damping parameters for GHG.Finally, in section V we apply these modifications to simulations of massless scalar fields in spherical symmetry.In section VI we give a brief summary.

A. Code bamps
All the simulations we present in this paper were carried out with the pseudospectral code bamps [27].Its grid is divided into three main patches: a cubed box at the center of the simulation domain, a transition shell in the shape of a cubed sphere, and a spherical shell cir- cling the domain up to the outer boundary (see Fig. 1).Each of these patches is itself divided into cells, inside of which the evolved fields are stored on Chebyshev-Gauss-Lobatto collocation points.hp-adaptive mesh refinement (AMR) controls the resolution that every portion of the fields is granted by locally adapting the number of collocation points.It does so both by adding or removing points in a given cell (p-refinement), and by subdividing or joining cells (hrefinement).In order to evaluate the refinement needed, AMR can assess the error produced during the simulation and also the smoothness of the evolved fields.Details of the refinement scheme are given in [28].
Neighboring cells communicate with each other through characteristic fields using the penalty method.At the outer boundary, we use the constraint preserving boundary conditions described in [29] and [27].Symmetries are handled by the cartoon method, using double cartoon for spherically symmetric runs.
In bamps we employ a standard free evolution approach.The constraints are only explicitly solved for the initial data but not at any timestep during the evolution.The scalar field initial data is generated by an elliptic solve integrated within bamps, which uses the hyperbolic relaxation method [30].Time evolution is carried out using a fourth order Runge-Kutta method.The timesteps are adjusted during the evolution in order to preserve the Courant-Friedrich-Lewy (CFL) condition.

B. Physical System
In geometric units the trace reversed Einstein field equations read where g ab is the 4d metric, R ab the Ricci tensor, T ab the energy momentum tensor and T = g ab T ab its trace.Time evolutions are enabled by a 3 + 1 splitting of the 4d metric g ab .This gives the 3d spatial metric γ ij and a normal unit vector n a = α −1 (1, −β i ), where α is the lapse and β i the shift.In these variables the line element becomes We denote 4d component indices with Latin letters starting from a and 3d spatial ones starting from i.
We consider a massless real scalar field φ minimally coupled to the Einstein field equations.The corresponding scalar field energy momentum tensor is given by

C. Initial Data
In this work we consider two families of initial data, one starting at a moment of time symmetry, the other with a predominantly incoming pulse.The moment of time symmetry data is constructed from a scalar field with a Gaussian profile of the form and with a vanishing derivative along n a n a ∇ a φ = 0 .
For moment of time symmetry initial data we use R 0 = 0 as in [10].Observe that (5) implies that our A is a half of the equivalent parameter η used in [10].
For the incoming initial data we use the same initial scalar field as in [12].Thus we have and with R 0 = 5 and the appropriate limits taken at the origin.
To solve the ADM constraints we make a conformal decomposition of the metric and solve the extended conformal thin sandwich equations (XCTS) [31,32].We choose the conformal metric to be flat, γij = δ ij , and take its time derivative to vanish, ∂ t γij = 0. Furthermore we impose maximal slicing on the initial data, so that both the trace of the extrinsic curvature and its time derivative vanish, K = γ ij K ij = 0, ∂ t K = 0.After fixing these variables the XCTS equations become a set of coupled elliptic partial differential equations for the conformal factor ψ, the shift β i and the lapse α.
At the outer boundary, ∂Ω, we impose Robin boundary conditions compatible with a 1/r decay for ψ, β i and α.Concretely, the boundary conditions are given by where s i = γij ∂ j r/L is the spatial normal to the outer boundary, with the normalization factor L such that γij s i s j = 1.
For incoming initial data we have to solve the full coupled set of the XCTS equations.For moment of time symmetry initial data on the other hand, we only have to solve one equation for ψ.Thanks to the choice in Eq. ( 6) the solution for shift and lapse are trivially given by β i = 0 and α = 1.The remaining equation for ψ has the form where π here is the mathematical constant not to be confused with the field π introduced later.We use the hyperbolic relaxation method [30] implemented in bamps to solve the XCTS equations or Eq. ( 12).

D. Evolution Equations
a. Scalar field evolution equations The equation of motion of a massless scalar field is given by the Klein Gordon equation ∇ a ∇ a φ = 0. We work under a first order reduction with π as the time reduction variable n a ∂ a φ, and χ i as the spatial reduction variable associated to the reduction constraint S i := ∂ i φ − χ i = 0.This yields the first order massless Klein Gordon evolution system where σ is a damping term and the spatial connection is denoted by (3) Γ i jk .
b. GHG evolution equations In bamps, the Einstein field equations are formulated with the first order reduction of the GHG formulation [33].The evolved variables are the metric components g ab , the time reduction variable Π ab corresponding to −n d ∂ d g ab and the spatial reduction variable Φ iab associated to the reduction constraint C iab := ∂ i g ab − Φ iab = 0.The evolution equations are The canonical constraint damping parameters are γ 1 = −1 and γ 4 = γ 5 = 1/2 as introduced in [27].We use αγ 0 = 2 for incoming data and αγ 0 = 4 for the moment of time symmetry data.In section IV we examine the role of γ 2 , directly analogous to the damping parameter of the scalar field σ, in relation to the spatial reduction constraint.The harmonic constraint corresponds to where H a is the gauge source function, which we are free to specify (see section II E).As pointed out in [33], {C a , ∂ t C a } = {0, 0} encodes the Hamiltonian and momentum constraints, thereby ensuring both that the coordinates are harmonic and that the Einstein field equations are satisfied.The system of equations ( 16)-( 18) is symmetric hyperbolic, ensuring a well-posed initial boundary value problem provided suitable initial and boundary values.

E. Gauge source functions
Harmonic coordinates satisfy The GHG formulation generalizes this gauge to a less restrictive condition with the addition of a source term H a to the right hand side of (20).This generalization is satisfied as long as the harmonic constraint (19), C a = 0, is satisfied.All the gauge freedom of GHG is contained in H a .When the harmonic constraint is satisfied, the expressions of Γ a in standard 3 + 1 formalism give the following evolution equations for the lapse and shift where d t = α n a ∂ a = ∂ t − β k ∂ k and the 3d contracted spatial connection is (3) Γ i = (3) Γ i jk γ jk .To guarantee symmetric hyperbolicity, the arbitrary function of spacetime H a cannot contain derivatives of the fields.In that sense the gauge freedom of GHG is rather limited.We assess here in detail the choice of Harmonic Damped Wave Gauge (HDWG) as well as two modifications of it used in studies of critical collapse.
a. HDWG As explained in [25], the harmonic condition (20) can lead to strong gauge dynamics, when the physical degrees of freedom are very dynamical.An approach to combat this effect is to reduce those gauge dynamics, with the time and spatial coordinates subject to slightly different considerations.For instance, to suppress gauge dynamics associated with the spatial coordinates, the shift can directly be damped through the gauge source function yielding a shift condition very similar to the Γ-driver one (57).The equations for harmonic slicing and a shift damped by a factor η in GHG are This choice of H a only damps the spatial gauge dynamics, but an extra term is necessary in order to damp the temporal gauge dynamics.Besides, the authors of [25] report a concern about the growth of α −1 √ γ in some simulations, where γ = det(γ ij ) is the spatial volume element. Since in order to exponentially suppress α −1 √ γ.This is the well known HDWG for GHG (see equation (A15) in [25]) with a canonical choice of coefficients being η L = 1 and η S = 2.Although the HDWG choice of H a was extremely successful yielding black hole simulations for the authors of [34], in critical collapse studies, we have found that it fell short to handle the spacetimes of interest.Some modifications of this function have shown to be more appropriate for critical collapse studies, as we present below.b.HDWG-α 2 H a In [13], equations ( 1) and ( 5), and more recently in [16], equation ( 4), we gave a variant of the HDWG gauge source function characterised by the powers of the lapse as These modifications take into account the collapse of the lapse in spacetimes of extreme curvature (see a description in section IV) so that relevant terms survive despite a vanishing lapse.In the shift evolution equation one can see that the damping of the spatial gauge dynamics takes place whether or not the lapse collapses.This success inspired one of the modifications we present below, the adjusted damping parameter (see IV). Large values of η S had to be used to deal with undesirable coordinate features for some families of data in [16].Empirically we found η L = 1 and η S = 2 to be a reliable choice.c.HDWG-ln(α) H a In [12] the authors introduce a new version of the HDWG gauge source function in order to further enhance the suppression of α −1 √ γ and to control the growth of α −1 .
where R(t) is a roll-on function and W (x i ) is a spatial weight function.R(t) allows to start simulations with maximal slicing and then move on smoothly to their modified HDWG whilst W (x i ) enforces pure harmonic gauge at the boundary.As we are only interested in testing the impact of the modified HDWG, in our tests we consider R(t) = W (x i ) = 1.The authors of [12] also mention that for reliable evolutions of critical collapse, the powers in these logarithms are higher than the ones typically needed in binary black hole evolutions.

III. DSS-COMPATIBLE COORDINATES
As is hinted in the previous subsection, choosing a gauge for critical collapse simulations is not trivial.In GHG, gauge source functions have been constructed by careful experimentation and used to control divergences and a collapsing lapse.Even then, code failures coinciding with undesirable coordinate features such as a shift with a large gradient, occur.Ultimately this may result in a coordinate singularity.In practice, since these coordinate features trigger large constraint violations they tend to result in a crash before unambiguously determining the presence of such a singularity.Nevertheless, our best hypothesis thus far is that at least a subset of our code crashes are due to these undesirable coordinate features.Therefore we look for a gauge choice that keeps the lapse and shift "well behaved" during the simulation by reducing any additional coordinate features.
A natural suggestion is to tailor coordinates to the symmetries of the studied system.In the case of DSS spacetimes, using DSS-adapted coordinates would have the key advantage of directly tackling the problem of coordinate singularities since, if one could compute through one complete period, subsequent periods ought to follow without extra gauge features, as long as numerical error were sufficiently well controlled.From the numerical point of view such adapted coordinates are therefore attractive because they should potentially reduce the computational cost to reach a desired error.In [26], the authors study coordinates based on the ADM equations and on the homothetic Killing vector describing the continuous self-similarity (CSS) or DSS present in critical collapse spacetimes.With a particular boundary prescription, they test two gauges adapted to self-similarity in spherical symmetry.They call these coordinates symmetry seeking.Unfortunately neither is a simple modification of popular dynamical gauge conditions now known to be well behaved when treating a range of generic 3+1 dimensional data.Therefore, with the same general motivation, we take a different approach.Starting from the expression of the assumed symmetry-DSS in this caseon the evolved variables, we derive a necessary condition for gauge choices to be compatible with DSS (CSS would be similar).This condition is not restricted to spherical symmetry and can be straightforwardly used with any 3+1 formulation of GR.

A. DSS-compatibility condition
a. Periodic rescaling By definition [5], there exist coordinates (T, x i ) in which the metric of a DSS spacetime displays a periodic conformal rescaling of the form where the conformal metric gab is periodic, with fixed period ∆ gab T + ∆, Such coordinates are referred to as DSS-adapted coordinates.Slow time is a time coordinate adapted to DSS.
It can be thought of as the logarithm of an ever decreasing spacetime scale.An example of slow time is the one computed as the logarithmic distance of proper time τ to the accumulation time τ * along a future directed timelike curve that terminates at the accumulation point.(In general the coordinate then needs to be suitably extended away from the observer).In (T, x i ) coordinates, the characteristic time scale is ∆, therefore it would be numerically beneficial to evolve the spacetime using T , avoiding then the need to resolve ever decreasing time scales.Obviously, combining (35) and (36) we obtain that, after a period ∆, the 4d metric is rescaled as This transformation of the metric can be translated to the 3+1 variables in (3) as The transformation after a period ∆ of other relevant quantities can be derived from these, giving b. Coordinate choice Let us briefly refer to ( 39)-( 44) holding at all times as (i).A consequence of these, in particular of ( 40)- (41), is that (ii) hold at all times too.Thus (ii) is necessary for (i), and in turn for coordinates to be adapted to DSS.It is not a sufficient condition.(ii) does not imply (i) because the time integration required in that step is also affected by the remaining degrees of freedom of GR.
In the 3 + 1 split, we choose coordinates by imposing evolution equations for the lapse and shift where the right hand-sides are functionals permitted to depend also upon derivatives of the metric.Given this gauge choice, (ii) implies that (iii) holds, with the arguments assumed to respect such a symmetry according to (39)- (41), namely (i) holds.(iii) is a necessary condition for (ii), because without (iii) it cannot be the case that (ii) holds for the given gauge choice.By transitivity, (iii) is also a necessary condition of (i).We call (iii) the DSS-compatibility necessary condition for an arbitrary choice of gauge evolution equations to render the associated coordinates DSS-adapted.
Since (ii) is not sufficient to have (i) DSS-adapted coordinates, neither is (iii).Note that (iii) is not even a sufficient condition for (ii) to hold, as without additionally assuming (i), (iii) alone could not yield (ii) One might hope that satisfying the DSS-compatible condition would at least remove undesirable gauge features emerging from the incompatibility between the gauge choice and the self-similarity.This in itself would also be a good step towards avoiding coordinate singularities, the best step being of course actually evolving in adapted coordinates as explained at the beginning of this section.
Evolving a spacetime in DSS-adapted coordinates is not needed to assess the presence of exact DSS, as this can be verified after the evolution through coordinate transformations from the evolved coordinates to constructed DSS-adapted ones (see for instance section V C).However, it would not only be desirable from the numerical point of view (potential efficiency and coordinate singularity avoidance), but also because the dynamical construction of DSS-adapted coordinates would be the least ambiguous demonstration of such a symmetry.
In practice we can now use conditions ( 49)-(50) to assess whether the different gauges that are commonly used in critical collapse studies may be compatible with the symmetry we expect at the threshold of collapse.We refer to gauge choices satisfying the DSS-compatibility condition (49) and (50) as DSS-compatible.If a gauge is not DSS-compatible, the associated coordinates will not be DSS-adapted.
B. Moving puncture gauge a. Bona-Masso slicing condition Most finite difference simulations using the BSSN [35][36][37] or conformal Z4 formulations [38][39][40] choose the Bona-Masso [41] family of slicing conditions Harmonic gauge corresponds to f (α) = 1.The use of ( 40) and ( 44) gives which satisfies the DSS-compatibility condition (49).More interesting is to consider 1 + log slicing, which corresponds to f (α) = 2 α , as it is the most popular choice.After a period ∆, we obtain which does not satisfy the condition (49).Similarly, shock avoiding slicing, corresponding to for a given constant κ, also fails to satisfy the DSScompatibility condition b. Γ-Driver shift condition The conformal nature of the BSSN and conformal Z4 formulations plays a role in the choice of shift condition.These formalisms are expressed in terms of a conformal spatial metric γij = χγ ij where the conformal factor is proportional to the determinant of the spatial metric χ = γ −1/3 .The Christoffel symbols associated to the conformal spatial metric Γi jk can be contracted with this metric to give (3) Γi = γjk(3) Γi jk .Using ( 39)-( 44) we obtain In terms of this connection, the Γ-driver condition reads as with µ S a constant.We find that this choice does satisfy the DSS-compatibility condition (50) We conclude that the most common gauge choices used in finite differencing codes will not provide coordinates adapted to DSS due to the slicing condition, with the exception of pure harmonic gauge.

C. Gauge in GHG
Since in GHG the gauge freedom is fixed through the gauge source function H a , the DSS-compatibility condition (49)-(50) can be rewritten in terms of H a .Using (39)-( 44) in ( 21)- (22) gives So requiring (49)-( 50) is equivalent to requiring The HDWG gauge source function does not satisfy the DSS-compatibility condition (61), since Similarly, neither variant of HDWG, HDWG-α 2 (29) nor HDWG-ln(α) (32), satisfies (61).We have instead that Eq. ( 29) and that Eq. ( 32) We conclude that none of these popular dynamical gauge conditions that have recently been used in studies of critical collapse satisfy the DSS-compatibility condition.

D. DSS-compatible gauge sources
Remaining in the context of GHG from now on, the DSS-compatibility condition (61) can be taken into account to construct gauge sources designed to satisfy it.As can be seen in ( 63), the variant HDWG-α 2 H a only fails to satisfy (61) in the factor proportional to η L , for a = 0.For a = i, the extra α −1 , which was initially added to counteract the collapse of the lapse in the shift evolution equation of (26), cancels out the factors e −2∆ that appear in γ ij after a ∆ period.However,  27), ( 33), ( 30) and (66) with fixed K = 10 and γ = 9.The coefficients ηL for HDWG and HDWG-α 2 all taken to be 1 because it is the canonical choice and it provides a fair comparison with HDWG-ln(α).We show the slicing condition corresponding to the DSS-compatible gauge source function with its ηL being 1 for the comparison, and also 4, a reliable value.
the same extra α −1 added to the first term of ( 29) with the same intention but regarding this time the slicing condition (21) fails to satisfy the DSS-compatibility condition and needs to be adjusted.
As mentioned in Sec.II E, the authors of [25] and [12] point out the importance of including natural logarithms of √ γα −1 in the slicing condition to prevent divergences.
Therefore we keep the logarithms, but suggest to balance out the e ±∆ factors inside the logarithm, in order to satisfy (61).This leads to such that This is far from a unique choice to fulfill the compatibility condition.One way to assess whether this choice is reasonable is to compare it with existing working gauge source functions.Fig. 2 shows, for specific values of the trace of the extrinsic curvature K and determinant of the spatial metric γ, the profiles of d t α as functions of the lapse α itself.In this figure we can see the different ways in which variations of HDWG deal with a near-vanishing lapse.For HDWG one might get the impression that d t α is always negative, but in fact it becomes positive, in a small range near α = 0 and the maximal value is also vanishingly small compared to the other gauges.That means that the HDWG gauge is counteracting the collapse of the lapse less strongly and also only at smaller α.The d t α corresponding to HDWG-ln(α) has a larger maximum and becomes positive at a larger α, intuitively resulting in pushing α towards larger values in an evolution, typically to the region in which d t α changes sign.In the unlikely case of the lapse somehow dropping to zero exactly it would also stay there.Very close to α = 0 it would be pushed away from zero, but initially only very slowly.The d t α corresponding to HDWG-α 2 instead diverges as α approaches zero, an effect of the lack of α 2 in (30).This gauge will therefore push very small α more aggressively to the equilibrium point with d t α = 0.Both modifications of HDWG successfully fight a vanishing lapse at least until it becomes too small to be dealt with.The DSS-compatible suggestion (66) keeps a factor α in the first term so it results in a very similar profile to the d t α from HDWG-ln(α).That can be better seen in Fig. 2 with the curve corresponding to the DSScompatible gauge and η L = 4. Coincidentally, that is also the value we found provided successful simulations.
In section V we present numerical critical collapse evolutions using this gauge with the choice of η L = 4 and η S = 6.Assessing how well these coordinates adapt to DSS can only be done once the self-similar phase is successfully reached.To ensure we can stably evolve our simulations up to that point, we next treat a second major obstruction, namely constraint violations.

IV. THE REDUCTION CONSTRAINT DAMPING SCHEME
Ensuring that the constraints of the system are well satisfied throughout the evolution guarantees that the numerically evolved data does truly represent the mathematical and physical problem we set out to evolve.Since numerical error cannot be completely avoided however, constraint violations will still occur.What is needed is a strategy to prevent them from dominating a simulation and potentially even causing the code to crash.
As suggested in [42], some systems of equations admit a damping scheme that reduces these violations by making the constraint surface (in solution space) an attractor.This was done explicitly for the Z4 and GHG systems in [33,43] by adding to the evolution equations terms containing multiples of the constraints and freely specifiable parameters.In ( 16)- (18) these are the terms with γ 0 , γ 1 and γ 2 .
As mentioned in section I, we are primarily concerned with violations of the reduction constraints, as these are not present in second order formulations, which yield a finer tuning of the threshold parameter and so are a point of suspicion.Consequently we restrain our analysis to a linear model for GHG whose only non-principal terms are those containing the reduction constraint, plus linear terms.

A. Mode analysis
Taking the first order GHG evolution equations ( 16)- (18) in vacuum, linearizing about a constraint satisfying background with perturbations that satisfy the harmonic constraint (19), working in the constant coefficient approximation and finally assuming that Π ab = Φ iab = 0 in the background gives where, here and for the rest of this section, we overload the notation so that g ab , and so forth, stand for the metric perturbation such that (69)-( 71) holds up to first order in the perturbation.The remaining coefficients α, β i and γ ij are the constant lapse, shift and spatial metric in the background.These simplifying assumptions suppress any coupling between different components, making the following analysis algebraically tractable.Here the analysis deviates from that done in [43] as the lapse and the shift are not taken to be 1 and 0 respectively but instead to be constant in the background.It differs from the analysis provided in [33] as we here take into account non-principal terms.Although it would be interesting to do so, dropping any of the remaining assumptions would make the computation very substantially more complicated and, since we are only trying to motivate a simple change to our constraint damping scheme, does not seem worthwhile here.
A mode analysis of the evolution system can tell us the rate at which the constraint violations grow in this approximation.Writing u = (g ab , Π ab , Φ iab ) T , the system of equations ( 69)-(71) can be written as Here A k correspond to the principal part matrices of GHG, whereas Bu corresponds to the subset of nonprincipal terms that survive linearization and our simplifying assumptions.Observe that the 4d indices ab can be omitted for the analysis.
Moving to the frequency domain, we make a mode Ansatz ũ = ũ0 e st+iω k x k where ũ is the Laplace-Fourier transform of u, ω k = |ω|ω k is an arbitrary 3d-vector in the frequency domain, and we write |ω| = γ ij ω i ω j .The system (72) then takes the form of the following eigenvalue problem where Using the unit frequency vector ωi (ω i ωi = 1) and its orthogonal 2d projection operator q ij (ω i q ij = 0) we can further simplify the system with a 2 + 1 decomposition where capital Latin letters A, B denote a projection by q ij , running over the 2d plane orthogonal to ωk .Our 2 + 1 notation denotes components in the direction of ωi as β ω = ωk β k and 2d projected ones as β A = q A k β k .Any other vector or covector is decomposed using an analogous notation.In this notation the mode variables in the eigenvalue problem (73) become ũ = (g, Π, Φω , ΦA ) T .The matrix M becomes The eigenvalues of M are The imaginary parts of these eigenvalues represent the speeds of the system's propagation modes and the real parts capture exponential decay or growth.The expected speed of light is captured by s 1 and s 2 for |ω| = 1.More interesting to our analysis is the role of the parameters γ 1 and γ 2 in the evolution of the constraints.We construct the propagation modes v = l • ũ from the left eigenvectors l of M (see Appendix A).In the generic case β ω ̸ = 0, the propagation modes, whose speeds contain damping terms, are where for the last equalities we exploited the orthogonality between q A i and ωi in |ω|ω i gab .This tells us that some modes of the reduction constraints propagate as v s3,4 with speed s 3,4 , whereas some propagate as v s5 with speed s 5 .We can then confirm by direct computation that Ciab is damped directly by a factor −αγ 2 .Similar results hold in the special cases β ω = 0 and γ 1 = 0.
Despite the simplistic nature of this treatment it helps us to understand that there is a clear relation between moments where the lapse collapses and the growth of the constraints, as we explain in the next section.
B. Adjusted damping a. Collapsing lapse The mode analysis shows that violations of the reduction constraints, C iab , are damped at a rate −αγ 2 .In the limit where the lapse tends to zero, constraint violations are no longer damped and can start to grow due to numerical round-off.It is an empirical fact that, with many popular gauge conditions, 3. The value of the lapse at the origin is plotted as a function of time in simulations of moment of time symmetry scalar field initial data.The orange line corresponds to a supercritical simulation that stopped and an apparent horizon was found at the time we observe the lapse reaches 0. The blue dashed line corresponds to a subcritical simulation where, despite a large peak of curvature exactly when the lapse approaches 0, the fields manage to eventually disperse.
the lapse α collapses to zero in the presence of an apparent horizon.Although this is a coordinate dependent feature, it is so consistent that a collapsing lapse is often used to detect the presence of a black hole.To the best of our knowledge, there is no mathematical proof in which this is demonstrated for GHG, but we do observe such behavior numerically.Intuitively, this means that spatial slices stretch and wrap around a singularity, causing a rapid growth of the radial metric components whilst proper time gradually freezes inside the horizon, but advances on the outside.This time freezing effect near a singularity can even be desirable for certain simulation purposes, for instance to steer clear of the singularity.
Although pure harmonic slicing is only "marginally singularity avoiding" in the terminology of Alcubierre [23], our experience with GHG in bamps is that every time we have confirmed the detection of an apparent horizon with either harmonic or one of the gauge source functions discussed above, the lapse had collapsed almost to zero, as Fig. 3 shows.Interestingly, we observe the same tendency when the curvature is extreme, even in the absence of an apparent horizon.For instance in evolutions of subcritical data close to the threshold of gravitational collapse, we observe a nearly vanishing lapse at the time where the curvature invariant peaks.A crucial example of this are our studies of critical phenomena in gravitational collapse.
b. Adjusted damping parameters Given that the rate of damping of C iab is suppressed as the lapse collapses, and that in critical phenomena simulations, close to the threshold, the lapse consistently collapses, we hypothesize that this vanishing lapse results in the large FIG. 4. The L2 norm of the sum of all constraints (the "constraint monitor") with and without the adjusted damping parameter; in the left for a subcritical amplitude; in the middle for an amplitude that was not classified without the adjusted damping of Eq. ( 79), but could be classified as subcritical with the adjusted damping; and a supercritical amplitude.The data used in these plots are for configurations 4 and 3 in table I.
constraint violations we observe close to the threshold.This may even result in code failures before field dispersion or trapped surface formation, thereby obstructing the threshold parameter search.This could account for at least a subset of the shortcomings of first order GHG in the tuning of the threshold parameter, as there is no reduction constraint to deal with in second order systems.To test our hypothesis we adjusted the parameter γ 2 in the GHG equations ( 16)-( 18) such that the damping terms in s 3,4 and s 5 decouple from the lapse with typically c γ2 = 2 or 4.
In this way, the damping of the reduction constraint violation speeds can resist a collapsing lapse Essentially with this change the damping timescale is given in terms of coordinate, rather than proper time.c.Damped Constraints As shown in Fig. 4, using the adjusted damping reduces the constraint violation of both subcritical and supercritical evolutions as desired.Crucially, we can see that it holds the constraints low enough for a previously unclassified spacetime to be classified as subcritical, enabling the continuation of the parameter search.In section V we now demonstrate that this improvement helps in the observation of critical phenomena using first order GHG.

V. CRITICAL COLLAPSE RESULTS
We now present critical collapse results for two oneparameter families of spherical initial data for the mass-less scalar field.We treat these spacetimes as a (very challenging) testbed for the code, and as such assume some familiarity with the standard notions of critical collapse.We refer the reader to [5] for a comprehensive overview.The first family is the moment of time symmetry data, Eqs. ( 5) and ( 6), the second incoming initial data , Eqs. ( 7) and (8).All bisection searches to the critical amplitude follow the same procedure.We classify data as subcritical if it evolves until the fields are fully dispersed.Although our code supports a robust method for spherical black hole excision [44], we are presently only interested in knowing whether or not an apparent horizon forms, so do not treat the subsequent black hole evolution carefully.The code therefore eventually fails for all supercritical data.In spherical symmetry, the presence of apparent horizons can be determined by zero-crossings of the expansion, which is algebraic in our evolved variables [45].We only classify failed evolutions as supercritical if they display such a negative expansion.Starting with a sufficiently wide initial window, we proceed to bisect, evolving data whose amplitude is in the middle of the regime.If the evolution is subcritical we update the lower bound and it if its supercritical we update the upper bound.We proceed with the bisection until we reach data that cannot be classified.In all plots the slow time T p , defined in (37), is computed from the proper time at the origin as and the accumulation time τ * following Eq.( 23) in [10].We computed the period ∆ using Eq. ( 24) in [10], Similarly, we computed the proper length x p from the x-coordinate as

A. Constraint violations
Critical phenomena of a massless scalar field minimally coupled to GR in spherical symmetry have been very well studied and reproduced.It was therefore surprising to us that our initial automated bisections of moment of time symmetry data with the standard first order GHG constraint damping setup with HDWG-ln(α) H a (32), apparently tuning up to 15 digits, displayed at most one and a half echoes.Closely looking at the data, large constraint violations were causing crashes and noisy data caused incidental negative values of the expansion.The 15 digit tuning was not trustworthy, and the solution space had drifted towards the subcritical side by classifying as supercritical what was simply noisy data without horizons.The last reasonable estimation would only be after a 6 digit tuning, already giving indeed the same echo and a half we observed with the initially wrongly-tuned 15 digits.
We obtained comparably poor results with the incoming data, only reliably tuning 7 digits.Those were our best results after increasing the damping parameters and forcing more refinement (to no avail).For this reason the tables and figures mention different refinement indicators and values of γ 2 for the two different families of data.With our initial gauge source function HDWG-α 2 , Eq. ( 29), we obtained similarly large constraint violations and unreliable results far from the 3 echoes that [10] found.
By instead choosing our adjusted damping parameter of Eq. (79), immune to the collapsing lapse that occurs as curvature approaches the supercritical side, we were able to control those constraint violations.As Fig. 4 shows, close to the time of highest curvature (t ≈ 5 for those simulations), the constraints are damped much more effectively when using the adjusted damping factor.This improvement of the constraints near criticality has allowed us to confidently classify data on both sides of the threshold.From the subcritical side, it avoids crashes caused by large constraint violations, allowing the simulation to last until the fields disperse, as the central plot in Fig. 4 shows.From the supercritical side, we observe clear and smooth indications of apparent horizons in the evolution of the expansion.
It should be stressed that the adjusted parameter does not provide an improvement of the constraints throughout the entire evolution, as can be seen at t ≈ 2.5 in Fig. 4 or at later times.Instead, it makes a clear improvement in the key case where the lapse collapses and the spacetime starts displaying critical phenomena, at t ≈ 5 in this case.
Table I provides an overview over all configurations discussed in the text.The table shows that the adjusted constraint damping parameter improved the estimation of the critical amplitude from 7 to 9 digits for the incoming data family (compare configurations 9 and 10) and from 6 to 11 digits for the moment of time symmetry family (compare configurations 3 and 4).As seen in Fig. 5, this has allowed us to observe up to 3 echoes.Fig. 6 shows 3 echoes for both the moment of time symmetry data and the incoming data.The former becomes approximately self-similar slightly earlier, as its first echo is stronger than the incoming data's first echo, whilst its last one is incomplete compared to the incoming data.In this plot the incoming data has been shifted to the right to verify that, indeed, different initial data displays the same echoing near criticality with a period of ∆ ≃ 3.44.We conclude that the adjusted damping was essential to see critical phenomena with GHG.

B. Coordinate behavior
Once our bisections reach spacetimes that are close enough to criticality, we can evolve those up to the self-similar phase and evaluate the utility of the DSScompatible gauge source function.
For this comparison we employ the same refinement strategy for both sets of data using the smoothness indicator for h-refinement (see [28] for details).This provides low constraint violations with all gauge combinations but does actually display a lower peak in the third echo for the moment of time symmetry data (compare the lower panel of Fig. 7 with the straight line in Fig. 6).The self-similar phase is well approached for all combinations once the adjusted damping parameter is used.As can be seen by comparing configuration 4 and 5 in table I the values of the critical regime do vary within a family with different h-refinement strategies after the 8th digit.
Table I shows the intervals that contain the critical amplitude A from the bisections.The values of A sub correspond to the curves in Fig. 7.These runs have been checked and are not always the raw output of the automated bisections.In particular, like in the case discussed in the previous subsection, the automated bisection of the moment of time symmetry data evolved with HDWG-α 2 gauge source function (configuration 2) gave a misleading 14 digit tuning that actually showed no improvement in the echoing and large constraint violations after the 12th digit.We therefore trust the bisection up to the 12th digit only.
It is expected that the use of different coordinates not only affect how close we can evolve to the threshold but it also can drift the solution space slightly because of different numerical error.For example in table I the critical intervals obtained with HDWG-α 2 (configuration 8)  On the left we show echoes of the highest subcritical amplitudes of two moment of time symmetry data bisections: one with the original damping parameter (configuration 3) and one with the adjusted damping parameter of Eq. (79) (configuration 4).The lack of reliable tuning beyond 6 digits is an obstacle to see critical phenomena with the original damping scheme.On the right we show the same for the incoming initial data family (configuration 10 and 9).and the DSS-compatible gauge (configuration 10) differ in their 10th digit but none is in principle more real.The true value could perhaps be more accurately obtained with careful convergence testing, but since round-off error may already contribute with this level of tuning, and we are interested here primarily in capturing the correct phenomenology, we make no attempt to do so.
We can use again Fig. 7 to assess from which point this difference in tuning stops being meaningful with respect to the observed critical phenomena.For both sets of data, the low level of tuning obtained with vanilla HDWG (of only 1 or 2 digits) had a clear impact on our capacity to observe critical phenomena as it does not give any echoes for either data type.The difference in tuning obtained with vanilla HDWG compared to the other gauge source functions is therefore a severe shortcoming.
Moving on to the other three gauge source functions, the lower panel of Fig. 7 displaying moment of time symmetry data shows that the better tuning obtained with the HDWG-α 2 variant does indeed give the highest last echo.There both the data evolved with HDWG-ln(α) and the DSS-compatible version behave in exactly the same way, so the difference in the tuning they provide, although not unphysical, is not significant.In the upper panel displaying incoming data, the best tuned data, also obtained with HDWG-α 2 H a , this time does not provide the best last echo.It comes third after the one obtained with HDWG-ln(α) H a and the DSS-compatible version.The improved tuning obtained with HDWG-α 2 is therefore not as meaningful as that obtained with the DSS-compatible gauge for this family of data.Again the constraint violations are all comparable so all critical in- tervals are valid, but it does seem that gauge differences affect how easily, after how many digits, we can evolve near criticality.
A sharp profile in the shift was particularly prominent when evolving vanilla HDWG with both sets of data.In order to examine coordinate behaviour more closely, we compare the last common amplitude run in bisections of all three different gauge source functions that approach criticality.For the incoming data this corresponds to evolutions of A = 0.108933281(456118).The evolution of this data with the DSS-compatible gauge shows a clear uninterrupted dispersion of the fields allowing to classify the spacetime as subcritical.In contrast, both adjustments to HDWG crash, classifying it as supercritical.A common feature to both of the HDWG-variant evolutions is the presence of what resembles a step function in the profile of the shift vector in space.The constraint violations are comparable enough to trust both classifications but we see that these sharp features in the shift seem to appear later in the bisection for this family of data with the DSS-compatible guage source function.This advantge does not seem to be the case for the last common amplitude in the moment of time symmetry data, A = 151675332(291050) which shows very similar profiles of the shift, none as smooth as would be desirable, for bisections with all three gauge source functions.This spacetime was classified as subcritical with HDWGln(α) H a whereas a horizon was found with the other two choices.
To summarize, the two HDWG-variants considered and the DSS-compatible gauge agree very well to to 9 digits for both sets of data with comparable critical results and constraint violations.With this level of tuning our results will be affected by round-off, so to tune further we would need to work with higher precision arithmetic.More work will be needed to completely avoid the gauge features that affect the bisections near criticality, but we nevertheless conclude that the DSS-compatible gauge source function maintains well behaved coordinates far into the critical regime.

C. DSS and gauge sources
As positive as it is to tackle the issues that were stopping our critical collapse bisections, the question of how well any of our gauges dynamically adapt to the approx- imate self-similarity of the evolved spacetimes remains.We compare directly with the HDWG-ln(α) non-DSScompatible gauge source function, because we built the DSS-compatible suggestion so as to approximately agree with it, as can be seen in Fig. 2, but also because they provide very comparable results for critical phenomena, as demonstrated in Fig. 7.
First, we examine the presence of DSS in the fields evolved.After the evolution we transform the evolved coordinates into constructed canonical DSS-adapted coordinates to verify the symmetry.These coordinates are slow time computed from proper time T p at the center, as defined in Eq. ( 37), and its spatial counterpart where x i p are proper lengths.In Fig. 8 we plot the scalar field of a near-critical spacetime evolved with HDWGln(α) and the DSS-compatible gauge sources, with respect to (T p , X i p ).This figure helps us also assess if any unwanted gauge features actually disrupt the selfsimilarity phase of the spacetime.This is luckily not the case for either gauge choice.It seems that any undesirable gauge features reported in the previous subsection V B are not sufficient to prevent reliable bisections, and not enough to disrupt the spacetimes that survive.The periodicity observed in these plots also confirms that we are close to the threshold, the only case where the similarity would be exact.
We proceed to examine whether our evolved coordinates are adapted to DSS.If our evolved time coordinate t was exactly T p , as defined in Eq. ( 37), then Eq. ( 40) would hold, which can be rewritten in the same way as Eq. ( 35) as where α is ∆-periodic in T p at fixed spatial similarity coordinate.Fig. 9 we plot the left hand side of Eq. ( 86) to check if indeed we find a periodic function of constant average as ln α(T p , X i p ) is.The periodicity is certainly noticeable for all gauge source functions, but the nonvanishing slope indicates that Eq. ( 86) is not satisfied for either gauge choice.This means that t ̸ = T p .
The DSS-adapted coordinates (T p , X i p ) are not unique.As it is explained in Appendix B, any coordinates (T, X i ) related to (T p , X i p ) by where f a T p , X i p is a function ∆-periodic in T p , will be themselves DSS-adapted.The existence of such a wide range of DSS-adapted coordinates could make us hope that even if our coordinate t does not exactly coincide with T p , it could still coincide with T in (87) and hence be DSS-adapted.However, a short calculation also provided in Appendix B shows that if t was DSS-adapted it would satisfy the equation where f must be a function ∆-periodic in T p .Fig. 9 shows that ln(α) + T p has some oscillations with a characteristic frequency, but it is not periodic in T p , which implies that t is not a DSS-adapted time coordinate.We can see more explicitly that t is not DSS-adapted at all.For our coordinate t to be DSS-adapted we should have that ln(α(t, x i )) + t = ln α(t, x i ) .
In Fig.  86) for the highest amplitude of bisections of incoming data on the top panel and moment of time symmetry data in the lower panel both with adjusted damping of Eq. ( 79) and both HDWG-variants Ha (configurations 10 and 8) and the DSS-compatible Ha, Eq. (65) (configuration 11).Despite the different data used, each curve seems identical in both panels, mostly determined by the choice of Ha.
Fig. 9: although the average through a constant value during the self-similar phase is much more noticeable than in Fig. 9, there is not such a clear periodicty.We can now conclude that the dynamical coordinates associated with neither gauge are adapted to DSS.Following the calculation in Section III we already knew that any non-DSS-compatible gauge choice could not yield DSSadapted coordinates, and that even a DSS-compatible one could not guarantee them.The gauge choice Eq. ( 65) is unfortunately not the needle in the haystack.
The choice we made for a DSS-compatible gauge source function is only a first attempt.Many other gauge source functions could be constructed that satisfy the DSScompatibility condition.As we already saw in Fig. 2 90) for the highest amplitude of bisections of incoming data on the top panel and moment of time symmetry data in the lower panel both with adjusted damping of Eq. ( 79) and both HDWG-variants Ha (configurations 10 and 8) and the DSS-compatible Ha, Eq. (65) (configuration 11).The top panel is zoomed out to show how the DSS phase stands out in the evolution.
qualitatively closest to HDWG-ln(α) gauge source function, which may suggest the reason that this variant performs the best of all the non-DSS-compatible choices.

VI. SUMMARY AND DISCUSSION
In this paper, using the classic example of scalar field collapse as a testbed [1], we have overcome two major obstacles that were stopping our critical collapse simulations of Brill waves in [16], namely constraint violations and undesirable coordinate features.
Taking into account the collapse of the lapse in nearcritical spacetimes, we modified the reduction constraint damping parameter of GHG, γ 2 , such that it keeps damping constraint violations near criticality.This can be thought of as changing the damping timescale from proper to coordinate time, and has allowed us to improve our critical parameter estimation and thus to clearly observe critical phenomena with two different families of real massless scalar fields in spherical symmetry minimally coupled to GR.Since pseudospectral codes typically employ first order systems of equations, GHG in our case, we suspected that the inefficient damping of violations of the reduction constraint of GHG in spacetimes where the lapse collapses was responsible for the poorer tuning provided by pseudospectral searches of critical collapse in comparison to those performed with finite difference codes and second order in space formulations of GR.With the improved damping scheme we can now confidently tune to a level competitive with bespoke spherical methods.We believe the bottleneck in tuning now to be the effect of round-off error.
Aiming to exploit the structure of a DSS spacetime to avoid coordinate singularities, we derived a necessary condition for a gauge choice to be compatible with DSS.The condition is appropriate also for CSS spacetimes.Using a specific gauge source function for GHG that satisfies the DSS-compatibility condition, we have reproduced the critical phenomena of a massless scalar field at the level of commonly used gauge choices.The spatial profiles of the fields suggests that this gauge source function respects DSS slightly better than the other gauge sources, particularly vanilla HDWG, which was not designed with critical collapse in mind (see [46] for discussion).However, at the level of tuning that can be achieved with double-precision arithmetic, we do not see much improvement in the observation of critical phenomena, constraint violation reduction nor in the avoidance of coordinate singularities with the DSS-compatible choice we made.This is more a practical than a principle issue.Every time we extend from one scale echo to the next, a larger class of coordinate singularities may occur if we do not satisfy the compatibility condition.Further work is however needed to come up with gauge source functions that give truly DSS-adapted coordinates.
The natural continuation of the present work is simply to apply our improvements to the more interesting setup of axisymmetry, both in the context of scalar fields and gravitational waves.We expect that an improved apparent horizon finder will be needed to classify spacetimes of the latter.By systematically treating each of the obstructions we encounter, we expect to be able to further investigate non-spherical spacetimes more accurately, with the ultimate aim of gaining a comprehensive understanding of the threshold of black hole formation in complete generality.corresponding to two sets of DSS-adapted coordinates x ā = t, x ī and x â = t, x î .
Let a spacetime enter the DSS phase and call echo 1 the event where the scalar field first reaches 0.6 and echo 2 the subsequent event where the field reaches again 0.6.In general we have t 2 , x i 2 = t 1 + δt, x i 1 + δx i with arbitrary δt and δx i whilst being in DSS-adapted coordinates implies δt = ∆ and δx i = 0.
Assuming that there is a functional relationship between the coordinates such that x ā = h ā (x a ) and x ā = f ā x â then t2 = h t t 2 , x i 2 = h t t 1 + δt, x i From Eq. (B2) we can conclude that f t t1 + ∆, x î 1 = f t t1 , x î 1 + ∆ and hence that the function must be ∆-periodic in t.Note that this is only the case because the DSS period is the same for both DSS-adapted coordinates, which is why it does not apply to the not-DSS-adapted coordinates as Eq.(B1) shows.Similarly, we have that x ī = fī t, x î is also ∆-periodic in t.From Eq. (B3) we have that where everything inside the square brackets, now referred to as Câ x b , is ∆-periodic in t.
As in Eq. ( 85), in these DSS-adapted coordinates we have ᾱ t, x ī = e − t α t, x ī , (B5) α t, x î = e − t α t, x î , (B6) where α t, x ī and α t, x î are ∆-periodic in t and t respectively.In terms of the conformal metrics gā b and gâ b, ∆-periodic in t and t respectively, we have that Eq. (B4) implies α t, x ī = −g āb ∇ āt ∇b t  where the right hand side should be ∆-periodic in t as Eq.(B7) shows.
For the coordinates in the main text, we computed slow time from proper time T p which corresponds to t and we would like to check whether our evolution coordinate t corresponds to another DSS-adapted coordinate here denoted as t.The function f in Eq. (89) corresponds to the right hand side of Eq. (B9).

2 4 FIG. 2 .
FIG.2.Time evolution equations of the lapse as functions of the lapse itself, for the various gauge source functions in Eqs.(27), (33), (30) and (66) with fixed K = 10 and γ = 9.The coefficients ηL for HDWG and HDWG-α 2 all taken to be 1 because it is the canonical choice and it provides a fair comparison with HDWG-ln(α).We show the slicing condition corresponding to the DSS-compatible gauge source function with its ηL being 1 for the comparison, and also 4, a reliable value.
FIG.3.The value of the lapse at the origin is plotted as a function of time in simulations of moment of time symmetry scalar field initial data.The orange line corresponds to a supercritical simulation that stopped and an apparent horizon was found at the time we observe the lapse reaches 0. The blue dashed line corresponds to a subcritical simulation where, despite a large peak of curvature exactly when the lapse approaches 0, the fields manage to eventually disperse.
FIG.5.On the left we show echoes of the highest subcritical amplitudes of two moment of time symmetry data bisections: one with the original damping parameter (configuration 3) and one with the adjusted damping parameter of Eq. (79) (configuration 4).The lack of reliable tuning beyond 6 digits is an obstacle to see critical phenomena with the original damping scheme.On the right we show the same for the incoming initial data family (configuration 10 and 9).

FIG. 6 .
FIG.6.Echoes of the highest subcritical amplitudes of incoming initial data (configuration 10) and moment of time symmetry data (configuration 4) and the adjusted damping of Eq. (79).The incoming data has been shifted to the right to enable a direct comparison with the moment of time symmetry data.

6 FIG. 8 .
FIG.8.The scalar field is plotted as a function of slow time and the similarity coordinate in the x-direction Xp as defined in (84).These correspond to the highest subcritical evolutions of incoming data with the adjusted damping of Eq. (79) and on the left with HDWG-ln(α) (configuration 10) and on the right with the DSS-compatible gauge source function (configuration 11).

TABLE I .
(18)ction intervals for various configurations.The critical amplitude A * is located in the interval [A sub , Asup].The identifier serves as reference for the configurations in the text.γ2 is the constraint damping factor for the reduction constraints, see Eqs. (16)-(18).Ha indicates the choice of gauge source function, which are introduced in Sec.II E and Sec.III D.