### 4.3 The BSSN formulation

The BSSN formulation is based on the 3+1 decomposition of Einstein’s field equations. Unlike the harmonic formulation, which has been motivated by the mathematical structure of the equations and the understanding of the Cauchy formulation in general relativity, this system has been mainly developed and improved based on its capability of numerically evolving spacetimes containing compact objects in a stable way. Interestingly, it turns out that in spite of the fact that the BSSN formulation is based on an entirely different motivation, mathematical questions like the well-posedness of its Cauchy problem can be answered, at least for most gauge conditions.

In the BSSN formulation, the three metric and the extrinsic curvature are decomposed according to

Here, and are the trace and the trace-less part, respectively, of the conformally-rescaled extrinsic curvature. The conformal factor is determined by the requirement for the conformal metric to have unit determinant. Aside from these variables one also evolves the lapse (), the shift () and its time derivative (), and the variable
In terms of the operator the BSSN evolution equations are
Here, quantities with a tilde refer to the conformal three metric , which is also used in order to raise and lower indices. In particular, and denote the covariant derivative and the Christoffel symbols, respectively, with respect to . Expressions with a superscript refer to their trace-less part with respect to the conformal metric. Next, the sum represents the Ricci tensor associated with the physical three metric , where
The term in Eq. (4.55) is set equal to the right-hand side of Eq. (4.59). The parameter in the latter equation modifies the evolution flow off the constraint surface by adding the momentum constraint to the evolution equation for the variable . This parameter was first introduced in [10] in order to compare the stability properties of the BSSN evolution equations with those of the ADM formulation.

The gauge conditions, which are imposed on the lapse and shift in Eqs. (4.52, 4.54, 4.55), were introduced in [52] and generalize the Bona–Massó condition [62] and the hyperbolic Gamma driver condition [11]. It is assumed that the functions , and are strictly positive and smooth in their arguments, and that and are smooth functions of their arguments. The choice

with a positive constant, corresponds to the evolution system used in many black-hole simulations based on slicing and the moving puncture technique (see, for instance, [423] and references therein). Finally, the source terms , and are defined in the following way: denoting by and the Ricci tensors belonging to the three-metric and the spacetime metric, respectively, and introducing the constraint variables
the source terms are defined as
For vacuum evolutions one sets , and . When matter fields are present, the Einstein field equations are equivalent to the evolution equations (4.524.59) setting , , and the constraints , and .

When comparing Cauchy evolutions in different spatial coordinates, it is very convenient to reformulate the BSSN system such that it is covariant with respect to spatial coordinate transformations. This is indeed possible; see [77, 82]. One way of achieving this is to fix a smooth background three-metric , similarly as in Section 4.1, and to replace the fields and by the scalar and vector fields

where and denote the determinants of and , and is the covariant derivative associated to the latter. If is flat and time-independent, the corresponding BSSN equations are obtained by replacing and in Eqs. (4.524.59, 4.60, 4.61, 4.644.66).

#### 4.3.1 The hyperbolicity of the BSSN evolution equations

In fact, the ADM formulation in the spatial harmonic gauge described in Section 4.2.3 and the BSSN formulation are based on some common ideas. In the covariant reformulation of BSSN just mentioned, the variable is just the quantity defined in Eq. (4.40), where is replaced by the conformal metric . Instead of requiring to vanish, which would convert the operator on the right-hand side of Eq. (4.60) into a quasilinear elliptic operator, one promotes this quantity to an independent field satisfying the evolution equation (4.59) (see also the discussion below Equation (2.18) in [390]). In this way, the -block of the evolution equations forms a wave system. However, this system is coupled through its principal terms to the evolution equations of the remaining variables, and so one needs to analyze the complete system. As follows from the discussion below, it is crucial to add the momentum constraint to Eq. (4.59) with an appropriate factor in order to obtain a hyperbolic system.

The hyperbolicity of the BSSN evolution equations was first analyzed in a systematic way in [373], where it was established that for fixed shift and densitized lapse,

the evolution system (4.53, 4.564.59) is strongly hyperbolic for and and symmetric hyperbolic for and . This was shown by introducing new variables and enlarging the system to a strongly or symmetric hyperbolic first-order one. In fact, similar first-order reductions were already obtained in [196, 188]. However, in [373] it was shown that the first-order enlargements are equivalent to the original system if the extra constraints associated to the definition of the new variables are satisfied, and that these extra constraints propagate independently of the BSSN constraints , and . This establishes the well-posedness of the Cauchy problem for the system (4.69, 4.53, 4.564.59) under the aforementioned conditions on and . Based on the same method, a symmetric hyperbolic first-order enlargement of the evolution equations (4.52, 4.53, 4.564.59) and fixed shift was obtained in [52] under the conditions and and used to construct boundary conditions for BSSN. First-order strongly-hyperbolic reductions for the full system (4.524.59) have also been recently analyzed in [82].

An alternative and efficient method for analyzing the system consists in reducing it to a first-order pseudodifferential system, as described in Section 3.1.5. This method has been applied in [308] to derive a strongly hyperbolic system very similar to BSSN with fixed, densitized lapse and fixed shift. This system is then shown to yield a well-posed Cauchy problem. In [52] the same method was applied to the evolution system (4.524.59). Linearizing and localizing, one obtains a first-order system of the form . The eigenvalues of are , , , , , , , where we have defined and . The system is weakly hyperbolic provided that

and it is strongly hyperbolic if, in addition, the parameter and the functions , , and can be chosen such that the functions
are bounded and smooth. In particular, this requires that the nominators converge to zero at least as fast as the denominators when , or , respectively. Since , the boundedness of requires that . For the standard choice , the conditions on the gauge parameters leading to strong hyperbolicity are, therefore, , and . Unfortunately, for the choice (4.62, 4.63) used in binary black-hole simulations these conditions reduce to
which is typically violated at some two-surface, since asymptotically, and while near black holes is small and positive. It is currently not known whether or not the Cauchy problem is well posed if the system is strongly hyperbolic everywhere except at points belonging to a set of zero measure, such as a two-surface. Although numerical simulations based on finite-difference discretizations with the standard choice (4.62, 4.63) show no apparent sign of instabilities near such surfaces, the well-posedness for the Cauchy problem for the BSSN system (4.524.59) with the choice (4.62, 4.63) for the gauge source functions remains an open problem when the condition (4.72) is violated. However, a well-posed problem could be formulated by modifying the choice for the functions and such that and are guaranteed to hold everywhere.

Yet a different approach to analyzing the hyperbolicity of BSSN has been given in [219, 220] based on a new definition of strongly and symmetric hyperbolicity for evolution systems, which are first order in time and second order in space. Based on this definition, it has been verified that the BSSN system (4.69, 4.53, 4.564.59) is strongly hyperbolic for and and symmetric hyperbolic for . (Note that this generalizes the original result in [373] where, in addition, was required.) The results in [220] also discuss more general 3+1 formulations, including the one in [308] and construct constraint-preserving boundary conditions. The relation between the different approaches to analyzing hyperbolicity of evolution systems, which are first order in time and second order in space, has been analyzed in [221].

Strong hyperbolicity for different versions of the gauge evolution equations (4.52, 4.54, 4.55), where the normal operator is sometimes replaced by , has been analyzed in [222]. See Table I in that reference for a comparison between the different versions and the conditions they are subject to in order to satisfy strong hyperbolicity. It should be noted that when and is replaced by , additional conditions restricting the magnitude of the shift appear in addition to and .

#### 4.3.2 Constraint propagation

As mentioned above, the BSSN evolution equations (4.524.59) are only equivalent to Einstein’s field equation if the constraints

are satisfied. Using the twice contracted Bianchi identities in their 3+1 decomposed form, Eqs. (4.45, 4.46), and assuming that the stress-energy tensor is divergence free, it is not difficult to show that the equations (4.524.59) imply the following evolution system for the constraint fields [52, 220]:
This is the constraint propagation system for BSSN, which describes the propagation of constraint violations, which are usually present in numerical simulations due to truncation and roundoff errors. There are at least three reasons for establishing the well-posedness of its Cauchy problem. The first reason is to show that the unique solution of the system (4.744.76) with zero initial data is the trivial solution. This implies that it is sufficient to solve the constraints at the initial time . Then, any smooth enough solution of the BSSN evolution equations with such data satisfies the constraint propagation system with , and , and it follows from the uniqueness property of this system that the constraints must hold everywhere and at each time. In this way, one obtains a solution to Einstein’s field equations. However, in numerical calculations, the initial constraints are not exactly satisfied due to numerical errors. This brings us to the second reason for having a well-posed problem at the level of the constraint propagation system; namely, the continuous dependence on the initial data. Indeed, the initial constraint violations give rise to constraint violating solutions; but, if these violations are governed by a well-posed evolution system, the norm of the constraint violations is controlled by those of the initial violations for each fixed time . In particular, the constraint violations must converge to zero if the initial constraint violations do. Since the initial constraint errors go to zero when resolution is increased (provided a stable numerical scheme is used to solve the constraints), this guarantees convergence to a constraint-satisfying solution. Finally, the third reason for establishing well-posedness for the constraint propagation system is the construction of constraint-preserving boundary conditions, which will be explained in detail in Section 6.

The hyperbolicity of the constraint propagation system (4.744.76) has been analyzed in [220, 52, 81, 80], and [315] and shown to be reducible to a symmetric hyperbolic first-order system for . Furthermore, there are no superluminal characteristic fields if . Because of finite speed of propagation, this means that BSSN with (which includes the standard choice ) does not possess superluminal constraint-violating modes. This is an important property, for it shows that constraint violations that originate inside black hole regions (which usually dominate the constraint errors due to high gradients at the punctures or stuffing of the black-hole singularities in the turducken approach [156, 81, 80]) cannot propagate to the exterior region.

In [353] a general result is derived, showing that under a mild assumption on the form of the constraints, strong hyperbolicity of the main evolution system implies strong hyperbolicity of the constraint propagation system, with the characteristic speeds of the latter being a subset of those of the former. The result does not hold in general if “strong” is replaced by “symmetric”, since there are known examples for which the main evolution system is symmetric hyperbolic, while the constraint propagation system is only strongly hyperbolic [108].