text:
# Consistent deformations method applied to a topological coupling of antisymmetric gauge fields in D=3

## Abstract

In this work we use the method of consistent deformations of the master equation by Barnich and Henneaux in order to prove that an abelian topological coupling between a zero and a two form fields in D=3 has no nonabelian generalization. We conclude that a topologically massive model involving the Kalb-Ramond two-form field does not admit a nonabelian generalization. The introduction of a connection-type one form field keeps the previous result.

, and thanks: Electronic address: carlos@fisica.ufc.br

PACS: 11.15.-q, 11.10.Ef, 11.10.Kk

Keywords: method of consistent deformations; nonabelian gauge theories; antisymmetric tensor gauge fields; BRST/anti-BRST symmetry; topological mass generation

Antisymmetric tensor gauge fields coupled to abelian or nonabelian gauge fields has been studied during the last decade. In particular, four-dimensional BF models has been extensively investigated in different contexts (see Ref.  and references therein). Besides supergravity and string theories , two form fields play an important role in topological mass generation mechanism to Kalb-Ramond fields . Indeed, a nonabelian theory involving an antisymmetric tensor field coupled to a gauge field appears as an alternative mechanism for generating vector bosons masses, similar to the theory of a heavy Higgs particle . Therefore analysis of nonabelian topological terms deserves some attention. It is worth to mention a generalization to a compact nonabelian gauge group of an abelian mechanism in the context of nonabelian quantum hair on black holes .

Using the technique of consistent deformation, Henneaux et al. , have proved that is not possible to generalize the four dimensional topological mass mechanism to its nonabelian counterpart with the same field contents and fulfilling the power-counting renormalization requirements. In this way, they put in more rigorous grounds the need to add an auxiliary field.

Let us begin by presenting a recently proposed abelian three-dimensional action with a topological term involving a two-form gauge field $`B_{\mu \nu }`$ and a scalar field $`\phi `$ :

$$S_{inv}^A=d^3x\left(\frac{1}{12}H_{\mu \nu \alpha }H^{\mu \nu \alpha }+\frac{1}{2}_\mu \phi ^\mu \phi +\frac{m}{2}ϵ^{\mu \nu \alpha }B_{\mu \nu }_\alpha \phi \right),$$
(1)
where $`H_{\mu \nu \alpha }`$ is the totally antisymmetric tensor $`H_{\mu \nu \alpha }=_\mu B_{\nu \alpha }+_\alpha B_{\mu \nu }+_\nu B_{\alpha \mu }`$. The action (1), is invariant under the transformation

$$\delta \phi =0,\delta B_{\mu \nu }=_{[\mu }\omega _{\nu ]}.$$
(2)
The model described by action (1) can be consistently obtained by dimensional reduction of a four-dimensional $`BF`$ model if we discard the Chern-Simons-like terms .

The purpose of the present note is to prove that there is no power counting renormalizable nonabelian generalization of the action (1). We begin making an analysis of the possibility to construct the nonabelian action with $`\phi \phi ^a`$ and $`B_{\mu \nu }B_{\mu \nu }^a`$, i.e, with the same field content, and the same number of local symmetries, by making use of the method of consistent deformation . As will be proved, there is a no-go theorem for this construction. The same occur with an introduction of a connection-type one form gauge field.

Let us now apply the consistent deformation method described in . We shall start therefore with the following invariant action

$$S_0^{}=d^3x\left(\frac{1}{12}H_{\mu \nu \alpha }^aH^{a\mu \nu \alpha }+\frac{1}{2}_\mu \phi ^a^\mu \phi ^a\right),$$
(3)
where now $`\phi ^a`$ and $`H_{\mu \nu \alpha }^a`$ are scalar fields and the abelian curvature tensor for a set of $`N`$ fields. All fields are valued in the Lie algebra $`𝒢`$ of some Lie group G. Since we are interested if the mass term can exist in a nonabelian extension of (1), the mass parameter will be considered as a deformation parameter. The action (3) is invariant under the transformations

$$\delta \phi ^a=0,\delta B_{\mu \nu }^a=_\mu \omega _\nu ^a_\nu \omega _\mu ^a.$$
(4)
Since the transformation of $`B_{\mu \nu }^a`$ is reducible, we introduce a set of ghosts $`(\eta _\mu ^a,\rho ^a)`$, where $`\eta _\mu ^a`$ is a ghost for the gauge transformation of $`B_{\mu \nu }^a`$, and $`\rho ^a`$ the ghost for ghost for taking into account this reducibility. For all fields of the model we introduce the corresponding antifields $`(B_{\mu \nu }^a,\phi ^a,\eta _\mu ^a,\rho ^a)`$. The antifields action reads

$$S_{ant}^{}=d^3x\left(\frac{1}{2}B^{\mu \nu a}_{[\mu }\eta _{\nu ]}^a+\eta ^{\mu a}_\mu \rho ^a\right).$$
(5)
The free action

$$S_0=S_0^{}+S_{ant}^{},$$
(6)
is solution of the master equation

$$(S_0,S_0)=0,$$
(7)
with

$$(S_0,S_0)=d^3x\left(\frac{\delta S_0}{\delta \phi ^a}\frac{\delta S_0}{\delta \phi ^a}+\frac{1}{2}\frac{\delta S_0}{\delta B^{a\mu \nu }}\frac{\delta S_0}{\delta B_{\mu \nu }^a}+\frac{\delta S_0}{\delta \eta ^{a\mu }}\frac{\delta S_0}{\delta \eta _\mu ^a}+\frac{\delta S_0}{\delta \rho ^a}\frac{\delta S_0}{\delta \rho ^a}\right)$$
(8)
The nilpotent BRST transformation $`s`$ on all fields and antifields is

$$\begin{array}{cc}s\phi ^a=0,\hfill & s\phi _a^{}=^2\phi ,\hfill \\ & \\ sB_{\mu \nu }^a=_\mu \eta _\nu ^a_\nu \eta _\mu ^a,\hfill & sB_a^{\mu \nu }=_\rho H_a^{\rho \mu \nu },\hfill \\ & \\ s\eta _\mu ^a=_\mu \rho ^a,\hfill & s\eta _a^\mu =_\rho B_a^{\rho \mu },\hfill \\ & \\ s\rho ^a=0,\hfill & s\rho _a^{}=_\mu \eta _a^\mu .\hfill \\ & \end{array}$$
(9)
We show in the table (1) below, the canonical dimension and the ghost number for all fields and antifields of the model.

Having the ghost number and dimension of all fields and antifields at hand, we are now able to solve our problem using the consistent deformation method. The action (6) will be deformed to a new action $`S`$ in powers of the deformation parameters:

$$S=S_0+\underset{i}{}g_iS_i+\underset{i,j}{}g_ig_jS_{ij}+\mathrm{},$$
(10)
where $`S_i,S_{ij}..`$ are local integrated polynomials with ghost number zero and dimension bounded by three, and $`g_i`$ are the deformed parameters with nonnegative mass dimension. The action (10) must satisfy the master equation

$$(S,S)=0.$$
(11)
Expanding the master equation (11) in powers of the deformation parameters, we have

$$(S_0,S_0)=0,$$
(12)
$$(S_0,S_i)=0,$$
(13)
$$2(S_0,S_{ij})+(S_i,S_j)=0.$$
(14)
The equation (12) is the master equation for the $`S_0`$, and it does not give any additional information. The equation (13) tell us that $`S_i`$ has to be a BRST invariant under (9). We must neglect BRST exacts, since this correspond to fields redefinitions. The last equation (14) is satisfied only if the antibracket $`(S_i,S_j)`$ is a trivial cocycle.

Let us now construct all $`S_i`$ solution of equation (13). First we focus our attention to terms that do not deform the gauge symmetry, i.e, terms constructed with the fields only. Due to trivial BRST transformation of $`\phi ^a`$, the all possible terms with this field are

$$S_1=d^3x\left(\alpha _a\phi _a\right),S_2=d^3x\left(\alpha _{ab}\phi _a\phi _b\right),$$
(15)
$$S_3=d^3x\left(\alpha _{abc}\phi _a\phi _b\phi _c\right),S_4=\left(d^3x\alpha _{abcd}\phi _a\phi _b\phi _c\phi _d\right),$$
(16)
$$S_5=d^3x\left(\alpha _{abcde}\phi _a\phi _b\phi _c\phi _d\phi _e\right),S_6=d^3x\left(\alpha _{abcdef}\phi _a\phi _b\phi _c\phi _d\phi _e\phi _f\right),$$
(17)
where $`\alpha ^{}s`$ are parameters. The most general invariant local integrable terms that can be constructed with $`B_{\mu \nu }^a`$ and $`\phi ^a`$ mixed are

$$S_7=d^3x\left(m_{ab}ϵ^{\mu \nu \alpha }H_{a\mu \nu \alpha }\phi _b\right),S_8=d^3x\left(m_{abc}ϵ^{\mu \nu \alpha }H_{a\mu \nu \alpha }\phi _b\phi _c\right)$$
(18)
$$S_9=d^3x\left(m_{abcd}ϵ^{\mu \nu \alpha }H_{a\mu \nu \alpha }\phi _b\phi _c\phi _d\right),$$
(19)
with $`m_{ab}`$ having dimension of mass, $`m_{abc}`$ of dimension $`1/2`$ and $`m_{abcd}`$ as a dimensionless parameter.

Observing the table (1), it is easy to see that it is impossible to construct invariant local integrated polynomials with dimension bounded by three with the antifields. This means that the algebra of the gauge symmetry is undeformed, i.e., we do not have a nonabelian generalization of the action (1), the only possibility being with an introduction of extra fields or non-renormalizable couplings.

Let us now introduce a set of abelian vectorial gauge field in order to implement the possible nonabelian generalization of (1). We take the mass dimension of all vector fields equal to one. Therefore, those fields assume a non-dynamical character. The BRST transformations are

$$sA_\mu ^a=_\mu c^a,sc^a=0,$$
(20)
where $`c^a`$ are the ghost for the abelian transformation of $`A_\mu ^a`$. We must add to the action (6) the corresponding antifield action

$$S_{ant}^{\prime \prime }=d^3xA_{a\mu }^{}^\mu c_a.$$
(21)
The new antifields have the following BRST transformations

$$sA_\mu ^a=0,sc^a=_\mu A^{a\mu }.$$
(22)
We show in the table (2) below the ghost number and dimension for the new fields (antifields).

The all possible invariant integrated local polynomials that can be constructed with all fields and antifields are

$$S_{10}=gd^3xf_{abc}\left(\phi _a^{}\phi _bc_c^\mu \phi _a\phi _bA_{c\mu }\right),$$
(23)
$$S_{11}=\mu _{ab}d^3x\left(A_{a\mu }^{}\eta _b^\mu c_a^{}\rho _b\right),$$
(24)
$$S_{12}=hd^3xk_{abc}\left(A_a^{}A_b^\mu c_c\frac{1}{2}c_a^{}c_bc_c\right),$$
(25)
where $`g,h`$ are dimensionless parameter, $`\mu `$ is a matrix with dimension 3/2, and $`f_{abc}(k_{abc})`$ are dimensionless parameters antisymmetric in its first(last) two indices. Now we perform the calculation of the antibrackets $`(S_i,S_j)`$, with $`i,j=1,2,\mathrm{}12`$, in order to fit the second order consistency condition. As we have already seen above, these antibrackets must be a BRST exact. The antibrackets $`(S_m,S_n)`$, for $`n,m=1,2,\mathrm{},9`$ are identically zero, due to absence of antifields in $`S_n,n=1,2,\mathrm{},9`$. The antibracket $`(S_{10},S_{10})`$ is

$`(S_{10},S_{10})`$ $`=`$ $`g^2{\displaystyle d^3xf_{abc}f_{ab^{}c^{}}\left(\phi _b^{}^{}\phi _bc_cc_c^{}+\phi _{[b}_\mu \phi _{b^{}]}A_c^{}^\mu c_c\right)}`$ (26)
$`{\displaystyle \frac{g^2}{2}}s\left({\displaystyle d^3xf_{abc}f_{ab^{}c^{}}\phi _b\phi _b^{}A_{c\mu }A_c^{}^\mu }\right),`$
where, $`\phi _{[b}_\mu \phi _{b^{}]}=\phi _b_\mu \phi _b^{}\phi _b^{}_\mu \phi _b`$. The first term in (26), is not a BRST trivial and it could jeopardize the nonabelian implementation. In order to circumvent this, we must have the identification $`hk_{abc}=gf_{abc}`$, and $`f_{abc}`$ being the structure constant of a Lie group. Therefore the $`S_{10}`$ and $`S_{12}`$ are replaced by the sum

$$S_{10}^{}=gd^3xf^{abc}\left(\phi ^a\phi ^bc^c^\mu \phi ^a\phi ^bA_\mu ^c+A^aA^{b\mu }c^c\frac{1}{2}c^ac^bc^c\right).$$
(27)
It is easy to see that now $`(S_{10}^{},S_{10}^{})`$ is BRST trivial

$$(S_{10}^{},S_{10}^{})=\frac{g^2}{2}s\left(d^3xf_{abc}f_{ab^{}c^{}}\phi _b\phi _b^{}A_{c\mu }A_c^{}^\mu \right).$$
(28)
The antibrackets $`(S_{10}^{},S_n)`$, with $`n=1,2,\mathrm{},6`$, give us constraint for the parameters $`\alpha `$: $`\alpha _a=\alpha _{abc}=\alpha _{abcde}=0`$, $`\alpha _{ab}=a_1\delta _{ab}`$, $`\alpha _{abcd}=a_2\delta _{ab}\delta _{cd}`$, $`\alpha _{abcdef}=a_3\delta _{ab}\delta _{cd}\delta _{ef}`$, i.e., only the terms $`\phi ^2=\phi _a\phi _a`$, $`(\phi ^2)^2`$ and $`(\phi ^2)^3`$ are permitted. The last antibrackets reads

$$(S_{11},S_{11})=0,$$

$`(S_{10}^{},S_{11})=g{\displaystyle d^3x}`$ $`f_{abc}\mu _{ab^{}}(\rho _b^{}\phi _b^{}\phi _c+\rho _b^{}A_{b\mu }^{}A_c^\mu \rho _b^{}c_b^{}c_c`$ (29)
$`\eta _b^{}^\mu _\mu \phi _b\phi _c\eta _b^{}^\mu A_{b\mu }^{}c_c),`$
$$(S_{10}^{},S_7)=gd^3xf_{abc}m_{b^{}a}\epsilon _{\mu \nu \alpha }H_b^{}^{\mu \nu \alpha }\phi _bc_c,$$

$$(S_{10}^{},S_8)=gd^3xf_{abc}(m_{b^{}c^{}a}+m_{b^{}ac^{}})\epsilon _{\mu \nu \alpha }H_b^{}^{\mu \nu \alpha }\phi _b\phi _c^{}c_c,$$

$$(S_{10}^{},S_9)=gd^3xf_{abc}(m_{b^{}c^{}d^{}a}+m_{b^{}c^{}ad^{}}+m_{b^{}ac^{}d^{}})\epsilon _{\mu \nu \alpha }H_b^{}^{\mu \nu \alpha }\phi _b\phi _c^{}\phi _d^{}c_c.$$

The last four antibrackets are not BRST trivial, representing thus an obstruction to the deformation of the master equation. The only way to remedy this is setting $`g=0`$, or setting $`S_7=S_8=S_9=S_{11}=0`$. In the case $`g=0`$ we have lost the deformation of the abelian algebra, i.e, we have a set of abelian fields not representing a nonabelian generalization of (1). In the case in which $`S_7=0`$, we have lost the mass generation of the model. We have thus proved that there are no nonabelian generalization of the action (1), even with an addition of an auxiliary vector gauge field.

It is interesting to remark that the introduction of an one form gauge connection $`A`$ is required to go further in the nonabelian generalization of our model (1), although our original abelian action (1) does not contain this field.

We can quote some works that built some kind of nonabelianization of the model under consideration. However, the work of Oda and Yahikozawa  as well as the work of Smailagic and Spallucci , have considered the connection one form $`A`$ as a flat background field (so imposing an extra constraint in the model). On the other hand, Del Cima et al. , have studied the finiteness of a three-dimensional extension of the BF model, called BFK model, which does not have kinetic term for the antisymmetric gauge field (this term violate the nonabelian gauge invariance which is straightforwardly implemented for vectorial fields).

Motivated by possible nonabelian topological mass generation for a Kalb-Ramond field in three dimensions, we have considered deformations of a model involving a topological coupling between a second rank antisymmetric tensor field and a scalar field, using the method of consistent deformations. But an obstruction arise, leading us to a no-go theorem, namely, if we require power-counting renormalizable couplings and the same field content, the nonabelian extension for the model is forbidden.

ACKNOWLEDGMENTS

We would like to thank Dr. O. S. Ventura for helpful discussions. This work was supported in part by Fundação Cearense de Amparo à Pesquisa-FUNCAP.


format: text

text:
# On free evolution of self gravitating, spherically symmetric waves. 11footnote 1Work supported by Antorchas, CONICOR, CONICET and Se.CyT, UNC

## 1 Introduction

In the past decade, due to the results of several numerical experiments, there has been a considerable understanding of the phenomena of the collapse of a spherically symmetric, selfgravitating scalar field. In particular universal scaling properties of the final mass of the black hole were discovered and then explained. For a review and references see  The present work is not aimed to continue that line, but rather to study instability problems which are recurrent in many numerical evolutions near black holes when those evolutions are completely free, that is when the constraints are not solved for, but remain being satisfied just because the evolution equations are solved for. There are few cases in which free evolution has been successful, some of them,  use light cone coordinates and in fact some of the equations being solved could be considered as constraints. Others, , uses a symmetric hyperbolic formulation and use the freedom still available in that setting to suppress the main instability found there and so reach a stable propagation, alas first order numerical dissipation is used in that scheme. Finally there are conformal space-time methods  where no instability seems to be present, nor should be very relevant, for only short time evolutions suffices to cover the whole future domain of dependence.

In this work we use Eddington-Finkelstein coordinates and are able to choose equations and variables so that the final system of equations is remarkably simple and manifestly symmetric hyperbolic. Such simplicity allows as to study the nature of the instability and proceed to avoid it.

In the next section we introduce the equations and discuss which evolution system, among all possible equivalent ones, is suitable for free evolution. It turns out that two out of the five characteristics of the system can be specified at will by modifying the evolution equations. This is done by adding terms proportional to the constraints. These characteristics can even become complex, signaling the failure of the evolution equations to be hyperbolic, and so to have a well posed initial value formulation. The change on the characteristics of the evolution system reflects on an identical change on the characteristics of the evolution equations satisfied by the constraint equations. Thus it is possible to understand that the change in boundary conditions induced by the change on the direction of the characteristics is just the one needed to ensure the correct preservation of the constraint equations under evolution. Among all possible evolution system we take a specific one which is both the simplest in terms of computations needed and also on dealing with constraint preservation. We then look at the residual gauge freedom of the Eddington-Finkelstein coordinates, that freedom is the only dynamics which is present in the vacuum case, and is the mode that can generate instabilities. After that we discuss the initial-boundary value problem for the specific evolution system chosen and determine which fields must be given in order to have a unique evolution, in particular we discuss a case in which the characteristics change at the inner boundary according to whether that boundary is inside or outside the horizon. We then discuss the instabilities of the system, due to the simplicity of the evolution system chosen this task is trivial and we can display explicitly the instability. We can decide in which cases there will be problems and the type of them. In particular near the stationary regime, that is when most of the scalar field has dissipated away or fallen into the black hole the system can be made stable using simple initial-boundary values for the gauge fields.

Having completed the analytic study of the equations in section 3 we discuss the numerical methods employed, look at the convergence tests and constraint preservation under evolution for several types of initial-boundary data sets.

In the next section we comment in several results obtained by the simulations, they do not pretend to be novel or interesting, but are done just to confirm the good properties of the evolution system chosen. They include ringing and tail studies, the relation mass-flux for boundary data and the relation total mass black hole mass for different types of collapse.

## 2 The equations

We consider the Eddington-Finkelstein coordinate system (), the time dependent spherically symmetric metric is

$$ds^2=(\alpha ^2+a^2\beta ^2)dt^2+2a^2\beta dtdr+a^2dr^2+r^2b^2d\mathrm{\Omega }^2$$

with $`a,b,\alpha `$ and $`\beta `$ being functions of $`r`$ and $`t`$, and $`d\mathrm{\Omega }^2`$ the metric of the unit-sphere. The areal locking condition ($`\dot{b}=0`$) and the condition that the vector $`_t_r`$ be null is achieved choosing the lapse and the shift as

$$\beta =\frac{raK^\theta _\theta }{1+raK^\theta _\theta }\alpha =\frac{a}{1+raK^\theta _\theta }.$$

where $`K^\theta _\theta `$ is the corresponding coordinate component of the extrinsic curvature of the constant time hypersurfaces. Moreover this choice corresponds to the choice $`b=1`$, and with it the metric becomes

$$ds^2=a^2((2\beta 1)dt+dr)(dt+dr)+r^2d\mathrm{\Omega }^2.$$

The full set of evolutions and constraint equations for the geometric variables $`(a,K^\theta {}_{\theta }{}^{},K^r{}_{r}{}^{})`$ are 222 Notice some corrections with respect to the equations in , also notice the difference on the choice of gravitational constant.

$`_ta`$ $`=`$ $`_r(a\beta )a^2(1\beta )K^r_r`$
$`_tK^\theta _\theta `$ $`=`$ $`\beta _rK^\theta {}_{\theta }{}^{}+{\displaystyle \frac{_r\beta }{ar}}+a(1\beta )K^\theta {}_{\theta }{}^{}(K^r{}_{r}{}^{}+2K^\theta {}_{\theta }{}^{})+{\displaystyle \frac{1\beta }{r^2}}(a{\displaystyle \frac{1}{a}})`$
$`_tK^r_r`$ $`=`$ $`\beta _rK^r{}_{r}{}^{}+{\displaystyle \frac{_r^2\beta }{a}}+{\displaystyle \frac{\beta 1}{a}}[{\displaystyle \frac{_r^2a}{a}}\left({\displaystyle \frac{_ra}{a}}\right)^2{\displaystyle \frac{2_ra}{ra}}]+{\displaystyle \frac{_r\beta _ra}{a^2}}`$
$`+a(1\beta )K^r{}_{r}{}^{}(K^r{}_{r}{}^{}+2K^\theta {}_{\theta }{}^{})+{\displaystyle \frac{8}{a}}\mathrm{\Phi }^2(\beta 1).`$
and

$`C_a`$ $`=`$ $`_ra+{\displaystyle \frac{1}{2r}}(a^3a)+{\displaystyle \frac{a^3r}{2}}K^\theta {}_{\theta }{}^{}(2K^r{}_{r}{}^{}+K^\theta {}_{\theta }{}^{})2ra(\mathrm{\Phi }^2+\mathrm{\Pi }^2)=0`$
$`C_{K^\theta _\theta }`$ $`=`$ $`_rK^\theta {}_{\theta }{}^{}+{\displaystyle \frac{K^\theta {}_{\theta }{}^{}K^r_r}{r}}{\displaystyle \frac{4}{a}}\mathrm{\Phi }\mathrm{\Pi }=0`$ (2)
respectively. While the wave equation for the massless scalar field becomes the first-order (in time) system

$`_t\mathrm{\Phi }`$ $`=`$ $`_r\left(\beta \mathrm{\Phi }+(1\beta )\mathrm{\Pi }\right)`$
$`_t\mathrm{\Pi }`$ $`=`$ $`{\displaystyle \frac{1}{r^2}}_r(r^2(\beta \mathrm{\Pi }+(1\beta )\mathrm{\Phi })`$ (3)
where $`\mathrm{\Pi }=_t\varphi `$, and $`\mathrm{\Phi }=_r\varphi `$.

In order to obtain first order evolution equations, we perform the change of variables $`(f,g,h,\chi _+,\chi _{})`$

$$f=K^\theta {}_{\theta }{}^{}ar+1g=\frac{a^2}{f}h=\frac{1}{af}(\frac{a}{2}_rf+_raa^2K^r{}_{r}{}^{}).$$

$$\chi _+=\frac{(\mathrm{\Pi }+\mathrm{\Phi })r}{2}\chi _{}=\frac{(\mathrm{\Pi }\mathrm{\Phi })r}{2}$$

with this change the shift becomes $`\beta ={\displaystyle \frac{f1}{f}}`$ and

$$ds^2=g\left((f2)dt+fdr\right)(dr+dt)+r^2d\mathrm{\Omega }^2.$$

Then the full set of evolution and constraint equations for the metric variables $`(f,g,h)`$ are

$`_tf`$ $`=`$ $`{\displaystyle \frac{f}{r}}+{\displaystyle \frac{g2+f}{r}}`$
$`_tg`$ $`=`$ $`{\displaystyle \frac{g}{f}}{\displaystyle \frac{f}{r}}+{\displaystyle \frac{f2}{f}}{\displaystyle \frac{g}{r}}+2gh{\displaystyle \frac{g(g2+f)}{rf}}`$ (4)
$`_th`$ $`=`$ $`{\displaystyle \frac{h}{r}}{\displaystyle \frac{1}{2fr}}{\displaystyle \frac{f}{r}}+{\displaystyle \frac{g2}{2fhr}}{\displaystyle \frac{g}{r}}{\displaystyle \frac{(gf)h}{fr}}{\displaystyle \frac{(g2+f)}{2fr^2}}+2({\displaystyle \frac{\chi _++\chi _{}}{fr}})^2`$
and

$`C_f`$ $`=`$ $`{\displaystyle \frac{f}{r}}2(f2)h+{\displaystyle \frac{f+g2}{r}}4{\displaystyle \frac{\chi _+^2}{fr}}=0`$
$`C_g`$ $`=`$ $`{\displaystyle \frac{g}{r}}{\displaystyle \frac{2g(hfr\chi _+^2+\chi _{}^2)}{fr}}=0`$ (5)
while the wave equation for the massless scalar field becomes

$`_t\chi _+`$ $`=`$ $`{\displaystyle \frac{\chi _+}{r}}{\displaystyle \frac{(f2)\chi _{}}{fr}}`$
$`_t\chi _{}`$ $`=`$ $`{\displaystyle \frac{f2}{f}}{\displaystyle \frac{\chi _{}}{r}}+{\displaystyle \frac{2\chi _{}}{f^2r}}(2r(f2)h+2gf)+{\displaystyle \frac{8(\chi _{})\chi _+^2}{f^3r}}{\displaystyle \frac{\chi _+}{r}}`$ (6)
where we have already used one of the constraint equations.

Clearly the evolution equations for the scalar field are symmetric hyperbolic as a subsystem and their characteristics are fixed, they correspond to the two null directions. On the other hand the set of evolution equations for the metric coefficients is not unique, for we can add to them terms proportional to the constraint equations and change some of their characteristics at will. To study this possibilities we add to the evolution equations for $`f`$, $`g`$, and $`h`$ respectively the following terms $`(K_{ff}1)C_f+K_{fg}C_g`$, $`(K_{gf}+\frac{g}{f})C_f+(K_{gg}\frac{f2}{f})C_g`$, $`(K_{hf}\frac{g2}{2frh})C_f+(K_{hg}+\frac{1}{2fr})C_g`$. We get a system with the following principal part matrix, 333Recall that the principal part matrix is the one form by the coefficients of the derivative terms appearing on the equations.

$$\left(\begin{array}{ccc}K_{ff}& K_{fg}& 0\\ K_{gf}& K_{gg}& 0\\ K_{hf}& K_{hg}& 1\end{array}\right)$$

The characteristics of this system are given by the one forms $`(1,1)`$, $`(1,\lambda _+)`$, and $`(1,\lambda _{})`$, where

$$\lambda _\pm =(K_{ff}+K_{gg}\pm \sqrt{(K_{ff}K_{gg})^2+4K_{fg}K_{gf}})/2.$$

The first one corresponds to the propagation of $`h`$ and it does so along the null incoming direction. It is the only fix characteristic. The other two can take any values, for instance if we take $`K_{fg}=K_{gf}=0`$ then the characteristics are $`(1,K_{ff})`$, and $`(1,K_{gg})`$, thus choosing these two values we can prescribe any propagation direction we please, they can be incoming, or outgoing time-like, null or even space-like. Even more, choosing $`K_{ff}=K_{gg}=0`$, $`K_{fg}=K_{gf}0`$ the characteristics become imaginary and so the system is not even hyperbolic! So we see that the different ways of writing Einstein’s equations are only equivalent in the sense that solutions of one are solutions of the others, but some of their evolution equations do not even have a well posed initial value problem. We have tried to evolve numerically a system as the last one (with $`K_{fg}=K_{gf}=1`$) and as expected it is completely unstable, giving an explosive blow up of higher wave number modes. To avoid the fact that no good boundary conditions are known for non-symmetric hyperbolic systems these equations were tested using periodic boundary conditions and the initial data did not satisfy the constraint equations. But the instability does not seems to be sensitive to the failure of the constraint equations to be fulfilled nor to the particular boundary values chosen for those instabilities appear locally at every point of the evolution hypersurface.

Thus, to start with, we must choose coefficients such that the system is hyperbolic. Since the coefficients $`K_{hf}`$, $`K_{hg}`$ do not alter the characteristics we shall take them to vanish, thus simplifying a bit the numerical computations. The sign of the characteristics at the boundary points determine whether or not values for the different fields must be given or not as boundary conditions. But not all these values are free, for they must be consistent with constraint propagation. In order for the constraints to remain enforced along evolution we must make sure that the resulting boundary conditions for the evolution equations obeyed for the constraints ensure that they have a unique solution, and it is the trivial one. It follows trivially from the form of the constraints that the characteristic matrix for the constraints evolution equations is:

$$\left(\begin{array}{ccc}K_{ff}& K_{fg}& \\ K_{gf}& K_{gg}& \end{array}\right)$$

So that the propagation directions for them coincide with the propagation directions of the pair $`(f,g)`$. This fact seems to imply that the freedom in choosing part of the pair $`(f,g)`$ at the boundaries is just the one needed for getting maximally dissipative boundary conditions for the constraint system, thus making sure they propagate correctly. Unfortunately even for this simple case a rigorous proof of this has proven be difficult 444The main difficulty arises when the $`K`$ coefficients depend on the metric components, and the resulting conditions are complicated to implement numerically. Thus it seems that the best one can do retaining simplicity is to choose coefficients so that the characteristics for the pair $`(f,g)`$ are incoming or tangent in both boundaries. This way no boundary condition is needed for the pair $`(f,g)`$ and therefore no boundary condition must be satisfied by the constraints, thus, there is no danger that the constraint would cease to be satisfied during evolution. The simplest case of all is when both are tangent, and this is the case for instance if we take all coefficients to vanish. With this choice the evolution system reduces to,

$`_tf`$ $`=`$ $`2(f2)h+{\displaystyle \frac{4\chi _+^2}{rf}}`$
$`_tg`$ $`=`$ $`2gh{\displaystyle \frac{2g\chi _+^2}{fr}}+{\displaystyle \frac{2g(f2)\chi _{}^2}{rf^2}}`$ (7)
$`_th`$ $`=`$ $`{\displaystyle \frac{h}{r}}{\displaystyle \frac{(g2)\chi _+^2}{f^2r^2}}+{\displaystyle \frac{4\chi _+\chi _{}}{f^2r^2}}+{\displaystyle \frac{g\chi _{}^2}{f^2r^2}}`$
And similarly the evolution equations for the constraints are just ordinary differential equations,

$`_tC_f`$ $`=`$ $`{\displaystyle \frac{2}{rf^2}}(hf^2r2\chi _+^2)C_f`$ (8)
$`_tC_g`$ $`=`$ $`{\displaystyle \frac{2g}{rf^3}}(f\chi _+^2f\chi _{}^2+\chi _{}^2)C_g+{\displaystyle \frac{2}{rf^2}}(hrf^2f\chi _+^22\chi _{}^2+f\chi _{}^2)C_g`$ (9)
Since they are homogeneous ordinary differential equations the constraints will remain to be satisfied along evolution if they are initially satisfied. No extra condition is need at (nor can be imposed on) the boundaries. This is the system we shall study in detail from now on. Note that this system not only is symmetric hyperbolic, but it is also diagonal. In fact we reach such a simple system when looking for variables which would give diagonal systems. In figure 1 we show all characteristics of the chosen system.

### 2.1 The Gauge Freedom

The chosen lapse and shift functions depend on the dynamical variables and therefore do not fix completely the gauge. To display this freedom we introduce the null coordinate $`u=(t+r)`$ 555We thanks M. Tiglio for pointing to us this way to proceed.. In the $`(u,r)`$ coordinates the metric becomes

$$ds^2=g((f2)du+2dr)du+r^2d\mathrm{\Omega }^2,$$

For Schwarzschild the choice $`g=1`$, $`f2=(12M/r)`$ gives the standard form of the metric. If we choose some other coordinate $`\stackrel{~}{u}=\stackrel{~}{u}(u)`$ the functional form remains the same if $`\stackrel{~}{g}=g/\frac{d\stackrel{~}{u}}{du}`$ and $`\frac{\stackrel{~}{f}2}{\stackrel{~}{}g}=\frac{f2}{g}`$.

The invariant ratio $`\frac{f2}{g}`$ is related to the mass function, indeed,

$$m=\frac{r}{2}(1+\frac{f2}{g})$$

is the usual form of the local mass for spherical symmetric space-times.

Another invariant is the surface where $`f`$ takes the value $`2`$, this can be seen to be an apparent horizon.

From the first constraint equation, (5), it follows that under this gauge transformation, in the vacuum case,

$$\stackrel{~}{h}=h+\frac{1}{2}\frac{\frac{d\stackrel{~}{u}}{du}}{\frac{d^2\stackrel{~}{u}}{du^2}}$$

To fix completely the gauge we must uniquely define the coordinate $`u`$, this can be done by choosing their values at the initial time slice and at the outer boundary of the integration region. We shall do this by choosing the values for $`h`$ at those points.

### 2.2 Initial Boundary Value Problem

We are interested in performing a time evolution between two fixed (areal) radiae, $`R_{in}`$, $`R_{out}`$, starting at a given initial surface, $`t=t_0`$. Thus we have to see what initial-boundary values can be prescribed.

With the standard choice of equations (6,7) it is clear that we need to give initial values for the fields $`(f,g,h,\chi _+,\chi _{})`$. Since they must satisfy the constraint not all initial data can be freely specified. We have taken a set of initial data sets for which the constraints are easily solved. We give an arbitrary initial value for $`\chi _+`$ and solve

$$\frac{m}{r}=\chi _{+}^{}{}_{}{}^{2},m|_{r=r_{in}}=m_0,$$

Then the following choice is a solution of both constraint equations,

$`\chi _{}`$ $`:=`$ $`0,`$
$`g`$ $`:=`$ $`1,`$ (10)
$`f`$ $`:=`$ $`1+2m/r`$
$`h`$ $`:=`$ $`{\displaystyle \frac{\chi _{+}^{}{}_{}{}^{2}}{rf}}`$
In our numerical simulations we have taken $`\chi _+`$ to be a Gaussian or a power of r times an exponential decay, so that the mass integration could be done exactly. The integration constant for $`m`$, the mass at $`R_{in}`$, can take any value, in general we have taken either 666 Due to the scaling properties of the system there is no loss of generality in choosing only just these two values. $`m|_{R_{in}}=1`$ or $`m|_{R_{in}}=0`$

We now discuss the boundary conditions (see figure (1). If the inner boundary is inside the apparent horizon (which is located at $`f=2`$) there is no incoming characteristic at it and so no boundary condition is needed (nor can be prescribed) there. In fact in this case this boundary is actually space-like. If the inner boundary is outside the horizon there is one characteristic incoming into the integration region, namely the one for $`\chi _{}`$, so we have to prescribe some value for it. To allow for an apparent horizon marching across the boundary the numerical code checks for the value of $`f`$ at the boundary and, depending whether its value is smaller or bigger than 2, imposes a boundary condition for $`\chi _{}`$. When needed we have normally used a null incoming radiation condition, ($`\chi _{}|_{R_{in}}=0`$) but it is clear than in that case there is no natural boundary condition which can mimic the physical collapse one is trying to simulate.

At the outer boundary, and provided the horizon is towards the inside of that boundary ($`f|_{R_{outer}}<2`$), there are two incoming characteristics, one corresponding to $`\chi _+`$ the other to $`h`$. We can prescribe arbitrary values for them there, and we have performed several runs prescribing either incoming scalar field radiation or pure gauge modes.

With these initial-boundary value conditions the problem is well posed and so, given smooth data we obtain a smooth solution valid for a finite time interval.

### 2.3 Analytic Instabilities

The evolution system is well posed, so there are no instabilities of the explosive type, that is those whose rate grows unbounded with the frequency, rather the expected instabilities have an exponent with bounded positive real part. These instabilities essentially affect the pair $`(f,g)`$ which controls the characteristics of the system, if they are present, and act for a long enough time they can ruin the hyperbolicity of the system –by making some of the propagation speeds to diverge– and so its well posedness. But much before that happens one of the characteristic speeds grows so much that the system becomes numerically unstable due to violations of the Courant-Friedrich stability condition.

For simplicity we analyze the instability for the vacuum case. In that case the equations are:

$`_tf`$ $`=`$ $`2(f2)h`$
$`_tg`$ $`=`$ $`2gh`$ (12)
$`_th`$ $`=`$ $`{\displaystyle \frac{h}{r}}`$
Given initial data for $`(f,g,h)`$ and boundary data at the outer boundary for $`h`$, we can first solve for $`h`$, $`h(t,r)=h_b(t+rR_{out})`$ and then for $`(f,g)`$. We have,

$`f(t,r)2`$ $`=`$ $`(f_0(r)2)e^{2_0^th_b(\stackrel{~}{t}+rR_{out})𝑑\stackrel{~}{t}}`$
$`g(t,r)`$ $`=`$ $`g_0(r)e^{2_0^th_b(\stackrel{~}{t}+rR_{out})𝑑\stackrel{~}{t}}`$
Thus we see that $`g`$ remains always positive, and the horizon –the radius at which $`f(t,r)=2`$ remains fixed. But the equal time surfaces tilt, and so if $`h`$ takes big positive values, $`f`$, and $`g`$ grow exponentially and we have an analytical instability, which in numerical simulations would result in a numerical one if no appropriate methods are used. If a scalar field is present, no matter how small, the tilt of the equal time surfaces would increase the propagation speed of the outgoing modes, $`v=\frac{f2}{f}`$. At points inside the horizon ($`f_02>0`$) and so there $`f2`$ grows exponentially. The propagation speed $`v`$ goes to $`1`$. At points outside the horizon ($`f_02<0`$) and so there $`f`$ eventually becomes negative. Notice that when $`f`$ vanish the equal time hypersurface becomes null and so the time evolution of all outgoing modes (in this case only the one from the scalar field) is ill posed in this gauge.

On the other hand if $`h`$ becomes large but negative, $`f2`$ and now the hypersurface of constant $`r`$ approach null surfaces. In that case the evolution still can be done, but the loss of accuracy is very important, even for moderate values of $`h`$.

We need therefore to keep $`h`$ moderately small during evolution and propagating inwards without any important residual part due to the numerical method. This can be achieved in many cases by just prescribing the initial-boundary conditions for $`h`$ so that it starts small. There are nevertheless situations for which $`h`$ becomes big and the simulation runs into the instabilities, this is so, for instance, if one takes boundary data for the incoming component of the scalar field with a very big amplitude (final masses of order $`10`$). In that case the $`h`$ generated along the evolution of the scalar field is important and the systems becomes unstable. On the other hand, for situations describing the final stages of collapse, namely after most of the scalar field has fallen into the black hole and the geometry has settled in a near stationary regime the instabilities do not show up. We have tested that with evolutions lasting for several thousand masses.

## 3 Numerical Methods

To perform the simulations we use a standard fourth order method of lines. The space derivatives discretization is done with standard fourth order (we have also used second order discretizations with similar results) discretized derivatives which are one sided at the boundaries. We have used symmetric and non-symmetric, but full fourth order, ones. The symmetric ones, , are guaranteed to be stable, but are only second or third order accurate at the last points of the grid near the boundary. This accounts, for instance, for spikes at the error on the constraint propagation check, some of these errors sometimes propagate inwards, but stay bounded and diminish as the grid spacing is decreased. The non-symmetric discretizations used were full fourth order (but with bigger errors near the boundary) and they give smaller errors on constraints evolution near the boundaries. But there is no guaranty that the scheme would be stable. Most of the calculations, including all the results and test given in this paper, were done with the symmetric ones. Their symmetry implies they are norm preserving and so they do not have any dissipation. Fourth order centered differences approximations introduce hight frequency modes with backwards propagations speeds (group velocities) which are 5/3 bigger than the exact speeds. These modes become important when analyzing the tail behavior of the scalar field, for there the values of the scalar field become so small that the contribution from them becomes important. The higher than light propagation of them implies that boundary effects take place before they should. These high frequency modes become smaller as the grid size is decreased, thus tail decay studies are possible if enough grid points are included, but this should be very costly in higher dimensions. For these studies it seems better to use centered difference schemes of second order, for they have a maximum group velocity of one.

Time discretization is done using the standard fourth order Runge-Kutta method.

We have done three types of simulations: The first type of simulation has no-incoming boundary conditions ($`\chi _+=h=0`$ at the outer boundary) and an initial data satisfying (10). One starts with a solution representing an inner black hole (usually of unit mass) and certain amount of infalling scalar field. The evolution should result in a bigger black hole and some amount of scalar field going outwards. Thus a ringing and latter a power law tail should be observed in the outgoing component of the scalar field, $`\chi _{}`$.

The second type are vacuum solutions representing pure gauge dynamics. They have Schwarzschild or Minkowski initial data and the dynamics is created by allowing incoming gauge ($`h|_{R_{out}}=h^b(t)`$) to enter the outer boundary. The evolution should have some dynamics and latter on, if the incoming gauge has a finite duration, should relax into a static black hole in a different gauge, through the simulation the mass should stay constant, as should the horizon position.

The third type of data starts with initial data corresponding to a Schwarzschild black hole of unit mass and latter on the evolution an incoming scalar field mode is injected into the system from the outer boundary during a finite amount of time ($`\chi _+|_{R_{out}}=\chi _+^b(t)`$. The evolution should create a bigger black hole and some outgoing radiation should be seen.

We shall discuss the results for these three type of simulations in the next section.

### 3.1 Convergence Test

We have perform convergence tests for the three types of simulations described above. In general it is difficult to analyst the convergence order in detail due to the fact that the space discrete derivative we use is not of the same order near the boundary (third or second order) rather than fourth. This gives rise to Q values different than the corresponding to fourth order (Q=16). Typically after boundary effects become important (that is when the fields reach the boundary) we get values over 12, but some times they drop to about 4 when the bulk of the solution is near the boundary. In figure: 3, and 3 we give values for the Q factors obtained using the $`L^2`$ and $`L^{\mathrm{}}`$ norms respectively, for a run with homogeneous boundary conditions and a Gaussian peak of the ingoing component of the scalar field midway in the integration region. This run has an initial black hole mass $`M_{Initial}=1`$, a total mass, $`M_{Total}=8.90`$ and a final black hole mass of $`M_{Final}=8.67`$ so it is well into the nonlinear regime. The peak observed at the end of the integration period is due to boundary effects, after that the values of the different Q’s go to values above 12, but never return to 16.

### 3.2 Constraint Evolution

The constraint equations are not solved for in any step of our simulations, and although they are satisfied along evolution for the exact equations here we have to check this happens at the discrete level. The figures (99) show the $`L^{\mathrm{}}`$ and $`L^2`$ norm of the constraint expressions along evolution of different type of data set used. For type I data we had: $`M_{Initial}=1`$, $`M_{Total}=1.15`$, $`M_{Final}=1.12`$. For the type II data we introduced at the outer boundary $`10`$ cycles of gauge mode of frequency $`2\pi `$, starting at $`t=10`$. The black hole mass was of course constant and had the value $`1`$. For the type III data we introduced at the outer boundary $`10`$ cycles of the $`\chi _+`$ mode with frequency $`2\pi `$, starting at $`t=10`$. We had: $`M_{Initial}=1`$, $`M_{Total}=2.6`$, $`M_{Final}=2.4`$.

It is clear that the constraint quantities remain bounded along evolution and that they go to zero as the grid size diminishes. Most of the contribution to the $`L^{\mathrm{}}`$ norm comes from the boundary, where the derivative operator used is only second or third order. These peaks near de boundary can be reduced to the same level than the rest if one uses a discrete derivative operators which are fourth order at all grid points.

If one looks at the ratios $`C_f/\frac{df}{dr}`$, $`C_g/\frac{dg}{dr}`$ which in some sense measure better the failure of the constraint to be satisfied one sees that those ratios do not change so much along evolution. For $`C_f`$ the change is of only about 7%, and are of the order of $`10^4`$. For $`C_g`$, since $`g`$ is chosen initially to be zero, the change at the beginning is rather big, but after a transient it gets to a plateau of the order of $`10^4`$.

## 4 Results

### 4.1 Decay of the scalar field, ringing and power law decay

We reproduce standard results on ringing and tail (power law decay) for the first type of simulation. In figure (11) it is plotted the log of the absolute value of $`\chi _{}`$ at $`r=25.5`$ with respect to time. It is a run where the outer boundary is at $`R_{outer}=6000`$, the space grid spacing was $`h/M=0.0234`$, where we have used the initial mass to make the ratio. footnoteIn most of the runs $`k/h2`$, this does not violate the Courant-Friedrichs-Lewi condition because we use Runge-Kutta for time evolution. The initial scalar field is given by a Gaussian $`\chi _+`$ pulse centered at $`r=10`$, giving: Initial black hole mass $`M_{Initial}=1`$, total mass $`M_{Total}=2.56`$, final black hole mass $`M_{Final}=2.45`$. At the beginning the characteristic ringing is seen and then a power law decay of the type $`t^\alpha `$ (this can be seen best plotting the local decay power $`t{\displaystyle \frac{d\mathrm{ln}(\chi _{})}{dt}}`$ as defined in , figure (11)). The power law decay computed $`t^{4.01}`$ (at $`t=7200`$), agrees quite well from the one expected for linear perturbation theory, namely $`t^4`$. See . In the plot we also show a lower resolution ($`h/M=0.0468`$) computation where it is clear that an instability starts to appear and become important at about half the time. This instability diminishes when the grid resolution is incremented, but still seems to be present at longer times. At the resolution used a feature departing form the power law decay is seeing after $`t=7200`$, this corresponds to an higher frequency mode propagating backwards at a speed $`1.66`$ the speed of light. These high frequency modes are common on all centered fourth order methods. Their amplitude decreases with increased precision. The time of arrival of that feature implies it originates near the inner boundary probably due to imperfect ingoing boundary conditions. After that feature, we again see (not shown in the figure) a even bigger feature, this one propagates at the correct speed of the problem, but its presence at that time implies that it is generated by the one traveling at higher speed when it enters the region near the outer boundary.

The possibility of giving boundary conditions which would automatically satisfy the constraint equations allow as to perform several interesting runs. In particular we made a couple of runs of type III with

$$\chi _+|_{R_{out}}=A_+(1\mathrm{cos}(\omega _+t))t[10,20]$$

with $`A_+=1/2\pi `$, and in one case $`\omega =2\pi `$, in the other $`\omega =2\pi /10`$. For these two boundary data we have,

$$_{10}^{20}\chi _+^2𝑑t=15A_+^2$$

Since the energy flux is given by

$$\frac{dm}{dt}=\frac{1}{g}(\chi _+^2\frac{(f2)^2}{f^2}\chi _{}^2),$$

and $`\frac{dg}{dt}\frac{2g}{fr}\chi _+^2`$, $`\frac{df}{dt}\frac{4}{fr}\chi _+^2`$, $`\frac{d\chi _{}}{dt}\frac{1}{r}\chi _+`$. We see that to a 1% accuracy the flux of both solutions should coincide and so the black hole masses should be similar. Even more, the ringing of both solutions should be very similar, for it mostly depends on the geometry. In figure (12) we see the mass of one of the simulations and the mass difference between both simulations multiplied by 5. We see that most of the difference is bounded to the region where the wave of scalar field is moving, leaving behind the same geometry. In figure (13) we see one of the ringings (the value of $`\chi _{}`$ at $`r=50`$) and the difference between the ringing of both solutions augmented by a factor 10. Again we see that most of the difference is at the moment where the two different wave fronts pass the point and then on the precise location of the maxima of the ringing.

### 4.2 Relation between the initial mass and the black hole mass

We have tested the relation between the mass gap given by the total initial mass minus the final black hole mass as a function of the total initial mass for type I and type III initial-boundary data sets. In figures (15, 15) we display the results. In both cases we start with a black hole of mass unit, and put our inner boundary inside its horizon.

## 5 Conclusions

It is difficult to see whether some of the ideas presented in this work can be extended to the full three dimensional case, where many more variables are present and where there is no preferred center in which to anchor a gauge as the one used here. But the model is so simple that perhaps the observations made in solving it can shed some light into the more difficult problem. Among the observations we have the following: a) If the gauge prescription does not fix completely the gauge it is expected there would be gauge modes propagating. In this case there was just one gauge freedom left (the value of $`h`$ at the initial surface and at the outer boundary), nevertheless there were three equations which acquired a non trivial propagation. b) Gauge modes, by being non-physical, can have singular behavior although the geometry given by the solution can be regular. It is clear in the model studied that the different components of the metric grow exponentially while the solution is still Schwarzschild. Thus, it seems to be necessary to isolate and keep under control the potentially unstable gauge modes. In our case this was done by choosing variables for which it was clear where the instability was, and so, choosing conveniently the initial-boundary data we could keep it small for many cases, although not always. c) The choice of gauge modes can also complicate the boundary value problem. Ideally some of the fields should acquire values at the boundary so that the constraint propagate correctly. In the model it can be seen that if part of the data for the pair $`(f,g)`$ must be given at the boundary (according to the values chosen for the coefficients $`(K_{ff},K_{fg},K_{gg},K_{gf})`$) that part must satisfy some evolution equation along the boundary, so actually there is no freedom left. But other gauge quantities, like $`h`$, can take any value. Thus, if the variables are not chosen appropriately it could become very cumbersome to find the correct boundary conditions which would, at the same time keep the constraint equations being satisfied and the gauge instabilities under control. Note that the procedure we used, getting some equations which are intrinsic to the boundary for some of the variables, is similar to the one used in , the only case where well posedness of the initial-boundary value problem in full general relativity has been asserted.

## 6 Acknowledgments

It is a pleasure to acknowledge several and important discussions with S. Frittelli, H.O. Kreiss, P. Laguna, and M. Tiglio.