Table of Contents Author Guidelines Submit a Manuscript
Abstract and Applied Analysis
VolumeΒ 2008, Article IDΒ 742040, 21 pages
http://dx.doi.org/10.1155/2008/742040
Research Article

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., [19]) 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., [1215]). 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., [35, 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]: πœ• 0 β„Ž 𝑖 𝑗 = βˆ’ 2 π‘Ž π‘˜ 𝑖 𝑗 + 2 β„Ž 𝑙 ( 𝑖 πœ• 𝑗 ) 𝑏 𝑙 , πœ• ( 2 . 1 ) 0 π‘˜ 𝑖 𝑗  𝑅 = π‘Ž 𝑖 𝑗 + ξ‚€ π‘˜ 𝑙 𝑙  π‘˜ 𝑖 𝑗 βˆ’ 2 π‘˜ 𝑖 𝑙 π‘˜ 𝑙 𝑗 ξ‚„ + π‘˜ 𝑖 𝑙 πœ• 𝑗 𝑏 𝑙 + π‘˜ 𝑙 𝑗 πœ• 𝑖 𝑏 𝑙 βˆ’ 𝐷 𝑖 𝐷 𝑗 𝑅 π‘Ž , ( 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 0 = βˆ’ π‘Ž 2 + 𝑏 𝑖 𝑏 𝑗 β„Ž 𝑖 𝑗 , 𝑔 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, πœ• 0 ∢ = ( πœ• 𝑑 βˆ’ 𝑏 𝑠 πœ• 𝑠 ) 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 𝑑 = π‘₯ 0 , π‘₯ 1 , π‘₯ 2 , π‘₯ 3 are specified. Furthermore, 𝑅 𝑖 𝑗 are the components of the spatial Ricci tensor 𝑅 𝑖 𝑗 = 1 2 β„Ž 𝑝 π‘ž ξ‚€ πœ• 𝑝 πœ• 𝑗 β„Ž 𝑖 π‘ž + πœ• 𝑖 πœ• 𝑝 β„Ž π‘ž 𝑗 βˆ’ πœ• 𝑝 πœ• π‘ž β„Ž 𝑖 𝑗 βˆ’ πœ• 𝑖 πœ• 𝑗 β„Ž 𝑝 π‘ž  + β„Ž 𝑝 π‘ž β„Ž π‘Ÿ 𝑠 ξ‚€ Ξ“ 𝑖 𝑝 𝑠 Ξ“ π‘ž 𝑗 π‘Ÿ βˆ’ Ξ“ 𝑝 π‘ž 𝑠 Ξ“ 𝑖 𝑗 π‘Ÿ  , ( 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 πœ• 0 π‘˜ 𝑖 𝑖 = π‘Ž π‘˜ 𝑙 𝑖 π‘˜ 𝑙 𝑖 βˆ’ 𝐷 𝑙 𝐷 𝑙 π‘Ž . ( 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 πœ• 0 1 π‘˜ = 3 π‘Ž π‘˜ 2 + π‘Ž 𝐴 𝑙 π‘š 𝐴 𝑙 π‘š βˆ’ 𝐷 𝑙 𝐷 𝑙 π‘Ž . ( 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 ) π‘Ž π‘˜ 2 which implies that ∫ π‘˜ β‰₯ [ ( 1 / 3 ) 𝑑 0 π‘Ž ( 𝜏 ) 𝑑 𝜏 + 1 / π‘˜ ( 0 ) ] βˆ’ 1 , 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 π‘˜ = 0 .

Alternatively, we can use harmonic slicing [2, 28] (a particular case of the Bona-Massó family of π‘˜ -driving slicing conditions ( πœ• 𝑑 βˆ’ 𝑏 𝑙 𝐷 𝑙 ) π‘Ž = βˆ’ π‘Ž 2 𝑓 ( π‘Ž ) π‘˜ , 𝑓 ( π‘Ž ) > 0 [31, 32]), which corresponds to setting πœ• 0 π‘Ž = βˆ’ π‘Ž 2 π‘˜ . ( 2 . 1 0 )

The equation on 𝐴 is obtained from (2.2), (2.8), and (2.1) as πœ• 0 𝐴 𝑖 𝑗 = π‘Ž 𝑅 𝑖 𝑗 + 1 3 π‘Ž π‘˜ 𝐴 𝑖 𝑗 βˆ’ 2 π‘Ž 𝐴 𝑖 𝑙 𝐴 𝑙 𝑗 + 2 9 π‘Ž π‘˜ 2 β„Ž 𝑖 𝑗 βˆ’ 1 3 π‘Ž 𝐴 𝑙 π‘š 𝐴 𝑙 π‘š β„Ž 𝑖 𝑗 + 𝐴 𝑖 𝑙 πœ• 𝑗 𝑏 𝑙 + 𝐴 𝑗 𝑙 πœ• 𝑖 𝑏 𝑙 βˆ’ 𝐷 𝑖 𝐷 𝑗 1 π‘Ž + 3 β„Ž 𝑖 𝑗 𝐷 𝑙 𝐷 𝑙 π‘Ž . ( 2 . 1 1 ) 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 / 1 2 ) l n ( d e t ( β„Ž 𝑖 𝑗 ) ) and the conformal metric β„Ž 𝑖 𝑗 = 𝑒 βˆ’ 4 πœ‘ β„Ž 𝑖 𝑗 , β„Ž 𝑖 𝑗 = 𝑒 4 πœ‘ β„Ž 𝑖 𝑗 . Using the Leibnitz formula for differentiating the determinant of a matrix πœ• d e t ( β„Ž 𝑖 𝑗 ) = d e t ( β„Ž 𝑖 𝑗 ) β„Ž 𝑙 π‘š πœ• β„Ž π‘š 𝑙 , ( 2 . 1 2 ) one finds that the derivative πœ• β„Ž 𝑖 𝑗 = 𝑒 βˆ’ 4 πœ‘  πœ• β„Ž 𝑖 𝑗 βˆ’ 1 3 β„Ž 𝑖 𝑗 β„Ž 𝑙 π‘š πœ• β„Ž 𝑙 π‘š ξ‚„ ( 2 . 1 3 ) of the conformal metric is trace-free. By applying operator πœ• 0 on the definition of πœ‘ and using (2.12) and (2.1), we get the second equation of our system πœ• 0 1 πœ‘ = βˆ’ 6 1 π‘Ž π‘˜ + 6 πœ• 𝑙 𝑏 𝑙 . ( 2 . 1 4 ) Now using (2.13) and (2.1), we obtain the third equation πœ• 0 β„Ž i j 𝐴 = βˆ’ 2 π‘Ž 𝑖 𝑗 β„Ž + 2 𝑙 ( 𝑖 πœ• 𝑗 ) 𝑏 𝑙 βˆ’ 2 3 β„Ž 𝑖 𝑗 πœ• 𝑙 𝑏 𝑙 , ( 2 . 1 5 ) 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 πœ• 0 = πœ• 𝑑 βˆ’ 𝑏 𝑠 πœ• 𝑠 .

The remaining two equations can be obtained from (2.11) which can be rewritten in terms of 𝐴 as πœ• 0 𝐴 𝑖 𝑗 = π‘Ž 𝑒 βˆ’ 4 πœ‘ 𝑅 𝑖 𝑗 ξ‚€ π‘˜ 𝐴 + π‘Ž 𝑖 𝑗 𝐴 βˆ’ 2 𝑖 𝑙 𝐴 𝑙 𝑗 + 2 9 π‘˜ 2 β„Ž 𝑖 𝑗 βˆ’ 1 3 𝐴 𝑙 π‘š 𝐴 𝑙 π‘š β„Ž 𝑖 𝑗  + 𝐴 𝑖 𝑙 πœ• 𝑗 𝑏 𝑙 + 𝐴 𝑗 𝑙 πœ• 𝑖 𝑏 𝑙 βˆ’ 2 3 𝐴 𝑖 𝑗 πœ• 𝑙 𝑏 𝑙 βˆ’ 𝑒 βˆ’ 4 πœ‘ 𝐷 𝑖 𝐷 𝑗 π‘Ž + 𝑒 βˆ’ 4 πœ‘ 1 3 β„Ž 𝑖 𝑗 𝐷 𝑙 𝐷 𝑙 π‘Ž . ( 2 . 1 6 ) The Ricci tensor in terms of the conformal metric becomes [2, 28] 𝑅 𝑖 𝑗 = 1 2 β„Ž 𝑝 π‘ž ξ‚€ πœ• 𝑝 πœ• 𝑗 β„Ž 𝑖 π‘ž + πœ• 𝑖 πœ• 𝑝 β„Ž π‘ž 𝑗 βˆ’ πœ• 𝑝 πœ• π‘ž β„Ž 𝑖 𝑗 βˆ’ πœ• 𝑖 πœ• 𝑗 β„Ž 𝑝 π‘ž  𝐷 βˆ’ 2 𝑖 𝐷 𝑗 β„Ž πœ‘ βˆ’ 2 𝑖 𝑗 β„Ž 𝑝 π‘ž 𝐷 𝑝 𝐷 π‘ž πœ‘ + β„Ž 𝑝 π‘ž β„Ž π‘Ÿ 𝑠 ξ‚€ Ξ“ 𝑖 𝑝 𝑠 Ξ“ π‘ž 𝑗 π‘Ÿ βˆ’ Ξ“ 𝑝 π‘ž 𝑠 Ξ“ 𝑖 𝑗 π‘Ÿ  + 4 πœ• 𝑖 πœ‘ πœ• 𝑗 β„Ž πœ‘ βˆ’ 4 𝑖 𝑗 β„Ž 𝑝 π‘ž πœ• 𝑝 πœ‘ πœ• π‘ž πœ‘ . ( 2 . 1 7 ) Here Ξ“ 𝑖 𝑗 π‘˜ = ( πœ• 𝑖 β„Ž π‘˜ 𝑗 + πœ• 𝑗 β„Ž 𝑖 π‘˜ βˆ’ πœ• π‘˜ β„Ž 𝑗 𝑖 ) / 2 and 𝐷 𝑖 𝑣 𝑗 = πœ• 𝑖 𝑣 𝑗 βˆ’ β„Ž 𝑝 π‘ž Ξ“ 𝑖 𝑗 𝑝 𝑣 π‘ž are the covariant derivatives associated with the conformal metric. The first line in (2.17) can be rewritten to obtain 𝑅 𝑖 𝑗 1 = βˆ’ 2 β„Ž 𝑝 π‘ž πœ• 𝑝 πœ• π‘ž β„Ž 𝑖 𝑗 + πœ• ( 𝑖 β„Ž 𝑝 π‘ž Ξ“ | | 𝑝 π‘ž | | 𝑗 ) + Ξ“ 𝑝 π‘ž ( 𝑖 πœ• 𝑗 ) β„Ž 𝑝 π‘ž 𝐷 βˆ’ 2 𝑖 𝐷 𝑗 β„Ž πœ‘ βˆ’ 2 𝑖 𝑗 β„Ž 𝑝 π‘ž 𝐷 𝑝 𝐷 π‘ž πœ‘ + β„Ž 𝑝 π‘ž β„Ž π‘Ÿ 𝑠 ξ‚€ Ξ“ 𝑖 𝑝 𝑠 Ξ“ π‘ž 𝑗 π‘Ÿ βˆ’ Ξ“ 𝑝 π‘ž 𝑠 Ξ“ 𝑖 𝑗 π‘Ÿ  + 4 πœ• 𝑖 πœ‘ πœ• 𝑗 β„Ž πœ‘ βˆ’ 4 𝑖 𝑗 β„Ž 𝑝 π‘ž πœ• 𝑝 πœ‘ πœ• π‘ž πœ‘ . ( 2 . 1 8 ) This suggests that we introduce a new variable Ξ“ 𝑗 = β„Ž 𝑝 π‘ž Ξ“ 𝑝 π‘ž 𝑗 = β„Ž 𝑝 π‘ž πœ• 𝑝 β„Ž π‘ž 𝑗 . ( 2 . 1 9 ) Substituting (2.18) in (2.16), we derive the fourth evolution equation πœ• 0 𝐴 𝑖 𝑗 1 = βˆ’ 2 π‘Ž 𝑒 βˆ’ 4 πœ‘ β„Ž 𝑝 π‘ž πœ• 𝑝 πœ• π‘ž β„Ž 𝑖 𝑗 + π‘Ž 𝑒 βˆ’ 4 πœ‘ πœ• ( 𝑖 Ξ“ 𝑗 ) βˆ’ 2 π‘Ž 𝑒 βˆ’ 4 πœ‘ 𝐷 𝑖 𝐷 𝑗 πœ‘ βˆ’ 2 π‘Ž 𝑒 βˆ’ 4 πœ‘ β„Ž 𝑖 𝑗 β„Ž 𝑝 π‘ž 𝐷 𝑝 𝐷 π‘ž πœ‘ βˆ’ 𝑒 βˆ’ 4 πœ‘ 𝐷 𝑖 𝐷 𝑗 1 π‘Ž + 3 𝑒 βˆ’ 4 πœ‘ β„Ž 𝑖 𝑗 𝐷 𝑙 𝐷 𝑙 π‘Ž + π‘Š 𝑖 𝑗 , ( 2 . 2 0 ) where π‘Š 𝑖 𝑗 = π‘Ž 𝑒 βˆ’ 4 πœ‘ Ξ“ 𝑝 π‘ž ( 𝑖 πœ• 𝑗 ) β„Ž 𝑝 π‘ž + π‘Ž 𝑒 βˆ’ 4 πœ‘ β„Ž 𝑝 π‘ž β„Ž π‘Ÿ 𝑠 ξ‚€ Ξ“ 𝑖 𝑝 𝑠 Ξ“ π‘ž 𝑗 π‘Ÿ βˆ’ Ξ“ 𝑝 π‘ž 𝑠 Ξ“ 𝑖 𝑗 π‘Ÿ  + 4 π‘Ž 𝑒 βˆ’ 4 πœ‘ πœ• 𝑖 πœ‘ πœ• 𝑗 πœ‘ βˆ’ 4 π‘Ž 𝑒 βˆ’ 4 πœ‘ β„Ž 𝑖 𝑗 β„Ž 𝑝 π‘ž πœ• 𝑝 πœ‘ πœ• π‘ž πœ‘ ξ‚€ π‘˜ 𝐴 + π‘Ž 𝑖 𝑗 𝐴 βˆ’ 2 𝑖 𝑙 𝐴 𝑙 𝑗 + 2 9 π‘˜ 2 β„Ž 𝑖 𝑗 βˆ’ 1 3 𝐴 𝑙 π‘š 𝐴 𝑙 π‘š β„Ž 𝑖 𝑗  + 𝐴 𝑖 𝑙 πœ• 𝑗 𝑏 𝑙 + 𝐴 𝑗 𝑙 πœ• 𝑖 𝑏 𝑙 βˆ’ 2 3 𝐴 𝑖 𝑗 πœ• 𝑙 𝑏 𝑙 . ( 2 . 2 1 )

The evolution equation for Ξ“ 𝑗 is obtained by applying the operator πœ• 0 to (2.19) and using the momentum constraint πœ• 0 Ξ“ 𝑖 = βˆ’ 2 π‘Ž πœ• 𝑝 𝐴 𝑝 𝑖 β„Ž βˆ’ 2 𝑝 π‘ž ξ‚€ πœ• 𝑝 π‘Ž  𝐴 π‘ž 𝑖 𝐴 + 2 π‘Ž 𝑝 π‘ž ξ‚€ πœ• 𝑝 β„Ž π‘ž 𝑖  + Ξ“ 𝑙 πœ• 𝑖 𝑏 𝑙 + 1 3 πœ• 𝑖 πœ• 𝑙 𝑏 𝑙 + β„Ž 𝑙 𝑖 πœ• 𝑠 πœ• 𝑠 𝑏 𝑙 . ( 2 . 2 2 ) Next we notice that β„Ž 𝑝 π‘ž 𝐷 𝑝 𝐴 𝑖 π‘ž = πœ• 𝑝 𝐴 𝑝 𝑖 βˆ’ Ξ“ 𝑠 𝐴 𝑠 𝑖 + 6 ( πœ• 𝑠 𝐴 πœ‘ ) 𝑠 𝑖 , and thus (2.4) takes the form πœ• 𝑝 𝐴 𝑝 𝑖 βˆ’ 2 3 πœ• 𝑖 π‘˜ βˆ’ Ξ“ 𝑠 𝐴 𝑠 𝑖 + 6 ( πœ• 𝑠 𝐴 πœ‘ ) 𝑠 𝑖 = 0 . ( 2 . 2 3 ) Solving this equation for πœ• 𝑝 𝐴 𝑝 𝑖 and substituting the result in (2.22), we derive the fifth equation of the BSSN system πœ• 0 Ξ“ 𝑖 4 = βˆ’ 3 π‘Ž πœ• 𝑖 π‘˜ + 𝑆 𝑖 , ( 2 . 2 4 ) where 𝑆 𝑖 = βˆ’ 2 π‘Ž Ξ“ 𝑠 𝐴 𝑠 𝑖 + 1 2 π‘Ž ( πœ• 𝑠 𝐴 πœ‘ ) 𝑠 𝑖 β„Ž βˆ’ 2 𝑝 π‘ž ξ‚€ πœ• 𝑝 π‘Ž  𝐴 π‘ž 𝑗 𝐴 + 2 π‘Ž 𝑝 π‘ž ξ‚€ πœ• 𝑝 β„Ž π‘ž 𝑗  + Ξ“ 𝑙 πœ• 𝑗 𝑏 𝑙 + 1 3 πœ• 𝑖 πœ• 𝑙 𝑏 𝑙 + β„Ž 𝑙 𝑗 πœ• 𝑠 πœ• 𝑠 𝑏 𝑙 . ( 2 . 2 5 ) 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 d e t ( β„Ž 𝑖 𝑗 ) = 1 + 𝛾 𝑙 𝑙 1 , πœ‘ = 𝛾 1 2 𝑙 𝑙 , 𝑒 βˆ’ 4 πœ‘ 1 = 1 βˆ’ 3 𝛾 𝑙 𝑙 , 𝑒 4 πœ‘ 1 = 1 + 3 𝛾 𝑙 𝑙 , β„Ž 𝑖 𝑗 = 𝛿 𝑖 𝑗 + 𝛾 𝑖 𝑗 βˆ’ 1 3 𝛿 𝑖 𝑗 𝛾 𝑙 𝑙 = ∢ 𝛿 𝑖 𝑗 + 𝛾 𝑖 𝑗 , π‘˜ = πœ… = ∢ πœ… 𝑙 𝑙 , 𝐴 𝑖 𝑗 = 𝐴 𝑖 𝑗 = πœ… 𝑖 𝑗 βˆ’ 1 3 𝛿 𝑖 𝑗 πœ… , Ξ“ 𝑖 = ∢ Ξ“ 𝑖 = πœ• 𝑙 𝛾 𝑖 𝑙 . ( 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 πœ• 𝑑 1 πœ‘ = βˆ’ 6 1 πœ… + 6 πœ• 𝑙 𝛽 𝑙 , πœ• ( 3 . 2 ) 𝑑 πœ• 𝛼 = βˆ’ πœ… , ( 3 . 3 ) 𝑑 πœ… = βˆ’ πœ• 𝑙 πœ• 𝑙 πœ• 𝛼 , ( 3 . 4 ) 𝑑 𝛾 𝑖 𝑗 = βˆ’ 2 𝐴 𝑖 𝑗 + 2 πœ• ( 𝑖 𝛽 𝑗 ) βˆ’ 2 3 𝛿 𝑖 𝑗 πœ• 𝑙 𝛽 𝑙 , πœ• ( 3 . 5 ) 𝑑 𝐴 𝑖 𝑗 1 = βˆ’ 2 πœ• 𝑙 πœ• 𝑙 𝛾 𝑖 𝑗 + πœ• ( 𝑖 Ξ“ 𝑗 ) βˆ’ 2 πœ• 𝑖 πœ• 𝑗 πœ‘ βˆ’ 2 𝛿 𝑖 𝑗 πœ• 𝑙 πœ• 𝑙 πœ‘ βˆ’ πœ• 𝑖 πœ• 𝑗 1 𝛼 + 3 𝛿 𝑖 𝑗 πœ• 𝑙 πœ• 𝑙 πœ• 𝛼 , ( 3 . 6 ) 𝑑 Ξ“ 𝑖 4 = βˆ’ 3 πœ• 𝑖 1 πœ… + 3 πœ• 𝑖 πœ• 𝑝 𝛽 𝑝 + πœ• 𝑝 πœ• 𝑝 𝛽 𝑖 . ( 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 ) 𝑙 𝐴 𝑖 𝑙 βˆ’ 2 3 πœ• 𝑖 πœ… = 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 . 1 0 )

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 πœ• 𝑙 πœ• 𝑙 πœ‘  = πœ• 𝑖 ξ‚€ πœ• 𝑙 𝐴 𝑖 𝑙 βˆ’ 2 3 πœ• 𝑖 πœ…  , πœ• 𝑑 ξ‚€ πœ• 𝑙 Ξ“ 𝑙 βˆ’ 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 𝑀 𝑖 = πœ• 𝑙 𝐴 𝑖 𝑙 βˆ’ 2 3 πœ• 𝑖 πœ… . ( 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 πœ• 𝑑 𝑀 𝑖 1 = βˆ’ 2 πœ• 𝑙 πœ• 𝑙 πœ• π‘š 𝛾 𝑖 π‘š + 1 2 πœ• 𝑖 πœ• 𝑙 Ξ“ 𝑙 + 1 2 πœ• 𝑙 πœ• 𝑙 Ξ“ 𝑖 βˆ’ 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, 1618, 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 𝐿 2 norms naturally associated with these scalar products are β€– 𝑒 β€– 2 = ( 𝑒 𝑖 , 𝑒 𝑖 ) and β€– 𝜌 β€– 2 = ( 𝜌 𝑖 𝑗 , 𝜌 𝑖 𝑗 ) . We introduce the energy function for the constraint quantity β€– β€– πœ• πœ– = 𝑑 𝑀 β€– β€– 2 + β€– β€– πœ• 𝑙 𝑀 β€– β€– 2 . ( 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 ξ‚€ πœ• 𝑀 πœ• 𝑛 𝑖 πœ•  ξ‚€ 𝑑 𝑀 𝑖  ≀ 0 o n πœ• Ξ© . ( 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 . 1 0 ) Similarly, the gradient of a scalar field πœ“ becomes πœ• 𝑖 πœ“ = 𝑛 𝑖 πœ• πœ• 𝑛 πœ“ + π‘š 𝑖 πœ• πœ• π‘š πœ“ + 𝑙 𝑖 πœ• πœ• 𝑙 πœ“ . ( 4 . 1 1 ) Next we note that at any point of the boundary, a symmetric trace free matrix is spanned by 𝑛 ( 𝑖 π‘š 𝑗 ) , 𝑛 ( 𝑖 𝑙 𝑗 ) , 𝑙 ( 𝑖 π‘š 𝑗 ) , 𝑙 𝑖 𝑙 𝑗 βˆ’ π‘š 𝑖 π‘š 𝑗 , 2 𝑛 𝑖 𝑛 𝑗 βˆ’ 𝑙 𝑖 𝑙 𝑗 βˆ’ π‘š 𝑖 π‘š 𝑗 . ( 4 . 1 2 ) Introducing the scalar functions 𝐴 1 = 2 𝐴 𝑖 𝑗 𝑛 ( 𝑖 π‘š 𝑗 ) , 𝐴 2 = 2 𝐴 𝑖 𝑗 𝑛 ( 𝑖 𝑙 𝑗 ) , 𝐴 3 = 2 𝐴 𝑖 𝑗 𝑙 ( 𝑖 π‘š 𝑗 ) , 1 𝐴 4 = 2 𝐴 𝑖 𝑗 ξ‚€ 𝑙 𝑖 𝑙 𝑗 βˆ’ π‘š 𝑖 π‘š 𝑗  1 , 𝐴 5 = 6 𝐴 𝑖 𝑗 ξ‚€ 2 𝑛 𝑖 𝑛 𝑗 βˆ’ 𝑙 𝑖 𝑙 𝑗 βˆ’ π‘š 𝑖 π‘š 𝑗  , ( 4 . 1 3 ) we rewrite 𝐴 as follows: 𝐴 𝑖 𝑗 ξ‚€ 𝑛 = 𝐴 1 ( 𝑖 π‘š 𝑗 )  ξ‚€ 𝑛 + 𝐴 2 ( 𝑖 𝑙 𝑗 )  ξ‚€ 𝑙 + 𝐴 3 ( 𝑖 π‘š 𝑗 )  ξ‚€ 𝑙 + 𝐴 4 𝑖 𝑙 𝑗 βˆ’ π‘š 𝑖 π‘š 𝑗  ξ‚€ + 𝐴 5 2 𝑛 𝑖 𝑛 𝑗 βˆ’ 𝑙 𝑖 𝑙 𝑗 βˆ’ π‘š 𝑖 π‘š 𝑗  ( 4 . 1 4 ) Substituting (4.14) into (4.2) and using (4.10) and (4.11) to replace partial derivatives with directional derivatives, we obtain 𝑀 𝑖 =  1 2 πœ• 1 πœ• π‘š 𝐴 1 + 2 πœ• πœ• πœ• 𝑙 𝐴 2 + 2 2 πœ• 𝑛 𝐴 5 βˆ’ 3 πœ• πœ… ξ‚„ 𝑛 πœ• 𝑛 𝑖 +  1 2 πœ• 1 πœ• 𝑛 𝐴 1 + 2 πœ• πœ• πœ• 𝑙 𝐴 3 βˆ’ πœ• πœ• π‘š 𝐴 4 βˆ’ 2 πœ• π‘š 𝐴 5 βˆ’ 3 πœ• πœ… ξ‚„ π‘š πœ• π‘š 𝑖 +  1 2 πœ• 1 πœ• 𝑛 𝐴 2 + 2 πœ• πœ• πœ• π‘š 𝐴 3 + πœ• πœ• 𝑙 𝐴 4 βˆ’ 2 πœ• 𝑙 𝐴 5 βˆ’ 3 πœ• πœ… ξ‚„ 𝑙 πœ• 𝑙 𝑖 . ( 4 . 1 5 ) The last equation implies that ξ‚€ πœ• 𝑀 πœ• 𝑛 𝑖 πœ•  ξ‚€ 𝑑 𝑀 𝑖  = πœ•  1 πœ• 𝑛 2 πœ• 1 πœ• π‘š 𝐴 1 + 2 πœ• πœ• πœ• 𝑙 𝐴 2 + 2 2 πœ• 𝑛 𝐴 5 βˆ’ 3 πœ• πœ… ξ‚„ πœ• 𝑛 Γ— πœ• 𝑑  1 2 πœ• 1 πœ• π‘š 𝐴 1 + 2 πœ• πœ• πœ• 𝑙 𝐴 2 + 2 2 πœ• 𝑛 𝐴 5 βˆ’ 3 πœ• πœ… ξ‚„ + πœ• πœ• 𝑛  1 πœ• 𝑛 2 πœ• 1 πœ• 𝑛 𝐴 1 + 2 πœ• πœ• πœ• 𝑙 𝐴 3 βˆ’ πœ• πœ• π‘š 𝐴 4 βˆ’ 2 πœ• π‘š 𝐴 5 βˆ’ 3 πœ• πœ… ξ‚„ πœ• π‘š Γ— πœ• 𝑑  1 2 πœ• 1 πœ• 𝑛 𝐴 1 + 2 πœ• πœ• πœ• 𝑙 𝐴 3 βˆ’ πœ• πœ• π‘š 𝐴 4 βˆ’ 2 πœ• π‘š 𝐴 5 βˆ’ 3 πœ• πœ… ξ‚„ + πœ• πœ• π‘š  1 πœ• 𝑛 2 πœ• 1 πœ• 𝑛 𝐴 2 + 2 πœ• πœ• πœ• π‘š 𝐴 3 + πœ• πœ• 𝑙 𝐴 4 βˆ’ 2 πœ• 𝑙 𝐴 5 βˆ’ 3 πœ• πœ… ξ‚„ πœ• 𝑙 Γ— πœ• 𝑑  1 2 πœ• 1 πœ• 𝑛 𝐴 2 + 2 πœ• πœ• πœ• π‘š 𝐴 3 + πœ• πœ• 𝑙 𝐴 4 βˆ’ 2 πœ• 𝑙 𝐴 5 βˆ’ 3 πœ• πœ… ξ‚„ . πœ• 𝑙 ( 4 . 1 6 ) 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 . 1 7 ) πœ• πœ• 𝑛 𝐴 1 = 0 , πœ• 𝑛 𝐴 2 = 0 , 𝐴 3 = 0 , 𝐴 4 = 0 , 𝐴 5 = 0 , πœ… = 0 . ( 4 . 1 8 ) 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 . 1 9 ) 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]: 1 2 πœ• 1 πœ• π‘š 𝐴 1 + 2 πœ• πœ• πœ• 𝑙 𝐴 2 + 2 2 πœ• 𝑛 𝐴 5 βˆ’ 3 πœ• 1 πœ• 𝑛 πœ… = 0 , 2 πœ• 1 πœ• 𝑛 𝐴 1 + 2 πœ• πœ• πœ• 𝑙 𝐴 3 βˆ’ πœ• πœ• π‘š 𝐴 4 βˆ’ 2 πœ• π‘š 𝐴 5 βˆ’ 3 πœ• 1 πœ• π‘š πœ… = 0 , 2 πœ• 1 πœ• 𝑛 𝐴 2 + 2 πœ• πœ• πœ• π‘š 𝐴 3 + πœ• πœ• 𝑙 𝐴 4 βˆ’ 2 πœ• 𝑙 𝐴 5 βˆ’ 3 πœ• πœ• 𝑙 πœ… = 0 . ( 4 . 2 0 ) 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 . 2 1 ) 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 ( πœ• 2 / πœ• 𝑛 2 ) derivatives) an evolution equation defined on the boundary 2 πœ• 2 πœ• 𝑑 2 ξ‚€ πœ• 𝐴 5 βˆ’ 2 πœ• 𝑙 2 + πœ• 2 πœ• π‘š 2  2 𝐴 5 βˆ’ 3 πœ• 2 πœ• 𝑑 2 4 πœ… + 3 ξ‚€ πœ• 2 πœ• 𝑙 2 + πœ• 2 πœ• π‘š 2  ξ‚€ πœ• πœ… = 2 πœ• 𝑙 2 βˆ’ πœ• 2 πœ• π‘š 2  πœ• 𝐴 4 + πœ• πœ• 𝑙 πœ• π‘š 𝐴 3 . ( 4 . 2 2 ) 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 . 2 3 ) 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 𝐿 2 norms naturally associated with these scalar products are β€– πœ‡ β€– 2 = ( πœ‡ , πœ‡ ) and β€– 𝑒 β€– 2 = ( 𝑒 𝑖 𝑗 π‘˜ , 𝑒 𝑖 𝑗 π‘˜ ) . Consider the system (4.3) describing the propagation of 𝐴 . We define the energy of the system to be πœ– 1 = 1 2 ξ‚€ β€– β€– πœ• 𝑑 𝐴 β€– β€– 2 + β€– β€– πœ• 𝑙 𝐴 β€– β€– 2  . ( 5 . 1 ) As in Section 4, by differentiating πœ– 1 in time, integrating terms with spatial derivatives by parts, and using (4.3), we obtain πœ• 𝑑 πœ– 1 = ξ€œ πœ• Ξ© ξ‚€ πœ• 𝐴 πœ• 𝑛 𝑖 𝑗 πœ•  ξ‚€ 𝑑 𝐴 𝑖 𝑗  𝑑 𝜎 . ( 5 . 2 ) Since the right-hand side of (5.2) is zero if either (4.17) or (4.18) holds, we conclude that πœ– 1 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 [3639] 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 ) ∈ 𝐿 2 ( Ξ© ) and also πœ• 𝑙 𝐴 𝑖 𝑗 ( 0 ) and πœ• 𝑙 πœ… ( 0 ) ∈ 𝐿 2 ( Ξ© ) , then there exists unique solutions 𝐴 and πœ… to (4.3) and (4.4). These solutions obey the following energy estimates: s u p 0 ≀ 𝑑 ≀ 𝑇  β€– β€– πœ• 𝑑 𝐴 𝑖 𝑗 β€– β€– 𝐿 2 ( Ξ© ) + β€– β€– πœ• 𝑙 𝐴 𝑖 𝑗 β€– β€– 𝐿 2 ( Ξ© ) ξ‚„ ≀ β€– β€– πœ• 𝑑 𝐴 𝑖 𝑗 β€– β€– ( 0 ) 𝐿 2 ( Ξ© ) + β€– β€– πœ• 𝑙 𝐴 𝑖 𝑗 β€– β€– ( 0 ) 𝐿 2 ( Ξ© ) , s u p 0 ≀ 𝑑 ≀ 𝑇  β€– β€– πœ• 𝑑 πœ… β€– β€– 𝐿 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 πœ• ( 𝑝 𝛽 π‘ž ) βˆ’ 2 3 𝛿 𝑝 π‘ž πœ• 𝑠 𝛽 𝑠  ( 6 . 1 ) πœ• πœ‡ πœ‘ + ( 1 βˆ’ πœ‡ ) πœ• πœ• 𝑛 πœ‘ = πœ‡ πœ‘ ( 0 ) + ( 1 βˆ’ πœ‡ ) ξ€œ πœ• 𝑛 πœ‘ ( 0 ) + 𝑑 0 1 6 ξ‚€ πœ• πœ‡ βˆ’ ( 1 βˆ’ πœ‡ )  πœ• πœ• 𝑛 𝑙 𝛽 𝑙 , ( 6 . 2 ) πœ• πœ‡ 𝛼 + ( 1 βˆ’ πœ‡ ) πœ• πœ• 𝑛 𝛼 = πœ‡ 𝛼 ( 0 ) + ( 1 βˆ’ πœ‡ ) πœ• 𝑛 𝛼 ( 0 ) , ( 6 . 3 ) 𝑛 𝑖 Ξ“ 𝑖 = 𝑛 𝑖 Ξ“ 𝑖 ∫ ( 0 ) + 𝑑 0  βˆ’ 4 3 πœ• 1 πœ• 𝑛 πœ… + 3 πœ• πœ• πœ• 𝑛 𝑝 𝛽 𝑝 + πœ• 𝑝 πœ• 𝑝 𝑛 𝑖 𝛽 𝑖 ξ‚„ , 𝜏 𝑖 Ξ“ 𝑖 = 𝜏 𝑖 Ξ“ 𝑖 ∫ ( 0 ) + 𝑑 0  βˆ’ 4 3 πœ• 1 πœ• 𝜏 πœ… + 3 πœ• πœ• πœ• 𝜏 𝑝 𝛽 𝑝 + πœ• 𝑝 πœ• 𝑝 𝜏 𝑖 𝛽 𝑖 ξ‚„ . ( 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 𝑙 ( 𝑝 π‘š π‘ž ) 𝑙 ( 𝑖 π‘š 𝑗 ) + 1 2 ξ‚€ 𝑙 𝑝 𝑙 π‘ž βˆ’ π‘š 𝑝 π‘š π‘ž 𝑙  ξ‚€ 𝑖 𝑙 𝑗 βˆ’ π‘š 𝑖 π‘š 𝑗  + 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  βˆ’ 4 3 πœ• πœ• πœ• 𝜏 πœ• πœ• 𝑛 πœ… + ξ‚€ 1 πœ• 𝑛 3 πœ• πœ• πœ• 𝜏 𝑝 𝛽 𝑝 + πœ• 𝑝 πœ• 𝑝 𝜏 𝑖 𝛽 𝑖 .  ξ‚„ ( 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 πœ• 𝑑 𝐴 𝑖 𝑗 1 ( 0 ) = βˆ’ 2 πœ• 𝑙 πœ• 𝑙 𝛾 𝑖 𝑗 ( 0 ) + πœ• ( 𝑖 Ξ“ 𝑗 ) ( 0 ) βˆ’ 2 πœ• 𝑖 πœ• 𝑗 πœ‘ ( 0 ) βˆ’ 2 𝛿 𝑖 𝑗 πœ• 𝑙 πœ• 𝑙 πœ‘ ( 0 ) βˆ’ πœ• 𝑖 πœ• 𝑗 1 𝛼 ( 0 ) + 3 𝛿 𝑖 𝑗 πœ• 𝑙 πœ• 𝑙 πœ• 𝛼 ( 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  βˆ’ 1 6 1 πœ… + 6 πœ• 𝑙 𝛽 𝑙 ξ‚„ , ξ€œ 𝛼 = 𝛼 ( 0 ) βˆ’ 𝑑 0 𝛾 πœ… , 𝑖 𝑗 = 𝛾 𝑖 𝑗 ξ€œ ( 0 ) + 𝑑 0  βˆ’ 2 𝐴 𝑖 𝑗 + 2 πœ• ( 𝑖 𝛽 𝑗 ) βˆ’ 2 3 𝛿 𝑖 𝑗 πœ• 𝑙 𝛽 𝑙 ξ‚„ , Ξ“ 𝑖 = Ξ“ 𝑖 ξ€œ ( 0 ) + 𝑑 0  βˆ’ 4 3 πœ• 𝑖 1 πœ… + 3 πœ• 𝑖 πœ• 𝑝 𝛽 𝑝 + πœ• 𝑝 πœ• 𝑝 𝛽 𝑖 ξ‚„ . ( 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 . 1 0 ) 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 ) βˆ’ 𝑑 0 2 πœ• 𝑗 ξ‚€ πœ• 𝑙 𝐴 𝑗 𝑙 βˆ’ 2 3 πœ• 𝑗 π‘˜  , πœ• 𝑙 Ξ“ 𝑙 βˆ’ 8 πœ• 𝑙 πœ• 𝑙 πœ‘ = πœ• 𝑙 Ξ“ 𝑙 ( 0 ) βˆ’ 8 πœ• 𝑙 πœ• 𝑙 πœ‘ ( 0 ) . ( 6 . 1 1 ) 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: 1 3 πœ• 𝑖 πœ• 𝑝 𝛽 𝑝 + πœ• 𝑝 πœ• 𝑝 𝛽 𝑖 βˆ’ 4 3 πœ• 𝑖 πœ… = 0 . ( 6 . 1 2 ) The boundary conditions on 𝛽 can be taken, for example, to be 𝛽 𝑖 π‘š 𝑖 = 0 , 𝛽 𝑖 𝑙 𝑖 πœ• = 0 , 𝛽 πœ• 𝑛 𝑖 𝑛 𝑖 = 0 , ( 6 . 1 3 ) or πœ• 𝛽 πœ• 𝑛 𝑖 π‘š 𝑖 πœ• = 0 , 𝛽 πœ• 𝑛 𝑖 𝑙 𝑖 = 0 , 𝛽 𝑖 𝑛 𝑖 = 0 . ( 6 . 1 4 ) 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 𝑑 1 𝛽 = 3 πœ• 𝑖 πœ• 𝑝 𝛽 𝑝 + πœ• 𝑝 πœ• 𝑝 𝛽 𝑖 βˆ’ 4 3 πœ• 𝑖 πœ… , ( 6 . 1 5 ) 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 . 1 6 ) By reasoning similar to the above, an energy estimate for β€– 𝛽 β€– 𝐻 1 ( Ξ© ) in terms of β€– πœ• 𝑖 πœ… β€– 𝐿 2 ( Ξ© ) can be obtained using techniques for hyperbolic equations (cf. [3739, 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 46. 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 46. 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 ) 𝑑 𝐴 𝑖 𝑗 1 = βˆ’ 2 πœ• 𝑙 πœ• 𝑙 𝛾 𝑖 𝑗 + πœ• ( 𝑖 Ξ“ 𝑗 ) βˆ’ 8 πœ• 𝑖 πœ• 𝑗 πœ• πœ‘ , ( 7 . 5 ) 𝑑 Ξ“ 𝑖 4 = βˆ’ 3 πœ• 𝑖 πœ… . ( 7 . 6 ) Taking the scalar product of (7.3) with πœ… , integrating the result by parts, and using (7.2) to replace πœ… with πœ• 𝑑 πœ‘ , we obtain 1 2 πœ• 𝑑  β€– β€– πœ… β€– β€– 2 β€– β€– πœ• + 3 6 𝑙 πœ‘ β€– β€– 2 ξ‚„ ξ€œ = βˆ’ 6 πœ• Ξ© ξ‚€ πœ• πœ‘  πœ• 𝑛 πœ… . ( 7 . 7 ) Also, we conclude from (7.2) that 3 6 β€– πœ• 𝑑 πœ‘ β€– 2 = β€– πœ… β€– 2 .

From (7.2) and (7.6), we observe that πœ• 𝑑 ( Ξ“ 𝑖 βˆ’ 8 πœ• 𝑖 πœ‘ ) = 0 , which implies that πœ• 𝑑 β€– β€– Ξ“ 𝑖 βˆ’ 8 πœ• 𝑖 πœ‘ β€– β€– 2 ξ‚€ πœ• = 0 , 𝑑 ( Ξ“ 𝑖 βˆ’ 8 πœ• 𝑖 πœ‘ ) , πœ• 𝑙 𝛾 𝑙 𝑖  = 0 . ( 7 . 8 )

Next we rewrite the right-hand side of (7.5) in the divergence form: πœ• 𝑑 𝐴 𝑖 𝑗 = πœ• 𝑙  βˆ’ 1 2 πœ• 𝑙 𝛾 𝑖 𝑗 + 𝛿 𝑙 ( 𝑖 ( Ξ“ 𝑗 ) βˆ’ 8 πœ• 𝑗 ) ξ‚„ πœ‘ ) . ( 7 . 9 ) Then we take the scalar product with 𝐴 𝑖 𝑗 , and integrate by parts to obtain 1 2 πœ• 𝑑 β€– β€– 𝐴 β€– β€– 2 + ξ‚€ βˆ’ 1 2 πœ• 𝑙 𝛾 𝑖 𝑗 + 𝛿 𝑙 ( 𝑖 ( Ξ“ 𝑗 ) βˆ’ 8 πœ• 𝑗 ) πœ‘ ) , πœ• 𝑙 𝐴 𝑖 𝑗  = ∫ πœ• Ξ©  βˆ’ 1 2 πœ• 𝛾 πœ• 𝑛 𝑖 𝑗 + 𝑛 ( 𝑖 ( Ξ“ 𝑗 ) βˆ’ 8 πœ• 𝑗 ) ξ‚„ 𝐴 πœ‘ ) 𝑖 𝑗 . ( 7 . 1 0 ) Replacing 𝐴 𝑖 𝑗 with (7.4) in the second term and using πœ• 𝑑 [ 𝛿 𝑙 ( 𝑖 ( Ξ“ 𝑗 ) βˆ’ 8 πœ• 𝑗 ) πœ‘ ) ] = 0 , we replace the last identity with 1 2 πœ• 𝑑  β€– β€– 𝐴 β€– β€– 2 + β€– β€– β€– 1 2 πœ• 𝑙 𝛾 𝑖 𝑗 βˆ’ 𝛿 𝑙 ( 𝑖 ( Ξ“ 𝑗 ) βˆ’ 8 πœ• 𝑗 ) β€– β€– β€– πœ‘ ) 2 ξ‚„ = ∫ πœ• Ξ©  βˆ’ 1 2 πœ• 𝛾 πœ• 𝑛 𝑖 𝑗 + 𝑛 ( 𝑖 ( Ξ“ 𝑗 ) βˆ’ 8 πœ• 𝑗 ) ξ‚„ 𝐴 πœ‘ ) 𝑖 𝑗 ( 7 . 1 1 )

From (7.7)–(7.11), we observe that the energy β€– β€– πœ… β€– β€– πœ– = 2 β€– β€– πœ• + 3 6 𝑙 πœ‘ β€– β€– 2 + β€– β€– Ξ“ 𝑙 βˆ’ 8 πœ• 𝑙 πœ‘ β€– β€– 2 + β€– β€– 𝐴 β€– β€– 2 + β€– β€– β€– 1 2 πœ• 𝑙 𝛾 𝑗 𝑖 βˆ’ 𝛿 𝑙 ( 𝑖 ( Ξ“ 𝑗 ) βˆ’ 8 πœ• 𝑗 ) β€– β€– β€– πœ‘ ) 2 ( 7 . 1 2 ) has its growth determined by the boundary terms πœ• 𝑑 ξ€œ πœ– = βˆ’ 6 πœ• Ξ© ξ‚€ πœ• πœ‘  ξ€œ πœ• 𝑛 πœ… βˆ’ πœ• Ξ© ξ‚€ πœ• 𝛾 πœ• 𝑛 𝑖 𝑗  𝐴 𝑖 𝑗 ξ€œ + 2 πœ• Ξ© 𝑛 ( 𝑖 ( Ξ“ 𝑗 ) βˆ’ 8 πœ• 𝑗 ) πœ‘ ) 𝐴 𝑖 𝑗 . ( 7 . 1 3 )

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 1 πœ• 𝑛 πœ‘ ( 0 ) πœ… βˆ’ 2 πœ• πœ• πœ• 𝑛 𝛾 3 ( 0 ) 𝐴 3 βˆ’ 2 πœ• πœ• 𝑛 𝛾 4 ( 0 ) 𝐴 4 βˆ’ 6 ξ‚€ 𝑛 πœ• 𝑛 𝛾 5 ( 0 ) 𝐴 5 + 4 𝑖 Ξ“ 𝑖 πœ• ( 0 ) βˆ’ 8  ξ‚„ , πœ• 𝑛 πœ‘ ( 0 ) 𝐴 5 ( 7 . 1 4 ) where 𝛾 1 𝛾 5 are coefficients of the decomposition 𝛾 𝑖 𝑗 ξ‚€ 𝑛 = 𝛾 1 ( 𝑖 π‘š 𝑗 )  ξ‚€ 𝑛 + 𝛾 2 ( 𝑖 𝑙 𝑗 )  ξ‚€ 𝑙 + 𝛾 3 ( 𝑖 π‘š 𝑗 )  ξ‚€ 𝑙 + 𝛾 4 𝑖 𝑙 𝑗 βˆ’ π‘š 𝑖 π‘š 𝑗  ξ‚€ + 𝛾 5 2 𝑛 𝑖 𝑛 𝑗 βˆ’ 𝑙 𝑖 𝑙 𝑗 βˆ’ π‘š 𝑖 π‘š 𝑗  . ( 7 . 1 5 ) 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 πœ‘ πœ• 0 π‘˜ = βˆ’ 6 π‘Ž πœ• 𝑝 πœ• 𝑝 πœ• πœ‘ + 𝐹 , ( 7 . 1 6 ) 0 1 πœ‘ = βˆ’ 6 1 π‘Ž π‘˜ + 6 πœ• 𝑙 𝑏 𝑙 , πœ• ( 7 . 1 7 ) 0 β„Ž 𝑖 𝑗 𝐴 = βˆ’ 2 π‘Ž 𝑖 𝑗 β„Ž + 2 𝑙 ( 𝑖 πœ• 𝑗 ) 𝑏 𝑙 βˆ’ 2 3 β„Ž 𝑖 𝑗 πœ• 𝑙 𝑏 𝑙 , 𝑒 ( 7 . 1 8 ) 4 πœ‘ πœ• 0 𝐴 𝑖 𝑗 = βˆ’ π‘Ž πœ• 𝑝  1 2 πœ• 𝑝 β„Ž 𝑖 𝑗 βˆ’ β„Ž 𝑝 ( 𝑖 ( Ξ“ 𝑗 ) βˆ’ 8 πœ• 𝑗 ) ξ‚„ πœ‘ ) + 𝐺 𝑖 𝑗 πœ• , ( 7 . 1 9 ) 0 Ξ“ 𝑖 4 = βˆ’ 3 π‘Ž πœ• 𝑖 π‘˜ + 𝑆 𝑖 , ( 7 . 2 0 ) where β„Ž 𝐹 = 6 π‘Ž 𝑝 π‘ž Ξ“ 𝑙 𝑝 π‘ž πœ• 𝑙 ξ‚€ πœ• πœ‘ βˆ’ 4 8 π‘Ž 𝑝 πœ‘ πœ•  ξ‚€ 𝑝 πœ‘  βˆ’ 1 4 𝑒 6 πœ‘ ξ‚€ πœ• 𝑝 πœ‘ πœ•  ξ‚€ 𝑝 𝑄  βˆ’ 𝑒 6 πœ‘ ξ‚€ 𝐷 𝑝 𝐷 𝑝 𝑄  + 1 3 𝑒 4 πœ‘ π‘Ž π‘˜ 2 + 𝑒 4 πœ‘ π‘Ž 𝐴 𝑝 π‘ž 𝐴 𝑝 π‘ž , 𝐺 𝑖 𝑗 ξ‚€ πœ• = βˆ’ π‘Ž 𝑝 β„Ž 𝑝 ( 𝑖 Ξ“  ξ‚€ 𝑗 ) βˆ’ 8 πœ• 𝑗 ) πœ‘  + 8 π‘Ž Ξ“ 𝑝 𝑖 𝑗 πœ• 𝑝 ξ‚€ πœ• πœ‘ βˆ’ 1 2 π‘Ž 𝑖 πœ‘ πœ•  ξ‚€ 𝑗 πœ‘  β„Ž + 4 π‘Ž 𝑖 𝑗 ξ‚€ πœ• 𝑝 πœ‘ πœ•  ξ‚€ 𝑝 πœ‘  βˆ’ 8 𝑒 6 πœ‘ ξ‚€ πœ• ( 𝑖 πœ‘ πœ•  ξ‚€ 𝑗 ) 𝑄  βˆ’ 𝑒 6 πœ‘ ξ‚€ 𝐷 𝑖 𝐷 𝑗 𝑄  + 8 3 𝑒 6 πœ‘ β„Ž 𝑖 𝑗 ξ‚€ πœ• 𝑝 πœ‘ πœ•  ξ‚€ 𝑝 𝑄  + 1 3 𝑒 6 πœ‘ β„Ž 𝑖 𝑗 ξ‚€ 𝐷 𝑝 𝐷 𝑝 𝑄  + 𝑒 4 πœ‘ π‘Š 𝑖 𝑗 . ( 7 . 2 1 ) Taking the scalar product of (7.16) with π‘˜ and integrating by parts on the right-hand side, we obtain ξ€œ Ξ© ξ‚€ πœ• 0 π‘˜  π‘˜ 𝑒 4 πœ‘ ξ€œ = 6 Ξ© ξ‚€ πœ• 𝑝 πœ‘ πœ•  ξ‚€ π‘ž  β„Ž ( π‘Ž π‘˜ ) 𝑝 π‘ž ξ€œ βˆ’ 6 πœ• Ξ© ξ‚€ 𝑛 𝑝 πœ• 𝑝 πœ‘  ξ€œ π‘Ž π‘˜ + Ξ©  6 ξ‚€ πœ• 𝑝 πœ‘  π‘Ž π‘˜ πœ• π‘ž β„Ž 𝑝 π‘ž ξ‚„ , + 𝐹 π‘˜ ( 7 . 2 2 ) 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: 1 2 πœ• 0  β€– β€– π‘˜ β€– β€– 2 βˆ— β€– β€– πœ• + 3 6 𝑙 πœ‘ β€– β€– 2 ξ‚„ ξ€œ = βˆ’ 6 πœ• Ξ© ξ‚€ 𝑛 𝑝 πœ• 𝑝 πœ‘  ξ€œ π‘Ž π‘˜ + Ξ© 𝐻 , ( 7 . 2 3 ) where β€– π‘˜ β€– 2 βˆ— = ∫ Ξ© π‘˜ 2 𝑒 4 πœ‘ , and ξ‚€ πœ• 𝐻 = 6 𝑝 πœ‘  [ πœ• π‘ž πœ• 𝑠 𝑏 𝑠 ξ‚€ πœ• + 6 π‘ž 𝑏 𝑠 πœ•  ξ‚€ 𝑠 πœ‘  ] β„Ž 𝑝 π‘ž ξ‚€ πœ• + 3 6 𝑝 πœ‘ πœ•  ξ‚€ π‘ž πœ‘  Γ— ξ‚€ π‘Ž 𝐴 𝑝 π‘ž βˆ’ πœ• ( 𝑝 𝑏 π‘ž ) + 1 3 β„Ž 𝑝 π‘ž πœ• 𝑠 𝑏 𝑠  + π‘˜ 2 𝑒 4 πœ‘ ξ‚€ βˆ’ 1 3 1 π‘Ž π‘˜ + 3 πœ• 𝑠 𝑏 𝑠  ξ‚€ πœ• + 6 𝑝 πœ‘  π‘Ž π‘˜ πœ• π‘ž β„Ž 𝑝 π‘ž + 𝐹 π‘˜ . ( 7 . 2 4 ) Also, from (7.17), it follows that 1 2 πœ• 0 β€– β€– πœ‘ β€– β€– 2 = ξ€œ Ξ©  βˆ’ 1 6 1 π‘Ž π‘˜ + 6 πœ• 𝑙 𝑏 𝑙 ξ‚„ πœ‘ . ( 7 . 2 5 ) Next, we notice from (7.16) and (7.20) that πœ• 0 ξ‚€ Ξ“ 𝑖 βˆ’ 8 πœ• 𝑖 πœ‘  ξ‚€ πœ• = 8 π‘Ž 𝑖 πœ‘  4 π‘˜ + 3 𝑒 6 πœ‘ ξ‚€ πœ• 𝑖 𝑄  4 π‘˜ βˆ’ 3 πœ• 𝑖 πœ• 𝑠 𝑏 𝑠 ξ‚€ πœ• βˆ’ 8 𝑖 𝑏 𝑠 πœ•  ξ‚€ 𝑠 πœ‘  + 𝑆 𝑖 . ( 7 . 2 6 ) Therefore, 1 2 πœ• 0 β€– β€– Ξ“ 𝑖 βˆ’ 8 πœ• 𝑖 πœ‘ β€– β€– 2 = ξ€œ Ξ© 𝐽 , ( 7 . 2 7 ) where  ξ‚€ πœ• 𝐽 = 8 π‘Ž 𝑝 πœ‘  4 π‘˜ + 3 𝑒 6 πœ‘ ξ‚€ πœ• 𝑝 𝑄  4 π‘˜ βˆ’ 3 πœ• 𝑝 πœ• 𝑠 𝑏 𝑠 ξ‚€ πœ• βˆ’ 8 𝑝 𝑏 𝑠 πœ•  ξ‚€ 𝑠 πœ‘  + 𝑆 𝑝 Ξ“ ξ‚„ ξ‚€ π‘ž βˆ’ 8 πœ• π‘ž πœ‘  β„Ž 𝑝 π‘ž + 1 2 ξ‚€ Ξ“ 𝑝 βˆ’ 8 πœ• 𝑝 πœ‘ Ξ“  ξ‚€ π‘ž βˆ’ 8 πœ• π‘ž πœ‘ 𝐴   2 π‘Ž 𝑝 π‘ž βˆ’ πœ• ( 𝑝 𝑏 π‘ž ) + 2 3 β„Ž 𝑝 π‘ž πœ• 𝑠 𝑏 𝑠 ξ‚„ . ( 7 . 2 8 ) Finally, taking scalar product of (7.19) with 𝐴 , integrating by parts in the right-hand side, using (7.18) to replace 𝐴 with πœ• 0 β„Ž , and regrouping, we derive our last energy identity: 1 2 πœ• 0  β€– β€– 𝐴 β€– β€– 2 βˆ— + β€– β€– β€– 1 2 πœ• 𝑙 β„Ž 𝑖 𝑗 βˆ’ β„Ž 𝑙 ( 𝑖 ξ‚€ Ξ“ 𝑗 ) βˆ’ 8 πœ• 𝑗 ) πœ‘  β€– β€– β€– 2 ξ‚„ ∫ = βˆ’ πœ• Ξ©  1 2 ξ‚€ 𝑛 𝑝 πœ• 𝑝 β„Ž 𝑖 𝑗  βˆ’ 𝑛 ( 𝑖 ξ‚€ Ξ“ 𝑗 ) βˆ’ 8 πœ• 𝑗 ) πœ‘ π‘Ž 𝐴  ξ‚„ 𝑖 𝑗 + ∫ Ξ© 𝐾 , ( 7 . 2 9 ) where β€– 𝐴 β€– 2 βˆ— = ∫ Ξ© 𝐴 𝑖 𝑗 𝐴 𝑖 𝑗 𝑒 4 πœ‘ , and  1 𝐾 = 2 πœ• 𝑝 β„Ž 𝑖 𝑗 βˆ’ β„Ž 𝑝 ( 𝑖 ξ‚€ Ξ“ 𝑗 ) βˆ’ 8 πœ• 𝑗 ) πœ‘ Γ— πœ•  ξ‚„   π‘ž ξ‚€ β„Ž 𝑠 ( π‘š πœ• 𝑛 ) 𝑏 𝑠 βˆ’ 1 3 β„Ž π‘š 𝑛 πœ• 𝑠 𝑏 𝑠  + ξ‚€ πœ• π‘ž 𝑏 𝑠 πœ•  ξ‚€ 𝑠 β„Ž π‘š 𝑛  βˆ’ ξ‚€ 𝐴 βˆ’ 2 π‘Ž π‘ž ( π‘š + β„Ž 𝑙 π‘ž πœ• ( π‘š 𝑏 𝑙 + β„Ž 𝑠 ( π‘š πœ• | | π‘ž | | 𝑏 𝑠 βˆ’ 2 3 β„Ž π‘ž ( π‘š πœ• | | 𝑠 | | 𝑏 𝑠 Ξ“  ξ‚€ 𝑛 ) βˆ’ 8 πœ• 𝑛 ) πœ‘  βˆ’ β„Ž π‘ž ( π‘š ξ‚€ ξ‚€ πœ• 8 π‘Ž 𝑛 ) πœ‘  4 π‘˜ + 3 𝑒 6 πœ‘ ξ‚€ πœ• 𝑛 ) 𝑄  4 π‘˜ βˆ’ 3 πœ• 𝑛 ) πœ• 𝑠 𝑏 𝑠 ξ‚€ πœ• βˆ’ 8 𝑛 ) 𝑏 𝑠  Γ— ξ‚€ πœ• 𝑠 πœ‘  + 𝑆 𝑛 ) β„Ž  ξ‚„ 𝑝 π‘ž β„Ž 𝑖 π‘š β„Ž 𝑗 𝑛 𝐴 + π‘Ž π‘š 𝑛 πœ• π‘ž ξ‚€ β„Ž 𝑝 π‘ž β„Ž 𝑖 π‘š β„Ž 𝑗 𝑛   + 𝐺 𝑖 𝑗 𝐴 𝑖 𝑗 . ( 7 . 3 0 ) 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

  1. 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
  2. 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
  3. 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
  4. 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
  5. 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
  6. 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
  7. 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
  8. 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
  9. 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
  10. 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
  11. 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
  12. 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
  13. 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
  14. 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
  15. 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
  16. 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
  17. 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
  18. 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
  19. 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
  20. 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
  21. 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
  22. 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
  23. 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
  24. 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
  25. 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
  26. 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
  27. 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
  28. M. Alcubierre, B. Brügmann, P. Diener et al., β€œGauge conditions for long-term numerical black hole evolutions without excision,” Phys