Type II critical collapse on a single fixed grid: a gauge-driven ingoing boundary method

We develop a numerical method suitable for gravitational collapse based on Cauchy evolution with an ingoing characteristic boundary. Unlike similar methods proposed recently (Ripley; Bieri, Garfinkle&Yau 2019/20), the numerical grid remains fixed during the evolution and no points need to be removed or added. Increasing coordinate refinement of the central region as the field collapses is achieved solely through the choice of spatial gauge and particularly its boundary condition. We apply this method to study critical collapse of a massless scalar field in spherical symmetry using maximal slicing and isotropic coordinates. Known results on mass scaling, discrete self-similarity and universality of the critical solution (Choptuik 1993) are reproduced using this considerably simpler numerical method.


Introduction
Critical phenomena in gravitational collapse are one of the most remarkable discoveries made through numerical methods applied to Einstein's field equations of general relativity. Since Choptuik's groundbreaking study of a massless scalar field coupled to the Einstein equations in spherical symmetry [1], similar phenomena have been discovered for a variety of matter models and even in vacuum, see [2] for a review article. Briefly, the idea is to choose a smooth one-parameter family of initial data such that in the future Cauchy development of such data, a black hole forms for large parameter values and the field disperses to flat spacetime for small parameter values. We are interested in the threshold between these two final states and the associated critical solution. In what has been termed Type II critical collapse, the black hole mass becomes infinitesimally small as the threshold, obeying a universal scaling law, and the critical solution is discretely self-similar and universal, i.e. independent of the particular one-parameter family of initial data chosen. (There is also Type I critical collapse in certain models, where the black hole mass is finite at the threshold and the critical solution is stationary or time-periodic.) It is this discrete self-similarity of near-critical evolutions that makes the problem so hard numerically: the solution repeats itself on smaller and smaller spatial scales, in shorter and shorter time intervals. Choptuik [1] implemented an adaptive mesh refinement algorithm [3] in order to be able to resolve the increasingly smaller length scales of the solution.
Alternate methods to tackle the same problems have subsequently been developed, e.g. formulations in double null coordinates with [4] and without [5] adaptive mesh refinement, although in the latter case grid points had to be added during the evolution in order to maintain accuracy.
Recently Bieri, Garfinkle and Yau [6] proposed a general method for Cauchy evolution in numerical relativity whereby the boundary of the finite spatial computational domain is expanded along a spacelike direction at each time step. Additional initial data must be specified on this surface. The advantage is that with such a setup, no outer boundary conditions need to be imposed because all the constant-time slices lie within the domain of dependence of the initial slice and the additional "tilted" spacelike surface. This proposal thus avoids the long-standing problem of imposing boundary conditions along a finite timelike surface in general relativity [7]. Other alternatives to this problem include Cauchy-characteristic matching [8], evolution on hyperboloidal slices compactified towards future null infinity [9,10,11] and the regular conformal field equations [12,13].
Related to Bieri et al. 's scheme is the "excision method" proposed by Ripley [14], whereby the computational domain is excised along a surface that is spacelike or tangent to an ingoing characteristic of the boundary on the initial slice. Again, no boundary conditions need to be impose at the outer boundary because in this case all characteristics leave the computational domain. This method appears to be well suited to gravitational collapse problems. A disadvantage is that grid points are lost during the evolution due to the excision procedure so that one will have to add grid points in the interior in order to maintain accuracy.
The method developed in the present paper is similar to Ripley's in that the outer boundary of the spatial computational domain is an ingoing characteristic. However, no grid points are excised; instead merely the spatial coordinates are changed as time proceeds.
In the Arnowitt-Deser-Misner (ADM) formulation of general relativity [15], the time vector field ∂/∂t is decomposed as 1 where n a = −α∇ a t is the unit future-directed timelike normal, α is the lapse function and β a the shift vector (which is spatial, n a β a = 0). This means that a point with spatial coordinates x i on the spatial slice at time t correponds to the point with spatial coordinates x i − β i dt on the slice at time t + dt if we drag it along the timelike normal. Consider now a shift vector field of the form where c is a constant w.r.t. the spatial coordinates x i . Hence identified points on the spatial slices will change coordinates according to as time increases from t to t + dt, so if we choose c < 0 then the coordinates "zoom in" isotropically towards the origin. The significance of (2) is that it is a homogeneous solution to the spatial isotropic gauge condition in spherical symmetry, Eq. (11) below, where β r = rβ so (2) corresponds to β = c = const. The value of the constant c will be fixed by the boundary condition on the shift in the isotropic gauge condition. For a suitable value of this (in general time-dependent) constant, the outer boundary can be made an ingoing characteristic (or spacelike) so that no boundary conditions on the evolved fields are needed. We supplement the isotropic spatial gauge condition with a maximal slicing condition. The advantage of such a slicing as compared with the polar slicing used by Choptuik [1] is that the coordinates remain regular at the apparent horizon when it forms, which allows for a more accurate determination of its location and mass.
This article is organised as follows. In Sect. 2 we set up our model problem of a massless scalar field in spherical symmetry, and we state the gauge conditions and their boundary conditions appropriate for our ingoing boundary method. In Sect. 3 we provide details on the numerical methods we use to solve the field equations. The numerical results are contained in Sect. 4. We set up two families of initial data, describe our method to tune to the critical parameter and to choose an appropriate outer boundary radius, and we present results on the mass scaling, discrete self-similarity and universality of the critical solution. In Sect. 5 we summarise, discuss potential challenges of the method and future applications.
2 Formulation of the model

Choice of gauge and variables
In spherical symmetry and isotropic coordinates, the spacetime metric takes the form ds 2 = −(α 2 −ψ 4 r 2 β 2 )dt 2 +2rβψ 4 dt dr +ψ 4 dr 2 + r 2 (dθ 2 + sin 2 θ dϕ 2 ) . (4) We impose maximal slicing and hence the extrinsic curvature has only one independent component in spherical symmetry: For reasons discussed shortly, we define a rescaled quantitŷ The massless scalar field Φ itself does not enter the equations but only its first derivativesξ where here and in the following a dash denotes a partial derivative w.r.t. r and a dot w.r.t. t. The fundamental variables α, β, ψ,K r r ,ξ andΠ depend on t and r only, and the powers of r in their definitions have been chosen so that they are all even functions of r with finite nonzero limits at the origin r = 0.

Field equations
The relevant components of the Einstein equations R ab = κ∇ a ∇ b Φ, where R ab is the spacetime Ricci tensor and κ = 8π in geometric units, are the momentum constraint rK r r + 5K r r + κπξ = 0 (8) and the Hamiltonian constraint When linearising (9) about a given background solution ψ, the coefficient of the undifferentiated term proportional to ψ is manifestly negative because of the negative powers of ψ in (9). If this was not the case then non-unique oscillatory solutions might exist; see [16,17] for further discussion of this issue. This is the reason for the choice of the powers of ψ in (6) and (7). The maximal slicing condition implies the following equation for the lapse, and preservation of the isotropic form of the metric (4) yields The equation of motion for the scalar field ∇ a ∇ a Φ = 0 reduces to the pair of first-order equationṡξ There are redundant evolution equations for ψ andK r r that can be used to monitor the accuracy of the code; the first will also be needed to specify boundary conditions:

Boundary conditions
A crucial feature of our method is the choice of gauge boundary conditions. We want the outer boundary to be ingoing null or spacelike, which corresponds to setting β .
with ν 1, where . = means equality at the outer boundary r = r max . For the results presented in Sect. 4 we will always choose ν = 1 corresponding to the boundary being null, although we will briefly discuss making ν a timedependent function in Sect. 5.
Since there are no ingoing characteristics at the outer boundary with this choice, the evolution equations (12) and (13) for the scalar field do not require any boundary conditions. We specify Dirichlet boundary conditions on ψ for the Hamiltonian constraint (9) by evolving (14) at the outer boundary. The momentum constraint (8) does not require a boundary condition as this is already fixed by demanding the solution to be regular at the origin.
What remains to be specified is an outer boundary condition on the lapse α for the maximal slicing condition (10). Freezing the lapse to its flat value α = 1 is not a good idea since this will lead to unacceptably large slice stretching as the physical size of the grid shrinks and the lapse collapses in the centre as the singularity is approached. Instead we simply advect the lapse along the shift at the outer boundary, as in the first terms of all the evolution equations: Another way of phrasing this is to extrapolate (in time) the value of the lapse at the outer boundary r = r max on the slice at time t + dt from its value at the identified radius on the slice at time t, which according to (3) is at (r + rβdt) r=rmax < r max (note β < 0 at r = r max ).

Evolution scheme
Given data at time t, we first evolve the scalar field variablesξ andΠ to the next timestep t + ∆t using (12) and (13). At the advanced time the radial ordinary differential equations (ODEs) (8)- (11) are solved in this order for K r r , ψ, α and β (notice they form a hierarchy). Dirichlet boundary values for ψ and α are supplied by evolving (14) and (17) at the outer boundary along with the other evolution equations, and the boundary condition for β is (16).

Discretisation
We use a fixed non-uniform radial grid at points is a cubic map from numerical to physical coordinates. Here r eff r max can be thought of as an "effective" radius the grid would have if the same resolution as close to the origin was used all the way to the outer boundary. We typically choose r eff ≈ 1 2 r max . It should be noted that a non-uniform grid is not essential for our method to work, it just saves computational resources since the distribution of grid points is better adapted to the features of the solution, which has its largest gradients close to the origin. We could just as well take r eff = r max corresponding to a uniform grid. With respect to the numerical coordinate x, the grid is equidistant and staggered at the origin: We use N = 500 grid points for the simulations presented in Sect. 4. The grid remains unchanged during the evolution. We use centred fourth-order finite differences to discretise the equations in r. Near the origin the finite-difference stencils are modified according to the known (even) r-parity of all the evolved variables. Near the outer boundary (fourth-order) backward finite differences are used.

ODE solvers
Following the method of lines, the evolution equations are integrated forward in time using a standard fourth-order Runge-Kutta method. Sixth-order Kreiss-Oliger dissipation [18] is added to the right-hand sides of the evolution equations in order to maintain numerical stability (a small coefficient ≈ 0.1 is found to be sufficient).
The radial ODEs are solved using a direct band-diagonal solver at each substep of the Runge-Kutta method.
Since the size of the metric functions α, β and ψ changes drastically during the evolution, it is important to adapt the size of the time step ∆t in order not to violate the Courant-Friedrichs-Lewy (CFL) condition for numerical stability. At each time step, we compute the characteristic speeds of the scalar wave and set the time step size according to .

Termination criteria
We terminate a simulation when either a black hole forms (i.e. the evolution is supercritical ) or the field disperses to flat spacetime (i.e. the evolution is subcritical ). Formation of a black hole is detected by looking for an apparent horizon (outermost marginally outer trapped surface). This is an r = const surface whose outgoing null expansion vanishes, where is the areal radius and a is an outward-pointing radial null vector. In our variables (21) is equivalent to The radius r AH of the apparent horizon is the largest zero of this equation, and the associated mass is It is this mass computed from the apparent horizon that will enter the scaling law in Sect. 4.3. Assuming cosmic censorship holds, formation of an apparent horizon implies the existence of an event horizon containing the apparent horizon in its interior. We consider an evolution to be subcritical if the maximum (w.r.t. r) of the scalar curvature drops below some fraction (typically 5%) of its maximum value attained during the evolution. Fig. 1 Penrose diagram of a typical near-critical spacetime. Shown are the initial spatial slice at t = 0 and a number of subsequent spatial slices and their ingoing null boundaries for two different initial boundary radii r 1 and r 2 as discussed in the main text. On the last slice the apparent horizon (AH) forms, which at later times converges to the event horizon (EH) of the black hole. Spacetime is close to the Type II critical solution roughly in the shaded region.

Initial data and bisection
We consider two very different families of initial data for the scalar field: (i) data that would be exactly ingoing in a flat metric (ψ = α = 1, β = 0), and (ii) data that are centred at the origin and initially at rest, We fix the parameters σ = 1 and (for the ingoing family) r 0 = 10, and we take the amplitude A as the critical parameter. For large values of A the solution forms a black hole whereas for small values it disperses. We use the bisection method to find an approximation to the critical amplitude A * .

Choosing the outer boundary radius
A typical Penrose diagram of a supercritical evolution close to the critical point is shown in Fig. 1. It becomes obvious from this diagram that the success of our method will depend on a good choice of the radius r max of the outer boundary on the initial spatial slice. If this is taken to be too large, r max = r 2 in Fig. 1, then despite the fact that the outer boundary is an ingoing characteristic, the apparent horizon forms at a very small radius compared to the radius of the outer boundary. We terminate the bisection scheme if the radius of the apparent horizon in the supercritical evolutions gets too small, say r AH < 0.01 r max , and start over with a smaller value of r max .
If on the other hand the initial boundary radius is chosen too small, r max = r 1 in Fig. 1, then the field escapes from the spatial domain before the apparent horizon forms. In a numerical evolution of this type we observe that the bulk of the scalar field moves out of the domain but the scalar curvature (25) remains large, unlike in a subcritical evolution. If this happens, we terminate the bisection scheme and repeat it with a larger value of r max .
Essentially this adds an outer bisection loop (for r max ) to the inner one (for A). In practice, one does not have to repeat the A-bisection all the way from the start because one can use a somewhat smaller A-interval of the previous r max -iteration as the initial interval for the A-bisection at the improved value of r max .
Using this procedure we determine r To provide some idea of how the physical size of the grid changes during a simulation, we plot in Fig. 2 the areal radius R max of the outer boundary as a function of the number of time steps n for a near-critical evolution. The exponential decrease of R max with n is well adapted to the expected discrete self-similarity of the critical solution, which repeats itself on smaller and smaller scales. Also shown in Fig. 2

Mass scaling
In Fig. 3 we plot the apparent horizon mass M vs. the distance A − A * to the critical amplitude for a series of supercritical evolutions. In a doublelogarithmic plot this forms a straight line with a periodic wiggle: as first observed numerically in [1,19] and predicted from a perturbative analysis of the critical solution in [19,20]. According to this analysis, the period of the function Ψ is related to the echoing exponent ∆ (cf. Sect. 4.4) via ∆ = 2 γ. The values of the mass scaling exponent γ and the echoing exponent ∆ obtained from a fit to our numerical data are shown in Table 1 and are in good agreement with the predicted values. It should be noted that ∆ can be determined more accurately from the echoing behaviour of the near-critical solution (Sect. 4.4).

Discrete self-similarity and universality of the critical solution
Variables that are scale invariant display discrete self-similarity in near-critical evolutions. One such scale-invariant variable for the scalar field is (30) Fig. 4 Discrete self-similarity: the scale-invariant variable X (i) (τ, ρ) for a near-critical evolution of the ingoing family is plotted as a function of ρ at two different times τ = −1.8 (left) and τ = −3.1 (right) as a solid curve with dots at every tenth grid point. In the same plots, we also show X (i) (τ − ∆, ρ − ∆) as a function of ρ with plus symbols (+) at every tenth grid point. The echoing exponent is taken to be ∆ = 3.44.
Discrete self-similarity is best described in terms of logarithmic coordinates where T * 0 is the accumulation time of the critical solution. The conjecture, first discovered numerically in [1], is that for the critical solution (indicated by the star), any scale-invariant variable such as X (30) obeys where ∆ is the echoing exponent. In Fig. 4 we plot X (i) (ρ, τ ) for a near-critical evolution of the ingoing family as a function of ρ at two different times τ , and we overlay X (i) (ρ − ∆, τ − ∆) in the same plots, using ∆ = 3.44. The accumulation time T * 0 has been determined by minimising the norm of the difference between both functions at one fixed time τ . The fact that the curves nearly coincide provides strong support of the echoing property (32). We can also see in Fig. 4 that the solution is well resolved numerically both at the original time τ and at the time of the echo τ − ∆, when the spatial scale has shrunk by a factor e ∆ ≈ 31.
Finally we investigate if the critical solution is universal, i.e. independent of the particular one-parameter family of initial data. In Fig. 5 we again plot X (i) (ρ, τ ) for a near-critical evolution of the ingoing family as a function of ρ at two different times τ , but this time we overlay X (ii) (ρ − δ, τ − δ) for a near-critical evolution of the centred family, where δ is an overall familydependent scale chosen such that the norm of the difference between the two solutions is minimal at one fixed τ . The fact that the curves nearly coincide also at a different time τ with the same constant offset δ strongly supports the conjecture that the critical solution is universal, as already argued in [1]. Fig. 5 Universality: the scale-invariant variable X (i) (τ, ρ) for a near-critical evolution of the ingoing family is plotted as a function of ρ at two different times τ = −1.8 (left) and τ = −3.1 (right) as a solid curve with dots at every tenth grid point. In the same plots, we also show X (ii) (τ − δ, ρ − δ) for a near-critical evolution of the centred family with plus symbols (+) at every tenth grid point. The same (family-dependent) constant offset δ = 0.245 is used in both plots.

Discussion
We presented a numerical method for gravitational collapse based on Cauchy evolution with an ingoing null boundary. The method is similar in spirit to the excision method of Ripley [14] but differs in that no grid points are removed from the computational domain; rather, the grid remains fixed and only the coordinates are adapted along with the evolution. This is achieved by adding a linear term to the shift vector that causes the coordinates to "zoom in" isotropically towards the centre. This linear term is a homogeneous solution to the isotropic spatial gauge condition. Another important ingredient is the treatment of the lapse function. We propose to use an advection equation for the lapse along the shift vector at the outer boundary in order to provide boundary values for the slicing condition (in our case, maximal slicing). This corresponds to interpolating the lapse from the previous time step and minimises the amount of slice stretching as the lapse collapses towards zero in the high curvature region in the centre.
We worked out the method in detail for the model problem of a massless scalar field coupled to the Einstein equations in spherical symmetry. Known results on critical behaviour [1] are reproduced: the mass-scaling relation including its fine structure [19,20], the discrete self-similarity (echoing) of the critical solution and its universality among different families of initial data. This demonstrates that the method is well suited to studying critical phenomena in gravitational collapse, while being considerably simpler than more commonly used methods that typically employ adaptive mesh refinement.
A price one has to pay for the simplicity of the method is that the outer boundary radius r max of the initial data slice needs to be chosen carefully so that a sufficiently large region of spacetime where the evolution is close to the critical solution can be explored. We optimised r max using an outer bisection loop depending on the outcome of the standard inner bisection along the critical parameter. One might wonder if this makes the method overly computationally expensive. Certainly in spherical symmetry this is not the case as a single evolution takes less than five minutes on a laptop even close to the critical point. Furthermore one does not have to restart the bisection for the critical parameter from the beginning for each r max iteration; instead, a smaller interval from the previous bisection can be used as an improved initial guess. Whether the method is competitive in axisymmetry or without any spacetime symmetries remains to be seen.
We have tried to alleviate the need for fine-tuning r max by equipping the algorithm with a control system similar to the one described in [21]: make ν in (16) time dependent and steer it so that a typical feature of the solution such as the minimum of the outgoing expansion (21) remains approximately at a constant coordinate radius. The larger the value of ν, the stronger the magnifying effect. For this to work, r max must be chosen somewhat larger than its optimal value for a null boundary, and ν must be taken somewhat larger than 1 initially, so that the control system has enough room to do its job. While performing reasonably well at early times, we have found such a control system to be ineffective in halting the rapid escape of the scalar field from the domain that often occurs just before an apparent horizon forms if the initial r max was chosen too small or the control system kept ν too large for too long a time. One should note that ν must not get smaller than 1, otherwise the boundary becomes timelike and boundary conditions for the evolved fields are needed.
This has been used in much numerical work, including the first study of critical behaviour in vacuum axisymmetric gravitational collapse by Abrahams and Evans [22], as well as e.g. [17,23,24,25]). The quasi-isotropic gauge condition admits homogeneous solutions analogous to the isotropic gauge condition in spherical symmetry, and our method can be carried over with very few modifications. Work along these lines is in progress. It is conceivable that our method can be made to work with other classes of spatial gauge conditions as well. Any elliptic shift condition such as the minimal strain or minimal distortion conditions [26] requires boundary conditions, and the freedom in choosing the boundary data can be used to make the outer boundary an ingoing null surface. Evolutionary shift conditions such as the hyperbolic Gamma-driver condition employed in some of the first successful binary black hole merger simulations [27] require initial conditions, and they could also be modified by adding lower-order terms, which could be used to a similar effect. These are interesting questions for further research.
Finally it should be stressed that this ingoing boundary method or the related method of Ripley [14] are not limited to studying critical collapse.
One can also start with a standard Cauchy evolution with timelike boundary (where of course boundary conditions must be imposed) and switch to the ingoing boundary method at a certain time. Combinations with the outgoing boundary method of Bieri et al. [6] are also possible.