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

Authors: Alexander M. Alekseenko
Journal: Abstract and Applied Analysis
 (2008)

---

## 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.


---

## Body

## 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 Section2, 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 the3+1 Arnowitt-Deser-Misner (ADM) decomposition
[29, 30]:
(2.1)∂0hij=−2akij+2hl(i∂j)bl,(2.2)∂0kij=a[Rij+(kll)kij−2kilkjl]+kil∂jbl+klj∂ibl−DiDja,(2.3)Rii+(kii)2−kijkij=0,(2.4)Djkij−Dikjj=0.
Here a denotes the
lapse, the bi are the components of the shift vector b, and hij are the components of the spatial metric h. The components of the 4-dimensional metric g are given by
(2.5)g00=−a2+bibjhij,g0i=bi,gij=hij.
Here hij denotes the matrix inverse to hij. Indices on all other quantities are raised and
traces taken with respect to the spatial metric. Also, ∂0:=(∂t−bs∂s) is the convective derivative and Di is the covariant derivative operator associated with the spatial metric. The extrinsic curvature kij is defined by (2.1). We assume that global Cartesian coordinates t=x0, x1, x2, x3 are specified. Furthermore, Rij are the components of the spatial Ricci tensor
(2.6)Rij=12hpq(∂p∂jhiq+∂i∂phqj−∂p∂qhij−∂i∂jhpq)+hpqhrs(ΓipsΓqjr−ΓpqsΓijr),
where Γijk are the spatial Christoffel symbols defined by Γijk=(∂ihkj+∂jhik−∂khji)/2.The operatorRij 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 kij, in general. However, a simplified equation can be derived for the trace of the extrinsic curvature k=kii. Contracting (2.2), and using (2.1) and (2.3), we find
(2.7)∂0kii=aklikli−DlDla.
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 k=kii and the trace-free part of the extrinsic curvature Aij=kij−(1/3)hijk as new variables. Then (2.7) yields
(2.8)∂0k=13ak2+aAlm⁡Alm−DlDla.
Unless the lapse function a 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 ∂tk≥(1/3)ak2 which implies that k≥[(1/3)∫0ta(τ)dτ+1/k(0)]−1, or that the solution k 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)DlDla=aklmklm;
with this condition, (2.8) reduces
to ∂0k=0.Alternatively, we can use harmonic slicing [2, 28] (a particular case of the
Bona-Massó family of k-driving slicing conditions (∂t−blDl)a=−a2f(a)k, f(a)>0 [31, 32]), which corresponds to setting
(2.10)∂0a=−a2k.The equation onA is obtained from (2.2), (2.8), and (2.1) as
(2.11)∂0Aij=aRij+13akAij−2aAilAjl+29ak2hij−13aAlmAlmhij+Ail∂jbl+Ajl∂ibl−DiDja+13hijDlDla.
To proceed with the derivation, we need a splitting for the spatial metric h compatible with the splitting of kij into k and Aij. In the BSSN formulation, the desired splitting is achieved by introducing the conformal factor φ=(1/12)ln⁡(det⁡(hij)) and the conformal metric hij=e−4φhij, hij=e4φhij. Using the Leibnitz formula for differentiating the
determinant of a matrix
(2.12)∂det⁡(hij)=det⁡(hij)hlm∂hml,
one finds that the derivative
(2.13)∂hij=e−4φ[∂hij−13hijhlm∂hlm]
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
(2.14)∂0φ=−16ak+16∂lbl.
Now using (2.13) and (2.1), we obtain the third equation
(2.15)∂0hij=−2aAij+2hl(i∂j)bl−23hij∂lbl,
where Aij=e−4φAij and bi=e−4φbi are the conformal analogs of the variables A and b. Beginning with the last equation, indices are lowered and raised with the conformal metric hij and its inverse hij=e4φhij. In this case, bs=bs, and it is easy to redefine ∂0=∂t−bs∂s.The remaining two equations can be obtained from (2.11) which can be rewritten in terms of A as
(2.16)∂0Aij=ae−4φRij+a(kAij−2AilAjl+29k2hij−13Alm⁡Alm⁡hij)+Ail∂jbl+Ajl∂ibl−23Aij∂lbl−e−4φDiDja+e−4φ13hijDlDla.
The Ricci tensor in terms of the conformal metric becomes [2, 28]
(2.17)Rij=12hpq(∂p∂jhiq+∂i∂phqj−∂p∂qhij−∂i∂jhpq)−2DiDjφ−2hijhpqDpDqφ+hpqhrs(ΓipsΓqjr−ΓpqsΓijr)+4∂iφ∂jφ−4hijhpq∂pφ∂qφ.
Here Γijk=(∂ihkj+∂jhik−∂khji)/2 and Divj=∂ivj−hpqΓijpvq are the covariant derivatives associated with the conformal metric. The first line in (2.17) can be rewritten to obtain
(2.18)Rij=−12hpq∂p∂qhij+∂(ihpqΓ|pq|j)+Γpq(i∂j)hpq−2DiDjφ−2hijhpqDpDqφ+hpqhrs(ΓipsΓqjr−ΓpqsΓijr)+4∂iφ∂jφ−4hijhpq∂pφ∂qφ.
This suggests that we introduce a new variable
(2.19)Γj=hpqΓpqj=hpq∂phqj.
Substituting (2.18) in (2.16), we derive the fourth evolution equation
(2.20)∂0Aij=−12ae−4φhpq∂p∂qhij+ae−4φ∂(iΓj)−2ae−4φDiDjφ−2ae−4φhijhpqDpDqφ−e−4φDiDja+13e−4φhijDlDla+Wij,
where
(2.21)Wij=ae−4φΓpq(i∂j)hpq+ae−4φhpqhrs(ΓipsΓqjr−ΓpqsΓijr)+4ae−4φ∂iφ∂jφ−4ae−4φhijhpq∂pφ∂qφ+a(kAij−2AilAjl+29k2hij−13Alm⁡Alm⁡hij)+Ail∂jbl+Ajl∂ibl−23Aij∂lbl.The evolution equation forΓj is obtained by applying the operator ∂0 to (2.19) and using the momentum constraint
(2.22)∂0Γi=−2a∂pApi−2hpq(∂pa)Aqi+2aApq(∂phqi)+Γl∂ibl+13∂i∂lbl+hli∂s∂sbl.
Next we notice that hpqDpAiq=∂pApi−ΓsAis+6(∂sφ)Ais, and thus (2.4) takes the form
(2.23)∂pApi−23∂ik−ΓsAis+6(∂sφ)Ais=0.
Solving this equation for ∂pApi and substituting the result in (2.22), we derive the fifth equation of the BSSN system
(2.24)∂0Γi=−43a∂ik+Si,
where
(2.25)Si=−2aΓsAis+12a(∂sφ)Ais−2hpq(∂pa)Aqj+2aApq(∂phqj)+Γl∂jbl+13∂i∂lbl+hlj∂s∂sbl.
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 bi or a shift that follows from the gamma-freezing condition ∂tΓi=0 (cf. [28]).


## 3. Linearization around Minkowski Space

Minkowski spacetime in Cartesian coordinates is represented by the trivial solution to the ADM system:hij=δij, kij=0, a=1, bi=0. Consider perturbations of the ADM variables around
the Minkowski spacetime: hij=δij+γij, kij=κij, a=1+α, and bi=βi, with the γij, κij, α, and βi 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
(3.1)det⁡(hij)=1+γll,φ=112γll,e−4φ=1−13γll,e4φ=1+13γll,hij=δij+γij−13δijγll=:δij+γij,k=κ=:κll,Aij=Aij=κij−13δijκ,Γi=:Γi=∂lγil.
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 φ, γij, κ, Aij, and Γi, we derive the linearization of the BSSN system
(3.2)∂tφ=−16κ+16∂lβl,(3.3)∂tα=−κ,(3.4)∂tκ=−∂l∂lα,(3.5)∂tγij=−2Aij+2∂(iβj)−23δij∂lβl,(3.6)∂tAij=−12∂l∂lγij+∂(iΓj)−2∂i∂jφ−2δij∂l∂lφ−∂i∂jα+13δij∂l∂lα,(3.7)∂tΓi=−43∂iκ+13∂i∂pβp+∂p∂pβi.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(3.8)∂l∂jγlj−8∂l∂lφ=0,∂lΓl−8∂l∂lφ=0,(3.9)∂lAil−23∂iκ=0.
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)Γi=∂lγil.The linearized problem then consists of determiningφ, α, κ, γ, A, Γ 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), A(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), A(0), Γ(0) satisfy the constraint equations (3.8)–(3.10). The linearized Hamiltonian and momentum constraints in the ADM system are ∂l∂iγli−∂l∂lγii=0 and ∂lκil−∂iκll=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
(4.1)∂t(∂l∂jγlj−8∂l∂lφ)=∂i(∂lAil−23∂iκ),∂t(∂lΓl−8∂l∂lφ)=0.
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
(4.2)Mi=∂lAil−23∂iκ.
Equation (3.9) is satisfied if
and only if Mi=0.Now let us derive an equation for the evolution ofM. By differentiating (3.6) with respect to time and using (3.5) to replace ∂tγij in the result, we obtain
(4.3)∂t2Aij=∂l∂lAij.
Terms in φ, α, κ, β, Γ cancel in view of (3.2)–(3.4), (3.7).Similarly, by differentiating (3.4) and substituting (3.3) for ∂tα, one derives
(4.4)∂t2κ=∂l∂lκ.
Differentiating (4.2) twice with respect to time and using (4.3) and (4.4), we obtain
(4.5)∂t2Mi=∂l∂lMi.
The initial value M(0) can be determined from A(0) and κ(0) using (4.2) and
must be zero for the physical initial data. The initial values ∂tM(0) can be calculated by differentiating (4.2) in time and substituting (3.4) and (3.6) for ∂tκ and ∂tAij. Thus we have
(4.6)∂tMi=−12∂l∂l∂mγim+12∂i∂lΓl+12∂l∂lΓi−4∂i∂l∂lφ.
It can be verified by substitution that if γ(0), Γ(0), and φ(0) satisfy (3.8) and (3.10) then ∂tMi(0)=0.We must now select the boundary conditions onM that imply the trivial evolution of (4.5). However, the boundary conditions on M are expected not to be given freely but rather be determined by the boundary conditions (and the data) on the variables A and κ. Similarly, the initial data M(0) and ∂tM(0) is determined by A(0), κ(0), γ(0), Γ(0), and φ(0). Our goal then is to construct the boundary conditions on the main variables A and κ in such a way as to guarantee homogeneous data for M. 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(vi,ui)=∫Ωviuidx and (ρij,σij)=∫Ωρijσijdx for the spaces of vector fields and matrix fields on Ω, respectively. The L2 norms naturally associated with these scalar products are ∥u∥2=(ui,ui) and ∥ρ∥2=(ρij,ρij). We introduce the energy function for the constraint quantity
(4.7)ϵ=∥∂tM∥2+∥∂lM∥2.
If we prove that the energy ϵ remains zero at all times, then, in view of the trivial initial data, we will have ∥M∥≡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)∂tϵ=∫∂Ω(∂∂nMi)(∂tMi)dσ.
The energy ϵ is not increasing if
(4.9)(∂∂nMi)(∂tMi)≤0on∂Ω.
The desired boundary conditions on A 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 ni be the unit vector perpendicular to the face. Let mi and li complement ni to form an orthonormal triple. For example, mi can be chosen to be any unit vector parallel to the boundary and li be the cross product of ni and mi, so that li=(n×m)i=εijk⁢njmk. At a face of a flat boundary, the divergence of a
vector field can be expressed in terms of the directional derivatives along vectors n, m, and l as
(4.10)∂ivi=ni∂∂nvi+mi∂∂mvi+li∂∂lvi.
Similarly, the gradient of a scalar field ψ becomes
(4.11)∂iψ=ni∂∂nψ+mi∂∂mψ+li∂∂lψ.
Next we note that at any point of the boundary, a symmetric trace free matrix is spanned by
(4.12)n(imj),n(ilj),l(imj),lilj−mimj,2ninj−lilj−mimj.
Introducing the scalar functions
(4.13)A1=2Aijn(imj),A2=2Aijn(ilj),A3=2Aijl(imj),A4=12Aij(lilj−mimj),A5=16Aij(2ninj−lilj−mimj),
we rewrite A as follows:
(4.14)Aij=A1(n(imj))+A2(n(ilj))+A3(l(imj))+A4(lilj−mimj)+A5(2ninj−lilj−mimj)
Substituting (4.14) into (4.2) and using (4.10) and (4.11) to replace partial derivatives with directional derivatives, we obtain
(4.15)Mi=[12∂∂mA1+12∂∂lA2+2∂∂nA5−23∂∂nκ]ni+[12∂∂nA1+12∂∂lA3−∂∂mA4−∂∂mA5−23∂∂mκ]mi+[12∂∂nA2+12∂∂mA3+∂∂lA4−∂∂lA5−23∂∂lκ]li.
The last equation implies that
(4.16)(∂∂nMi)(∂tMi)=∂∂n[12∂∂mA1+12∂∂lA2+2∂∂nA5−23∂∂nκ]×∂t[12∂∂mA1+12∂∂lA2+2∂∂nA5−23∂∂nκ]+∂∂n[12∂∂nA1+12∂∂lA3−∂∂mA4−∂∂mA5−23∂∂mκ]×∂t[12∂∂nA1+12∂∂lA3−∂∂mA4−∂∂mA5−23∂∂mκ]+∂∂n[12∂∂nA2+12∂∂mA3+∂∂lA4−∂∂lA5−23∂∂lκ]×∂t[12∂∂nA2+12∂∂mA3+∂∂lA4−∂∂lA5−23∂∂lκ].
Either of the following two sets of boundary conditions imply that ((∂/∂n)Mi)(∂tMi)=0 on ∂Ω:
(4.17)A1=0,A2=0,∂∂nA3=0,∂∂nA4=0,∂∂nA5=0,∂∂nκ=0;(4.18)∂∂nA1=0,∂∂nA2=0,A3=0,A4=0,A5=0,κ=0.
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 A3, A4, A5, and κ drop out as a direct consequence of (4.17). Also (4.3) can be solved for the second normal derivatives of A1 and A2 in terms of the temporal and tangential derivatives using the following representation of the Laplace operator:
(4.19)∂j∂j=∂∂n∂∂n+∂∂m∂∂m+∂∂l∂∂l.
It can be easily noticed that (4.17) implies that all tangential and temporal derivatives of A1 and A2 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, Mi|∂Ω=0 is equivalent to the set of the following differential boundary conditions that can be used in a numerical method [13, 16]:
(4.20)12∂∂mA1+12∂∂lA2+2∂∂nA5−23∂∂nκ=0,12∂∂nA1+12∂∂lA3−∂∂mA4−∂∂mA5−23∂∂mκ=0,12∂∂nA2+12∂∂mA3+∂∂lA4−∂∂lA5−23∂∂lκ=0.
In this case, we can prescribe the Dirichlet data for A3 and A4 and the Neumann data for κ. Then (4.20) contains normal derivatives of A1, A2, and A5 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 A 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(4.21)∂∂nMini=0,Mili=0,Mimi=0.
By applying ∂/∂m to the second equation in (4.20), ∂/∂l 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/∂n2) derivatives) an evolution equation defined on the boundary
(4.22)2∂2∂t2A5−(∂2∂l2+∂2∂m2)A5−23∂2∂t2κ+43(∂2∂l2+∂2∂m2)κ=(∂2∂l2−∂2∂m2)A4+∂∂l∂∂mA3.
One can use (4.22) to construct
nonhomogeneous mixed Dirichlet and Neumann data that imply (4.21). Specifically, we may prescribe Dirichlet data for A3, A4, and κ that is compatible at corners but, otherwise, arbitrary. Then we can use (4.22) to evolve A5 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 A5 are determined, we can use the second and third equations of (4.20) to determine nonhomogeneous Neumann data for A1 and A2. Similarly, one can give the Dirichlet data for A3, A4, and A5 and use (4.22) to determine the Dirichlet data for κ and the Neumann data for A1 and A2. The corresponding nonhomogeneous algebraic conditions on Aij then read
(4.23)2∂∂nAijn(imj)=∂∂nA1,2∂∂nAijn(ilj)=∂∂nA2,2Aijl(imj)=A3,12Aij(lilj−mimj)=A4,16Aij(2ninj−lilj−mimj)=A5.
If desired, energy estimates may be obtained on the solutions A 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 ofA 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 A(0) and ∂tA(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(μ,ν)=∫Ωμνdxdydz and (uijk,vijk)=∫Ωuijkvijkdxdydz for the spaces of scalar fields and triple-indexed fields on Ω. The L2 norms naturally associated with these scalar products are ∥μ∥2=(μ,μ) and ∥u∥2=(uijk,uijk). Consider the system (4.3) describing the propagation of A. We define the energy of the system to be
(5.1)ϵ1=12(∥∂tA∥2+∥∂lA∥2).
As in Section 4, by differentiating ϵ1 in time, integrating terms with spatial derivatives by parts, and using (4.3), we obtain
(5.2)∂tϵ1=∫∂Ω(∂∂nAij)(∂tAij)dσ.
Since the right-hand side of (5.2) is zero if either (4.17) or (4.18) holds, we conclude that ϵ1 and, therefore, A 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), A(0), κ(0), φ(0), α(0), and Γ(0) be given, and let ∂tA(0) and ∂tκ(0) be determined by (3.6) and (3.4), respectively. If ∂tA(0) and ∂tκ(0)∈L2(Ω) and also ∂lAij(0) and ∂lκ(0)∈L2(Ω), then there exists unique solutions A and κ to (4.3) and (4.4). These solutions obey the following energy estimates:
(5.3)sup⁡0≤t≤T[∥∂tAij∥L2(Ω)+∥∂lAij∥L2(Ω)]≤∥∂tAij(0)∥L2(Ω)+∥∂lAij(0)∥L2(Ω),sup⁡0≤t≤T[∥∂tκ∥L2(Ω)+∥∂lκ∥L2(Ω)]≤∥∂tκ(0)∥L2(Ω)+∥∂lκ(0)∥L2(Ω).


## 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 A and κ be taken in either form (4.17) or (4.18). According to Theorem 5.1, the matrix A 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 Nijpq is defined below (6.4)
(6.1)Nijpqγpq+(δipδjq−Nijpq)∂∂nγpq=Nijpqγpq(0)+(δipδjq−Nijpq)∂∂nγpq(0)+∫0t(Nijpq+(δipδjq−Nijpq)∂∂n)(2∂(pβq)−23δpq∂sβs)(6.2)μφ+(1−μ)∂∂nφ=μφ(0)+(1−μ)∂∂nφ(0)+∫0t16(μ−(1−μ)∂∂n)∂lβl,(6.3)μα+(1−μ)∂∂nα=μα(0)+(1−μ)∂∂nα(0),(6.4)niΓi=niΓi(0)+∫0t[−43∂∂nκ+13∂∂n∂pβp+∂p∂pniβi],τiΓi=τiΓi(0)+∫0t[−43∂∂τκ+13∂∂τ∂pβp+∂p∂pτiβi].
The projection operator Nijpq corresponding to (4.17) is given by
(6.5)Nijpq=2n(pmq)n(imj)+2n(plq)n(ilj)
and the one corresponding to (4.18) is
(6.6)Nijpq=2l(pmq)l(imj)+12(lplq−mpmq)(lilj−mimj)+16(2npnq−lplq−mpmq)(2ninj−lilj−mimj).
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), τi can be either vector li or mi. 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
(6.7)τi∂∂nΓi=τi∂∂nΓi(0)+∫0t[−43∂∂τ∂∂nκ+∂∂n(13∂∂τ∂pβp+∂p∂pτiβi)].Theorem 6.1.
LetA 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 A(0), κ(0), ∂tA(0), and ∂tκ(0) and the initial data γ(0), Γ(0), α(0), and φ(0) be such that
(6.8)∂tAij(0)=−12∂l∂lγij(0)+∂(iΓj)(0)−2∂i∂jφ(0)−2δij∂l∂lφ(0)−∂i∂jα(0)+13δij∂l∂lα(0),∂tκ(0)=−∂l∂lα(0).
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
(6.9)φ=φ(0)+∫0t[−16κ+16∂lβl],α=α(0)−∫0tκ,γij=γij(0)+∫0t[−2Aij+2∂(iβj)−23δij∂lβl],Γi=Γi(0)+∫0t[−43∂iκ+13∂i∂pβp+∂p∂pβi].
Moreover, if A 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
(6.10)∂tAij=∂tAij(0)+∫0t∂l∂lAij,∂tκ=∂tκ(0)+∫0t∂l∂lκ,
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
(6.11)∂l∂jγlj−8∂l∂lφ=∂l∂jγlj(0)−8∂l∂lφ(0)−∫0t2∂j(∂lAjl−23∂jk),∂lΓl−8∂l∂lφ=∂lΓl(0)−8∂l∂lφ(0).
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 ∂tΓi=0 yielding an elliptic equation for β that needs to be solved at each time step:
(6.12)13∂i∂pβp+∂p∂pβi−43∂iκ=0.
The boundary conditions on β can be taken, for example, to be
(6.13)βimi=0,βili=0,∂∂nβini=0,
or
(6.14)∂∂nβimi=0,∂∂nβili=0,βini=0.
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 βi and integrating by parts, we can show that ∥βi∥H1(Ω) is bounded by ∥∂iκ∥L2(Ω) 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
(6.15)∂t2β=13∂i∂pβp+∂p∂pβi−43∂iκ,
which corresponds to a dynamic gamma-freezing condition ∂tΓi=∂t2βi. In this case, in addition to (6.4) and (6.13), one can
consider either of two sets of radiative boundary conditions
(6.16)βimi=0,βili=0,(∂tβi+∂∂nβi)ni=0,(∂tβi+∂∂nβi)mi=0,(∂tβi+∂∂nβi)li=0,βini=0.
By reasoning similar to the above, an energy estimate for ∥β∥H1(Ω) in terms of ∥∂iκ∥L2(Ω) can be obtained using techniques for hyperbolic equations (cf. [37–39, 42]).After the boundary conditions forβ, A, and κ are chosen, one can set Γi−Γ(0)i=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
(7.1)a=e6φQ,
which results in significant simplifications in our equations. The result, however, is expected to carry over to the harmonic slicing and the k-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 Q≡1 reduces to α=6φ. Assuming zero-shift perturbation βi=0 for further simplicity, (3.2), (3.4)–(3.7) are restated as
(7.2)∂tφ=−16κ,(7.3)∂tκ=−6∂l∂lφ,(7.4)∂tγij=−2Aij,(7.5)∂tAij=−12∂l∂lγij+∂(iΓj)−8∂i∂jφ,(7.6)∂tΓi=−43∂iκ.
Taking the scalar product of (7.3) with κ, integrating the result by parts, and using (7.2) to replace κ with ∂tφ, we obtain
(7.7)12∂t[∥κ∥2+36∥∂lφ∥2]=−6∫∂Ω(∂∂nφ)κ.
Also, we conclude from (7.2) that 36∥∂tφ∥2=∥κ∥2.From (7.2) and (7.6), we observe that ∂t(Γi−8∂iφ)=0, which implies that
(7.8)∂t∥Γi−8∂iφ∥2=0,(∂t(Γi−8∂iφ),∂lγli)=0.Next we rewrite the right-hand side of (7.5) in the divergence form:
(7.9)∂tAij=∂l[−12∂lγij+δl(i(Γj)−8∂j)φ)].
Then we take the scalar product with Aij,  and integrate by parts to obtain
(7.10)12∂t∥A∥2+(−12∂lγij+δl(i(Γj)−8∂j)φ),∂lAij)=∫∂Ω[−12∂∂nγij+n(i(Γj)−8∂j)φ)]Aij.
Replacing Aij with (7.4) in the second term and using ∂t[δl(i(Γj)−8∂j)φ)]=0, we replace the last identity with
(7.11)12∂t[∥A∥2+∥12∂lγij−δl(i(Γj)−8∂j)φ)∥2]=∫∂Ω[−12∂∂nγij+n(i(Γj)−8∂j)φ)]AijFrom (7.7)–(7.11), we observe that the energy
(7.12)ϵ=∥κ∥2+36∥∂lφ∥2+∥Γl−8∂lφ∥2+∥A∥2+∥12∂lγji−δl(i(Γj)−8∂j)φ)∥2
has its growth determined by the boundary terms
(7.13)∂tϵ=−6∫∂Ω(∂∂nφ)κ−∫∂Ω(∂∂nγij)Aij+2∫∂Ωn(i(Γj)−8∂j)φ)Aij.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
(7.14)∂tϵ=∫∂Ω[−6∂∂nφ(0)κ−12∂∂nγ3(0)A3−2∂∂nγ4(0)A4−6∂∂nγ5(0)A5+4(niΓi(0)−8∂∂nφ(0))A5],
where γ1 – γ5 are coefficients of the decomposition
(7.15)γij=γ1(n(imj))+γ2(n(ilj))+γ3(l(imj))+γ4(lilj−mimj)+γ5(2ninj−lilj−mimj).
Energy (7.12) is conserved if the
initial data is chosen as to satisfy (∂/∂n)φ(0)=0, γ1(0)=γ2(0)=0, and (∂/∂n)γ3(0)=(∂/∂n)γ4(0)=(∂/∂n)γ5(0)=0 at the boundary. Notice that Γi(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, (∂/∂n)γ1(0)=(∂/∂n)γ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
(7.16)e4φ∂0k=−6a∂p∂pφ+F,(7.17)∂0φ=−16ak+16∂lbl,(7.18)∂0hij=−2aAij+2hl(i∂j)bl−23hij∂lbl,(7.19)e4φ∂0Aij=−a∂p[12∂phij−hp(i(Γj)−8∂j)φ)]+Gij,(7.20)∂0Γi=−43a∂ik+Si,
where
(7.21)F=6ahpqΓpql∂lφ−48a(∂pφ)(∂pφ)−14e6φ(∂pφ)(∂pQ)−e6φ(DpDpQ)+13e4φak2+e4φaApqApq,Gij=−a(∂php(i)(Γj)−8∂j)φ)+8aΓijp∂pφ−12a(∂iφ)(∂jφ)+4ahij(∂pφ)(∂pφ)−8e6φ(∂(iφ)(∂j)Q)−e6φ(DiDjQ)+83e6φhij(∂pφ)(∂pQ)+13e6φhij(DpDpQ)+e4φWij.
Taking the scalar product of (7.16) with k and integrating by parts on the right-hand side, we obtain
(7.22)∫Ω(∂0k)ke4φ=6∫Ω(∂pφ)(∂q(ak))hpq−6∫∂Ω(np∂pφ)ak+∫Ω[6(∂pφ)ak∂qhpq+Fk],
where ni is the outward normal vector to the boundary in the sense of the conformal metric h. Substituting (7.17) for ak in the second term and regrouping, we obtain our first energy identity:
(7.23)12∂0[∥k∥∗2+36∥∂lφ∥2]=−6∫∂Ω(np∂pφ)ak+∫ΩH,
where ∥k∥∗2=∫Ωk2e4φ, and
(7.24)H=6(∂pφ)[∂q∂sbs+6(∂qbs)(∂sφ)]hpq+36(∂pφ)(∂qφ)×(aApq−∂(pbq)+13hpq∂sbs)+k2e4φ(−13ak+13∂sbs)+6(∂pφ)ak∂qhpq+Fk.
Also, from (7.17), it follows
that
(7.25)12∂0∥φ∥2=∫Ω[−16ak+16∂lbl]φ.
Next, we notice from (7.16) and (7.20) that(7.26)∂0(Γi−8∂iφ)=8a(∂iφ)k+43e6φ(∂iQ)k−43∂i∂sbs−8(∂ibs)(∂sφ)+Si.
Therefore,
(7.27)12∂0∥Γi−8∂iφ∥2=∫ΩJ,
where
(7.28)J=[8a(∂pφ)k+43e6φ(∂pQ)k−43∂p∂sbs−8(∂pbs)(∂sφ)+Sp](Γq−8∂qφ)hpq+12(Γp−8∂pφ)(Γq−8∂qφ)[2aApq−∂(pbq)+23hpq∂sbs].
Finally, taking scalar product of (7.19) with A, integrating by parts in the right-hand side, using
(7.18) to replace A with ∂0h, and regrouping, we derive our last energy identity:
(7.29)12∂0[∥A∥∗2+∥12∂lhij−hl(i(Γj)−8∂j)φ)∥2]=−∫∂Ω[12(np∂phij)−n(i(Γj)−8∂j)φ)]aAij+∫ΩK,
where ∥A∥∗2=∫ΩAijAije4φ, and
(7.30)K=[12∂phij−hp(i(Γj)−8∂j)φ)]×{[∂q(hs(m∂n)bs−13hmn∂sbs)+(∂qbs)(∂shmn)−(−2aAq(m+hlq∂(mbl+hs(m∂|q|bs−23hq(m∂|s|bs)(Γn)−8∂n)φ)−hq(m(8a(∂n)φ)k+43e6φ(∂n)Q)k−43∂n)∂sbs−8(∂n)bs)×(∂sφ)+Sn))]hpqhimhjn+aAmn∂q(hpqhimhjn)}+GijAij.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 φ, k, h, A, Γ, ∂φ, and ∂h, but not derivatives of k, A, 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 bs are prescribed. Otherwise, the terms ∂i∂jbs 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 Section7 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.



---

*Source: 742040-2008-05-26.xml*