text:
# 1 Introduction

## 1 Introduction

The Stefan problem concerns solid-liquid or liquid-vapor phase transitions when moving unknown beforehand surface $`S`$ of phase transition is formed (see, e.g., ). In fact, a formulation of the Stefan problem was given for the first time by G. Lame and B.P. Clapeiron in 1831 for a particular case of equal temperature of liquid and crystalline phases . In 1889 J. Stefan published four papers devoted to the subject (in particular, to the description of soil freezing) in which the problem was formulated in a general form . According to it, for the interphase the following condition

$$K_{sol}\frac{T(x_S+0,t)}{x}K_{liq}\frac{T(x_S0,t)}{x}=L\rho _{sol}V_S,$$
(1)
defining Stefan’s problem, has been suggested. Here $`V_S=d\xi _S/dt`$ is the velocity of the boundary surface $`S`$, $`K_{sol}`$ and $`K_{liq}`$ are thermal conductivities of material for solid and liquid phases, $`L`$ and $`\rho _{sol}`$ are the melting heat and density, correspondingly. Condition (1) has a clear physical meaning. Indeed, according to the Fourier law, heat flow $`j`$ is proportional to the temperature gradient,

$$𝐣=KgradT.$$

Therefore, the left-hand side of (1) is the heat absorbed in the unit of area per the unit of time. The expression in the right-hand side is the heat connected with freezing or melting of material crossed per the unit of time by the the unit of area.

A complete mathematical formulation of the Stefan problem includes, besides (1), a condition of continuity on the surface $`S`$ separating solid and liquid phases

$$T|_S=T^{},$$
(2)
where $`T^{}`$ denotes temperature of the phase transition, and the energy conservation law

$$\rho C\frac{T}{t}=div𝐣+q(𝐱,t).$$

Here $`q(𝐱,t)`$ represents the power of external heat sources, $`C`$ is the specific heat. In the original Stefan papers $`q(𝐱,t)0,`$ so that the whole heat transfer has been considered to be a consequence of the temperature gradient inside the medium.

If one also specifies initial and boundary conditions, the Stefan problem can be solved more often approximately, but sometimes exactly. Particular examples of suitable boundary conditions are considered below.

Relations (1) and (2) are usually used in numerical algorithms explicitly. Another approach was suggested by A.N. Tikhonov and A.A. Samarskii in 1953 . According to it, conditions (1) and (2) are included themselves into the energy conservation equation to obtain generalized formulation of the Stefan problem in the form

$$(\rho C+L\delta (TT^{}))\left(\frac{T}{t}+𝐯gradT\right)=div(KgradT)+q(𝐱,t),$$
(3)
where the term $`L\delta (TT^{})T/t`$ describes the additional heat input expended on the phase transformation, $`𝐯gradT`$ takes into account possible temperature change due to convection (hereafter we ignore it for simplicity). The main idea of this approach is quite clear, too. Namely, it is suggested to treat the heat of fusion $`L`$ as an additional component of the thermal capacity $`\rho C`$ which gives contribution only at the point of phase transition.

Lately Samarskii and his followers have turned this idea into effective numerical algorithms (see, e.g. ). But even in those papers equation (3) is considered only as a corollary of the condition (1). For example, it was derived in  by substituting expression $`L\delta (TT^{})T/t`$ instead of the term $`L\delta (x\xi _S(t))V_S`$, which is assumed to be included in the heat equation to account for the heat absorption on 2-dimensional interface $`S`$.

The purpose of this paper is to show that the condition (3) supplies us with a more powerful description of phase transitions, that may be used even in the case when (1) and (2) are not applicable.

## 2 Heuristic arguments

As it was mention above, the possibility of solving the classical Stefan problem by making use of condition (3) has been demonstrated by Samarskii and his co-authors. Therefore, we only consider an example when (3) is applicable and (1), (2) are not. To this end let us study the following problem:

$$(\rho C+L\delta (TT^{}))\frac{T}{t}div(kgradT)+q(t),$$
(4)
$$T(𝐱,0)=T_0<T^{},$$

where all parameters of (4) are suggested to be independent of $`𝐱`$. Due to the spatial uniformity, it is evident that the condition

$$gradT=0$$

holds on the solutions of (4). In this case Eq. (4) is reduced to an ordinary differential one

$$(\rho C+L\delta (TT^{}))\frac{dT}{dt}=q(t)$$
(5)
with the initial condition

$$T(0)=T_0.$$

Integrating both sides of (5) over $`t`$ just near the phase transition temperature $`T^{}`$, one obtains

$$_{T^{}0}^{T^{}+0}(\rho C+L\delta (TT^{}))𝑑T=_t^{t+\delta t}q(t)𝑑t,$$
(6)
where $`\delta t`$ is a time necessary for the phase transition. It is evident from (6) that

$$\delta t\frac{L}{Q},$$
(7)
where $`Q`$ is the maximum value of $`q(t)`$ in the interval $`(t,t+\delta t)`$. The inequality (7) means that the phase transition at a fixed spatial point lasts a finite, distinct from zero, time.

This simple example shows something completely different from the Stefan description of the phase transition. Let us examine it carefully.

1. First of all, instead of gradual warming (or cooling) up the pattern due to the influence of one of its boundaries, here we have an uniformly heated layer. Therefore, creation of 2-D surface $`S(y,z)`$ separating the solid and liquid phases in $`x`$ is evidently impossible due to a total equivalence of all spatial points $`x`$.
2. One can also expect that finiteness of the phase transition time, $`\delta t`$, forces all points within a spatial layer of nonzero thickness to be at the same temperature $`T^{}`$. This is expected even in the case when the power deposition $`q(x,t`$), unlike in the example considered, is spatially irregular111 Indeed, let material has just reached the temperature $`T^{}`$ at some point $`𝐱`$ and now begins receiving its portion of the heat necessary for melting. Then another adjacent point $`𝐱`$ $`+`$ $`\mathrm{\Delta }𝐱`$, which attained the melting temperature merely a little earlier, can be still in the state of heat receiving and, therefore, must have the same temperature $`𝐓^{}`$ (see Fig.1)..

If we consider $`\delta `$-function in (3) as a limit of a bounded function $`D(TT^{})`$ localized in the vicinity of $`T=T^{}`$, then the possibility to obtain the solution shown in Fig. 1(a), or analytically,

$$\frac{T}{t}0,gradT0,$$

follows from the indefiniteness

$$D(TT^{})\frac{T}{t}\mathrm{}\mathrm{\hspace{0.33em}0},$$

springing up in the left-hand side of (3). Clearly, that this indefiniteness can take a finite value and compensate in a space region with nonzero thickness the spatially distributed source $`q(𝐱,t)`$ which contributes to the right-hand side of (3). This, of course, is no more true if the external sources are absent and heat enters the pattern only through its boundary.

In a general case, one can expect existence of two jumps for spatial derivatives of temperature on the boundaries $`S`$ of the volume $`V_T^{}`$ with $`T=T^{}`$, instead of one for the classical Stefan problem, but the condition (1) is hardly met for any of them (see Fig. 1(a), where intersection of the boundary $`S`$ by (x,y)-plane in points 1, 2, 3 and 4 is seen). Indeed, to prove the existence of two jumps — one from the side of the solid and another from the side of the melted phase — it is sufficient only to show that the spatial derivative on the surface $`S`$, taken externally, is not equal to zero. The co-ordinates of the boundary $`\stackrel{}{\xi }(t)`$ can be found as a solution of an equation

$$T(𝐱,t)T^{}=0,$$

where $`T(𝐱,t)`$ is the solution of the heat equation (3) outside the volume $`V_T^{}`$. Taking the total temporal derivative, one obtains

$$\frac{T}{t}+gradT\frac{d\stackrel{}{\xi }}{dt}=0.$$

Thus $`gradT=0`$ automatically implies $`T/t=0.`$ It is evident that such conditions are impossible if the external sources are not adjusted specially to stabilize the temperature in the infinitesimal layers adjacent to the volume $`V_T^{}`$ just before and just after the phase transition.

## 3 Beam induced phase transitions

To verify the conclusions which we have just come to, let us study numerically the dynamics of phase transition induced by a short powerful ion beam in solids. At present this technology is really used for modification of surface layers to create new materials with unique physical and chemical properties (see, e.g. ). The process is underlain by the equation for heat transfer which we discussed in the previous sections:

$$\rho (T)c(T)\frac{T}{t}=\frac{}{x}\left(k(T)\frac{T}{x}\right)+q.$$
(8)
The initial and boundary conditions could be taken in the form:

$$T(x,0)=T_0,\frac{T(0,t)}{x}=\frac{T(l_0,t)}{x}=0.$$

Let us consider, for definiteness, an iron pattern, which thermal properties are described in popular reference books, and choose dimensionless (DL for brevity) variables

$$T:=T/T_0,x:=x/l_0,t:=t/\tau $$

as follows:

$`T_0=293K,`$ $`l_0=10^5m(\mathrm{the}\mathrm{pattern}\mathrm{thickness}),`$

$`\tau =\mathrm{3\hspace{0.33em}10}^7s(\mathrm{duration}\mathrm{of}\mathrm{ion}\mathrm{beam}\mathrm{pulse}\mathrm{from}\mathrm{an}\mathrm{accelerator}).`$

For DL power deposition $`q`$ we take a simple model, shown in Fig. 2 and 3,

with analytical representation

$$q(x,t)=Qq_1(x)q_2(t),$$

where

$$q_i(z)=\frac{1}{1+\mathrm{exp}\mu _i(zz_i)}$$

and $`Q`$ describe the total DL energy brought into the pattern (here $`Q=59.44,x_1=0.07,t_1=1,\mu _i=100`$). For simplicity, we neglect in (8) a small difference between physical parameters for the solid and liquid phases.

Now, using of the general idea due to Tikhonov and Samarskii , we assume an expression:

$$\rho (T)c(T)=1+\lambda \delta (TT^{},\mathrm{\Delta })$$

for DL specific heat, where $`\lambda `$ denotes the DL heat of fusion and $`\delta (TT^{},\mathrm{\Delta })`$ is an approximate $`\delta `$-function, smoothed with the help of the Gaussian distribution of width $`\mathrm{\Delta }`$ (see Fig. 4)222There were other methods of smoothing in original papers by Samarskii et al. They used regularization on the space grid..

Now Eq. (8) can be solved numerically on the space-time grid $`x`$ and $`t`$ with steps $`h_x`$ and $`h_t`$, within intervals $`x(0,1),t(0,t_{max})`$:

$$x_j=h_xj,j=0,\mathrm{},n_x,h_x=1/n_x,$$

$$t_k=h_tk,k=0,\mathrm{},n_t,h_t=t_{max}/n_t,$$

where $`n_x`$ and $`n_t`$ are numbers of partitions.

The following difference scheme with weights $`\gamma `$ was implemented (see  for details):

$$e_j^k\frac{T_j^{k+1}T_j^k}{h_t}=k_0\left[\gamma \frac{T_{j+1}^{k+1}2T_j^{k+1}+T_{j1}^{k+1}}{h_x^2}+(1\gamma )\frac{T_{j+1}^k2T_j^k+T_{j1}^k}{h_x^2}\right]+q_j^{k+\frac{1}{2}},$$
(9)
where

$$T_j^k=T(x_j,t_k),e_j^k=\rho (T_j^k)c(T_j^k),q_j^{k+\frac{1}{2}}=q(x_j,t_k+\frac{h_t}{2}),$$

and the upper index numerates different moments of time (time “levels”), the lower one specifies a set of spatial co-ordinates. The scheme is absolutely convergent at $`\gamma =0.5`$ and possesses the second-order accuracy for both variables.

From initial condition $`T(x,0)=T_0`$, values $`T_j^0`$ $`(j=0,\mathrm{},n_x)`$ on a zero time level are known. The boundary conditions

$$\frac{T_1^kT_1^k}{2h_x}=\frac{T_{n_x+1}^kT_{n_x1}^k}{2h_x}=0,k=1,\mathrm{},n_t$$

allow one to introduce symmetric points $`x_1=h_x`$ and $`x_{n_x+1}=1+h_x`$ with appropriate values $`T_1^k=T_1^k`$ and $`T_{n_x+1}^k=T_{n_x1}^k`$ respectively. So, we can use the Eq. (9) in points $`x_0`$ and $`x_{n_x}`$. Using initial and boundary conditions, we obtain a system of $`n_x`$ linear algebraical equations with the same number of variables. Thus, under the accepted approximation, we reduced the partial differential equation (8) to system (9) of linear algebraic equations. The matrix of this system is tridiagonal and after its solution333Recursive relations for determining the solution of algebraic problem (9) comprise the well-known sweep method, called also forward-backward or Thomas algorithm .  we obtain value $`T_j^1`$ $`(j=0,\mathrm{},n_x)`$ at the first time level. Repeating this process, values $`T_j^k`$ on every time level $`k`$ are computed.

The result of straightforward verification of the Stefan condition (1) is shown in Fig. 5, where the function

$$\varphi (t)=k\left(\frac{T}{x}|_{T_A}\frac{T}{x}|_{T_B}\right)\lambda \frac{d\xi }{dt}$$
(10)
is depicted.

The left and right points, in which the spatial derivatives of temperature were taken in (10), are shown in Fig. 4. They define a spatial layer which nearly the whole fusion energy is absorbed within. From Fig. 5 one can see that condition (1) is satisfied indeed, but only after a characteristic relaxation time $`𝐭_\mathrm{𝟏}`$ has elapsed. The physical meaning of $`t_1`$ is clear from Fig. 6. Namely, it corresponds to the transition from a rapid to slow motion of the exterior interphase surface. In the case when boundary motion is rapid, the heat necessary for fusion is brought into the melting layer directly from the external source $`q(x,t).`$ The slow motion corresponds to the ordinary Stefan mode when the process is controlled mainly by the heat entered into the layer through its boundary.

It is also seen from Fig. 6 that transition to the Stefan mode takes place earlier than the external source to be totally turned off:

$$t_1<\tau .$$

Time $`t_2`$ shown in Fig. 6 denotes a moment when the thickness of the melted material begins to diminish due to heat escape into a more cooler solid phase.

Fig. 7 and 8 also confirm the conclusions which we have come to in the previous section. Formation of the “tableland” (whose height corresponds to the fusion temperature) for spatial temperature distribution is distinctly seen in Fig. 7. The narrow strip restricted by two dashed lines in Figs. 7 and 8 exhibits the width of the smoothed $`\delta `$-function. We believe that existence of two breaks for the spatial derivative is masked in Fig. 7 with this $`\delta `$-function smearing. Fig. 8 demonstrates a temperature evolution for two divorced spatial points. One can make sure that the above mentioned time interval corresponding to the same temperature at the different spatial points really exists. It is evident that such a behavior of temperature has nothing to do with the traditional description in the framework of (1) and (2).

Fig. 9 shows a time-dependence of the interphase coordinate. Numbers 1 and 2 denote the regions where verification of the Stefan condition (1) is impossible due to $`𝚫`$-instability. It means that small variations of fusion temperature value, $`T^{}`$, lead to a drastic change of interphase position (see dotted lines in Fig. 9).

## 4 Track formation in solids

The next example demonstrating the preference for the $`\delta `$-function approach is connected with the problem of track formation in solids. In fact, at present nobody knows with certainty the main mechanism responsible for these processes. Furthermore, it seems like the universal model explaining all of them does not exist and different materials behave differently under heavy ion attack. Here we assume the so-called thermal spike model based on the following system of two coupled nonlinear differential equations (see, e.g.  and references therein):

$$\rho C_e(T_e)\frac{T_e}{t}=\frac{1}{r}\frac{}{r}\left[rK_e(T_e)\frac{T_e}{r}\right]g(T_eT_i)+q(r,t),$$
(11)
$$\rho C_i(T_i)\frac{T_i}{t}=\frac{1}{r}\frac{}{r}\left[rK_i(T_i)\frac{T_i}{r}\right]+g(T_eT_i),$$
(12)
where $`T_e`$ and $`T_i`$ are electrons and lattice temperatures, respectively, $`C_e`$$`C_i`$ and $`K_e,K_i`$ specific heat and thermal conductivity for the electronic system and lattice, $`\rho `$ is the material density, $`g`$ the electron-atom coupling, $`q(r,t)`$ the power brought on the electronic system, $`r`$ the radius in cylindrical geometry with the ion path as the axis. One can see that electrons receive their energy directly from the external source $`q(r,t)`$ which takes into account ion energy loss in electron gas. The characteristic duration of source activity is usually in the range $`10^{15}5\times 10^{15}`$ s. According to (12), atoms are heated due to electron-atom coupling represented by the term $`g(T_eT_i)`$. Nuclear interaction of atoms with the projectile ion is relatively small and, therefore, can be neglected. It is clear that coupling is the most effective at the beginning of the relaxation process when $`T_eT_i`$ and $`g(T_eT_i)gT_e`$.

The initial conditions can be chosen in a form

$$T_e(r,0)=T_i(r,0)=T_0,$$

and the boundary ones444 One should take into account that there is no heat transfer at the center of track. can be taken as

$$\left(\frac{T_e}{r}\right)_{r=0}=\left(\frac{T_i}{r}\right)_{r=0}=0,T_e(r_{max},t)=T_i(r_{max},t)=T_0,$$

where $`r_{max}`$ was taken of order $`10^5cm`$.

The thermal spike model explains track formation as a structural transition of lattice due to its warming-up and melting with subsequent quenching. These processes are usually accompanied with disorder creation in the lattice. Indeed, rapid quenching leads to a “conservation” of atoms’ random places that were in the melted material just before cooling. For amorphous materials, which are characterized by high disorder of atoms’ positions and small values of thermal conductivity, quenching, quite the contrary, leads to putting atoms’ places in order. But in either case, structural modifications are observed in the microscope as an ion trace in solid.

Besides thermal spike, one may assume the ion spike as well, when the track is formed due to the electrostatic repulsion of ionized atoms. The main reason justifying our utilization of system (11), (12) is an agreement of nuclear track radii, calculated in this framework, with the experimental data . The total formulation of the model includes many physical details, such as a description of the source $`q(x,t)`$, and is outside the scope of this publication. Here we only touch some problems concerning the main topic of the paper.

A numerical algorithm similar to that described above has been elaborated for numerical solving system (11), (12). The radial distribution of the lattice temperature $`T_i`$ around the path of $`\mathrm{𝐏𝐛}`$ in amorphous $`\mathrm{𝐆𝐞}`$ at kinetic energy of impinging ions of about $`110`$ MeV is shown in Fig. 10 (for two different moments of time). One can see the typical “tablelands” similar to those discussed in the previous chapter and which could not be obtained in frames of the classical Stefan approach. Our calculations show that the “tableland” exists here only during a transitory time $`t_1`$, when the material is under a strong exposure of the source $`g(T_eT_i)gT_e`$. It is shorter than $`\tau =\varrho C_e/g10^{12}`$ s (see Fig. 10), where $`\tau `$ is electron-atom relaxation time555The formula for $`\tau `$ estimation follows from Eq. (11)., characterizing duration of source activity in (12) (compare this conclusion with data presented in Fig. 6).

It is interesting to note that there is a real, met in the nature, “regularization” of $`\delta `$-function analogous to that implemented in this paper. For materials with a complex molecular structure (high temperature superconductors, biological molecules, alloys etc.), the melting temperature is not fixed but, instead, smeared within a characteristic interval where atom bonds of different type are gradually destroyed with temperature increase. In this case the only possible approach to the problem should be based on the condition (3). An approximate $`\delta `$-function, analogous to that shown in Fig. 4, can be extracted here directly from the experiment. In , a model based on the smeared $`\delta `$-function approach and Eqs. (11), (12) was used for computation of effective electron-atom relaxation time $`\tau `$ in a high temperature superconductor. The established $`\tau `$ turned out to be in a good agrement with experimentally observed values.

## 5 Conclusion

To the best of our knowledge, the peculiarities of phase transition dynamics, we discussed in this paper, have never been considered explicitly in mathematical physics. This fact may be explained partially by the necessity to use very powerful spatially distributed external sources of heat, in order the above mentioned effects to be urgent. Such sources were hardly available for industrial applications not long ago. However, the examples which given above are likely evidences of the fact that such sources, “interfering” in the thermal conductive processes, are integral parts of all most recent ion beam technologies. The numerical investigations, which have been undertaken, show that the $`\delta `$-function approach to phase transitions is a suitable instrument to tackle these problems, though the authors of this idea have never used it in such a context.

Within networks of scientific papers devoted to the interphase motion problem, one can distinguish in retrospect the following logical order: formulation of the Stefan problem (Lame and Clapeiron; Stefan) $``$ application of the $`\delta `$-function approach for equivalent representation of it (Tikhonov and Samarskii) $``$ numerical implementations of the idea (Samarskii with co-authors; this review) $``$ description of materials with interval distributed fusion temperature (a natural physical interpretation of the previous step). Here we have shown that it is more expedient to turn over this order and take its last element as the basis for solving both the classical Stefan problem and a more general one including spatially distributed sources. In other words, both of these solutions could be considered as an idealized limiting case $`\mathrm{\Delta }0`$ of a natural physical point of view that none of phase transitions take place at the exactly defined value of fusion temperature. A peculiarity of a new, found in this paper, solution to Eq. (3) is its “tableland” behavior seen in Figs. 1(a), 7, 8 and 10. Such solutions could never be obtained in the framework of the classical Stefan formulation (compare with Fig. 1(b)).

We would like to express our gratitude to professor M. Toulemonde from CIRIL (Caen) for interest and useful discussion. The authors are grateful to doctors E. Airjan, I. Amirkhanov, I.N. Goncharov and T.P. Puzynina from JINR (Dubna) for stimulating discussions.

This investigation has been supported in part by the Russian Foundation for Basic Research, project No. 02-01-00606.


format: text

text:
# Einstein boundary conditions for the 3+1 Einstein equations

## Abstract

In the 3+1 framework of the Einstein equations for the case of vanishing shift vector and arbitrary lapse, we calculate explicitly the four boundary equations arising from the vanishing of the projection of the Einstein tensor along the normal to the boundary surface of the initial-boundary value problem. Such conditions take the form of evolution equations along (as opposed to across) the boundary for certain components of the extrinsic curvature and for certain space-derivatives of the three-metric. We argue that, in general, such boundary conditions do not follow necessarily from the evolution equations and the initial data, but need to be imposed on the boundary values of the fundamental variables. Using the Einstein-Christoffel formulation, which is strongly hyperbolic, we show how three of the boundary equations up to linear combinations should be used to prescribe the values of some incoming characteristic fields. Additionally, we show that the fourth one imposes conditions on some outgoing fields.

The advent of sensitive gravitational wave detectors such as LIGO lig  and LISA lis  motivates the need for accurate templates of gravitational wave emission from powerful sources such as the collisions of black holes with other black holes or with neutron stars Finn and Thorne (2000). Even though the early and late stages of such merging events can be modelled analytically by perturbative methods Price and Pullin (1994); Gleiser et al. (1996), the actual merger phase is highly non-linear and can only be expected to be modelled accurately by numerical simulations. Most numerical simulations of dynamical black-hole spacetimes break down well before the relevant information can be extracted, for reasons that are yet to be pinned down. Numerical analysis on the side, there are many fundamental issues in the Einstein equations that may be relevant to the stability of any numerical implementation of the binary black hole merger, among them: intrinsic ill-posedness of the related Cauchy problem, constraint violations, poor physical choices of binary black hole data and poor choices of boundary conditions. The solution to long term numerical instabilities of binary black hole mergers is likely to involve an orchestrated control of all these factors, in addition to the other elements of general numerical analysis that apply to time-dependent problems with boundaries. In this work, we point out that even before thinking about a numerical integration of the Einstein equations in a region with timelike boundaries, the possible implications of the Einstein equations to the boundaries need to be taken into consideration, as a matter of principle. To date, most numerical simulations with boundaries have largely disregarded this issue, as have most mathematical studies of the mixed initial-boundary value problem of the Einstein equations, with exceptions explicitly noted in the following.

Consider the 3+1 formulation of the vacuum Einstein equations with vanishing shift vector, in which the metric has the form

$$ds^2=\alpha ^2dt^2+\gamma _{ij}dx^idx^j$$
(1)
where $`\alpha `$ is the lapse function and $`\gamma _{ij}`$ is the Riemannian metric of the spatial three-slices of fixed value of $`t`$. The Einstein equations $`G_{ab}=0`$ take the ADM form York (1979)

$`\dot{\gamma }_{ij}`$ $`=`$ $`2\alpha K_{ij},`$ (2)
$`\dot{K}_{ij}`$ $`=`$ $`\alpha (R_{ij}2K_{il}K^l{}_{j}{}^{}+KK_{ij})D_iD_j\alpha ,`$ (3)
with the constraints

$`RK_{ij}K^{ij}+K^2`$ $`=`$ $`0,`$ (4)
$`D_jK^j{}_{i}{}^{}D_iK`$ $`=`$ $`0.`$ (5)
Here an overdot denotes a partial derivative with respect to the time coordinate $`/t`$, indices are raised with the inverse metric $`\gamma ^{ij}`$, $`D_i`$ is the covariant three-derivative consistent with $`\gamma _{ij}`$ , $`R_{ij}`$ is the Ricci curvature tensor of $`\gamma _{ij}`$, $`R`$ its Ricci scalar and $`K\gamma ^{ij}K_{ij}`$. Because there are no time derivatives of the lapse in any of the equations, the lapse function is completely free. Because a time derivative of the constraints turns out to be a linear combination of the constraints themselves and their space derivatives by virtue of the evolution equations, they will vanish at subsequent times in the domain of dependence of the initial slice if vanishing initially (how stable the propagation of the zero values of the constraints is to small perturbations does not concern us at this time Frittelli (1997)).

Strictly speaking, the initial-value problem consists of finding a spatially periodic or square integrable solution to the Einstein equations in the half-spacetime $`t0`$. In practice, one looks for a spatially non-periodic solution in $`t0,x^i𝒱`$, where $`𝒱`$ is a box in $`R^3`$. There arises the problem of specifying data on the timelike boundaries of the region of interest. The problem is rooted in the fact that in order to specify such data one needs to know the values of the solution on the boundary for all time, which, of course, one does not know. In an initial value problem without constraints, the boundary values may or not be constrained by the choice of initial data. It seems to us that if the problem, in addition, has constraints, the boundary values will know it, and it should be natural to expect consequences to the boundary values from the presence of the constraints. The point has been raised in a precise sense by Stewart Stewart (1998) and further pursued in a similar direction by Calabrese et al Calabrese et al. (2002) and by Szilágyi and Winicour Szilágyi and Winicour .

The existence of at least some such consequences is intuitively predictable from the point of view that the 3+1 equations derive from a 4-d covariant formulation, which is more fundamental, that is: $`G_{ab}=0`$. The 3+1 split is simply a choice of linearly independent combinations of the ten components of the Einstein tensor. The choice is precisely that which removes any time derivatives of second order from four of the equations, and is accomplished by projecting the Einstein tensor along the direction normal to the slicing $`G_{ab}n^b=0`$.

One can easily convince oneself (or else refer to Frittelli and Gómez (2003)) that the projection of the Einstein tensor along the direction normal to a surface given by a fixed value of any of the coordinates $`x^i`$ removes all second derivatives with respect to that coordinate $`x^i`$. Let’s say we look at our boundary at a fixed value of $`x`$ and write down the projection of the Einstein tensor of the metric (1) along the normal vector $`e^b=(0,\gamma ^{xi})`$, which may be transformed to unit length by dividing by $`\sqrt{\gamma ^{xx}}`$, but that is clearly not necessary. The four components of the projection of the Einstein tensor along $`e^b`$ are

$`G_{tb}e^b`$ $`=`$ $`{\displaystyle \frac{1}{2}}\gamma ^{ix}((\mathrm{ln}\gamma ),_{it}\gamma ^{kl}\dot{\gamma }_{ik,l})KD^x\alpha +K_k^xD^k\alpha `$ (6a)
$`+\alpha \left(\gamma ^{kl}\mathrm{\Gamma }_{kl}^jK_j^x+\gamma ^{ix}\mathrm{\Gamma }_{ik}^jK_j^k\right)=0`$
$`G_{yb}e^b`$ $`=`$ $`{\displaystyle \frac{\dot{K}_y^x}{\alpha }}+KK_y^x+R_y^x{\displaystyle \frac{1}{\alpha }}D^xD_y\alpha =0`$ (6b)
$`G_{zb}e^b`$ $`=`$ $`{\displaystyle \frac{\dot{K}_z^x}{\alpha }}+KK_z^x+R_z^x{\displaystyle \frac{1}{\alpha }}D^xD_z\alpha =0`$ (6c)
$`G_{xb}e^b`$ $`=`$ $`{\displaystyle \frac{\dot{K}\dot{K}_x^x}{\alpha }}{\displaystyle \frac{1}{2}}(R+K^{ij}K_{ij}+K^2)+KK_x^x`$ (6d)
$`+R_x^x+{\displaystyle \frac{1}{\alpha }}\left(D^jD_j\alpha D^xD_x\alpha \right)=0`$
Here the time derivative of the components of the extrinsic curvature is applied after raising an index, that is: $`\dot{K}_j^i(\gamma ^{ik}K_{kj}),_t`$. The reader can verify that $`R_y^x`$ and $`R_z^x`$ do not involve second derivatives with respect to $`x`$ of any of the variables and that the combination $`R_x^x\frac{1}{2}R`$ doesn’t either. In the framework of the ADM equations (2), even though, to our knowledge, there is no clear notion of which kind of boundary conditions would be appropriate for a unique solution, one could think of (6a) as an evolution equation for a particular combination of the first space-derivatives of the metric, of Eqs. (6b) and (6c) as evolution equations for $`K_y^x`$ and $`K_z^x`$ respectively, and of Eq. (6d) as an evolution equation for the combination $`K_y^y+K_z^z`$. Or one could freely combine them linearly to obtain evolution equations for combinations of these fundamental variables. In any case, these are evolution equations along the boundary itself which can be viewed as a sort of mixed Neumann-Dirichlet boundary conditions Szilágyi et al. (2000).

The situation clarifies in the case that one considers a first-order reduction of the ADM equations. In that case, the space-derivatives of the metric are adopted as additional variables $`q_{ijk}\gamma _{ij,k}`$. This has the effect of eliminating all appearances of $`x`$derivatives of the fundamental variables from Eqs (6). However, a straight reduction in this fashion is not well posed in the sense that it is not strongly hyperbolic. If one is to consider a first-order reduction, one is much better off adopting one that is strongly hyperbolic in the sense that there exists a complete set of null eigenvectors of the principal symbol of the resulting system of equations.

In a generic strongly hyperbolic system of PDE’s, the null eigenvectors constitute characteristic fields that propagate with the characteristic speeds of the system, and the notion of consistent boundary conditions can be defined unambiguously, as follows Gustaffson et al. (1995). Characteristic fields that are outgoing at the boundary are determined by the initial values and cannot be prescribed freely at the boundary. On the other hand, characteristic fields that are incoming at the boundary must be prescribed in order for there to be a unique solution, and, furthermore, can be prescribed arbitrarily (we do not concern ourselves with stability for the moment).

In the case that there exist constraints, it seems to us that this notion of consistent boundary conditions may still be used but must be adapted to the particularities of a constrained system. To wit: if the incoming characteristic fields are to lead to a solution that satisfies the constraints inside the region of interest, then they must “know” about the constraints and, in general, may not be prescribed as freely as in the case where there are no constraints. In the case of the Einstein equations, we have shown that there are four boundary equations implied by the ten Einstein equations themselves, which are completely equivalent to linear combinations of the evolution equations and constraints that have no second derivatives across the boundary. Moreover, in a sense, the boundary equations are to the boundary surface what the constraints are to the initial slice. To see this consider the Bianchi identities $`_aG^{ab}=0`$, which take the following form in a generic coordinate system $`(t,x^i)`$ adapted to spacelike slices:

$$_xG^{xb}=_tG^{tb}_yG^{yb}_zG^{zb}{}_{}{}^{4}\mathrm{\Gamma }_{}^{a}{}_{ac}{}^{}G_{}^{cb}{}_{}{}^{4}\mathrm{\Gamma }_{}^{b}{}_{ac}{}^{}G_{}^{ac}$$
(7)
where $`{}_{}{}^{4}\mathrm{\Gamma }_{}^{c}_{ab}`$ are the Christoffel symbols of the four-metric $`g_{ab}`$. The right-hand side obviously contains no derivatives with respect to $`x`$ higher than second order. Therefore, $`_xG^{xb}`$ must be of second order in $`_x`$ at most, which implies that $`G^{xb}`$ involve, at most, first-order derivatives with respect to $`x`$. Of course, the reasoning runs as well for the $`y`$derivatives, the $`z`$derivatives, and the $`t`$derivatives, so that $`G^{ay}`$ has not second $`y`$derivatives, $`G^{az}`$ has no second $`z`$derivatives and $`G^{tb}`$ have no second time derivatives. In fact, traditionally this reasoning has been implicitly used to argue that the constraints $`G^{tb}=0`$ are the only conditions needed at the initial slice to specify a solution to the Einstein equations via the initial-value problem Weinberg (1971).

We conclude that $`G_{ab}e^bg_{ac}G^{cx}=0`$ contain as much knowledge about the solution of the Einstein equations in the interior as the boundary values may be expected to “know”. For this reason, we think that imposing the boundary conditions (6) on the incoming characteristic fields is necessary in order to evolve a solution to the evolution equations from constrained initial data that will satisfy the four constraints at a later time slice as well.

To illustrate the point, we now consider the Einstein-Christoffel formulation of the 3+1 equations Anderson and York (1999). By defining first-order variables as the following 18 linearly independent combinations of the space-derivatives of the metric:

$$f_{kij}\mathrm{\Gamma }_{(ij)k}+\gamma _{ki}\gamma ^{lm}\mathrm{\Gamma }_{[lj]m}+\gamma _{kj}\gamma ^{lm}\mathrm{\Gamma }_{[li]m}$$
(8)
where $`\mathrm{\Gamma }^k_{ij}`$ are the Christoffel symbols of $`\gamma _{ij}`$, and choosing the lapse function as $`\alpha Q\sqrt{\gamma }`$ with $`Q`$ assumed arbitrarily prescribed a priori, the 3+1 equations can be put in the following equivalent form Kidder et al. (2000):

$`\dot{\gamma }_{ij}`$ $`=`$ $`2\alpha K_{ij}`$ (9a)
$`\dot{K}_{ij}+\alpha \gamma ^{kl}_lf_{kij}`$ $`=`$ $`\alpha \{\gamma ^{kl}(K_{kl}K_{ij}2K_{ki}K_{lj})+\gamma ^{kl}\gamma ^{mn}(4f_{kmi}f_{[ln]j}+4f_{km[n}f_{l]ij}f_{ikm}f_{jln}`$ (9b)
$`+8f_{(ij)k}f_{[ln]m}+4f_{km(i}f_{j)ln}8f_{kli}f_{mnj}+20f_{kl(i}f_{j)mn}13f_{ikl}f_{jmn})`$
$`_i_j\mathrm{ln}Q_i\mathrm{ln}Q_j\mathrm{ln}Q+2\gamma _{ij}\gamma ^{kl}\gamma ^{mn}(f_{kmn}_l\mathrm{ln}Qf_{kml}_n\mathrm{ln}Q)`$
$`+\gamma ^{kl}[(2f_{(ij)k}f_{kij})_l\mathrm{ln}Q+4f_{kl(i}_{j)}\mathrm{ln}Q3(f_{ikl}_j\mathrm{ln}Q+f_{jkl}_i\mathrm{ln}Q)]\}`$
$`\dot{f}_{kij}+\alpha _kK_{ij}`$ $`=`$ $`\alpha \{\gamma ^{mn}[4K_{k(i}f_{j)mn}4f_{mn(i}K_{j)k}+K_{ij}(2f_{mnk}3f_{kmn})]`$ (9c)
$`+2\gamma ^{mn}\gamma ^{pq}[K_{mp}(\gamma _{k(i}f_{j)qn}2f_{qn(i}\gamma _{j)k})+\gamma _{k(i}K_{j)m}(8f_{npq}6f_{pqn})`$
$`+K_{mn}(4f_{pq(i}\gamma _{j)k}5\gamma _{k(i}f_{j)pq})]K_{ij}_k\mathrm{ln}Q`$
$`+2\gamma ^{mn}(K_{m(i}\gamma _{j)k}_n\mathrm{ln}QK_{mn}\gamma _{k(i}_{j)}\mathrm{ln}Q)\}`$
with the constraints

$`𝒞`$ $``$ $`\gamma ^{ij}\gamma ^{kl}\{2(_kf_{ijl}_if_{jkl})+K_{ik}K_{jl}K_{ij}K_{kl}`$ (10a)
$`+\gamma ^{mn}[f_{ikm}(5f_{jln}6f_{ljn})+13f_{ikl}f_{jmn}`$
$`+f_{ijk}(8f_{mln}20f_{lmn}]\}=0`$
$`𝒞_i`$ $``$ $`\gamma ^{kl}\{\gamma ^{mn}[K_{ik}(3f_{lmn}2f_{mnl})K_{km}f_{iln}]`$ (10b)
$`+_iK_{kl}_kK_{il}\}=0`$
$`𝒞_{kij}`$ $``$ $`_k\gamma _{ij}2f_{kij}+4\gamma ^{lm}(f_{lm(i}\gamma _{j)k}\gamma _{k(i}f_{j)lm})=0`$
to be imposed only on the initial data. Here $`𝒞`$ and $`𝒞_i`$ are the Hamiltonian and momentum constraints, respectively, and $`𝒞_{ijk}`$ represent the definition of the additional 18 first-order variables $`f_{kij}`$ needed to reduce the ADM equations to full first-order form (they represent (8) inverted for $`\gamma _{ij,k}`$ in terms of $`f_{kij}`$). No other constraints are needed to pick a solution of the Einstein equations out of the larger set of solutions of the Einstein-Christoffel equations.

In this language, the boundary conditions (6b-6d) require the derivatives of the lapse $`\alpha `$ and the components $`R_y^x,R_z^x`$ and $`R_x^x(1/2)R`$ expressed in terms of $`f_{kij}`$, the point being that there is a way to express their occurrence so that no $`x`$derivatives of any $`f_{kij}`$ appears. As they stand, they remain as evolution equations for the components $`K_y^x,K_z^x`$ and $`K_y^y+K_z^z`$ respectively. However, (6a) takes the following form

$`G_{tb}e^b`$ $`=`$ $`\dot{f}^x{}_{m}{}^{}{}_{}{}^{m}\dot{f}^m{}_{m}{}^{}{}_{}{}^{x}+\alpha \{2K_l^x(4f^l{}_{m}{}^{}{}_{}{}^{m}3f^m{}_{m}{}^{}{}_{}{}^{l})`$ (11)
$`+K_l^mf^x{}_{m}{}^{}{}_{}{}^{l}K^x\alpha +K^{xl}_l\alpha \}=0`$
where $`\dot{f}^k{}_{i}{}^{}{}_{}{}^{j}_t(\gamma ^{km}f_{mil}\gamma ^{jl})`$. This is an evolution equation for a particular linear combination of the components of $`f_{kij}`$.

The next step is to figure out which variables these four equations should prescribe at the boundary. We start by recalling the characteristic fields of the problem with respect to the unit vector $`\xi ^i\gamma ^{xi}/\sqrt{\gamma ^{xx}}`$ which is normal to the boundary $`x=x_0`$ for the region $`xx_0`$. In this case, we have that

$${}_{}{}^{+}U_{i}^{j}K_i^j+\frac{f^x{}_{i}{}^{}^j}{\sqrt{\gamma ^{xx}}}$$
(12)
are six characteristic fields that are outgoing at the boundary, and

$${}_{}{}^{}U_{i}^{j}K_i^j\frac{f^x{}_{i}{}^{}^j}{\sqrt{\gamma ^{xx}}}$$
(13)
are another six characteristic fields that are incoming at the boundary. The remaining characteristic fields are the six components of $`f^y{}_{i}{}^{}^j`$, the six components of $`f^z{}_{i}{}^{}^j`$ and the six components of $`\gamma _{ij}`$, all 18 of which travel upwards along the boundary (they have vanishing characteristic speeds).

Inverting to find the fundamental fields in terms of the characteristic fields we have

$`K_i^j`$ $`=`$ $`{\displaystyle \frac{{}_{}{}^{+}U_{i}^{j}+{}_{}{}^{}U_{i}^{j}}{2}}`$ (14)
$`f^x{}_{i}{}^{}^j`$ $`=`$ $`{\displaystyle \frac{{}_{}{}^{+}U_{i}^{j}{}_{}{}^{}U_{i}^{j}}{2\sqrt{\gamma ^{xx}}}}`$ (15)
One can immediately see that the boundary conditions (6b) and (6c) turn into evolution equations for $`{}_{}{}^{}U_{y}^{x}`$ and $`{}_{}{}^{}U_{z}^{x}`$ assuming the outgoing fields $`{}_{}{}^{+}U_{y}^{x}`$ and $`{}_{}{}^{+}U_{z}^{x}`$ are known, respectively. In fact, these equations yield the values of $`{}_{}{}^{}U_{y}^{x}`$ and $`{}_{}{}^{}U_{z}^{x}`$ linearly in terms of the outgoing fields and in terms of sources known from the previous time step. So the transverse components $`a=y,z`$ of the projection $`G_{ab}e^b=0`$ are quite welcome and are actually necessary and consistent with the initial-boundary value problem of the Einstein equations.

Turning now to the longitudinal component $`a=x`$ of $`G_{ab}e^b=0`$, which is Eq. (6d), we see that, because $`KK_x^x=K_y^y+K_z^z`$, the equation takes the form

$$({}_{}{}^{+}\dot{U}_{y}^{y}+{}_{}{}^{+}\dot{U}_{z}^{z})+({}_{}{}^{}\dot{U}_{y}^{y}+{}_{}{}^{}\dot{U}_{z}^{z})+\mathrm{}=0.$$
(16)
On the other hand, because

$$\dot{f}^x{}_{m}{}^{}{}_{}{}^{m}\dot{f}^m{}_{m}{}^{}{}_{}{}^{x}=\dot{f}^x{}_{y}{}^{}{}_{}{}^{y}+\dot{f}^x{}_{z}{}^{}{}_{}{}^{z}\dot{f}^y{}_{y}{}^{}{}_{}{}^{x}\dot{f}^z{}_{z}{}^{}^x$$
(17)
then the time component $`a=t`$ of $`G_{ab}e^b=0`$, which is Eq. (11), takes the form

$$({}_{}{}^{+}\dot{U}_{y}^{y}+{}_{}{}^{+}\dot{U}_{z}^{z})({}_{}{}^{}\dot{U}_{y}^{y}+{}_{}{}^{}\dot{U}_{z}^{z})\dot{f}^y{}_{y}{}^{}{}_{}{}^{x}\dot{f}^z{}_{z}{}^{}{}_{}{}^{x}+\mathrm{}=0.$$
(18)
Both equations thus involve the same incoming and outgoing characteristic fields: $`{}_{}{}^{}U_{y}^{y}+{}_{}{}^{}U_{z}^{z}`$ and $`{}_{}{}^{+}U_{y}^{y}+{}_{}{}^{+}U_{z}^{z}`$. This means that they determine both the outgoing and incoming fields $`{}_{}{}^{+}U_{y}^{y}+{}_{}{}^{+}U_{z}^{z}`$ and $`{}_{}{}^{}U_{y}^{y}+{}_{}{}^{}U_{z}^{z}`$. One way (by no means the only one) to make this point explicit is the following. Subtracting the two equations with appropriate factors one obtains an evolution equation for $`{}_{}{}^{}U_{y}^{y}+{}_{}{}^{}U_{z}^{z}`$, which is welcome as a prescription of boundary values that are consistent with the Einstein equations:

$$2\alpha G_{xb}e^b\frac{2}{\sqrt{\gamma ^{xx}}}G_{tb}e^b={}_{}{}^{}\dot{U}_{y}^{y}+{}_{}{}^{}\dot{U}_{z}^{z}+\mathrm{}=0,$$
(19)
which also involves time derivatives of some of the fields that travel along the boundary, but this is fine. In fact, we have shown already Frittelli and Gómez (2003) that at least in the case of the Einstein-Christoffel formulation with spherical symmetry, this resulting boundary evolution equation is being used in numerical integration and is referred to as a constraint-preserving boundary condition Calabrese et al. (2002).

So far, the projection of the Einstein equations in the direction normal to the boundary is seen as providing prescriptions for the boundary values that are necessary in the light of the initial-boundary value problem of this strongly hyperbolic formulation. Any solution to the ten Einstein equations must satisfy such boundary equations. Therefore, failure to impose such conditions on the boundary would result in a solution to the evolution equations with constrained initial data which would not solve all four of the constraints at a later time slice.

There remains the question of the fourth boundary equation implied by adding $`G_{tb}e^b=0`$ and $`G_{xb}e^b=0`$ with appropriate factors in order to eliminate the occurrence of the incoming fields. This equation is linearly independent of the other three, therefore it will not be satisfied as a consequence of imposing the other three. One can see by inspection that this is an evolution equation for the outgoing field $`{}_{}{}^{+}U_{y}^{y}+{}_{}{}^{+}U_{z}^{z}`$:

$$2\alpha G_{xb}e^b+\frac{2}{\sqrt{\gamma ^{xx}}}G_{tb}e^b={}_{}{}^{+}\dot{U}_{y}^{y}+{}_{}{}^{+}\dot{U}_{z}^{z}+\mathrm{}=0,$$
(20)
which involves also the time derivatives of some of the fields that travel along the boundary, but it does not involve the time derivative of any incoming characteristic fields. This equation is not “necessary” in the sense that the outgoing fields are already determined by propagation from their initial values. However, this equation must be consistent with the outgoing fields determined by the initial values, otherwise one would not be solving the full set of Einstein equations. Therefore, it is actually necessary for this equation to be satisfied. In fact, in principle, one may simply use this equation for those outgoing characteristic fields involved, with the confidence that this boundary prescription is consistent with the initial values, so long as the initial values satisfy the constraints with sufficient accuracy. The actual numerical implementation of this boundary condition is a subject for subsequent study and lies beyond our current scope.

In summary, in the case of the Einstein-Christoffel formulation, the four “Einstein boundary conditions” defined as $`G_{ab}e^b=0`$ up to linear combinations constrain the values of three of the six incoming characteristic fields. The other three incoming characteristic fields are free, since there are no other available conditions on the boundary. In fact, the 18 first-order constraints $`𝒞_{ijk}`$ –which are needed in order to pick an actual solution of the Einstein equations out of the larger set of solutions of the Einstein-Christoffel equations– are preserved trivially by the evolution. This can be seen easily by taking a time derivative of $`𝒞_{ijk}`$ and using the evolution equations to eliminate the occurrence of time derivatives of the fundamental variables in the right-hand side. The reader may verify that this leads to $`\dot{𝒞}_{ijk}=F_{ijk}(C,C_i,C_{ijk})`$, that is: ordinary differential equations for the propagation of $`𝒞_{ijk}`$. This means that the 18 first-order constraints are trivially satisfied on the boundary by virtue of the initial data and the evolution equations, and, consequently, they may not prescribe conditions on the incoming fields. There thus remain three incoming fields that can be specified in any arbitrary fashion for a solution to the Einstein equations. They may, for instance, be “frozen” by setting their time derivatives to zero, as is done by some authors Kidder et al. (2000, 2001).

We have thus obtained a complete set of boundary conditions for the Einstein-Christoffel formulation from the projection of the Einstein tensor normally to the boundary and the fact that the 18 constraints that realize the first-order reduction are preserved by ordinary differential equations. The incoming fields that are prescribed by the three boundary equations carry, in a sense, the vanishing of $`G_{ab}e^b`$ into the interior, where the evolution equations hold. Since $`G_{ab}e^b`$ are linear combinations of the evolution equations and the constraints, it follows that the constraints will be satisfied in the interior. Thus the boundary conditions preserve the constraints.

The argument makes no use of the auxiliary system of propagation of the constraints implied by the evolution equations. Still, the resulting boundary conditions appear to be completely consistent with constraint propagation in the following sense. In the three-dimensional case that interests us here, if the evolution equations are satisfied, the constraint functions propagate according to a strongly hyperbolic system Kidder et al. (2001). In this auxiliary system, three characteristic constraints are incoming, three are outgoing and the rest are “static”. In the context of this auxiliary system, the three incoming constraints must be set to zero along the boundary. One may thus infer that there will be three conditions on the boundary values of the fundamental fields that are necessary and sufficient to preserve the constraints. Since we have three Einstein boundary conditions and no other conditions on the boundary, the Einstein boundary conditions and the three conditions that arise from constraint propagation must be equivalent in some definite sense. We argue that the equivalence is through linear combinations of the evolution equations. In the spherically symmetric case we have shown that this is precisely the case Frittelli and Gómez (2003). Work is currently in progress to demonstrate explicitly the equivalence of the two sets of conditions in the three dimensional case and will appear elsewhere. Interestingly, since constraint propagation is at the basis of Calabrese et al’s “constraint-preserving” boundary conditions Calabrese et al. , it would thus appear that projecting the Einstein equations normally to the boundary would shortcut the procedure leading to “constraint-preserving” boundary conditions, providing an exact nonlinear version of the conditions and a geometric basis for them at the same time.

We emphasize that there is a very subtle issue of consistency between the constrained initial data traveling towards the boundary, and the conditions $`G_{ab}e^b=0`$ on the boundary that apply to the outgoing fields. The issue turns critical in numerical simulations because of the inaccuracies involved in satisfying the initial constraints. Clearly, the outgoing fields propagating out from the initial data will carry and perhaps increase the effects of the initial errors, and there is a good chance that they will fail to satisfy (20), the odds of failure increasing with time. This issue may have large consequences for achieving long time numerical evolution. This observation applies especially to numerical simulations implementing black hole excision Shoemaker et al. . As far as we know, no consistency issues in strongly hyperbolic systems with constraints have been addressed in the reference literature on finite-difference methods such as Gustaffson et al. (1995).

Interestingly, the other constrained system of equations in physics which is oftentimes used as a model for the Einstein equations –that is, the Maxwell equations– does not have this issue. When the Maxwell equations $`^aF_{ab}=J_b`$ are written in terms of a vector potential $`F_{ab}=_{[a}A_{b]}`$, the other half of the system $`_{[a}F_{bc]}=0`$ is identically satisfied and the time component $`A_0=\varphi `$ of the vector potential may be considered arbitrary (just as the lapse function is arbitrary insofar as the 3+1 equations are concerned). It is a straightforward exercise to put the equations in “ADM form” by defining a variable $`K_i\dot{A}_i`$, and to show that the projection $`e^b^a_{[a}A_{b]}=e^bJ_b`$ for $`e^b=\delta _x^b`$ turns out to be exactly one of the evolution equations. Exactly the same thing happens in the Lorenz gauge. The underlying reason being, of course, that there is only one constraint $`Cn^b^a_{[a}A_{b]}n^bJ_b`$ with $`n^b=\delta _0^b`$ which is preserved by the evolution by means of the trivial equation $`\dot{C}=0`$.

A separate issue not treated here is whether or not the “Einstein boundary conditions”, that is: $`G_{ab}e^b=0`$ up to linear combinations, would lead to a well-posed initial-boundary value problem for the Einstein-Christoffel formulation in the sense that the boundary conditions preserve the well-posedness of the initial value problem. In this respect, Calabrese et al Calabrese et al.  show how to obtain a subset of the constraint-preserving boundary conditions that lead to an overall well-posed initial-boundary value problem –including constraint propagation– for the linearized “generalized” Einstein-Christoffel system. The main difference between the Einstein-Christoffel system and its modified counterpart (aside from the fact that the variable $`f_{kij}`$ is defined differently) is that the constraints of the modified system propagate in a symmetric hyperbolic manner, a difference that does not affect the well-posed character, so one may expect that a subset of well-posed boundary conditions can be found in the Einstein-Christoffel case as well. Furthermore, the “Einstein boundary conditions” clearly exist and can be written down for any formulation of the Einstein equations, including those that are symmetric hyperbolic and that have symmetric hyperbolic constraint propagation, by simply expressing Eqs. (6) in terms of the fundamental variables of the problem. There are preliminary indications that, in such cases, linear combinations of $`G_{ab}e^b=0`$ exist which are well posed in the sense of Calabrese et al. .

But most of the numerical relativity effort is currently driven by formulations that are not strongly hyperbolic (see, for instance, Baumgarte and Shapiro (2003); Shibata and Nakamura (1995); Laguna and Shoemaker (2002); Alcubierre et al. (2001); New et al. (2000); Miller and Suen ; Baker et al. (2002)) in which case sensible boundary conditions are also needed, whether or not they are well-posed. For those formulations, the constraint-preserving boundary conditions of Calabrese et al may not even be defined, yet the projections $`G_{ab}e^b=0`$ evaluated on the boundary are always available, are easy to write in exact non-linear form, and provide physical boundary conditions that are clearly necessary in order to enforce the constraints, the question of well-posedness becoming irrelevant.

There is the question of why should numerical simulations use these boundary equations at all if they can do just fine by simply “pushing the boundaries out” –even though doing so increases the computational requirements considerably. In this respect, we point out the existence of the following conjecture based on the fact, observed by many authors, that the run time of numerical simulations is noticeably shortened by growing constraint violations. The conjecture, yet to be proven rigorously, is that the run time might be extended by efficiently controlling the constraint violations. Since we have shown that the boundaries are intimately connected with constraint propagation, if the conjecture turns out to be true then the “Einstein boundary conditions” may have a role to play in extending the run time of numerical simulations.

###### Acknowledgements.

We gratefully acknowledge support by the NSF under grants No. PHY-0070624 to Duquesne University and No. PHY-0135390 to Carnegie Mellon University.