Evolutions of Gowdy, Brill and Teukolsky initial data on a smooth lattice

TL;DR

This paper demonstrates that a lattice method for computational general relativity can effectively evolve initial data for various space-times, matching results from existing methods and handling axisymmetric instabilities and boundary conditions well.

Contribution

The paper introduces a lattice-based approach for evolving initial data in general relativity, showing its effectiveness for Brill, Teukolsky, and Gowdy space-times compared to traditional methods.

Findings

Lattice method matches results from conventional techniques.

Effective handling of axisymmetric instabilities.

Good performance in wave passage through boundaries.

Abstract

Numerical results, based on a lattice method for computational general relativity, will be presented for Cauchy evolution of initial data for the Brill, Teukolsky and polarised Gowdy space-times. The simple objective of this paper is to demonstrate that the lattice method can, at least for these space-times, match results obtained from contemporary methods. Some of the issues addressed in this paper include the handling of axisymmetric instabilities (in the Brill space-time) and an implementation of a Sommerfeld radiation condition for the Brill and Teukolsky space-times. It will be shown that the lattice method performs particularly well in regard to the passage of the waves through the outer boundary. Questions concerning multiple black-holes, mesh refinement and long term stability will not be discussed here but may form the basis of future work.

Equations

$$d{s}^2=\left(g_{\alpha \beta }-\frac{1}{3}{R}_{\alpha \mu \beta \nu }{x}^{\mu }{x}^{\nu }-\frac{1}{6}{R}_{\alpha \mu \beta \nu ,\gamma }{x}^{\mu }{x}^{\nu }{x}^{\gamma }+\cdots \right)d{x}^{\alpha }d{x}^{\beta }$$

$$d{s}^2=\left(g_{\alpha \beta }-\frac{1}{3}{R}_{\alpha \mu \beta \nu }{x}^{\mu }{x}^{\nu }\right)d{x}^{\alpha }d{x}^{\beta }$$

$$L_{pq}^2=g_{\alpha \beta }\Delta x_{pq}^{\alpha }\Delta x_{pq}^{\beta }-\frac{1}{3}{R}_{\alpha \mu \beta \nu }{x}_p^{\alpha }{x}_p^{\beta }{x}_q^{\mu }{x}_q^{\nu }+\mathcal{O}\left(R{L}^5\right)$$

$$\frac{d{L}_{oq}}{dt}=[v_{\mu }(N{n}^{\mu }){]}_o^q$$

$$\frac{d{L}_{oq}}{dt}=-N{K}_{\alpha \beta }{v}_{oq}^{\alpha }{v}_{oq}^{\beta }{L}_{oq}+\mathcal{O}\left(L^2\right)$$

$$\frac{d{L}_{oq}}{dt}=-\frac{1}{2}\left({\left(N{K}_{\alpha \beta }\right)}_{q\overline{q}}{v}_{qo\overline{q}}^{\alpha }{v}_{qo\overline{q}}^{\beta }+{\left(N{K}_{\alpha \beta }\right)}_{o\overline{o}}{v}_{oq\overline{o}}^{\alpha }{v}_{oq\overline{o}}^{\beta }\right)L_{oq}+\mathcal{O}\left(L^2\right)$$

$$\frac{d{\tilde{L}}_{oq}}{dt}$$

$$\frac{d{x}_p^{\alpha }}{dt}$$

$$\frac{d{K}_{xx}}{dt}$$

$$\frac{d{K}_{yy}}{dt}$$

$$\frac{d{K}_{zz}}{dt}$$

$$\frac{d{K}_{xy}}{dt}$$

$$\frac{d{K}_{xz}}{dt}$$

$$\frac{d{K}_{yz}}{dt}$$

$$R_{xyxy},R_{xyxz},R_{xyyz},R_{xzxz},R_{xzyz},R_{yzyz}$$

$$R_{txxy},R_{tyxy},R_{tzxy},R_{txxz},R_{tyxz},R_{tzxz},R_{tyyz},R_{tzyz}$$

$$R_{txtx},R_{tyty},R_{tztz},R_{txty},R_{txtz},R_{tytz}$$

$$R_{xyxy,t}$$

$$R_{xyxz,t}$$

$$R_{xyyz,t}$$

$$R_{xzxz,t}$$

$$R_{xzyz,t}$$

$$R_{yzyz,t}$$

$$R_{txxy,t}$$

$$R_{tyxy,t}$$

$$R_{tzxy,t}$$

$$R_{txxz,t}$$

$$R_{tyxz,t}$$

$$R_{tzxz,t}$$

$$R_{tyyz,t}$$

$$R_{tzyz,t}$$

$$R_{\alpha \beta \mu \nu ,\gamma }=R_{\alpha \beta \mu \nu †\gamma }-m^{\lambda }{}_{\alpha \gamma }{R}_{\lambda \beta \mu \nu }-m^{\lambda }{}_{\beta \gamma }{R}_{\alpha \lambda \mu \nu }-m^{\lambda }{}_{\mu \gamma }{R}_{\alpha \beta \lambda \nu }-m^{\lambda }{}_{\nu \gamma }{R}_{\alpha \beta \mu \lambda }$$

$$0$$

Full text

1 Introduction

With the recent successful detection of gravitational waves, and the reasonable expectation of more to follow, there will soon be a wealth of new information about the universe allowing ever more detailed questions to be asked. But the computational methods that have served us well for today’s questions may well prove to be inadequate for the questions that arise in the near future. So it seems that there is good reason to continue to develop new approaches to computational general relativity. One such approach, known as smooth lattice general relativity, will be described in this paper. As its name suggests it is based on a lattice and it employs a metric that is locally smooth.

The smooth lattice method differs from traditional numerical methods in computational general relativity in a number of important aspects. The space-time manifold consists of a large collection of overlapping computational cells with local Riemann normal coordinates used in each cell. The computational cells are a set of vertices and legs that define small subsets of the manifold. The use of local Riemann normal coordinates in each each cell not only reduces the complexity of the evolution equations but it also explicitly incorporates the Einstein equivalence principle into the formalism. The lattice method provides an elegant separation between the topological properties of the space-time (by specifying combinatoric data such as the connections between cells, vertices etc.) and the metric properties (by specifying data such as leg-lengths, curvature components etc. within each cell). A key element of the lattice method is that it uses the second Bianchi identities to evolve the Riemann curvatures. More details of the lattice method will be given later in section (3).

Previous applications of the lattice method includes the Schwarzschild [brewin:2010-02], Oppenheimer-Snyder [brewin:2009-05] and Kasner [brewin:2014-01] space-times. Though these were important tests of the lattice method, they lacked some of the more challenging aspects expected in full 3-dimensional computational general relativity, in particular the presence of gravitational waves and their interactions with the outer boundaries on a finite computational grid. In this paper evolutions of a smooth lattice with zero shift for the Gowdy [gowdy:1971-01], Brill [brill:1959-01] and Teukolsky [teukolsky:1982-01] spacetimes will be presented. The objective is not to explore any new features of these space-times but rather to use them as examples of the smooth lattice method.

The boundaries in the Gowdy space-time will be handled using standard periodic boundary conditions while the Brill and Teukolsky space-times will require an out outgoing radiation condition. The Brill space-time adds the extra complexity of the numerical instabilities that arise from the use of a lattice adapted to the axisymmetry. These issues will be addressed in the following sections.

This class of space-times has been studied extensively by other authors. See [berger:1993-01, hern:1998-01, garfinkle:2004-01] for the Gowdy space-time, [eppley:1979-01, Choptuik:2003-01, alcubierre:2000-04, miyama:1981-01] for Brill waves and [baumgarte:1998-01, abrahams:1998-01] for Teukolsky waves.

The structure of this paper is as follows. The notation used in this paper will be defined in the following section. Sections (3,4) provide a broad summary of the smooth lattice method including details of the evolution equations on a typical lattice. The specific details of the lattice, the construction of the initial data and the evolution equations for each of the three spacetimes are given sections (5,6,7). This is followed by a short discussion on the use of the Einstein toolkit [loffler:2011-01] before the results are presented in section (9). Most of the algebraic calculations are deferred to the appendices (The transition matrices–Complete evolution equations).

2 Notation

Throughout this paper Greek letters will denote space-time indices while spatial indices will be denoted by just three Latin letters, $i,j$ and $k$. The remaining Latin letters will serve as vertex labels. One small exception to these rules will be noted in Appendix (Cartan structure equations) where Latin indices will be used (extensively) to record frame components for differential forms.

The coordinates for a typical Riemann normal frame will be denoted by either $(t,x,y,z)$ or $x^{\mu}$ while globally defined coordinates will be denoted by the addition of a tilde such as $({\tilde{t}},{\tilde{x}},{\tilde{y}},{\tilde{z}})$ or ${{\tilde{x}}}^{\mu}$. A tilde will also be used to denote tensor components in the global frame, e.g., ${\tilde{T}_{xy}}$ would be the ${\tilde{x}}{\tilde{y}}$ component of the tensor $T$ in the global coordinate frame. Note that the global coordinates are not an essential part of the smooth lattice method. They appear in this paper solely to assist in setting the initial data and also when comparing the evolved data against the exact solution or against data obtained by other numerical means (e.g., a finite difference code).

A key element of the smooth lattice method is that it employs many local Riemann normal frames. This introduces a minor bookkeeping issue – if a tensor is defined across two frames, how should its components in each frame be recorded? Let $\overline{a}$ and $\overline{b}$ be the Riemann normal frames associated with the pair of vertices $a$ and $b$. Consider a vector $v$ defined over this pair of frames. Then the components, in the frame $\overline{b}$, of the vector $v$ at vertex $a$ will be denoted by $v^{\alpha}{}_{a\overline{b}}$ while $v^{\alpha}{}_{a\overline{a}}$ denotes the components, in $\overline{a}$, of $v$ at $a$. Similar notation will be used for other tensors, for example $R^{\alpha}{}_{\beta p\overline{q}}$ would denote the components of the Ricci tensor at the vertex $p$ in the frame $\overline{q}$.

It is customary to denote the Cauchy time parameter by the symbol $t$. However, that symbol is reserved for the time coordinate of a typical local Riemann normal frame and thus some other symbol is required, for example ${\tilde{t}}$ with a corresponding time derivative operator $d/d{\tilde{t}}$. The proliferation of tildes that would follow from this choice can be avoided with the following convention – replace $d/d{\tilde{t}}$ with $d/dt$ and take the $d/dt$ to be the time derivative operator associated with the Cauchy time parameter ${\tilde{t}}$. This convention applies only to the operator $d/dt$, thus a (partial) time derivative such as $v^{\mu}{}_{,t}$ should be understood as a derivative with respect to the Riemann normal coordinate $t$.

The signature for the metric, Riemann and Ricci tensors follows that of Misner, Thorne and Wheeler [mtw:1973-01].

3 Smooth lattices

A smooth lattice is a discrete entity endowed with sufficient structure to allow it to be used as a useful approximation to a smooth geometry (which in the context of computational general relativity is taken to be a solution of the Einstein equations). The typical elements of a smooth lattice are combinatoric data such as vertices, legs, etc. and geometric data such as a coordinates, the Riemann and metric tensors and any other geometric data needed to make the approximation to the smooth geometry meaningful.

An $n$-dimensional smooth lattice can be considered as a generalisation of an $n$-dimensional piecewise linear manifold. The later are constructed by gluing together a collection of flat $n$-simplices in such a way as to ensure that the resulting object is an $n$-dimensional manifold, that the points common to any pair of $n$-simplices form sub-spaces of dimension $n-1$ or less and that the metric is continuous across the interface between every pair of connected $n$-simplices.

In a smooth lattice the cells need not be simplices, they are required to overlap with their neighbours and the curvature may be non-zero throughout each cell. The picture to bear in mind is that the cells of a smooth lattice are akin to the collection of coordinate charts that one would normally use to cover a manifold. The overlap between each pair of charts is non-trivial and allows for coordinate transformations between neighbouring charts. So too for the smooth lattice – each pair of neighbouring cells overlap to the extent that a well defined transition function can be constructed. This is an essential element of the smooth lattice formalism – it is used extensively when computing various source terms in the equations that control the evolution of the lattice (see appendix (The transition matrices) for further details). Another important feature of the smooth lattice is that each cell of the lattice need not be flat. The intention here is to better allow the smooth lattice to approximate smooth geometries than could otherwise be achieved using piecewise flat simplices (compare the approximation of a sphere by spherical triangles as opposed to flat triangles). The smooth lattice should also provide smoothly varying estimates for various quantities (for example the geodesic length of a leg) in the overlap region between a pair of cells. The use of the adjective smooth in the name smooth lattice is intended to capture the idea that all quantities on the lattice should vary smoothly (as best as possible) across the lattice.

Denote the smooth geometry by $(g,M)$ where $g$ is the metric on the $n$-dimensional manifold $M$. A smooth lattice representation of $(g,M)$ can be constructed in a number of stages, in particular, choose a set of cells $M_{i},\>i=1,2,3,\cdots$ that cover $M$, add the vertices and legs and finally add the metric data to the lattice.

The cells $M_{i},\>i=1,2,3,\cdots$ must be chosen so that each point in $M$ is contained in at least one $M_{i}$ and each point in each $M_{i}$ should also be a point in $M$. Now decorate $M$ by introducing a set of vertices $V$ and a set of legs $L$ as follows. Add one or more vertices to each cell and in each cell label one of these as the central vertex for that cell (which will later serve as the origin of a set of coordinates local to the cell). Thus each cell will contain one central vertex as well as other vertices (which are also the central vertices of other neighbouring cells). The legs $L$ of the lattice are chosen as the geodesics that connects the central vertices between pairs of neighbouring cells. Paths other than geodesics could be used but since the geodesic is defined intrinsically by the underlying smooth geometry it is a natural choice. There is, however, the issue of the uniqueness of the geodesic – if the curvature is too large or the vertices too far apart there may not exist a unique geodesic joining the pair of vertices. This problem can be overcome by a suitable choice of cells – in regions where the curvature is large the cells should be small and closely packed while in other regions, where the curvature is weak, the cells can be well spaced out. It is well known that such a construction is always possible (in the absence of curvature singularities).

The next step in the construction is to assign metric data to the cells. In each cell $M_{i}$, expand the metric around the central vertex in terms of a local set of Riemann normal coordinates $x^{\alpha}$ (see [willmore:1996-01, synge:1960-01, brewin:2009-03]), that is

[TABLE]

The coefficients $g_{\alpha\beta}$, $R_{\alpha\mu\beta\nu}$ etc. can be obtained by projecting their corresponding quantities from the smooth metric onto a local orthonormal basis on the central vertex.

At this stage the lattice is an exact copy of the original smooth geometry but with additional structure (the vertices, legs, cells, coordinates etc.). The approximation is introduced by truncating the series expansion for the metric at some finite order. The lattice will then no longer be an exact copy of the original smooth metric and should be considered an entity in its own right and will be denoted by $(g,M,V,L)$. The original smooth geometry will now be denoted by $({\tilde{g}},{\tilde{M}})$.

For the space-times considered in this paper the metric in each cell will be taken as

[TABLE]

where $g_{\alpha\beta}={\rm diag}(-1,1,1,1)$. This form of the metric will lead to estimates for the geodesic lengths that differ from that given by $({\tilde{g}},{\tilde{M}})$. By inspection of the (3.1) and (3.2) it is should be clear that for a typical leg $(p,q)$ in ${\tilde{M}}$ and $M$, the geodesic lengths, using the two metrics ${\tilde{g}}$ and $g$, will differ by a term of order $\mathcal{O}\left(RL^{5}\right)$ where $R$ and $L$ are estimates of the largest curvatures and lengths in any of the cells that contain this leg.

If $(p,q)$ is a leg in the smooth lattice then the (squared) geodesic length can be estimated (see [brewin:2009-03, synge:1960-01]) on the smooth lattice using

[TABLE]

where where $\Delta x^{\alpha}_{pq}=x^{\alpha}_{q}-x^{\alpha}_{p}$. Of course other sources of truncation errors will arise as part of the numerical evolution of the lattice data so this $\mathcal{O}\left(RL^{5}\right)$ truncation is the best that can expected at this level of approximation. To obtain higher order approximations would require not only retaining more terms in the series expansion for the metric but would also require the cells to overlap beyond nearest neighbours.

Imagine for the moment that the truncation errors on the right hand side of (3.3) where discarded. This leaves one equation that links the vertex coordinates, the leg lengths and the curvatures. It might be thought that given sufficiently many leg-lengths that the curvatures and coordinates could be computed by solving (3.3). Past experience shows that even though the equations can be solved (in some cases) the resulting evolution of the lattice did not converge to the continuum space-time. It was found that correct evolutions could be obtained by evolving either the leg-lengths and curvatures or equally by evolving the coordinates and the curvatures. Both approaches will be discussed in more detail in section (4.2).

3.1 Continuous time smooth lattices

The construction of the smooth lattice as described above would naturally lead, for the case of computational general relativity, to a structure that is discrete in both space and time. There is, however, an alternative picture in which the lattice evolves smoothly in time while retaining its discrete spatial structure. This allows for a fairly simple construction of a Cauchy initial value problem on such a lattice (as described later in the following section. For the remainder of this paper, the smooth lattice, its coordinates, leg lengths and Riemann curvatures should be considered to evolve smoothly with time.

4 Cauchy evolution of a smooth lattice

Suppose that the spacetime $({\tilde{g}},{\tilde{M}})$ can be foliated by a one parameter family of spatial hypersurfaces $\Sigma({\tilde{t}})$ (i.e., each $\Sigma({\tilde{t}})$ is a Cauchy surface in $({\tilde{g}},{\tilde{M}})$). Each element of this family could be represented by a lattice with 3-dimensional computational cells denoted by $\Sigma_{i}$. The 4-dimensional computational cells $M_{i}$ of $M$ will be taken as the space-time volume swept out by the corresponding $\Sigma_{i}$ for an infinitesimal increment in the Cauchy time parameter ${\tilde{t}}$. Thus a single $M_{i}$ is a 4-dimensional cylinder, with a 3-dimensional base $\Sigma_{i}$, that connects a pair of infinitesimally close Cauchy surfaces while the set of all $M_{i},i=1,2,3,\cdots$ fills out the space-time region between that pair of Cauchy surfaces.

The dynamical variables on a smooth lattice can be chosen to include the Riemann curvatures on the central vertex and either the (squared) leg-lengths or the Riemann normal coordinates for each vertex in each cell. In either case, the addition of the extrinsic curvatures (at the central vertex) allows the full set of evolution equations for the lattice to be given in first order form.

4.1 Lapse and shift

In the standard formulation of the Cauchy initial value problem for general relativity the lapse function and shift vectors can be freely specified at each point in the space-time. This naturally carries over to the smooth lattice by allowing the lapse function and shift vector to be freely specified on the central vertex of each cell.

In computational general relativity it is usually the case that once the lapse function and shift vector have been fully specified then there are no remaining coordinate freedoms. This is not exactly true on a smooth lattice – each cell carries its own local set of coordinates and specifying the lapse and shift at one point in that cell is not sufficient to properly constrain the coordinates on the remaining vertices. What remains is the freedom to orient the coordinate axes within each cell. Thus using boosts and spatial rotations the $t$-axis can be aligned with the world-line of the central vertex (for the case of zero shift) while the spatial axes can be given some preferred alignment with some of the remaining vertices of the cell.aaaThis picture changes slightly if the coordinates are evolved, see the comment at the end of section (4.2). This is a choice that depends on the structure of the cells and possibly on any symmetries that might exist in the space-time.

In each of the space-times considered in this paper the shift vector will be set equal to zero (i.e., the world-lines of the vertices will be normal to the Cauchy surfaces) while the lapse function will be given as a function on the set of central vertices.

4.2 Evolving the legs and coordinates

The only legs that will be evolved in a cell are those that are directly connected to the central vertex. There are two reasons for making this choice. First, legs that are not tied to the central vertex are likely to incur a larger truncation error than legs closer to the central vertex (such as those tied to that vertex). Second, there is no contribution to the leg-length from the Riemann tensor for legs directly connected to the central vertex thus avoiding any issues of accounting for time derivatives of such terms.

Consider a typical cell with central vertex $o$ and let $q$ be any of its vertices. A standard result from differential geometry, known as the first variation of arc-length [berger-marcel:2003, chavel:2006-01, hicks:1965-01], states that for a one-parameter family of geodesics, the arc-length $L_{oq}$ will evolve according to

[TABLE]

where $v^{\alpha}$ is the (forward pointing) unit tangent vector to the geodesic, $n^{\mu}$ is the (future pointing) unit tangent vector to the vertex world-line and $N$ is the lapse function. For a short leg, where the lapse and extrinsic curvatures are approximately constant across the leg, this result can be estimated by [brewin:2009-04]bbbThis paper contains a number of small errors that do not effect the final results. A corrected version can be found on ArXiv:0903.5365

[TABLE]

Since $N$ and $K_{\alpha\beta}$ are defined on the vertices there is an ambiguity in attempting to apply this equation to any leg – each leg is defined by two vertices so which vertex should supply the required values? As there is no clear reason to prefer one vertex over the other it seems reasonable to take the average from both vertices, that iscccThis result can also be obtained directly from (4.1) as shown in Appendix (Evolution of $L_{oq}$).

[TABLE]

A simple generalisation of this result can be obtained by noting that any 3-geodesic within a Cauchy surface can be arbitrarily approximated by a large sequence of short 4-geodesics of the space-time. The arc-length for each short 4-geodesic is subject to the above evolution equation and thus, on summing over all contributions to the path and taking a suitable limit, it follows that

[TABLE]

where $s$ is the proper distance along the path and ${\tilde{L}}_{oq}=\int_{o}^{q}ds$ is the arc-length of the 3-geodesic.

Using this equation to evolve the leg-lengths requires a re-appraisal of how the legs of the lattice are interpreted. In the standard formulation [brewin:2014-01], the legs of the lattice are geodesics in space-time (and will appear as chords connecting the vertices) whereas in this alternative interpretation the geodesics now lie entirely within a Cauchy surface.

The evolution equation (4.4) is suitable for simple lattices, such as the Gowdy lattice, where information about $N$ and $K_{\alpha\beta}$ can be deduced along the entire path. In all other cases, such as the Brill and Teukolsky lattices, the former evolution equation (4.3) must be used.

As the leg-lengths evolve, so too must the Riemann normal coordinates. So it is natural to ask: What are the appropriate evolution equations for the $x^{\alpha}$? A simple calculation, as detailed in [brewin:2014-01], shows that for any vertex $p$ in a cell

[TABLE]

A short independent derivation of this equation can also be found in Appendix (Evolution of $x^{i}$). Note that in choosing to evolve the coordinates, the freedom to adapt the coordinates to the lattice, as described in section (4.1), can only be imposed either on the initial Cauchy surface or at future times by applying suitable rotations.

4.3 Evolving the extrinsic curvatures

In [brewin:2014-01] the evolution equations for the extrinsic curvatures where given for the particular case of a unit lapse. The method employed in that paper can be easily repeated for the more general case of a non-constant lapse. The results are as followsdddThese equations can also be obtained directly by projecting the Arnowitt, Deser and Misner (ADM) 3+1 equations [mtw:1973-01], with zero shift, onto a local orthonormal frame.

[TABLE]

These equations apply at the central vertex where, in the Riemann normal frame of this vertex, $n^{\alpha}=\delta^{\alpha}_{t}$ and where the covariant derivatives $N_{;\alpha\beta}$ coincides with the partial derivatives $N_{,\alpha\beta}$.

4.4 Evolving the Riemann curvatures

In 4-dimensions there are 20 algebraically independent components of the Riemann tensor at any one point and in each cell these are taken to be

[TABLE]

Of these, the first 14 will be evolved while the remaining 6 will be set by applying the vacuum Einstein equations (see section (4.5)).

The evolution equations for the Riemann curvatures are based upon the second Bianchi identities. At the origin of the local frame (i.e., the central vertex) the connection vanishes and thus these equations take the simple form

[TABLE]

There is, however, a small bump in the road in using these equations to evolve the curvatures – the only data immediately available are the point values for the curvatures in each cell and thus some process must be applied to estimate the partial derivatives in each cell. It is possible to use a finite difference approximation using data from neighbouring cells but in doing so a proper account must be made of the different orientations of the neighbouring frames. This is clearly true for the spatial derivatives where neighbouring frames may differ by boosts and rotations. It is also true for the time derivatives due to progression of boosts needed to keep the world-line of the origin of the local frame normal to the Cauchy surfaces. Thus $R_{tzxy,x}$, for example, will consist not only of the raw partial derivatives (i.e., taking the raw data from neighbouring frames without regard for coordinate transformations) but also of terms that account for the boosts and rotations between neighbouring frames. The details are spelt out in full, for the particular class of lattices used in this paper, in Appendix (Source terms) leading to expression such as

[TABLE]

in which the $R_{\alpha\beta\mu\nu{{\dagger}}\gamma}$ are the raw partial derivatives of $R_{\alpha\beta\mu\nu}$ and the $m^{\alpha}{}_{\beta\gamma}$ are geometrical data built solely from the structure of the lattice (i.e., they depend only on the leg-lengths and Riemann normal coordinates). This result is very much like the usual definition of a covariant derivative. This does of course lead to a significant increase in the number of terms in each equation. The full set of equations (for a zero shift) can be found in Appendix (Complete evolution equations).

4.5 The vacuum Einstein equations

The second Bianchi identities provide no information about the time derivatives of the Riemann components such as $R_{txtx}$. Consequently such components can not be evolved but rather must be determined algebraically by applying the (vacuum) Einstein equations. Thus the 6 curvatures $R_{txtx},R_{txty}\cdots R_{tytz}$ are obtained from

[TABLE]

4.6 Constraint equations

The constraints consist not only of the four standard Hamiltonian and momentum constraints, which on a lattice take the form

[TABLE]

but also the extra constraints that arise from allowing the Riemann curvatures to be evolved. These constraints follow from the second Bianchi identities, namely

[TABLE]

Note that $R_{txyz}$ is not one of the 20 chosen $R_{\alpha\beta\mu\nu}$ but it can be computed directly using $R_{txyz}=R_{tyxz}-R_{tzxy}$.

5 Gowdy polarised cosmologies

Polarised Gowdy cosmologies on $T^{3}\times R$ are a class of solutions of the vacuum Einstein equations that posses two linearly independent spatial Killing vectors. The metric, in coordinates adapted to the symmetries, is commonly written in the form [new-kc:1998-01, alcubierre:2004-02]

[TABLE]

where $P$ and $\lambda$ are functions of $({\tilde{t}},{\tilde{z}})$ and where $\partial/\partial{\tilde{x}}$ and $\partial/\partial{\tilde{y}}$ are the two Killing vectors. Each of the spatial coordinates $({\tilde{x}},{\tilde{y}},{\tilde{z}})$ are required to be periodic (to respect the $T^{3}$ topology). The functions $P$ and $\lambda$ used in this paper are those given by New-Watt etal [new-kc:1998-01], namely,

[TABLE]

with ${\tilde{z}}$ restricted to $\left[-0.5,0.5\right]$. The domain for ${\tilde{x}}$ and ${\tilde{y}}$ can be chosen as any finite interval, e.g., $\left[0,1\right]$.

The metric is singular only at ${\tilde{t}}=0$ and consequently initial data should be set at some other time (e.g., at ${\tilde{t}}=1$ as described below). The Gowdy initial data will be evolved away from the ${\tilde{t}}=0$ singularity.

5.1 A Gowdy lattice

A lattice that represents the spatial part of this metric is rather easy to construct. Start by discretising the ${\tilde{z}}$ axis into a finite number of points labelled from [math] to $N_{z}$ with the point labelled [math] identified with that labelled $N_{z}$ (i.e., two labels for a single point). These points will soon be identified as the vertices of the lattice. Note that there are no legs at this stage, these will be added later. Now use the Killing vectors $\partial/\partial{\tilde{x}}$ and $\partial/\partial{\tilde{y}}$ to drag the discretised ${\tilde{z}}$ axis along the ${\tilde{x}}$ and ${\tilde{y}}$ axis. The legs of the lattice can now be constructed as the space-time geodesics that connect pairs of points (now taken as vertices of the lattice). This leads to the simple lattice shown in figure (1) consisting of $N_{z}$ computational cells labelled from [math] to $N_{z}$ with cell [math] identified with cell $N_{z}$. This lattice contains three classes of legs, one for each of the three coordinate axes, namely, $L_{xx},L_{yy}$ and $L_{zz}$. Other data that must be carried by the lattice include the extrinsic curvatures, $K_{\alpha\beta}$, the Riemann curvatures, $R_{\alpha\beta\mu\nu}$ and the lapse function $N$.

Consider a typical computational cell, as shown in figure (1), and ask the question: How should the Riemann normal frame be constructed? Let $\partial_{\alpha}$ be the unit basis vectors for the Riemann normal frame. Now choose the origin of the Riemann normal frame to be (permanently) attached to the central vertex. Next, use boosts to ensure that $\partial_{t}$ is normal to the Cauchy surface, then use rotations to ensure that the vertices of $L_{zz}$ lie on the $z$-axis and also for the vertices of $L_{xx}$ to lie in the $xz$-plane. Given the symmetries of the Gowdy space-time it is no hard to appreciate that the $(t,x,y,z)$ coordinates of the seven vertices of the cell $M_{p}$ will be of the following form

[TABLE]

where the time coordinate is given by $2t=-K_{\alpha\beta}x^{\alpha}x^{\beta}$ (see [brewin:2010-03]).

Note that this construction also ensures that the Riemann normal axes are aligned with their Gowdy counterparts (as a consequence of the Gowdy metric being diagonal).

5.2 Initial data

A straightforward computation on the Gowdy metric reveals that there are three non-trivial extrinsic curvatures, ${\tilde{K}}_{xx},{\tilde{K}}_{yy}$ and ${\tilde{K}}_{zz}$ and five non-trivial Riemann curvatures, ${\tilde{R}}_{xyxy},{\tilde{R}}_{xzxz},{\tilde{R}}_{yzyz},{\tilde{R}}_{txxz}$ and ${\tilde{R}}_{tyyz}$. The lattice values for the extrinsic and Riemann curvatures, $K_{\alpha\beta}$ and $R_{\alpha\beta\mu\nu}$, were computed by projecting their counterparts, ${\tilde{K}}_{\alpha\beta}$ and ${\tilde{R}}_{\alpha\beta\mu\nu}$, onto the local Riemann normal frame. This provides not only a way to identify the non-trivial components on the lattice but also a simple way to assign the initial data.

The leg-lengths $L_{xx},L_{yy}$ and $L_{zz}$ were set as follows. The $L_{xx}$ were computed as the length of the geodesic connecting $(1,0,0,{\tilde{z}})$ to $(1,\delta{\tilde{x}},0,{\tilde{z}})$ with $\delta{\tilde{x}}=0.0001$. A similar approach was used to compute the $L_{yy}$ this time using the points $(1,0,0,{\tilde{z}})$ and $(1,0,\delta{\tilde{y}},{\tilde{z}})$ with $\delta{\tilde{y}}=\delta{\tilde{x}}=0.0001$. A common value for $L_{zz}$ was chosen for all cells, namely

[TABLE]

This in turn required the ${\tilde{z}}$ coordinate to be unequally spaced from cell to cell. Starting with ${\tilde{z}}_{0}=-0.5$ the successive ${\tilde{z}}_{p}$ for $p=1,2,3\cdots N_{z}-1$ where found by treating the equation

[TABLE]

as a non-linear equation for ${\tilde{z}}_{p}$ given ${\tilde{z}}_{p-1}$.

5.3 Evolution equations

The evolution equations for $L_{xx},L_{yy}$ and $L_{zz}$ follow directly from equation (4.4) by making appropriate use of the symmetries built into the Gowdy lattice, in particular that the legs are aligned to the coordinate axes and thus $v^{\alpha}_{ox}=(0,1,0,0)$, $v^{\alpha}_{oy}=(0,0,1,0)$ and $v^{\alpha}_{oz}=(0,0,0,1)$ while rotational symmetry ensures that the integrand in (4.4) is constant along the $x$ and $y$ axes. This leads to the following evolution equations for $L_{xx},L_{yy}$ and $L_{zz}$ in cell $p$,

[TABLE]

and where $s$ is the arc-length along the leg connecting successive cells (i.e., along the ${\tilde{z}}$-axis of the lattice) and where the limits $(p,p+1)$ are understood to denote the corresponding vertices.

The evolution equations for the extrinsic and Riemann curvatures can be constructed in at least two ways. In the first approach the evolution equations for the ${\tilde{K}}_{\alpha\beta}$ and ${\tilde{R}}_{\alpha\beta\mu\nu}$ can be projected onto the the local Riemann normal frame. The second approach is to impose the known symmetries on the the complete set of equations given in Appendix (Complete evolution equations). Both approaches lead to the following set of equations for the extrinsic curvatures,

[TABLE]

and for the Riemann curvatures,

[TABLE]

where

[TABLE]

5.4 The lapse function

The lapse function can be freely chosen across the lattice either by way of an explicit function (e.g. $N=1$) or by evolving the lapse along with other lattice data. This second choice will taken in this paper where three different methods for evolving the lapse will be used, namely

[TABLE]

where ${\rm TrK}=K_{xx}+K_{yy}+K_{zz}$. The $1+\log$ and harmonic lapse equations are standard gauge choices and need no explanation while the third equation, as its name suggests, is designed to track the exact solution. This exact lapse equation can be obtained as follows. First note that for the exact solution $N^{2}={\tilde{g}}_{zz}$. Then use $d{\tilde{g}}_{zz}/dt=-2N{\tilde{K}}_{zz}$ to obtain $dN/dt=-{\tilde{K}}_{zz}$ whereupon the result follows by noting that ${\tilde{K}}_{zz}={\tilde{g}}_{zz}K_{zz}=N^{2}K_{zz}$.

Many other choices are of course possible but those just given stand out as they allow for a direct comparison with either the exact solution (5.1–5) or with the results from the Cactus code.

Initial values for the lapse will be discussed later in section (9.1).

5.5 Constraints

The only constraints that survive under the symmetries inherent in the Gowdy space-time are (4.34,4.37,4.38) and can be written as

[TABLE]

where $R_{xyxy{{\dagger}}z}$, $m^{x}{}_{zx}$ and $m^{y}{}_{zy}$ are given by (5.20,5.21). Note also that trivial factors have been cleared from the first two equations. This set of constraints were not imposed during the evolution but were instead used as a quality control on the evolved data (see section (9.1)).

5.6 Numerical dissipation

It was found that for some choices of the lapse function, most notably the $1+\log$ choice, the addition of some numerical dissipation could significantly prolong the evolution.

The particular form of numerical dissipation used here is based upon the familiar Kreiss-Oliger approach in which an additional term is added to the right hand side of selected evolution equations, in our case, the evolution equations for the extrinsic and Riemann curvatures. In each case the modified evolution equation in cell $p$ was of the form

[TABLE]

where $\epsilon$ is a small number (in the results described below $\epsilon=0.8$). The first term on the right hand side is the right hand side of the evolution equations (5.10-5.3) while the second term is a naive approximation to $\epsilon L_{zz}^{5}d^{6}Y/ds^{6}$. The important point is that the dissipation scales as $\mathcal{O}\left(L_{zz}^{5}\right)$ and thus will vanish in the limit as $L_{zz}\rightarrow 0$.

6 Brill waves

Brill waves [brill:1959-01] are time and axisymmetric solutions of the vacuum Einstein equations generated by initial data of the form

[TABLE]

in which $({\tilde{\rho}},{\tilde{\phi}},{\tilde{z}})$ are cylindrical polar coordinates and where $\psi({\tilde{\rho}},{\tilde{z}})$ and $q({\tilde{\rho}},{\tilde{z}})$ are a class of functions subject to the conditions of asymptotic flatness, the vacuum Einstein equations and reflection symmetry across both ${\tilde{z}}=0$ and ${\tilde{\rho}}=0$. The reflection symmetry across ${\tilde{\rho}}=0$ follows from the condition that the data be well behaved at ${\tilde{\rho}}=0$. However, the condition that the data be reflection symmetric across ${\tilde{z}}=0$ has no physical basis and is introduced only to reduce the bulk of the numerics (i.e., the data can be evolved in the quarter plane (${\tilde{\rho}}>0,{\tilde{z}}>0$) rather than the half plane (${\tilde{\rho}}>0,|{\tilde{z}}|<\infty$)).

Brill showed that the initial data will have a finite ADM mass when the functions $q$ and $\psi$ behave as $q=\mathcal{O}\left({\tilde{r}}^{-2}\right)$ and $\psi=1+\mathcal{O}\left({\tilde{r}}^{-1}\right)$ as ${\tilde{r}}\rightarrow\infty$ where ${\tilde{r}}^{2}={\tilde{\rho}}^{2}+{\tilde{z}}^{2}$. He also showed that for the initial data to be well behaved near the ${\tilde{\rho}}=0$ coordinate singularity, $q$ must behave like $q=\mathcal{O}\left({\tilde{\rho}}^{2}\right)$ as ${\tilde{\rho}}\rightarrow 0$ which can also be expressed as

[TABLE]

while the reflection symmetric conditions on $q$ and $\psi$ requires

[TABLE]

The condition that $\psi=1+\mathcal{O}\left({\tilde{r}}^{-1}\right)$ as ${\tilde{r}}\rightarrow\infty$ was implemented using a standard mixed outer boundary condition,

[TABLE]

Finally, the vacuum Einstein equations requires $\psi$ to be a solution of the Hamiltonian constraint which in this case takes the form

[TABLE]

where $\nabla^{2}$ is the (flat space) Laplacian in the cylindrical coordinates $({\tilde{\rho}},{\tilde{\phi}},{\tilde{z}})$. The three momentum constraints provide no new information as they are identically satisfied for any choice of $q$ and $\psi$.

6.1 Eppley Initial data

The function $q({\tilde{\rho}},{\tilde{z}})$ was chosen as per Eppley [eppley:1977-01], namely

[TABLE]

with $n=5$ (any $n\geq 4$ would be sufficient to satisfy $q=\mathcal{O}\left({\tilde{\rho}}^{-2}\right)$). The parameter $a$ governs the wave amplitude with $a=0.01$ in the results presented below. Even though this is a weak amplitude it is sufficient to test the lattice method.

The Hamiltonian constraint (6.6), subject to the boundary conditions (6.4–6.5), was solved for $\psi$ using standard second order centred finite differences (including on the boundaries). The grid comprised $2048\times 2048$ equally spaced points covering the rectangle bounded by ${\tilde{\rho}}={\tilde{z}}=0$ and ${\tilde{\rho}}=z=20$. The finite difference equations were solved (with a maximum residual of approximately $10^{-13}$) using a full multigrid code. The full Brill 3-metric was then constructed using the reflection symmetry across $z=0$ and the rotational symmetry around the $z$-axis.

Since the Brill initial data is axisymmetric it is sufficient to use a 2-dimensional lattice on which to record the initial data for the lattice. An example of such a lattice is shown in figure (2). Each cell contains legs that are (at ${\tilde{t}}=0$) aligned to the Brill $({\tilde{\rho}},{\tilde{z}})$ axes as well as a set of diagonal legs. A full 3-dimensional lattice could be constructed by rotating this 2-dimensional lattice around the symmetry axis (as indicated in figure (2)). In our computer code the right portion of lattice covered the domain bounded by ${\tilde{\rho}}={\tilde{z}}=0$, ${\tilde{z}}=\pm 5$ and ${\tilde{\rho}}=5$ while the left portion was obtained by reflection symmetry across ${\tilde{\rho}}=0$. This places the symmetry axis mid-way from left to right across the lattice (this is the blue axis shown in figure (2)).

Each cell of the lattice contains 9 vertices $o,a,b,\cdots,h$ plus one additional vertex $p$ connected just to the central vertex $o$. The purpose of the extra vertex $p$ is that the collection of all such vertices defines the image of the 2-dimensional lattice under the action of the rotational symmetry. Figure (2) shows two such additional lattices in which each yellow leg has vertices of the form $(o,p)$.

In each cell the local Riemann normal coordinates $(t,x,y,z)$ were chosen as follows

[TABLE]

for some set of numbers $x_{a},z_{a},\cdots y_{p}$ and where the labels $o,a,b,\cdots,h$ follow the pattern shown in figure (4).

The leg-lengths and Riemann normal coordinates were set by first distributing the $N_{x}\times N_{z}$ vertices as equally spaced points in the $({\tilde{\rho}},{\tilde{z}})$ domain, $(-5,-5)$ to $(5,5)$, and then integrating the geodesic equations as a two-point boundary value problem for each leg in each cell.

The remaining initial data on the lattice consists of the non-zero components of the Riemann and extrinsic curvatures along with either the leg-lengths or the vertex coordinateseeeThe choice depends on which evolution scheme is used – evolving the leg-lengths or evolving the coordinates.. Given the symmetries of the Brill metric it is not hard to see that the there are only 4 non-trivial extrinsic curvatures, $K_{xx}$, $K_{yy}$, $K_{zz}$ and $K_{xz}$ and 8 non-trivial extrinsic curvatures, $R_{xyxy}$, $R_{yzyz}$, $R_{xzxz}$, $R_{xyyz}$, $R_{txxz}$, $R_{tzxz}$, $R_{tyxy}$ and $R_{tyyz}$. Each of these 12 curvatures were given initial values by projecting their counterparts from the Brill metric (extended to 3+1 form using a unit lapse and setting $d\psi/dt=dq/dt=0$ at ${\tilde{t}}=0$) onto the local orthonormal frame.

6.2 Evolution equations

The initial data just described has only 12 non-trivial components for the Riemann and extrinsic curvatures. It is easy to see that that this situation is preserved by the evolution equations. For example, equation (4.9) shows that $dK_{xy}/dt=0$ for this particular set of initial data. Thus all of the symmetries in the initial data will be preserved throughout the evolution (e.g., $K_{xy}$ will remain zero for all time). This leads to the following set of evolution equations for the 4 extrinsic curvatures

[TABLE]

while the evolution equations for the 8 Riemann curvatures are

[TABLE]

where $m^{x}{}_{yy}$, $m^{x}{}_{zx}$ and $m^{x}{}_{zz}$ are solutions of

[TABLE]

where $v^{\alpha}_{ab\overline{c}}=x^{\alpha}_{b\overline{c}}-x^{\alpha}_{a\overline{c}}$. The equations for $m^{x}{}_{yy}$, $m^{x}{}_{zx}$ and $m^{x}{}_{zz}$ were obtained by a simple application of equation (.16) to the $xz$-plane (leading to equations (6.25) and (6.26)) and the $yz$-plane (leading to equation (6.24)).

The final set of evolution equations required are those for the leg-lengths or the vertex coordinates. In contrast to the Gowdy lattice it was decided to evolve the vertex coordinates. There are two reasons for doing so. First, the above evolution equations for the $R_{\alpha\beta\mu\nu}$ refer directly to the vertex coordinates and second, solving the coupled set of non-linear equation (3.3) for the the vertex coordinates involves not only extra work but was observed to lead to asymmetric evolutions (i.e., the evolved data failed to be reflection symmetric across the symmetry axis). This loss of symmetry was attributed to the algorithm [brewin:2014-01] used to solve these equationsfffThe algorithm in [brewin:2014-01] computes the coordinates one by one visiting the vertices in a clock wise order. But for two cells on either side of the symmetry axis, one cell should be processed clockwise and the other anti-clockwise..

6.3 Numerical dissipation

Other authors [Choptuik:2003-01, garfinkle:2001-01] have noted that the singular behaviour of the evolution equations on the symmetry axis can cause numerical instabilities to develop along the symmetry axis. This problem can be avoided by either using a fully 3-dimensional formulation (which is computationally expensive) or mitigated by introducing numerical dissipation. Similar instability problems were expected on the 2-dimensional axisymmetric lattice. By direct experiment it was found that good damping of the numerical instabilities could be obtained by applying a Kreiss-Oliger dissipation to the evolution equations. The standard practice is to weight the dissipation term by powers of the discretisation scale (i.e., powers of $L$) to ensure that the dissipation terms do not dominate the truncation errors inherent in the numerical integrator. For a 4th-order Runge-Kutta integrator (as used here) this would require a dissipation term of order $\mathcal{O}\left(L^{6}\right)$ which would be the case for a 6th-order derivative term (as used in the Gowdy lattice (5.6)). However, on this simple Brill lattice, where cells interact only by nearest neighbours, the best that can be done is to use a 2nd-derivative dissipation term. The choice used in the results given below was

[TABLE]

where $\epsilon$ is a small number and the first term on the right hand side is time derivative without dissipation while the second term is a crude estimate of $\mathcal{O}\left(L^{2}\right)\nabla^{2}Y$ on the cell (the subscripts correspond to the vertices displayed in figure (4)). The dissipation was applied only to the Riemann curvatures as no significant gains were noted when the dissipation was also applied to the extrinsic curvatures. In the results presented below $\epsilon=1.0$ (this was the smallest value of $\epsilon$ that allowed the evolution to remain stable to at least $t=10$).

6.4 Inner boundary conditions

Figure (2) show three copies of the 2-dimensional lattice sharing the common symmetry axis. Away from the symmetry axis the three copies of the lattice provide sufficient data to estimate $y$ derivatives of data on the lattice. However, this construction clearly fails at the symmetry axis. One consequence of this can be seen in equation (6.24) which, when expressed in terms of the coordinates and leg-lengths, leads to $m^{x}{}_{yy}\approx-(1/L_{yy})(dL_{yy}/dx)$ where $x$ is the proper distance measured along the $x$-axis. This shows that $m^{x}{}_{yy}$ is singular on the symmetry axis (where $L_{yy}=0$). The upshot is that any $y$ derivative, on this choice of lattice, will by singular on the symmetry axis (e.g., all of the $y$ derivatives in equations (4.13–4.26)).

One approach to dealing with this problem is to return to equations (4.13–4.26) and make direct use of the rotational symmetry to express all of the $y$ derivatives in terms of the (manifestly non-singular) $x$ derivatives on the symmetry axis. As an example, let $V_{\alpha\beta}$ be the components of a tensor $V$ on the lattice. Now consider a copy of the lattice rotated by $\pi/2$ about the symmetry axis. Denote the components of $V$ on the second lattice by $V^{\prime}_{\alpha\beta}$. Then $V^{\prime}_{\alpha\beta}=V_{\alpha\beta}$ by rotational symmetry. However, on the symmetry axis the coordinates for both lattices are related by $x^{\prime}=y$, $y^{\prime}=-x$ and $z^{\prime}=z$ thus the usual tensor transformation law would give $V^{\prime}_{xy}=-V_{yx}$. But $V^{\prime}_{xy}=V_{xy}$ and thus $V_{xy}=-V_{yx}$ on the symmetry axis. Now suppose $V_{\alpha\beta}=W_{\alpha,\beta}$ for some tensor $W$. It follows that $W_{x,y}=-W_{y,x}$ on the symmetry axis. This idea can be applied to any tensor on the lattice in particular to the derivatives of $R_{\alpha\beta\mu\nu}$.

It is also possible to gain information about the curvature components by considering a rotation of $\pi$ rather than $\pi/2$. Following the steps described above, the result is that any component of a tensor with an odd number of $x$ indices will be anti-symmetric across the symmetry axis while the remaining components will be symmetric. This shows immediately that $K_{xz}$, $R_{xyyz}$, $R_{tyxy}$ and $R_{tzxz}$ must vanish on the symmetry axis.

The upshot is that the evolution equations (4.13–4.26) can be reduced, on the symmetry axis, to just 5 non-zero equations

[TABLE]

Though these equations are non-singular there remains a numerical problem with cells near the symmetry axis – their proximity to the symmetry axis can lead to instabilities in the evolution.

A better approach, described in more detail below, is to excise a strip of cells containing the symmetry axis (as shown in figure (2)) and to interpolate from outside the strip to recover the time derivatives of the Riemann curvatures within the strip. This, along with numerical dissipation, proved to be crucial in obtaining stable evolutions.

The interpolation near the symmetry axis was implemented as follows. The cells of the 2-dimensional lattice where indexed by rows and columns aligned to the ${\tilde{x}}$ and ${\tilde{z}}$ axes. Each cell was given an index pair such as $(i,j)$ with $i$ denoting the number of columns from the ${\tilde{x}}=0$ axis (i.e., the symmetry axis) and $j$ the number of rows from the ${\tilde{z}}=0$ axis. The interpolation used data from the cells $i=3,4,5,6,7$, for a given $j$, to supply data for the cells with $i=-2,-1,0,1,2$, for the same $j$. In each case the interpolation was tailored to respect the known symmetry of the data across the symmetry axis. Thus for $dR_{xyxy}/dt$, which is symmetric across ${\tilde{x}}=0$, a polynomial of the form $y(x)=a_{0}+a_{2}x^{2}+\cdots a_{8}x^{8}$ was used. For anti-symmetric data the polynomial was of the form $y(x)=a_{1}x+a_{3}x^{3}+\cdots a_{9}x^{9}$. The five coefficients $a_{0},a_{2},\cdots a_{8}$ and $a_{1},a_{3},\cdots a_{9}$ were determined using trivial variations of standard methods for polynomial interpolation. The choice of interpolation indices $i=3,4,5,6,7$, which correspond to the light blue strip in figure (2), was found by trial and error as it gave stable evolutions (in conjunction with the numerical dissipation) without being overly expensive.

There is a simple variation on this interpolation scheme in which the data from the symmetry axis (i.e., equations (6.28–6.4)) is included in the data used to build the polynomial. Thus data on the cells $i=0,3,4,5,6,7$ would be used to build data for cells $i=-2,-1,1,2$. The evolutions that resulted form this construction were highly unstable and crashed at approximately $t=4.7$.

6.5 Outer boundary conditions

The outer boundary of the lattice is defined to be a skin of cells one cell deep on the outer edges of the lattice (as indicated by the orange region in figure (2)). In each of the boundary cells the Riemann and extrinsic curvatures were evolved by way of an outgoing radiation boundary condition of the form

[TABLE]

where $f$ is one of the Riemann and extrinsic curvatures and $n$ is the outward pointing unit normal to the cell (at the central vertex). The ${\tilde{x}}^{i}$ are constants set equal to the Brill the coordinates $({\tilde{\rho}},{\tilde{z}})$ of the central vertex at $t=0$. Finally, ${\tilde{r}}=({\tilde{\rho}}^{2}+{\tilde{z}}^{2})^{1/2}$. The leg-lengths and Riemann normal coordinates in each cell were not evolved but rather copied across from the nearest inward neighbouring cell.

This is an extremely simplistic set of boundary conditions (particularly so for the leg-lengths and coordinates). It was chosen simply to get a numerical scheme up and running. The surprise it that it works very well (as discussed below in section (9.2)).

6.6 Constraints

Only five of the ten constraints (4.34–4.43) survive once the axisymmetry of the Brill space-time is imposed. The surviving equations are (4.34,4.35,4.37,4.38,4.40) and can be written in the form

[TABLE]

where some simple numerical factors have been factored out.

7 Teukolsky linearised waves

The results for the Gowdy and Brill spacetimes are promising but a proper test of the smooth lattice method requires that it be applied to truly 3-dimensional data, i.e., initial data devoid of any symmetries such as the Teukolsky linearised waves [teukolsky:1982-01] described by the metric

[TABLE]

where

[TABLE]

and where $Q(x)$ is an arbitrary function of $x$. Note that this form of the metric differs slightly from that given by Teukolsky. Here the function $F$ has been expressed as an explicit combination of ingoing and outgoing waves (thus ensuring time symmetric initial data). Note also that the derivatives of $F$ are taken with respect to ${\tilde{r}}$ rather than $x$ as used by Teukolsky. Consequently, the signs of the odd-derivatives of $F$ in the expressions for $A$, $B$ and $C$ have been flipped.

Following Baumgarte and Shapiro [baumgarte:1998-01], the function $Q(x)$ was chosen to be

[TABLE]

as this produces initial data describing a compact wave centred on the origin with a wave amplitude controlled by the parameter $a$.

Note that the metric (7.1) is not an exact solution of the vacuum Einstein equations but rather a solution of the linearised equations in the sense that $G_{ab}(g)=\mathcal{O}\left(a^{2}\right)$.

This form of the metric requires some care when setting the initial data near ${\tilde{r}}=0$ (where the coordinates are singular). A better choice is to express the metric in standard Cartesian coordinates. At the moment of time symmetry, ${\tilde{t}}=0$, the Cartesian components, $h_{ij}$, of the 3-metric are given by

[TABLE]

where ${\tilde{r}}=({\tilde{x}}^{2}+{\tilde{y}}^{2}+{\tilde{z}}^{2})^{1/2}$.

The 3-dimensional lattice was built by a simple generalisation of the 2-dimensional lattice used for the Brill waves. The grid was built from a set of $N_{x}\times N_{y}\times N_{z}$ equally spaced points in a the 3-dimensional volume bounded by $|x|=|y|=|z|=5$. The points were then identified as the vertices of the lattice while on each of the ${\tilde{x}}{\tilde{y}}$, ${\tilde{x}}{\tilde{z}}$ and ${\tilde{y}}{\tilde{z}}$ planes, legs were added in exactly the same pattern as for the 2-dimensional Brill lattice, recall figure (2). Consequently many of the ideas discussed in regard to the Brill lattice carry over to the this lattice. Initial data for the coordinates and leg-lengths were assigned by integrating the geodesic equations as two-point boundary problems for each leg of the lattice (this was time consuming but only needed to be done once). The outer boundary conditions were exactly as per equation (6.33) but on this occasion applied to all six faces of the lattice. Geodesic slicing was used (i.e., zero shift and unit lapse) and as there are no symmetries, the full set of evolution equations (4.6–4.11) and (4.13–4.26) were used (see also Appendix (Complete evolution equations)). The implementation of the numerical dissipation is in this case slightly different to that for the 2-dimensional lattice. The appropriate version of (6.27) for the 3-dimensional lattice is

[TABLE]

where the sum on the right hand side includes contributions from the 6 immediate neighbouring cells. The term in the second set of brackets in this expression is an approximation to $\mathcal{O}\left(L^{2}\right)\nabla^{2}Y$ and thus will converge to zero on successively refined lattices.

Since the Teukolsky space-time carries no symmetries it follows that none of the constraints (4.34–4.43) will be trivially satisfied throughout the evolution. Including results for all 10 constraints is somewhat of an overkill so results will be presented (in section (9.3)) for just the Hamiltonian constraint, namely,

[TABLE]

8 Cactus

The combination of the open source code Cactus [goodale:2002-01] and the Einstein Toolkit [loffler:2011-01] (collectively referred to here as the Cactus code) provide a well understood framework for computational general relativity. The Cactus code was used largely out of the box but with some simple extensions for setting the initial data for the Brill and Teukolsky space-times. A new thorn was written for the Brill space-time to set the initial data from the discretised metric provided by the same multigrid code used to set the lattice initial data. For the Teukolsky metric the EinsteinInitialData/Exact thorn was extended to include the exact 3-metric given in equations (7.7–7.12). These changes were made to ensure that the lattice and Cactus evolutions were based on exactly the same initial data.

The Cactus initial data were built over the same domain as used in the corresponding lattice initial data. The initial data were integrated using the standard BSSN and ADM thorns. The BSSN thorn used a 4th order Runge-Kutta integrator and artificial dissipation was applied to all dynamical variables with a dissipation parameter equal to $0.1$. The ADM integrations used a two-step iterated Crank-Nicholson scheme without artificial dissipation. The time step in each case was chosen to ensure a Courant factor of $1/8$.

The Cactus code does not provide values for the components of either the 3 or 4 dimensional Riemann tensor. However the spatial components, such as ${\tilde{R}}_{xyxy}$, can be reconstructed from the 3 dimensional components of the Ricci tensor and metric using a combination of the Gauss-Codazzi equations

[TABLE]

and the equation

[TABLE]

where ${\tilde{h}}_{\alpha\beta}$ is the 3-metric, ${\tilde{R}}_{\alpha\beta}$ is the 3-Ricci tensor and ${\tilde{R}}={\tilde{h}}^{\alpha\beta}{\tilde{R}}_{\alpha\beta}$.

Since the Cactus and lattice data are expressed in different frames some post-processing of the data is required before the two sets of data can be compared. There are two aspects to this, first, mapping points between the respective spaces (e.g., given a point in the Cactus coordinates what is the corresponding point in the lattice?) and second, comparing the data at those shared points. Recall that when constructing the initial data for the Brill and Teukolsky lattices, the vertices of the lattice were taken as the uniformly distributed grid points in the Brill and Teukolsky coordinates. This correspondence is preserved throughout the evolution by the zero shift condition. This is not the case for the Gowdy space-time where the initial data was constructed on an unequally spaced grid (see section (5.1)) while in contrast the Cactus code uses an equally spaced grid. In this case the conversion of tensor components, such as ${\tilde{R}}_{\alpha\beta\mu\nu}$, from the Cactus data into a form suitable for comparison with the lattice data entails two steps, first, the tensor is projected onto a local orthonormal frame, second, the radial ${\tilde{z}}$ coordinate is converted to a radial proper distance ${\tilde{s}}$. Since the Gowdy metric is diagonal the projection onto the coordinate aligned orthonormal frame is trivial, for example $R_{xyxy}={\tilde{h}}_{xx}{\tilde{h}}_{yy}{\tilde{R}}_{xyxy}$, while the proper distance between successive grid points can be computed by

[TABLE]

where the limits $(i,i+1)$ are understood to represent the corresponding grid points. The integral was estimated by a cubic polynomial based on the grid points $(i-1,i,i+1,i+2)$.

9 Results

The evolution equations for the Brill and Teukolsky lattices were integrated using a 4th order Runge-Kutta routine with a fixed time step $\delta t$ chosen to satisfy a Courant condition of the form $\delta t<C\delta L$ where $\delta L$ is the shortest leg-length on the lattice and where $C$ is a Courant factor with $0<C<1$. The same integration scheme was used for the Gowdy lattice apart from one small change where the Courant condition was based upon $N\delta t<C{\rm min}(L_{zz})$ where $N$ is the largest lapse on the lattice. This Courant condition uses the shortest $L_{zz}$ for the simple reason that the evolution equations (5.7,5.8) for $L_{xx}$ and $L_{yy}$ admit a re-scaling of $L_{xx}$ and $L_{yy}$ and thus their values can not influence $\delta t$.

A trial and error method was first used to find any time step that yielded a stable evolution (despite the cost). This allowed a more informed judgement to made by careful examination of the history of the leg-lengths. Thus for the Gowdy lattices the time step was chosen as $\delta t=0.0512/N_{z}$ corresponding to a Courant factor of $1/20$, while for the Brill and Teukolsky lattices the time step, with $C=1/8$, was set by $\delta t=1.25/(N_{z}-1)$.

9.1 Gowdy

There are two obvious tests that can be applied to the lattice data, first, a comparison against the exact data and, second, a comparison against numerical results generated by the Cactus code. Other tests that can be applied include basic convergence tests as well as observing the behaviour of the constraints.

The initial data for the lapse was chosen according to the comparison being made. The comparisons with the Cactus data were based on a unit lapse, $N=1$, while the comparisons with the exact solution used initial values taken from the exact solution, $N=e^{\lambda}/4$ at ${\tilde{t}}=1$.

The dissipation parameter $\epsilon$ (see equation (5.6)) was set equal to $0.8$ (which was found by trial and error as the smallest value that ensured good stability for the $1+\log$ lapse). The integral in equation (5.9) was estimated using a 4th order interpolation built from 5 cells centred on this leg.

Selected results can be seen in figures (5–9) and show that the lattice method works well with excellent agreement against the exact and numerical solutions. Note that since the lattice expands by factors of order 100, the $L_{zz}$ have been uniformly scaled to squeeze the lattice into the range $[-0.5,0.5]$. Figure (5) shows a comparison of the original and scaled data. Figures (8,9) show the behaviour of selected constraints as well as basic convergence tests.

9.2 Brill

The results for the Brill initial data are shown in figures (10–13). In all cases the dissipation parameter $\epsilon$ for the lattice was set equal to $1.0$ (except as noted in figure (13)). The Cactus BSSN data was computed on a full 3-dimensional grid and thus there is no reason to expect any instabilities on the symmetry axis. This allows a much small dissipation parameter, $\epsilon=0.1$, to be used for the BSSN evolutions. The Cactus ADM thorn does not appear to support any form of Kreiss-Oliger numerical dissipation.

The expected behaviour for the Brill wave is that the curvature will be propagated away from the symmetry axis with the wave hitting the edges of the outer boundary by about $t=5$ followed by the four corners by about $t=7$ and will completely cross the boundary by about $t=10$. As the wave moves across the grid it should leave zero curvature in its wake (though the extrinsic curvatures need not return to zero).

The results for all three methods at $t=5$ are shown in figure (10) where it is clear that though there is some good agreement in the propagation of the main the wave there are also some notable differences. The ADM method shows a series of parallel waves propagating in from the outer boundary towards the symmetry axis (such waves will later be referred to as boundary waves, these waves are particular evident in movies from $t=0$ to $t=10$) while the BSSN data shows a non-propagating bump close to the origin. In contrast the lattice data shows a smooth behaviour in the wave with no apparent boundary waves nor any sign of a bump. By $t=10$ (see figure (11)) the ADM data shows not only the boundary waves but also reflected waves from the outer boundary. Similar reflected waves can also be seen in the BSSN results though with a significantly smaller amplitude. The bump in the BSSN data has remained in place and has grown in amplitude. The lattice data shows no signs of reflection but there is a very small bump that correlates with the wings of the BSSN bump.

It is reasonable to ask why the three methods should give such different results in the region behind the main wave. The smooth profile in the lattice data might be due to the large dissipation parameter compared to that used in the ADM and BSSN data. The boundary waves in the ADM data are clearly associated with the boundary conditions while the cause of the bump in the BSSN data is not so easy to identify from these plots. A more detailed analysis will be given later when discussing the Teukolsky data where similar behaviour was observed.

The effects of changing dissipation parameter on the evolution of the lattice data is shown in figure (13). This shows clearly how crucial the numerical dissipation is in controlling the instabilities. The figure also shows that despite the significant dissipation ($\epsilon=1.0$) required to suppress the axis instability, the broad features of the main wave are largely unaffected.

Figure (12) shows the behaviour of the constraints $C_{1}$ (6.34) and $C_{4}$ (6.6) over the period $t=0$ to $t=10$. The remaining three constraints are not shown as they show much the same behaviour. Each plot contains four curves corresponding to different lattices scales, $N_{z}=101$ (red), $N_{z}=201$ (blue), $N_{z}=401$ (green) and $N_{z}=801$ (black). These show that the constraints appear to decrease as $N_{z}$ is increased. It also appears that the constraints settle to a non-zero value as $t$ increases. This could be due to truncation errors inherent in the solution of the Hamiltonian equation (6.6) coupled with the interpolation to the lattice (though this claim was not tested). The two bumps in the left figure, one just after $t=5$ and one close to $t=10$ are most likely due to reflections from the outer boundary (this too was not tested).

9.3 Teukolsky

The Teukolsky data is specified on a full 3-dimensional grid/lattice and is thus not susceptible to the axis instability seen in the Brill data. This allows for a much smaller dissipation parameter to be used for the lattice, ADM and BSSN codes, in this case $\epsilon=0.1$.

The results for the Teukolsky initial data are shown figures (14–18) and bear some similarities with the results for the Brill initial data. However, in this case the boundary and reflected waves appear to be much less noticeable while the bump in the BSSN data is still present and is more pronounced than in the Brill wave data.

The plots in figure (17) show that the bump in the BSSN data is a numerical artefact. The figure shows that as the spatial resolution is decreased (i.e., increasing $N_{z}$) the amplitude of the bump, at $t=5$, decreases. The figure also shows that the amplitude of the bump grows with time. No attempt was made to determine the source of the bump.

In order to better understand the influence of the outer boundary condition on the evolution it was decided to run the lattice, ADM and BSSN codes on two different sets of initial data, each with the same spatial resolution but with one grid twice the size of the other (i.e., one grid had boundaries at $\pm 5$ and the other at $\pm 10$). The influence of the outer boundary condition on the evolution was then be measured by comparing the evolution on the common region. The results are shown in figure (18). The right panel shows the evolution of $R_{xyxy}$ on the lattice on both grids with $N_{z}=101$ for the red curve and $N_{z}=201$ for the blue curve. Notice how the red curve lies entirely on top of the blue curve even as the wave passes through the $\pm 5$ boundary. The left panel shows the difference in $R_{xyxy}$ between the two grids for the lattice data (red curve) and for the BSSN data (green curve, using $N_{z}=100$ and $N_{z}=200$). This shows clearly that the boundary waves for both methods are present well before the main wave hits the boundary. It also shows that the amplitude for the BSSN data is much larger than for the lattice data. Note also that the boundary waves do not propagate very far into the grid (in stark contrast to the ADM Brill waves). By $t=10$ the main wave has left the smaller grid and the data in the left panel describes a mix of waves dominated by the reflected waves. This figure also shows that the BSSN data contains a long wavelength mode while the waves in the lattice data are much smaller in amplitude and are dominated by high frequency modes (which are rapidly suppressed by the numerical dissipation).

The evolution of the Hamiltonian constraint (7.14) is shown in figure (16). The linear growth in the constraint for the BSSN data is due solely to the growth of the BSSN bump at the origin. The sharp rise in the constraint for the lattice data for $N_{z}=201$ is due to the onset of a small instability in the lattice near the origin. This can also be seen in the small bump in the lower right plot of figure (18). This instability can be suppressed by increasing the dissipation parameter but at the expense of compromising the quality of the evolution. The source of this instability is thought to be due to the residual extrinsic curvatures driving the lattice vertices in different directions leading to distorted computational cells that break the near-planar assumptions built into the derivation of equations (.16). This is an important issue for the viability of the lattice method and will be explored in more detail in subsequent work.

10 Discussion

The passage of the waves through the outer boundaries appear to be better handled by the lattice method than both the ADM and BSSN methods. This is particularly true for the Brill waves but less so for the Teukolsky waves. It is reasonable to ask if this is a generic feature of the lattice method and if so, then which features of the lattice method gives rise to this result? An argument can be made that this behaviour may well be germane to the lattice method. The basis of the argument is the simple observation that in any small region of space-time covered by Riemann normal coordinates the first order coupled evolution equations for the Riemann curvatures (4.13–4.26) can be de-coupled to second order equations in which the principle part is the wave operatorgggThis is shown in detail in section 4.3 and 4.4 of [brewin:2010-03] but note that the author failed to explicitly state that all computations were for the principle part of the equations.. That is, for each Riemann component such as $R_{xyxy}$,

[TABLE]

where the term $\mathcal{O}\left(R^{2}\right)$ is a collection of terms quadratic in the $R_{\alpha\beta\mu\nu}$. The natural outgoing boundary condition for this wave equation is the Sommerfeld condition as per equation (6.33). Thus it is not surprising that the lattice method works as well as it does. This result is a direct consequence of the use of Riemann normal coordinates. In a generic set of coordinates the principle part would not be the wave operator.

As encouraging as the results may appear to be there remain many questions about the method. How does it behave for long term integrations? What are its stability properties? How can it be extended to higher order methods? How can mesh refinement be implemented? How well does it work on purely tetrahedral meshes? How well does it work for non-unit lapse functions? How can black holes be incorporated into a lattice (punctures or trapped surfaces?) and how would these holes move through the lattice? How can energy flux, ADM mass and other asymptotic quantities be computed on a lattice?

These are all important question and must answered before the lattice method can be considered for serious work in computational general relativity. These questions will be addressed in later papers.

The transition matrices

The transition matrices play a central role in the computation of the derivatives such as $R_{xyxy,z}$. They are used to import data from neighbouring cells so that the vertices of a chosen cell are populated with data expressed in the frame of that cell. A finite difference estimate can then be made for the required partial derivatives.

The purpose of this appendix is to extend the approach given in [brewin:2014-01]. In that paper particular attention was paid to the form of the transition matrix for a cubic lattice. It was argued that, with sufficient refinement of the lattice, the transition matrices should vary smoothly across the lattice and should converge to the identity matrix in the continuum limithhhBoth of these conditions apply to cubic lattices but need not apply for other lattices..

The particular feature of the cubic lattice that makes it attractive for our purposes is that it is easily sub-divided in a manner that preserves its original structure. This allows a whole family of cubic lattices to be constructed, with arbitrarily small cells, and thus it is easy to investigate the continuum limit of the lattice.

For a vertex $p$ with neighbour $q$ the transition matrixiiiThere is one such matrix for each pair $(p,q)$. In this paper the transition matrix will be denoted by $M$ rather than $M(p,q)$ as used in [brewin:2014-01]. $M$ allows data such as $v^{\alpha}_{q\overline{q}}$ to be imported from $\overline{q}$ to $\overline{p}$ via

[TABLE]

When constructing a frame within a cell there is considerable freedom in locating the origin and orientation of the coordinate axes. A simple and natural choice is to locate the origin on the central vertex and to align the coordinate axes with various sub-spaces of the cell (e.g., align the $x$-axis to the leg $(0,1)$, the $y$-axis to the plane spanned by the legs $(0,1)$ and $(0,2)$ etc.).

Without further information about the relationship of one cell to another little can be said about the corresponding transition matrices. However, for the cubic lattice it is not hard to see that the frames for a typical pair of cells can be chosen so that the transition matrix will be of the form

[TABLE]

where $m^{\alpha}{}_{\beta}=\mathcal{O}\left(L\right)$ are determined from the data in the pair of cells (i.e., the coordinates and leg-lengths). This form of $M$ ensures that it converges to the identity matrix in the continuum limit (e.g., by successive refinements of the cubic lattice). Note that the $m^{\alpha}{}_{\beta}$ must be subject to a constraint since the resulting transition matrix must preserve scalar products. That is, for any pair of vectors $u$ and $v$,

[TABLE]

which leads immediately to

[TABLE]

This shows that the $m_{\alpha\beta}$ define a skew-symmetric $4\times 4$ matrix determined by just six independent entries (corresponding to three boosts and three rotations).

The $m^{\alpha}{}_{\beta}$ were computed in [brewin:2014-01] by applying (.1) to a specially chosen set of vectors. A different approach will be taken in this paper, one that will be seen to be more in the spirt of Cartan’s method of local frames (see Appendix Cartan structure equations).

First recall that the lattice is assumed to be a discrete approximation to some possibly unknown smooth geometry. Thus it is reasonable to requite that the $m^{\alpha}{}_{\beta}$ should also be smooth functions across the lattice. This allows the $m_{\alpha\beta}$ to be expanded as a Taylor series based on the vertex $p$. That is

[TABLE]

for some set of coefficients $m_{\alpha\beta\gamma}$.

Now consider a closed path such as that defined by the four vertices $o,a,b,c$ in figure (4). Clearly

[TABLE]

where $v^{\alpha}_{pq\overline{r}}$ are defined by $v^{\alpha}_{pq\overline{r}}=x^{\alpha}_{q\overline{r}}-x^{\alpha}_{p\overline{r}}$ and $x^{\alpha}_{q\overline{r}}$ are the coordinates of vertex $q$ in the frame $\overline{r}$. However, the vector joining vertices $a$ to $b$ can also be expressed in terms of the frame $\overline{a}$. Likewise, the vector joining $b$ to $c$ can be expressed in terms of the frame $\overline{c}$. Using the transformation law given by (.1) leads to

[TABLE]

Substituting this pair of equations into (.6) leads to

[TABLE]

This construction can be applied to each of the 6 coordinate planes leading to 24 equations for the 24 unknowns $m^{\alpha}{}_{\beta\gamma}$. In the cases of a lattice that evolves continuously in time it is possible (see Appendix The time components of $m^{\alpha}{}_{\beta\gamma}$) to solve these equations for 15 of the $m^{\alpha}{}_{\beta\gamma}$ in terms of the extrinsic curvatures $K^{i}{}_{j}$ and the lapse function $N$. This leaves just 9 equations (based on the spatial coordinate planes) for the 9 remaining $m^{\alpha}{}_{\beta\gamma}$.

Though it is possible to use the above equations (.9) to directly compute the $m^{\alpha}{}_{\beta\gamma}$ doing so might introduce a systematic bias due to the asymmetric arrangement of the legs relative to the central vertex. An improved set of equations can be obtained simply by adding together the equations that would arise from each of the four tiles of figure (4) attached to the central vertex $\overline{o}$. This leads to the following set of equations

[TABLE]

Now since each $v^{\alpha}_{pq\overline{r}}=\mathcal{O}\left(L\right)$ it follows that the right hand side of (.10) is $\mathcal{O}\left(L^{2}\right)$ and thus

[TABLE]

which allows the terms $v^{\alpha}_{fh\overline{g}}$ and $v^{\alpha}_{df\overline{e}}$ on the right hand side of .10 to be replaced by their counterparts leading to

[TABLE]

Finally note that

[TABLE]

and therefore

[TABLE]

These are the equations that were used in the computer code to compute the $m^{\alpha}{}_{\beta\gamma}$.

Cartan structure equations

Equations (.4) and (.16) bear a striking similarity to the Cartan structure equationsjjjLatin indices will be used in this appendix to denote frame components (rather than spatial indices). This follows standard notation for differential forms.

[TABLE]

in which $\omega^{i}$ are the basis 1-forms, $\omega^{i}{}_{j}$ are the connection 1-forms and where the metric is given by $g=g_{ij}\omega^{i}\omega^{j}$ with $g_{ij}={\rm diag}(-1,1,1,1)$.

The purpose of this appendix is to show how equations (.4) and (.16) can be obtained from the Cartan structure equations (.1) and (.2).

To start the ball rolling, note that equations (.4) and (.1) agree upon choosing $m^{i}{}_{j}=\omega^{i}{}_{j}$. Showing that the remaining pair of equations (.16) and (.2) agree requires a bit more work. Start by integrating (.2) over the tile $R$ defined by the vertices $b,d,f,h$ in figure (4)

[TABLE]

where $\omega^{i}{}_{j}$ has been expanded as $\omega^{i}{}_{jk}\omega^{k}$. This equation can be re-written using Stoke’s theorem as

[TABLE]

The path integral on the left can be split into four pieces, one the four edges of the tile. On each edge set $\omega^{i}=dx^{i}$ where $x^{i}$ are the local Riemann normal coordinates appropriate to the edge (e.g., along the edge $(b,d)$ use the coordinates of frame $\overline{c}$). Thus

[TABLE]

where $v^{i}_{pq\overline{r}}=x^{i}_{q\overline{r}}-x^{i}_{p\overline{r}}$. The area integral on the right hand side of (.4) can be estimated to leading order in the length scale $L$ by approximating $\omega^{i}{}_{jk}$ by its value at the vertex $o$. Thus

[TABLE]

and noting that the integrand on the right is just the area 2-form for the tile leads to the estimate

[TABLE]

The integrated form of the Cartan equation (.4) can now be re-written as

[TABLE]

which agrees (apart from the Greek/Latin indices), to leading order in $L$, with (.16) provided $m^{i}{}_{jk}=\omega^{i}{}_{jk\overline{o}}$.

Source terms

A lattice would normally consist of a finite number of local frames, one for each central vertex. But there is nothing to stop the construction of a local frame at every point in the lattice. The new frames could be introduced by any rule but for a smooth lattice it is reasonable to require that the frames vary smoothly across the lattice. This will certainly be the case when the transition matrices are of the form

[TABLE]

The addition of these extra frames makes it easier to discuss differentiation on the lattice.

Consider a cell $p$ and some point $q$ within that cell. Let $v^{\alpha}$ be the components of a typical vector at $q$ expressed in the local frame of $q$, that is $v^{\alpha}_{q}=v^{\alpha}_{q\overline{q}}$. The components of the vector in the frame $\overline{p}$ would then be given by $M^{\alpha}{}_{\beta q\overline{p}}v^{\alpha}$. This allows the derivatives of $v^{\alpha}$ at $p$ and in $\overline{p}$ to be computed as follows

[TABLE]

At this point there is a slight problem with the notation. The last term on the right hand side above is a derivative of $v^{\alpha}$ formed from the raw point values of the $v^{\alpha}$. That derivative takes no account of the transition matrices and thus is not the partial derivative (indeed the partial derivative is the term on the left hand side). To emphasise this distinction the following notation will be used. Define a new derivative operator ${{\dagger}}$ bykkkBut note that mixed ${{\dagger}}$ derivatives need not commute.

[TABLE]

Then the equation (.4) can be written as

[TABLE]

where it is understood that all terms are evaluated at $p$ and in $\overline{p}$. By following a similar line of reasoning it is not hard to see that, for example,

[TABLE]

As a consistency check it is rather easy to see that applying this notation to $0=g_{\alpha\beta;\gamma}=g_{\alpha\beta,\gamma}$ leads directly to equation (.4). To see that this is so first note that $g_{\alpha\beta q\overline{q}}={\rm diag}(-1,1,1,1)$ at every vertex $q$ and thus the derivatives $g_{\alpha\beta{{\dagger}}\gamma}$ are zero everywhere. This leads immediately to equation (.4).

It should be noted that the hessian of lapse $N_{,ij}$ could be computed entirely from data within a single frame or by sharing data, such as $N_{,i}$, between neighbouring frames. In the later case some care must be taken when computing terms like $N_{{{\dagger}}x{{\dagger}}y}$ since the colon derivatives need not commute.lllFor example, $N_{,x,y}=N_{{{\dagger}}x{{\dagger}}y}-m^{y}{}_{xy}N_{,y}-m^{z}{}_{xy}N_{,z}$ while $N_{,y,x}=N_{{{\dagger}}y{{\dagger}}x}-m^{x}{}_{yx}N_{,x}-m^{z}{}_{yx}N_{,z}$ and as $N_{,x,y}=N_{,y,x}$ it follows that $N_{{{\dagger}}x{{\dagger}}y}-N_{{{\dagger}}y{{\dagger}}x}=m^{y}{}_{xy}N_{,y}+m^{z}{}_{xy}N_{,z}-m^{x}{}_{yx}N_{,x}-m^{z}{}_{yx}N_{,z}$ which in general will not be zero.

The time components of $m^{\alpha}{}_{\beta\gamma}$

In a lattice that is discrete in both space and time there would be 24 distinct $m^{\alpha}{}_{\beta\gamma}$ in each computational cell. However, in the case of a continuous time lattice with a zero shift vector at each central vertex, 15 of the 24 $m^{\alpha}{}_{\beta\gamma}$ can be expressed in terms of the lapse function $N$ and the extrinsic curvature $K_{ij}$, namely

[TABLE]

The key to this computation will be the application of (.9) to two carefully chosen tiles, in particular a time-like tile (generated by the evolution of a spatial leg) and a spatial tile (where all of the vertices lie in one Cauchy surface).

Showing that $m^{t}{}_{ij}=m^{t}{}_{ji}$

Consider a spatial tile in which all of the vertices of the tile lie within one Cauchy surface, Thus the $t$ component of the various $v^{\alpha}$ in (.9) are zero. This leads immediately to

[TABLE]

where the implied sum over $j$ includes only the spatial terms (since $v^{t}=0$). Since this equation must be true for all choices of $v^{i}_{oa\overline{o}}v^{j}_{oc\overline{o}}$ it follows that

[TABLE]

Showing that $m^{t}{}_{it}=(\log N)_{,i}$

Consider now the time-like tile generated by the leg $(oa)$ as it evolves between a pair of nearby Cauchy surfaces (as indicated by vertices $(o,a,b,c)$ in figure (4)). The two time-like edges $(oc)$ and $(ab)$ are tangent to the world-lines normal to the Cauchy surface while the space-like edges $(oa)$ and $(bc)$ are the two instances of the leg $(ab)$, one at time $t$ the other at $t+\delta t$. Since the shift vector is assumed to vanish at each central vertex, it follows that

[TABLE]

Likewise, for the spatial edges the $v^{\alpha}$ will have a zero $t$ component and thus will be of the form

[TABLE]

for some choice of $v^{i}_{oa\overline{o}}$ and $v^{i}_{cb\overline{c}}$. With this choice for the $v^{\alpha}$ and noting that $N_{ab\overline{a}}=N_{oc\overline{o}}+\mathcal{O}\left(L\right)$, the $t$ component of equation (.9) is given by

[TABLE]

Noting that $m^{t}{}_{ti}=0$ and estimating the left hand side by $N_{,io\overline{o}}v^{i}_{oa\overline{o}}\delta t$ leads to

[TABLE]

and since the $v^{i}_{oa\overline{o}}$ are arbitrary, it follows that

[TABLE]

in which it is understood that all terms are evaluated at $o$ in the frame $\overline{o}$.

Showing that $m^{t}{}_{ij}=-K_{ij}$

This computation follows on directly from the previous computation. This time our attention is on the spatial terms of equation (.9), namely

[TABLE]

Now recall that $v^{\alpha}_{pq\overline{r}}$ is defined by $v^{\alpha}_{pq\overline{r}}=x^{\alpha}_{q\overline{r}}-x^{\alpha}_{p\overline{r}}$ and as $x^{\alpha}_{o\overline{o}}=x^{\alpha}_{c\overline{c}}=0$ it follows that

[TABLE]

and on taking a limit as $\delta t\rightarrow 0$ leads immediately to the evolution equations

[TABLE]

for the coordinates $x^{i}_{a\overline{o}}(t)$. Now take $d/dt$ of $g_{ij}x^{i}_{a\overline{o}}x^{j}_{a\overline{o}}$ and use equation (4.5) to obtain

[TABLE]

which when combined with the above result leads to

[TABLE]

and as the first term on right vanishes due to $m_{ijt}=-m_{jit}$ the above can be further simplified to

[TABLE]

But from (.5), $m_{itj}=-m_{tij}=-m_{tji}=m_{jti}$, and as the $x^{i}_{a\overline{o}}$ are arbitrary (since the vertex $a$ can be chosen anywhere in the cell) the previous equation can only be true provided

[TABLE]

or equally

[TABLE]

Showing that $m^{i}{}_{jt}=0$

The next task is to show that $m^{i}{}_{jt}=0$. This is rather easy to do. Having just shown that $m^{t}{}_{ij}=-K_{ij}$ means that equation (.15) can also be written as

[TABLE]

which when compared with (4.5) shows that

[TABLE]

for any choice of $x^{j}_{a\overline{o}}N_{o}$. This in turn requires $m^{i}{}_{jt}=0$.

Evolution of $x^{i}$

Our aim here is to obtain evolution equations for the spatial coordinates $x^{i}(t)$ of each vertex in a computational cell.

To begin, consider two points $p$ and $q$ chosen arbitrarily in a typical cell. Equation (.15) can be applied to this pair of points leading to

[TABLE]

Now combine this pair by contracting (.1) with $x^{j}_{q\overline{o}}$ and (.2) with $x^{j}_{p\overline{o}}$ while noting that $m_{ijt}=-m_{jit}$ to obtain

[TABLE]

After shuffling terms across the equals sign this can also be re-written as

[TABLE]

This equation must be true for all choices of $(p,q)$. As the bracketed term on the left hand side depends only on $p$, that term must match the only $p$ dependent term on the right hand side, namely the $x^{j}_{p\overline{o}}$. Thus it follows that

[TABLE]

for some scalar $\alpha$. But upon setting $p=q$ in (.4) it follows that

[TABLE]

which when applied to (.5) leads to

[TABLE]

and thus $\alpha=0$. This leads immediately to

[TABLE]

with a similar result for the point $q$. Since the point $p$ is arbitrary it follow that this result holds for any point in the computational cell.

Evolution of $L_{oq}$

Equation (4.3) can be obtained from (4.1) as follows. Let $(o,q)$ be a typical leg connected to the central vertex of some cell. Our first step is to express the various vectors at $o$ and $q$ in terms of the local frames $\overline{o}$ and $\overline{q}$. Since the shift vector is assumed to be zero across the lattice it is follows that the unit normals take the simple form

[TABLE]

while

[TABLE]

which follows directly from the definition of Riemann normal coordinates $x^{\alpha}$. Recall that $x^{\alpha}_{a\overline{b}}$ are the Riemann normal coordinates of the vertex $a$ in the frame $\overline{b}$. Note also that the forward pointing unit tangent vectors $v^{\alpha}_{o\overline{o}}$ and $v^{\alpha}_{q\overline{q}}$ are given by

[TABLE]

Now substitute the above equations (.1–.6) into (4.1) to obtain

[TABLE]

where $t$ is the Riemann normal time coordinate. However, as shown in [brewin:2010-03],

[TABLE]

which using (.3–.4) can also be written as

[TABLE]

and thus

[TABLE]

which leads immediately to equation (4.3).

Complete evolution equations

The following are the complete set of evolution equations for the 14 Riemann curvatures for the particular case of a zero shift vector. These were obtained by applying the process outlined in appendix (Source terms) to the second Bianchi identities (4.13–4.26).

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]