# An axisymmetric generalized harmonic evolution code

###### Abstract

We describe the first axisymmetric numerical code based on the generalized harmonic formulation of the Einstein equations, which is regular at the axis. We test the code by investigating gravitational collapse of distributions of complex scalar field in a Kaluza-Klein spacetime. One of the key issues of the harmonic formulation is the choice of the gauge source functions, and we conclude that a damped wave gauge is remarkably robust in this case. Our preliminary study indicates that evolution of regular initial data leads to formation both of black holes with spherical and cylindrical horizon topologies. Intriguingly, we find evidence that near threshold for black hole formation the number of outcomes proliferates. Specifically, the collapsing matter splits into individual pulses, two of which travel in the opposite directions along the compact dimension and one which is ejected radially from the axis. Depending on the initial conditions, a curvature singularity develops inside the pulses.

^{†}

^{†}preprint: AEI-2009-108

## I Introduction

In general, a detailed investigation of fully nonlinear gravitational dynamics is impossible by other than numerical means. Luckily, the numerical methods have recently reached the level of maturity that finally allows addressing many long-standing puzzles. Perhaps the most remarkable is the progress achieved in solving a general-relativistic two-body problem—the coalescence of black holes FP2 (); FP1 (); FP3 (); RIT (); nasa (); AEI (); CaltechCornell (). Driven by the gravitational wave detection prospects, the problem of the collision of two black holes or neutron stars continues to be the central front where an overwhelming majority of numerical relativity research is done. Fortunately, the computational methods employed there are portable and—as demonstrated below—can readily be applied on other problems of interest.

The success of the numerical simulations was backed up by a
parallel development of the software and the hardware, which provided
the necessary computational resources. The rapid hardware evolution combined with the
persisting regularity problems in axial symmetry eventually led to
a direct transition from highly symmetrical spherical configurations
to fully general 3D situations without any symmetries at all,
essentially bypassing the intermediate axisymmetric case.
However, here we argue that important theoretical and practical
reasons exist to explore axisymmetry better, and
we describe a new regular numerical code that, we believe, will be capable of
achieving this. ^{1}^{1}1See Garfinkle_axi (); Choptuik_axi (); Rinne_axi_1 (); Rinne_axi_2 ()
for alternative approaches.

Before describing any concrete setup we would point out one important possible
use of an axisymmetric code, specifically, that
it can be regarded as an efficient “calibration tool” for more
general 3D codes. ^{2}^{2}2and as a probe of reliability of the
cartoon methods cartoon_method (); FP1 () used to effectively
simulate axisymmetric spacetimes in 3D. Indeed, we expect that an
intrinsically axisymmetric code applied to, say, the head-on collision of
two black holes would be capable of following the evolution and the resulting gravitational radiation
more accurately compared to Cartesian 3D numerical implementations, both because of explicit use of
the symmetry and since higher numerical resolution can be employed
for given hardware resources.

Several interesting open problems arise in axisymmetric gravitational collapse situations. In particular, it remains unclear whether or not the weak cosmic censorship is violated in collapse of prolate Brill waves BrillWaves (); Rinne_axi_1 (). An independent observation of universality in critical collapse of gravitational waves AbrahamsEvans () is pending, as well as further investigation of the non-spherical unstable mode that apparently shows up at threshold for black hole formation in axisymmetric collapse of a scalar field crit_collapse_axi (). A basic problem of mathematical relativity concerning the stability of black holes with respect to nonlinear axisymmetric perturbations, can be equivalently addressed in a collapse situation by computing how fast a newly formed black hole radiates away higher multipole moments.

In addition, we shall mention other, rarely cited in the context of numerical relativity, axisymmetric systems which are of great interest in theoretical work on higher-dimensional gravity. One of the most basic motivations for studying higher-dimensional spacetimes relies on the observation that the Einstein equations, describing classical General Relativity (GR), are independent of the spacetime dimension. Nevertheless, certain properties of the solutions to these equations vary dramatically with the dimension. A striking example is that axisymmetric black holes in dimensions greater than four do not necessarily have spherical horizons, but also admit horizons of toroidal topology BlackRing (). Moreover, and in sharp contrast to their four-dimensional counterparts, higher dimensional black holes do not respect the Kerr limit MP (), and they are unstable in certain range of parameters GL (). One of the fundamental unresolved puzzles HM (); Choptuik_BS () is whether or not the instability leads to fragmentation of the horizon and exposition of the inner singularity, hence violating the cosmic censorship hypothesis.

Since it is unlikely that any of these these problems can be addressed analytically in a systematic manner, we turn to computational methods. However, solving the Einstein equations numerically is notoriously difficult and depends crucially on the way these equations are formulated and evolved. In this paper we focus on the generalized harmonic (GH) formulation Friedrich85 (); Friedrich96 (); Garfinkle2001 () that has recently gained popularity because of its great success in the simulations of black hole binaries FP2 (); FP3 (); CaltechCornell (); FP1 (); Lindblom_etal ().

In a nutshell, the GH approach is a way to write the field equations such that the resulting system is manifestly hyperbolic, taking the form of a set of quasilinear wave equations for the metric components. As the name suggests, the GH method generalizes the harmonic approach, which achieves strong hyperbolicity by choosing the harmonic spacetime coordinates. In the GH approach much of the coordinate freedom is regained through the introduction of certain gauge source functions, which also maintains the desirable property of strong hyperbolicity of the field equations. In fact, the source functions can be thought of as representing the coordinate freedom of the Einstein equations, and when constructing solutions of the equations, via an initial value approach, for example, they must be completely specified in some fashion. Choosing the source functions in a controlled way is a key issue of the GH technique and after testing several recent prescription, we conclude that a damped wave gauge LS () is remarkably robust in collapse situations.

Applications of the GH approach in spherically-symmetric situations were studied in SorkinChoptuik () and here we employ the method in axisymmetry. In both cases it is natural to use coordinates in which the symmetries of the spacetime are explicit. However, these coordinates, are formally singular: at the origin in spherical symmetry and on the axis in axial symmetry. Thus, the field equations have to be regularized in numerical implementations; here we describe a regularization procedure that is compatible with the GH formulation.

Our numerical implementation of the GH system is a free evolution code that advances initial data by solving a set of wave equations. In addition, there are also constraint equations that must be satisfied during the evolution. Although the constraints are consistently preserved in the GH approach in the continuum limit, in numerical computations at finite resolution, constraint violations generically develop. In order to maintain stability these deviations must be damped and we discuss an effective method that achieves this.

Having in mind implications for higher-dimensional GR we
test our new code by studying gravitational
collapse of a complex scalar field
in a -dimensional Kaluza-Klein (KK) spacetime. This background,
which has a single compact extra dimension curled into a circle, is a
classical example of a higher-dimensional compactified spacetime that
in certain limits can appear four dimensional, for example when the
size of the compact dimension is small.^{3}^{3}3In this paper,
however, the size of the KK circle is arbitrary.
Assuming spherical symmetry in the infinite -dimensional portion of the space
makes the problem 2+1 that depends on time and two spatial
coordinates: one in the radial direction, and one along the KK
circle.

We perform a series of numerical simulations where the initial distribution of the scalar matter is freely specified and the outcome of the evolution depends on the “strength” of the initial data as well as on its topology. The weak data correspond to the dispersion of relatively dilute pulses, while a typical strong data configuration leads to black hole formation or nearly does so. Our preliminary results indicate a wide range of the black hole forming scenarios, including how many holes form and of what topology. For instance, a static distribution of matter centered at the axis and localized in the KK direction with the energy-density above certain threshold collapses to form a black hole with a quasi-cylindrical horizon, smeared along the extra-dimension. Data with the energy-density below this threshold evolve to form a quasi-spherical horizon centered around the initial matter distribution. By further reducing the density we find that resulting black holes become progressively smaller, and at some critical density a pulse of matter is emitted radially away from the axis and a curvature singularity develops inside it. We find that for slightly lower initial densities the evolving matter splits into several pulses and two of them individually collapse to form the singularities moving apart along the KK dimension. Finally, when the initial matter distribution becomes dilute below a certain limit, no black holes or curvature singularities are created.

In the next section we describe the class of effectively 2+1 dimensional models where our code can currently be applied. In Sec. III we present the basic formulas of the GH formulation and discuss the constraints and a method to damp their violations. We describe an axis regularization in Sec. III.1 and boundary conditions in Sec. III.2. Although, in this work we integrate full -dimensional equations, in the Appendix we describe an alternative approach—one that uses a dimensional reduction on the symmetry and integrates the reduced 2+1 equations—and compare its performance with ours. Coordinate conditions and the initial data problem are formulated in Secs. III.3 and III.4, respectively. We use several diagnostics to probe the spacetimes that we construct, including computation of asymptotic measurables and apparent horizons, and describe that in Sec. III.5. After elaborating on our numerical algorithm in Sec. IV we test its performance in Sec. V, giving detailed accounts of various aspects, such as specific coordinate choices, constraint damping, numerical dissipation, and convergence. We conclude in Sec. V.4, outline possibilities of improvement, and discuss some future prospects.

## Ii The Setup

We consider a -dimensional spacetime that possesses the isometry group. We will further assume that the corresponding Killing vectors are orthogonal to the closed hypersurface they generate. The symmetry reduces the problem to effectively 2+1, which depends on time, , and two spatial dimensions that we denote by and . The spatial coordinates can be either infinite or finite. For instance, taking and assuming asymptotic flatness correspond to the usual four-dimensional axially-symmetric situation without angular momentum, while setting and assuming periodic and infinite can describe dynamics in a five-dimensional Kaluza-Klein background.

The most general -dimensional line element with these isometries can be written as

(1) | |||||

Here is the -dimensional metric, is the metric on a unit n-sphere, , running over , and the -metric and scalar are functions of and , alone.

We take a complex, minimally coupled to gravity scalar field to represent the matter of the theory and write the total action of the system as

(2) | |||||

where is the -dimensional Newton constant.

In the next section we will describe our strategy to solving the equations derived from this action. We will focus on asymptotically flat spacetimes times a circle (having the topology ) but remark that asymptotically de Sitter (dS) and anti-de Sitter (AdS) spacetime are also included in our model (2) when the potential of the scalar field satisfies with positive and negative , respectively.

## Iii Generalized-harmonic formulation and its reduction to 2D case

In order to numerically solve the Einstein equations derived from (2), we use the generalized harmonic formulation. To make the description as self-contained as possible, we summarize below basic facts regarding the approach, (more details can be found in e.g. FP1 (); Lindblom_etal (); SorkinChoptuik ()) and adapt it to the 2+1 situation of interest.

We begin by noting that whenever isometries are present one could perform a Kaluza-Klein reduction on them, and in our case (2) the reduction yields lower-dimensional, 2+1 Einstein equations coupled to the scalar, which is related to the size of the -sphere, and the matter. Initially, we had indeed performed such a reduction and coded the reduced equations; see the Appendix A for details and comparison of the methods. However, after experimenting with the reduced and the full -dimensional versions of the equations, we found that numerical solution of the latter is generically more stable. Therefore, in what follows we adopt the unreduced approach.

The Einstein equations on a -dimensional spacetime obtained by varying the action (2) can be written in the form

(3) |

where is the Ricci tensor, is the energy-momentum tensor of the matter with trace . Hereafter we use units in which the -dimensional Newton constant satisfies .

The Ricci tensor that appears in the left-hand-side of (3) contains various second derivatives of the metric components : these second derivatives collectively constitute the principal part of , viewed as an operator on . This principal part can be decomposed into a term , plus mixed derivatives of the form . Without the mixed derivatives, (3) would represent manifestly (and strongly) hyperbolic wave equations for the Friedrich96 (). One can view the GH formulation of general relativity as a particular method that eliminates the mixed second derivatives appearing in (3); see Friedrich85 (); Garfinkle2001 (); FP1 (); FP3 (); Lindblom_etal ().

One requires that the coordinates satisfy

(4) |

where are arbitrary “gauge source functions” which are to be viewed as specified quantities, and are the contracted Christoffel symbols. One defines the GH constraint

(5) |

which clearly must vanish provided (4) holds, and then modifies the Einstein equations as follows:

(6) |

This last equation can be written more explicitly as

(7) | |||||

Now, provided that are functions of the coordinates and the metric only, but not of the metric derivatives—namely —the field equations (7) form a manifestly hyperbolic system. The source functions are arbitrary at this stage and their specification is equivalent to choosing the coordinate system (“fixing the gauge”). Determining an effective prescription for the source functions is thus crucial for the efficacy of the GH approach, and our strategies for fixing the are discussed in Sec. III.3.

After the coordinates have been chosen, we integrate the equations forward in time. Consistency of the scheme requires that the GH constraint (5) be preserved in time. The contracted Bianchi identities guarantee that this is indeed the case, since, using those identities, one can show FP1 (); Lindblom_etal () that itself satisfies a wave equation,

(8) |

Thus, assuming that the evolution is generated from an initial hypersurface on which , and constraint preserving boundary conditions are used during the evolution, (8) guarantees that for all future (or past) times.

Although the GH constraint is preserved at the continuum level, in numerical calculations, where equations are discretized on a lattice the constraint cannot be expected to hold exactly. It appears that numerical solutions of (7) can admit “constraint violating modes”, with the result that the desired continuum solution is not obtained in the limit of vanishing mesh size. However, an effective way of preventing the development of such modes in numerical calculations exists: one adds terms to the field equations that are explicitly designed to damp constraint violations (see e.g. KST ()). Following the approach of Pretorius FP1 (); FP3 () that builds on earlier works gm_sys (); Gundlachetal () we define the constraint damping terms

(9) |

and solve the modified equations of the form

(10) |

Here, is the unit time-like vector normal to the hypersurfaces, that can be written as

(11) |

and is an adjustable parameter that controls the damping timescale. Specifically, it is shown in Gundlachetal () that small constraint perturbations about Minkowski background decay exponentially with a characteristic timescale of order . We note that the constraint damping term contains only first derivatives of the metric and hence does not affect the principal (hyperbolic) part of the equations.

### iii.1 Regularization of the axis,

Having described the GH formulation, we now specialize to the symmetric case. We note first that in our coordinates (1) adapted to the symmetry the line element of the flat spacetime becomes

(12) |

and that in this case the source function (4) does not vanish but becomes

(13) |

where are angular coordinates of the sphere’s line element . Since near the axis a general spacetime is locally flat, the radial component of the source function is generically singular at , diverging as . To regularize this radial component, we thus subtract the singular background contribution by transforming

(14) |

and prescribe gauge conditions using the regular sources.

Invariance of the line element (1) under the reflection in our case implies that the metric components and are odd functions of , while and are even in . The GH constraint (4) then dictates that , regularized via (14) is an odd function of , while and are even in .

Moreover, the requirement that the surface area of an -sphere must
vanish at the axis^{4}^{4}4that is, that the radial and areal
coordinates coincide at the axis, to avoid a conical singularity
there. implies . We note that
this is an extra condition on , which thus has to satisfy both this
relation, as well as the constraint that it have vanishing radial
derivative at —specifically that .
Therefore, at we essentially have three conditions on the two
fields and . In the continuum, and given regular initial
data, the evolution equations will preserve regularity, however, in a
finite-differencing numerical code this will be true only up to
discretization errors. As a general rule-of-thumb, the number of
boundary conditions should be equal to the number of evolved variables
in order to avoid regularity problems and divergences of a numerical
implementation.

An elegant way to deal with this regularity issue involves definition
of a new variable, , ^{5}^{5}5We note that a similar
variable was introduced in lambda_ref (), also for the purpose
of regularization.:

(15) |

At the axis one then has . Therefore, after changing variables from to by using in all equations, and imposing at the axis, one ends up with a system where there is no over-constraining due to the demand of regularity at . Crucially, we note that the hyperbolicity of the GH system is not affected by the change of variables.

While we note that a more straightforward regularization method that
maintains as a fundamental dynamical variable and employs
analytical Taylor-series expansion of the equations in the vicinity of
can also be used in simulations, its reliability degrades in the
strong field regime, where the regularization that uses remains consistently
accurate.^{6}^{6}6 See, however,
SorkinChoptuik (); SorkinOren () for accurate simulations of the
strong gravity regime in spherical symmetry where Taylor-series
approach is used.

### iii.2 The field equations and Kaluza-Klein boundary conditions

With the metric ansatz (1) and the regularized source function (14), our equations become 8 equations for 8 variables: six components of the 3-metric , and real and complex scalars and correspondingly. Schematically the system can be written as

(16) | |||||

(17) | |||||

(18) |

Here ellipses denote terms that may contain the metric and it derivatives and/or the source functions, in various combinations. These equations are to be evolved forward in time starting from the initial () time slice, where values for the fields and their first time derivatives must be prescribed.

In order to completely specify the problem we have to provide boundary conditions which the above equations are subject to. In this paper we will be interested in a -dimensional Kaluza-Klein spacetime of topology , where the direction is considered periodic with asymptotic length namely that , and for the future use we also define half-period . Asymptotically, this spacetime becomes Minkowski times the compact circle and the corresponding boundary conditions are

(19) |

### iii.3 Coordinate conditions

As we have already mentioned, fixing the coordinates in the GH approach amounts to specifying the source functions . The choice that we find to perform best in our case is a variant of the damped-wave gauge (DWG) condition proposed recently in LS () (see also PretoriusChoptuik_coll ())

(20) |

where is the spatial metric whose
determinant has the factor
removed in accordance with the regularization (14), and
and are free parameters.^{7}^{7}7We note that the
spatial part of DWG is essentially a version of the popular
-driver condition LS (); gamma_dr (). Since the gauge function (20) depends only on
the metric, we can simply set without destroying
hyperbolicity of our system. Below we refer to this approach as the
algebraic DWG condition.

An alternative method that preserves the hyperbolicity was originally devised by Pretorius for binary black hole simulations FP2 (); FP3 () and has also proven useful in the studies of the gravitational collapse of scalar field in spherical symmetry SorkinChoptuik (). This strategy elevates the status of the to independent dynamical variables that satisfy time-dependent partial differential equations. The evolution equations for the are designed so that the ADM kinematic variables—lapse and shift which (implicitly) result from the time development—have certain desirable properties. For example, the equation for is tailored in an attempt to keep the value of the lapse function of order unity everywhere—including near the surfaces of the black holes—during the evolution.

One specific prescription for achieving this type of control evolves the gauge source functions according to

(21) |

where is the covariant wave operator, and
and are adjustable constants.^{8}^{8}8In certain situations it is
convenient to assume that and are given functions of
space and time rather than mere constants. For example, one might
require that the gauge driver is switched on gradually in time, or
that it be active only in certain regions, e.g. in the vicinity of a
black hole, and that its effect vanish asymptotically. Thus the
temporal source function satisfies a wave equation similar to those
that govern the metric components in the system
(10). The first term on the right-hand-side
of (III.3) is designed to “drive” to a value that
results in a lapse that is approximately . The second,
“frictional” term tends to confine to this value. For the
case of the spatial coordinates, Pretorius found that the simplest
choice of spatially harmonic gauge, was sufficient in
simulations of binary black hole collisions. A slight generalization
of this technique was considered in Scheel_etal () where instead
of using , the spatial components of the source functions are
evolved according to

(22) |

where is an additional parameter.

Variants of gauge drivers were further investigated in the recent Lindblom_etal_gauge (); SorkinChoptuik (); LS (). Specifically in LS () the following hyperbolic first-order drivers were proposed

(23) |

such that all time-independent solutions of this system satisfy , where are certain predetermined target gauge functions, for instance (20), and the parameters and are freely specified.

In Sec. V.1 we compare the performance of these strategies and argue that in collapse situations the algebraic DWG approach is the most robust of all. Namely, it does not require an extensive fine-tuning of parameters, and it is long-term stable.

### iii.4 Initial data

We now consider specification of initial data, which are values for the fields and their first time derivatives at . For simplicity we restrict attention to time-symmetric initial conditions for the metric and .

Given this assumption, initial data for the scalar field reduces to the specification of , which we take to have the form of a generalized Gaussian,

(24) |

where is a complex amplitude and and are real adjustable parameters, supplemented by the choice

(25) |

that satisfies , where is another parameter.

The momentum constraint is satisfied for our initial data, and writing the initial metric as

(26) |

one can show that the Hamiltonian constraint becomes an elliptic equation for

(27) |

that is subject to the boundary conditions, and .

We will assume initial harmonic coordinates, which implies, see (4), that the lapse can be determined in terms of as

(28) |

After substituting this into equation (III.4) we solve for and initialize our basic variables and their derivatives:

(29) |

We finally note that in order to use the DWG conditions correctly, we have to smooth out the transition from initially harmonic to later DWG coordinates. As described in Sec. V.1 this can be done stably by multiplying in (20) by a time-dependent factor that gradually grows from 0 to 1. For the driver version of DWG we also initialize at .

### iii.5 Spacetime diagnostics

We employ several diagnostics in order to probe the geometry of the spacetimes we construct.

#### iii.5.1 Asymptotic charges: Mass and tension.

Far away from an isolated system a natural radial coordinate is well defined by comparison with the flat background, and the charges that characterize the solution can be found from the asymptotic radial behavior of the metric functions. In asymptotically flat spacetimes with no angular momentum there is exactly one charge: the ADM mass. However, in our system with the extra compact dimension, the tension, associated with varying the length of the compact direction, can also be defined. A derivation of this result can be found in KolSorkinPiran1 (); HO1 (), but here we note why the appearance of an additional charge could be expected. Asymptotically the metric becomes -independent since from the lower-dimensional perspective the -dependent Kaluza-Klein modes are massive (with the discrete masses ) and decay exponentially . In this situation asymptotic Kaluza-Klein reduction is possible with the result that behaves effectively as a scalar field that carries an (unconserved) charge. Designating the asymptotic fall-off of the metric functions and as

(30) |

the mass and tension of the solutions are defined as KolSorkinPiran1 (); HO1 ()

(31) |

where is the surface area of a unit -sphere, and is the asymptotic length of the KK circle.

#### iii.5.2 Apparent horizons

It is well known that a process that concentrates sufficient mass-energy within a small enough volume can lead to the formation of a black hole. In numerical calculations based on a space-plus-time split, black hole formation is often inferred by the appearance of apparent horizons. An apparent horizon is defined as the outermost of the marginally trapped surfaces, on which future-directed null geodesics have zero divergence. Specifically, in our case we will be searching for simply connected horizons of two categories, see Fig. 1: (i) Quasi-spherical horizons of topology , in which case the surface describing the horizon can be taken as

(32) |

where

(33) | |||||

where is the point where the horizon is centered, and (ii) Horizons of topology smeared along the direction, which do not intersect the axis. In this case a convenient parametrization is

(34) |

In either case we define an outward-pointing spacelike unit normal to the surface ,

(35) |

The vanishing of the divergence, , of the outgoing null rays defined by can be expressed as

(36) |

Substituting the expressions (32) or (34) into (36) yields ordinary second order differential equations for or correspondingly. The equations could be solved by “shooting”, subject to appropriate boundary conditions, however, here we instead use the following point-wise relaxation method, that is more suitable for parallel numerical implementations. We start with supplying an initial guess for the entire function and iterate the parabolic equation

(37) |

in unphysical “time” until a solution is found. This equation implies that depending on whether is inside or outside the horizon at any given moment, it will expand or shrink in the next instant. In numerical simulations the hope is that the initial will “flow” to the apparent horizon in finite time. We find that in our case this indeed happens when the initial guess is “reasonably close” to the final solution.

## Iv Numerical approach

Here we describe our strategy for the numerical solution of the GH system (with a scalar matter source) in Kaluza-Klein spacetime.

### iv.1 The numerical grid and the algorithm

We cover the -- space by a discrete lattice denoted by
, where and
are integers and , , and define the
grid spacings in the temporal and two spatial directions,
respectively. As described in the next section, the spatial domain is
compactified, and hence a grid of finite size extends
from the axis to spatial infinity. Approximations to the dynamical
fields, collectively denoted here by , are evaluated at each grid
point, yielding the discrete unknowns . In the
interior of the domain, the GH equations and the gauge-driver
equations are discretized using^{9}^{9}9Here stands for any of
the mesh-spacings and that define our
numerical grid. finite difference approximations
(FDAs), which replace continuous derivatives with the discrete
counterparts, examples of which are given in Table 1. As
in FP1 (); FP2 (); SorkinChoptuik () our scheme directly integrates the
second-order-in-time equations.

Centered derivatives | |
---|---|

( | |

( | |

( | |

( | |

Backward derivatives | |

Following discretization, we thus obtain finite difference equations at every mesh point for each dynamical variable. Denoting any single such equation as

(38) |

we then iteratively solve the entire system of algebraic equations as follows.

First, we note that for those variables that are governed by equations of motion that are second order in time, our discretization of the equations of motion results in a three level scheme which couples advanced-time unknowns at to known values at retarded times and . In order to determine the advanced-time values for such variables, we employ a point-wise Newton-Gauss-Seidel scheme: starting with a guess for (typically, we take ) we update the unknown using

(39) |

Here, is the residual of the finite-difference equation (38), evaluated using the current approximation to , and the diagonal Jacobian element is defined by

(40) |

In the cases where we used gauge drivers (III.3) we found that an iteration based on the Crank-Nicholson discretization scheme of the corresponding first order equations performed well. Specifically, writing any such equation schematically as , we update using

(41) |

We iterate (39) and (41) over all equations until the total residual norm, see (49), falls below a desired threshold.

In order to inhibit high-frequency^{10}^{10}10“High-frequency” refers
to modes having a wavelength of order of the mesh spacings, and . instabilities which often plague FDA
equations, our scheme incorporates explicit numerical dissipation of
the Kreiss-Oliger (KO) type KO (). Following FP1 (), at the
interior grid points, , and for each dynamical variable we apply a low-pass KO
filter by making the replacement

(42) | |||||

at both the and time-levels before updating the unknowns. Here is a positive parameter satisfying that controls the amount of dissipation. An extension of the dissipation to the boundaries FP1 () was also tried, but this has not resulted in any significant improvement of the performance of the code.

### iv.2 Coordinates and boundary conditions

While the physical, asymptotically flat (times a circle, in our case) spacetime extends to spatial infinity, in a numerical code one can only use grids of finite size. Instead of using a standard strategy that deals with this issue by introducing an outer boundary at some finite radius where approximate boundary conditions are imposed, we adopt another technique—proven successful in previous work in numerical relativity, see e.g. FP2 (); SorkinChoptuik (); Choptuik_BS ()—and compactify the spatial domain. We find that compactifying the radial direction and imposing the (exact) Dirichlet conditions (19) at the edge of the domain works well, provided that we use sufficient dissipation. In particular, it is known that due to the loss of resolution near the compactified outer boundary (assuming a fixed mesh spacing in the compactified coordinates), outgoing waves generated by the dynamics in the interior will be partially reflected as they propagate toward the edge of the computational domain, and these reflections will then tend to corrupt the interior solution. By adding sufficient dissipation one can damp the waves in the outer region, attenuating any unphysical influx of radiation. This enables one to use the compactification meaningfully.

The results presented in this paper were obtained using the compactification of the form

(43) |

where the compactified ranges from to for values of the original radial coordinate . In practice, we define a uniform grid in the compactified and use chain-rule to replace derivatives with respect to with the derivatives with respect to in all dynamical equations. The asymptotic boundary conditions (19) at are then imposed exactly: , and .

We have previously described the boundary (regularity) conditions at in Sec. III.1. Denoting by , for , the advanced-time value at the axis for any of the variables, ,, and that have vanishing derivative at , we use the update , which is based on a second-order backwards difference approximation (see Table 1) of . For the quantities and , which are odd in as , we simply use .

For simplicity, in this paper we consider only configurations that have reflection symmetry about which together with periodicity in z, implies

(44) |

We update points at and using backward difference approximation similar to that in Tab. 1. In particular for , we use at , and at .

### iv.3 The elliptic equation of initial data

It turns out that generally the equation (III.4) for the conformal factor is ill-posed, namely that the linear equation governing small perturbations about a solution of (III.4) does not admit a unique solution for given boundary conditions York (). Therefore, an attempt to solve (III.4) using standard relaxation methods will generically fail, as the relaxation is not guaranteed to converge. However, a method circumventing this difficulty exists York (), and this is through a rescaling that transforms the Hamiltonian constraint (III.4) into

(45) | |||

By choosing the power such that terms in this equation which are proportional to have only non positive ’s one renders the problem well posed York (). In the case of a free scalar field having the potential , which we assume here, this implies that , see a related discussion in SorkinPiran ().

Note that the physical data are , not , and therefore we solve equation (IV.3) in an iterative manner where we start with an initial guess and at each iteration update that is then used in (IV.3) in order to solve for . Most of the results obtained in this paper were generated using and . Interestingly, it turns out that provided the initial guess is not too distant from the solution, the original equation (III.4) is numerically stable without the rescaling. This happens, for instance, in the weak field regime where the guess relaxes to the solution without any trouble.

### iv.4 Excision

We use excision to dynamically exclude from the computational domain a region interior to the apparent horizon that would eventually contain the black hole singularity. This approach relies on the observation that in spacetimes that satisfy the null energy condition and assuming the cosmic censorship holds, the apparent horizon is contained within the event horizon. This ensures that the excluded region is causally disconnected from the rest of the domain (see thornburg_excision () and the references therein for further discussion). Operationally, once an apparent horizon is found, we introduce an excision surface, , contained within the apparent horizon, and such that all characteristics at are pointing inwards. This specific characteristic structure eliminates the need for boundary conditions at : rather, advanced-time unknowns located on the excision surface are computed using finite difference approximations to the interior evolution equations, but where centered difference formulas are replaced with the appropriate one-sided expressions given in Table 1.

Currently our apparent horizon finder is only capable of locating horizons of the shapes depicted in Fig. 1. In this case the radius of the excision surface typically satisfies , where is the coordinate radius.

## V Performance of the code

In this section we investigate the performance of the code in series of simulations that evolve regular initial distribution of complex scalar field in five and six dimensional Kaluza-Klein spacetime. The code uses pamr/amrd infrastructure pamr_amrd () where our suitably interfaced numerical routines are called by the amrd driver. All our results are generated using an initial scalar field profile of the Gaussian form (24), with fixed values and , so that the scalar pulse is always initially centered at the axis, see Fig. 2.

In addition, we set and use , unless otherwise specified. We choose the initial frequency in (25) to be , and the asymptotic size of the KK circle, . The scalar field potential used here is taken to be of the form , where without much loss of generality we set .

Because we mostly use a time-explicit finite difference scheme, we expect restrictions on the ratio (the Courant factor) that can be used while maintaining numerical stability. In the results discussed below we chose in the weak-field regime, and for evolving the strong data. Taking larger usually leads to amplification and dominance of numerical errors near the axis, that shows up as diverging high frequency oscillations. In most cases we use uniform grids of the same size in both spatial directions, , and our lowest and highest resolution simulations have and , respectively. The highest-resolution runs generally required for stability.

As expected, the outcome of the collapse depends on the strength of the initial data. Low density distributions describe weakly gravitating scalar pulses which completely disperse in all instances. Strong data generate spacetimes in which black holes form, or almost form. For fixed and a single free parameter that controls the strength of the pulse—and hence, the outcome of the evolution—is the initial amplitude . In this case “strong initial data” will refer to situations when an increase of the initial amplitude by less than leads to formation of a curvature singularity, and the other data will be called “weak”. The critical amplitude is at the threshold for the singularity formation.

For each spacetime that we construct, the asymptotic mass and tension are computed using (31), where the constants and are found by fitting the metric components and with the functions of the form in the asymptotic region. The errors in the constants are determined by the fitting uncertainty which in our case is (larger for weaker data). We find that our initial data sets usually have , and it follows from (31) that , within the numerical accuracy of our method.

A useful quantity that illustrates energy-momentum distribution is the density . Figure 3 depicts at several moments during the 5D evolution of weak initial data defined by , using algebraic damped wave gauge (DWG), see (20). Figure 3 shows that in the early stages of the evolution the pulse (initially localized around the center ) undergoes a gravitational collapse that is accompanied by a growth of the central energy-density. At a later time, however, the pulse bounces off while forming several shells, and disperses. The distribution of the energy density is anisotropic inside the shells, with a higher concentration of energy occurring along the axis, , and at the equator, . The matter waves that travel along the compact KK circle collide to form a typical interference picture (see Fig. 3) and the pattern becomes increasingly complicated in the course of time, when more and more waves undergo interaction. As expected in the KK background the -dependence of solutions vanishes in the asymptotic region, since the -dependent modes are massive and fall off exponentially fast.

In order to obtain more quantitative insight into the process, we plot in Fig. 4 the evolution of the logarithm of at the center. During the highly dynamical early epoch, lasting until , the field collapses, bounces off the center, spreads along the -direction, and starts dispersing to infinity. In the first stage of the dispersion, lasting until , the decay of follows a power-law, and in the late times the decay is exponential with a characteristic time scale of a few hundreds of . The maximal scalar curvature is attained during the initial collapse phase; in the shown simulation the curvature reaches the magnitude of order of (in units of ). The discrete spectrum of the high-frequency oscillations consists of the normal frequencies defined by the size of the KK circle. We observe that the higher frequency modes decay faster, and the late-time spectrum is dominated by a few lowest frequency modes, see also bottom right panel of Fig. 3.

When the initial amplitude of the scalar field increases, the maximal curvature achieved during the collapse grows, and when surpasses a certain threshold, the curvature diverges, signaling the appearance of a singularity; see Fig. 5. Currently we are able to estimate the critical amplitude with a relatively low precession of approximately 1 part in 800. In five dimensions the amplitude is between and such that initial data with completely disperses, and the data gives rise to curvature singularities and black holes.

A covariant way to illustrate the distribution of matter
and the geometry of the spacetime is provided by the Ricci^{11}^{11}11In this and other figures Ricci scalar is
measured in units of . scalar, and
Fig. 6 shows its evolution computed in the
evolution of our strongest subcritical data set,
. The initial pulse collapses,
bounces off and forms several shells, that start dispersing after
. While the anisotropy of the energy distribution inside
the shells is small for weak data, it is obvious in the strong field regime.
Some matter is ejected along the equator, , and two
distinctive pulses, moving in the opposite directions, form at the axis.
The pulses collide and interact, which results in smearing of the energy-density along the axis,
see bottom right panel of Fig. 6. The distribution is not
stationary, rather the matter is continuously leaking to infinity,
and the late-time behavior is qualitatively similar to the weak field case,
shown in Fig. 3.

It turns out that for higher initial amplitudes, , the pulses that form at the axis are able to individually collapse and develop curvature singularities, indicated by diverging scalar curvature and energy-momentum density while the lapse remains finite (of order one) everywhere. Presently, our horizon-finder is not designed to locate the moving apparent horizons that may arise in this case around the singularities, and it would be very interesting to verify whether or not such engulfing horizons indeed form. The time when the curvature singularities appear depends on the strength of the initial data in the manner that the closer to threshold we are the later into the evolution the curvature diverges. For instance, in the case of the pulses, collapse and become singular as they cross the circle and collide at , see Fig 7. However, for , the curvature inside each pulse blows up earlier, when they reach .

Intriguingly, we observe that for even stronger initial data formation of the curvature singularity ensues differently. Specifically, the pulse of matter which is usually emitted during the initial collapse-bounce stage outwards along the equator, is now seen to also be able to collapse and develop a singularity. Snapshots of the process are shown in Fig. 8 that was obtained in the evolution of the data defined by . Since all our attempts to detect an apparent horizon, engulfing the singular region and the center, have failed we believe that the horizon, if it forms in this case, must be localized around the moving curvature singularity.

In Fig. 9 we plot the Ricci scalar along the equator at , shortly before the appearance of a singularity at causes the code to break down. Figure 9 shows that the maximal value of the Ricci scalar along this slice is determined by the numerical resolution: the finer meshes we use, the larger the values attained. This signals emergence of a curvature rather than a coordinate singularity.

Figure 10 shows snapshots of the energy-density distribution obtained in a supercritical simulation initiated by . The resulting spacetime has the total mass and the tension . The initial stages of this evolution are qualitatively similar to those of weaker data, however, after the first bounce off the center (at ) the pulse recollapses and a black hole forms as indicated by the appearance of an apparent horizon. The energy-density and scalar curvature both blow up in the vicinity of in finite time, while the lapse and shift remain finite.

In order to locate apparent horizons we solve (37) iteratively starting with some initial guess for the entire function that parametrizes the horizon. In the case of we are searching for a horizon that has spherical topology. We use the parametrization (32) with and solve (37) initialized by . Setting the target accuracy to and using grid points to represent the horizon, we are able to solve the equation in iterations. The resulting horizon is shown in bottom panels in Fig. 10.

After an apparent horizon is found we use excision to remove a region inside the horizon that contains curvature singularity. However, it turns out that regardless of our specific gauge choice the lapse keeps evolving inside the horizon, and continuously decreases until it is reaching the magnitudes of order near the excision boundary. In this situation truncation errors in quantities near occasionally cause the computation of non-positive values for the lapse, which immediately leads to code failure. Unfortunately, this happens when the apparent horizon is still fairly dynamical and continues to change its shape and size. Therefore, we are currently not able to determine the stationary black hole state in the end of the evolution.

Having described the low mass configurations we will now briefly discuss more massive solutions of certain type. One initial configuration that we have evolved consisted of a nearly uniform along distribution of the matter—achieved by setting in (24)—with the initial amplitude of . The Ricci scalar and the energy-density computed in this evolution are shown in Fig. 11 together with the resulting apparent horizon that appears at and has a cylindrical topology. In this simulation we used the algebraic DWG condition, and in this case we were not able to locate a surface on which all characteristics will point inwards. Therefore, no excision was employed in this evolution, and the eventual failure of the code was caused by collapse of the lapse at the axis.

While the total mass of this spacetime, , is determined by the asymptotic fall-off (31), the mass of the “black string” (a black hole with the smeared horizon) can be estimated from the average size of its horizon. The last moment of the evolution before the code had crashed is shown in bottom panels of Fig. 11. The horizon is nearly uniform along and is located at , which yields the mass of , where is the uncompactified average areal radius of the horizon. This value is below the critical mass, , needed for stability of the extended solutions of this type GL (). Therefore, at this stage the system is probably far from approaching stationarity. Since the total available mass is higher than several scenarios for subsequent dynamics can be imagined. If the black hole accretes enough matter and increases its mass above critical, it can reach the uniform black-string end-state. Otherwise, the evolution would probably have to proceed via forming a localized black hole in the first place. This black hole then may or may not accrete additional matter, and as a result either to grow and become a black string or to settle down to a stationary solution of spherical topology. While further investigation of the process is clearly needed, Fig. 11 shows a developing progressive localization of energy-density and curvature around , indicating that formation of a localized black hole first is more probable in this specific case.

We turn now to a detailed description of the performance of the code. Since we have implemented a free evolution scheme, we can assess the convergence of our numerical solutions by monitoring discrete versions of the Hamiltonian and momentum constraints, which are defined by contracting the Einstein equations with the unit normal vector to the hypersurfaces, i.e. , where is the Einstein tensor. One way of doing this involves evaluation of the following -norm of finite-differencing variables

(46) |

In the next section, we compute the -norms of the constraints at each time step and examine their behavior as a function of various parameters of the problem.

### v.1 Coordinate conditions

In the weak-gravity regime where an initial pulse of matter completely disperses to infinity we find clear advantage of the algebraic DWG conditions. They are robust, almost do not require fine-tuning, and are extremely stable. Additionally, as described in the next section, the constraints remain well preserved in this case, such that only very small damping (controlled by the parameter , see (9)) needs to be added to the equations.