An axisymmetric generalized harmonic evolution code

An axisymmetric generalized harmonic evolution code

Evgeny Sorkin Max Planck Institute for Gravitational Physics (Albert Einstein Institute) Am Muehlenberg 1, D-14476, Golm, Germany

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. 111See 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. 222and 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.333In 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


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


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


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


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


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


This last equation can be written more explicitly as


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,


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


and solve the modified equations of the form


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


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


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


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


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 axis444that 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, 555We note that a similar variable was introduced in lambda_ref (), also for the purpose of regularization.:


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.666 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


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


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


where is the spatial metric whose determinant has the factor removed in accordance with the regularization (14), and and are free parameters.777We 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


where is the covariant wave operator, and and are adjustable constants.888In 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


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


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,


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


that satisfies , where is another parameter.

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


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


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


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


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


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


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




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

Figure 1: We are locating apparent horizons of two types: one, having a spherical topology, is conveniently parametrized by (plotted here as centered at ), and another, having a cylindrical topology, is parametrized by .

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


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


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


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 using999Here 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
Table 1: Examples of second order finite-differencing approximation to derivatives calculated at a grid point that we use in our discretization scheme in the interior of the domain (centered stencil), and at the excision boundary (backward stencil) of the numerical grid.

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


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


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


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


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-frequency101010“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


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


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


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


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.

Figure 2: We choose the initial distribution of the scalar field as a generalized Gaussian (24) and assume that the initial time derivative is given by (25). We mostly consider localized initial pulses (such as the one shown here), that are obtained by setting in (24). The radial direction is compactified, and the hypersurfaces and identified in a Kaluza-Klein spacetime under consideration.

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.

Figure 3: Snapshots of the evolution of energy-momentum density in the simulation of the initial data defined by in 5D, using algebraic DWG. The horizontal axis represents the compactified radial direction, and the vertical axis is along the periodic -direction. The lapse of time is measured in units of the total mass. The initial collapse of the pulse around the center is accompanied by a growth of the amplitude of by about an order of magnitude. The pulse then bounces off and several outgoing shells form. The shells propagate along the periodic -direction and interact, creating a typical interference pattern. As expected, the -dependence vanishes in the asymptotic regions. In the later stages of the evolution is small, and the frequency spectrum of the oscillations along is dominated by a few lowest eigenfrequencies associated with the compact KK circle.

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.

Figure 4: The logarithm of the central energy-density, , in a five-dimensional weak field simulation defined by that uses algebraic DWG. The resulting spacetime has the mass and the tension is zero within the numerical accuracy. In the early stages of the evolution, the pulse collapses and grows. Subsequently, the energy density decreases following a power-law dependence during the period lasting from until . After that decays exponentially . The high-frequency oscillations are associated with the eigenmodes of the compact KK circle.

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.

Figure 5: A log-log plot of the Ricci scalar (in units of ) as a function of the proper time, , both evaluated at the center, in several simulations in . The first peak in the curves corresponds to the initial collapse phase, and for stronger data the Ricci scalar diverges in finite time, signaling a curvature singularity (case ). In subcritical cases, the curvature decreases after the initial growth, and the secondary peaks correspond to collisions of the shells traveling along the KK circle and arriving back at .

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 Ricci111111In 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.

Figure 6: The Ricci scalar at several instants of the evolution of the initial data characterized by , currently our strongest subcritical data set. The numerical resolution used in this simulation is . After the initial collapse stage an expanding shell of matter forms. The energy-density distribution inside the shell quickly becomes anisotropic with most of the energy localized within individual pulses traveling along the axis and inside the pulse emitted radially along the equator, . The pulses at the axis undergo interactions and spread the energy-density along . In late stages of the evolution, all the matter is radiated away to infinity.
Figure 7: Snapshots of the Ricci scalar at several instants of the evolution of the initial data characterized by at the resolution of . The matter distributes anisoropically inside the bouncing shells such that two localized pulses, traveling along the axis, form. When the pulses collide at a curvature singularity appears.
Figure 8: The Ricci scalar at several instants during the evolution defined by . A curvature singularity forms along the equator, inside the lump of matter emitted in the radial direction during the initial collapse-bounce process. This simulation uses the resolution of .

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 9: The Ricci scalar (in units of ) along the equator at , shortly before the evolution defined by develops a singularity at . The maximal value that the Ricci scalar attains in this evolution depends on the resolution (here ), being larger on finer numerical meshes, which indicates the geometrical nature of the emerging singularity.
Figure 10: The distribution of energy-density in the evolution of the initial data, defined by in a simulation that uses the resolution of . A strong gravitational interaction in this case precludes any significant amount of matter from escaping to infinity and leads to formation of a black hole soon after the first bounce. The black hole is centered around and is signaled by the appearance of an apparent horizon (designated by a thick dashed line) and an explosive growth of the scalar curvature and the energy-density.

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.

Figure 11: Distributions of the energy-density (on the right) and the Ricci scalar (on the left) computed in the evolution of a nearly uniform along initial data defined by , and , that uses the resolution of . A strong gravitational interaction leads to a rapid formation of a black string, signaled by the appearance of an apparent horizon of cylindrical topology (shown as a thick dashed line). Note, that the matter and the curvature inside the horizon are strongly localized around the center .

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


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.

Since at , we assume harmonic conditions , choosing the source functions according to (