format: text

text:
# A cure for unstable numerical evolutions of single black holes: adjusting the standard ADM equations

## I Introduction

In numerical relativity, Einstein’s equations are usually solved as an initial value problem. That is, the spacetime is foliated with spacelike hypersurfaces. These hypersurfaces, or slices, are characterized by their intrinsic geometry (spatial metric $`g_{ij}`$) and extrinsic curvature $`K_{ij}`$. Subsequent slices in a given foliation are connected via the lapse function $`\alpha `$ and shift vector $`\beta ^i`$ . Under this framework, the components of Einstein’s equation naturally separate into constraint and evolution equations for the dynamical variables $`g_{ij}`$ and $`K_{ij}`$. Thus, a typical procedure to construct a spacetime consists of first specifying Cauchy data $`(g_{ij},K_{ij})`$ that satisfy the constraints on the initial slice and then applying the evolution equations to update these data into the next slice. This procedure is known as free or unconstrained evolution.

In order to propagate $`(g_{ij},K_{ij})`$ using the evolution equations, one must provide in addition a prescription for choosing the lapse function and shift vector. One generally regards the lapse-shift prescription as a choice of coordinates. At the continuum level, given a lapse-shift choice and a specific set of Cauchy data, the evolution equations yield not only a unique spacetime metric expressed in a given coordinate system but also evolved data $`(g_{ij},K_{ij})`$ that continues to satisfy the constraints.

Numerical approximations are likely to complicate the picture described above. For instance, because of truncation errors, the data $`(g_{ij},K_{ij})`$ do not satisfy the continuum Einstein constraints but rather their discrete approximations. Even if consistency between the numerical approximations of the evolution and constraint equations is achieved, the evolved data $`(g_{ij},K_{ij})`$ will at best satisfy the constraints up to truncation errors. Of course, this is what one should expect. However, the presence of numerical errors could also trigger fast growing modes that render the numerical evolution unstable. Characterizing and controlling these growing modes has been and continues to be one of the most difficult and demanding tasks in numerical relativity.

The origin of these destabilizing modes is multiple. Unstable modes could be purely numerical (e.g., the Courant instability) or they could be present already at the continuum level. An example of continuum instabilities are the so called constraint-violating modes . These are solutions of the Einstein evolution equations that do not satisfy the Hamiltonian and momentum constraints. The general perception is that this class of solutions consists mostly of rapidly growing solutions, although there is no formal proof that this is indeed the case. In summary, even if one were to start an evolution with “perfect” initial data that satisfies the constraints, truncation errors could cause the numerical solution to drift into a constraint-violating, and perhaps unbound, solution.

Various groups have accumulated numerical evidence supporting as a possible cause of instabilities the particular 3+1 form used to recast Einstein’s equation as an initial value problem. These observations have partially motivated the development of hyperbolic formulations of Einstein equations . Hyperbolic formulations have the advantage that mathematical tools exist to prove existence and uniqueness of the solutions as well as well-posedness, thus providing important information for the implementation of stable discretization algorithms. If excision is used to handle the singularities, numerical simulations of black hole spacetimes require the imposition not only of outgoing, radiative boundary conditions far away from the holes but also of conditions on the excision boundary that respect the propagation of physical modes to flow down into the hole. Because hyperbolic formulations also yield information regarding the characteristics of all the field variables in the system, in principle this knowledge could be extremely useful when applying these boundary conditions.

Unfortunately, so far, none of the numerical efforts based on hyperbolic formulations of Einstein equations have been able to demonstrate their clear superiority. Currently, the formulation originally developed by Shibata and Nakamura , and later re-introduced by Baugarte and Shapiro , (BSSN) seems to be the least prone to instabilities. Interestingly enough, the BSSN formulation is not explicitly hyperbolic. Just recently, a numerical implementation of the Einstein-Christoffel hyperbolic system  has been applied to single black hole spacetime and produced three-dimensional evolutions with stability properties comparable to those using the BSSN formulation . The common lore these days is, however, that the standard Arnowitt-Deser-Misner (ADM)  formulation is the one which most easily suffers from instabilities . On the other hand, the ADM formulation has the advantage of containing a minimal number of equations, making it attractive for numerical studies.

The main objective of this paper is to show that, given a choice of lapse function and shift vector, it is possible to obtain long-term stable and convergent simulations using the standard ADM formulation if “appropriate” terms involving the constraints are added to the evolution equations. Of course, here the key ingredient is to understand what constitutes an appropriate choice of constraint terms. We show that at least in spherical symmetry, it is possible to obtain a definite prescription for adding these constraint-terms that only depends on the choice of gauge or coordinate conditions, namely the lapse and shift. The idea of adding constraint-terms to the evolution equations (adjusted ADM formulations) is not new . Just recently, Yoneda and Shinkai  studied in full detail the propagation of the constraints in the family of adjusted ADM systems.

In order to help identify what constitutes a suitable choice of constraint-terms, we introduce a definition of quasi well-posedness. The idea of quasi well-posedness is simply to make choices of lapse and shift such that, when the ADM system is enlarged with the addition of new variables (typically, but not always, quantities related to the first spatial derivatives of the three metric), one obtains a strongly hyperbolic system. Once a strongly hyperbolic system is obtained, one can make use of the standard results that apply to those systems regarding the existence, uniqueness and well-posedness of the solutions. The essence of quasi well-posedness is then that the properties of the strongly hyperbolic system are to some degree inherited by the original ADM system. The advantage of this procedure is that by looking at the strongly hyperbolic system, one is also able to gain insight about the propagation of the constraints and characteristics of the original system, thus improving the chances of obtaining stable numerical evolutions.

In Sec. II, we introduce the concept of quasi well-posedness. The explicit form of the equations and initial data used in our numerical evolutions are described in Sec. III. The prescription to adjust the standard ADM equations is presented in Sec. IV. We investigate in Sec. V the implications of the presence of gauge modes. Next, in Sec. VI, we discuss the quasi well-posedness properties of the adjusted ADM system under three lapse and shift conditions. For each of these conditions, we present numerical results showing the stability properties of the simulations. We end in Sec. VII with concluding remarks.

## II Quasi Well-Posedness

To understand the general idea of quasi well-posedness, let us consider the initial value problem for the wave equation $`_{tt}\phi =_{xx}\phi `$, ignoring boundaries. By adding a new variable, $`\pi =_t\phi `$, this equation takes the form

$$_tu=P_{xx}u+Qu,$$
(1)
where

$$u=\left(\begin{array}{c}\phi \\ \pi \end{array}\right),P=\left(\begin{array}{cc}0& 0\\ 1& 0\end{array}\right)\text{and}Q=\left(\begin{array}{cc}0& 1\\ 0& 0\end{array}\right).$$

One can view the system of equations (1) as the analogue of the standard ADM equations. Now, adding one more variable, $`\xi =_x\phi `$, we rewrite (1) as a first order system in both time and space. That is,

$$_t\stackrel{~}{u}=A_x\stackrel{~}{u}+B\stackrel{~}{u},$$
(2)
with

$$\stackrel{~}{u}=\left(\begin{array}{c}\phi \\ \pi \\ \xi \end{array}\right),A=\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 1\\ 0& 1& 0\end{array}\right)\text{and}B=\left(\begin{array}{ccc}0& 1& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right).$$

Notice that one needs to supplement the system (2) with the the constraint $`C=\xi _x\phi `$. For this simple case, one can find estimates for the solution $`u`$ from (1) directly, but the idea here is not to do so but rather to use what is known about the first order system (2). The matrix $`A`$ has the following set of eigenvalues and eigenvectors:

$$\begin{array}{ccc}\lambda _1\hfill & =1\hfill & \text{ with eigenvector}\stackrel{}{e}_1=[0,1,1]\hfill \\ \lambda _2\hfill & =0\hfill & \text{ with eigenvector}\stackrel{}{e}_2=[1,0,0]\hfill \\ \lambda _3\hfill & =1\hfill & \text{ with eigenvector}\stackrel{}{e}_3=[0,1,1].\hfill \end{array}$$
(3)
Given that the matrix $`A`$ has real eigenvalues and a complete system of eigenvectors, the system of equations (2) is strongly hyperbolic ***Following , a system is called symmetric hyperbolic if the matrix $`A`$ is Hermitian. It is called strictly hyperbolic if the eigenvalues are real and distinct; it is called strongly hyperbolic if the eigenvalues are real and there exists a complete system of eigenvectors; and, finally, it is called weakly hyperbolic if the eigenvalues are real. and thus well-posed. In this case, the matrix $`A`$ is also Hermitian (symmetric hyperbolicity). and has real and distinct eigenvalues (strictly hyperbolic)

Since the system (2) is well-posed, we can establish existence and uniqueness of the solution to (2), together with the usual condition:

$$||\stackrel{~}{u}(t,)||^2De^{at}||\stackrel{~}{u}(0,)||^2,$$
(4)
with $`a`$ and $`D`$ independent of the initial data and $`||||`$ the $`L_2`$ norm. Furthermore, because the vector $`\stackrel{~}{u}`$ contains $`u`$ as part of its components,

$$||u(t,)||^2||\stackrel{~}{u}(t,)||^2.$$

It follows then from equation (4) that the solution of the system (1) satisfies the condition

$$||u(t,)||^2De^{at}(||u(0,)||^2+||_xu(0,)||^2).$$
(5)
We should remark that, although (5) was derived using $`L_2`$ norms, the solution could in general be bounded by the initial data in some other Sobolev norm. The importance of the condition (5) and similar bounds is that it guarantees, among other things, that the solutions to (1) will not have unbounded high frequency modes. In numerical calculations however, this is not enough. For practical purposes, even large power-law growths such as $`t^4`$ are likely to be extremely problematic .

The idea is then to investigate quasi well-posedness in the standard ADM formulation. That is, given a choice of lapse and shift, one enlarges the standard ADM system introducing new variables that render the system at least strongly hyperbolic. Existence, uniqueness and well-posedness for the ADM equations will then follow from the corresponding usual properties of the strongly hyperbolic system. In some instances, to obtain such hyperbolic formulations it is not enough just to add the spatial derivatives of the metric as new variables, but one also has to add the constraints to the evolution equations, or even take combinations of these evolution equations. But the general idea is the same: add the constraints and/or take the same combination in the ADM equations as one does in the strongly hyperbolic system. It will be in this sense that we will refer, in what follows, to the quasi well-posedness of some ADM systems.

It is important to mention that for this analysis we have neglected the effects from boundary conditions. For the case under consideration (single black hole evolutions via excision of the singularity) boundary conditions at the excised region do not affect quasi well-posedness since, as we shall see, all the characteristics of the system are outgoing into the hole. At the outer edge of the computational domain, we use the exact analytic solution as boundary condition. Our numerical results show that this does not influence the properties drawn from quasi well-posedness for the class of gauge choices and added constraint-terms we considered.

## III The $`\dot{g}`$$`\dot{K}`$ Equations

Although the idea of quasi well-posedness can be applied to more general systems, throughout this paper we will concentrate on the ADM evolution equations ($`\dot{g}\dot{K}`$ equations) obtained from the vanishing of the Ricci tensor, i.e.

$`_tg_{ij}_\beta g_{ij}`$ $`=`$ $`2\alpha K_{ij},`$ (6)
$`_tK_{ij}_\beta K_{ij}`$ $`=`$ $`_i_j\alpha +\alpha (R_{ij}+KK_{ij}2K_{ik}K_j^k)`$ (7)
As customary in numerical relativity, we will focus on free evolutions; that is, we will not enforce the constraints

$`2H`$ $``$ $`R+K^2K_{ij}K^{ij}=0`$ (8)
$`M^i`$ $``$ $`_j(K^{ij}g^{ij}K)=0,`$ (9)
except at the initial data. Finally, to facilitate our analysis, we will only consider a single black hole in spherical coordinates. The most general time-dependent spherically symmetric metric in a 3+1 form is:

$`ds^2=(\alpha ^2+a^2\beta ^2)dt^2+2a^2\beta dtdr+a^2dr^2+b^2d\mathrm{\Omega }^2,`$ (10)
where all functions are assumed to depend only on $`r`$ and $`t`$; $`d\mathrm{\Omega }^2d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2`$, and $`\beta \beta ^r`$. In these spherical coordinates, the spatial metric, extrinsic curvature and Ricci tensor are diagonal:

$`g_{ij}`$ $`=`$ $`\text{diag}(a^2,b^2,b^2\mathrm{sin}^2\theta )`$ (11)
$`K_j^i`$ $`=`$ $`\text{diag}(K_a,K_b,K_b)`$ (12)
$`R_j^i`$ $`=`$ $`\text{diag}(R_a,R_b,R_b).`$ (13)
The $`\dot{g}`$$`\dot{K}`$ equations take the form

$`\left(_t\beta _r\right)a`$ $`=`$ $`\alpha aK_a+a_r\beta `$ (14)
$`\left(_t\beta _r\right)b`$ $`=`$ $`\alpha bK_b`$ (15)
$`\left(_t\beta _r\right)K_a`$ $`=`$ $`{\displaystyle \frac{1}{a^2}}\left(_r^2\alpha {\displaystyle \frac{1}{a}}_ra_r\alpha \right)+\alpha \left[R_a+\left(K_a+2K_b\right)K_a\right]`$ (16)
$`\left(_t\beta _r\right)K_b`$ $`=`$ $`{\displaystyle \frac{1}{ba^2}}_rb_r\alpha +\alpha \left[R_b+\left(K_a+2K_b\right)K_b\right],`$ (17)
and the constraints

$`H`$ $``$ $`{\displaystyle \frac{R_a}{2}}+R_b+2K_aK_b+K_b^2=0`$ (18)
$`M`$ $``$ $`_rK_b+(K_bK_a){\displaystyle \frac{_rb}{b}}=0,`$ (19)
where the components of the Ricci tensor are given by

$`R_a`$ $`=`$ $`{\displaystyle \frac{2}{a^3b}}\left(a_r^2b+_ra_rb\right)`$ (20)
$`R_b`$ $`=`$ $`{\displaystyle \frac{1}{a^3b^2}}\left[ba_r^2b+b_ra_rb+a^3a\left(_rb\right)^2\right].`$ (21)
Our numerical results consist of free-evolutions of analytic, single black hole initial data $`(a,b,K_a,K_b)`$ in a computational domain $`r_err_o`$, where $`r_o`$ denotes the location of the outer boundary far from the hole and $`r_e`$ the excision boundary inside the black hole horizon. The numerical code used solves the $`\dot{g}\dot{K}`$ equations in the interior of the computational domain by the method of lines. That is, the $`\dot{g}\dot{K}`$ system has the generic form

$$_tu=\beta _ru+\rho ,$$
(22)
with $`u=u(a,b,K_a,K_b)`$ and $`\rho `$ given by the r.h.s. of equations (1417). The starting point is then to approximate the derivatives appearing in $`\rho `$ by second order accurate, centered finite differences. On the other hand, the the spatial derivative in the “advection” term $`\beta _ru`$ is approximated by an upwind discretization  of the form:

$$_ru_i\mathrm{\Delta }_ru_i\frac{(u_{i+1}u_{i1})}{2\mathrm{\Delta }r}\pm q\frac{(u_{i1}3u_i+3u_{i\pm 1}u_{i\pm 2})}{3\mathrm{\Delta }r}$$
(23)
where $`q0`$ and the choice of sign is given by $`\pm =\text{sign}(\beta )`$. As we shall later see, for the black hole spacetime metrics under consideration, the shift vector is always non-negative. Thus, the discretization (23) is only needed with the upper sign. The above upwind discretization has truncation error:

$$\tau _i\mathrm{\Delta }_ru_i_ru_i=\frac{1}{6}(12q)\mathrm{\Delta }r^2_r^3u_i\frac{q}{6}\mathrm{\Delta }r^3_r^4u_i.$$
(24)
Thus, the discretization is second order accurate for any choice of the adjustable parameter $`q`$, except for $`q=0.5`$ when the truncation errors are of $`O(\mathrm{\Delta }r^3)`$. In our simulations, we set $`q=0.5`$. This value yields both stable evolutions and minimizes the amount of dissipation introduced by the discretization. After the discretization of the spatial differential operators, Eq. (22) becomes a set of ordinary differential equations for the interior grid points $`\{r_i\}_{i=1,\mathrm{},N1}`$. The grid values $`r_0`$ and $`r_N`$ denote the locations of the inner and outer boundaries, respectively. The temporal updating of these equations is carried out via an iterated Crank Nicholson method . Notice that at the outer-most grid point $`r_{N1}`$, the upwind scheme (23) requires “ghost” values at $`r_{N+1}`$. These values are computed using second or third order accurate extrapolations.

Next is the implementation of boundary conditions. As mentioned before, there are two boundaries in the problem, one far from the black hole at $`r_N=r_o`$ and a second at the point of excision, $`r_0=r_e`$. The boundary condition used at $`r_o`$ is that all field variables take the values provided by the exact analytic solution. It is important to notice that in previous efforts in numerical evolutions of spherically symmetric black hole spacetimes, the stability of the evolutions depended not only in the location of the outer boundary but also on the imposition of “outgoing” boundary conditions. One of these outgoing conditions is constructed by assuming that the solution error $`u\overline{u}`$ is of the form $`ϵ(tr)/r^n`$, with $`\overline{u}`$ the analytic exact solution. This condition yields

$$_tu=_ru+_r\overline{u}+n\frac{(u\overline{u})}{r},$$
(25)
to be applied at the nodal grid point $`r_N=r_o`$ The outgoing condition (25) is then added to the system (22) and handled with the same method of lines used for the interior points. Obviously, there is no reason to believe that the solution error should behave as $`ϵ(tr)/r^n`$ for general choices of lapse and shift. Another approach to construct outer boundary conditions that has been used is to blend  the numerical solution to the analytic one beyond a certain radius, $`r`$.

Finally, at the inner or excision boundary, the working assumption is that all of the fields have outgoing (into the hole) characteristics. We will show that for the gauge or coordinate conditions under consideration, this is indeed the situation. Therefore, there is no need to impose any boundary condition, and one can just apply at $`r_e`$ the same system of equations used in the interior of the computational domain. However, finite difference discretization would require “ghost” values at $`r_1`$. We construct these values by extrapolation. Another possibility, which is becoming popular in three-dimensional simulations , is instead to extrapolate the r.h.s. of Eq. (22) to $`r_0=r_e`$. Either approach was stable for the cases considered in the present work.

An important requirement when performing numerical simulations of black holes in which the singularities are excised from the computational domain is to use coordinates that are regular and penetrate the horizon. Foliations that penetrate the horizon facilitate the task of removing (i.e. excising) a region containing the black hole singularity while preserving the causal structure of the spacetime exterior to the event horizon. Thus the numerical results we present here consist of numerical evolutions that for infinite resolution correspond to the solution of a single black hole expressed in ingoing-Eddington-Finkelstein (iEF) and also Painlevé-Gullstrand (PG) . coordinates. The iEF coordinates coincide with Kerr-Schild coordinates in the case of zero angular momentum . Recently, Martel and Poisson  have showed that the PG and iEF coordinates are members of the same one-parameter family of coordinate systems.

In iEF coordinates, the line element (10) takes the explicit form

$$ds^2=\left(1\frac{2m}{r}\right)dt^2+\frac{4m}{r}dtdr+\left(1+\frac{2m}{r}\right)dr^2+r^2d\mathrm{\Omega }^2.$$
(26)
The ADM variables in these coordinates are given by:

$`\alpha `$ $`=`$ $`\left(1+{\displaystyle \frac{2m}{r}}\right)^{1/2}`$ (27)
$`\beta `$ $`=`$ $`{\displaystyle \frac{2m}{r}}\left(1+{\displaystyle \frac{2m}{r}}\right)^1`$ (28)
$`a`$ $`=`$ $`\left(1+{\displaystyle \frac{2m}{r}}\right)^{1/2}`$ (29)
$`b`$ $`=`$ $`r`$ (30)
$`K_a`$ $`=`$ $`{\displaystyle \frac{2m}{r^3}}(r+m)\left(1+{\displaystyle \frac{2m}{r}}\right)^{3/2}`$ (31)
$`K_b`$ $`=`$ $`{\displaystyle \frac{2m}{r^2}}\left(1+{\displaystyle \frac{2m}{r}}\right)^{1/2}.`$ (32)
Similarly, for PG coordinates

$$ds^2=\left(1\frac{2m}{r}\right)dt^2+2\sqrt{\frac{2m}{r}}dtdr+dr^2+r^2d\mathrm{\Omega }^2,$$
(33)
and

$`\alpha `$ $`=`$ $`1`$ (34)
$`\beta `$ $`=`$ $`\sqrt{{\displaystyle \frac{2m}{r}}}`$ (35)
$`a`$ $`=`$ $`1`$ (36)
$`b`$ $`=`$ $`r`$ (37)
$`K_a`$ $`=`$ $`{\displaystyle \frac{\beta }{2r}}`$ (38)
$`K_b`$ $`=`$ $`{\displaystyle \frac{\beta }{r}}.`$ (39)
The geometrical interpretation of the iEF coordinate system is that in addition to having a timelike killing vector, the combination of timelike and radial tangent vectors $`\stackrel{}{}_t\stackrel{}{}_r`$ remains null. In terms of the spacetime metric, this condition is stated as $`g_{tt}2g_{tr}+g_{rr}=0`$, or similarly in terms of 3+1 metric functions as $`\alpha =a(1\beta )`$. On the other hand, the PG coordinate system can be viewed as that anchored to a family of freely moving observers (time-like) starting at infinity with vanishing velocity .

Finally, in numerical evolutions of spherically symmetric spacetimes, it is useful to monitor not only the Hamiltonian and momentum constraints but to pay attention also to the mass function

$$(r,t)=\frac{b}{2}(1_\mu b^\mu b).$$
(40)
In vacuum, this $``$ is the gauge invariant definition of mass.

## IV Adjusted ADM Systems

Each of the equations in the $`\dot{g}`$$`\dot{K}`$ or standard ADM system (14-17) has the form

$$_tu_{(n)}\beta _ru_{(n)}=B_{(n)}u_{(n)}+C_{(n)}$$
(41)
(no summation over $`n`$) with

$`u_{(n)}`$ $`=`$ $`(a,b,K_a,K_b)`$ (42)
$`B_{(n)}`$ $`=`$ $`\text{diag}(\alpha K_a+_r\beta ,\alpha K_b,\mathrm{\hspace{0.17em}2}\alpha K_b,\alpha K_a)`$ (43)
The $`C_{(n)}`$ contain the remaining terms that cannot be written in the form $`B_{(n)}u_{(n)}`$, with $`B_{(n)}`$ independent of $`u_{(n)}`$. If one views each equation as independent, that is $`\beta `$, $`B_{(n)}`$ and $`C_{(n)}`$ as given and independent of $`u_{(n)}`$, Eqs. (41) admit rapidly growing solutions if $`B_{(n)}>0`$, exponentially growing if $`B_{(n)}`$ and $`\beta `$ are non-negative constants. Obviously, when the equations (41) are considered as a coupled system of equations, one cannot guarantee the above conclusion. However, we have been able to track the problems observed in our numerical evolutions to those terms for which $`B_{(n)}>0`$.

In iEF and PG coordinates, the second and fourth components of $`B_{(n)}`$ are non-positive. The first component of $`B_{(n)}`$ vanishes for PG coordinates and is positive for iEF. Finally, the third component of $`B_{(n)}`$ is always non-negative, for both iEF and PG coordinates. Our numerical experiments indicate that the origin of the instabilities is due to the term involving $`B_{(n=3)}=2\alpha K_b`$. The term in iEF coordinates involving $`B_{(n=1)}`$ with the “wrong” sign did not seem to affect the evolutions (see results below).

The objective is then to find a way to “change the sign” of $`B_{(n=3)}`$. Fortunately, the combination $`K_bK_a`$ also appears in the Hamiltonian constraint (18). One can then add to the evolution equation for $`K_a`$ a term of the form $`\mu \alpha H`$. As a consequence, $`B_{(n=3)}=2(1\mu )\alpha K_b`$. Therefore, in principle, any value of $`\mu 1`$ should produce stable evolutions. Our numerical simulations show that $`\mu =2`$ is an optimal value to reach quickly the time independent solution. However, the same experiments indicate also that even values of $`\mu 0.5`$ yield stable evolutions. The reason for this is likely because the simple analysis above does not take into consideration the non-linear coupling in the equations.

## V Linear Gauge Modes

A potential source of instabilities in numerical simulations is the presence of gauge modes. Gauge modes arise because a prescription of the lapse and shift is not sufficient for complete gauge fixing. If the lapse and shift functions are for instance determined from algebraic expressions or differential equations, there is still remaining freedom left to perform coordinate transformation that take us from $`(g_{ij},K_{ij})`$ to $`(\stackrel{~}{g}_{ij},\stackrel{~}{K}_{ij})`$, leaving the lapse-shift prescription invariant. In other words, the data $`(g_{ij},K_{ij})`$ will be unique up to coordinate transformations that leave the lapse-shift prescription invariant. One can then encounter the situation in which the transformed pair $`(\stackrel{~}{g}_{ij},\stackrel{~}{K}_{ij})`$ continues to satisfy the constraints, but it possesses unbounded growths. In principle, because $`(\stackrel{~}{g}_{ij},\stackrel{~}{K}_{ij})`$ are allowed solutions to the evolution equations, one could through numerical convergence and monitoring of the constraints single out these modes. In practice the situation is not that simple. Gauge conditions that allow the rapid growth of grid functions are likely to trigger numerical instabilities.

To investigate these modes, let us consider an infinitesimal coordinate transformation,

$$x^\mu x^\mu +\xi ^\mu ,$$
(44)
This transformation induces a linear perturbation of the spacetime metric,

$$g_{\mu \nu }g_{\mu \nu }+\delta g_{\mu \nu },$$
(45)
of the form

$$\delta g_{\mu \nu }=_\xi g_{\mu \nu }.$$
(46)
In spherical symmetry, the most general coordinate transformation is

$`t`$ $``$ $`t+\delta t(r,t)`$ (47)
$`r`$ $``$ $`r+\delta r(r,t),`$ (48)
or equivalently

$$\xi ^\mu =(\delta t,\delta r,0,0).$$
(49)
Given (49), the gauge induced perturbations of the metric takes the form

$`\delta g_{tt}`$ $`=`$ $`2\left(\alpha ^2+a^2\beta ^2\right)_t\delta t2a^2\beta _t\delta r_t\left(\alpha ^2+a^2\beta ^2\right)\delta t_r\left(\alpha ^2+a^2\beta ^2\right)\delta r`$ (50)
$`\delta g_{tr}`$ $`=`$ $`a^2\beta _t\delta ta^2_t\delta r\left(\alpha ^2+a^2\beta ^2\right)_r\delta ta^2\beta _r\delta r_t\left(a^2\beta \right)\delta t_r\left(a^2\beta \right)\delta r`$ (51)
$`\delta g_{rr}`$ $`=`$ $`2a^2\beta _r\delta t2a^2_r\delta r_t\left(a^2\right)\delta t_r\left(a^2\right)\delta r`$ (52)
$`\delta g_{\theta \theta }`$ $`=`$ $`_t\left(b^2\right)\delta t_r\left(b^2\right)\delta r`$ (53)
$`\delta g_{\varphi \varphi }`$ $`=`$ $`\delta g_{\theta \theta }\mathrm{sin}^2\theta ,`$ (54)
with all remaining components vanishing. On the other hand, from (10), the components of the perturbed metric $`\delta g_{\mu \nu }`$ above are also given by:

$`\delta g_{tt}`$ $`=`$ $`2\alpha \delta \alpha +2a\delta a\beta ^2+2a^2\beta \delta \beta `$ (55)
$`\delta g_{tr}`$ $`=`$ $`2a\delta a\beta +a^2\delta \beta `$ (56)
$`\delta g_{rr}`$ $`=`$ $`2a\delta a`$ (57)
$`\delta g_{\theta \theta }`$ $`=`$ $`2b\delta b.`$ (58)
In studying these gauge perturbations, we consider a number of different prescriptions for choosing the lapse and shift. These impose restrictions on the components of $`\delta g_{\mu \nu }`$ to linear order, but are not sufficient for complete gauge fixing.

One piece of evidence usually given to argue for the highly unstable properties of the standard ADM formulation is its inability to simulate single black hole spacetimes. Specifically, the problem consists of using a known analytic black hole solution of Einstein’s equation to construct initial data, set boundary conditions and specify the lapse and shift. Given this input, the output is the numerical evolution of the spatial metric $`g_{ij}`$ and extrinsic curvature $`K_{ij}`$. This simple setup does not yield stable evolutions.

Choosing the lapse and shift to be given by the exact analytic solutions implies that in the coordinate transformations above, $`\delta \alpha =\delta \beta =0`$. We shall call this the “exact lapse + exact shift” (EL+ES) condition, It is clear from (5054) that the EL+ES choice does not fix the coordinates uniquely, but rather results in a non-trivial system of PDE’s to be satisfied by $`\delta t`$ and $`\delta r`$. The solutions to this PDE system form an equivalence class of gauges all satisfying the same condition on the lapse and shift. At the analytic level, one selects a unique member of this equivalence class by imposing boundary and initial conditions. Numerically however, because of truncation errors, this is only possible if the evolution is stable.

Substitution of $`\delta \alpha =\delta \beta =0`$ in (55-57) yields $`\delta g_{tr}\beta \delta g_{rr}=\delta g_{tt}\beta ^2\delta g_{rr}=0`$. Using (5052), these conditions can be rewritten in terms of the gauge perturbation $`\xi ^\mu `$ as

$$_t\xi =A_r\xi +B\xi $$
(59)
where

$$\xi =\left(\begin{array}{c}\delta t\\ \delta r\end{array}\right)\text{and}A=\left(\begin{array}{cc}\beta & 0\\ \frac{\alpha ^2}{a^2}& \beta \end{array}\right).$$

The matrix $`A`$ has a degenerate eigenvalue $`\lambda =\beta `$ and corresponding eigenvector $`\stackrel{}{e}=[0,\mathrm{\hspace{0.17em}1}]`$. Therefore, the system is only weakly hyperbolic and thus ill-posed. It is then not possible to guarantee the absence of rapidly growing gauge modes. Gauge instabilities, by themselves, do not violate the constraints; however, they are likely to trigger numerical instabilities and thereby couple to constraint-violating modes. Note that this conclusion does not depend on the use of the standard ADM formulation, but it applies to any initial value formulation of Einstein’s equation. This is likely the reason why it has not been possible to produce hyperbolic formulations of Einstein equations using the EL+ES lapse-shift condition. Other groups have reported numerical instabilities associated with an EL+ES prescription. (See, for example, .) Our result, which provides some analytic insight into the source of these instabilities, extends the analysis of Ref.  by including perturbations of a general spherically symmetric spacetime (10) and by making no assumptions about the form of the coupling between gauge and constraint-violating perturbations.

Finally, the inability to guarantee the absence of unbound gauge modes does not necessarily imply that it is impossible to design an evolution scheme that is long-term stable and convergent. As we shall see, by adjusting the standard ADM system with constraint terms, stable evolutions are possible even in the presence of these gauge modes.

## VI Stable adjusted ADM systems in 1d

We consider next a series of lapse-shift choices. For each choice, we investigate (1) the quasi well-posedness of the resulting system of adjusted ADM evolution equations, (2) the propagation of the constraints and (3) the convergence and stability of numerical simulations.

### A Exact Lapse + Area Locking

We can take advantage of the assumed spherical symmetry of the problem and “lock” the area of constant-$`r`$ surfaces. That is, we exploit the lapse-shift freedom and set $`_tg_{\theta \theta }=0`$, or equivalently $`_tb=0t`$. From Eq. (15), this yields

$$0=\beta _rb+\alpha bK_b,$$
(60)
which can be seen as an algebraic equation to solve for $`\beta `$ or $`\alpha `$. We will use (60) as

$$\beta =\frac{\alpha bK_b}{_rb},$$
(61)
with $`b`$ determined by the initial data. In our case, for both iEF and PG coordinates, $`b=r`$. In addition, we choose an exact lapse, i.e. an arbitrary but a priori specified function of spacetime. Here again, since the goal is to reproduce numerically the analytic solution, we set the lapse to that given by the iEF or PG solutions.

This exact lapse, area locking (EL+AL) gauge condition has previously been investigated in Ref. . However, the implementation of locking the area was done at the numerical level. That is, condition (61) was not explicitly used. Instead, during the temporal updating of grid functions, a correction to the shift was introduced to keep the area locked. With this numerical area-locking and with blending outer boundary conditions, stable simulations were reported in Ref..  for computational domains with $`r_o11m`$.

Interestingly enough, this EL+AL choice of lapse and shift yields an ADM system of equations already first order in space and quasi-linear; namely,

$$_tu=A_ru+Bu$$
(62)
with

$$u=\left(\begin{array}{c}a\\ K_a\\ K_b\end{array}\right)\text{and}A=\left(\begin{array}{ccc}\beta & 0& \alpha ar\\ & & \\ \frac{r_r\alpha +2(1\mu )\alpha }{a^3r}& \beta & 0\\ & & \\ \frac{\alpha }{a^3r}& 0& \beta \end{array}\right),$$

where $`\beta `$ is given by (61) and we have set $`b=r`$ to simplify notation. In the numerical evolutions, however, we do not set $`b=r`$. The numerical code includes the evolution equation for $`b`$. We do not explicitly write the matrix $`B`$ since it is not necessary for the analysis below. However, we must emphasize that, in order to achieve stability, a value of $`\mu 1`$ is in principle needed (see Sec. IV).

Although this EL+AL case is not representative of the general ADM equations, where second spatial derivatives do appear, we will use it to introduce the main techniques and ideas regarding quasi well-posedness. We notice first that the matrix $`A`$ has eigenvalues

$`\lambda _1`$ $`=`$ $`\beta `$ (63)
$`\lambda _2`$ $`=`$ $`\beta +{\displaystyle \frac{\alpha }{a}}`$ (64)
$`\lambda _3`$ $`=`$ $`\beta {\displaystyle \frac{\alpha }{a}},`$ (65)
and corresponding eigenvectors

$`\stackrel{}{e}_1`$ $`=`$ $`[0,1,0]`$ (66)
$`\stackrel{}{e}_2`$ $`=`$ $`[a^2r,{\displaystyle \frac{r_r\alpha +2(1\mu )\alpha }{\alpha }},1]`$ (67)
$`\stackrel{}{e}_3`$ $`=`$ $`[a^2r,{\displaystyle \frac{r_r\alpha +2(1\mu )\alpha }{\alpha }},1].`$ (68)
Because all of the eigenvalues are real and distinct, the system (62) is strictly hyperbolic independent of the addition of the constraint term. The EL+AL system of equations is then an example of a hyperbolic system that, unless suitable constraint-terms are added, is subject to developing rapidly growing solutions (see Ref.  for another example). It is also important to notice that $`\lambda _1`$ represents a characteristic speed corresponding to propagation along the timelike normal to the foliation. Similarly, $`\lambda _2`$ and $`\lambda _3`$ are characteristic speeds along the light cone.

Let’s now consider the particular case of initial data and lapse function constructed from the single black hole solution in iEF or PG coordinates. The eigenvalues are given in these coordinates by

$`\lambda _1`$ $`=`$ $`{\displaystyle \frac{2m}{r+2m}}=\sqrt{{\displaystyle \frac{2m}{r}}}`$ (69)
$`\lambda _2`$ $`=`$ $`1=\sqrt{{\displaystyle \frac{2m}{r}}}+1`$ (70)
$`\lambda _3`$ $`=`$ $`{\displaystyle \frac{2mr}{2m+r}}=\sqrt{{\displaystyle \frac{2m}{r}}}1,`$ (71)
where the first and second equalities are for iEF and PG coordinates, respectively. Since the excision boundary $`r_e`$ is by construction inside the black hole horizon (i.e. $`r_er_h2m`$), we have that all the eigenvalues are positive there. By looking at the principal part of Eq. (62), it is easy to see that non-negative eigenvalues imply propagation of field variables in the direction of decreasing $`r`$-coordinate. Therefore, at $`r_e`$ all the fields propagate out of the computational domain into the hole singularity, thus boundary conditions are not required. On the other hand, at the outer boundary $`r_o`$, one has $`\lambda _3<0`$, with the remaining eigenvalues still non-negative. At $`r_o`$, one has then two ingoing modes ($`\lambda _1`$ and $`\lambda _2`$) and only one outgoing ($`\lambda _3`$). Given this information, it is perfectly possible to impose a condition suppressing modes entering the computational domain. However, as we mentioned before, we choose not to do so and impose at $`r_o`$ the analytic, exact solutions. The direct consequence will be the appearance of a pulse at the outer boundary of the computational domain due to discontinuities in the truncation errors between the outermost evolved point and the boundary point $`r_o`$. This pulse propagates in the direction of the black hole and leaves the computational domain through the excision boundary.

The next step is to analyze the effect that the EL+AL choice has on the propagation of the constraints. At the continuum level, for arbitrary choice of lapse, shift and initial data ($`g_{ij},K_{ij}`$) satisfying the constraints (8) and (9), the evolution equations (6) and (7) guarantee, ignoring for the moment boundary conditions, that the evolved data will continue to satisfy the constraints. If one now takes into consideration boundary conditions, it is important to keep in mind that boundary data ($`g_{ij},K_{ij}`$) must satisfy the constraints. By looking at the way the constraints propagate, i.e. their characteristics, one gains insight about the allowed boundary conditions consistent with the constraints. Another important aspect of well-posedness in the propagation of the constraints is that it guarantees that there will be no unbounded high frequency growth appearing in the constraints if they are not exactly satisfied at the initial slice (for example, due to numerical errors). This well-posedness for the constraint propagation is a non-trivial property, not possible to prove for generic formulation of Einstein’s equations .

What we look for are evolution equations for the constraints, equations that would hold if the system (62) is satisfied. They can be found by taking time derivatives in both sides of equations (18) and (19), replacing the time derivatives of the metric by the r.h.s. of equations (1417), and finally expressing the metric and its spatial derivatives in terms of the constraints and their spatial derivatives. Following this procedure, it is not too difficult to show that

$$_tv=P_rv+Qv,$$
(72)
where now

$$v=\left(\begin{array}{c}H\\ M\end{array}\right)\text{and}P=\left(\begin{array}{cc}\beta & \frac{4\alpha }{a^2}\\ & \\ \frac{\alpha }{4}& \beta \end{array}\right).$$
(73)
The matrix $`P`$ has eigenvalues $`\overline{\lambda }_1=\lambda _2,\overline{\lambda }_2=\lambda _3`$ , with $`\lambda _2`$ and $`\lambda _3`$ given by (64) and (65), respectively. This implies that the system (72) is also strictly hyperbolic with characteristic speeds along the light cone. Notice that the constraints at $`r_e`$ propagate out of the computational domain into the hole singularity, consistent with the outgoing propagation of the field variables at $`r_e`$, namely the tilting, into the black hole, of the light cone. At the outer boundary $`r_o`$, there is an ingoing mode ($`\overline{\lambda }_2>0`$). Therefore, one has, as expected, to be careful to provide data at $`r_o`$ consistent with this entering mode. Since we are imposing at the outer boundary the analytic iEF and PG solutions, the data at $`r_o`$ already satisfy the constraints. However, as mentioned before, choosing the lapse and shift does not completely fix the gauge freedom, thus one still has to be careful handling the gauge modes described in Sec. V.

Figure 1 shows the $`L_2`$ norm of the Hamiltonian constraint as a function of time. The initial data is given by the iEF analytic solution (26), and the lapse and shift are chosen from the EL+AL conditions. Similar results were obtained with PG coordinates. The computational domain extends from $`r_e=1m`$ to $`r_o=40m`$. We tried larger and smaller values for $`r_o`$. However, the stability of the simulations was not affected by the location of the outer boundary. We use an upwind parameter $`q=0.5`$ and a constraint-term parameter $`\mu =2`$. We show runs for resolutions of $`\mathrm{\Delta }r=m/5,m/10,m/20,m/40`$ with $`\mathrm{\Delta }t=0.25\mathrm{\Delta }r`$. The run with $`\mathrm{\Delta }r=m/5`$ has a resolution similar to those used in three-dimensional simulations of black hole collisions.

Figure 2 shows the $`L_2`$ norm of the Hamiltonian constraint (top) and the $`L_2`$ norm of the mass (40) error (bottom) taken at time $`t=200m`$ as a function of resolution. The convergence rate implied by the Hamiltonian constraint is 2.18 and by the mass error 1.7. The convergence rate from the Hamiltonian constraint is larger than second order because we used third order accurate discretizations of the advection term as well as third order accurate extrapolations at the excision. On the other hand, the convergence rate obtained from the mass error is less than second order because the mass function is proportional to $`g^{rr}`$, a quantity difficult to handle numerically near the singularity. The reason for using the Hamiltonian constraint and mass to monitor accuracies and convergence is because of the gauge invariant nature of these quantities. Finally for reference, we show in Fig. 3 the $`L_2`$ norm of the solution Hamiltonian constraint as a function of time for runs with resolution $`\mathrm{\Delta }r=m/10`$ varying the parameter $`\mu `$. Different lines correspond to values of $`\mu =0.0,\mathrm{\hspace{0.17em}0.025},\mathrm{\hspace{0.17em}0.5},`$ and 2 in order of stability improvement. It is clear from this figure the dramatic effect that the added constraint term has on the stability of the simulations.

### B Ingoing Null + Area Locking

The ingoing null + area locking (IN+AL) recipe to specify the lapse and shift consists of, in addition to locking the areal coordinate, imposing the condition that the vector $`_t_r`$ remains null throughout the evolution. This null condition is only compatible with the iEF case since by construction the iEF solution is based on ingoing null observers. An analogous (ingoing timelike) condition can be obtained for the PG case. In terms of the spacetime metric, the ingoing null condition is stated as $`g_{tt}2g_{tr}+g_{rr}=0`$, or similarly in terms of 3+1 metric functions in (10) as

$$\alpha =a(1\beta ).$$
(74)
Conditions (74) and (60) yield

$$\alpha =\frac{a}{arK_b+1}\text{and}\beta =\frac{arK_b}{arK_b+1},$$
(75)
where we have set $`b=r`$ since by construction $`b`$ remains locked to $`r`$. This prescription for the lapse and shift has been successfully applied in the past .

In order to make the ADM equations in the IN+AL gauge a first order in space system, we need to introduce two new variables: $`w_ra`$ and $`y_rK_b`$. Even after this choice is made, there is no unique way of writing the resulting equations as a quasi-linear system. The reason is the ambiguity one encounters when dealing with the terms involving $`_ra`$. One has the choice to either keep it as $`_ra`$, substitute it with $`w`$, or a combination of both. Either choice changes the principal part of the equation. However, in our case, we need only to find a choice that yields a strongly hyperbolic system. It turns out that the simplest choice of replacing $`_ra`$ by $`w`$ everywhere yields a system that is well-posed. That is, the resulting ADM equations have the form

$$_tu=A_ru+Bu,$$
(76)
with

$$u=\left(\begin{array}{c}a\\ K_a\\ K_b\\ y\\ w\end{array}\right)$$
(77)
and

$$A=\frac{1}{z^2}\left(\begin{array}{ccccc}0& 0& 0& 0& 0\\ & & & & \\ 0& zarK_b& 0& r& a^2\\ & & & & \\ 0& 0& 0& 0& 0\\ & & & & \\ 0& zaK_b& 0& z+a^2r^2K_b^2& a^1K_b\\ & & & & \\ 0& za^2& 0& a^2r& arK_b(1+z)\end{array}\right),$$
(78)
where we have introduced $`z1+arK_b`$ to simplify notation. The corresponding eigenvalues and eigenvectors are

$`\lambda _1`$ $`=`$ $`0\text{ with eigenvectors}\stackrel{}{e}_1=[1,0,0,0,0],\stackrel{}{e}_2=[0,0,1,0,0],`$ (79)
$`\lambda _2`$ $`=`$ $`1\text{ with eigenvectors}\stackrel{}{e}_3=[0,1,0,0,za^2],\stackrel{}{e}_4=[0,0,0,1,a^2r],`$ (81)
$`\lambda _3`$ $`=`$ $`{\displaystyle \frac{arK_b1}{z}}\text{ with eigenvector}\stackrel{}{e}_5=[0,1,0,aK_b,a^2].`$ (83)
Notice that all the eigenvectors are independent, and, thus, the system is strongly hyperbolic. Also notice that the eigenvalues $`\lambda _2`$ and $`\lambda _3`$ are again the characteristic speeds ($`\beta \pm \alpha /a`$) along the light cone. Furthermore, in iEF coordinates, $`\lambda _3=(2mr)/(2m+r)`$. Therefore, one encounters a similar situation to that of the EL+AL case; namely, at $`r_e`$ all the eigenvalues are non-negative, and at $`r_o`$ one has that only $`\lambda _3<0`$. The existence, uniqueness, and well-posedness for IN+AL follows then as with EL+AL case.

Regarding the constraints, their evolution is also described by a strongly hyperbolic system with characteristic speeds along the light cone. It is important to stress that our ADM equations already imply this, i.e. it does not have any relation to making the evolution equations first order in space. In fact, the principal part of the evolution equations for the constraints is exactly the same as in the EL+AL gauge, but now the lapse and shift are given by (75). Thus, the analysis and conclusions also follow as in that case.

An interesting aspect of the IN+AL choice is that it is possible to find the general solution to Einstein’s equations. We start by defining $`fraK_b`$. We then use the momentum constraint (19) to eliminate $`K_a`$ from the other equations. The outcome is that we need only to solve three of the four equations (14), (16-18). We choose to work with the Hamiltonian constraint, (18), and Eqs. (14) and (17). The resulting system of equations reads

$`0`$ $`=`$ $`a(f^2+a^21)+2raf_rf+2r(1f^2)_ra`$ (84)
$`(_t_r)a`$ $`=`$ $`{\displaystyle \frac{a(f^2+a^21)}{2r(1+f)^2}}`$ (85)
$`(_t_r)f`$ $`=`$ $`{\displaystyle \frac{(f^2+a^21)}{r(1+f)}}.`$ (86)
The general solution of these equations is

$`f`$ $`=`$ $`{\displaystyle \frac{2m}{r}}C\left(1{\displaystyle \frac{2m}{r}}\right)`$ (87)
$`a^2`$ $`=`$ $`(1+C)(1+f)`$ (88)
where $`C=C(t+r)`$ and $`m`$ a constant. In terms of the 3+1 variables, the solutions reads

$$a^2=(1+C)(1+f),b^2=r^2,K_a=_r\left(\frac{f}{a}\right),K_b=\frac{f}{ar},\alpha =\frac{a}{1+f},\beta =\frac{f}{1+f},$$
(89)
with the line element (10) explicitly given by

$`ds^2`$ $`=`$ $`(1+C)^2\left(1{\displaystyle \frac{2m}{r}}\right)dt^2+\mathrm{\hspace{0.17em}2}(1+C)\left[{\displaystyle \frac{2m}{r}}\left(1{\displaystyle \frac{2m}{r}}\right)C\right]dtdr`$ (91)
$`+(1+C)\left[1+{\displaystyle \frac{2m}{r}}\left(1{\displaystyle \frac{2m}{r}}\right)C\right]dr^2+r^2d\mathrm{\Omega }^2.`$
By setting $`C=0`$ one recovers the iEF solution (26). Also, it is not difficult to show from the gauge invariant definition of mass (40) that the parameter $`m`$ is indeed the mass of the black hole. An important property of the general solution (91) is that it shows explicitly the residual gauge freedom associated with the IN+AL choice of lapse and shift. We have found explicitly the equivalence-class of solutions gauge related to the iEF solutions that satisfy the IN+AL lapse-shift condition.

We repeat the same type of numerical experiment as with the EL+AL case, same parameter values, resolutions, boundary conditions and initial data. The lapse and shift are however constructed in this case from (75). Figure 4 shows the $`L_2`$ norm of the Hamiltonian constraint as a function of time for different resolutions. Figure 5 shows the $`L_2`$ norms of the Hamiltonian constraint (top) and mass error (bottom) for different resolutions. The convergence rate are similar to those in EL+AL, namely 2.19 from the Hamiltonian constraint and 1.7 from the mass error. Finally, Fig. 6 shows the $`L_2`$ norm of the Hamiltonian constraint as a function of time for a resolution of $`\mathrm{\Delta }r=m/10`$. Different lines correspond to values of $`\mu =0.0,\mathrm{\hspace{0.17em}0.025},\mathrm{\hspace{0.17em}0.5},`$ and 2 in order of stability improvement. It is clear from this results that the stability behavior of the system of equations under the IN+AL gauge choice follows closely that of EL+AL.

### C Exact Lapse + Exact Shift

Finally, we consider the case in which the lapse function and shift vector are prescribed by the exact analytic solutions. To date, it has not been possible to obtain with the standard ADM system long-term stable and convergent numerical evolutions under the EL+ES choice. Before we present results of EL+ES evolutions using the adjusted ADM system of equations, let us investigate its quasi well-posedness properties. To make the evolution equations first order in space, we add a new variable $`y=_rb`$. Once more, the resulting set of equations have the form

$$_tu=A_ru+Bu,$$
(92)
with

$$u=\left(\begin{array}{c}a\\ b\\ K_a\\ K_b\\ y\end{array}\right)$$
(93)
and

$$A=\left(\begin{array}{ccccc}\beta & 0& 0& 0& 0\\ & & & & \\ 0& 0& 0& 0& 0\\ & & & & \\ \frac{b_r\alpha +2(1\mu )\alpha y}{a^3b}& 0& \beta & 0& \frac{2(1\mu )\alpha }{a^2b}\\ & & & & \\ \frac{\alpha y}{ba^3}& 0& 0& \beta & \frac{\alpha }{a^2b}\\ & & & & \\ 0& 0& 0& \alpha b& \beta \end{array}\right).$$
(94)
The matrix $`A`$ has eigenvalues

$`\lambda _1`$ $`=`$ $`\beta `$ (95)
$`\lambda _2`$ $`=`$ $`\beta `$ (96)
$`\lambda _3`$ $`=`$ $`0`$ (97)
$`\lambda _4`$ $`=`$ $`\beta +\alpha /a`$ (98)
$`\lambda _5`$ $`=`$ $`\beta \alpha /a.`$ (99)
There are two fields that propagate with characteristic speed ($`\beta `$) along the normal to the hypersurfaces in the foliation, one field with zero speed and two other fields with characteristic speeds ($`\beta \pm \alpha /a`$) along the light cone. It is not difficult to show that the eigenvectors corresponding to the eigenvalue $`\lambda _1=\lambda _2=\beta `$ are not distinct, whatever the value of $`\mu `$. Therefore, the system of equations is only weakly hyperbolic and thus not quasi well-posed, as we have already seen by considering the EL+ES gauge perturbations by themselves in Sec. V.

As in the previous two cases, the the constraints propagate according to

$$_tv=P_rv+Qv,$$
(100)
where now

$$v=\left(\begin{array}{c}H\\ M\end{array}\right)$$
(101)
and

$$P=\left(\begin{array}{cc}\beta & \frac{4\alpha }{a^2}\\ & \\ \frac{\alpha }{4}& \beta \end{array}\right)$$

This is, we have the same situation as before, namely propagation of the constraints along with characteristic speeds along the light cone.

We have carried out the same numerical experiments as for the previous two cases. However, for the iEF solution, the simulations with computational domains for $`r_o40m`$ crashed. We have been able to identify two possible sources behind this problem. One of them is the pulse originated at the outer boundary due to discontinuities of the truncation error. This pulse propagates in the direction of the black hole (i.e. decreasing $`r`$-coordinate) with the characteristic speed $`\beta `$. As the pulse travels its amplitude grows in time. This effect is shown in Fig. 7. Here we plot the solution error for the metric function $`a`$ as a function of space and time for a resolution $`\mathrm{\Delta }r=m/10`$ and outer boundary located at $`r_o=40m`$. The simulation stops because this pulse-error increases to the point that the metric function $`a`$ becomes negative. For small computational domains ($`r_o<40m`$), the crossing time of this pulse is short enough and does not allow a catastrophic growth of the pulse. The pulse is able to leave the computational domain through the excision boundary without crashing the simulation. Since the initial amplitude of this pulse is $`O(\mathrm{\Delta }r^2)`$, i.e. the accuracy of the discretization, in principle one could find a resolution small enough such that the growth of the pulse does not affect the life of the simulation. However, accessing those fine resolutions in three-dimensional simulations is likely to be impractical. The second, and perhaps more severe, source of the problem is the presence of the zero velocity mode. Zero velocity modes were in principle also allowed in the case of IN+AL, see Eq. (79). However, the adjusted IN+AL system yields only numerical solutions of the general form (91), which clearly does not contain a zero velocity mode. If $`\mu =0`$, these zero velocity modes are not suppressed and eventually terminate the simulation. This catastrophic effect induced by zero velocity modes has been previously noticed by Alcubierre et.al .

Figure 8 shows the $`L_2`$ norm of the Hamiltonian constraint as a function of time with resolution $`\mathrm{\Delta }r=m/10`$ for different locations of the outer boundary. For computational domains with approximately $`r_o40m`$, the zero velocity mode is still present, but is damped eventually. The reasons why this mode stops growing remain unclear. Nonetheless, there is a strong indication that this behavior is connected to the particular choice of coordinates used to set the exact lapse and shift. If one, instead of the iEF solution, sets the lapse and shift from the PG solution, the outcome of the simulations is completely different. EL+ES simulations with PG lapse and shift are long term stable and convergent for arbitrary sizes of the computational domain as long as the system of equations was adjusted with $`\mu 1`$. Figure 9 shows the $`L_2`$ norm of the solution error for the metric function $`a`$ as a function of time for different resolutions with $`r_o=40m`$ for EL+ES in PG coordinates.

## VII Conclusions

We have demonstrated that, at least for the case of single black hole spacetimes in spherical symmetry, it is possible to obtain long-term stable and convergent numerical simulations using the standard ADM system of equations if the equations are adjusted by introducing terms involving the constraints. Results were presented for three choices of lapse and shift. In addition, we introduced the concept of quasi well-posedness, which appears to be useful in characterizing the properties of the system of evolution equations. We are currently investigating the extension of this approach to three-dimensional evolutions.

###### Acknowledgements.

This work was supported by NSF PHY9800973, NSF PHY0090091 and by the Eberly research funds of the Pennsylvania State University. We wish to thank D. Arnold, A. Ashtekar, L. Lehner, R. Matzner, R. Price, O. Reula, H. Shinkai and B. Whiting for comments and helpful discussions.


text:
# References

BATALIN-TYUTIN QUANTISATION OF THE SPINNING PARTICLE MODEL
Subir Ghosh

Physics and Applied Mathematics Unit,

Indian Statistical Institute,

203 B. T. Road, Calcutta 700035,

India.

Abstract:
The Spinning Particle Model for anyon is analysed in the Batalin-Tyutin scheme of quantisation in extended phase space. Here additional degrees of freedom are introduced in the phase space such that all the constraints in the theory are rendered First Class that is commuting in the sense of Poisson Brackets. Thus the theory can be studied without introducing the Dirac Brackets which appear in the presence of non-commuting or Second Class constraints. In the present case the Dirac Brackets make the configuration space of the anyon non-canonical and also being dynamical variable dependent, poses problems for the quantisation programme. We show that previously obtained results, (e.g. gyromagnetic ratio of anyon being 2), are recovered in the Batalin-Tyutin variable independent sector in the extended space. The Batalin-Tyutin variable contributions are significant and are computable in a straightforward manner. The latter can be understood as manifestations of the non-commutative space-time in the enlarged phase space.

I: Introduction

Theoretical models for anyons , excitations of arbitrary spin and statistics in 2+1-dimensions, have been proposed from different perspectives, the most popular ones being the point charge- Chern-Simons gauge field interacting theory  and the later Symplectic Models (SM)  and the Spinning Particle Models (SPM) . The SM and SPM can be related to the anyon field equation proposed in . It has been questioned  whether the former interacting theory is a minimal description of the anyon. The latter models are free from complications of this type and both of them agree regarding results, (such as rigidity of the angular momentum and the gyromagnetic ratio $`(g)`$ being equal to $`2`$), for free anyons as well as anyons interacting with abelian external gauge field, respectively.

However, the SM or SPM’s are also interesting from another point of view. They have a close connection with the non-commutative space-time theories, which have created a lot of interest lately .

In SM , the symplectic structure has been posited, which gives rise to the non- commutativity in the particle configuration space. On the other hand, in SPM , the above mentioned feature is a result of the constraint structure of the model. The set of Second Class Constraints (SCC), with non-commuting Poisson Brackets (PB), (in the sense of Dirac classification scheme ), induces non-trivial Dirac Brackets (DB) among the coordinate variables. Since we will concentrate on SPM  in the present work, these comments will be elaborated.

Let us now impose a bit of caution on the methodology of evaluation of the $`g`$ value of anyon in SM and SPM. One compares the Hamiltonian (in the former ) and equations of motion (in the latter ) with the corresponding expressions obtained for a charged particle with spin in backgroung electromagnetic field in conventional (commuting) configuration space. The validity of this matching procedure can be questioned on the grounds that the nature of the configuration spaces in the two systems that are being compared are qualitatively different. But more importantly, the quantisation procedure runs into trouble due to operator ordering ambiguities in the DB’s in SP and SPM. This complication arises since in theories with non-linear constraints, the DBs can involve dynamical variables which is true for the present case. When these DBs are elevated to commutation relations in the process of quantisation, the above mentioned feature creates problems.

Finally, (as has been pointed out before ), the induced non-vanishing DB algebra $`[A_a(𝐱,t),A_b(𝐲,t)]`$ between supposedly ”external” electromagnetic field variables makes us wonder as to how far the treatment of the gauge fields being non-dynamical is justified. The problems related to the non-canonical nature of the coordinate system have also been discussed in .

Now we come to the motivation of the present work. All of the above problems can be naturally addressed in the Batalin-Tyutin (BT) quantisation scheme , (which is a particular type of the more general quantisation method ), where additional phase space degrees of freedom, (BT variables), are introduced in such a way that existing SCC’s are converted into (abelian) First Class Constraints (FCC, see ), thereby enhancing the gauge symmetry of the extended system. Thus the problematic DB’s can be avoided. One can resort to Dirac quantisation of FCC system by requiring that the physical states obey the FCCs, (i.e. $`FCCphysicalstate=0`$). Or else one can work in a convenient gauge where the problematic DBs do not appear. Later on we will use the first alternative.

Apart from the above reasons, we will find that by its own right, the problem of BT extension of SPM has several interesting features, some of them not addressed in the literature. Firstly, the SCC system in SPM can be expressed  in a covariant (and neat) form which however is reducible that is the SCCs are not independent. The BT formulation has been developed for irreducible or independent set of SCCs only. This forces us to choose an independent set of SCC thereby losing manifest covariance. As un unsolved problem, it will be interesting to study the BT extension of a generic reducible SCC system.

Let us put the present work in its proper perspective. We have developed a framework for treating the Spinning Particle Model as an FCC system in a BT-extended phase space where the background space-time is normal, (that is commuting), so that comparison between results obtained for anyons and analogous models in conventional space-time do not pose any ambiguity. (In particular, we will recover the result that for anyons $`g=2`$ but will indicate the presence of correction terms as well.) Also the canonical quantisation programme is rigorous since DBs do not appear and so operator ordering problems are absent in the canonical commutation relations.

Lastly we note that the present BT extension is classical in the sense that operator ordering problems in the extended system have not been addressed .

The paper is organised as follows: in section II we briefly reproduce the skeleton of BT formalism. This will also help us to fix the notations. We quickly introduce the SPM in section III for completeness. Section IV is devoted to developing the BT extension of SPM. This constitutes the main body of our work. Section V deals with the quantisation procedure and recovering of previous results. In section VI we conclude with a discussion and some future lines of work.

II: Batalin-Tyutin Formulation

In this section we state the main results of BT prescription  to be used. The basic idea behind the scheme is to introduce additional phase space variables (BT variables) $`\varphi _a^\alpha `$, besides the existing physical degrees of freedom $`(q,p)`$, such that all the constraints in the extended system are reduced to FCCs. This means that one has to modify the original constraints and Hamiltonian accordingly by putting BT-extension terms in them. The way to achieve this at the classical level has been provied in . Let us consider a set of constraints $`(\mathrm{\Theta }_\alpha ^a,\mathrm{\Psi }_ł)`$ and a Hamiltonian operator $`H`$ with the following Poisson Bracket (PB) relations,

$$\{\mathrm{\Theta }_\alpha ^a(q,p),\mathrm{\Theta }_\beta ^b(q,p)\}\mathrm{\Delta }_{\alpha \beta }^{ab}(q,p)0;\{\mathrm{\Theta }_\alpha ^a(q,p),\mathrm{\Psi }_\beta ^l(q,p)\}0$$

$$\{\mathrm{\Psi }_l(q,p),\mathrm{\Psi }_n(q,p)\}0;\{\mathrm{\Psi }_l(q,p),H(q,p)\}0.$$
(1)
In the above $`(q,p)`$ are referred to as physical variables and ”$``$” means that the equality holds on the constraint surface and $`\{p_a,x_b\}=g_{ab};g_{ab}=diag.(1,1,1)`$. Clearly $`\mathrm{\Theta }_\alpha ^a`$ and $`\mathrm{\Psi }_l`$ are SCC and FCC  respectively. The latter are responsible for gauge invariance and the former can be used as operator identities provided one uses the DBs  as defined below,

$$\{A(q,p),B(q,p)\}_{DB}=\{A,B\}\{A,\mathrm{\Theta }_\alpha ^a\}\mathrm{\Delta }_{ab}^{\alpha \beta }\{\mathrm{\Theta }_\beta ^b,B\},\mathrm{\Delta }_{\alpha \beta }^{ab}\mathrm{\Delta }_{bc}^{\beta \gamma }=\delta _c^a\delta _\gamma ^\alpha .$$
(2)
However, in systems with non-linear SCCs, in general the DBs can become dynamical variable dependent  due to the $`\{A,\mathrm{\Theta }_\alpha ^a\}`$ and $`\mathrm{\Delta }_{ab}^{\alpha \beta }`$ terms, leading to problems for the quantisation programme. To cure this type of pathology, BT formalism is a systematic framework where one introduces the BT variables $`\varphi _a^\alpha `$, obeying

$$\{\varphi _a^\alpha ,\varphi _b^\beta \}=\omega _{ab}^{\alpha \beta }=\omega _{ba}^{\beta \alpha },$$
(3)
where $`\omega _{ab}^{\alpha \beta }`$ is a constant (or at most a c-number funtion) matrix, with the aim of modifying the SCC $`\mathrm{\Theta }_\alpha ^a(q,p)`$ to $`\stackrel{~}{\mathrm{\Theta }}_\alpha ^a(q,p,\varphi _a^\alpha )`$ such that,

$$\{\stackrel{~}{\mathrm{\Theta }}_\alpha ^a(q,p,\varphi ),\stackrel{~}{\mathrm{\Theta }}_\beta ^b(q,p,\varphi )\}=0;\stackrel{~}{\mathrm{\Theta }}_\alpha ^a(q,p,\varphi )=\mathrm{\Theta }_\alpha ^a(q,p)+\mathrm{\Sigma }_{n=1}^{\mathrm{}}\stackrel{~}{\mathrm{\Theta }}_\alpha ^{a(n)}(q,p,\varphi );\stackrel{~}{\mathrm{\Theta }}^{a(n)}O(\varphi ^n)$$
(4)
The explicit terms in the above expansion are ,

$$\stackrel{~}{\mathrm{\Theta }}_\alpha ^{a(1)}=X_{\alpha \beta }^{ab}\varphi _b^\beta ;\mathrm{\Delta }_{\alpha \beta }^{ab}+X_{\alpha \gamma }^{ac}\omega _{cd}^{\gamma \delta }X_{\beta \delta }^{bd}=0$$
(5)
$$\stackrel{~}{\mathrm{\Theta }}_\alpha ^{a(n+1)}=\frac{1}{n+2}\varphi _d^\delta \omega _{\delta \gamma }^{dc}X_{cb}^{\gamma \beta }B_{\beta \alpha }^{ba(n)};n1$$
(6)
$$B_{\beta \alpha }^{ba(1)}=\{\stackrel{~}{\theta }_\beta ^{b(0)},\stackrel{~}{\theta }_\alpha ^{a(1)}\}_{(q,p)}\{\stackrel{~}{\theta }_\alpha ^{a(0)},\stackrel{~}{\theta }_\beta ^{b(1)}\}_{(q,p)}$$
(7)
$$B_{\beta \alpha }^{ba(n)}=\mathrm{\Sigma }_{m=0}^n\{\stackrel{~}{\theta }_\beta ^{b(nm)},\stackrel{~}{\theta }_\alpha ^{a(m)}\}_{(q,p)}+\mathrm{\Sigma }_{m=0}^n\{\stackrel{~}{\theta }_\beta ^{b(nm)},\stackrel{~}{\theta }_\alpha ^{a(m+2)}\}_{(\varphi )};n2$$
(8)
In the above, we have defined,

$$X_{\alpha \beta }^{ab}X_{bc}^{\beta \gamma }=\omega _{\alpha \beta }^{ab}\omega _{bc}^{\beta \gamma }=\delta _\alpha ^\gamma \delta _c^a.$$
(9)
A very useful idea is to introduce the Improved Variable $`\stackrel{~}{f}(q,p)`$  corresponding to each $`f(q,p)`$,

$$\stackrel{~}{f}(q,p,\varphi )f(\stackrel{~}{q},\stackrel{~}{p})=f(q,p)+\mathrm{\Sigma }_{n=1}^{\mathrm{}}\stackrel{~}{f}(q,p,\varphi )^{(n)};\stackrel{~}{f}^{(1)}=\varphi _c^\beta \omega _{\beta \gamma }^{cb}X_{bd}^{\gamma \delta }\{\theta _\delta ^a,f\}_{(q,p)}$$
(10)
$$\stackrel{~}{f}^{(n+1)}=\frac{1}{n+1}\varphi _c^\beta \omega _{\beta \gamma }^{cb}X_{bd}^{\gamma \delta }G(f)_\delta ^{d(n)};n1$$
(11)
$$G(f)_\beta ^{b(n)}=\mathrm{\Sigma }_{m=0}^n\{\stackrel{~}{\theta }_\beta ^{b(nm)},\stackrel{~}{f}^{(m)}\}_{(q,p)}+\mathrm{\Sigma }_{m=0}^{(n2)}\{\stackrel{~}{\theta }_\beta ^{b(nm)},\stackrel{~}{f}^{(m+2)}\}_{(\varphi )}+\{\stackrel{~}{\theta }_\beta ^{b(n+1)},\stackrel{~}{f}^{(1)}\}_{(\varphi )}$$
(12)
which have the property $`\{\stackrel{~}{\mathrm{\Theta }}_\alpha ^a(q,p,\varphi ),\stackrel{~}{f}(q,p,\varphi )\}=0`$. It can be proved that extensions of the original FCC $`\mathrm{\Psi }_l`$ and Hamiltonian $`H`$ are simply,

$$\stackrel{~}{\mathrm{\Psi }}_l=\mathrm{\Psi }(\stackrel{~}{q},\stackrel{~}{p});\stackrel{~}{H}=H(\stackrel{~}{q},\stackrel{~}{p}).$$
(13)
One can also reexpress the converted SCCs as $`\stackrel{~}{\mathrm{\Theta }}_\alpha ^a\mathrm{\Theta }_\alpha ^a(\stackrel{~}{q},\stackrel{~}{p})`$. The following identification theorem,

$$\{\stackrel{~}{A},\stackrel{~}{B}\}=\stackrel{~}{\{A,B\}_{DB}};\{\stackrel{~}{A},\stackrel{~}{B}\}_{\varphi =0}=\{A,B\}_{DB};\stackrel{~}{0}=0,$$
(14)
will play a crucial role in our later application. Hence the outcome of the BT extension is the closed system of FCCs with the FC Hamiltonian given below,

$$\{\stackrel{~}{\mathrm{\Theta }}_\alpha ^a,\stackrel{~}{\mathrm{\Theta }}_\beta ^b\}=\{\stackrel{~}{\mathrm{\Theta }}_\alpha ^a,\stackrel{~}{\mathrm{\Psi }}_l\}=\{\stackrel{~}{\mathrm{\Theta }}_\alpha ^a,\stackrel{~}{H}\}=0;\{\stackrel{~}{\mathrm{\Psi }}_l,\stackrel{~}{\mathrm{\Psi }}_n\}0;\{\stackrel{~}{\mathrm{\Psi }}_l,\stackrel{~}{H}\}0.$$
(15)
We will see that due to the non-linearity in the SCCs, the extensions in the improved variables, (and subsequently in the FCCs and FC Hamiltonian), turn out to be infinite series. This type of situation has been encountered before .

III: Spinning Particle Model revisited

The SPM proposed by us  where an anyon interacts with a $`U(1)`$ gauge field in 2+1-dimensions is given by the Lagrangian,

$$L=(m^2U^aU_a+\frac{j^2}{2}\sigma ^{ab}\sigma _{ab}+mjϵ^{abc}U_a\sigma _{bc})^{\frac{1}{2}}+eU_aA^a,$$
(16)
where

$$U^a=\frac{dx^a}{d\tau };\sigma ^{ab}=\frac{1}{2}ϵ^{abc}\sigma _c=\mathrm{\Lambda }_c^a\frac{d\mathrm{\Lambda }^{cb}}{d\tau }:\mathrm{\Lambda }_c^a\mathrm{\Lambda }^{cb}=\mathrm{\Lambda }_c^a\mathrm{\Lambda }^{bc}=g^{ab}.$$

Here $`(x^a,\mathrm{\Lambda }^{ab})`$ is a Poincare group element and also a set of dynamical variables of the theory. The canonical momenta are defined in the following way ,

$$p^a=\frac{L}{U_a}\pi ^aeA^a,S^{ab}=\frac{L}{\sigma _{ab}}\frac{1}{2}ϵ^{abc}S_c.$$
(17)
The phase space algebra is,

$$\{x_a,x_b\}=0;\{p_a,x_b\}=g_{ab};\{\pi _a,\pi _b\}=eF_{ab}=e(_aA_b_bA_a),$$
(18)
$$\{S^a,S^b\}=ϵ^{abc}S_c;\{S^a,\mathrm{\Lambda }^{0b}\}=ϵ^{abc}\mathrm{\Lambda }_c^0;\{\mathrm{\Lambda }^{0a},\mathrm{\Lambda }^{0b}\}=0.$$
(19)
In the free theory, one encounters the following set of FCC and SCC respectively,

$$\mathrm{\Psi }_1\pi _a\pi ^am^20;\mathrm{\Psi }_2\pi .S\frac{mj}{2}0,$$

$$\mathrm{\Theta }_1^a\mathrm{\Lambda }^{0a}\frac{p^a}{m}0;\theta _2^aS^{ab}p_b=ϵ^{abc}p_bS_c0.$$

In the free theory, $`\mathrm{\Psi }_1`$ and $`\mathrm{\Psi }_2`$ are respectivly the mass-shell condition and the Pauli-Lubanski relation.

In the interacting theory, One obtains to $`O(e)`$ the following set of FCCs and SCCs ,

$$\mathrm{\Psi }_1\pi _a\pi ^am^2+\frac{2eF_{ab}\pi ^a}{\pi ^2}\mathrm{\Theta }_2^b0;\mathrm{\Psi }_2\pi .S\frac{mj}{2}+\frac{ejF_{ab}\pi ^a}{2m\pi ^2}\mathrm{\Theta }_2^b0,$$
(20)
$$\mathrm{\Theta }_1^a\mathrm{\Lambda }^{0a}\frac{\pi ^a}{m}0;\theta _2^aS^{ab}\pi _b=ϵ^{abc}\pi _bS_c0.$$
(21)
Actually one can check that $`S^aS_a=\frac{j^2}{4}`$ but this is not independent of $`\mathrm{\Psi }_l`$ given above (20). To verify the constraint algebra after computing the PBs, one has to invoke the relation $`S^a=\frac{j}{2m}\pi ^a`$, which is consistent with the SCC $`\mathrm{\Theta }_2^a`$ and the normalisation agrees with the free theory. The above relation is valid also at the level of DBs .

Note that in the above set of SCC, $`\mathrm{\Theta }_1^a`$ has been imposed from outside such that the angular coordinates are properly restricted. It is easy to ascertain that the SCCs $`\mathrm{\Theta }_\alpha ^a`$ form a reducible set, since

$$\pi ^a\mathrm{\Theta }_{1a}=\frac{1}{2}m\mathrm{\Theta }^{1a}\mathrm{\Theta }_{1a};\pi ^a\mathrm{\Theta }_{2a}=0.$$

(With $`e=0`$, these features are true for the free theory stated above.) However, one works with this reducible set because the manifestly covariant structure simplifies calculations and one can invert  the SCC algebra matrix perturbativly to get the DBs. To $`O(e)`$ we obtain the DB between two generic field as ,

$$\{A,B\}_{DB}=\{A,B\}+\frac{eF_{ab}}{m^2}\{A,\mathrm{\Theta }_2^a\}\{\mathrm{\Theta }_2^b,B\}$$

$$\frac{1}{2}ϵ^{abc}S_c\{A,\mathrm{\Theta }_{1a}\}\{\mathrm{\Theta }_{1b},B\}+\frac{1}{2m}\{A,\mathrm{\Theta }_2^a\}\{\mathrm{\Theta }_{1a},B\}$$

$$\frac{1}{2m}\{A,\mathrm{\Theta }_1^a\}\{\mathrm{\Theta }_{2a},B\}.$$
(22)
In particular the following DB

$$\{x^a,x^b\}_{DB}=\frac{1}{2m^2}ϵ^{abc}S_c+O(e),$$
(23)
gives rise to the non-commutative space-time, which in turn generates the arbitrary spin contribution in the angular momentum. Note that the above algebra (23) is non-trivial in the free theory as well. Fixing the gauge condition $`x_0=\tau `$ as the proper time and using $`\mathrm{\Psi }_1`$ one ends up with the Hamiltonian,

$$H=(m^2\pi ^i\pi _i)^{\frac{1}{2}}eA_0,$$
(24)
where the $`\mathrm{\Theta }_2^a`$ dependent term has been dropped since we used DBs . To $`O(e)`$ the 2+1-dimensional analogue of the 3+1-dimensional Bargmann-Michel- Telegdi equation  is,

$$\dot{S}^a=\{H,S^a\}_{DB}=\frac{e}{m}F^{ab}S_b.$$
(25)
Comparison with the original equation  reveals that $`g=2`$ for the particle.

As we have mentioned before, the comparison has been carried out between the present system in non-commutative configuration space and the original Bargmann-Michel-Telegdi equation in normal coordinate system.

The connection between our result and that of SM  is very direct but subtle. Remember that in SM , a modified expression, (with an $`O(eF)`$ term), for the ”mass-shell” condition is used. (In our analysis, an analogous $`\mathrm{\Theta }_2^a`$-dependent term in (24) was ”strongly” put to zero since it was proportional to the SCCs and we use DBs.) The analysis in SM is for FCC system, according to Dirac, where one demands that the FCCs annihilate the physical states. (In our BT extended theory, we will also follow this route.) Now in , one solves perturbativly, for small $`e`$ the non-canonical symplectic algebra, (which is equivalent to our DBs), in terms of a set of canonical phase space variables and rewrites the Hamiltonian in terms of these new variables and finally compares this expression with the Hamiltonian of a charged particle in ordinary space time. Again the $`g=2`$ result  is reproduced. Obviously this derivation  also suffers from the same conceptual drawback as the previous one . With this background, we now move on to the BT extension of SPM.

IV:Batalin-Tyutin extension of Spinning Particle Model

This section comprises of the main body of our work where we introduce the BT machinary  in order to take into account the SCCs of the theory but at the same time avoid using a non-canonical coordinates. As mentioned before, we now use a smaller set of SCCs which are irreducible. Also note that in the extended space the constraints will be modified and we can no longer use the original constraints to simplfy the SCC algebra. Choosing $`\mathrm{\Theta }_\alpha ^i`$ as the irreducible set of SCCs, the PB algebra is,

$$\{\mathrm{\Theta }_\alpha ^i,\mathrm{\Theta }_\beta ^j\}=\mathrm{\Delta }_{\alpha \beta }^{ij}$$
(26)
$$\mathrm{\Delta }_{11}^{ij}=\frac{eF^{ij}}{m^2}$$
(27)
$$\mathrm{\Delta }_{12}^{ij}=[(\pi .\mathrm{\Lambda })g^{ij}\pi ^i\mathrm{\Lambda }^{0j}+\frac{e}{m}ϵ^{jbc}F_b^iS_c]=\mathrm{\Delta }_{21}^{ji};\pi .\mathrm{\Lambda }=\pi _a\mathrm{\Lambda }^{0a}$$
(28)
$$\mathrm{\Delta }_{22}^{ij}=[\pi _0(\pi .S)ϵ^{ij}+e(F^{ij}S^2+S_bF^{bi}S^j+S^iF^{jb}S_b];S^2=S_aS^a.$$
(29)
Next, following (5), we need to compute $`X_{\alpha \beta }^{ij}`$ as the whole calculational scheme rests on this quantity and its inverse. This we do perturbativly by first considering the free theory, (i.e. $`e=0`$), where,

$$\mathrm{\Theta }_1^i\mathrm{\Lambda }^{0i}\frac{p^i}{m};\mathrm{\Theta }_2^iϵ^{ibc}p_bS_c$$
(30)
$$\{\mathrm{\Theta }_\alpha ^i,\mathrm{\Theta }_\beta ^j\}=\mathrm{\Delta }_{\alpha \beta }^{ij}=\left(\begin{array}{cc}0& ((p.\mathrm{\Lambda })g^{ij}p^i\mathrm{\Lambda }^{0j})\\ ((p.\mathrm{\Lambda })g^{ij}p^j\mathrm{\Lambda }^{0i})& p_0(p.S)ϵ^{ij}\end{array}\right).$$
(31)
For the free theory we propose,

$$X_{\alpha \beta }^{ij}_{e=0}x_{\alpha \beta }^{ij}=\left(\begin{array}{cc}0& ((p.\mathrm{\Lambda })g^{ij}p^i\mathrm{\Lambda }^{0j})\\ g^{ij}& \frac{1}{2}p_0(p.S)ϵ^{ij}\end{array}\right).$$
(32)
One can check that $`x_{\alpha \beta }^{ij}`$ satisfies (5) for the free theory, provided we choose $`\omega _{ab}^{\alpha \beta }=ϵ^{\alpha \beta }g_{ab},ϵ^{12}=1`$. Since there are some artitrariness involved in $`x_{\alpha \beta }^{ij}`$ and $`\omega _{ab}^{\alpha \beta }`$, their choices are dictated by convenience. The inverse is defined below,

$$x_{ij}^{\alpha \beta }x_{\beta \gamma }^{jk}=\delta _\gamma ^\alpha g_k^i,$$

$$x_{jk}^{\gamma \delta }=\left(\begin{array}{cc}\frac{p_0(p.S)}{2(p.\mathrm{\Lambda })}(ϵ_{jk}\frac{ϵ_{jl}p_l\mathrm{\Lambda }_{0k}}{p_0\mathrm{\Lambda }_{00}})& g_{jk}\\ \frac{1}{(p.\mathrm{\Lambda })}(g_{jk}+\frac{p_j\mathrm{\Lambda }_{0k}}{p_0\mathrm{\Lambda }_{00}})& 0\end{array}\right).$$
(33)
This is an exact result. Now, for the interacting theory, we find to $`O(e)`$,

$$X_{\alpha \beta }^{ij}=x_{\alpha \beta }^{ij}(p_a\pi _a)+ey_{\alpha \beta }^{ij}$$
(34)
$$y_{11}^{ij}=\frac{F^{il}\mathrm{\Lambda }_{0l}p^j}{m^2p_0\mathrm{\Lambda }_{00}(p^k\mathrm{\Lambda }_{0k})};y_{22}^{ij}=\frac{1}{2}(S_0)^2F^{ij}$$
(35)
$$y_{21}^{ij}=\frac{S_0ϵ^{rl}F^{0r}S^l}{p_0(p.S)}g^{ij}$$
(36)
$$y_{12}^{ij}=\frac{1}{m}F_b^iϵ^{jbc}S_c+\frac{F^{il}\mathrm{\Lambda }_{0l}p_kϵ^{jk}p_0(p.S)}{2m^2p_0\mathrm{\Lambda }_{00}(p^k\mathrm{\Lambda }_{0k}}+\frac{S_0ϵ^{rl}F^{0r}S^l}{p_0(p.S)}((p.\mathrm{\Lambda })g^{ij}p^i\mathrm{\Lambda }^{0j})$$
(37)
The inverse, to $`O(e)`$, is of the form,

$$X_{ij}^{\alpha \beta }=x_{ij}^{\alpha \beta }(p_a\pi _a)ex_{ik}^{\alpha \mu }(p)y_{\mu \nu }^{kl}(p)x_{lj}^{\nu \beta }(p).$$
(38)
In the present work the explicit form of $`X_{ij}^{\alpha \beta }`$ will not be utilised. From the relation,

$$\stackrel{~}{\mathrm{\Theta }}_\alpha ^{i(1)}=X_{\alpha \beta }^{ij}\varphi _j^\beta ,$$
(39)
the explicit expressions for $`\stackrel{~}{\mathrm{\Theta }}_\alpha ^{i(1)}`$ (i.e. one BT variable extension term) are,

$$\stackrel{~}{\mathrm{\Theta }}_1^{i(1)}=(g^{ik}(\pi .\mathrm{\Lambda })\pi ^i\mathrm{\Lambda }^{0k})\varphi _k^2+e(y_{12}^{ik}\varphi _k^2+y_{11}^{ik}\varphi _k^1),$$
(40)
$$\stackrel{~}{\mathrm{\Theta }}_2^{i(1)}=\frac{1}{2}\pi _0(\pi .S)ϵ^{ik}\varphi _k^2+\varphi ^{1i}+e(y_{21}^{ik}\varphi _k^1+y_{22}^{ik}\varphi _k^2),$$
(41)
We emphasise that the series for $`\stackrel{~}{\mathrm{\Theta }}_1^i`$ and $`\stackrel{~}{\theta }_2^i`$ do not terminate and the higher order terms in $`\varphi _i^\alpha `$ can be derived by a straightforward but extremely tedious calculation. However, to check whether our target of converting SCCs to FCCs has really been achieved (to $`O(e)`$), one can convince oneself that,

$$\{(\mathrm{\Theta }_\alpha ^i+\stackrel{~}{\mathrm{\Theta }}_\alpha ^{i(1)}),(\mathrm{\Theta }_\beta ^j+\stackrel{~}{\mathrm{\Theta }}_\beta ^{j(1)})\}=0+O(\varphi ).$$

To check the cancellation of $`O(\varphi )`$ terms, $`\mathrm{\Theta }_\alpha ^{i(2)}`$ terms are required.

Now we move on to the $`oneBTextension`$ of the physical degrees of freedom, i.e the Improved Variables. The BT extensions of the BT variables themselves will vanish. Let us start by constructing the extensions for $`\pi ^a`$. From the generic expression given in (10), we have,

$$\stackrel{~}{\pi }^{a(1)}=\varphi _c^\beta \omega _{\beta \gamma }^{cb}X_{bd}^{\gamma \delta }\{\mathrm{\Theta }_\delta ^d,\pi ^a\}_{(q,p)}$$

$$=\frac{e}{m}(\varphi ^{2j}x_{ji}^{11}\varphi ^{1j}x_{ji}^{21})F^{ia}e\varphi _i^2ϵ^{ide}S_eF_d^a$$

$$=\frac{e}{m}F^{ia}[(\varphi ^{2j}\frac{p_0(p.S)}{2(p.\mathrm{\Lambda })}(ϵ_{ji}\frac{ϵ_{jl}p_l\mathrm{\Lambda }_i^0}{p_0\mathrm{\Lambda }_{00}})$$

$$+\varphi ^{ij}\frac{1}{(p.\mathrm{\Lambda })}(g_{ij}+\frac{p_j\mathrm{\Lambda }_{0i}}{p_0\mathrm{\Lambda }_{00})})]e\varphi ^2_iϵ^{ide}S_eF_d^a.$$
(42)
In a straightforward manner, one can compute the rest of the Improved Variables, from the structures given below,

$$\stackrel{~}{x}^{a(1)}=\varphi _c^\beta \omega _{\beta \gamma }^{cb}X_{bd}^{\gamma \delta }\{\mathrm{\Theta }_\delta ^d,x^a\}_{(q,p)}$$

$$=\frac{g^{ai}}{m}(\varphi ^{2j}x_{ji}^{11}\varphi ^{1j}x_{ji}^{21})+ϵ^{aic}S_c(\varphi ^{2j}X_{ji}^{12}\varphi ^{1j}X_{ji}^{22}),$$
(43)
$$\stackrel{~}{S}^{a(1)}=\varphi _c^\beta \omega _{\beta \gamma }^{cb}X_{bd}^{\gamma \delta }\{\mathrm{\Theta }_\delta ^d,S^a\}_{(q,p)}$$

$$=ϵ^{aic}\mathrm{\Lambda }_{0c}(\varphi ^{1j}X_{ji}^{21}+\varphi ^{2j}X_{ji}^{11})(g^{ai}(\pi .S)S^i\pi ^a)(\varphi ^{2j}X_{ji}^{12}\varphi ^{1j}X_{ji}^{22})$$
(44)
$$\stackrel{~}{\mathrm{\Lambda }}^{0a(1)}=\varphi _c^\beta \omega _{\beta \gamma }^{cb}X_{bd}^{\gamma \delta }\{\mathrm{\Theta }_\delta ^d,\mathrm{\Lambda }^{0a}\}_{(q,p)}$$

$$=(g^{ai}(\pi .\mathrm{\Lambda })\mathrm{\Lambda }^{0i}\pi ^a)(\varphi ^{1j}X_{ji}^{22}\varphi ^{2j}X_{ji}^{12}).$$
(45)
These Improved Variables also comprise of infinite sequences of higher order $`\varphi _i^\alpha `$ terms. As a non-trivial consistency check, we have tested the validity of the assertion that, to $`O(\varphi )`$, $`\stackrel{~}{\mathrm{\Theta }}_\alpha ^a\mathrm{\Theta }_\alpha ^a(\stackrel{~}{q},\stackrel{~}{p})`$ holds. To examine the $`two\varphi `$-term, $`O(\varphi ^2)`$-terms in $`\stackrel{~}{\mathrm{\Theta }}_\alpha ^i`$ and $`\stackrel{~}{f}(q,p)`$ are required. In the next section we will make use of these results to redetermine the $`g`$-value for the anyon in the extended phase space. However, as we have stressed before, the main achievement is that a consistent framework has been provided wherein quantisation of the SPM is unambiguous and derivation of the previous results are more transparant.

V: Application - $`g`$ in extended space

Let us start by recovering the ”mass shell” condition $`\stackrel{~}{\mathrm{\Psi }}_1`$ in BT extended space, which is simply given by

$$\stackrel{~}{\mathrm{\Psi }}_1\stackrel{~}{\pi }_a\stackrel{~}{\pi }^am^2+\frac{2e\stackrel{~}{F}_{ab}\stackrel{~}{\pi }^a}{\stackrel{~}{\pi }^2}\stackrel{~}{\mathrm{\Theta }}_2^b0.$$
(46)
Remembering that all the Improved Variables are of the form

$$\stackrel{~}{𝒜}(q,p,\varphi )=𝒜(q,p)+O(\varphi )+..;\stackrel{~}{F}=F(\stackrel{~}{x})=F(x)+(F)\varphi +..,$$

and neglecting $`F`$ terms, it is clear that $`\varphi `$-independent polynomial expressions of relevant operators will remain intact and so up to non-$`\varphi `$ terms, the previous relations survive. Hence proceeding in the same way as before, in the extended space, one can compare the Schrodinger equation to derive $`g=2+O(\varphi )`$ for anyons. However, more work is to be done to ascertain whether the $`\varphi `$-terms can contribute to $`g`$ even in the non-relativistic limit $`p_0m>>𝐩`$ used in . Notice that the Pauli-Lubanski relation is modified to,

$$\stackrel{~}{\mathrm{\Psi }}_2\stackrel{~}{\pi }.\stackrel{~}{S}\frac{mj}{2}+\frac{ej\stackrel{~}{F}_{ab}\stackrel{~}{\pi }^a}{2m\stackrel{~}{\pi }^2}\stackrel{~}{\mathrm{\Theta }}_2^b0.$$

This discussion corresponds to the SM model .

Returning to our SPM model , we follow the discussion in Section III and introduce the gauge fixing condition for $`\stackrel{~}{\mathrm{\Psi }}_1`$ to be $`\stackrel{~}{x}_0\tau =0`$ and obtain the Hamiltonian as,

$$\stackrel{~}{H}=(m^2\stackrel{~}{\pi }^i\stackrel{~}{\pi }_i+\stackrel{~}{\mathrm{\Theta }}_2^aterm)^{\frac{1}{2}}e\stackrel{~}{A}_0.$$
(47)
We immedietly notice that the Hamiltonian in extended space is not a rational function of the phase space variables and also that the $`\mathrm{\Theta }_2^a`$ term can not be dropped. But the theorem in (14) comes to the rescue. The $`\mathrm{\Theta }_2^a`$ can contribute to the equations of motion of the generic form $`\dot{\stackrel{~}{𝒜}}=\{\stackrel{~}{H},\stackrel{~}{𝒜}\}`$ by terms proportional to $`\stackrel{~}{\mathrm{\Theta }}_2^a`$ only since since PBs in extended space with $`\stackrel{~}{\mathrm{\Theta }}_\alpha ^i`$ vanish. So long as we are concerned with computing PBs, the theorem (14) can be applied without going to the details of the explicit structure of the respective operators in the PB. Hence utilising (14), we can simply transform the BMT equation of motion for the spin variable , in terms of the Improved Variables and get,

$$\dot{\stackrel{~}{S}}^a=\frac{e}{m}\stackrel{~}{F}^{ab}\stackrel{~}{S}_b\frac{e}{m}F^{ab}\stackrel{~}{S}_b.$$
(48)
This is our cherished Bargmann-Michel-Telegdi equation for anyon in BT- extended space. Obviously in the the BT extended space, the expression for $`g`$ is modified to $`g=2+O(\varphi )`$ for anyons, with the additional terms coming from the BT- variable dependent terms.

To facilitate a comparison with the results obtained in the physical variable sector, (without introducing BT variables), we can rewrite the Improved Variables by their physical counterpart and the BT extension terms. Since the phase space PB algebra (18,19, 3) is trivially known, we can compute $`\dot{S}^a`$ from (51) in a straightforward way. Obviously it will be of the form,

$$\dot{\stackrel{~}{S}}^a\frac{e}{m}F^{ab}\stackrel{~}{S}_b+O(\varphi ).$$
(49)
In a restricted way, we can eliminate the BT variables in terms of the physical variables if we adopt the Dirac quantisation prescription for systems with FCCs only, where the physical sector of Hilbert space is required to satisfy $`FCCPhysicalState=0`$ . In the present case, this reduces to

$$\stackrel{~}{\mathrm{\Theta }}_\alpha ^iPh.St.=\stackrel{~}{\mathrm{\Psi }}_lPh.St.=0.$$
(50)
Since we are interested in substituting the BT variables in the expression of the Hamiltonian in (47), we can solve the constraints for $`e=0`$ as the terms involving BT variables are already of $`O(e)`$. This leads to the set of equations,

$$\stackrel{~}{\mathrm{\Theta }}_1^i_{e=0}=\mathrm{\Theta }_1^i_{e=0}(g^{ij}(p.\mathrm{\Lambda })p^i\mathrm{\Lambda }^{oj})\varphi _j^2=0,$$

$$\stackrel{~}{\mathrm{\Theta }}_2^i_{e=0}=\mathrm{\Theta }_2^i_{e=0}+\varphi ^{1i}+\frac{1}{2}m^2S_0ϵ^{ik}\varphi _k^2=0.$$
(51)
From the solutions of the above equations, we find,

$$\varphi ^{1i}=\mathrm{\Theta }_2^i\frac{1}{2}m^2S_0ϵ^{ij}\varphi _j^2,$$

$$\varphi ^{2i}=\frac{\mathrm{\Theta }_{1j}}{(p.\mathrm{\Lambda })}(g^{ij}+\frac{p^i\mathrm{\Lambda }^{0j}}{p_0\mathrm{\Lambda }_{00}})$$
(52)
Notice that the BT-variables in this restricted derivation, being proportional to $`\mathrm{\Theta }_\alpha ^i`$, are non-vanishing even in the free theory since in the extended sector $`\mathrm{\Theta }_\alpha ^i`$ are no longer the constraints. It should be emphasised that this is not the whole story since $`\stackrel{~}{\mathrm{\Theta }}_\alpha ^i`$ consists of an infinite number of terms. Also there are indications that higher order terms will survive even if, on top of our small $`e`$ restriction, we incorporate the non-relativistic limit, (that is $`p_0m>>𝐩`$, as has been done in ). On the other hand, putting these expressions in the BMT equation shows the leading term in $`g=2`$ survives but there can be non-trivial velocity dependent corrections. Putting the BT variables back in (49) we get the BMT equation completely in the physical sector.

VI: Conclusion

To conclude, in the present work we have formulated an extension of the Spinning Particle Model for anyon  along the lines of Batalin-Tyutin quantisation scheme . The reason behind working in the extended phase space lies in the presence of non-linear Second Class Constraints in the model, which induce a non-canonical structure in the particle coordinate. Specifically, the position coordinates become non-commuting which creates problem for the quantisation programme. Also this makes the comparison between results obtained here and in conventional space-time models, difficult. The gyromagnetic ratio of anyon was obtained to be 2  by invoking precisly this type of matching.

To avoid this type of non-trivial structure in the space-time, the Batalin-Tyutin formalism  is adopted where extra BT-variables are introduced in the phase space in such a way that all the constraints become First Class in the extended space and the problematic Dirac Brackets can be avoided. Hence the commuting space-time structure is kept intact. However, the price to pay is that the extensions of the constraints and relevant variables turn out to be infinite serieses, (even to lowest order in $`e`$, the electromagnetic coupling), with higher powers of BT variables. One has to be very careful in taking the non-relativistic limit  and a priory it is difficult to determine whether the serises will terminate or not.

In the present work we have computed explicitly the one-$`\varphi `$ extensions of all the constraints and degrees of freedom. From the nature of the extensions, it is clear the BT-variable independent relations remain the same as the original ones. Hence one might say that the BT terms are effects of the non-commuting space-time algebra. However, we note that these results are partial in the sense that higher BT-variable terms may also contribute to this order of accuracy. One of our immediate goals in this area is to compute explicitly effect of the BT variable terms in observable quantities, e.g. $`g`$.

From another point of view, this work is significant since it may provide a mapping between theories in non-commuting space-time on one hand, and theories in conventional configuration space with extra spin degrees of freedom. Subsequently, Batalin-Tyutin extension can be introduced and one can check if results of simple non-commutative models  are reproduced.

Another interesting problem of the formal kind is to develop the Batalin-Tyutin scheme for reducible Second Class Constraint systems, such as the Spinning Particle Model in manifestly covariant form. One has to be careful in introducing the BT variables since reducibility in the system will be reflected in the number of these degrees of freedom.