Abstract and Applied Analysis

VolumeΒ 2008, Article IDΒ 742040, 21 pages

http://dx.doi.org/10.1155/2008/742040

## Constraint-Preserving Boundary Conditions for the Linearized Baumgarte-Shapiro-Shibata-Nakamura Formulation

Department of Mathematics, California State University Northridge, 18111 Nordhoff street, Northridge, CA 91330, USA

Received 3 August 2007; Accepted 19 February 2008

Academic Editor: NorimichiΒ Hirano

Copyright Β© 2008 Alexander M. Alekseenko. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

We derive two sets of explicit homogeneous algebraic constraint-preserving boundary conditions for the second-order in time reduction of the linearized Baumgarte-Shapiro-Shibata-Nakamura (BSSN) system. Our second-order reduction involves components of the linearized extrinsic curvature only. An initial-boundary value problem for the original linearized BSSN system is formulated and the existence of the solution is proved using the properties of the reduced system. A treatment is proposed for the full nonlinear BSSN system to construct constraint-preserving boundary conditions without invoking the second order in time reduction. Energy estimates on the principal part of the BSSN system (which is first order in temporal and second order in spatial derivatives) are obtained. Generalizations to the case of nonhomogeneous boundary data are proposed.

#### 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 and and the Neumann data for . Then (4.20) contains normal derivatives of , , and and, therefore, produces mixed conditions on these quantities. The difficulty using this condition, though, is that it contains tangential derivatives of the unknown fields and it is not obvious if it leads to a well-posed evolution of and . In general, well-posedness of the differential boundary conditions similar to (4.20) can be established by employing the techniques of the Laplace-Fourier transform developed in [26, 34, 35]. However, in some cases, classical energy estimates still can be obtained by the standard energy techniques as is suggested in the following example.

Consider a combination of the Neumann and Dirichlet conditions By applying to the second equation in (4.20), to the third equation of (4.20), and subtracting the results from the normal derivative of the first equation, we derive (using (4.3) and (4.4) to eliminate derivatives) an evolution equation defined on the boundary One can use (4.22) to construct nonhomogeneous mixed Dirichlet and Neumann data that imply (4.21). Specifically, we may prescribe Dirichlet data for , , and that is compatible at corners but, otherwise, arbitrary. Then we can use (4.22) to evolve from the initial data along the boundary. Since (4.22) is essentially a wave equation defined on the boundary surface, it admits a well-posed numerical formulation assuming that additional compatibility conditions are specified at corners in order to guarantee the smoothness of the solution (cf. [4]). Once the values of are determined, we can use the second and third equations of (4.20) to determine nonhomogeneous Neumann data for and . Similarly, one can give the Dirichlet data for , , and and use (4.22) to determine the Dirichlet data for and the Neumann data for and . The corresponding nonhomogeneous algebraic conditions on then read If desired, energy estimates may be obtained on the solutions and by a two step argument using the standard techniques for hyperbolic equations. The boundary conditions given by (4.23) are analogous to the conditions introduced in [4] for the Einstein-Christoffel formulation and in [18] for the first-order reduction of the BSSN formulation.

#### 5. Evolution of and ; Second Order in Time Reduction

We will argue now that the boundary conditions (4.17) and (4.18) lead to a well-posed problem for the linearized BSSN system. We assume that the initial values and are determined from (3.1) and (3.6), respectively, and satisfy the constraint equations. Also, we suppose that the conditions (4.17) or (4.18) are specified at the domain's boundary.

We introduce scalar products and for the spaces of scalar fields and triple-indexed fields on . The norms naturally associated with these scalar products are and . Consider the system (4.3) describing the propagation of . We define the energy of the system to be As in Section 4, by differentiating in time, integrating terms with spatial derivatives by parts, and using (4.3), we obtain Since the right-hand side of (5.2) is zero if either (4.17) or (4.18) holds, we conclude that and, therefore, remain bounded.

Similarly consider (4.4). The boundary conditions (4.17) and (4.18) imply that the trivial Neumann or Dirichlet data are given for , respectively. In a similar fashion, we can show that a solution to (4.4) must stay bounded as well.

We take advantage of the simplicity of the homogeneous data case and invoke the energy methods developed in [36–39] to prove the existence of the solution to (4.3) and (4.4). By employing the standard first-order reduction for each component of (4.3) and (4.4), we can
reduce these equations to six uncoupled first-order symmetric hyperbolic
systems. Moreover, both conditions (4.17) and (4.18) yield maximally nonnegative boundary conditions for the resulting first-order systems. Furthermore, we can invoke the results of [39, 40] for each individual first-order system to formulate the following theorem.
Theorem 5.1. *Let the initial data , , , , , and be given, and let and be determined by (3.6) and (3.4), respectively. If and and also and , then there exists unique solutions and to (4.3) and (4.4). These solutions obey the following energy estimates:
*

#### 6. The Initial Boundary Value Problem for the Linearized BSSN: Prescribed and Gamma-Freezing Shift

We recall that conditions (4.20) may also be formally imposed in a numerical method. However, these conditions contain tangential derivatives of the unknown fields, therefore their well-posedness is not easy to establish (see, however, results of [16, 23, 26, 27] for a relevant treatment of the generalized harmonic formulation of Einstein's equations). Also, the conditions (4.23) are formally well-posed, but deriving the complete estimates on their solutions is difficult. We therefore restrict our attention to the case of either condition (4.17) or (4.18).

Let system (3.2)–(3.7) be provided with relevant initial data and boundary conditions for and be taken in either form (4.17) or (4.18). According to Theorem 5.1, the matrix can be computed from (4.3), and can be determined from (4.4). Let us assume that the shift perturbation is known. Integration of (3.5), (3.2), (3.3), and (3.7) with respect to time yields the boundary conditions on the variables , , , and , where the projection operator is defined below (6.4)
The projection operator corresponding to (4.17) is given by
and the one corresponding to (4.18) is
In (6.2) and (6.3), in the case of (4.17), we have ; in the case of (4.18) . Furthermore, in (6.4), can be either vector or . If the Dirichlet data is given on , and are coupled through an integral equation (6.4). If the Neumann data is specified for , for example, in (4.17), the last equation of (6.4) can be replaced by
Theorem 6.1. *Let and be the solutions of (4.3) and (4.4) given by Theorem 5.1 for either boundary condition (4.17) or (4.18). Let also the initial data , , , and and the initial data , , , and be such that
**
Then a solution to (3.2)–(3.7) satisfying either boundary condition (4.17) or (4.18) and also boundary conditions
(6.1)–(6.4) is given by
**
Moreover, if and satisfy the constraint equation (3.9), then , , and defined by (6.9) satisfy constraints (3.8) and (3.10) provided that they satisfy them at the initial time.**Proof. *Equations (3.2), (3.3), (3.5), and (3.7) are verified by substitution. Replacing , , , and in (3.6) and (3.4) by their expressions from (6.9) and using (6.8), we obtain
which is a consequence of (4.3) and (4.4). Substituting (6.9) into (6.1)–(6.4), we verify the boundary conditions.

Now consider constraints (3.8). Replacing , , and using (6.9), we obtain
It follows that (3.8) is met as
long as it is satisfied initially and (3.9) is true. Constraint (3.10) can be proven similarly.

The situation is similar when is to be determined from the gamma-freezing condition yielding an elliptic equation for that needs to be solved at each time step: The boundary conditions on can be taken, for example, to be or A complete well-posedness proof requires a careful consideration of the coupling between (6.12) and the main system (3.2)–(3.7). However, taking a scalar product of (6.12) and and integrating by parts, we can show that is bounded by using the standard techniques for elliptic equations.

For computational purposes is beneficial to replace the elliptic equation (6.12) with a time-dependent hyperbolic equation [41] to have which corresponds to a dynamic gamma-freezing condition . In this case, in addition to (6.4) and (6.13), one can consider either of two sets of radiative boundary conditions By reasoning similar to the above, an energy estimate for in terms of can be obtained using techniques for hyperbolic equations (cf. [37–39, 42]).

After the boundary conditions for , , and are chosen, one can set on the boundary. Then the conditions for the rest of the variables follow by integration with respect to time as in the previous examples.

#### 7. Energy Estimates for the BSSN System with Boundaries

We use differentiation in time to propose boundary conditions (4.17) and (4.18) (and the associated conditions (6.1)–(6.4)) in Sections 4–6. Our derivation, however, relies strongly on the linearization assumption and does not extend to the nonlinear case. The difficulty appears to be the extra derivatives of the inverse metric resulting from differentiation of (2.20) with respect to time. These extra terms contaminate the principal part and break the similarity with the linear case. In other words, we could have proposed boundary conditions for the nonlinear case on the basis of (4.17) (or (4.18)), but we will not be able to prove the well-posedness of the new conditions by repeating the argument of Sections 4–6. In this section, we try to correct this flaw by establishing an energy estimate without reduction to second order in time. We use approach proposed by Gundlach and Martin-Garcia [18, 22] for unbounded domain and extend their proof to the case of bounded domains.

Following [18], we use the densitized lapse which results in significant simplifications in our equations. The result, however, is expected to carry over to the harmonic slicing and the -driving slicing as well.

We illustrate the method in the linear case first. Under the densitized lapse assumption, (3.3) is replaced with the linearized densitized lapse condition. The latter in the special case reduces to . Assuming zero-shift perturbation for further simplicity, (3.2), (3.4)–(3.7) are restated as Taking the scalar product of (7.3) with , integrating the result by parts, and using (7.2) to replace with , we obtain Also, we conclude from (7.2) that .

From (7.2) and (7.6), we observe that , which implies that

Next we rewrite the right-hand side of (7.5) in the divergence form: Then we take the scalar product with , and integrate by parts to obtain Replacing with (7.4) in the second term and using , we replace the last identity with

From (7.7)–(7.11), we observe that the energy has its growth determined by the boundary terms

Expression (7.13) can be used to propose examples of new energy- (and constraint-) preserving boundary conditions for the linearized BSSN formulation. However, here we will use (7.13) to discuss the validity of conditions (4.17) and (4.18) proposed in Section 4. In the case of (4.17) and the associated conditions (6.1), (6.2), and (6.4), expression (7.13) reduces to where – are coefficients of the decomposition Energy (7.12) is conserved if the initial data is chosen as to satisfy , , and at the boundary. Notice that cannot be given freely but is expected to be subject to constraint (3.10).

Similarly, condition (4.18) is both constraint- and energy-preserving for (7.2)–(7.6) if , , and on .

We are now ready to establish the nonlinear analog of (7.13). Substituting (7.1) into the BSSN equations and distributing and expanding covariant derivatives, we rewrite (2.8), (2.14), (2.15), (2.20), and (2.24) as where Taking the scalar product of (7.16) with and integrating by parts on the right-hand side, we obtain where is the outward normal vector to the boundary in the sense of the conformal metric . Substituting (7.17) for in the second term and regrouping, we obtain our first energy identity: where , and Also, from (7.17), it follows that Next, we notice from (7.16) and (7.20) that Therefore, where Finally, taking scalar product of (7.19) with , integrating by parts in the right-hand side, using (7.18) to replace with , and regrouping, we derive our last energy identity: where , and Notice that boundary terms in (7.23) and (7.29) differ from those in the linear case in the spatial metric only. Notice also that the right-hand sides of (7.23)–(7.29) are combinations of , , , , , , and , but not derivatives of , , and , implying that (7.23)–(7.29) are a closed estimate and may be proposed for proving local well-posedness of the initial-boundary value problem for the BSSN system. In this derivation, we assumed that contravariant components of the conformal shift vector are prescribed. Otherwise, the terms have to be expanded to ensure that they do not break the hyperbolicity of the equations.

#### 8. The Concluding Remarks

We have summarized the methods available for the construction of constraint-preserving boundary conditions in numerical general relativity and applied them to the BSSN system. In the simplest case when the system is linearized around the flat space and the boundary data is homogeneous, we proved the well-posedness of the proposed conditions. The obtained results will be extended in future to prove the well-posedness of the nonhomogeneous boundary conditions that are used in actual simulations of Einstein equations. In particular, estimates of Section 7 will be extended to include the frozen coefficient regime and the linearization around an arbitrary background. Another development lies in the construction of well-posed Sommerfeld-type constraint-preserving boundary conditions for the BSSN formulation. A more thorough study of the propagation of the constraint quantities needs to be done and methods that allow for driving the perturbed solution to a constraint-satisfying solution are developed in the BSSN formulation. We expect that further study of the BSSN equations and generalization of the results presented in this paper will provide answers to these important questions.

#### Acknowledgments

The author thanks D. Arnold for helpful suggestions at the earlier stages of this paper, O. Sarbach and G. Nagy for helping to correct an error in the first version of the paper, and L. Lindblom and M. Scheel for fruitful discussions and interest in this work. The author also thanks H. Pfeiffer for reading the second draft of the manuscript and suggesting helpful comments. The author expresses his most cordial gratitude to the emeritus professor at CSUN, Dr. L.L. Foster, for her help in the preparation of the final version of this paper.

#### References

- A. Anderson and J. W. York, Jr., βFixing Einstein's equations,β
*Physical Review Letters*, vol. 82, no. 22, pp. 4384β4387, 1999. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH Β· View at MathSciNet - T. W. Baumgarte and S. L. Shapiro, βNumerical integration of Einstein's field equations,β
*Physical Review D*, vol. 59, no. 2, Article ID 024007, p. 7, 1999. View at Publisher Β· View at Google Scholar Β· View at MathSciNet - G. Calabrese, L. Lehner, and M. Tiglio, βConstraint-preserving boundary conditions in numerical relativity,β
*Physical Review D*, vol. 65, no. 10, Article ID 104031, p. 13, 2002. View at Publisher Β· View at Google Scholar Β· View at MathSciNet - G. Calabrese, J. Pullin, O. A. Reula, O. Sarbach, and M. Tiglio, βWell posed constraint-preserving boundary conditions for the linearized Einstein equations,β
*Communications in Mathematical Physics*, vol. 240, no. 1-2, pp. 377β395, 2003. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH Β· View at MathSciNet - H. Friedrich and G. Nagy, βThe initial boundary value problem for Einstein's vacuum field equation,β
*Communications in Mathematical Physics*, vol. 201, no. 3, pp. 619β655, 1999. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH Β· View at MathSciNet - S. Frittelli and O. A. Reula, βFirst-order symmetric hyperbolic Einstein equations with arbitrary fixed gauge,β
*Physical Review Letters*, vol. 76, no. 25, pp. 4667β4670, 1996. View at Publisher Β· View at Google Scholar - L. E. Kidder, M. A. Scheel, and S. A. Teukolsky, βExtending the lifetime of 3D black hole computations with a new hyperbolic system of evolution equations,β
*Physical Review D*, vol. 64, no. 6, Article ID 064017, p. 13, 2001. View at Publisher Β· View at Google Scholar Β· View at MathSciNet - M. Shibata and T. Nakamura, βEvolution of three-dimensional gravitational waves: harmonic slicing case,β
*Physical Review D*, vol. 52, no. 10, pp. 5428β5444, 1995. View at Publisher Β· View at Google Scholar Β· View at MathSciNet - B. Szilágyi and J. Winicour, βWell-posed initial-boundary evolution in general relativity,β
*Physical Review D*, vol. 68, no. 4, Article ID 041501, p. 5, 2003. View at Publisher Β· View at Google Scholar Β· View at MathSciNet - O. Brodbeck, S. Frittelli, P. Hübner, and O. A. Reula, βEinstein's equations with asymptotically stable constraint propagation,β
*Journal of Mathematical Physics*, vol. 40, no. 2, pp. 909β923, 1999. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH Β· View at MathSciNet - M. A. Scheel, L. E. Kidder, L. Lindblom, H. P. Pfeiffer, and S. A. Teukolsky, βToward stable 3D numerical evolutions of black-hole spacetimes,β
*Physical Review D*, vol. 66, no. 12, Article ID 124005, p. 4, 2002. View at Publisher Β· View at Google Scholar Β· View at MathSciNet - M. Holst, L. Lindblom, R. Owen, H. P. Pfeiffer, M. A. Scheel, and L. E. Kidder, βOptimal constraint projection for hyperbolic evolution systems,β
*Physical Review D*, vol. 70, no. 8, Article ID 084017, p. 17, 2004. View at Publisher Β· View at Google Scholar Β· View at MathSciNet - L. Lindblom, M. A. Scheel, L. E. Kidder, H. P. Pfeiffer, D. Shoemaker, and S. A. Teukolsky, βControlling the growth of constraints in hyperbolic evolution systems,β
*Physical Review D*, vol. 69, no. 12, Article ID 124025, p. 14, 2004. View at Publisher Β· View at Google Scholar Β· View at MathSciNet - M. Alcubierre, G. Allen, B. Brügmann, E. Seidel, and W.-M. Suen, βTowards an understanding of the stability properties of the $3+1$ evolution equations in general relativity,β
*Physical Review D*, vol. 62, no. 12, Article ID 124011, 15 pages, 2000. View at Publisher Β· View at Google Scholar Β· View at MathSciNet - B. Kelly, P. Laguna, K. Lockitch et al., βCure for unstable numerical evolutions of single black holes: adjusting the standard ADM equations in the spherically symmetric case,β
*Physical Review D*, vol. 64, no. 8, Article ID 084013, 14 pages, 2001. View at Publisher Β· View at Google Scholar - L. E. Kidder, L. Lindblom, M. A. Scheel, L. T. Buchman, and H. P. Pfeiffer, βBoundary conditions for the Einstein evolution system,β
*Physical Review D*, vol. 71, Article ID 064020, 22 pages, 2005. View at Publisher Β· View at Google Scholar - D. N. Arnold and N. Tarfulea, βBoundary conditions for the Einstein-Christoffel formulation of Einstein's equations,β in
*Proceedings of the 6th Mississippi State—UBA Conference on Differential Equations and Computational Simulations*, vol. 15 of*Electronic Journal of Differential Equations Conference*, pp. 11β27, Southwest Texas State University, San Marcos, Tex, USA, 2007. View at Zentralblatt MATH Β· View at MathSciNet - C. Gundlach and J. M. Martín-García, βSymmetric hyperbolicity and consistent boundary conditions for second-order Einstein equations,β
*Physical Review D*, vol. 70, no. 4, Article ID 044032, p. 16, 2004. View at Publisher Β· View at Google Scholar Β· View at MathSciNet - S. Frittelli and R. Gómez, βEinstein boundary conditions for the $3+1$ Einstein equations,β
*Physical Review D*, vol. 68, no. 4, Article ID 044014, p. 6, 2003. View at Publisher Β· View at Google Scholar Β· View at MathSciNet - O. Sarbach and M. Tiglio, βBoundary conditions for Einstein's field equations: mathematical and numerical analysis,β
*Journal of Hyperbolic Differential Equations*, vol. 2, no. 4, pp. 839β883, 2005. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH Β· View at MathSciNet - N. Tarfulea,
*Constraint-preserving boundary conditions for Einstein's equations*, Ph.D. dissertation, University of Minnesota, Minneapolis, Minn, USA, 2004. View at Zentralblatt MATH Β· View at MathSciNet - C. Gundlach and J. M. Martín-García, βSymmetric hyperbolic form of systems of second-order evolution equations subject to constraints,β
*Physical Review D*, vol. 70, no. 4, Article ID 044031, p. 14, 2004. View at Publisher Β· View at Google Scholar Β· View at MathSciNet - O. Rinne, βStable radiation-controlling boundary conditions for the generalized harmonic Einstein equations,β
*Classical and Quantum Gravity*, vol. 23, no. 22, pp. 6275β6300, 2006. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH Β· View at MathSciNet - O. A. Reula and O. Sarbach, βA model problem for the initial-boundary value formulation of Einstein's field equations,β
*Journal of Hyperbolic Differential Equations*, vol. 2, no. 2, pp. 397β435, 2005. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH Β· View at MathSciNet - O. Sarbach, βAbsorbing boundary conditions for Einstein's field equations,β
*Journal of Physics*, vol. 91, no. 1, Article ID 012005, 15 pages, 2007. View at Publisher Β· View at Google Scholar - H.-O. Kreiss and J. Winicour, βProblems which are well-posed in a generalized sense with applications to the Einstein equations,β
*Classical and Quantum Gravity*, vol. 23, no. 16, pp. S405βS420, 2006. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH Β· View at MathSciNet - H.-O. Kreiss, O. A. Reula, O. Sarbach, and J. Winicour, βWell-posed initial-boundary value problem for the harmonic Einstein equations using energy estimates,β
*Classical and Quantum Gravity*, vol. 24, no. 23, pp. 5973β5984, 2007. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH Β· View at MathSciNet - M. Alcubierre, B. Brügmann, P. Diener et al., βGauge conditions for long-term numerical black hole evolutions without excision,β
*Phys*