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

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 3+1 Arnowitt-Deser-Misner (ADM) decomposition [29, 30]: 𝜕𝜕₀ℎ_𝑖𝑗=−2𝑎𝑘_𝑖𝑗+2ℎ_𝑙(𝑖𝜕_𝑗)𝑏^𝑙,(2.1)₀𝑘_𝑖𝑗=𝑎𝑅_𝑖𝑗+𝑘𝑅^𝑙𝑙𝑘^𝑖𝑗−2𝑘_𝑖𝑙𝑘^𝑙𝑗+𝑘^𝑖𝑙𝜕_𝑗𝑏^𝑙+𝑘_𝑙𝑗𝜕_𝑖𝑏^𝑙−𝐷_𝑖𝐷_𝑗𝑎,(2.2)_𝑖𝑖+𝐷𝑘^𝑖𝑖²−𝑘_𝑖𝑗𝑘^𝑖𝑗=0,(2.3)_𝑗𝑘_𝑖𝑗−𝐷_𝑖𝑘^𝑗𝑗=0.(2.4) 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 𝑔₀₀=−𝑎²+𝑏_𝑖𝑏_𝑗ℎ^𝑖𝑗,𝑔_0𝑖=𝑏_𝑖,𝑔_𝑖𝑗=ℎ_𝑖𝑗.(2.5) 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 𝑅_𝑖𝑗=12ℎ^𝑝𝑞𝜕^𝑝𝜕_𝑗ℎ_𝑖𝑞+𝜕_𝑖𝜕_𝑝ℎ_𝑞𝑗−𝜕_𝑝𝜕_𝑞ℎ_𝑖𝑗−𝜕_𝑖𝜕_𝑗ℎ_𝑝𝑞+ℎ^𝑝𝑞ℎ^𝑟𝑠Γ^𝑖𝑝𝑠Γ_𝑞𝑗𝑟−Γ_𝑝𝑞𝑠Γ_𝑖𝑗𝑟,(2.6) where Γ_𝑖𝑗𝑘 are the spatial Christoffel symbols defined by Γ_𝑖𝑗𝑘=(𝜕_𝑖ℎ_𝑘𝑗+𝜕_𝑗ℎ_𝑖𝑘−𝜕_𝑘ℎ_𝑗𝑖)/2.

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 𝜕₀𝑘^𝑖𝑖=𝑎𝑘^𝑙𝑖𝑘_𝑙𝑖−𝐷^𝑙𝐷_𝑙𝑎.(2.7) 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 𝐴_𝑖𝑗=𝑘_𝑖𝑗−(1/3)ℎ_𝑖𝑗𝑘 as new variables. Then (2.7) yields 𝜕₀𝑘=13𝑎𝑘²+𝑎𝐴^𝑙𝑚𝐴_𝑙𝑚−𝐷^𝑙𝐷_𝑙𝑎.(2.8) 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 𝜕_𝑡𝑘≥(1/3)𝑎𝑘² which implies that 𝑘≥[(1/3)∫^𝑡0𝑎(𝜏)𝑑𝜏+1/𝑘(0)]⁻¹, 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] 𝐷^𝑙𝐷_𝑙𝑎=𝑎𝑘^𝑙𝑚𝑘_𝑙𝑚;(2.9) with this condition, (2.8) reduces to 𝜕₀𝑘=0.

Alternatively, we can use harmonic slicing [2, 28] (a particular case of the Bona-Massó family of 𝑘-driving slicing conditions (𝜕_𝑡−𝑏^𝑙𝐷_𝑙)𝑎=−𝑎²𝑓(𝑎)𝑘, 𝑓(𝑎)>0 [31, 32]), which corresponds to setting 𝜕₀𝑎=−𝑎²𝑘.(2.10)

The equation on 𝐴 is obtained from (2.2), (2.8), and (2.1) as 𝜕₀𝐴_𝑖𝑗=𝑎𝑅_𝑖𝑗+13𝑎𝑘𝐴^𝑖𝑗−2𝑎𝐴_𝑖𝑙𝐴^𝑙𝑗+29𝑎𝑘²ℎ_𝑖𝑗−13𝑎𝐴^𝑙𝑚𝐴_𝑙𝑚ℎ_𝑖𝑗+𝐴_𝑖𝑙𝜕_𝑗𝑏^𝑙+𝐴_𝑗𝑙𝜕_𝑖𝑏^𝑙−𝐷_𝑖𝐷_𝑗𝑎+13ℎ^𝑖𝑗𝐷^𝑙𝐷_𝑙𝑎.(2.11) 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 𝜑=(1/12)ln(det(ℎ_𝑖𝑗)) and the conformal metric ℎ_𝑖𝑗=𝑒^−4𝜑ℎ_𝑖𝑗, ℎ^𝑖𝑗=𝑒^4𝜑ℎ^𝑖𝑗. Using the Leibnitz formula for differentiating the determinant of a matrix 𝜕det(ℎ_𝑖𝑗)=det(ℎ_𝑖𝑗)ℎ^𝑙𝑚𝜕ℎ_𝑚𝑙,(2.12) one finds that the derivative 𝜕ℎ_𝑖𝑗=𝑒^−4𝜑𝜕ℎ^𝑖𝑗−13ℎ^𝑖𝑗ℎ^𝑙𝑚𝜕ℎ_𝑙𝑚(2.13) 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 𝜕₀𝜑=−16𝑎𝑘+16𝜕^𝑙𝑏^𝑙.(2.14) Now using (2.13) and (2.1), we obtain the third equation 𝜕₀ℎ_ij=−2𝑎𝐴_𝑖𝑗+2ℎ_𝑙(𝑖𝜕_𝑗)𝑏^𝑙−23ℎ^𝑖𝑗𝜕_𝑙𝑏^𝑙,(2.15) where 𝐴_𝑖𝑗=𝑒^−4𝜑𝐴_𝑖𝑗 and 𝑏_𝑖=𝑒^−4𝜑𝑏_𝑖 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 ℎ^𝑖𝑗=𝑒^4𝜑ℎ^𝑖𝑗. 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 𝜕₀𝐴_𝑖𝑗=𝑎𝑒^−4𝜑𝑅_𝑖𝑗+𝑎𝑘𝐴_𝑖𝑗−2𝐴_𝑖𝑙𝐴^𝑙𝑗+29𝑘²ℎ_𝑖𝑗−13𝐴^𝑙𝑚+𝐴_𝑖𝑙𝐴_𝑙𝑚ℎ_𝑖𝑗𝜕_𝑗𝑏^𝑙+𝐴_𝑗𝑙𝜕_𝑖𝑏^𝑙−23𝐴^𝑖𝑗𝜕_𝑙𝑏^𝑙−𝑒^−4𝜑𝐷_𝑖𝐷_𝑗𝑎+𝑒^−4𝜑13ℎ^𝑖𝑗𝐷^𝑙𝐷_𝑙𝑎.(2.16) The Ricci tensor in terms of the conformal metric becomes [2, 28] 𝑅_𝑖𝑗=12ℎ^𝑝𝑞𝜕^𝑝+ℎ_𝑝𝑞𝜕_𝑗ℎ_𝑖𝑞+𝜕_𝑖𝜕_𝑝ℎ_𝑞𝑗−𝜕_𝑝𝜕_𝑞ℎ_𝑖𝑗−𝜕_𝑖𝜕_𝑗ℎ_𝑝𝑞−2𝐷^𝑖𝐷_𝑗𝜑−2ℎ_𝑖𝑗ℎ^𝑝𝑞𝐷_𝑝𝐷_𝑞𝜑ℎ^𝑟𝑠Γ^𝑖𝑝𝑠Γ_𝑞𝑗𝑟−Γ_𝑝𝑞𝑠Γ_𝑖𝑗𝑟+4𝜕^𝑖𝜑𝜕_𝑗𝜑−4ℎ_𝑖𝑗ℎ^𝑝𝑞𝜕_𝑝𝜑𝜕_𝑞𝜑.(2.17) Here Γ_𝑖𝑗𝑘=(𝜕_𝑖ℎ_𝑘𝑗+𝜕_𝑗ℎ_𝑖𝑘−𝜕_𝑘ℎ_𝑗𝑖)/2 and 𝐷_𝑖𝑣_𝑗=𝜕_𝑖𝑣_𝑗−ℎ^𝑝𝑞Γ_𝑖𝑗𝑝𝑣_𝑞 are the covariant derivatives associated with the conformal metric. The first line in (2.17) can be rewritten to obtain 𝑅_𝑖𝑗=−12ℎ^𝑝𝑞+ℎ_𝑝𝑞𝜕_𝑝𝜕_𝑞ℎ_𝑖𝑗+𝜕_(𝑖ℎ^𝑝𝑞Γ||_𝑝𝑞||_𝑗)+Γ_{𝑝𝑞(𝑖}𝜕_𝑗)ℎ^𝑝𝑞−2𝐷_𝑖𝐷_𝑗𝜑−2ℎ_𝑖𝑗ℎ^𝑝𝑞𝐷_𝑝𝐷_𝑞𝜑ℎ^𝑟𝑠Γ^𝑖𝑝𝑠Γ_𝑞𝑗𝑟−Γ_𝑝𝑞𝑠Γ_𝑖𝑗𝑟+4𝜕^𝑖𝜑𝜕_𝑗𝜑−4ℎ_𝑖𝑗ℎ^𝑝𝑞𝜕_𝑝𝜑𝜕_𝑞𝜑.(2.18) This suggests that we introduce a new variable Γ_𝑗=ℎ^𝑝𝑞Γ_𝑝𝑞𝑗=ℎ^𝑝𝑞𝜕_𝑝ℎ_𝑞𝑗.(2.19) Substituting (2.18) in (2.16), we derive the fourth evolution equation 𝜕₀𝐴_𝑖𝑗=−12𝑎𝑒^−4𝜑−2𝑎𝑒_−4𝜑ℎ^𝑝𝑞𝜕_𝑝𝜕_𝑞ℎ_𝑖𝑗+𝑎𝑒^−4𝜑𝜕_(𝑖Γ_𝑗)−2𝑎𝑒^−4𝜑𝐷_𝑖𝐷_𝑗𝜑ℎ_𝑖𝑗ℎ^𝑝𝑞𝐷_𝑝𝐷_𝑞𝜑−𝑒^−4𝜑𝐷_𝑖𝐷_𝑗𝑎+13𝑒^−4𝜑ℎ_𝑖𝑗𝐷^𝑙𝐷_𝑙𝑎+𝑊_𝑖𝑗,(2.20) where 𝑊_𝑖𝑗=+𝑎𝑒4𝑎𝑒^−4𝜑Γ_{𝑝𝑞(𝑖−4𝜑}𝜕_𝑗)ℎ^𝑝𝑞+𝑎𝑒^−4𝜑ℎ^𝑝𝑞ℎ^𝑟𝑠Γ^𝑖𝑝𝑠Γ_𝑞𝑗𝑟−Γ_𝑝𝑞𝑠Γ_𝑖𝑗𝑟+𝑎𝑘𝐴_𝑖𝑗𝜕_𝑖𝜑𝜕_𝑗𝜑−4𝑎𝑒^−4𝜑ℎ_𝑖𝑗ℎ^𝑝𝑞𝜕_𝑝𝜑𝜕_𝑞𝜑−2𝐴_𝑖𝑙𝐴^𝑙𝑗+29𝑘²ℎ_𝑖𝑗−13𝐴^𝑙𝑚+𝐴_𝑖𝑙𝐴_𝑙𝑚ℎ_𝑖𝑗𝜕_𝑗𝑏^𝑙+𝐴_𝑗𝑙𝜕_𝑖𝑏^𝑙−23𝐴^𝑖𝑗𝜕_𝑙𝑏^𝑙.(2.21)

The evolution equation for Γ_𝑗 is obtained by applying the operator 𝜕₀ to (2.19) and using the momentum constraint 𝜕₀Γ_𝑖=−2𝑎𝜕^𝑝𝐴_𝑝𝑖−2ℎ^𝑝𝑞𝜕^𝑝𝑎𝐴_𝑞𝑖+2𝑎𝐴^𝑝𝑞𝜕^𝑝ℎ_𝑞𝑖+Γ^𝑙𝜕_𝑖𝑏^𝑙+13𝜕^𝑖𝜕_𝑙𝑏^𝑙+ℎ_𝑙𝑖𝜕^𝑠𝜕_𝑠𝑏^𝑙.(2.22) Next we notice that ℎ^𝑝𝑞𝐷_𝑝𝐴_𝑖𝑞=𝜕^𝑝𝐴_𝑝𝑖−Γ_𝑠𝐴^𝑠𝑖+6(𝜕_𝑠𝜑)𝐴^𝑠𝑖, and thus (2.4) takes the form 𝜕^𝑝𝐴_𝑝𝑖−23𝜕^𝑖𝑘−Γ_𝑠𝐴^𝑠𝑖+6(𝜕_𝑠𝜑)𝐴^𝑠𝑖=0.(2.23) Solving this equation for 𝜕^𝑝𝐴_𝑝𝑖 and substituting the result in (2.22), we derive the fifth equation of the BSSN system 𝜕₀Γ_𝑖=−43𝑎𝜕^𝑖𝑘+𝑆_𝑖,(2.24) where 𝑆_𝑖=−2𝑎Γ_𝑠𝐴^𝑠𝑖+12𝑎(𝜕_𝑠𝜑)𝐴^𝑠𝑖−2ℎ^𝑝𝑞𝜕^𝑝𝑎𝐴_𝑞𝑗+2𝑎𝐴^𝑝𝑞𝜕^𝑝ℎ_𝑞𝑗+Γ^𝑙𝜕_𝑗𝑏^𝑙+13𝜕^𝑖𝜕_𝑙𝑏^𝑙+ℎ_𝑙𝑗𝜕^𝑠𝜕_𝑠𝑏^𝑙.(2.25) 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 𝜕_𝑡Γ_𝑖=0 (cf. [28]).

3. Linearization around Minkowski Space

Minkowski spacetime in Cartesian coordinates is represented by the trivial solution to the ADM system: ℎ_𝑖𝑗=𝛿_𝑖𝑗, 𝑘_𝑖𝑗=0, 𝑎=1, 𝑏_𝑖=0. Consider perturbations of the ADM variables around the Minkowski spacetime: ℎ_𝑖𝑗=𝛿_𝑖𝑗+𝛾_𝑖𝑗, 𝑘_𝑖𝑗=𝜅_𝑖𝑗, 𝑎=1+𝛼, 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 det(ℎ_𝑖𝑗)=1+𝛾^𝑙𝑙,𝜑=112𝛾^𝑙𝑙,𝑒^−4𝜑=1−13𝛾^𝑙𝑙,𝑒^4𝜑=1+13𝛾^𝑙𝑙ℎ_𝑖𝑗,=𝛿_𝑖𝑗+𝛾_𝑖𝑗−13𝛿^𝑖𝑗𝐴_𝑖𝑗𝛾^𝑙𝑙=∶𝛿_𝑖𝑗+𝛾_𝑖𝑗,𝑘=𝜅=∶𝜅^𝑙𝑙,=𝐴_𝑖𝑗=𝜅_𝑖𝑗−13𝛿^𝑖𝑗𝜅,Γ_𝑖=∶Γ_𝑖=𝜕^𝑙𝛾_𝑖𝑙.(3.1) 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 𝜕_𝑡𝜑=−16𝜅+61𝜕^𝑙𝜕𝛽_𝑙,(3.2)_𝑡𝜕𝛼=−𝜅,(3.3)_𝑡𝜕𝜅=−𝜕^𝑙𝜕_𝑙𝛼,(3.4)_𝑡𝛾_𝑖𝑗=−2𝐴_𝑖𝑗+2𝜕_(𝑖𝛽_𝑗)−23𝛿^𝑖𝑗𝜕𝜕^𝑙𝛽_𝑙,(3.5)_𝑡𝐴_𝑖𝑗=−12𝜕^𝑙𝜕_𝑙𝛾_𝑖𝑗+𝜕_(𝑖Γ_𝑗)−2𝜕_𝑖𝜕_𝑗𝜑−2𝛿_𝑖𝑗𝜕^𝑙𝜕_𝑙𝜑−𝜕_𝑖𝜕_𝑗𝛼+13𝛿^𝑖𝑗𝜕𝜕^𝑙𝜕_𝑙𝛼,(3.6)_𝑡Γ_𝑖=−43𝜕^𝑖𝜅+13𝜕^𝑖𝜕^𝑝𝛽_𝑝+𝜕^𝑝𝜕_𝑝𝛽_𝑖.(3.7)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 𝜕^𝑙𝜕^𝑗𝛾_𝑙𝑗−8𝜕^𝑙𝜕_𝑙𝜕𝜑=0,𝜕^𝑙Γ_𝑙−8𝜕^𝑙𝜕_𝑙𝜑=0,(3.8)_𝑙𝐴_𝑖𝑙−23𝜕^𝑖𝜅=0.(3.9) 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 Γ_𝑖=𝜕^𝑙𝛾_𝑖𝑙.(3.10)

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 𝜑(0), 𝜅(0), 𝛾(0), 𝐴(0), Γ(0) can be determined from 𝛾(0) and 𝜅(0) using (3.1). It can be verified that if 𝛾(0) and 𝜅(0) satisfy the linearized Hamiltonian and momentum constraints of the ADM system, then 𝜑(0), 𝜅(0), 𝛾(0), 𝐴(0), Γ(0) satisfy the constraint equations (3.8)–(3.10). The linearized Hamiltonian and momentum constraints in the ADM system are 𝜕^𝑙𝜕^𝑖𝛾_𝑙𝑖−𝜕^𝑙𝜕_𝑙𝛾^𝑖𝑖=0 and 𝜕^𝑙𝜅_𝑖𝑙−𝜕_𝑖𝜅^𝑙𝑙=0.

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 𝜕_𝑡𝜕^𝑙𝜕^𝑗𝛾_𝑙𝑗−8𝜕^𝑙𝜕_𝑙𝜑=𝜕^𝑖𝜕^𝑙𝐴_𝑖𝑙−23𝜕^𝑖𝜕_𝑡𝜅,𝜕^𝑙Γ_𝑙−8𝜕^𝑙𝜕_𝑙𝜑=0.(4.1) 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 𝑀_𝑖=𝜕^𝑙𝐴_𝑖𝑙−23𝜕^𝑖𝜅.(4.2) Equation (3.9) is satisfied if and only if 𝑀_𝑖=0.

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 𝜕^2𝑡𝐴_𝑖𝑗=𝜕^𝑙𝜕_𝑙𝐴_𝑖𝑗.(4.3) Terms in 𝜑, 𝛼, 𝜅, 𝛽, Γ cancel in view of (3.2)–(3.4), (3.7).

Similarly, by differentiating (3.4) and substituting (3.3) for 𝜕_𝑡𝛼, one derives 𝜕^2𝑡𝜅=𝜕^𝑙𝜕_𝑙𝜅.(4.4) Differentiating (4.2) twice with respect to time and using (4.3) and (4.4), we obtain 𝜕^2𝑡𝑀_𝑖=𝜕^𝑙𝜕_𝑙𝑀_𝑖.(4.5) The initial value 𝑀(0) can be determined from 𝐴(0) and 𝜅(0) using (4.2) and must be zero for the physical initial data. The initial values 𝜕_𝑡𝑀(0) can be calculated by differentiating (4.2) in time and substituting (3.4) and (3.6) for 𝜕_𝑡𝜅 and 𝜕_𝑡𝐴_𝑖𝑗. Thus we have 𝜕_𝑡𝑀_𝑖=−12𝜕^𝑙𝜕_𝑙𝜕^𝑚𝛾_𝑖𝑚+12𝜕^𝑖𝜕^𝑙Γ_𝑙+12𝜕^𝑙𝜕_𝑙Γ_𝑖−4𝜕_𝑖𝜕^𝑙𝜕_𝑙𝜑.(4.6) It can be verified by substitution that if 𝛾(0), Γ(0), and 𝜑(0) satisfy (3.8) and (3.10) then 𝜕_𝑡𝑀_𝑖(0)=0.

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 𝑀(0) and 𝜕_𝑡𝑀(0) is determined by 𝐴(0), 𝜅(0), 𝛾(0), Γ(0), and 𝜑(0). 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 𝜖=‖‖𝜕_𝑡𝑀‖‖²+‖‖𝜕_𝑙𝑀‖‖².(4.7) If we prove that the energy 𝜖 remains zero at all times, then, in view of the trivial initial data, we will have ‖𝑀‖≡0 (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 𝜕_𝑡𝜖=_𝜕Ω𝜕𝜕𝑛𝑀^𝑖𝜕^𝑡𝑀^𝑖𝑑𝜎.(4.8) The energy 𝜖 is not increasing if 𝜕𝜕𝑛𝑀^𝑖𝜕^𝑡𝑀^𝑖≤0on𝜕Ω.(4.9) 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 𝜕^𝑖𝑣_𝑖=𝑛^𝑖𝜕𝜕𝑛𝑣^𝑖+𝑚^𝑖𝜕𝜕𝑚𝑣^𝑖+𝑙^𝑖𝜕𝜕𝑙𝑣^𝑖.(4.10) Similarly, the gradient of a scalar field 𝜓 becomes 𝜕_𝑖𝜓=𝑛_𝑖𝜕𝜕𝑛𝜓+𝑚^𝑖𝜕𝜕𝑚𝜓+𝑙^𝑖𝜕𝜕𝑙𝜓.(4.11) Next we note that at any point of the boundary, a symmetric trace free matrix is spanned by 𝑛_(𝑖𝑚_𝑗),𝑛_(𝑖𝑙_𝑗),𝑙_(𝑖𝑚_𝑗),𝑙_𝑖𝑙_𝑗−𝑚_𝑖𝑚_𝑗,2𝑛_𝑖𝑛_𝑗−𝑙_𝑖𝑙_𝑗−𝑚_𝑖𝑚_𝑗.(4.12) Introducing the scalar functions 𝐴𝐴41==21𝐴^𝑖𝑗𝑛_(𝑖𝑚_𝑗),𝐴2=2𝐴^𝑖𝑗𝑛_(𝑖𝑙_𝑗),𝐴3=2𝐴^𝑖𝑗𝑙_(𝑖𝑚_𝑗),2𝐴^𝑖𝑗𝑙^𝑖𝑙_𝑗−𝑚_𝑖𝑚_𝑗,𝐴5=16𝐴^𝑖𝑗2𝑛^𝑖𝑛_𝑗−𝑙_𝑖𝑙_𝑗−𝑚_𝑖𝑚_𝑗,(4.13) we rewrite 𝐴 as follows: 𝐴_𝑖𝑗=𝐴1𝑛_(𝑖𝑚_𝑗)+𝐴2𝑛^(𝑖𝑙_𝑗)+𝐴3𝑙^(𝑖𝑚_𝑗)+𝐴4𝑙^𝑖𝑙_𝑗−𝑚_𝑖𝑚_𝑗+𝐴52𝑛^𝑖𝑛_𝑗−𝑙_𝑖𝑙_𝑗−𝑚_𝑖𝑚_𝑗(4.14) Substituting (4.14) into (4.2) and using (4.10) and (4.11) to replace partial derivatives with directional derivatives, we obtain 𝑀_𝑖=12𝜕𝜕𝑚𝐴1+21𝜕𝜕𝑙𝐴2+2𝜕𝜕𝑛𝐴5−32𝜕𝜕𝑛𝜅𝑛^𝑖+12𝜕𝜕𝑛𝐴1+21𝜕𝜕𝑙𝐴3−𝜕𝜕𝑚𝐴4−𝜕𝜕𝑚𝐴5−23𝜕𝜕𝑚𝜅𝑚^𝑖+12𝜕𝜕𝑛𝐴2+21𝜕𝜕𝑚𝐴3+𝜕𝜕𝑙𝐴4−𝜕𝜕𝑙𝐴5−23𝜕𝜕𝑙𝜅𝑙^𝑖.(4.15) The last equation implies that 𝜕𝜕𝑛𝑀^𝑖𝜕^𝑡𝑀^𝑖=×𝜕𝜕𝑛𝜕_𝑡12𝜕𝜕𝑚𝐴1+21𝜕𝜕𝑙𝐴2+2𝜕𝜕𝑛𝐴5−32𝜕𝜕𝑛𝜅+𝜕21𝜕𝜕𝑚𝐴1+21𝜕𝜕𝑙𝐴2+2𝜕𝜕𝑛𝐴5−32𝜕𝜕𝑛𝜅×𝜕𝜕_𝑡𝑛12𝜕𝜕𝑛𝐴1+21𝜕𝜕𝑙𝐴3−𝜕𝜕𝑚𝐴4−𝜕𝜕𝑚𝐴5−23𝜕𝜕𝑚𝜅+𝜕21𝜕𝜕𝑛𝐴1+21𝜕𝜕𝑙𝐴3−𝜕𝜕𝑚𝐴4−𝜕𝜕𝑚𝐴5−23𝜕𝜕𝑚𝜅×𝜕𝜕_𝑡𝑛12𝜕𝜕𝑛𝐴2+21𝜕𝜕𝑚𝐴3+𝜕𝜕𝑙𝐴4−𝜕𝜕𝑙𝐴5−23𝜕𝜕𝑙𝜅21𝜕𝜕𝑛𝐴2+21𝜕𝜕𝑚𝐴3+𝜕𝜕𝑙𝐴4−𝜕𝜕𝑙𝐴5−23𝜕𝜕𝑙𝜅.(4.16) Either of the following two sets of boundary conditions imply that ((𝜕/𝜕𝑛)𝑀^𝑖)(𝜕_𝑡𝑀_𝑖)=0 on 𝜕Ω: 𝐴1=0,𝐴2=0,𝜕𝜕𝜕𝑛𝐴3=0,𝜕𝜕𝑛𝐴4=0,𝜕𝜕𝑛𝐴5=0,𝜕𝜕𝑛𝜅=0;(4.17)𝜕𝑛𝐴1=0,𝜕𝜕𝑛𝐴2=0,𝐴3=0,𝐴4=0,𝐴5=0,𝜅=0.(4.18) 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 𝐴3, 𝐴4, 𝐴5, and 𝜅 drop out as a direct consequence of (4.17). Also (4.3) can be solved for the second normal derivatives of 𝐴1 and 𝐴2 in terms of the temporal and tangential derivatives using the following representation of the Laplace operator: 𝜕^𝑗𝜕_𝑗=𝜕𝜕𝑛𝜕𝜕𝑛+𝜕𝜕𝑚𝜕𝜕𝑚+𝜕𝜕𝑙𝜕𝜕𝑙.(4.19) It can be easily noticed that (4.17) implies that all tangential and temporal derivatives of 𝐴1 and 𝐴2 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, 𝑀_𝑖|_𝜕Ω=0 is equivalent to the set of the following differential boundary conditions that can be used in a numerical method [13, 16]: 112𝜕𝜕𝑚𝐴1+12𝜕𝜕𝑙𝐴2+2𝜕𝜕𝑛𝐴5−32𝜕𝜕𝑛𝜅=0,21𝜕𝜕𝑛𝐴1+21𝜕𝜕𝑙𝐴3−𝜕𝜕𝑚𝐴4−𝜕𝜕𝑚𝐴5−23𝜕𝜕𝑚𝜅=0,2𝜕𝜕𝑛𝐴2+21𝜕𝜕𝑚𝐴3+𝜕𝜕𝑙𝐴4−𝜕𝜕𝑙𝐴5−23𝜕𝜕𝑙𝜅=0.(4.20) In this case, we can prescribe the Dirichlet data for 𝐴3 and 𝐴4 and the Neumann data for 𝜅. Then (4.20) contains normal derivatives of 𝐴1, 𝐴2, and 𝐴5 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 𝜕𝜕𝑛𝑀^𝑖𝑛^𝑖=0,𝑀_𝑖𝑙^𝑖=0,𝑀_𝑖𝑚^𝑖=0.(4.21) 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 2𝜕²𝜕𝑡²𝐴5−𝜕²𝜕𝑙²+𝜕²𝜕𝑚²𝐴5−32𝜕𝜕𝑡²²𝜅+43𝜕𝜕𝑙²²+𝜕²𝜕𝑚²𝜅=𝜕𝜕𝑙²²−𝜕²𝜕𝑚²𝐴4+𝜕𝜕𝑙𝜕𝜕𝑚𝐴3.(4.22) One can use (4.22) to construct nonhomogeneous mixed Dirichlet and Neumann data that imply (4.21). Specifically, we may prescribe Dirichlet data for 𝐴3, 𝐴4, and 𝜅 that is compatible at corners but, otherwise, arbitrary. Then we can use (4.22) to evolve 𝐴5 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 𝐴5 are determined, we can use the second and third equations of (4.20) to determine nonhomogeneous Neumann data for 𝐴1 and 𝐴2. Similarly, one can give the Dirichlet data for 𝐴3, 𝐴4, and 𝐴5 and use (4.22) to determine the Dirichlet data for 𝜅 and the Neumann data for 𝐴1 and 𝐴2. The corresponding nonhomogeneous algebraic conditions on 𝐴_𝑖𝑗 then read 2𝜕𝜕𝑛𝐴^𝑖𝑗𝑛_(𝑖𝑚_𝑗)=𝜕2𝜕𝜕𝑛𝐴1,𝜕𝑛𝐴^𝑖𝑗𝑛_(𝑖𝑙_𝑗)=𝜕2𝐴_𝑖𝑗𝜕𝑛𝐴2,1𝑙_(𝑖𝑚_𝑗)=𝐴3,2𝐴^𝑖𝑗𝑙^𝑖1𝑙_𝑗−𝑚_𝑖𝑚_𝑗=𝐴4,6𝐴^𝑖𝑗2𝑛^𝑖𝑛_𝑗−𝑙_𝑖𝑙_𝑗−𝑚_𝑖𝑚_𝑗=𝐴5.(4.23) 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 𝐴(0) and 𝜕_𝑡𝐴(0) 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 𝜖₁=12‖‖𝜕^𝑡𝐴‖‖²+‖‖𝜕_𝑙𝐴‖‖².(5.1) As in Section 4, by differentiating 𝜖₁ in time, integrating terms with spatial derivatives by parts, and using (4.3), we obtain 𝜕_𝑡𝜖₁=_𝜕Ω𝜕𝜕𝑛𝐴^𝑖𝑗𝜕^𝑡𝐴^𝑖𝑗𝑑𝜎.(5.2) 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 𝛾(0), 𝐴(0), 𝜅(0), 𝜑(0), 𝛼(0), and Γ(0) be given, and let 𝜕_𝑡𝐴(0) and 𝜕_𝑡𝜅(0) be determined by (3.6) and (3.4), respectively. If 𝜕_𝑡𝐴(0) and 𝜕_𝑡𝜅(0)∈𝐿₂(Ω) and also 𝜕_𝑙𝐴_𝑖𝑗(0) and 𝜕_𝑙𝜅(0)∈𝐿₂(Ω), then there exists unique solutions 𝐴 and 𝜅 to (4.3) and (4.4). These solutions obey the following energy estimates: sup_{0≤𝑡≤𝑇}su‖‖p_{0≤𝑡≤𝑇}𝜕^𝑡𝐴_𝑖𝑗‖‖_𝐿2(Ω)+‖‖𝜕_𝑙𝐴_𝑖𝑗‖‖_𝐿2(Ω)≤‖‖𝜕^𝑡𝐴_𝑖𝑗(0)‖‖_𝐿2(Ω)+‖‖𝜕_𝑙𝐴_𝑖𝑗(0)‖‖_𝐿2(Ω),‖‖𝜕^𝑡𝜅‖‖_𝐿2(Ω)+‖‖𝜕_𝑙𝜅‖‖_𝐿2(Ω)≤‖‖𝜕^𝑡𝜅(0)‖‖_𝐿2(Ω)+‖‖𝜕_𝑙𝜅(0)‖‖_𝐿2(Ω).(5.3)

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) 𝑁^{𝑝𝑞𝑖𝑗}𝛾_𝑝𝑞+𝛿^𝑝𝑖𝛿^𝑞𝑗−𝑁^{𝑝𝑞𝑖𝑗}𝜕𝜕𝑛𝛾^𝑝𝑞=𝑁^{𝑝𝑞𝑖𝑗}𝛾_𝑝𝑞(0)+𝛿^𝑝𝑖𝛿^𝑞𝑗−𝑁^{𝑝𝑞𝑖𝑗}𝜕𝜕𝑛𝛾^𝑝𝑞(0)+^𝑡0𝑁^{𝑝𝑞𝑖𝑗}+𝛿^𝑝𝑖𝛿^𝑞𝑗−𝑁^{𝑝𝑞𝑖𝑗}𝜕𝜕𝑛2𝜕^(𝑝𝛽_𝑞)−23𝛿^𝑝𝑞𝜕^𝑠𝛽_𝑠(6.1)𝜇𝜑+(1−𝜇)𝜕𝜕𝑛𝜑=𝜇𝜑(0)+(1−𝜇)𝜕𝜕𝑛𝜑(0)+^𝑡016𝜇−(1−𝜇)𝜕𝜕𝑛𝜕^𝑙𝛽_𝑙,(6.2)𝜇𝛼+(1−𝜇)𝜕𝜕𝑛𝛼=𝜇𝛼(0)+(1−𝜇)𝜕𝜕𝑛𝛼(0),(6.3)𝑛^𝑖Γ_𝑖=𝑛^𝑖Γ_𝑖(0)+∫^𝑡0−43𝜕𝜕𝑛𝜅+13𝜕𝜕𝑛𝜕^𝑝𝜏_𝑖𝛽_𝑝+𝜕^𝑝𝜕_𝑝𝑛^𝑖𝛽_𝑖,Γ_𝑖=𝜏^𝑖Γ_𝑖(0)+∫^𝑡0−43𝜕𝜕𝜏𝜅+13𝜕𝜕𝜏𝜕^𝑝𝛽_𝑝+𝜕^𝑝𝜕_𝑝𝜏^𝑖𝛽_𝑖.(6.4) The projection operator 𝑁^{𝑝𝑞𝑖𝑗} corresponding to (4.17) is given by 𝑁^{𝑝𝑞𝑖𝑗}=2𝑛^(𝑝𝑚^𝑞)𝑛_(𝑖𝑚_𝑗)+2𝑛^(𝑝𝑙^𝑞)𝑛_(𝑖𝑙_𝑗)(6.5) and the one corresponding to (4.18) is 𝑁^{𝑝𝑞𝑖𝑗}=2𝑙^(𝑝𝑚^𝑞)𝑙_(𝑖𝑚_𝑗)+12𝑙^𝑝+1𝑙^𝑞−𝑚^𝑝𝑚^𝑞𝑙^𝑖𝑙_𝑗−𝑚_𝑖𝑚_𝑗62𝑛^𝑝𝑛^𝑞−𝑙^𝑝𝑙^𝑞−𝑚^𝑝𝑚^𝑞2𝑛^𝑖𝑛_𝑗−𝑙_𝑖𝑙_𝑗−𝑚_𝑖𝑚_𝑗.(6.6) In (6.2) and (6.3), in the case of (4.17), we have 𝜇=0; in the case of (4.18) 𝜇=1. 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 𝜏^𝑖𝜕𝜕𝑛Γ^𝑖=𝜏^𝑖𝜕𝜕𝑛Γ^𝑖(0)+^𝑡0−43𝜕𝜕𝜏𝜕𝜕𝑛𝜅+𝜕𝜕𝑛31𝜕𝜕𝜏𝜕^𝑝𝛽_𝑝+𝜕^𝑝𝜕_𝑝𝜏^𝑖𝛽_𝑖.(6.7)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 𝐴(0), 𝜅(0), 𝜕_𝑡𝐴(0), and 𝜕_𝑡𝜅(0) and the initial data 𝛾(0), Γ(0), 𝛼(0), and 𝜑(0) be such that 𝜕_𝑡𝐴_𝑖𝑗(0)=−12𝜕^𝑙−2𝛿_𝑖𝑗𝜕_𝑙𝛾_𝑖𝑗(0)+𝜕_(𝑖Γ_𝑗)(0)−2𝜕_𝑖𝜕_𝑗𝜑(0)𝜕^𝑙𝜕_𝑙𝜑(0)−𝜕_𝑖𝜕_𝑗𝛼(0)+13𝛿^𝑖𝑗𝜕𝜕^𝑙𝜕_𝑙𝛼(0),_𝑡𝜅(0)=−𝜕^𝑙𝜕_𝑙𝛼(0).(6.8) 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 𝜑=𝜑(0)+^𝑡0−16𝜅+61𝜕^𝑙𝛼=𝛼(0)−_𝑡0𝛽_𝑙,𝛾𝜅,_𝑖𝑗=𝛾_𝑖𝑗(0)+^𝑡0−2𝐴^𝑖𝑗+2𝜕_(𝑖𝛽_𝑗)−23𝛿^𝑖𝑗Γ_𝑖𝜕^𝑙𝛽_𝑙,=Γ_𝑖(0)+^𝑡0−43𝜕^𝑖𝜅+13𝜕^𝑖𝜕^𝑝𝛽_𝑝+𝜕^𝑝𝜕_𝑝𝛽_𝑖.(6.9) 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 𝜕_𝑡𝐴_𝑖𝑗=𝜕_𝑡𝐴_𝑖𝑗(0)+^𝑡0𝜕^𝑙𝜕_𝑙𝐴_𝑖𝑗,𝜕_𝑡𝜅=𝜕_𝑡𝜅(0)+^𝑡0𝜕^𝑙𝜕_𝑙𝜅,(6.10) 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 𝜕^𝑙𝜕^𝑗𝛾_𝑙𝑗−8𝜕^𝑙𝜕_𝑙𝜑=𝜕^𝑙𝜕^𝑗𝛾_𝑙𝑗(0)−8𝜕^𝑙𝜕_𝑙𝜑(0)−^𝑡02𝜕^𝑗𝜕^𝑙𝐴_𝑗𝑙−23𝜕^𝑗𝜕_𝑙𝑘,Γ_𝑙−8𝜕^𝑙𝜕_𝑙𝜑=𝜕^𝑙Γ_𝑙(0)−8𝜕^𝑙𝜕_𝑙𝜑(0).(6.11) 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 𝜕_𝑡Γ_𝑖=0 yielding an elliptic equation for 𝛽 that needs to be solved at each time step: 13𝜕^𝑖𝜕^𝑝𝛽_𝑝+𝜕^𝑝𝜕_𝑝𝛽_𝑖−43𝜕^𝑖𝜅=0.(6.12) The boundary conditions on 𝛽 can be taken, for example, to be 𝛽_𝑖𝑚^𝑖=0,𝛽_𝑖𝑙^𝑖=0,𝜕𝜕𝑛𝛽^𝑖𝑛^𝑖=0,(6.13) or 𝜕𝜕𝑛𝛽^𝑖𝑚^𝑖=0,𝜕𝜕𝑛𝛽^𝑖𝑙^𝑖=0,𝛽_𝑖𝑛^𝑖=0.(6.14) 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 ‖𝛽_𝑖‖_𝐻1(Ω) is bounded by ‖𝜕_𝑖𝜅‖_𝐿2(Ω) 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 𝜕^2𝑡𝛽=13𝜕^𝑖𝜕^𝑝𝛽_𝑝+𝜕^𝑝𝜕_𝑝𝛽_𝑖−43𝜕^𝑖𝜅,(6.15) which corresponds to a dynamic gamma-freezing condition 𝜕_𝑡Γ_𝑖=𝜕^2𝑡𝛽_𝑖. In this case, in addition to (6.4) and (6.13), one can consider either of two sets of radiative boundary conditions 𝛽_𝑖𝑚^𝑖=0,𝛽_𝑖𝑙^𝑖=0,𝜕_𝑡𝛽_𝑖+𝜕𝜕𝑛𝛽^𝑖𝑛^𝑖𝜕=0,_𝑡𝛽_𝑖+𝜕𝜕𝑛𝛽^𝑖𝑚^𝑖=0,𝜕_𝑡𝛽_𝑖+𝜕𝜕𝑛𝛽^𝑖𝑙^𝑖=0,𝛽_𝑖𝑛^𝑖=0.(6.16) By reasoning similar to the above, an energy estimate for ‖𝛽‖_𝐻1(Ω) in terms of ‖𝜕_𝑖𝜅‖_𝐿2(Ω) can be obtained using techniques for hyperbolic equations (cf. [37–39, 42]).

After the boundary conditions for 𝛽, 𝐴, and 𝜅 are chosen, one can set Γ_𝑖−Γ(0)_𝑖=0 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 𝑎=𝑒^6𝜑𝑄,(7.1) 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 𝑄≡1 reduces to 𝛼=6𝜑. Assuming zero-shift perturbation 𝛽_𝑖=0 for further simplicity, (3.2), (3.4)–(3.7) are restated as 𝜕_𝑡𝜑=−1𝜕6𝜅,(7.2)_𝑡𝜕𝜅=−6𝜕^𝑙𝜕_𝑙𝜑,(7.3)_𝑡𝜕𝛾_𝑖𝑗=−2𝐴_𝑖𝑗,(7.4)_𝑡𝐴_𝑖𝑗=−12𝜕^𝑙𝜕𝜕_𝑙𝛾_𝑖𝑗+𝜕_(𝑖Γ_𝑗)−8𝜕_𝑖𝜕_𝑗𝜑,(7.5)_𝑡Γ_𝑖=−43𝜕^𝑖𝜅.(7.6) Taking the scalar product of (7.3) with 𝜅, integrating the result by parts, and using (7.2) to replace 𝜅 with 𝜕_𝑡𝜑, we obtain 12𝜕^𝑡‖‖𝜅‖‖²+36‖‖𝜕_𝑙𝜑‖‖²=−6_𝜕Ω𝜕𝜕𝑛𝜑𝜅.(7.7) Also, we conclude from (7.2) that 36‖𝜕_𝑡𝜑‖²=‖𝜅‖².

From (7.2) and (7.6), we observe that 𝜕_𝑡(Γ_𝑖−8𝜕_𝑖𝜑)=0, which implies that 𝜕_𝑡‖‖Γ_𝑖−8𝜕_𝑖𝜑‖‖²=0,𝜕_𝑡(Γ_𝑖−8𝜕_𝑖𝜑),𝜕^𝑙𝛾_𝑙𝑖=0.(7.8)

Next we rewrite the right-hand side of (7.5) in the divergence form: 𝜕_𝑡𝐴_𝑖𝑗=𝜕^𝑙−12𝜕^𝑙𝛾_𝑖𝑗+𝛿_𝑙(𝑖(Γ_𝑗)−8𝜕_𝑗)𝜑).(7.9) Then we take the scalar product with 𝐴_𝑖𝑗, and integrate by parts to obtain 12𝜕^𝑡‖‖𝐴‖‖²+−12𝜕^𝑙𝛾_𝑖𝑗+𝛿_𝑙(𝑖(Γ_𝑗)−8𝜕_𝑗)𝜑),𝜕_𝑙𝐴_𝑖𝑗=∫^𝜕Ω−12𝜕𝜕𝑛𝛾^𝑖𝑗+𝑛_(𝑖(Γ_𝑗)−8𝜕_𝑗)𝜑)𝐴^𝑖𝑗.(7.10) Replacing 𝐴_𝑖𝑗 with (7.4) in the second term and using 𝜕_𝑡[𝛿_𝑙(𝑖(Γ_𝑗)−8𝜕_𝑗)𝜑)]=0, we replace the last identity with 12𝜕^𝑡‖‖𝐴‖‖²+‖‖‖12𝜕^𝑙𝛾_𝑖𝑗−𝛿_𝑙(𝑖(Γ_𝑗)−8𝜕_𝑗)𝜑)‖‖‖²=∫^𝜕Ω−12𝜕𝜕𝑛𝛾^𝑖𝑗+𝑛_(𝑖(Γ_𝑗)−8𝜕_𝑗)𝜑)𝐴^𝑖𝑗(7.11)

From (7.7)–(7.11), we observe that the energy 𝜖=‖‖𝜅‖‖²+36‖‖𝜕_𝑙𝜑‖‖²+‖‖Γ_𝑙−8𝜕_𝑙𝜑‖‖²+‖‖𝐴‖‖²+‖‖‖12𝜕^𝑙𝛾_𝑗𝑖−𝛿_𝑙(𝑖(Γ_𝑗)−8𝜕_𝑗)𝜑)‖‖‖²(7.12) has its growth determined by the boundary terms 𝜕_𝑡𝜖=−6_𝜕Ω𝜕𝜕𝑛𝜑𝜅−_𝜕Ω𝜕𝜕𝑛𝛾^𝑖𝑗𝐴^𝑖𝑗+2_𝜕Ω𝑛_(𝑖(Γ_𝑗)−8𝜕_𝑗)𝜑)𝐴^𝑖𝑗.(7.13)

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 𝜕_𝑡𝜖=∫_𝜕Ω−6𝜕𝜕𝑛𝜑(0)𝜅−12𝜕𝜕𝑛𝛾3(0)𝐴3−2𝜕𝜕𝑛𝛾4(0)𝐴4−6𝜕𝜕𝑛𝛾5(0)𝐴5+4𝑛^𝑖Γ_𝑖(0)−8𝜕𝜕𝑛𝜑(0)𝐴(57.,14) where 𝛾1 – 𝛾5 are coefficients of the decomposition 𝛾_𝑖𝑗=𝛾1𝑛_(𝑖𝑚_𝑗)+𝛾2𝑛^(𝑖𝑙_𝑗)+𝛾3𝑙^(𝑖𝑚_𝑗)+𝛾4𝑙^𝑖𝑙_𝑗−𝑚_𝑖𝑚_𝑗+𝛾52𝑛^𝑖𝑛_𝑗−𝑙_𝑖𝑙_𝑗−𝑚_𝑖𝑚_𝑗.(7.15) Energy (7.12) is conserved if the initial data is chosen as to satisfy (𝜕/𝜕𝑛)𝜑(0)=0, 𝛾1(0)=𝛾2(0)=0, and (𝜕/𝜕𝑛)𝛾3(0)=(𝜕/𝜕𝑛)𝛾4(0)=(𝜕/𝜕𝑛)𝛾5(0)=0 at the boundary. Notice that Γ_𝑖(0) 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 𝜑=0, (𝜕/𝜕𝑛)𝛾1(0)=(𝜕/𝜕𝑛)𝛾2(0)=0, and 𝛾3(0)=𝛾4(0)=𝛾5(0)=0 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 𝑒𝜕^4𝜑𝜕₀𝑘=−6𝑎𝜕^𝑝𝜕_𝑝𝜑+𝐹,(7.16)₀𝜑=−16𝑎𝑘+16𝜕^𝑙𝜕𝑏^𝑙,(7.17)₀ℎ_𝑖𝑗=−2𝑎𝐴_𝑖𝑗+2ℎ_𝑙(𝑖𝜕_𝑗)𝑏^𝑙−23ℎ^𝑖𝑗𝑒𝜕_𝑙𝑏^𝑙,(7.18)_4𝜑𝜕₀𝐴_𝑖𝑗=−𝑎𝜕^𝑝12𝜕^𝑝𝜕ℎ_𝑖𝑗−ℎ_𝑝(𝑖(Γ_𝑗)−8𝜕_𝑗)𝜑)+𝐺_𝑖𝑗,(7.19)₀Γ_𝑖=−43𝑎𝜕^𝑖𝑘+𝑆_𝑖,(7.20) where 𝐹=−6𝑎𝑒ℎ^𝑝𝑞_6𝜑Γ^𝑙𝑝𝑞𝜕_𝑙𝜑−48𝑎𝜕^𝑝𝜑𝜕_𝑝𝜑−14𝑒^6𝜑𝜕^𝑝𝜑𝜕_𝑝𝑄𝐺_𝑖𝑗𝐷^𝑝𝐷_𝑝𝑄+13𝑒^4𝜑𝑎𝑘²+𝑒^4𝜑𝑎𝐴^𝑝𝑞𝐴_𝑝𝑞,=+−𝑎4𝑎ℎ𝜕^𝑝_𝑖𝑗ℎ_𝑝(𝑖Γ^𝑗)−8𝜕_𝑗)𝜑+8𝑎Γ^𝑝𝑖𝑗𝜕_𝑝𝜑−12𝑎𝜕_𝑖𝜑𝜕_𝑗𝜑+8𝜕^𝑝𝜑𝜕_𝑝𝜑−8𝑒^6𝜑𝜕^(𝑖𝜑𝜕_𝑗)𝑄−𝑒^6𝜑𝐷^𝑖𝐷_𝑗𝑄3𝑒^6𝜑ℎ_𝑖𝑗𝜕^𝑝𝜑𝜕_𝑝𝑄+13𝑒^6𝜑ℎ_𝑖𝑗𝐷^𝑝𝐷_𝑝𝑄+𝑒^4𝜑𝑊_𝑖𝑗.(7.21) Taking the scalar product of (7.16) with 𝑘 and integrating by parts on the right-hand side, we obtain _Ω𝜕⁰𝑘𝑘𝑒^4𝜑=6_Ω𝜕^𝑝𝜑𝜕_𝑞(𝑎𝑘)ℎ^𝑝𝑞−6_𝜕Ω𝑛^𝑝𝜕_𝑝𝜑𝑎𝑘+_Ω6𝜕^𝑝𝜑𝑎𝑘𝜕_𝑞ℎ^𝑝𝑞+𝐹𝑘,(7.22) 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: 12𝜕⁰‖‖𝑘‖‖^2∗+36‖‖𝜕_𝑙𝜑‖‖²=−6_𝜕Ω𝑛^𝑝𝜕_𝑝𝜑𝑎𝑘+_Ω𝐻,(7.23) where ‖𝑘‖^2∗=∫_Ω𝑘²𝑒^4𝜑, and 𝐻=×6𝜕𝑎_𝑝𝐴𝜑_𝑝𝑞[𝜕_𝑞𝜕_𝑠𝑏^𝑠+6𝜕_𝑞𝑏^𝑠𝜕^𝑠𝜑]ℎ^𝑝𝑞+36𝜕_𝑝𝜑𝜕_𝑞𝜑−𝜕^(𝑝𝑏^𝑞)+13ℎ^𝑝𝑞𝜕_𝑠𝑏^𝑠+𝑘²𝑒^4𝜑−13𝑎𝑘+13𝜕^𝑠𝑏^𝑠+6𝜕^𝑝𝜑𝑎𝑘𝜕_𝑞ℎ^𝑝𝑞+𝐹𝑘.(7.24) Also, from (7.17), it follows that 12𝜕⁰‖‖𝜑‖‖²=_Ω−16𝑎𝑘+16𝜕^𝑙𝑏^𝑙𝜑.(7.25) Next, we notice from (7.16) and (7.20) that𝜕₀Γ^𝑖−8𝜕_𝑖𝜑=8𝑎𝜕_𝑖𝜑𝑘+43𝑒^6𝜑𝜕^𝑖𝑄𝑘−43𝜕^𝑖𝜕_𝑠𝑏^𝑠−8𝜕_𝑖𝑏^𝑠𝜕^𝑠𝜑+𝑆_𝑖.(7.26) Therefore, 12𝜕⁰‖‖Γ_𝑖−8𝜕_𝑖𝜑‖‖²=_Ω𝐽,(7.27) where 𝐽=8𝑎𝜕_𝑝𝜑𝑘+43𝑒^6𝜑𝜕^𝑝𝑄𝑘−43𝜕^𝑝𝜕_𝑠𝑏^𝑠−8𝜕_𝑝𝑏^𝑠𝜕^𝑠𝜑+𝑆_𝑝Γ^𝑞−8𝜕_𝑞𝜑ℎ^𝑝𝑞+12Γ^𝑝−8𝜕_𝑝𝜑Γ_𝑞−8𝜕_𝑞𝜑2𝑎𝐴^𝑝𝑞−𝜕^(𝑝𝑏^𝑞)+23ℎ^𝑝𝑞𝜕_𝑠𝑏^𝑠.(7.28) 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: 12𝜕⁰‖‖𝐴‖‖^2∗+‖‖‖12𝜕^𝑙ℎ_𝑖𝑗−ℎ_𝑙(𝑖Γ^𝑗)−8𝜕_𝑗)𝜑‖‖‖²=−∫^𝜕Ω12𝑛^𝑝𝜕_𝑝ℎ_𝑖𝑗−𝑛^(𝑖Γ^𝑗)−8𝜕_𝑗)𝜑𝑎𝐴^𝑖𝑗+∫_Ω𝐾,(7.29) where ‖𝐴‖^2∗=∫_Ω𝐴_𝑖𝑗𝐴^𝑖𝑗𝑒^4𝜑, and 𝐾=12𝜕^𝑝×𝜕ℎ_𝑖𝑗𝑞−ℎ_𝑝(𝑖Γ^𝑗)−8𝜕_𝑗)𝜑−−ℎ2^𝑠(𝑚𝑎𝜕𝐴^𝑛)𝑏^𝑠_𝑞(𝑚−13ℎ^𝑚𝑛𝜕_𝑠𝑏^𝑠+𝜕^𝑞𝑏^𝑠𝜕^𝑠ℎ_𝑚𝑛+ℎ_𝑙𝑞𝜕_(𝑚𝑏^𝑙+ℎ_𝑠(𝑚𝜕||_𝑞||𝑏^𝑠−23ℎ^𝑞(𝑚−ℎ_𝑞(𝑚𝜕||_𝑠||𝑏^𝑠Γ^𝑛)−8𝜕_𝑛)𝜑8×𝑎𝜕𝜕^𝑛)_𝑠𝜑𝑘+34𝑒^6𝜑𝜕^𝑛)𝑄𝑘−43𝜕^𝑛)𝜕_𝑠𝑏^𝑠−8𝜕_𝑛)𝑏^𝑠𝜑+𝑆_𝑛)ℎ^𝑝𝑞ℎ^𝑖𝑚ℎ^𝑗𝑛+𝑎𝐴_𝑚𝑛𝜕_𝑞ℎ^𝑝𝑞ℎ^𝑖𝑚ℎ^𝑗𝑛+𝐺^𝑖𝑗𝐴^𝑖𝑗.(7.30)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.