1. Introduction

The widely used treatment of Einstein's equations in numerical relativity is to cast them to the form of a nonlinear hyperbolic system with constraints (e.g., [1–9]) and to solve by employing sophisticated discretization techniques. In the course of solution, the constraint part is either monitored, or explicitly imposed. It was observed that the solution of the evolution part with no constraints produces a violation which grows rapidly breaking computations in a short time [7, 10, 11]. An attempt to control constraint violation, by projecting the solution, or by incorporating constraint quantities in the evolution equations, results in a longer life time of calculations (e.g., [12–15]). It was found in [13] that exponentially growing constraint violating solutions converge to unstable solutions of the dynamic equations, which suggests that the constraint violation is closely related to loss of stability in the system.

An exact solution to evolution equations in the entire space has a property that it satisfies constraint equations automatically as long as it satisfies them initially. However, in numerical simulations, because of the roundoff and truncation errors, one cannot hope for automatic constraint compliance. Instead, care must be taken to ensure that the inserted perturbations are small, and remain small during the evolution.

The behavior of the solution can be improved significantly by introducing special sets of boundary data, the so-called constraint-preserving boundary conditions, or conditions that imply trivial evolution of constraints [13, 16]. Several sets of such data were proposed for various first-order formulations of Einstein's equations (e.g., [3–5, 9, 13, 16, 17, 18]). An approach not involving first-order reduction has been proposed in [19] where boundary conditions were constructed by projecting Einstein's equations on time-like boundaries. These conditions are typically written as a system of partial differential equations restricted to the boundary, and in cases when the equations are time-dependent and decouple from the bulk system, the equations may be integrated in time to produce regular Dirichlet data that is compatible with constraints [4, 18, 20].

In this work, two sets of well-posed homogeneous algebraic constraint-preserving boundary conditions for the linearized Baumgarte-Shapiro-Shibata-Nakamura formulation [2, 8] are constructed. As being common, our derivation starts from considering the evolution equations for constraint quantities and looks for sets of data for the variables of the main system that guarantee zero Dirichlet data for the constraint quantities. The procedure is similar to the procedure found in [4] but (a) does not employ reduction to first order, and (b) does not involve integration of equations in time along the boundary. Instead, following [17, 21], we rewrite the equations in a special form to find well-posed constraint-preserving boundary conditions by direct inspection. The approach can be generalized to produce boundary conditions of the evolving type (see, [4]) and the differential type [13, 16]. To further justify the proposed conditions, we derive an energy estimate for the nonlinear BSSN system with boundaries extending the results of [18, 22].

Recently, a significant progress has been achieved in establishing the well-posedness of the constraint-preserving boundary conditions for the generalized harmonic formulation of the Einstein equations [13, 16, 23]. Methods have been developed to study the well-posedness of the higher-order differential conditions [24, 25] and Sommerfeld-type constraint-preserving conditions [23, 26, 27]. This work is intended to prepare the ground to formulate similar results for the BSSN formulation.

In Section 2, we summarize the derivation of the BSSN formulation and consider the choices of lapse and shift most commonly found in numerical relativity. In Section 3, the linearized BSSN equations are introduced. In Section 4, the constraint-preserving boundary conditions are derived for the second order in time reduction of the linearized BSSN system. In Sections 5 and 6, we extend these results and define the initial-boundary value problem for the original linearized BSSN system. In Section 7, we obtain an energy estimate for the nonlinear BSSN system without invoking the second order in time reduction and discuss its applications to the validation of the proposed boundary conditions.

2. The Trace-Free Decomposition of the Arnowitt-Deser-Misner System

To point out some facts about the nature and properties of the BSSN formulation (see [2, 8] and also, some special cases in [14, 28]), let us briefly recall the derivation in the case of vacuum fields where the right-hand side of Einstein's equation is zero. A reader not interested in the BSSN derivation may proceed to Section 3 where the linearized system is given.

The derivation starts from the Arnowitt-Deser-Misner (ADM) decomposition [29, 30]: Here denotes the lapse, the are the components of the shift vector , and are the components of the spatial metric . The components of the 4-dimensional metric are given by Here denotes the matrix inverse to . Indices on all other quantities are raised and traces taken with respect to the spatial metric. Also, is the convective derivative and is the covariant derivative operator associated with the spatial metric. The extrinsic curvature is defined by (2.1). We assume that global Cartesian coordinates , , , are specified. Furthermore, are the components of the spatial Ricci tensor where are the spatial Christoffel symbols defined by .

The operator in (2.2) contains second-order spatial derivatives of unknown fields and is very difficult to analyze. As a result, it is onerous to judge properties of (2.2) and properties of , in general. However, a simplified equation can be derived for the trace of the extrinsic curvature . Contracting (2.2), and using (2.1) and (2.3), we find The simplicity of (2.7) suggests that the evolution of the trace of the extrinsic curvature be separated from the system. Specifically, we introduce the trace of the extrinsic curvature and the trace-free part of the extrinsic curvature as new variables. Then (2.7) yields Unless the lapse function is chosen with care, (2.8) is expected to be unstable. For example, for a spatially independent lapse and zero-shift vector, (2.8) yields an estimate which implies that , or that the solution is unbounded in a finite time. (This is a well-known example of a coordinate singularity.) A well-posed choice of the lapse function is given by the maximal slicing condition [28] with this condition, (2.8) reduces to .

Alternatively, we can use harmonic slicing [2, 28] (a particular case of the Bona-Massó family of -driving slicing conditions , [31, 32]), which corresponds to setting

The equation on is obtained from (2.2), (2.8), and (2.1) as To proceed with the derivation, we need a splitting for the spatial metric compatible with the splitting of into and . In the BSSN formulation, the desired splitting is achieved by introducing the conformal factor and the conformal metric , . Using the Leibnitz formula for differentiating the determinant of a matrix one finds that the derivative of the conformal metric is trace-free. By applying operator on the definition of and using (2.12) and (2.1), we get the second equation of our system Now using (2.13) and (2.1), we obtain the third equation where and are the conformal analogs of the variables and . Beginning with the last equation, indices are lowered and raised with the conformal metric and its inverse . In this case, , and it is easy to redefine .

The remaining two equations can be obtained from (2.11) which can be rewritten in terms of as The Ricci tensor in terms of the conformal metric becomes [2, 28] Here and are the covariant derivatives associated with the conformal metric. The first line in (2.17) can be rewritten to obtain This suggests that we introduce a new variable Substituting (2.18) in (2.16), we derive the fourth evolution equation where

The evolution equation for is obtained by applying the operator to (2.19) and using the momentum constraint Next we notice that , and thus (2.4) takes the form Solving this equation for and substituting the result in (2.22), we derive the fifth equation of the BSSN system where Equations (2.8), (2.14), (2.15), (2.20), and (2.24) constitute the core of the BSSN formulation. These equations are usually supplemented by one or more equations describing the choice of the lapse and shift functions. Thus the lapse and shift are not given a priori but dynamically determined from the metric and other quantities. In this work, we will assume the harmonic lapse condition (2.10). Further, we consider either a prescribed shift or a shift that follows from the gamma-freezing condition (cf. [28]).

3. Linearization around Minkowski Space

Minkowski spacetime in Cartesian coordinates is represented by the trivial solution to the ADM system: , , , . Consider perturbations of the ADM variables around the Minkowski spacetime: , , , and , with the , , , and being small. Substituting these expressions into the definitions of the BSSN variables and neglecting terms of second- and higher-order in perturbations we arrive at Substituting the linearized quantities in (2.14), (2.15), (2.8), (2.20), (2.24), and (2.10), and ignoring the terms which are second- and higher-order in , , , , and , we derive the linearization of the BSSN system Notice that the linearized harmonic lapse condition is included in this system in the form of (3.3).

Linearization of the Hamiltonian and the momentum constraint equations yields correspondingly Since it can be written both in terms of and , the Hamiltonian constraint appears in two different forms. Also, introducing the new variable entails an artificial constraint

The linearized problem then consists of determining , , , , , from (3.2)–(3.7), given the initial data and an admissible boundary data. The constraint equations (3.8)-(3.9) may or may not be imposed during the process of solution. The initial data , , , , can be determined from and using (3.1). It can be verified that if and satisfy the linearized Hamiltonian and momentum constraints of the ADM system, then , , , , satisfy the constraint equations (3.8)–(3.10). The linearized Hamiltonian and momentum constraints in the ADM system are and .

4. Constraint-Preserving Boundary Conditions

The BSSN system is a constrained evolution system in the following sense: (3.2)–(3.7) that are called the evolution equations contain both temporal and spatial derivatives; (3.8)–(3.10) which are called the constraint equations involve spatial derivatives only. It has been assumed for a long time that for properly selected boundary data, a solution to (3.2)–(3.7) satisfies the constraints automatically if it satisfies them initially. Examples of such data, however, were constructed only recently [18] for a first-order reduction of the BSSN system. In this paper, we propose a set of constraint-preserving boundary conditions for the BSSN system in second order.

First, we notice that for a solution of (3.2)–(3.7), the constraints (3.8) and (3.10) are consequences of (3.9), so we can focus on just the last one. In view of (3.2), (3.6), and (3.7), the time derivative of (3.8) becomes These equations state that both parts of (3.8) are satisfied as long as they are satisfied initially and (3.9) is true. Similarly, if (3.9) holds, then the time derivative of (3.10) is zero in view of (3.5) and (3.7). Thus (3.10) remains zero provided that it is zero initially.

We will now construct boundary conditions for the system (3.2)–(3.7) that preserve (3.9). We introduce the new variable Equation (3.9) is satisfied if and only if .

Now let us derive an equation for the evolution of . By differentiating (3.6) with respect to time and using (3.5) to replace in the result, we obtain Terms in , , , , cancel in view of (3.2)–(3.4), (3.7).

Similarly, by differentiating (3.4) and substituting (3.3) for , one derives Differentiating (4.2) twice with respect to time and using (4.3) and (4.4), we obtain The initial value can be determined from and using (4.2) and must be zero for the physical initial data. The initial values can be calculated by differentiating (4.2) in time and substituting (3.4) and (3.6) for and . Thus we have It can be verified by substitution that if , , and satisfy (3.8) and (3.10) then .

We must now select the boundary conditions on that imply the trivial evolution of (4.5). However, the boundary conditions on are expected not to be given freely but rather be determined by the boundary conditions (and the data) on the variables and . Similarly, the initial data and is determined by , , , , and . Our goal then is to construct the boundary conditions on the main variables and in such a way as to guarantee homogeneous data for . The standard approach is to study the relationship between the main and the constraint variables (cf. [3, 4, 9, 16–18, 21, 24]) by considering both the definition (4.2) and the evolution equations (4.3), (4.4).

We introduce the scalar products and for the spaces of vector fields and matrix fields on , respectively. The norms naturally associated with these scalar products are and . We introduce the energy function for the constraint quantity If we prove that the energy remains zero at all times, then, in view of the trivial initial data, we will have (cf. [33]).

Differentiating in time and using Green's first identity component-wise to transfer the spatial derivative in the second term and using (4.5), we obtain The energy is not increasing if The desired boundary conditions on and will follow immediately if we rewrite (4.9) in terms of the main variables.

We assume that the boundary is a combination of arbitrarily oriented planes and consider any of its faces. Let be the unit vector perpendicular to the face. Let and complement to form an orthonormal triple. For example, can be chosen to be any unit vector parallel to the boundary and be the cross product of and , so that . At a face of a flat boundary, the divergence of a vector field can be expressed in terms of the directional derivatives along vectors , , and as Similarly, the gradient of a scalar field becomes Next we note that at any point of the boundary, a symmetric trace free matrix is spanned by Introducing the scalar functions we rewrite as follows: Substituting (4.14) into (4.2) and using (4.10) and (4.11) to replace partial derivatives with directional derivatives, we obtain The last equation implies that Either of the following two sets of boundary conditions imply that on : Let us prove the constraint-preserving character of (4.17). For (4.17), imply that the second multiplier of the first term in (4.16) is zero. Thus the first term drops out. Now consider the first multipliers in the second and third terms of (4.16). By commuting partial derivatives, terms in , , , and drop out as a direct consequence of (4.17). Also (4.3) can be solved for the second normal derivatives of and in terms of the temporal and tangential derivatives using the following representation of the Laplace operator: It can be easily noticed that (4.17) implies that all tangential and temporal derivatives of and vanish at the face of the boundary.

The condition (4.18), in turn, eliminates the first multiplier of the first term and the second multipliers in the second and third terms.

More examples of constraint-preserving boundary conditions can be proposed by the inspection of (4.16). For example, is equivalent to the set of the following differential boundary conditions that can be used in a numerical method [13, 16]: In this case, we can prescribe the Dirichlet data for