Black hole critical behaviour with the generalized BSSN formulation

# Black hole critical behaviour with the generalized BSSN formulation

Arman Akbarian Department of Physics and Astronomy, University of British Columbia, Vancouver BC, V6T 1Z1 Canada    Matthew W. Choptuik CIFAR Cosmology and Gravity Program
Department of Physics and Astronomy, University of British Columbia, Vancouver BC, V6T 1Z1 Canada
###### Abstract

The development of hyperbolic formulations of Einstein’s equations has revolutionized our ability to perform long-time, stable, accurate numerical simulations of strong field gravitational phenomena. However, hyperbolic methods have seen relatively little application in one area of interest, type II critical collapse, where the challenges for a numerical code are particularly severe. Using the critical collapse of a massless scalar field in spherical symmetry as a test case, we study a generalization of the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation due to Brown that is suited for use with curvilinear coordinates. We adopt standard dynamical gauge choices, including 1+log slicing and a shift that is either zero or evolved by a Gamma-driver condition. With both choices of shift we are able to evolve sufficiently close to the black hole threshold to (1) unambiguously identify the discrete self-similarity of the critical solution, (2) determine an echoing exponent consistent with previous calculations, and (3) measure a mass scaling exponent, also in accord with prior computations. Our results can be viewed as an encouraging first step towards the use of hyperbolic formulations in more generic type II scenarios, including the as yet unresolved problem of critical collapse of axisymmetric gravitational waves.

###### pacs:
04.25.dc, 04.40.-b, 04.40.Dg

## I Introduction

In this paper we investigate the application of the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation of Einstein’s equations Shibata and Nakamura (1995); Baumgarte and Shapiro (1999), as well as the dynamical coordinate choices typically associated with it, within the context of critical gravitational collapse. The BSSN formulation is a recasting of the standard 3+1 Arnowitt-Deser-Misner (ADM) Arnowitt et al. ((1962) equations that is known to be strongly hyperbolic Sarbach et al. (2002); Yoneda and Shinkai (2002) and suitable for numerical studies. It has been widely used in numerical relativity and provides a robust and stable evolution for the spacetime geometry. Most notably, various implementations of this formulation have allowed successful computation of dynamical spacetimes describing binaries of gravitationally-compact objects Campanelli et al. (2006); Baker et al. (2006); Choptuik et al. ((2015). The standard gauge choices in BSSN—namely the 1+log slicing condition Bona et al. (1997) and the Gamma-driver shift condition Alcubierre et al. (2003)—are partial differential equations (PDEs) of evolutionary type. Furthermore, the BSSN approach results in a set of so-called free evolution equations, meaning that the Hamiltonian and momentum constraints are only solved at the initial time. Thus, once initial data has been determined, one only has to solve time-dependent PDEs in order to compute the geometric variables in the BSSN scheme. In particular, during the evolution there is no need to solve any elliptic equations, which in general could arise either from the constraints or from coordinate conditions. This is advantageous since it can be quite challenging to implement efficient numerical elliptic solvers.

In addition to the BSSN approach, the numerical relativity community has adopted the generalized harmonic (GH) Pretorius (2005a) formulation of Einstein’s equations, which is also strongly hyperbolic and has performed very well in simulations of compact binaries Pretorius (2005b); Choptuik et al. ((2015). Like BSSN, the GH formulation is of evolutionary type so that all of the metric components satisfy time-dependent PDEs. It too uses dynamical coordinate choices: in this case one needs to provide a prescription for the evolution of the harmonic functions defined by .

Despite the tremendous success of these hyperbolic formulations in evolving strongly gravitating spacetimes containing black holes and neutron stars, they have not seen widespread use in another area of strong gravity physics typically studied via numerical relativity, namely critical phenomena in gravitational collapse. First reported in Choptuik (1993) and briefly reviewed below, critical phenomena emerge at the threshold of black hole formation and present significant challenges for thorough and accurate computational treatment. The original observation of critical behaviour as well as many of the subsequent studies were restricted to spherical symmetry (for a review, see Gundlach (1998); Gundlach and Martín-García (2007)) and there is a clear need to extend the work to more generic cases. In this respect the BSSN and GH formulations would appear to be attractive frameworks. However, it is not yet clear if these hyperbolic formulations, in conjunction with the standard dynamical gauge choices that have been developed, will allow the critical regime to be probed without the development of coordinate pathologies. Particularly notable in this regard is an implementation of the GH formulation that was employed by Sorkin and Choptuik Sorkin and Choptuik (2010) to study the critical collapse of a massless scalar field in spherical symmetry. Despite extensive experimentation with a variety of coordinate conditions, the code that was developed was not able to calculate near-critical spacetimes: coordinate singularities invariably formed once the critical regime was approached. A natural question that then arises is whether the BSSN formulation (including the standard dynamical gauge choices used with it) is similarly problematic or if it provides an effective framework to study critical phenomena.

Here we begin the task of addressing this question by revisiting the model of spherically symmetric massless scalar collapse. We use a generalization of the BSSN formulation due to Brown Brown (2009) that is well suited for use with curvilinear coordinates. The choice of a massless scalar field as the matter source has the great advantage that the nature of the critical solution is very well known Gundlach (1995, 1997); Hod and Piran (1997); Garfinkle and Duncan (1998); Martín-García and Gundlach (1999, 2003), making it straightforward for us to determine if and when our approach has been successful. We note that although the calculations described below are restricted to spherical symmetry our ultimate goal is to develop an evolutionary scheme—including gauge choices—that can be applied to a variety of critical phenomena studies in axial symmetry and ultimately generic cases.

We now briefly review the main concepts and features of black hole critical phenomena that are most pertinent to the work in this paper. Full details and pointers to the extensive literature on the subject may be found in review articles Gundlach (1998); Gundlach and Martín-García (2007).

Critical phenomena in gravitational collapse can be described as a phase transition, analogous to that in a thermodynamical system. Under certain assumptions, a matter source coupled to the Einstein gravitational field will evolve to one of two distinct final phases. On the one hand, weak initial data will eventually disperse to infinity leaving flat spacetime as the end state. On the other hand, sufficiently strong data will develop significant self gravitation and then collapse, resulting in a final phase which contains a black hole. Quite generically, remarkable behaviour emerges at and near the transition between these phases, and this behaviour is precisely what we mean by the critical phenomena in the system under consideration.

It transpires that there are two broad classes of critical phenomena that can be distinguished by the behaviour of the black hole mass at threshold. The class of interest here, known as type II, is characterized by infinitesimal mass at the transition. Further, the black hole mass, , satisfies a scaling law:

 MBH∼|p−p⋆|γ, (1)

where is an arbitrary parameter that controls the strength of the matter source at the initial time, is the parameter value at threshold and the mass scaling exponent, , is a constant that is independent of the choice of the initial data. Type II behaviour is also characterized by the emergence of a unique solution at threshold which is generically self-similar. In some cases, including the massless scalar field, the self-similarity is discrete. Specifically, in spherically symmetric critical collapse with discrete self-similarity (DSS), as we find

 Z⋆(ρ+Δ,τ+Δ)∼Z⋆(ρ,τ), (2)

where represents some scale-invariant component (function) of the critical solution. Here and are logarithmically rescaled values of the areal radius, , and polar time, , respectively, and is the accumulation time at which the central singularity associated with the DSS solution forms. has been normalized so that it measures proper time at the origin. As with , the echoing (rescaling) exponent, , is a universal constant for a specific matter source; i.e. it is independent of the form of the initial data.

Another feature of type II collapse, intimately related to the self-similarity of the critical solution, is that the curvature can become arbitrarily large: in the limit of infinite fine-tuning, , a naked singularity forms. Furthermore, the echoing behaviour (2) results in the development of fine structure in the solution around the center of the scaling symmetry. Observing this structure and measuring the echoing exponent associated with it requires a code that can reliably evolve solutions very close to the critical spacetime and that provides sufficient numerical resolution in the vicinity of the accumulation point .

As mentioned above, most studies of critical phenomena have assumed spherical symmetry. This is particularly so for the case of type II behaviour where the resolution demands dictated by the self-similarity of the critical solutions makes multi-dimensional work extremely computationally intensive. As far as we know, the only work in spherical symmetry to have used a purely evolutionary approach based on the BSSN or GH forms of the Einstein equations is Sorkin and Choptuik (2010) which, as we have noted, was not successful in isolating the critical solution.111However, see Pretorius and Choptuik (2000) for an investigation of type II behaviour in the collapse of a scalar field in 2+1 AdS spacetime that employs an ad hoc free evolution scheme. In axisymmetry there have been two investigations of type II collapse of massless scalar fields Choptuik et al. (2003, 2004), and several of type II collapse of pure gravitational waves (vacuum) Abrahams and Evans (1993); Alcubierre et al. (2000); Garfinkle and Duncan (2001); Rinne (2008); Sorkin (2011). Of these, only Alcubierre et al.’s Alcubierre et al. (2000) and Sorkin’s Sorkin (2011) calculations of vacuum collapse adopted hyperbolic formalisms, and only the scalar field calculations—which employed a modified ADM formulation and partially constrained evolution—were able to completely resolve the critical behaviour, including the discrete self-similarity of the critical solutions. In the fully three-space dimension (3D) context there have also been a few studies of type II collapse to date. Perhaps most notable is the recent work of Healy and Laguna Healy and Laguna (2014) which used a massless scalar field as a matter source and the BSSN formulation with standard dynamical gauge choices. The authors were able to observe the mass scaling (1) with a measured consistent with calculations in spherical symmetry. However, they were not able to conclusively see the discrete self-similarity of the critical solution; in particular they could not accurately measure the echoing exponent, . This shortcoming was attributed to a lack of computational resources rather than a breakdown of the underlying methodology, including the coordinate conditions that were adopted. Finally, there have been two attempts to probe the black hole threshold for the collapse of pure gravitational waves in 3D Lara (2006); Hilditch et al. (2013). Both employed a BSSN approach with, for the most part, standard dynamical gauge choices. In both cases problems with the gauge apparently precluded calculation near the critical point (although resolution limitations may also have been an issue) and neither the mass scaling nor the echoing exponent could be be estimated in either study.

We can thus summarize the state of the art in the use of hyperbolic formulations for the study of type II critical collapse as follows: to our knowledge there has been no implementation of a fully evolutionary scheme, based on either BSSN or GH, that has allowed for evolution sufficiently close to a precisely critical solution to allow the unambiguous identification of discrete self-similarity (or continuous self-symmetry for that matter). Again, and particularly in light of the experience of Sorkin and Choptuik (2010), the key aim of this paper is to investigate the extent to which it is possible to use a BSSN scheme to fully resolve type II solutions. A major concern here is the appropriate choice of coordinate conditions, not least since dynamical gauge choices can be prone to the development of gauge shocks and other types of coordinate singularities Alcubierre (1997, 2003). Such pathologies could, in principle, prevent a numerical solver from evolving the spacetime in or near the critical regime.

Now, as Garfinkle and Gundlach have discussed in detail Garfinkle and Gundlach (1999), an ideal coordinate system for numerical studies of type II collapse is one which adapts itself to the self-similarity: for the DSS case this means that the metric coefficients and relevant matter variables are exactly periodic in the coordinates in the fashion given by (2). Clearly, if the coordinate system is adapted, then other than at the naked singularity—which is inaccessible via finite-precision calculations—it should remain non-singular during a numerical evolution. One can then argue that ensuring that the numerical scheme has adequate resolution will be the key to successful simulation of the critical behaviour. At the same time, it is also clear that there will be coordinate systems which do not necessarily adapt but which nonetheless remain non-singular during critical collapse, at least over some range of scales, and which are therefore potentially useful for numerical calculations. We will see below that there is strong evidence that the coordinate systems we have used belong to the latter class, and weaker evidence that they do adapt to the self-similarity.

Another potential source of problems, which is not specific to hyperbolic formulations, relates to our restriction to spherical symmetry. As is well known, the singular points of curvilinear coordinate systems, in our case, can sometimes require special treatment to ensure that numerical solutions remain regular there. In critical collapse the highly dynamical nature of the solution near might naturally be expected to exacerbate problems with regularity. In the work described below we have paid special attention to the ability of our approach to both fully resolve the near-critical configuration and maintain regularity of the solution at the origin.

The remainder of this paper is organized as follows: in Sec. II we review the generalized BSSN formulation and display the equations of motion for our model system. Sec. III expands the discussion of the issue of regularity at the coordinate singularity point, describes the numerical approach we have adopted, and provides details concerning the various tests and diagnostics we have used to validate our implementation. In Sec. IV we present results computed using two distinct choices for the shift vector and provide conclusive evidence that the generalized BSSN formulation is capable of evolving in the critical regime in both cases. Sec. V contains some brief concluding remarks, and further details concerning the BSSN formalism in spherical symmetry and the scalar field equations of motion are included in App. A and App. B, respectively. We adopt units where the gravitational constant and the speed of light are both unity: .

## Ii Equations of Motion

The dynamical system we intend to study in the critical collapse regime is a real, massless scalar field, , self gravitating via Einstein’s equations,

 Gμν=8πTμν. (3)

Here, is the energy-momentum tensor associated with the minimally coupled :

 Tμν=∇μ∇νΨ−12gμν∇ηΨ∇ηΨ, (4)

and the evolution of the scalar field is given by

 ∇μ∇μΨ=0. (5)

The time-development of the geometry is then given by recasting Einstein’s equations as an evolution system based on the usual 3+1 expression for the spacetime metric:

 ds2=−α2dt2+γij(dxi+βidt)(dxj+βjdt). (6)

Here, the 3-metric components, , are viewed as the fundamental dynamical geometrical variables and the lapse function, , and shift vector, , which encode the coordinate freedom of general relativity, must in general be prescribed independently of the equations of motion.

### ii.1 Generalized BSSN

We now summarize the BSSN formulation of Einstein’s equations and describe how it can be adapted to curvilinear coordinates. Readers interested in additional details are directed to Baumgarte and Shapiro ((2010) for a more pedagogical discussion.

In the standard ADM formulation Arnowitt et al. ((1962); York ((1979), the dynamical Einstein equations are rewritten as evolution equations for the 3-metric and the extrinsic curvature . The first difference between the BSSN formulation and the ADM decomposition is the conformal re-scaling of the ADM dynamical variables:

 γij=e4ϕ~γij, (7)
 Kij=e4ϕ~Aij+13γijK, (8)

where is the conformal factor, is the conformal metric, is the conformally rescaled trace-free part of the extrinsic curvature and is the trace of the extrinsic curvature. Here by fixing the trace of , and the determinant of the conformal metric, the set of primary ADM dynamical variables transforms to the new set:

 {γij,Kij}→{ϕ,~γij,K,~Aij}, (9)

in the BSSN formulation.

In the original BSSN approach, the conformal metric is taken to have determinant . However this choice is only suitable when we adopt coordinates in which the determinant of the flat-space metric reduces to unity. This is the case, of course, for Cartesian coordinates but is not so for general curvilinear systems. For instance, the flat 3-metric in spherical coordinates:

 ds2=dr2+r2dθ2+r2sin2θdϕ2, (10)

has determinant . Recently, Brown Brown (2009) has resolved this issue by introducing a covariant version of the BSSN equations—the so-called generalized BSSN formulation, which we will hereafter refer to as G-BSSN—in which the primary dynamical variables are tensors so that the formulation can be adapted to non-Cartesian coordinate systems. In G-BSSN we no longer assume that the conformal 3-metric has determinant one. Rather, becomes a true scalar and for its dynamics to be determined a prescription for the time evolution of the determinant of must be given. In the following this will be done by requiring that the determinant be constant in time.

Another main difference between the ADM decomposition and BSSN is that the mixed spatial derivative terms occurring in the 3-Ricci tensor are eliminated through the definition of a new quantity, :

 ~Γk≡~γij~Γkij, (11)

which becomes an additional, independent dynamical variable. Note that is not a vector as it is coordinate dependent. To extend this redefinition so that it is well suited for all coordinate choices, in G-BSSN we define

 ~Λk≡~γij(~Γkij−˚Γkij)=~Γk−˚Γkij~γij, (12)

where denotes the Christoffel symbols associated with the flat metric. This definition makes this so-called conformal connection, , a true vector and it becomes a primary dynamical variable in G-BSSN.

We now summarize the G-BSSN equations, referring the reader to Alcubierre and Mendez (2011) for more details, including a full derivation. We begin by defining , the time derivative operator acting normally to the slices:

 ∂⊥≡∂t−L→β, (13)

where denotes the Lie derivative along . We then have

 ∂⊥ϕ=−16αK+σ16~Dkβk, (14)
 ∂⊥~γij=−2α~Aij−σ23~Aij~Dkβk, (15)
 ∂⊥K=−γijDjDiα+α(~Aij~Aij+13K2)+4π(ρ+S), (16)
 ∂⊥~Aij = e−4ϕ[−DiDjα+α(Rij−8πSij)]TF (17) + α(K~Aij−2~Ail~Alj)−σ23~Aij~Dkβk,
 ∂⊥~Γi = −2~Aij∂jα+~γlj∂j∂lβi (18) + 2α(~Γijk~Akj−23~γij∂jK+6~Aij∂jϕ−8π~γijSj) + σ3[2~Γi~Dkβk+~γli∂l(~Dkβk)].

Here, a superscript TF denotes the trace-free part (with respect to the 3-metric ) of a tensor, and is the covariant derivative associated with the conformal metric . Additionally, the quantity is an adjustable parameter that is discussed below and typically is either 0 or 1. Note that all the Lie derivatives in the G-BSSN equations operate on true tensors and vectors of weight 0. For instance,

 L→β~Aij=βk∂k~Aij+~Aik∂jβk+~Akj∂iβk. (19)

Furthermore, in G-BSSN, rather than evolving (18), the redefined conformal connection, , is evolved via

 ∂t~Λk=∂t~Γk−˚Γkij∂t~γij, (20)

where the time derivative is eliminated using (15). In equation (17), denotes the 3-Ricci tensor associated with and can be written as the sum

 Rij=Rϕij+~Rij, (21)

where is given by

 Rϕij=−2~Di~Djϕ−2~γij~Dk~Dkϕ+4~Diϕ~Djϕ−4~γij~Dkϕ~Dkϕ, (22)

and is the 3-Ricci tensor associated with the conformal metric:

 ~Rij= − 12~γlm∂m∂l~γij+~γk(i∂j)~Γk+~Γk~Γ(ij)k (23) + ~γlm(2~Γkl(i~Γj)km+~Γkim~Γklj).

The matter fields , , and are defined by

 ρ=nμnνTμν, (24)
 S=γijSij, (25)
 Si=−γijnμTμj, (26)
 Sij=γiμγiνTμν, (27)

where is the unit normal vector to the slices.

As mentioned previously, we need to prescribe dynamics for the determinant of to have a complete set of equations of motion for the G-BSSN dynamical variables. One approach is to fix the determinant to its initial value by demanding that

 ∂t~γ=0. (28)

This is the so-called Lagrangian option and is associated with the choice in the equations. Another option is to define the determinant to be constant along the normal direction to the time slices, which can be implemented by requiring . This is usually referred to as the Lorentzian option, and is associated with the choice . Here we choose (28), i. e. .

Note that in the G-BSSN equations the divergence of the shift vector,

 ~Dkβk=1√~γ∂k(√~γβk), (29)

no longer reduces to since the determinant of the conformal metric is not necessarily 1, but by virtue of the choice (28) is equal to that of the initial background flat metric in the chosen curvilinear coordinates.

As usual, when setting initial data for any given evolution of the coupled Einstein-matter equations we must solve the Hamiltonian and momentum constraints. In terms of the G-BSSN variables these are

 H ≡ ~γij~Di~Djeϕ−eϕ8~R+e5ϕ8~Aij~Aij (30) − e5ϕ12K2+2πe5ϕρ=0,
 Mi≡~Dj(e6ϕ~Aji)−23e6ϕ~DiK−8πe6ϕSi=0. (31)

### ii.2 G-BSSN in spherical symmetry and gauge choices

In spherical symmetry a generic form of the conformal metric is given by

 ~γij=⎛⎜ ⎜⎝~γrr(t,r)000r2~γθθ(t,r)000r2~γθθ(t,r)sin2θ⎞⎟ ⎟⎠. (32)

Similarly, a suitable ansatz for the traceless extrinsic curvature is

 ~Aij=⎛⎜ ⎜⎝~Arr(t,r)000r2~Aθθ(t,r)000r2~Aθθ(t,r)sin2θ⎞⎟ ⎟⎠. (33)

The shift vector and have only radial components:

 βi=[β(t,r),0,0], (34)
 ~Λi=[~Λ(t,r),0,0]. (35)

Given (32-35), the G-BSSN equations become a set of first order evolution equations for the 7 primary variables

 {ϕ(t,r),~γrr(t,r),~γθθ(t,r),K(t,r), ~Arr(t,r),~Aθθ(t,r),~Λ(t,r)}.

These are coupled to the evolution equation (5) for the scalar field and constrained by the initial conditions (3031). The explicit expressions for the full set of equations of motion are given in App. A.

To fix the time slicing we implement a non-advective222 The terminology non-advective derives from the absence of an “advective” term, , on the left hand side of equations (36,38). we note that we also experimented with the advective versions of the equations. the results were very similar to those for the non-advective case; in particular, near-critical solutions exhibiting echoing and scaling could also be obtained. version of the 1+log slicing condition:333 The reader can easily check that in the case of zero shift, the lapse choice given by (36) combined with (14) implies . in Cartesian coordinates , so this last equation gives , where the function is time independent. the choice then yields an algebraic expression for the lapse, , which is the origin of the terminology “1+log slicing”.

 ∂tα=−2αk. (36)

for the spatial coordinates we either choose a zero shift:

 βi=0, (37)

or use what we will term the gamma-driver condition:

 ∂tβi=μ~Λi−ηβi. (38)

Here, and are adjustable parameters which we set to and , where is the total mass of the system measured at infinity (see Sec. III.4.1). We emphasize that (38) is not the usual Gamma-driver equation used in the standard BSSN approach:

 ∂tβi=μ~Γi−ηβi, (39)

but since it is a natural extension of the above to the G-BSSN case we have opted to use the same nomenclature. In the rest of this paper, we frequently refer to the shift vector evolved via (38) as . Explicitly, in spherical symmetry is defined by

 ∂tβG(t,r)=μ~Λ(t,r)−ηβG(t,r). (40)

## Iii Numerics

We use a second order finite differencing method to discretize equations (14-17) and (20). Further, the equations of motion are transformed to a compactified radial coordinate that we denote by and which is defined in terms of the original coordinate by

 r=e~r−eδ+R∞R∞−~r−R∞R∞−δ, (41)

where and are parameters with typical values and . It is straightforward to verify the following: (1) the radial domain maps to the computational domain , (2) the derivative decreases toward the origin (), so that a uniform grid on is a non-uniform grid on with approximately times more resolution close to the origin relative to the outer portion of the solution domain, (), where the support of the scalar field is initially concentrated, (3) the parameter can be used to adjust the resolution near the origin; specifically, decreasing increases the resolution near . For notational simplicity, however, in the following we omit the explicit dependence of the fields on and denote the spacetime dependence of any dynamical variable as previously: .

We use a finite difference grid that is uniform in and analytically transform all -derivative terms in the equations of motion to their -coordinate counterparts prior to finite-differencing.

We also developed a Maple-based toolkit Akbarian (2014) that automates the process of discretizing an arbitrary derivative expression. This toolkit handles boundary conditions and generates a point-wise Newton-Gauss-Seidel solver in the form of Fortran routines for a given set of time dependent or elliptic PDEs . The calculations in this paper were all carried out using this infrastructure.

### iii.1 Initialization

The matter content is set by initializing the scalar field to a localized Gaussian shell:

 Ψ(0,r)=pexp(−(r−r0)2σ2r), (42)

where , and are parameters. Note that here is the non-compactified radial coordinate which is related to the compactified coordinate via (41). A typical initial profile for the scalar field in our calculations has , , and of order . We use the overall amplitude factor as the tuning parameter to find critical solutions. We initialize the conformal metric (32) to the flat metric in spherical symmetry,

 ~γrr(0,r)=~γθθ(0,r)=1, (43)

and initialize the lapse function to unity,

 α(0,r)=1. (44)

We also demand that the initial data be time-symmetric,

 ~Arr(0,r)=~Aθθ(0,r)=K(0,r)=0, (45)
 β(0,r)=~Λ(0,r)=0, (46)
 ∂tΨ(t,r)|t=0=0, (47)

which means that the momentum constraint (31) is trivially satisfied. This leaves the Hamiltonian constraint (30) which is solved as a two-point boundary value problem for the conformal factor at the initial time,

 ψ(r)≡eϕ(0,r). (48)

The outer boundary condition for ,

 ψ(r)|r=∞=1, (49)

follows from asymptotic flatness, while at we have

 ∂rψ(r)|r=0=0 (50)

since must be an even function in for regularity at the origin.

### iii.2 Boundary Conditions

Due to the fact that the metric has to be conformally flat at the origin we have

 ~γrr(t,0)=~γθθ(t,0). (51)

Further, since we are using the Lagrangian choice, , the determinant of must at all times be equal to its value at the initial time, so

 ~γrr~γ2θθ=1. (52)

From these two results we have

 ~γrr(t,0)=~γθθ(t,0)=1. (53)

Using (53) and (15) it is then easy to see that we must also have

 ~Arr(t,0)=~Aθθ(t,0)=0. (54)

As is usual when working in spherical coordinates, many of the boundary conditions that must be applied at follow from the demand that the solution be regular there. Essentially, the various dynamical variables must have either even or odd “parity” with respect to expansion in as . Variables with even parity, typically scalars or diagonal components of rank-2 tensors, must have vanishing radial derivative at , while odd parity functions, typically radial components of vectors, will themselves vanish at the origin.

Applying these considerations to our set of unknowns we find

 ∂r~γrr(t,r)|r=0=∂r~γθθ(t,r)|r=0=0, (55)
 ∂r~Arr(t,r)|r=0=∂r~Aθθ(t,r)|r=0=0, (56)
 β(t,0)=~Λ(t,0)=0, (57)
 ∂rK(t,r)|r=0=∂rϕ(t,r)|r=0=∂rΨ(t,r)|r=0=0. (58)

We use equations (53,54,57) to fix the values of the functions at the origin and a forward finite-differencing of (58) to update , and at . Further, we apply a forward finite-differencing of (55,56) to update the values of the function at the grid point next to the origin. The 1+log condition (36) can be used directly at . Again, we emphasize that all of the -derivative terms of the boundary conditions described above are analytically transformed to the numerical coordinate, , before the equations are finite-differenced.

Since we are using compactified coordinates, all the variables are set to their flat spacetime values at the outer boundary :

 ~γrr=~γθθ=eϕ=α=1at:(t,∞), (59)
 ~Arr=~Aθθ=K=~Λ=β=Ψ=0at:(t,∞). (60)

Here, we emphasize that spatial infinity, , corresponds to the finite compactified (computational) coordinate point .

### iii.3 Evolution scheme and regularity

We implemented a fully implicit, Crank-Nicolson Crank and Nicolson (1996) finite differencing scheme to evolve the system of G-BSSN equations. The precise form of the continuum equations used is given in App. A and all derivatives, both temporal and spatial, were approximated using second-order-accurate finite difference expressions.

During an evolution the correct limiting behaviour of the spatial metric components must be maintained near to ensure a regular solution. For example, the limiting values of the conformal metric components and are given by

 ~γrr(t,r)=1+O(r2), (61)
 ~γθθ(t,r)=1+O(r2). (62)

If the discrete approximations of the metric functions do not satisfy these conditions, then irregularity will manifest itself in the divergence of various expressions such as the Ricci tensor component (95)

 Rrr=2~γrr−~γθθr2~γθθ+⋯, (63)

which should converge to a finite value at the origin if conditions (61,62) hold.

One approach to resolve potential regularity issues is to regularize the equations Ruiz et al. (2008); Alcubierre and Mendez (2011); Sorkin and Choptuik (2010), by redefining the primary evolution variables, so that the equations become manifestly regular at the origin. Another approach is to use implicit or partially implicit methods Cordero-Carrión and Cerda-Duran (2012). As recently shown by Montero and Cordeo-Carrion Montero and Cordero-Carrión (2012), such schemes can yield stable evolution without need for explicit regularization. Baumgarte et al. Baumgarte et al. (2013) also adopted a similar approach—using a partially implicit scheme without regularization—in an implementation of the G-BSSN formulation in spherical polar coordinates.

As mentioned, our implementation is fully implicit and we have also found that our generalized BSSN equations can be evolved without any need for regularization at the origin, even in strong gravity scenarios where the spacetime metric has significant deviations from flatness near the origin.

That said, we also experimented with other techniques aimed at improving regularity. For example, using the constraint equation (52) and the fact that is trace-free,

 ~Arr~γrr+2~Aθθ~γθθ=0, (64)

we can compute and in terms of and , respectively, rather than evolving them. However, when we did this we found no significant improvement in regularity relative to the original scheme.

Finally, to ensure our solutions remain smooth on the scale of the mesh we use fourth order Kreiss-Oliger dissipation Kreiss and Oliger ((1973) in the numerical solution updates.

### iii.4 Tests

This section documents various tests we have made to validate the correctness of our numerical solver as well as the consistency of the finite-differencing method used to evolve the system of G-BSSN equations. We use a variety of diagnostic tools, including monitoring of the constraint equations, convergence tests of the primary dynamical variables, and a direct computation to check if the metric and matter fields calculated via the G-BSSN formulation satisfy the covariant form of Einstein’s equations. All of the calculations were performed using the 1+log slicing condition (36) and either or where satisfies the Gamma-driver condition (40).

#### iii.4.1 Constraints and conserved quantities

We monitor the evolution of the constraint equations (30,31) during a strongly-gravitating evolution where the nonlinearities of the equations are most pronounced. As demonstrated in Fig. 1 (a), at resolutions typical of those used in our study, the Hamiltonian constraint is well preserved during such an evolution and, importantly, the deviations from conservation converge to zero at second order in the mesh spacing as expected.

The total mass-content of the spacetime seen at spatial infinity (the ADM mass) is a conserved quantity. Here, using the G-BSSN variables the Misner-Sharp mass function is given by

 M(r)=r~γ1/2θθe2ϕ2[1−~γθθ~γrr(1+r∂r~γθθ2~γθθ+2r∂reϕeϕ)2]. (65)

The total mass, , can be evaluated at the outer boundary,

 MADM≡M(r=∞). (66)

The deviation of the total mass from its time average is plotted in Fig. 1(b); as the resolution of the numerical grid increases the variations converge to zero in a second order fashion.

#### iii.4.2 Convergence test

As mentioned in Sec. III.1 and Sec. III.3, we implemented our code using second order finite differencing of all spatial and temporal derivatives. Denoting any continuum solution component by , where is the spatial coordinate, and a discrete approximation to it computed at finite difference resolution, , by , to leading order in we expect

 qh(t,X)=q(t,X)+h2e2[q](t,X)+…. (67)

Fixing initial data, we perform a sequence of calculations with resolutions , and and then compute a convergence factor, , defined by

 C(t;q)=||qh(t,X)−qh/2(t,X)||2||qh/2(t,X)−qh/4(t,X)||2, (68)

where is the norm, i.e. the root mean square (RMS) value. It is straightforward to argue from (67) that, for sufficiently small , should approach 4 if the solution is converging at second order. The values of the convergence factor for a selection of dynamical variables are plotted for a strong-data evolution in Fig. 1(c) and provide clear evidence that the solution is second-order convergent throughout the time evolution.

#### iii.4.3 Direct validation via Einstein’s equations

A direct method to test the fidelity of our numerical solver involves the evaluation of a residual based on the covariant form of Einstein’s equations. We start with a reconstruction of the four-dimensional metric in spherical symmetry,

 ds2=(−α2+β2a2)dt2+2a2βdtdr+a2dr2+r2b2dΩ2, (69)

using the primary G-BSSN variables, and . In particular, and are simply given by

 a(t,r)=e4ϕ(t,r)~γrr(t,r), (70)
 b(t,r)=e4ϕ(t,r)~γθθ(t,r). (71)

We then check to see if the metric (69) satisfies the covariant Einstein equations (3) to the expected level of truncation error. Specifically, defining the residual

 Eμν≡Gμν−8πTμν, (72)

and replacing all derivatives in with second order finite differences, we expect to converge to zero as as .444Although it is not crucial for the usefulness of this test, we discretize the using a difference scheme that is distinct from the one used in the main code. Precisely this behaviour is shown in Fig. 1(d). This is a particularly robust test of our implementation since the non-trivial components of the covariant Einstein equations are quite complicated and, superficially at least, algebraically independent of the BSSN equations. For instance, the component of the residual (72) is given by

 Etr = 2βrα2(∂raa−2∂rbb+∂rαα) (73) + 2βα2(−∂2rbb+∂ra∂rbab−(∂rΨ)22+∂rα∂rbαb) + 2α2(∂t∂rbb+∂rΨ∂tΨ2−∂ta∂rbab−∂tb∂rααb) + 2rα2(−∂taa+∂tbb)

and depends on all of the dynamical variables of the system. The observed convergence of the residual is only plausible if (1) our G-BSSN equations (14-18) have been correctly derived from the covariant Einstein equations, (2) we have discretized the geometric and matter equations properly, and (3) we have solved the full set of discretized equations correctly.

### iii.5 Finding black hole threshold solutions

The strength of the initial data can be set by adjusting the amplitude of the scalar field, , in (42). For weak enough initial data (small enough ), the matter shell will reach the origin and then disperse, with the final state being a flat spacetime geometry. Sufficiently strong initial data, on the other hand (large enough ), results in a matter concentration in the vicinity of the origin which is sufficiently self-gravitating that a black hole forms. Using a binary search, we can find the threshold initial data, defined by , for which results in dispersal while yields black hole formation. At any stage of the calculation, the binary search is defined by two “bracketing” values, and , such that evolutions with and result in dispersal and black hole formation, respectively. It is convenient to define the amount of parameter tuning that has occurred by the dimensionless quantity

 δp≡ph−plpl. (74)

The dispersal case can be detected easily as the scalar field leaves the vicinity of the origin and the geometry approaches flat spacetime. To detect black hole formation, we use an apparent horizon finder to locate a surface on which the divergence of the outgoing null rays vanishes. We first define the divergence function

 Θ=qμν∇μkν, (75)

where is the induced metric on the constant surface. In spherical symmetry with metric (69) we have

 qμν=diag(0,0,r2b2,r2b2sin2θ), (76)

where is the null outgoing vector given by

 kμ=1√2[aβ−α,a,0,0]. (77)

Therefore, (75) becomes

 Θ=√2rb(rα∂t(b)+(1a−βα)∂r(rb)). (78)

The formation of an apparent horizon555Technically a marginally trapped surface—the apparent horizon being the outermost of these. is signaled by the value of the function crossing zero at some radius and, modulo cosmic censorship, implies that the spacetime contains a black hole. We note that since the focus of our work was on the critical (threshold) solution we made no effort to continue evolutions beyond the detection of trapped surfaces.

## Iv Results

In this section we describe results from two sets of numerical experiments to study the efficacy of the G-BSSN formulation in the context of critical collapse. Again, our calculations use the standard 1+log slicing condition for the lapse, and a shift which is either zero or determined from the Gamma-driver condition.

### iv.1 Zero shift

We first perform a collection of numerical experiments where the shift vector is set to zero. As described in Sec. III.5, in principle we can find the black hole threshold solution using a binary search algorithm which at any stage is defined by two values and , with , and where corresponds to dispersal (weak data) while corresponds to black hole formation (strong data).

As discussed in the introduction, the massless scalar collapse model has a very well-known critical solution, and we summarize the features most relevant to our study here. The threshold configuration is discretely self-similar with an echoing exponent measured from the first calculations to be Choptuik (1993). Following the original studies, Gundlach Gundlach (1995) showed that the construction of the precisely discretely self-similar spacetime could be posed as an eigenvalue problem, the solution of which led to the more accurate value . This estimate was subsequently improved by Martin-Garcia and Gundlach to  Martín-García and Gundlach (2003).

The original calculations determined a value for the mass-scaling exponent Choptuik (1993); further work based on perturbation theory gave  Gundlach (1997); Martín-García and Gundlach (1999). Here it is important to note that, as pointed out independently by Gundlach Gundlach (1997) and Hod and Piran Hod and Piran (1997), the simple power law scaling (1) gets modified for discretely self-similar critical solutions to

 lnM=γln|p−p⋆|+c+f(γln|p−p⋆|+c), (79)

where is a universal function with period and is a constant depending on the initial data. This results in the superposition of a periodic “wiggle” in the otherwise linear scaling of as a function of .

Finally, Garfinkle and Duncan Garfinkle and Duncan (1998) pointed out that near-critical scaling is seen in physical quantities other than the mass and, dependent on the quantity, in the subcritical as well as supercritical regime. In particular they argued that in subcritical evolutions the maximum central value, , of the four-dimensional Ricci scalar, , defined by

 Rmax≡maxtR(0,t), (80)

should satisfy the scaling

 Rmax∼|p−p⋆|−2γ, (81)

where the factor in the scaling exponent can be deduced from the fact that the curvature has units of . For the discretely self-similar case this scaling law is also modulated by a wiggle with period , which for the massless scalar field is about 4.61.

Using initial data given by (42) we tune so that it is close to the critical value: typically this involves reducing the value of defined by (74) so that it is about , which is a few orders of magnitude larger than machine precision. Our implementation includes code that actively monitors the dynamical variables for any indications of coordinate singularities or other pathologies which could cause the numerical solver to fail. Provided that such pathologies do not develop, we expect to observe features characteristic of critical collapse—discrete self-similarity and mass scaling in particular—to emerge as .

One way the discrete self-similarity of the critical solution is manifested is as a sequence of “echoes”—oscillations of the scalar field near the origin such that after each oscillation the profile of the scalar field is repeated but on a scale smaller than that of the preceding echo (see Eq. (2)). The oscillations are similarly periodic in the logarithmic time scale , where is the proper time measured at the origin,

 T(t)≡∫t0α(^t,0)d^t, (82)

and is the accumulation time at which the singularity forms (always at ). Furthermore, viewed at the origin, the oscillations of the scalar field occur at a fixed amplitude of about (with our units and conventions for the Einstein’s equations). As shown in Fig. 2, when we tune the initial data to the critical value, the central value of scalar field exhibits oscillatory behaviour and the amplitude is close to the expected value. The anticipated periodicity in logarithmic time is also apparent with a measured , in agreement with previous results. We thus have strong evidence that the evolution has indeed approached the critical regime and that the measured oscillations are true echos rather than numerical artifacts.

Evidence that our code correctly captures the expected critical scaling behaviour (81) of is presented in Fig. 3. We find , consistent with previous calculations. We note that we can measure scaling from our computations up to (or ). However, in Fig. 3 we have excluded the last few values closest to the critical point from both the plot and the linear fit: specifically, we truncate the fit at . The rationale for this is that we use the largest subcritical value of as an approximation to the critical value rather than, for example, implementing a multi-parameter fit that includes as one of the parameters. Our estimate of thus has an intrinsic error of and by fitting to data with we render the error in the estimate essentially irrelevant. We note that consistent with the early observations of the robustness of mass scaling in the model Choptuik (1993), measuring the exponent can be achieved by moderate tuning, in this case , (i.e. ). However, to be able to observe the echoing exponent (the oscillations around the fitting line, for example) we need to tune much closer to the critical value.

The echoing behaviour of the critical solution is also reflected in the geometry of spacetime and the matter variables other than the scalar field. Fig. 4 shows the radial profile of the G-BSSN variable at an instant close to the accumulation time . As seen in this plot, fine structure develops in the function in the near-critical regime. Observe that from the definition (7) and the choice (28), the scalar is the ratio of the determinant of the 3-metric, , to the determinant of the flat metric, :

 ϕ=112ln(γ/˚γ). (83)

The radial matter density, , is a convenient diagnostic quantity for viewing near-critical evolution. Snapshots of this function from a typical marginally subcritical calculation are shown in Fig. 5: the echoing behavior is clearly evident in the sequence. The number of echos is dependent on the degree to which the solution has been tuned to criticality. In this case, where , we expect and see about 4 echos (last frame of the figure). Here we note that each of the echos in corresponds to half of one of the scalar field oscillations shown in Fig. 2 (where the inset shows about full cycles).

Fig. 6(a) plots the central matter density for a marginally supercritical calculation. In accord with the self-similar nature of the near-critical solution, the central density grows exponentially with time. Fig. 6(b) is a snapshot of the extrinsic curvature at the critical time while Fig. 6(c) shows the dynamics of the central value of the lapse function and compares it with from the calculations performed with described in the next section. Fig. 6(d) displays the profile of the lapse at the critical time.

We note that we have not fully resolved the critical behaviour in the sense of having tuned to the limit of machine precision, , which would capture roughly 2 additional echos in the threshold solution (one full echo in the scalar field). In principle, by setting sufficiently large we could almost certainly do so since there are no indications that our method would break down at higher resolution and closer to criticality. However, we estimate that the required compute time for a complete critical search would increase from weeks to several months and we have thus not done this. Ultimately, a more effective approach to enhancing the resolution would be to incorporate a technique such as adaptive mesh refinement into our solver.

The results displayed in Figs. 2–6 provide strong evidence that the coordinate system consisting of 1+log lapse and zero shift remains non-singular in the critical regime, at least for the range of scales probed for . Additionally, the approximate periodicity in that can be seen, for example, in (Fig. 4) and (Fig. 5) suggests that the coordinates may be adapting to the self-similarity. Whether or not this is actually the case is a matter requiring further study.

### iv.2 Gamma-driver shift

We now briefly report on experiments similar to those of the previous section but where the shift was evolved with the Gamma-driver condition (38). A principal observation is that this gauge also facilitates near-critical evolutions with results similar to the choice. In particular, we are again able to observe all of the characteristics of the black hole threshold solution.

The gauge condition (38) acts as a damping factor for the conformal connection, , and we would therefore expect to observe a significant change in the profile of at threshold relative to the zero-shift case. This expectation is borne out by the comparison illustrated in Figs. 7 (a) and (b). When , diverges as close to the origin while it appears to have finite amplitude for . We find that the shift develops very sharp oscillations near the origin; some typical behaviour can be seen in the plot of shown in  Fig. 7(c). We believe that these oscillations are genuine and our expectation is that will diverge in the precise critical limit. Further, we observe that the oscillations can create numerical artifacts and generally require higher resolution relative to the case, as well as dissipation, to be controlled. Indeed, when using the Gamma-driver condition we find that Kreiss/Oliger dissipation is crucial to suppress unresolved high frequency oscillations close to the origin. Fig. 7 (d) shows the growth in the norm of the extrinsic curvature during a subcritical evolution. The norm of does not exhibit any significant difference for the two choices of the shift.

As was the case for the calculations, the results shown in Figs. 6 and 7 strongly suggest that the combination of 1+log slicing and Gamma driver shift provides a coordinate system which is adequate for computing the near-critical solution. In addition, the approximate periodicity seen in Figs. 6(b), 6(c), 7(a) and 7(b) suggest that this gauge may also be adapting to the self-symmetry.

## V Conclusion

We have described a numerical code that implements a generalized BSSN formulation adapted to spherical symmetry. Using standard dynamical coordinate choices, including 1+log slicing and a shift which either vanished or satisfied a Gamma-driver condition, we focused specifically on the applicability of the formulation and the gauge choices to studies of type II critical phenomena. As a test of the approach we revisited the model of massless scalar collapse, where the properties of the critical solution are very well known from previous work. For both choices of the shift, we found that our code was able to generate evolutions that were very close to criticality so that, in particular, we could observe the expected discrete self-similarity of the critical solution. To our knowledge, this is the first fully evolutionary implementation of a hyperbolic formulation of Einstein’s equations that has been able to unequivocally resolve discrete self-similarity in type II collapse. Furthermore, measured properties from near-critical solutions, including the mass-scaling and echoing exponents, are in agreement with previous work. Our results strongly suggest that the G-BSSN formulation, in conjunction with standard dynamical coordinate conditions, is capable of evolving the spacetime near criticality without the development of coordinate pathologies. There is also some evidence that both gauges adapt to the self-similarity, but we have not yet studied this issue in any detail.

We found that certain of the primary G-BSSN variables diverge as the critical solution is approached: this is only to be expected since the precisely critical solution contains a naked singularity. Dealing with such solution features in a stable and accurate manner presents a challenge for any code and in our case we found that a combination of a non-uniform grid and Kreiss/Oliger dissipation was crucial. Our use of a time-implicit evolution scheme may have also been important although we did not experiment with that aspect of our implementation. However, we suspect that the implicit time-stepping helped maintain regularity of the solutions near , as other researchers have found.

Given the success of the G-BSSN approach, it is natural to consider its generalization and application to settings with less symmetry, but where curvilinear coordinates are still adopted. In particular, one axisymmetric problem that has yet to be resolved is the collapse of pure gravity waves. This scenario arguably provides the most fundamental critical phenomena in gravity as the behaviour must be intrinsic to the Einstein equations, rather than being dependent on some matter source. Critical collapse of gravitational waves–with mass scaling and echoing—was observed by Abrahams and Evans over 20 years ago Abrahams and Evans (1993). However, their original results have proven very difficult to reproduce (or refute) Alcubierre et al. (2000); Garfinkle and Duncan (2001); Rinne (2008); Sorkin (2011); Hilditch et al. (2013). We refer the reader to the recent paper by Hilditch et al. Hilditch et al. (2013) for detailed discussions concerning some apparent inconsistencies among the follow-up studies, as well as the challenges and complications involved in evolving various types of nonlinear gravitational waves. We are currently extending the methodology described above to the axisymmetric case with plans to use the resulting code to study vacuum collapse. Results from this undertaking will be reported in a future paper.

###### Acknowledgements.
This research was supported by NSERC, CIFAR and by a 4FY scholarship to AA from UBC. The authors thank William G. Unruh for insightful comments and a thorough reading of a draft of this paper.

## Appendix A BSSN in Spherical Symmetry

In this appendix, we provide the explicit expressions of the G-BSSN evolution equations in spherical symmetry.

The evolution equations (14-15) for and the components of the conformal metric simplify to

 ∂tϕ=16αK+β∂rϕ+σ16B, (84)
 ∂t~γrr=−2α~Arr+β∂r~γrr+2~γ