text:
# Vistas in numerical relativity

## 1 Introduction

The Laser Interferometric Gravitational-wave Observatories LIGO/VIRGO  is a broad band detector, targeting gravitational radiation from compact sources of a few solar masses . Notable sources are binary coalescence of neutron stars and black holes, as well as emissions from black hole-torus systems as recently proposed . Their gravitational wave emissions provide a record of these strongly interacting catastrophic events and could contain most of the total energy released. What is the first source that LIGO/VIRGO may detect? Our knowledge of the gravitational wave-forms may be the determining factor in answering this question. For this reason, numerical simulations of general relativity (numerical relativity ) are receiving much attention in efforts towards matched filtering in searches for binary coalescence involving neutron stars and black holes .

### 1.1 Some astrophysical problems for numerical relativity

There exists a broad spectrum of candidate sources of gravitational waves. Many sources have been proposed, for instance the coalescence of binaries of neutron stars and black holes , supernovae (see ), and rapidly spinning neutron stars .

The gravitational wave-forms produced by neutron star-neutron star coalescence in the inspiral phase are well-understood with post-Newtonian expansion technique (see  for a recent review). Black hole-black hole coalescence is very promising because of the expected larger amplitude signal. However, their event rate is highly uncertain . Their wave-forms in the merger phase, believed to be relevant for LIGO/VIRGO detection, is not well understood. This is left as a challenge for numerical relativity (see, e.g., ).

A new model was recently proposed for gravitational radiation from at torus around a black hole . One particular feature is that it is associated with long gamma-ray bursts, whose intrinsic durations are about 20s on average. These energetic events are observed at a rate of about 2 per day and, corrected for beaming, at a rate of about 1 per year within a distance of 100Mpc. An important problem is the nature of the torus’ mass-quadrupole moment, which may be determined by numerical relativity.

### 1.2 Numerical relativity: integration on thin ice

Long-time integrations for the purpose of matched filtering requires accurate conservation of the elliptic constraints representing energy and momentum conservation. This is a natural requirement, which becomes apparent in numerical experiments .

The intrinsic hyperbolic structure of wave-motion permits well-posed initial value problems, as follows from the energy method. In the continuum limit, conservation of energy and momentum is exact and drops out of the “energy balance sheet”. Consequently, well-posedness for general relativity reduces to the anticipated results for its hyperbolic structure.

The nonlinear nature of the Einstein equations tends to introduce numerically a departure from exact conservation of the elliptic constraints. And the initial value problem for non-hyperbolic equations is often ill-posed, as for instance the Laplace and backward heat equation. The familiar Hadamard counter example to well-posedness is given by the solution

$`u_n={\displaystyle \frac{1}{n}}\mathrm{sin}nx\mathrm{sinh}ny`$ (1)
of Laplace’s equation on the upper half-plane:

$`u_{xx}+u_{yy}=0(\mathrm{}<x<\mathrm{},y0)`$ (2)
subject to the Cauchy data

$`u(x,0)=0,u_y(x,0)={\displaystyle \frac{1}{n}}\mathrm{sin}nx.`$ (3)
The solution blows up in the face of large $`n`$, even though the initial data approach zero in the norm of continuously differentiable functions. Numerically, this tends to result in ill-posedness: rapidly growing errors, regardless of the accuracy of initial conditions. The backward heat equation serves to illustrate similar ill-posedness in the presence of a first-order time-derivative.

The observed instabilities associated with constraint violations  suggests a link between variables with time-derivatives in the energy-momentum constraints and ill-posedness. Well-posed numerical relativity requires these constraints to be evolved in a consistent manner. A key test is the ‘dynamical’ evolution of a Schwarzschild black hole, where the dynamics derives from a singularity avoiding foliation of spacetime (see, e.g., ). In the absence of a covariant separation of the hyperbolic-elliptic structure of general relativity, we shall in this lecture discuss a recent proposal for an advanced-retarded evolution of the Einstein equations for exact preservation of the constraints under dynamical evolution.

## 2 The Einstein equations, spacetime foliation and conservation laws

The Einstein equations describe the structure of a curved spacetime manifold $`M`$ with four-covariant metric $`g_{ab}`$ in response to a stress-energy tensor $`T_{ab}`$. They are

$`R_{ab}{\displaystyle \frac{1}{2}}g_{ab}R=8\pi T_{ab}`$ (4)
as an equation for the Ricci tensor $`R_{bd}=R_{bcd}^c`$ and its scalar curvature $`R=R_c^c`$. These are expressions in terms of the Riemann tensor $`R_{bcd}^a`$. The left hand-side is commonly denoted by the Einstein tensor $`G_{ab}`$, with the property that $`^aG_{ab}0`$ (the Bianchi identity) is consistent with energy-momentum conservation $`_aT^{ab}=0`$. We are at liberty to include a cosmological constant term $`\mathrm{\Lambda }g_{ab}`$ on the right hand-side of (4). The equations of motion (4) derive from the Hilbert action

$`𝒮={\displaystyle _M}\sqrt{g}Rd^4x.`$ (5)
Translation invariance is a symmetry in this action (5). By Noether’s theorem, this leads to four conservation laws of energy and momentum. The associated gauge group is the Lorentz group SO(3,1,$`R`$) of boosts and rotations. These conservation laws become explicit on spacelike hypersurfaces $`\mathrm{\Sigma }_t`$: $`t=`$const., where $`t`$ denotes a timelike coordinate. Combined, the surfaces $`\mathrm{\Sigma }_t`$ provide a foliation of spacetime.

### 2.1 Spacetime foliation in spacelike hypersurfaces $`\mathrm{\Sigma }_t`$

Hypersurfaces $`\mathrm{\Sigma }_t`$ of constant time come with two vectors:

$`𝒩_a=g_{ta},n_a=_at/\sqrt{_at^at},`$ (6)
where $`n_a`$ denotes the unit normal ($`n^2=1)`$ to $`\mathrm{\Sigma }_t`$. (The vector $`𝒩_a`$ is commonly denoted by $`t_a`$ .) Generally, the covariant vectors $`𝒩_a`$ and $`n_a`$ are independent. Marching from one hypersurface to the next brings along a variation $`dt`$, along with the covariant displacement

$`ds_a=𝒩_adt.`$ (7)
The displacement $`ds_{(a)}`$ expresses $`𝒩_a`$ as a “flow of time.” It can be expressed in terms of orthogonal projections on $`n_a`$ and $`\mathrm{\Sigma }_t`$, in terms of the lapse function $`N`$ and shift functions $`N_a`$,

$`𝒩^a=Nn^a+N^a.`$ (8)
Here $`N=𝒩_an^a`$ and $`N_a=h_a^b𝒩_b`$, expressed in the metric

$`h_{ab}=g_{ab}+n_an_b`$ (9)
as the orthogonal projection of $`g_{ab}`$ onto $`\mathrm{\Sigma }_t`$. Note that $`ds^2=𝒩^2dt^2=g_{tt}dt^2`$ as the square of (7), so that $`g_{tt}=N^2+N_cN^c`$. With $`n_a=(n_t,0,0,0)`$, it follows that

$`g_{ab}=\left(\begin{array}{cc}N^cN_cN^2& N_j\\ N_i& h_{ij}\end{array}\right),`$ (12)
where $`i,j`$ refer to the spatial coordinates $`x^i`$ of $`(t,x^i)`$. The lapse function satisfies $`\sqrt{g}=N\sqrt{h}`$. The four degrees of freedom in the five functions $`(N,N_a)`$ are algebraically equivalent to $`𝒩_a`$. An equivalent expression for the line-element, in so-called 3+1 form , is

$`ds^2=\alpha ^2dt^2+\gamma _{ij}(dx^i+\beta ^i)(dx^j+\beta ^j),`$ (13)
where $`\alpha =N`$ is referred to as the redshift factor and $`\gamma _{ij}\beta ^j=g_{it}`$.

### 2.2 Conservation of energy and momentum

Coordinate invariance introduces a certain degeneracy in the Einstein equations. There are no second-order time-derivatives of $`g_{ab}`$ in the components $`G_{nb}=G_{ab}n^a`$ of the Einstein tensor $`G_{ab}=R_{ab}(1/2)g_{ab}R`$. Consequently, the expression $`G_{nb}`$ forms entirely out of Cauchy data on $`\mathrm{\Sigma }_t`$ (data and their first time-derivatives). The embedding of $`\mathrm{\Sigma }_t`$ in four-dimensional spacetime is expressed in terms of the symmetric extrinsic curvature tensor $`K_{ab}`$. If $`\stackrel{~}{n}_b`$ denotes a unit tangent to a geodesic orthogonal to $`\mathrm{\Sigma }_t`$, then

$`K_{ab}=_a\stackrel{~}{n}_b={\displaystyle \frac{1}{2}}_nh_{ab}.`$ (14)
Thus, $`K_{ab}`$ represents a time-like derivative of the metric in $`\mathrm{\Sigma }_t`$, which is a velocity of $`h_{ab}`$. We use here the sign convention in ; $`K_{ab}`$ with opposite sign is commonly used in numerical relativity. We then have

$`\begin{array}{cc}\hfill D_bK_a^bD_aK& =8\pi T_{an},\hfill \\ \hfill ^{(3)}R+K^2K_{ab}K^{ab}& =16\pi T_{nn}\hfill \end{array}`$ (17)
a consequence of the projection $`{}_{}{}^{(3)}R_{bcd}^{a}`$ of the four-dimensional Riemann tensor $`R_{bcd}^a`$ on $`\mathrm{\Sigma }_t`$ (see  for detailed calculations). Here, $`K=K_c^c`$ denotes the trace of $`K`$. These expressions (17) are, respectively, the conservation laws of linear momentum and energy. These equations are elliptic in the spatial coordinates internal to $`\mathrm{\Sigma }_t`$.

## 3 Two marching methods for hyperbolic formulations

A practical frame-work for numerical relativity consists of marching data from one space-like hypersurface $`\mathrm{\Sigma }_t`$: $`t=`$const. to the next. The hypersurface $`\mathrm{\Sigma }_t`$ is generally dynamical. A slicing of spacetime either comes before or after a choice of dynamical variables. Chosen before, the variables live in $`\mathrm{\Sigma }_t`$ as projections of the underlying four-covariant metric. Chosen after, one continues to work with four-covariant metric parametrized over the hypersurfaces $`\mathrm{\Sigma }_t`$. There is no dictum for the order of these choices, but they do give manifestly different formulations.

### 3.1 Slice first: the Hamiltonian approach

In the Hamiltonian approach, we consider first a choice of foliation of spacetime in spacelike hypersurface $`\mathrm{\Sigma }_t`$ of constant coordinate time $`t`$. Their dynamics is described by the projected metric $`h_{ab}`$ with canonical momentum $`\pi _{ab}`$, satisfying the Hamiltonian equations

$`\dot{h}_{ab}={\displaystyle \frac{\delta H}{\delta \pi ^{ab}}},\dot{\pi }^{ab}={\displaystyle \frac{\delta H}{\delta h_{ab}}},`$ (18)
where the dot denotes the Lie derivative $`_t`$ with respect to the vector field $`t^a=𝒩^a`$ of the flow of time (8). This Lie derivative reduces to differentiation with respect to $`t`$ in our coordinate system $`(t,x^i)`$.

The Hamiltonian equations derive from the Hilbert action (5). An excellent presentation is given in Appendix E in , which is briefly summarized here.

The Lagrangian density $`=R\sqrt{g}`$ can be expressed as a sum of the three-curvature $`{}_{}{}^{(3)}R`$ of $`\mathrm{\Sigma }_t`$ and a quadratic form of the extrinsic curvature tensor $`K_{ab}=_an_b`$, given by $`=\sqrt{h}N[^{(3)}R+K_{ab}K^{ab}K^2].`$ We have  $`K_{ab}=L_nh_{ab}/2=[\dot{h}_{ab}D_aN_bD_bN_a]/(2N),`$ where $`D_a=h_a^b_a`$ denotes the derivative internal to $`\mathrm{\Sigma }_t`$. It follows that

$`\pi ^{ab}={\displaystyle \frac{}{\dot{h}}}=\sqrt{h}(K^{ab}Kh^{ab}).`$ (19)
Coordinate invariance of the Einstein equations leaves $`g_{ta}`$ and, hence, $`(N,N_a)`$ freely specifyable. Tracing back, we indeed find no first-order time derivatives of $`g_{tb}`$ in the Lagrangian density $``$, whereby the associated canonical momenta vanish. The Hamiltonian density associated with the dynamical degrees of freedom reads therefore $`=\pi ^{ab}\dot{h}_{ab}.`$ The variational derivative of $`H=_{\mathrm{\Sigma }_t}\sqrt{h}d^3x`$ with respect to these gauge functions obtains the conservation laws of energy and momentum (17). The variational derivative with respect to the dynamical variables $`(h_{ab},\pi ^{ab})`$ obtains the ADM formulation

$`\dot{h}_{ab}2D_{(a}N_{b)}=2h^{1/2}N(\pi _{ab}{\displaystyle \frac{1}{2}}h_{ab}\pi )`$ (20)
and

$`\begin{array}{ccc}\hfill \dot{\pi }^{ab}2\pi ^{c(a}D_cN^{b)}& =\hfill & Nh^{1/2}(^{(3)}R\frac{1}{2}^{(3)}Rh^{ab})\hfill \\ & & +h^{1/2}(D^aD^bNh^{ab}D^cD_cN)\hfill \\ & & +\frac{1}{2}Nh^{1/2}h^{ab}(\pi :\pi \frac{1}{2}\pi ^2)\hfill \\ & & 2Nh^{1/2}(\pi ^{ac}\pi _c^b\frac{1}{2}\pi \pi ^{ab})\hfill \\ & & +h^{1/2}D_c(h^{1/2}N^c\pi ^{ab}).\hfill \end{array}`$ (26)
In numerical relativity, these Hamiltonian evolution equations are often considered in terms of the non-canonical pair $`(\gamma _{ij},K_{ij})`$, with $`\gamma _{ij}`$ as in (13) and $`K_{ij}`$, where $`i,j`$ refer to the spatial components in $`(t,x^i)`$. Thus, (20) and (26) become, using the vacuum case of (17),

$`\begin{array}{cc}(t_\beta )\gamma _{ij}\hfill & =2\alpha K_{ij},\hfill \\ (_t_\beta )K_{ij}\hfill & =D_iD_j\alpha +\alpha [^{(3)}R_{ij}2K_{ij}^2+KK_{ij}],\hfill \end{array}`$ (29)
where we use the definition of $`K_{ij}`$ with opposite sign of (14), following the convention in this context.

### 3.2 Hyperbolic formulations in the Hamiltonian approach

Hyperbolic systems of equations in the non-canonical variables $`(\gamma _{ij},K_{ij})`$, or closely related variables, have been derived in various forms by several groups, notably ; see  for a comprehensive review. This approach typically comprises constraints on the lapse function . For a recent comparison study between formulations in Hamiltonian variables and related hyperbolic formulations, see . Different formulations display various degrees of numerical stability , which appears not to depend significantly on the degree of hyperbolicity .

### 3.3 Dynamical conservation of constraints

The evolution equations for general relativity may formally be modified, such that the energy and momentum constraints become a stable manifold of physical solutions. This has recently been considered in a linearized treatment  and in Ashtekar’s formulation . This holds some promise in providing a unified treatment of the dynamical and the elliptic parts of general relativity. Numerical results on accuracy and stability are inconclusive at present .

### 3.4 Slice last: the four-covariant approach

General relativity can be written as nonlinear wave equations for the Riemann-Cartan connections in the tetrad formulation. This builds on Pirani’s arguments concerning the role of the Riemann tensor in gravitational waves  and on Yang-Mills formulations of general relativity, following Utiyama  and developed by Ashtekar and co-workers . Starting point is a divergence equation for the Riemann tensor with an anti-symmetric derivative of the stress-energy tensor as a source-term.

The interwoveness of wave motion and causal structure distinguishes gravity from other field theories. This becomes apparent in nonlinear wave equations for the connections on the curved spacetime manifold side-by-side with equations of structure for the evolution of the metric in the tangent bundle.

The tetrad approach  bears some relation to but is different from Ashtekar’s propram on nonperturbative quantum gravity. The original Ashtekar variables are SU(2,C) soldering forms and an associated complex connection in which the constraint equations become polynomial. The Riemann-Cartan variable is a real SO(3,1,R) connection. In Ashtekar’s variables, a real spacetime is recovered from the complex one by reality constraints. See Barbero for a translation of Ashtekar’s approach into SO(3,R) phase space with real connections . The main innovation in  is the incorporation of the Lorentz gauge condition (36) which obtains new hyperbolic evolution equations in four-covariant form (below).

Following Pirani, we take the view that gravitational wave-motion is contained in the Riemann tensor, $`R_{abcd}`$. It satisfies the Bianchi identity

$`3_{[e}R_{ab]cd}=_eR_{abcd}+_aR_{bacd}+_bR_{eacd}=0.`$ (30)
This gives rise to the homogeneous divergence equation $`^aR_{abcd}=0,`$ where $`R_{abcd}=(1/2)ϵ_{ab}^{ef}R_{efcd}`$ denotes its dual. Upon interaction with matter in accord with the Einstein equations, the Ricci tensor satisfies $`R_{ab}=8\pi [T_{ab}\frac{1}{2}g_{ab}T].`$ The Bianchi identity above also gives $`^dR_{abcd}=2_{[b}R_{a]c}`$, which obtains the inhomogeneous divergence equation

$`^aR_{abcd}=16\pi (_{[c}T_{d]b}{\displaystyle \frac{1}{2}}g_{b[d}_{c]}T).`$ (31)
The quantity on the right hand-side shall be referred to as $`16\pi \tau _{bcd}`$. This term is divergence free, $`^b\tau _{bcd}0`$, on account of the conservation law $`^aT_{ab}=0`$ and consistent with divergence-free condition $`^b^aR_{abcd}=0`$ on the left hand-side (31) (by anti-symmetry of the Riemann tensor in its first two indices).

Introduce the Riemann-Cartan connections $`\omega _{a\mu \nu }=(e_\mu )^c_a(e_\nu )_c`$ associated with a tetrad $`\{(e_\mu )^b\}_{\mu =1}^4`$. Then the above-mentioned homogeneous and inhomogeneous divergence equations take the form

$`\widehat{}^aR_{ab\mu \nu }=0,\widehat{}^aR_{ab\mu \nu }=16\tau _{b\mu \nu },`$ (32)
where the $`\omega _{a\mu \nu }`$ define a gauge-covariant derivative in accord with the Yang-Mills construction $`\widehat{}_a=_a+[\omega _a,].`$ The first of (32) gives rise to the representation $`R_{ab\mu \nu }=_a\omega _{b\mu \nu }_b\omega _{a\mu \nu }+[\omega _a,\omega _b]_{\mu \nu }.`$ The gauge covariant derivative satisfies the identity $`\widehat{}_a(e_\mu )^b0`$, which implies the equations of structure $`_{[a}(e_\mu )_{b])}=(e^\nu )_{[b}\omega _{a]\nu \mu }`$ – leaving $`_t(e_\mu )_t`$ undefined. Next, define $`\xi ^b=(_t)^b`$, and introduce the tetrad lapse functions

$`N_\mu =(e_\mu )_a\xi ^a`$ (33)
as freely specifiable functions. Thus, the equations of structure become a system of first-order ordinary differential equations

$`_t(e_\mu )_b+\omega _{t\mu }^\nu (e_\nu )_b=_bN_\mu +\omega _{b\mu }^\nu N_\nu .`$ (34)
The tetrad lapse functions are algebraically equivalent to the familiar Hamiltonian lapse, $`N`$, and shift functions, $`N_a`$, through (12):

$`g_{at}=N_\alpha (e^\alpha )_a=(N_qN^qN^2,N_p).`$ (35)
The term $`\omega _{b\mu \nu }N^\nu `$ on the right hand-side of (34) shows that the tetrad lapse functions introduce different transformations on each of the legs; the term $`\omega _{t\mu }^\nu (e_\nu )_b`$ on the left hand-side introduces a transformation with applies to all four legs simultaneously. It follows that it is the infinitesimal Lorentz transformations $`\omega _{t\mu }^{}`$ which provide the internal gauge transformations.

### 3.5 Hyperbolic equations in the four-covariant approach

An important aspect here is internal gauge fixing on the Lorentz group SO(3,1,$`R`$) associated with the choice of tetrad. To fix gauge, we propose using the Lorentz gauge

$`c_{\mu \nu }:=^a\omega _{a\mu \nu }=0.`$ (36)
This fixes unique evolution equations for the internal gauge, resulting in nonlinear wave equations for the connections $`\omega _{a\mu \nu }`$. These complement the equations of structure for the evolution of the tetrad legs, and together form a complete system of evolution equations. The Lorentz gauge (36) defines a choice of acceleration of the tetrad legs, through the infinitesimal Lorentz transformations $`\omega _{t\mu \nu }`$ mentioned above. In a different context of compact gauge groups and a metric with Euclidean signature, its geometric significance has been interpreted by Lewandowski et al. (1983) . These six constraints (36) can be given a hyperbolic implementation by application of the divergence technique

$`\widehat{}^a\{R_{ab\mu \nu }+g_{ab}c_{\mu \nu }\}=16\pi \tau _{b\mu \nu },`$ (37)
which preserves $`c_{\mu \nu }=0`$ is preserved in the future domain of dependence of the support of physical initial data. By explicit calculation, we have

$`\widehat{\mathrm{}}\omega _{a\mu \nu }R_a^c\omega _{c\mu \nu }[\omega ^c,_a\omega _c]_{\mu \nu }=16\pi \tau _{a\mu \nu },`$ (38)
where $`\widehat{\mathrm{}}`$ denotes the Yang-Mills wave operator $`\widehat{}^2`$. The Ricci tensor on the left hand-side may be understood in terms of $`T_{ab}`$ using the Einstein equations. The above provides the following covariant separation .
Theorem 1. Gravitational waves propagate on a curved spacetime manifold by nonlinear wave equations in a Lorentz gauge on the Riemann-Cartan connections. In response to this wave motion, the causal structure of the manifold evolves in the tangent bundle by the equations of structure. The Hamiltonian lapse and shift functions find their algebraic counterparts in the tetrad lapse functions $`N_\mu `$.

We remark that away from the matter source, the vacuum equations read

$`\widehat{\mathrm{}}\omega _{a\mu \nu }[\omega ^c,_a\omega _c]_{\mu \nu }=0.`$ (39)
This vacuum case has been considered numerically in one-dimensional Gowdy-tests  by using the underlying first-order system for the components of the Riemann tensor (37).

## 4 Hyperbolic Einstein-MHD equations

Relativistic hydrodynamics and magnetohydrodynamics has received considerable attention in the simulations of astrophysical jets . Recently, these efforts result in simulations of astrophysical jets around black holes  with extensions to flows around rotating black holes for a few dynamical time-scales . The latter shows a transition of accretion disk outflows towards a state of differential rotation in the vicinity of the black hole.

GRBs from rotating black holes are associated with a compact torus or disk, representing binary black hole-neutron star coalescence , failed-supernovae  or hypernovae . The torus or disk may well be magnetized with the remnant field of the progenitor star, i.e., the neutron star in the coalescence scenario or young massive star in the failed-supernova or hypernova scenario. The suggests simulating the creation of gravitational waves by high-density matter around a black hole in the approximation of ideal magnetohydrodynamics.

A hyperbolic formulation of ideal magnetohydrodynamics, used in the simulation of relativistic jets , is

$`\begin{array}{cc}\hfill _aT^{ab}& =0,\hfill \\ \hfill _a(u^{[a}h^{b]}+g^{ab}c)& =0,\hfill \\ \hfill _a(ru^a)& =0,\hfill \end{array}`$ (43)
expressing conservation of energy-momentum, Faraday’s equations in divergence form which conserves the constraint $`c=u^ch_c=0`$ , and conservation of baryon number. This is a single fluid description of an ideal, inviscid fluid with stress-energy tensor

$`T_{ab}=(r+\gamma P/(\gamma 1)+h^2)u^au^b+(P+h^2/2)g^{ab}h^ah^b,`$ (44)
where $`u^a`$ denotes the velocity four-vector, $`r`$ the comoving rest mass density and $`P`$ the pressure with polytropic equation of state $`P=Kr^\gamma `$ and polytropic index $`\gamma `$.

Perhaps Theorem 1 and (43) may serve as a starting point for hyperbolic Einstein-MHD equations for the purpose of numerical simulations.

## 5 Past and future gauge in numerical relativity

The Einstein equations pose six equations for dynamical evolution plus the four constraints of energy and momentum conservation. The latter involve the normal vector $`n^a`$ to the surfaces of foliation, i.e.: the projection operator onto these surfaces. The six dynamical degrees of freedom pertain to variables subject to second-order time-derivatives, while gauge-variables are subject to constraints on their first-order time-derivatives. In a discrete setting, the first live on three and the second on two hypersurfaces of constant time. This suggests to consider the ten degrees of freedom involved to be the six dynamical degrees of freedom supplemented with four gauge-variables from the past. The gain is exact conservation of energy-momentum, traded off against exact projections in the past.

Non-exact projections naturally permit an uncertainty between the three-metric and its canonical momentum within the underlying context of a four-covariant theory, i.e.: also in regards to the association with the hypersurface at hand. In the covariant approach of , this would thus reflect an uncertainty in the tetrad elements, which define the projection, and their connections. This points towards a potential connection to quantum gravity. Indeed, soon after this work was proposed , the author learned of a very interesting independent discussion on the problem of consistent discretizations in this context .

### 5.1 A discretized initial value problem

We illustrate our this approach on the vacuum Einstein equations

$`R_{ab}=0.`$ (45)
The Ricci tensor $`R_{ab}`$ is a second-order expression in the metric $`g_{ab}`$, whereby (45) defines a relationship between metric data $`(g_{ab}^{n1},g_{ab}^n,g_{ab}^{n+1})`$ on a triple of time-slices $`t_{n1}<t_n<t_{n+1}`$:

$`R_{ab}(g_{ab}^{n+1},g_{ab}^n,g_{ab}^{n1})=0.`$ (46)
Here, $`R_{bd}=R_{bcd}^a`$ derived from the Riemann tensor

$`R_{bcd}^a=_d\mathrm{\Gamma }_{bc}^a_c\mathrm{\Gamma }_{bd}^a+\mathrm{\Gamma }_{bc}^e\mathrm{\Gamma }_{ed}^a\mathrm{\Gamma }_{bd}^e\mathrm{\Gamma }_{ec}^a.`$ (47)
This expression (47) can be discretized by finite differencing on a triple of time-slices with preservation of the quasi-linear second-order structure of $`R_{ab}`$.

Algebraic gauge-fixing takes the form of specifying the components $`N_a=g_{ta}`$ in coordinates $`\{x^a\}_{a=1}^4`$ with $`t=x^1`$ time-like. A gauge-choice on a triple of time-slices amounts to a choice of $`(N_a^{n1},N_a^n,N_a^{n+1})`$. This gauge-choice in the metric arises explicitly in the Hamiltonian constraints of energy-momentum conservation. The components $`h_{ij}=g_{ij}`$, where $`i,j=2,3,4`$ refer to projections of the metric into the time-slice $`t=`$const., which describe the dynamical part of the metric. The combination $`(h_{ij},N_a)`$ reflects the mixed hyperbolic-elliptic structure in numerical relativity and (45) represents ten evolution equations in these variables on a triple of time-slices.

In algebraic gauge-fixing, we prescribe $`N_a^{n+1}`$ as a future gauge in computing $`h_{ij}^{n+1}`$ on a future hypersurface $`t=t_{n+1}`$ from data at present and past hypersurfaces $`t=t_{n1}`$ and $`t=t_n`$. Re-introducing $`N_a^{n1}`$ as dynamical gauge in the past gives closure, leaving $`h_{ij}^{n1}`$ fixed. This combination of ten degrees of freedom defines an advanced hyperbolic-retarded elliptic evolution of the metric. The paritioning of the metric in past and future variables as

$`g_{ab}=(h_{ij}^{n+1},N_a^{n1})=\left(\begin{array}{cccc}N_1^{n1}& N_2^{n1}& N_3^{n1}& N_4^{n1}\\ N_2^{n1}& h_{xx}^{n+1}& h_{xy}^{n+1}& h_{xz}^{n+1}\\ N_3^{n1}& h_{xy}^{n+1}& h_{yy}^{n+1}& h_{yz}^{n+1}\\ N_4^{n1}& h_{xz}^{n+1}& h_{zy}^{n+1}& h_{zz}^{n+1}\end{array}\right)`$ (52)
thus obtains ten dynamical variables in the ten equations

$`R_{ab}(h_{ij}^{n+1},N_a^{n1},\mathrm{})=0\text{at}t=t_n.`$ (53)
The dots refer to the remaining data $`(h_{ij}^{n1},h_{ij}^n,N_a^n,N_a^{n+1})`$, which are kept fixed while solving for $`(h_{ij}^{n+1},N_a^{n1})`$. Thus, (53) which takes into account all ten Einstein equations with no reduction of variables. Time-stepping by (53) evolves the metric into the future with dynamical gauge in the past, in an effort to satisfy energy-momentum conservation within the definition of the discretized Ricci tensor. Because (53) comprises derivatives of $`N_a`$ only to first-order in time, numerically though the data $`N_a^{n+1}`$ and $`N_a^{n1}`$, we anticipate that the evolution of $`N_a`$ is of first-order in the $`t`$discretization $`\mathrm{\Delta }t`$. This introduces non-exactness in $`h_{ij}^{n1}`$ as projections of $`g_{ab}`$ on $`t\mathrm{\Delta }t`$ to within the same order of accuracy. It may result in a first-order drift in the $`t`$labeling of the hypersurfaces – permitted by coordinate invariance.

### 5.2 A polarized Gowdy wave

The presented approach can be illustrated on a polarized Gowdy wave. Gowdy cosmologies have compact space-like hypersurfaces with two Killing vectors $`_\sigma `$ and $`_\delta `$. With cyclic boundary conditions, the space-like hypersurfaces are homeomorphic to the three-torus as a model universe collapsing into a singularity. The associated line-element is (see, e.g., )

$`ds^2=e^{(\tau \lambda )/2}\left(e^{2\tau }d\tau ^2+d\theta ^2\right)+d\mathrm{\Sigma }^2,`$ (54)
where $`\lambda =\lambda (\tau ,\theta )`$ and $`d\mathrm{\Sigma }`$ denotes the surface element in the space supported by the two Killing vectors. Polarized Gowdy waves form a special case, which permit a reduction to

$`d\mathrm{\Sigma }^2=e^\tau \left(e^Pd\sigma ^2+e^Pd\delta ^2\right).`$ (55)
Here $`P`$ satisfies a linear wave-equation $`P_{\tau \tau }=e^{2\tau }P_{\theta \theta }`$; a long wave-length solution is

$`P_0(\tau ,\theta )=Y_0(e^\tau )\mathrm{cos}\theta ,`$ (56)
where $`Y_0`$ is the Bessel function of the second kind of order zero. This leaves

$`\lambda (\tau ,\theta )={\displaystyle \frac{1}{2}}Y_0(e^\tau )Y_1(e^\tau )e^\tau \mathrm{cos}2\theta +{\displaystyle \frac{1}{2}}{\displaystyle _{e^\tau }^1}\left(Y_0^2(s)+Y_0^2(s)\right)s𝑑s.`$ (57)
A spectrally accurate numerical integration is described in .

The implicit equation (53) for the dynamical variables $`(h_{ij}^{n+1},N_a^{n1})`$ has been implemented numerically. We have done so by solving for the all ten components $`(h_{ij}^{n+1},N_a^{n1}`$) using Newton iterations on these variables. This procedure uses a numerical evaluation of the Jacobian

$`J_{AB}={\displaystyle \frac{R_A}{U_B}}`$ (58)
where the capital indices $`A,B=1,2,\mathrm{},10`$ refer to the labeling

$`\begin{array}{cc}\hfill R_A& =(R_{11},R_{22},R_{33},R_{44},R_{12},R_{13},R_{14},R_{23},R_{24},R_{34}),\hfill \\ \hfill U_B& =(N_1^{n1},h_{11}^{n+1},h_{22}^{n+1},h_{33}^{n+1},N_2^{n1},N_3^{n1},N_4^{n1},h_{23}^{n+1},h_{24}^{n+1},h_{34}^{n+1}).\hfill \end{array}`$ (61)
The Ricci tensor (47) has been implemented by second-order finite differencing, such that it remains quasi-linear in the second derivatives. In particular, the Christoffel symbols

$`\mathrm{\Gamma }_{ab}^c={\displaystyle \frac{1}{2}}g^{ce}(g_{cb,a}+g_{ac,b}g_{ab,c})`$ (62)
is obtained by symmetric finite-differencing on the metric components, and itself differentiated by the product rule following individual numerical differentiations of $`g^{ab}`$ and $`(g_{cb,a}+g_{ac,b}g_{ab,c})`$. The choice of future gauge $`N_a^{n+1}`$ is provided by the the components

$`g_{at}=(e^{(\tau \lambda )/2},0,0.0)`$ (63)
of the analytical line-element (54-57), which facilitates error analysis by direct comparison of the numerical results with the analytic solution. It will be appreciated that in principle other choices of $`N_a^{n+1}`$ can be made.

Fig. 1 shows numerical results for evolution of initial data on the interval $`0\tau 4`$. The results show that all Einstein equations in the form of $`R_{ab}=0`$ are satisfied with arbitrary precision, while the metric components are solved accurately to within one percent. The asymptotic behavior of the implicit corrections to the lapse functions are shown in Fig. 2. Note that these corrections are finite to first-order in $`\mathrm{\Delta }t`$: the corrections $`\delta N`$ on the past lapse satisfy

$`{\displaystyle \frac{\delta N}{N}}={\displaystyle \frac{N^{n+1}(\tau _n)N^{n1}(\tau _{n+2})}{N^n(\tau _{n+1})}}=O\left(\mathrm{\Delta }\tau \right).`$ (64)
This first-order dependence is a testimonial to the fact that the lapse function appears in the Einstein equations to first-order in time.
Theorem 2. A choice of gauge in the future and a dynamical gauge in the past obtains a discretization of general relativity consistent with energy and momentum conservation, permitting all ten Einstein equations to be solved in ten dynamical variables in the presence on non-exact projections of the four-covariant metric on the past hypersurface of constant time.

## 6 Summary and conclusion

Well-posed numerical relativity is a long-standing challenge in the calculation of wave-forms for astrophysical sources of gravitational radiation. A necessary condition for stable numerical relativity is accurate conservation of the energy and momentum constraints (“integration on thin ice”). This has been pursued by implementing these constraints dynamically . Here, we have explored a consistent discretization for exact conservation of energy-momentum constraints using a choice of gauge in the future and a dynamical gauge in the past. This permits integration of all ten Einstein equations, while allowing for in-exact projections of the four-covariant metric onto the surfaces of foliation of spacetime. The simulation of a nonlinear one-dimensional Gowdy wave by implicit time-stepping according to the ten discretized vacuum Einstein equations (53) serves to illustrate a numerical implementation.

A major open problem is obtaining sufficient conditions for stability. In this presented approach, it becomes of interest to consider novel definitions of the future gauge as a function of present gauge. We leave this as a suggested direction for future development. In light of the recent discussion by Gambini and Pullin (2002) , the question arises: is well-posed numerical relativity related to consistent discretization in quantum gravity?

Evolving eternal Schwarzschild black holes in 3+1 may serve as a test problem for these developments . More generally, it would be of interest to consider exact conservation of energy-momentum in higher dimensions, perhaps using a combination of any of the modern hyperbolic formulations and efficient elliptic solvers.

Acknowledgments. The author thanks the APCTP at Pohang University for hosting a very stimulating Winter School on black hole astrophysics and stimulating discussions with L. Wen. This research is supported by NASA Grant 5-7012 and an MIT C.E. Reed Award.

References