format: text

text:
# Hyperbolic slicings of spacetime: singularity avoidance and gauge shocks

## I Introduction

The choice of good coordinates plays a crucial role in finding solutions of the Einstein equations. This is particularly important in the case of numerical simulations of strongly gravitating systems, where a bad coordinate choice can easily lead to the formation of a coordinate singularity which stops the numerical simulation and severely limits the region of the spacetime covered by it. Coordinate singularities are not the only concern: the presence of physical singularities, such as those associated with black holes, can also have a deep impact in a numerical simulation, as unless special care is taken the time slices can march right onto the physical singularity very early during a simulation.

When considering dynamical evolutions of spacetime based on a 3+1 decomposition Arnowitt62 ; York79 , the coordinate choice naturally separates in two different aspects: the choice of a specific foliation of spacetime into spatial hyper-surfaces (also known as the “slicing”), associated with the lapse function $`\alpha `$, and the choice of the way in which the lines of constant spatial coordinates (the “time lines”) propagate from one hyper-surface to the next, associated with the shift vector $`\beta ^i`$. Here I will concentrate fully in the role played by the choice of a slicing condition and leave the choice of a shift vector for a later work.

In order to specify a foliation of spacetime, one needs to prescribe a way to calculate the lapse function $`\alpha `$, which measures the proper time interval between neighboring hyper-surfaces along their normal direction. There is, of course, an infinite number of ways in which one can choose the lapse function, but typically the different choices can be classified in the following way: 1) Prescribed slicings, where the lapse is specified as an a priori known function of space and time, 2) algebraic slicing conditions, where the lapse is specified algebraically as some function of the geometric variables (metric and extrinsic curvature) at each hyper-surface, 3) elliptic slicing conditions, where the lapse is obtained by solving an elliptic differential equation at every time step that typically enforces some geometric condition on the spatial hyper-surfaces, and 4) time derivative slicing conditions, where the time derivative of the lapse is specified as some algebraic function of the geometric quantities and the lapse is evolved as just another dynamical quantity (this last case is often included in the algebraic slicing class mentioned above).

An example of a prescribed slicing is the so-called “geodesic slicing”, where one simply takes $`\alpha =1`$. An example of an algebraic slicing condition is “harmonic slicing” where one takes instead $`\alpha =\sqrt{\gamma }`$, with $`\gamma `$ the determinant of the spatial metric. A well known elliptic slicing condition is the “maximal slicing” condition Smarr78b , which requires that the spatial volume elements remain constant during the evolution. Elliptic conditions are typically robust and well behaved, but have the drawback of being computationally expensive. Algebraic conditions are much easier to apply, but are difficult to analyze in a general case. Time derivative slicing conditions, on the other hand, have the advantage of being both easy to implement and, in the particular case when they lead to hyperbolic equations (as is the case of the Bona-Masso family Bona94b , see below), much easier to analyze as well. They also include many well known algebraic conditions as their integral form.

In this paper, I will only consider hyperbolic slicing conditions and study in which ways such conditions can lead to pathological slicings and how can those pathologies best be avoided. There are many different ways in which a foliation of spacetime can become pathological: The slices can hit a physical singularity, the slices can hit a coordinate singularity where the volume elements become zero (the normal lines focus), the slices can become non-smooth at some point or the slices can remain smooth but stop being space-like (they can become null at a point, for example). Of the different possible pathologies mentioned above, I will concentrate on two specific types: “focusing singularities” Bona88  defined as those for which the spatial volume elements vanish at a bounded rate, and “gauge shocks” Alcubierre97a  defined as solutions for which the lapse becomes discontinuous as a consequence of the crossing of the characteristic lines associated with the propagation of the gauge, and the time slices therefore become non-smooth.

This paper is organized as follows: In Sec. II, I introduce the Bona-Masso hyperbolic slicing condition. Focusing singularities are defined in Section III where I also find under which circumstances the Bona-Masso slicing condition is singularity avoiding. In Sec. IV, I introduce the idea of a gauge shock, derive a condition to avoid them and see what slicings obey that condition either exactly or approximately. I conclude in Sec. V.

## II The Bona-Masso family of hyperbolic slicing conditions

The Bona-Masso family of slicing conditions Bona94b  is a time evolution type of condition for which the lapse is chosen to satisfy the following evolution equation

$$\frac{d}{dt}\alpha \left(_t_\beta \right)\alpha =\alpha ^2f(\alpha )K,$$
(1)
with $`_\beta `$ the Lie derivative with respect to the shift vector $`\beta ^i`$, $`K`$ the trace of the extrinsic curvature and $`f(\alpha )`$ a positive but otherwise arbitrary function of $`\alpha `$. The reason why $`f(\alpha )`$ has to be positive will become clear below. Here we just mention the fact that the right hand side of condition (1) is the most general term one can construct that involves only first order spatial scalars Bona94b .

There are several things that are important to mention about the family of slicing conditions (1). First, we notice that even if this family was proposed in the context of the Bona-Masso hyperbolic reformulation of the Einstein equations Bona89 ; Bona92 ; Bona93 ; Bona94b ; Bona97a , it is in fact quite general and can be used successfully with any particular form of the 3+1 evolution equations. This was shown recently when such a condition was used together with the Baumgarte-Shapiro Shibata-Nakamura (BSSN) formulation Baumgarte99 ; Shibata95  to obtain long-term stable and accurate evolutions of black hole spacetimes Alcubierre01a ; Alcubierre02a . Also, condition (1) is a generalization of slicing conditions that have been used in evolution codes based on the Arnowitt-Deser-Misner (ADM) formulation Arnowitt62 ; York79  since the early 90’s Bernstein93a ; Anninos94c .

Second, condition (1) can also be trivially adapted to the case when, instead of the lapse $`\alpha `$, one evolves a densitized lapse of the form

$$Q:=\alpha \gamma ^{\sigma /2}$$
(2)
with $`\gamma `$ the determinant of the spatial metric $`\gamma _{ij}`$ and $`\sigma `$ a constant parameter. Such a densitized lapse (particularly the case $`\sigma =1`$) has recently been advocated in the context of hyperbolic reformulations of the Einstein equations (see for example Anderson99 ; Kidder01a ). In terms of $`Q`$, condition (1) becomes

$$\frac{d}{dt}Q:=\frac{Q^2}{\gamma ^{\sigma /2}}\left(f+\sigma \right)K,$$
(3)
where to calculate the Lie derivative with respect to $`\beta ^i`$ contained in the operator $`d/dt`$ one must use the fact that $`Q`$ is a density of weight $`\sigma `$.

Finally, it is also important to mention that the shift terms included through the Lie derivative in condition (1) are such that one is guaranteed to obtain precisely the same spacetime foliation regardless of the value of the shift vector. This would seem to be a natural requirement for any slicing condition. However, it is plausible that in a particular situation one would like to choose a slicing condition and a shift vector that are closely interrelated (see for example Balakrishna96a ). Indeed, generalizations of the Bona-Masso slicing condition that in the presence of a non-zero shift vector do not have the form (1) have already been used in the literature. For example, Refs. Alcubierre00a ; Yo02a  use $`_t\alpha =\alpha f(\alpha KD_i\beta ^i)`$ instead of (1) (with $`D_i`$ the spatial covariant derivative), and Ref. Alcubierre02a  simply uses $`_t\alpha =\alpha ^2f(KK_0)`$ (with $`K_0=K(t=0))`$ even with a non-zero shift vector. Which off these generalizations is best under different circumstances is an important question that I will leave for a future work.

With these comments in mind, let us go back to condition (1). Taking an extra time derivative we find

$$\frac{d^2}{dt^2}\alpha =\alpha ^2f\left[\frac{d}{dt}K\alpha (2f+\alpha f^{})K^2\right],$$
(4)
with $`f^{}:=df/d\alpha `$. From the ADM evolution equations one easily finds that, in vacuum,

$$\frac{d}{dt}K=D^2\alpha +\alpha K_{ij}K^{ij},$$
(5)
with $`D^2`$ the Laplace operator associated with the spatial metric. This last equation implies

$`{\displaystyle \frac{d^2}{dt^2}}\alpha \alpha ^2fD^2\alpha =`$
$`\alpha ^3f\left[K_{ij}K^{ij}\left(2f+\alpha f^{}\right)K^2\right].`$ (6)
Equation (6) shows that the lapse obeys a wave equation with a quadratic source term in $`K_{ij}`$. It is because of this that we say the slicing condition (1) is a hyperbolic slicing condition: it implies that the lapse evolves with a hyperbolic equation (but see the discussion on gauge shocks below for a more formal proof of hyperbolicity). The wave speed associated with equation (6) along a specific direction $`x^i`$ can be easily seen to be

$$v_g=\alpha \sqrt{f\gamma ^{ii}}.$$
(7)
Notice that this will only be real if $`f(\alpha )0`$, which explains why we asked for $`f(\alpha )`$ to be positive. In fact, $`f(\alpha )`$ must be strictly positive because if it was zero there would be no complete set of eigenvectors and we would not have a strongly hyperbolic system (see Sec. IV.2).

To see how the gauge speed $`v_g`$ is related to the speed of light consider for a moment a null world-line. It is not difficult to find that such a world-line will have a coordinate speed along direction $`x^i`$ given by

$$v_l=\alpha \sqrt{\gamma ^{ii}},$$
(8)
so the gauge speed (7) can be smaller or larger that the speed of light depending on the value of $`f`$.

Notice that having a gauge speed that is larger than the speed of light does not lead to any causality violations, as the superluminal speed is only related with the propagation of the coordinate system. One could argue that superluminal gauge speeds are not desirable as they would allow gauge effects to propagate out of black hole horizons, for example. Empirically, however, the most successful slicing conditions for the simulation of black hole spacetimes have been precisely those that have superluminal (even extremely large) gauge speeds: maximal slicing, which as an elliptic condition has, at least formally, an infinite speed of propagation, and the 1+log slicing condition which is a member of the Bona-Masso family that has superluminal gauge speeds whenever the lapse is small (see following section).

### II.1 Relating lapse and spatial volume elements

The ADM evolution equation for the spatial metric $`\gamma _{ij}`$ is given by

$$\frac{d}{dt}\gamma _{ij}=2\alpha K_{ij},$$
(9)
which implies the following evolution equation for the spatial volume elements $`\gamma ^{1/2}`$:

$$\frac{d}{dt}\gamma ^{1/2}=\alpha \gamma ^{1/2}K.$$
(10)
Comparing this with equation (1), and taking $`f=1`$, we can trivially solve for $`\alpha `$ in terms of $`\gamma ^{1/2}`$ to obtain (in the case of zero shift vector)

$$\alpha =h(x^i)\gamma ^{1/2},$$
(11)
with $`h(x^i)`$ a time independent function. It is very important to stress the fact that the previous relation holds only when moving along the normal direction to the hypersurfaces, and not when moving along the time lines which will differ from the normal direction for any non-zero shift vector. That is, we are relating the lapse to the volume elements as seen by the normal observers. In the following, whenever we relate the lapse to the volume elements, it should always be understood that we are referring to the volume elements associated with the normal observers.

It is not difficult to show (see following section) that equation (11) is equivalent to the condition

$$\mathrm{}t=g^{\mu \nu }\mathrm{\Gamma }_{\mu \nu }^0=0,$$
(12)
with $`g_{\mu \nu }`$ the spacetime metric. That is, $`f=1`$ corresponds to the case when the time coordinate is a harmonic function. Because of this the case $`f=1`$ is known as “harmonic slicing”. Notice also that in this case the gauge speed is identical to the physical speed of light, i.e. the gauge propagates along null lines. Harmonic slicing can be seen to be equivalent to having a time independent lapse density $`Q`$ of weight $`\sigma =1`$ (again, assuming the shift vanishes).

One can construct other well known families of slicing conditions by choosing different forms of $`f(\alpha )`$. For example, if we choose $`f(\alpha )=N`$, with $`N`$ a constant, we obtain what can be called the “generalized harmonic slicing condition”, which can also be easily integrated to give

$$\alpha =h(x^i)\gamma ^{N/2}.$$
(13)
And if we take $`f(\alpha )=N/\alpha `$ we obtain the “1+log” family Abrahams92a ; Anninos93b , which again can be integrated to find

$$\alpha =h(x^i)+\mathrm{ln}\left(\gamma ^{N/2}\right).$$
(14)
The 1+log family mimics maximal slicing and has strong singularity avoiding properties (see section III below). In particular, 1+log slicing with $`N=2`$ has been found empirically to be very robust when evolving black hole spacetimes Arbona99 ; Alcubierre01a ; Alcubierre02a . As mentioned before, one can easily see that the gauge speed associated with the 1+log family can become far larger than the speed of light as the lapse becomes smaller ($`v_g/v_l=\sqrt{N/\alpha }`$).

More generally, using equation (10) one can find that for arbitrary $`f(\alpha )`$ the following relation between $`\alpha `$ and $`\gamma ^{1/2}`$ holds

$$d\mathrm{ln}\gamma ^{1/2}=\frac{d\alpha }{\alpha f(\alpha )},$$
(15)
which implies

$$\gamma ^{1/2}=F(x^i)\mathrm{exp}\left\{\frac{d\alpha }{\alpha f(\alpha )}\right\}.$$
(16)
with $`F(x^i)`$ again a time independent function. This last expression will be the starting point when we discuss focusing singularities below.

### II.2 The foliation equation

Consider now a spacetime with metric $`g_{\mu \nu }`$, and assume that we have a foliation of this spacetime into spatial hypersurfaces. Such a foliation allows us to define a time function $`T`$ whose level sets correspond to the members of the foliation. One can show that the Bona-Masso slicing condition (1) can be written as a generalized wave equation for the time function $`T`$ in the following way

$$\left[g^{\mu \nu }+\left(1\frac{1}{f(\alpha )}\right)n^\mu n^\nu \right]T_{;\mu \nu }=0,$$
(17)
with $`n^\mu `$ the unit normal vector to the spatial hypersurfaces. I will call equation (17) the “foliation equation”. Notice that for harmonic slicing the above equation reduces to the simple wave equation (12).

The fact that the foliation equation above is equivalent to the Bona-Masso slicing condition (1) can be shown by choosing 3+1 coordinates $`\{t,x^i\}`$ adapted to the foliation. In this coordinate system we can take $`T=t`$, so equation (17) becomes

$$\left[g^{\mu \nu }+\left(1\frac{1}{f(\alpha )}\right)n^\mu n^\nu \right]\mathrm{\Gamma }_{\mu \nu }^0=0.$$
(18)
Using now the 3+1 expressions for the components of the 4-metric $`g_{\mu \nu }`$ we obtain the following expressions for the Christoffel symbols

$`\mathrm{\Gamma }_{00}^0`$ $`=`$ $`{\displaystyle \frac{1}{\alpha }}\left(_t\alpha +\beta ^i_i\alpha \beta ^i\beta ^jK_{ij}\right),`$ (19)
$`\mathrm{\Gamma }_{0i}^0`$ $`=`$ $`{\displaystyle \frac{1}{\alpha }}\left(_i\alpha \beta ^mK_{im}\right),`$ (20)
$`\mathrm{\Gamma }_{ij}^0`$ $`=`$ $`{\displaystyle \frac{1}{\alpha }}K_{ij}.`$ (21)
Substituting this into (18), and using the 3+1 expressions for the normal vector $`n^\mu `$, we find

$$_t\alpha \beta ^i_i\alpha +\alpha ^2f(\alpha )K=0,$$
(22)
which is precisely the Bona-Masso condition (1).

Notice that if we take $`f>1`$, one can always have a unit normal vector $`n^\mu `$ such that the coefficient of the $`_t^2T`$ term in the foliation equation changes sign, and the equation apparently becomes elliptic. The system is in fact still hyperbolic (see Sec. IV.2), and the change in signature just reflects the fact that for $`f>1`$ the characteristic cones can tilt beyond the time axis.

The foliation equation (17) will prove to be very important when we study the formation of gauge shocks.

## III Focusing singularities

A very important property of slicing conditions is that of “singularity avoidance”. Singularity avoidance refers to the property of certain slicing conditions of slowing down coordinate time, by making the lapse go to zero, when the spatial volume elements $`\sqrt{\gamma }`$ go to zero (this is known as the “collapse of the lapse”). Recent advances in black hole excision techniques Thornburg87 ; Thornburg93 ; Seidel92a ; Anninos94e  would seem to minimize the need for singularity avoidance in the choice of slicing conditions. One should remember, however, that singularity avoidance is not only important when one is interested in studying black hole spacetimes where real physical singularities are present, but is also needed in order to prevent the formation of coordinate singularities caused by the focusing of the normal lines in regions with strong gravitational fields.

Bona et. al. have shown Bona97a  that the slicing condition (1) can avoid so-called “focusing singularities” for some choices of the function $`f(\alpha )`$. Here I will extend their analysis and show explicitly what type of behavior $`f(\alpha )`$ must have as $`\alpha `$ approaches zero in order to avoid such singularities.

Following Bona97a , we define a focusing singularity as a place where the spatial volume elements $`\gamma ^{1/2}`$ vanish at a bounded rate. Let us assume that such a singularity occurs after a finite proper time $`\tau _s`$ away from our initial time slice (as measured by the normal observers). From the definition of the lapse we see that the elapsed coordinate time will then be

$$\mathrm{\Delta }t=_0^{\tau _s}\frac{d\tau }{\alpha }.$$
(23)
We will further characterize the singularity by the rate at which $`\gamma ^{1/2}`$ approaches zero as a function of proper time. We will say a singularity is of order $`m`$ if $`\gamma ^{1/2}`$ approaches zero as

$$\gamma ^{1/2}\left(\tau _s\tau \right)^m,$$
(24)
with $`m`$ some constant power. Notice that $`m`$ must be strictly positive for there to be a singularity at all, and it must be larger than or equal to $`1`$ for the singularity to be approached at a bounded rate.

As the volume elements $`\gamma ^{1/2}`$ approach zero, there are clearly three possible behaviors for the lapse: 1) $`\alpha `$ remains finite as $`\gamma ^{1/2}`$ vanishes, 2) $`\alpha `$ vanishes as $`\gamma ^{1/2}`$ vanishes, and 3) $`\alpha `$ vanishes before $`\gamma ^{1/2}`$ vanishes. Case 1 would clearly imply that coordinate time remains finite at the singularity, so the singularity would not be avoided. However, from equation (16) one can easily see that if the lapse remains always finite it is impossible for the volume elements to ever vanish (remember that $`f(\alpha )`$ is never allowed to be zero). We then conclude that case 1 can never happen, which implies that the Bona-Masso slicing condition (1) always causes the lapse to collapse when the volume elements approach zero, for any $`f(\alpha )>0`$. Case 3, on the other hand, implies that the time slices stop advancing a finite coordinate time before the singularity is reached (the time slices can in fact move back under certain conditions, see below). We will call such behavior “strong singularity avoidance”. Finally, case 2 corresponds to the case when the lapse becomes zero at the same time as the volume elements. Whether in such a case the singularity is reached after a finite or infinite coordinate time will depend on the speed at which $`\alpha `$ approaches zero at the singularity. We will say that a slicing is “marginally singularity avoiding” if the singularity is reached after an infinite coordinate time.

To study under which conditions we can have strong or marginal singularity avoidance we must say something about the form of the function $`f(\alpha )`$. From now on we will therefore assume that, as $`\alpha `$ approaches zero, the function $`f(\alpha )`$ behaves as

$$f(\alpha )=A\alpha ^n,$$
(25)
with both $`A`$ and $`n`$ constants and $`A>0`$. Such an assumption implies that

$`{\displaystyle \frac{d\alpha }{\alpha f(\alpha )}}`$ $`=`$ $`{\displaystyle \frac{1}{A}}{\displaystyle \frac{d\alpha }{\alpha ^{n+1}}}`$ (28)
$`=`$ $`\{\begin{array}{cc}\mathrm{ln}\alpha ^{1/A}\hfill & n=0\hfill \\ 1/(nA\alpha ^n)\hfill & n0\hfill \end{array}`$
Let us first consider the case $`n0`$. From equation (16) we now find that

$$\gamma ^{1/2}\mathrm{exp}\left(\frac{1}{nA\alpha ^n}\right).$$
(29)
As $`\alpha `$ approaches zero we have two separate cases depending on the sign of $`n`$:

$$\underset{\alpha 0}{lim}\gamma ^{1/2}=\{\begin{array}{cc}\text{finite}\hfill & n<0\hfill \\ 0\hfill & n>0\hfill \end{array}$$
(30)
Since for $`n<0`$ the volume elements remain finite as the lapse approaches zero we conclude that such a case corresponds to strong singularity avoidance. On the other hand, for $`n>0`$ both the lapse and the volume elements go to zero at the same time so we can at most have marginal singularity avoidance.

For the case $`n=0`$ we find, again using (16), that

$`\gamma ^{1/2}\alpha ^{1/A}`$
$`\alpha \gamma ^{A/2}.`$ (31)
It is then clear that in this case $`\alpha `$ and $`\gamma ^{1/2}`$ also vanish at the same time.

We now need to decide if the cases $`n>0`$ and $`n=0`$, for which the lapse and the volume elements become zero at the same time, reach the singularity in an infinite or a finite coordinate time. For this we need to study the behavior of $`\alpha `$ as a function of proper time $`\tau `$ as we approach the singularity. Starting from equation (23) for the elapsed coordinate time we find

$`\mathrm{\Delta }t`$ $`=`$ $`{\displaystyle _0^{\tau _s}}{\displaystyle \frac{d\tau }{\alpha }}`$ (32)
$`=`$ $`{\displaystyle _{\alpha _0}^0}{\displaystyle \frac{d\tau /d\alpha }{\alpha }}𝑑\alpha ,`$
where $`\alpha _0`$ is the initial lapse and where we are already assuming that we are interested in the case when $`\alpha `$ vanishes at $`\tau _s`$. Equation (32) implies that if $`d\tau /d\alpha `$ remains different from zero as we approach the singularity then $`\mathrm{\Delta }t`$ will diverge and we will have marginal singularity avoidance. On the other hand, if $`d\tau /d\alpha `$ vanishes at the singularity as $`\alpha ^p`$ (with $`p`$ some positive power) or faster, then the integral will converge and the singularity will be reached in a finite coordinate time.

To find the behavior of $`d\tau /d\alpha `$ as we approach the singularity, we notice that equation (24) implies

$$d\gamma ^{1/2}/d\tau m\left(\tau _s\tau \right)^{m1},$$
(33)
and

$$\frac{d\mathrm{ln}\gamma ^{1/2}}{d\tau }=\frac{m}{\left(\tau _s\tau \right)}.$$
(34)
From this, together with equation (15), we now find

$$\frac{d\alpha /d\tau }{\alpha f(\alpha )}=\frac{m}{\left(\tau _s\tau \right)},$$
(35)
which can be integrated to give

$$\tau =\tau _s\mathrm{exp}\left(\frac{1}{m}\frac{d\alpha }{\alpha f(\alpha )}\right).$$
(36)
If we now take $`f(\alpha )`$ given by (25) we finally find

$`\tau `$ $`=`$ $`\tau _s\mathrm{exp}\left({\displaystyle \frac{1}{mA}}{\displaystyle \frac{d\alpha }{m\alpha ^{n+1}}}\right)`$ (39)
$`=`$ $`\{\begin{array}{cc}\tau _s\alpha ^{1/mA}\hfill & n=0\hfill \\ \tau _s\mathrm{exp}\left[1/(nmA\alpha ^n)\right]\hfill & n>0\hfill \end{array}`$
The case $`n<0`$ is not of interest here since we already showed that for such a case the lapse will vanish before we reach the singularity.

Let us consider the case $`n>0`$ first. The derivative of $`\tau `$ with respect to $`\alpha `$ then turns out to be

$$\frac{d\tau }{d\alpha }=\frac{1}{mA\alpha ^{n+1}}\mathrm{exp}\left(\frac{1}{nmA\alpha ^n}\right),$$
(40)
from which it is easy to see that as $`\alpha `$ approaches zero, $`d\tau /d\alpha `$ also approaches zero faster than any power. As we have seen, this means that the singularity is reached in a finite coordinate time, so the case $`n>0`$ does not avoid the singularity (not even marginally).

For $`n=0`$ we have, on the other hand

$$\frac{d\tau }{d\alpha }=\frac{1}{mA}\alpha ^{1/mA1}.$$
(41)
We then see that

$$\underset{\alpha 0}{lim}\frac{d\tau }{d\alpha }=\{\begin{array}{cc}0& mA<1\hfill \\ 1& mA=1\hfill \\ \mathrm{}& mA>1\hfill \end{array}$$
(42)
The case $`mA<1`$ therefore reaches the singularity in a finite coordinate time, while the cases $`mA1`$ reach it in an infinite coordinate time and are therefore marginally singularity avoiding.

Our final result can be summarized as follows: If $`f(\alpha )`$ behaves as $`f=A\alpha ^n`$ for small $`\alpha `$ and we have a singularity of order $`m`$, then

1. For $`n<0`$ we have strong singularity avoidance.
2. For $`n=0`$ and $`mA1`$ we have marginal singularity avoidance.
3. For both $`n>0`$, and $`n=0`$ with $`mA<1`$, we do not have singularity avoidance, even though the lapse collapses to zero at the singularity.

In the particular case when we have a singularity of order $`m=1`$, then harmonic slicing ($`n=0`$, $`A=1`$) marks the boundary between avoiding and reaching the singularity.

As a final observation, let us consider again the case $`n<0`$, for which we have shown that we have strong singularity avoidance. Looking at our original slicing condition (1) we see that if $`n2`$ then, as the lapse approaches zero, one can not guarantee that $`_t\alpha `$ will also approach zero. The lapse can therefore easily become negative and the slices will not only avoid the singularity but can in fact back away from it. This type of behavior is probably not desirable, as one runs the risk of having the time slices stop being space-like (they advance in one region and move back in another). If we want to guarantee that we have strong singularity avoidance without the lapse becoming negative we must limit ourselves to the region $`2<n<0`$. Notice that the 1+log family corresponds to $`n=1`$, and is precisely in the middle of this range, which probably accounts for the fact that empirically it has been found to be a very good choice.

## IV Gauge shocks

In physics, one talks about “shock waves” as solutions to the hydrodynamic equations where very sharp density gradients propagate through a medium at speeds that are higher than the speed of sound in that medium. Mathematically, shocks are discontinuous solutions of a non-linear hyperbolic system of equations characterized by the fact that characteristic lines converge at the discontinuity. The discontinuity propagates with a speed called the “shock speed” that is somewhere in between the values of the characteristic speeds in the regions behind and in front of the shock. Usually, in order to completely determine the form and speed of a shock one needs to supplement the evolution equations with extra conditions coming from physical considerations known as “entropy-conditions”(see for example Lax73 ). Such entropy conditions are necessary because once a discontinuity forms, the possible mathematical extensions through it are no longer unique and one needs a mechanism to choose the physically allowed solutions.

It is well known that physical shocks (i.e. shocks in the geometry) do not appear in general relativity. Solutions of the Einstein equations normally called “shock fronts” refer to discontinuities in the curvature of spacetime present in the initial data that propagate with the speed of light. These type of solutions are not shocks in a mathematical sense but are called instead “contact discontinuities”. Here we will not consider discontinuities in the geometry, but rather solutions to our hyperbolic slicing conditions that start from smooth initial data and develop discontinuities later, when the characteristic lines associated with the gauge cross. Since this is the defining property of a shock, we will call those solutions “gauge shocks”. However, we must stress the fact that, in contrast to the case of hydrodynamics, once a gauge shock forms one should make no attempt to continue the solution any further since such a shock will indicate that our coordinate system has broken down, and there is no physical principle that can be used to extend the solution beyond this point. We will therefore not have shock waves as such (i.e. propagating shocks), but just shock formation.

Gauge shocks were first studied in References Alcubierre97a ; Alcubierre97b , where it was found that discontinuities in the lapse can easily develop starting from smooth initial data in a wide variety of cases. It was also shown how, in some particular cases, one can even predict the exact time when a gauge shock would form by just analyzing the initial data.

In Alcubierre97a  a particular condition on the function $`f(\alpha )`$ was found for which gauge shocks would not form:

$$1f\alpha f^{}/2=0.$$
(43)
This condition, however, was derived using a hyperbolic formulation of the 3+1 evolution equations (the Bona-Masso formulation). In the following sections I will re-derive the condition making no reference to the Einstein equations in any form.

### IV.1 Linear degeneracy and shocks

Consider a system of equations of the form

$$_tu_i+_xF_i=q_ii\{1,\mathrm{},N_u\},$$
(44)
where $`F_i`$ and $`q_i`$ are arbitrary, possibly non-linear, functions of the $`u`$’s but not their derivatives. Notice that the system above can be written also as

$$_tu_i+\underset{j}{}M_{ij}_xu_j=q_ii\{1,\mathrm{},N_u\},$$
(45)
with $`M_{ij}=F_i/u_j`$ the Jacobian matrix.

Let $`\lambda _i`$ be the eigenvalues of the Jacobian matrix $`M`$. The system of equations is called “hyperbolic” if all the $`\lambda _i`$ are real. Further, the system is said to be “strongly hyperbolic” if the matrix $`M`$ has a complete set of eigenvectors $`𝐞_i`$ . Let us assume that this is the case, we then define the eigenfields $`w_i`$ in the following way

$$𝐮=R𝐰𝐰=R^1𝐮,$$
(46)
where $`R`$ is the matrix of column eigenvectors $`𝐞_i`$. One can show that the matrix $`R`$ is such that

$$RMR^1=\mathrm{\Lambda },$$
(47)
with $`\mathrm{\Lambda }=\mathrm{diag}(\lambda _i)`$. The evolution equation for the eigenfields $`w_i`$ then turns out to be

$$_tw_i+\lambda _i_xw_i=q_i^{},$$
(48)
with $`q_i^{}`$ a function of the $`w^{}s`$ but not their derivatives. In terms of the eigenfields the system transforms into a series of coupled advection equations with characteristic speeds given by the eigenvalues $`\lambda _i`$.

If a given eigenvalue $`\lambda _i`$ is independent of its corresponding eigenfield $`w_i`$, then we say that the eigenfield is “linearly degenerate” Lax73 ; Leveque92 , that is

$$\frac{\lambda _i}{w_i}=\underset{j=1}{\overset{N_u}{}}\frac{\lambda _i}{u_j}\frac{u_j}{w_i}=_u\lambda _i𝐞_i=\mathrm{\hspace{0.17em}0}.$$
(49)
Linear degeneracy is a sufficient condition for there not to be shocks associated with a given eigenfield. One can understand this intuitively by noticing that linear degeneracy implies that the characteristic lines do not change in response to changes in the field propagating along them.

### IV.2 Avoiding gauge shocks

In order to study the effects of our slicing condition without having to worry about the evolution of the spacetime itself we will now assume that we have a known background spacetime with metric $`g_{\mu \nu }`$. In that spacetime, we will consider some initial spatial slice, and then construct a foliation according to our slicing condition. Let $`T`$ be the time function associated with our foliation (that is, each spatial hypersurface will correspond to $`T=\mathrm{constant}`$). In Sec. II.2 above we saw that for the Bona-Masso family of slicing conditions, the function $`T`$ will satisfy the following foliation equation

$$\left[g^{\mu \nu }+\left(1\frac{1}{f(\alpha )}\right)n^\mu n^\nu \right]_\mu _\nu T=0,$$
(50)
where $`n^\mu `$ is the unit normal vector to the hypersurfaces and where $`_\mu `$ denotes covariant differentiation with respect to the 4-metric $`g_{\mu \nu }`$. The unit normal vector $`n^\mu `$ can be easily constructed from the time function $`T`$ in the following way

$$n_\mu =\frac{_\mu T}{\left(_\nu T^\nu T\right)^{1/2}},$$
(51)
where the overall minus sign is there to guarantee that we have a future pointing normal vector.

Let us now calculate the increment in $`T`$ if we move a proper distance $`d\tau `$ along the normal direction

$$dT=\left(d\tau n^\mu \right)_\mu T=d\tau \left(_\mu T^\mu T\right)^{1/2}.$$
(52)
On the other hand, from the definition of $`\alpha `$ we have $`d\tau =\alpha dT`$. Comparing both expressions we find

$$\alpha =\left(_\mu T^\mu T\right)^{1/2}.$$
(53)
The normal vector then takes the form:

$$n_\mu =\alpha _\mu T.$$
(54)
Consider now a particular point $`𝒫`$ in spacetime. To study the evolution of $`T`$ close to that point we construct locally flat coordinates $`(t,x^i)`$, so the metric close to $`𝒫`$ becomes the flat metric $`\eta _{\mu \nu }`$ and the Christoffel symbols vanish. Equation (50) then reduces to

$$\left[\eta ^{\mu \nu }an^\mu n^\nu \right]_\mu _\nu T=0,$$
(55)
where we have defined $`a:=1/f1`$. Expanding this equation we find

$`\left(1+a(n^0)^2\right)_t^2T+\left(\delta ^{ij}an^in^j\right)_i_jT`$
$`2an^0n^i_i_tT=0.`$ (56)
Let us now define $`\mathrm{\Pi }:=_tT`$ and $`\mathrm{\Psi }_i:=_iT`$. Equation (56) can then be transformed into the system

$`_t\mathrm{\Pi }`$ $`=`$ $`{\displaystyle \underset{i}{}}{\displaystyle \frac{2an^0n^i}{1+a(n^0)^2}}_i\mathrm{\Pi }`$ (57)
$`+{\displaystyle \underset{ij}{}}{\displaystyle \frac{\delta ^{ij}an^in^j}{1+a(n^0)^2}}_i\mathrm{\Psi }_j,`$
$`_t\mathrm{\Psi }_i`$ $`=`$ $`_i\mathrm{\Pi }.`$ (58)
In our locally flat coordinate system, the contravariant components of the unit normal vector become

$$n^0=+\alpha \mathrm{\Pi },n^i=\alpha \mathrm{\Psi }_i,$$
(59)
and the lapse reduces to

$$\alpha =\left(\mathrm{\Pi }^2𝚿^2\right)^{1/2},$$
(60)
with $`𝚿^2=_i\mathrm{\Psi }_i^2`$. The system (57)-(58) then takes the following form

$`_t\mathrm{\Pi }`$ $`=`$ $`{\displaystyle \underset{i}{}}{\displaystyle \frac{2a\alpha ^2\mathrm{\Pi }\mathrm{\Psi }_i}{1+a\alpha ^2\mathrm{\Pi }^2}}_i\mathrm{\Pi }`$ (61)
$`+`$ $`{\displaystyle \underset{ij}{}}{\displaystyle \frac{\delta _{ij}a\alpha ^2\mathrm{\Psi }_i\mathrm{\Psi }_j}{1+a\alpha ^2\mathrm{\Pi }^2}}_i\mathrm{\Psi }_j,`$
$`_t\mathrm{\Psi }_i`$ $`=`$ $`_i\mathrm{\Pi }.`$ (62)
In order to see if our system of equations is hyperbolic we consider derivatives along a fixed spatial direction, say $`x`$, and neglect derivatives along different directions. It is evident that in this case the variables $`\mathrm{\Psi }_q`$, with $`qx`$, can be considered as fixed and we only need to analyze the sub-system $`\{\mathrm{\Pi },\mathrm{\Psi }_x\}`$. The Jacobian matrix $`M`$ for this reduced system becomes:

$$M=\left(\begin{array}{cc}\frac{2a\alpha ^2\mathrm{\Pi }\mathrm{\Psi }_x}{1+a\alpha ^2\mathrm{\Pi }^2}& \frac{1a\alpha ^2\mathrm{\Psi }_x^2}{1+a\alpha ^2\mathrm{\Pi }^2}\\ 1& 0\end{array}\right).$$
(63)
The eigenvalues of this matrix are easily found to be

$`\lambda _\pm `$ $`=`$ $`{\displaystyle \frac{1}{1+a\alpha ^2\mathrm{\Pi }^2}}\{a\alpha ^2\mathrm{\Pi }\mathrm{\Psi }_x`$ (64)
$`\pm [1+a\alpha ^2(\mathrm{\Pi }^2\mathrm{\Psi }_x^2)]^{1/2}\},`$
with corresponding eigenvectors

$$𝐞_\pm =(\lambda _\pm ,1).$$
(65)
Having found the eigenvalues and eigenvectors we can now ask if our system of equations is hyperbolic. Consider first the case $`f=0`$. In such a case one can show that the eigenvalues reduce to

$$\lambda _\pm =\mathrm{\Psi }_x/\mathrm{\Pi },$$
(66)
that is, both eigenvalues are equal and the eigenvectors (65) do not form a complete set, so the system is only weakly hyperbolic.

In the case when $`f<0`$, one can use the fact that $`\mathrm{\Pi }^2>𝚿^2\mathrm{\Psi }_x^2`$ (which will hold for spacelike slices) to show that the term inside the square root in equation (64) is always negative and the system is not hyperbolic.

The $`f>0`$ case turns out to be more difficult to analyze. If $`0<f1`$, then it is not difficult to prove that the term inside the square root is always positive. This means that we have two distinct real eigenvalues and a complete set of eigenvectors, and the system is therefore strongly hyperbolic. If $`f>1`$, on the other hand, the term inside the square root can become negative for sufficiently large $`\mathrm{\Psi }_q^2:=_{ix}\mathrm{\Psi }_i^2`$. It would then appear that in such a case we do not have a hyperbolic system. However, the fact that hyperbolicity depends on the size of $`\mathrm{\Psi }_q^2`$ indicates that this is only a coordinate problem. Indeed, if we reorient our spatial coordinates in a way such that $`n^\mu `$ only has components along the $`(t,x)`$ directions, then $`\mathrm{\Psi }_q^2`$ vanishes and the eigenvalues become real. The fact that for a different orientation of the spatial coordinates we can have complex eigenvalues is just an indication that for $`f>1`$ the characteristic cones can tilt beyond the $`(t,x)`$ coordinate plane. This is analogous to solving the hydrodynamic equations using a supersonic reference frame: the hydrodynamic equations become elliptic, but this is just a consequence of choosing a bad coordinate system. We then conclude that for $`f>0`$, an orientation of the coordinates always exists such that our system of equations is strongly hyperbolic.

Using the expression for the eigenvectors above, the condition for linear degeneracy, Eq. (49), takes the form

$$C_\pm :=\lambda _\pm \frac{\lambda _\pm }{\mathrm{\Pi }}+\frac{\lambda _\pm }{\mathrm{\Psi }_x}=0.$$
(67)
A straightforward calculation gives the following independent linear combinations of the previous conditions

$`C_++C_{}=0=\alpha ^2\mathrm{\Pi }\left[2a\left(1+a\right)+\alpha a^{}\right]`$
$`\left\{{\displaystyle \frac{\alpha ^2\mathrm{\Pi }^2\left(1+a+a\alpha ^2\mathrm{\Psi }_q^2\right)+3\alpha ^2\left(\mathrm{\Pi }^2\mathrm{\Psi }_q^2\right)3}{\left(1+a\alpha ^2\mathrm{\Pi }^2\right)^3}}\right\},`$ (68)
$`C_+C_{}=0=\alpha ^2\mathrm{\Psi }_x^2\left[2a\left(1+a\right)+\alpha a^{}\right]`$
$`\left\{{\displaystyle \frac{3\alpha ^2\mathrm{\Pi }^2\left(1+a+a\alpha ^2\mathrm{\Psi }_q^2\right)+\alpha ^2\left(\mathrm{\Pi }^2\mathrm{\Psi }_q^2\right)1}{\left(1+a\alpha ^2\mathrm{\Pi }^2\right)^3\sqrt{1+a\alpha ^2\mathrm{\Psi }_q^2}}}\right\},`$ (69)
where $`a^{}:=da/d\alpha `$. For arbitrary $`\alpha `$, $`\mathrm{\Pi }`$ and $`\mathrm{\Psi }_i`$, the only way in which both these conditions can hold is if we take

$$2a\left(1+a\right)+\alpha a^{}=0,$$
(70)
which turns out to be equivalent to

$$1f\alpha f^{}/2=0.$$
(71)
This is precisely the condition (43) found in reference Alcubierre97a . The main difference between the derivation above and the one of reference Alcubierre97a  is that in that reference the condition was derived using a hyperbolic formulation of the Einstein equations (the Bona-Masso formulation), while the new derivation is based purely on analyzing the slicing condition and makes no reference to the Einstein equations in any form, showing the generality of the condition.

As already shown in reference Alcubierre97a , condition (43) can be trivially solved to find

$$f=1+k/\alpha ^2,$$
(72)
with $`k>0`$ an arbitrary constant. For $`k=0`$ we recover harmonic slicing. On the other hand, if $`k0`$ we see that for small $`\alpha `$ we have $`f\alpha ^2`$, so the results of Sec. III imply that the slicing will be strongly singularity avoiding. However, we can also see that in this case we are in the regime for which the lapse can easily become negative. This means that the solution (72) has a serious drawback for any non-zero value of $`k`$, since it can allow the lapse to become negative as it collapses toward zero.

In the next sections we will see how one can still obtain useful slicing conditions by looking for approximate solutions to condition (43).

### IV.3 Zero order shock avoidance

In the previous section we found that if we want to guarantee that no shocks will form, then we must choose the function $`f(\alpha )`$ in a way that is incompatible with having a lapse that does not become negative when it collapses toward zero (except in the specific case of harmonic slicing). Here we will relax our requirements and look for approximate solutions of condition (43).

We start by assuming that the lapse is very close to $`1`$, that is

$$\alpha =1+ϵ,$$
(73)
with $`ϵ1`$. Notice that the limit above applies to situations that are close to flat space, but generally does not apply to strong field regions (like the region close to a black hole) where the lapse can be expected to be considerably smaller than 1. However, in such regions considerations about singularity avoidance are probably more important. Our aim is to find slicing conditions that can avoid singularities in strong field regions, and at the same time don’t have a tendency to generate shocks in weak field regions.

We can now expand $`f`$ in terms of $`ϵ`$ as

$`f`$ $`=`$ $`a_0+a_1ϵ+𝒪(ϵ^2)`$ (74)
$`=`$ $`a_0+a_1(\alpha 1)+𝒪(ϵ^2),`$
and look for solutions to (43) to lowest order in $`ϵ`$.

Substituting (74) into (43) we find

$$1a_0a_1/2+𝒪(ϵ)=0.$$
(75)
This means that if we want condition (43) to be satisfied to zero order in $`ϵ`$ we must have

$$a_1=2\left(1a_0\right),$$
(76)
which implies

$`f`$ $`=`$ $`a_0+2\left(1a_0\right)ϵ+𝒪(ϵ^2)`$ (77)
$`=`$ $`\left(3a_02\right)+2\left(1a_0\right)\alpha +𝒪(ϵ^2).`$
We must remember that (77) is just an expansion for small $`ϵ`$. Any form of the function $`f(\alpha )`$ that has the same expansion to first order in $`ϵ`$ will also satisfy condition (43) to zero order. One family of such functions emerges if we ask for $`f(\alpha )`$ to have the following form

$$f=\frac{p_0}{1+q_1ϵ}.$$
(78)
It is not difficult to show that for this to have an expansion of the form (77) we must ask for

$$p_0=a_0,q_1=2\left(a_01\right)/a_0,$$
(79)
which implies

$`f`$ $`=`$ $`{\displaystyle \frac{a_0^2}{a_0+2\left(a_01\right)ϵ}}`$ (80)
$`=`$ $`{\displaystyle \frac{a_0^2}{\left(2a_0\right)+2\left(a_01\right)\alpha }}.`$
Notice that if we take $`a_0=1`$, we recover harmonic slicing. But there is one other case that is of special interest: For $`a_0=2`$ the previous solution reduces to $`f=2/\alpha `$, which corresponds to a member of the 1+log family. The crucial observation here is that, as already mentioned in Sec. II.1, this specific member of the 1+log family is precisely the one that has been found empirically to be very robust in black hole simulations Arbona99 ; Alcubierre01a ; Alcubierre02a . The fact that it is the only member of the 1+log family that satisfies condition (43) even approximately means that one should in fact expect it to be particularly well behaved.

### IV.4 First order shock avoidance

We can now go one order higher in $`ϵ`$ to obtain other interesting forms of $`f`$. Taking now

$`f`$ $`=`$ $`a_0+a_1ϵ+a_2ϵ^2+𝒪(ϵ^3)`$ (81)
$`=`$ $`a_0+a_1(\alpha 1)+a_2(\alpha 1)^2+𝒪(ϵ^3),`$
we find, after substituting into condition (43), that

$$\left(1a_0a_1/2\right)\left(3a_1/2+a_2\right)ϵ+𝒪(ϵ^2)=0.$$
(82)
Condition (43) will be satisfied to first order if we take

$`a_1`$ $`=`$ $`2\left(1a_0\right),`$ (83)
$`a_2`$ $`=`$ $`3a_1/2=3\left(1a_0\right).`$ (84)
So the expansion of $`f`$ takes the form

$`f`$ $`=`$ $`a_0+2\left(1a_0\right)ϵ3\left(1a_0\right)ϵ^2+𝒪(ϵ^3)`$ (85)
$`=`$ $`\left(6a_05\right)+8\left(1a_0\right)\alpha `$
$`3\left(1a_0\right)\alpha ^2+𝒪(ϵ^3).`$
Just as before, we can now look for rational functions that have the above expansion. One such possibility is to ask for $`f`$ to have the form

$$f=\frac{p_0}{1+q_1ϵ+q_2ϵ^2}.$$
(86)
In order for $`f`$ to have the expansion (85) we must take

$`p_0`$ $`=`$ $`a_0,`$ (87)
$`q_1`$ $`=`$ $`2\left(a_01\right)/a_0,`$ (88)
$`q_2`$ $`=`$ $`\left(1a_0\right)\left[3a_0+\left(1a_0\right)\right]/a_0^2,`$ (89)
which means that $`f`$ takes the final form

$`f`$ $`=`$ $`{\displaystyle \frac{a_0^3}{a_0^22a_0\left(1a_0\right)ϵ+\left(1a_0\right)\left(4a_0\right)ϵ^2}}`$ (90)
$`=`$ $`{\displaystyle \frac{a_0^3}{\left(43a_0\right)+\alpha \left(1a_0\right)\left[\left(4a_0\right)\alpha 8\right]}}.`$
If we take $`a_0=1`$ we again recover harmonic slicing. Another interesting case is obtained by asking for

$$43a_0=0a_0=4/3,$$
(91)
since in that case we will have $`f\alpha ^1`$ for small $`\alpha `$, and as we have seen this implies good singularity avoidance. For such a choice the function $`f`$ reduces to

$$f=\frac{8}{3\alpha \left(3\alpha \right)}.$$
(92)
For small $`\alpha `$, this form of $`f`$ behaves as a member of the 1+log family. Moreover, it satisfies condition (43) to higher order than the usual choice $`f=2/\alpha `$. One could worry about the fact that for the choice (92) the function $`f`$ can become negative for $`\alpha >3`$. However, such a situation is unlikely to arise in practice since the initial lapse is always taken to be at most 1 throughout the region of interest, and later evolution usually makes it even smaller. Because of these facts, the slicing condition given by the choice (92) would seem to be a good candidate for a robust slicing condition when evolving systems with strong gravitational fields (including black holes). Whether this is true or not can only be settled by numerical experiments where this form of $`f`$ is tested against more traditional choices.

Notice also that, whereas in the case when $`f=2/\alpha `$ the asymptotic gauge speed in regions where $`\alpha 1`$ is $`\sqrt{2}1.41`$, in the case (92) the gauge speed in those regions is only $`\sqrt{4/3}1.15`$, which is much closer to the physical speed of light and might represent an extra advantage as gauge effects will not propagate much faster than physical effects.

## V Conclusion

I have considered the Bona-Masso hyperbolic family of slicing conditions and have studied under which circumstances such hyperbolic slicings can avoid two different types of singularities: focusing singularities, defined as those for which the spatial volume elements vanish at a bounded rate, and gauge shocks, defined as coordinate singularities for which the lapse becomes discontinuous as a consequence of the crossing of the characteristic lines associated with the propagation of the gauge.

In the case of focusing singularities, I have extended the analysis of Bona et. al. and shown that, depending on the form that the function $`f(\alpha )`$ defining the slicing takes in the limit of small $`\alpha `$, one can have three different types of behavior: the lapse vanishes before the spatial volume elements do (collapse of the lapse and strong singularity avoidance), the lapse vanishes at the same time as the spatial volume elements and the singularity is reached after an infinite coordinate time (collapse of the lapse and marginal singularity avoidance), or the lapse vanishes at the same time as the volume elements but the singularity is still reached in a finite coordinate time (collapse of the lapse but no singularity avoidance). Harmonic slicing falls into the marginal singularity avoiding case, whereas the more commonly used 1+log family falls into the strong singularity avoiding case.

For the case of gauge shocks I have re-derived, in a way that is independent of the Einstein equations, a condition on the function $`f(\alpha )`$ found previously that avoids them. This condition, unfortunately, is severely restrictive and implies that the lapse can easily become negative during evolutions. I have therefore studied different forms of the function $`f(\alpha )`$ that satisfy the condition only approximately. This study has shown that one specific member of the 1+log family that has previously been found empirically to be particularly robust, is in fact the only member of that family that satisfies the condition for shock avoidance even to lowest order. By asking for the shock avoidance condition to be satisfied to higher order, I have found a new form of the function $`f(\alpha )`$ that has the potential of being even more robust than the 1+log slicings for simulations of strongly gravitating systems.

As a final comment, it is important to mention that elliptic slicing conditions such as maximal slicing can easily avoid both focusing singularities (since the volume elements are not allowed to change) and gauge shocks (since the speed of propagation is infinite), so they should in principle be more robust that the slicings considered here. Elliptic conditions, however, are considerably more computationally expensive. Not only that, but they are typically much more restrictive. For example, maximal slices might not always exist (they typically do not exist in cosmological scenarios). Finally, elliptic equations require boundary conditions that might be difficult to impose in some situations, in particular when one has internal boundaries such as those associated with black hole excision. It is because of these reasons that one should consider alternatives to elliptic gauge conditions, such as those studied here.

###### Acknowledgements.

I wish to thank Bernd Brügmann, Dario Nuñez, Ed Seidel and Deirdre Shoemaker for many useful comments. I also thank the “Caltech Visitors Program for the Numerical Simulation of Gravitational Wave Sources” for their hospitality when completing the final version of this manuscript. This work was supported by CONACyT Mexico under the repatriation program.