{'fraction_non_alphanumeric': 0.03289425270439326, 'fraction_numerical': 0.012804474207079255, 'mean_word_length': 4.07070895522388, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An Einstein-Bianchi system for Smooth Lattice General Relativity. I. The Schwarzschild spacetime

The second Bianchi identity can be recast as an evolution equation for the Riemann curvatures. Here we will report on such a system for a vacuum static spherically symmetric spacetime. This is the first of two papers. In the following paper we will extend the ideas developed here to general vacuum spacetimes. In this paper we will demonstrate our ideas on a Schwarzschild spacetime and give detailed numerical results. For suitable choices of lapse function we find that the system gives excellent results with long term stability.


Introduction
Despite a slow start, hyperbolic formulations of the Einstein equations have in recent times become the system of choice for numerical relativity.
The confidence afforded to hyperbolic systems is borne out not just by the recent success in numerical relativity [1,2,3] but also from their strict mathematical underpinnings (which guarantees that future evolutions exist and that they depend smoothly on the initial data, these are key aspects of the theorems that demonstrate the stability of the system, see [4] for details).
One of the earlier hyperbolic formulations was given by Friedrich [5] in which he used the second Bianchi identities to evolve the Weyl curvatures in-situ with the metric. This idea has been extended by many other authors [6,7,8,9] and the resulting equations are commonly referred to as an Einstein-Bianchi system.
Yet despite their mathematical elegance and the virtues that this would bestow upon a numerical code there seems to be have been very few numerical applications employing an Einstein-Bianchi system (though see [10,11]).
In this paper we will report on a simple Einstein-Bianchi system adapted to a discrete lattice for static spherically symmetric spacetimes. We were lead to this formulation not by way of Friedrich's paper but rather as a direct extension of our own ideas developed in an earlier series of papers [12,13,14,15]. In those papers we used the spatial form of the second Bianchi identities (i.e., the second Bianchi identity for the 3-metric) to compute the 3-Riemann curvatures across a Cauchy surface. This device proved to be the key element in obtaining accurate and stable evolutions of the initial data.
Our longer term intention is to employ an Einstein-Bianchi system to evolve a 3-dimensional lattice. This will require not only evolution equations for the legs of the lattice, such as those given in [13], but also evolution equations for the curvatures. This will be the subject of the second paper in this series.
For the simple case presented here we find that the system works very well. The evolutions are stable, though this depends on the choice of the lapse function, see sections (4.1) and (6.1). We also find that the constraints are well behaved (they appear to grow linearly with time and converge to zero as the lattice is refined, see section (6.2)).
As this paper borrows heavily from two of our previous papers, which we refer to as Paper 1 and Paper 2, we will skip over many of the derivations and arguments assuming instead that the reader is familiar with the material in Paper 1 and Paper 2.

Spherically symmetric spacetimes
In this paper we will be constructing lattice approximations to the Schwarzschild spacetime in various slicings. In each case the continuum metric can be written in the form ds 2 = −N (r, t) 2 dt 2 + A(r, t) 2 dr 2 + B(r, t) 2 dΩ 2 (2.1) for some set of functions N (r, t), A(r, t), B(r, t) and where dΩ 2 = dθ 2 + sin 2 θdφ 2 is the metric of the unit 2-sphere. We have introduced this coordinate form of the continuum metric simply as a precursor to the introduction of the lattice. As we shall soon see, we will use the coordinate lines and their local tangent vectors as a scaffold on which to build the lattice, after which we will have no further need for the coordinates (indeed we could dispense with the coordinates altogether at the possible expense of the clarity of the exposition).
Consider a local orthonormal tetrad built from the future pointing unit normal t µ to a typical Cauchy surface and three unit vectors m µ x , m µ y and m µ z where m µ z is parallel to the radial axis (see figure (2)). These basis vectors are also tangent vectors to the coordinate axes of a local Riemann normal frame. We will use this tetrad to record the frame components of the extrinsic and Riemann curvatures on the lattice. Our notation, which we borrow from Paper 1, will be to use script characters to denote frame quantities, Also, to avoid an overflow of symbols, we will allow L xx and L zz to represent both the length of the corresponding leg as well the leg itself.
In this class of spacetimes, and on this tetrad, we know that the extrinsic curvature is diagonal and that a basis for the non-trivial Riemann curvatures is given by while the non-trivial vacuum Einstein equations yield Combining the above shows that we can express all of the non-trivial Riemann curvatures solely in terms of R xyxy and R xzxz , namely In obtaining these relations we used g µν = −t µ t ν + m µ x m ν x + m µ y m ν y + m µ z m ν z . Note that R xyxy and R xzxz are not independent for the simple equation R = 0 leads to 0 = R xyxy + 2R xzxz (2.11) We will use this equation as a check on our numerical integrations (see section (6.2) for more details).

The lattice
The symmetries in the Schwarzschild spacetime allows us to use a very simple ladder-like structure for the lattice, as indicated in figure (1). One way to imagine the construction of the lattice is to consider the coordinate mesh generated by setting t = constant and θ = π/2 in the coordinate form of the metric in (2.1). Then the rungs of the ladder are generated by small increments in φ leading to L xx ≈ B∆φ while the side rails would coincide with two radial curves (i.e., φ = constant) with L zz ≈ A∆r. Clearly, specifying all of the L xx and L zz is equivalent to specifying the metric components A(r, t) and B(r, t). Note that throughout this paper we treat the L xx and L zz as continuous functions of time.
We will label the nodes from 0 to n and on the few occasions where we need to discuss more than one leg at a time we will write (L xx ) i to denote an L xx leg at node i. In the same way (L zz ) i will denote the L zz that joins the nodes i and i + 1. Similar notation will be used for other data on the lattice.
The initial data (as described in section (5)) are constructed in a way that guarantees reflection symmetry at the throat (which is always tied to node 0).
In our computer code we extend our lattice a small way over the throat, by including the nodes -3 to -1, so that we can readily impose the reflection symmetries (by simply copying data across the throat, at no point do we independently evolve any of the data to the left of the throat).

The evolution equations
Our present task is to develop evolution equations for the leg-lengths, the extrinsic curvatures and, the principle innovation in this paper, evolution equations for the Riemann curvatures.
A simple derivation of the evolution equations for our lattice can be obtained from a general pair of equations developed in Paper 2. There it was shown that the first and second variations of arc lengths can be written in a form remarkably similar to the ADM equations, namely Note that in the following we will ignore the leading error terms O (L 3 ). Applying these equations to the two legs L xx and L zz of our spherically symmetric lattice leads immediately to The last part of the picture is to provide evolution equations for the Riemann curvatures, R xyxy and R xzxz . The basic idea is to rearrange the terms in the Bianchi identities to isolate the time derivatives while estimating the spatial derivatives from data imported from neighbouring cells. The calculations are straight-forward but a bit tedious to present here so we defer the full details to the Appendix. This leads to the following evolution equations The Riemann curvatures R xyxy and R xzxz would normally not be evolved but rather derived from the lattice data such as the leg lengths L xx and L zz . In Paper 1 we used (discrete versions of) the geodesic deviation equation and the spatial Bianchi identity * to compute the 3-dimensional Riemann curvatures 3 R xyxy and 3 R xzxz on the lattice. In raising the R xyxy and R xzxz to dynamical variables on the lattice we are forced to view equations (4.9,4.10) as constraints on the lattice data.
In section (6.2) we shall present discretised versions of these constraints which we will later use to check the quality of our numerical results.
The one remaining constraint is the standard momentum constraint (see Paper 1 for details)

Artificial viscosity
Our numerical experiments (which we will present shortly) showed that the future evolutions can be subject to high-frequency instabilities. This was seen to occur only in the cases where the lapse function was controlled by its own evolution equation (e.g., as in Harmonic slicing). For such cases we found that stability could be recovered with the addition of an artificial viscosity term to the evolution equations.
Let W be any one of the dynamical variables, L xx , L zz , K xx , K zz . Then the artificial viscosity is introduced by the addition of a simple dissipation term to the evolution equation for W . After some experimentation we settled on the following form where dW i /dt is the right hand side of the original evolution equation (4.3-4.6) and µ is a (small) constant. Other choices were tried but this form seemed to produce stable evolutions for the longest periods of time. Note that we do not add the dissipation terms to the evolution equations for the curvatures (doing so seemed to make no difference to the evolutions and had no effect in controlling the instabilities).
How should µ be chosen? We need to choose it large enough to ensure that the evolution is stable over a given time interval while also keeping it sufficiently small so as to not effect the large scale features of the numerical solution. By trial and error we found that setting µ = 0.08 worked well for evolutions to t = 100.0 using n = 2048 nodes. We also found that as the number of nodes was increased we had to make a proportionate increase in µ to maintain the same quality of the evolution over the same time interval.
The dissipation term is easily seen to be a finite difference approximation to µL zz d 2 W/dz 2 and thus it may appear to be like a Kreiss-Oliger term that vanishes in the continuum limit. However, since we are forced to set µ = O (n) and as L zz = O (1/n) we see that the term µL zz is approximately constant, say µ , and thus the dissipation term is actually of the form µ d 2 W/dz 2 . This is a standard dissipation term commonly used in hydrodynamic simulations and it does not vanish in the continuum limit.

Initial data
The initial data on the lattice are the L xx , L zz , K xx , K zz , R xyxy and R xzxz at each node of the lattice. Their time symmetric initial values were set by a combination of the Hamiltonian constraint, the geodesic deviation equation and the Bianchi identities. A full account of the choices made in coming to the equations described below can be found Paper 1. Here we will just quote the relevant equations simply to provide explicit details of how we constructed our initial data.
To ensure that the initial data is time symmetric we set K xx = 0 and K zz = 0.
The L zz were set according to the method of Bernstein, Hobill and Smarr [16] using n = 2048 on a grid of length 800.0m. The ADM mass, m, was set to be 1.0 and the L xx , R xyxy and R xzxz , for i = 1, 2, 3 . . . n, were set according At the reflection symmetric throat (i.e., at node 0) we set (L xx ) 0 = 0.01 and (R xyxy ) 0 = −2(R xzxz ) 0 = 0.25.

Results
In all of our results we used a 4-th order Runge-Kutta integrator with the time step set equal to 1/2 the smallest L zz on the lattice (which happens to be (L zz ) 0 ).

Slicing conditions
We ran our code for eight distinct slicing conditions, some were set by simple algebraic expressions while others involved differential operators.
We made four choices for the algebraic slicings, and three choice for the differential slicings, The eighth slicing condition was the simple case of geodesic slicing N = 1.
The algebraic slicings were introduced after our early explorations with the differential lapses, all of which developed high-frequency instabilities after a short time (well before t = 100). The algebraic slicings did not require any artificial viscosity and performed remarkably well, showing no signs of instabilities to at least t = 1000 (excluding the lapse (6.2) which hits the singularity at t ≈ 32). We have not run our codes beyond t = 1000 so we can not comment its stability for t > 1000.

Code tests and results
We subjected our code to many of the tests used in Paper 1, such as the time at which geodesic slicing hits the singularity, the rate at which the lapse at the throat collapses in maximal slicing and the constancy of L xx on the horizon. The results for these various slicings are shown in figures (4-9). All of the results are as expected. For the geodesic slicing the code crashes at approximately one time step short of the singularity. The familiar exponential collapse of the lapse for maximal slicing is evident in figure (6). In this case it is known that the lapse at the throat should behave as N ∼ β exp(αt) for t → ∞ with α = −(2/3) (3/2) ≈ −0.54433, see [17]. We estimated the slope of ln N vs t from our numerical data to be −0.54215 which agrees with the exact value to within 0.4 percent.
We also have a new test obtained by a simple combination of the evolution equations. From equations (4.7,4.8) we find that and as 0 = R xyxy + 2R xzxz on the initial slice (by construction, see (5.3)) we conclude that 0 = R xyxy + 2R xzxz for all time. This is not surprising, our evolution equations for the curvatures are based on the Bianchi identities and these are guaranteed to preserve the constraints. If we now set 0 = R xyxy + 2R xzxz in (4.8) and combine the result with (4.3) we find This gives us a new test of our code, that the quantity L 3 xx R xyxy should be constant throughout the evolution. Importantly this applies to all slicing conditions. In figure (9) we have plotted the fractional variations in L 3 xx R xyxy for two choices of slicings. We see that the errors for the 1 + log slicing are much larger than those for the algebraic slicing which we attribute to the use of an artificial viscosity. This last claim is easily checked by varying the artificial viscosity parameter µ. We find that the errors in L 3 xx R xyxy varies linearly with µ. Note that in obtaining equation (6.9) we have ignored the higher order error terms that would arise if we had carried through the O (L 3 ) truncation error from (4.1). Thus even if we set µ = 0 we can expect some variation of L 3 xx R xyxy over time (though this variation should vanish more rapidly than O (L 3 )).
We also have three constraint equations, namely the geodesic deviation equation (4.9), the 3-dimensional Bianchi identity (4.10) and the momentum constraint (4.11). The discrete form of these equations are whereR xzxz is the average of R xzxz across L zz while D/Dz andD/Dz are discrete derivative operators defined as follows. For a typical smooth function f (z) sampled at the grid points z i we define D f where we have introduced the superscripts +, o andto denote quantities at the grid points z i+1 , z i and z i−1 respectively. Note that the sample points z i are constructed from the lattice L zz by the recurrence relation z i+1 = z i + (L zz ) i with z 0 = 0. In this notation we haveR xyxy := (R + xyxy + R o xyxy )/2. Finally we note that the 3-curvatures can be computed from the 4-curvatures by way of the Gauss equation, Ideally we would like to see P = Q = M = 0 but in reality we expect P i , Q i and M i to be non-zero but small. This is indeed what we observe, see figure (9). We also computed a crude estimate of the rate of convergence (of Q, P and M to zero at a fixed time) by running our code twice, once with n = 2048 and once with n = 1024 and then forming suitable ratios of the constraints at the horizon. In this manner we estimated, in the absence of artificial viscosity, that P We also tried settingD/Dz := D/Dz andR xyxy := R o xyxy in the discrete Bianchi constraint but this lead to a reduction in the rate of convergence. The form of the discrete Bianchi constraint as given above (6.11) is readily seen [14] to be a second-order accurate estimate to the continuum Bianchi identity at the centre of the leg L zz .
One might ask why we have not included the Hamiltonian constraint in our code tests. The simple answer is that it is trivially satisfied by our discrete equations. This follows from the discussion surrounding equation (6.8) where we showed that 0 = R xyxy + 2R xzxz for all time. It follows that the Hamiltonian H := G µν t µ t µ will also vanish for all time. Note that this analysis was based on our discrete equations, not on the continuum equations. We did indeed check that our code maintained 0 = R xyxy +2R xzxz throughout the evolution.

A Bianchi identities
Here we will use the Bianchi identities to obtain evolution equations for the two curvatures R xyxy and R xzxz . We will follow the method given in our earlier paper [14] in which we used data imported from the neighbouring computational cells to estimate (by a finite difference approximation) the various derivatives required in the Bianchi identities. We will employ Riemann normal coordinates † , one for each computational cell, with the origin centred on the central vertex and the coordinate axes aligned with those described in section (2), see also figure (2). In these coordinates, the metric in a typical computational cell is given by where L is a typical length scale for the computational cell and g µν and R µανβ are constant throughout the computational cell. A convenient choice for g µν is diag(−1, 1, 1, 1) (such a choice can always be made by suitable gauge transformations within the class of Riemann normal frames). In this case the frame components R xyxy and R xzxz reduce to the coordinate components R xyxy and R xzxz respectively. A further advantage of using Riemann normal coordinates is that at the origin, where the connection vanishes, covariant derivatives reduce to partial derivatives.
The two Bianchi identities that we need are This pair of equations contains 4 spatial derivatives each of which we will estimate by a finite difference approximation. But in order to do so we must first have a sampling of the 4 curvatures at a cluster of points near and around the central vertex. Our simple ladder-like lattice, with its collection of computational cells along one radial axis, would allow us to compute only the z partial derivatives. For the x and y derivatives we will need to extend the lattice along the x and y axes. In short we need a truly 3 dimensional lattice. Fortunately this is rather easy to do for this spacetime. We can use the spherical symmetry of the Schwarzschild spacetime to clone copies of the ladder (by spherical rotations) so that a typical central vertex of the parent ladder-lattice becomes surrounded by 4 copies of itself. It has two further nearby vertices, fore and aft along the radial axis, that are themselves central vertices of neighbouring cells in the original ladder-like lattice. In figure (3) we display an xz slice of the cloned lattice.
We now need the coordinates of all six of the neighbouring vertices. This would require a solution of for the x µ i for given values for the L ij and R µανβ . However, as we are only going to use these coordinates to construct transformation matrices which will in turn multiply the Riemann curvatures, it is sufficient to solve (A.3) using a flat metric. Note that the above equations can only be used to compute (in fact estimate) the spatial coordinates of the vertices. For the time coordinates we can appeal to the smoothness of the underlying metric ‡ to argue that for each vertex t = O (L 2 ). The result is that the typical central vertex, with coordinates (0, 0, 0, 0), will have 6 neighbouring central vertices with coordinates as per Table (1  This accounts for the structure of our lattice but what values should we assign to the curvatures at the newly created vertices? Let (A) P Q denote the value of a quantity A at the vertex P in the local Riemann normal frame for vertex Q. Since our spacetime is spherically symmetric we can assert that Then the idea that we will import data from neighbouring cells can be expressed as (A) P Q = (U ) P Q (A P P ) where (U ) P Q is the transformation matrix, evaluated at P , from the Riemann normal frame of P to that of Q. This matrix will be composed of spatial rotations and boosts.
To get the correct estimates for the first partial derivatives we need only compute U to terms linear in the leg-lengths.
As an example, let us suppose we wished to compute v µ ,x for a spherically symmetric vector field v on the lattice. We start with (v) 10 = (U ) 10  The columns in the above matrices are labelled (t, x, y, z) from left to right and likewise for the rows. As we will latter be forming products of these matrices with the curvatures it is sufficient to compute these matrices as if we were working in flat spacetime. Thus to leading order in the lattice spacing we find § and thus In a similar manner we find For the rotations we use standard Euclidian trigonometry, for the boost we use the definition n µ i − n ν j = −K µ ν (x ν i − x ν j ) where n µ a is the future pointing unit normal to the Cauchy surface at the point a.
For the remaining two matrices, (U ) 50 and (U ) 60 , the job is quite simple, these matrices are built solely on boosts. This leads to Returning now to the construction of (v) 10 , we have We are now in a position to finally compute (v t ,x ) 00 , to wit Here we have written the truncation errors as O (L a ) with a > 0 for it is not clear, at this level of analysis, what the exact nature of this term is (save that it vanishes as L → 0). Since our spacetime is spherically symmetric we have and thus Similar calculations can be used to compute all of the spatial derivatives of v µ at the central vertex.
We can now return to the principle objective of this section -to compute the various partial derivatives of the curvatures. We proceed exactly as above but with a minor change in that we will no longer carry the truncation errors within the calculations. Thus we have for i = 1, 2, 3, 4, 5, 6 and (U µ ν ) i0 = g µρ g ντ (U ρ τ ) i0 with g µν = diag(−1, 1, 1, 1). And, as before, due to spherical symmetry.
Using the above expressions for the (U ) i0 and the following finite difference approximations We have dropped the 00 subscript as we no longer need to distinguish between the neighbouring frames. By spherical symmetry we have while from the vacuum Einstein equations we have Combining the last few equations leads to Substituting these into the above equations (A.4-A.7) and subsequently into the previous expressions for the Bianchi identities (A.1,A.2) leads to the following pair of equations Our job is almost complete, but we still have two tasks ahead of us i) to introduce a lapse function and ii) to account for the limited time interval over which a single Riemann normal frame can be used. The first task is rather easy, we simply make the coordinate substitution t → N t leading to and where we now have (g µν ) o = diag (−N 2 , 1, 1, 1). The lapse N can be freely chosen at each vertex of the lattice (but subject to the obvious constraint that N > 0). The second task is a bit more involved. We know that each Riemann normal frame is limited in both space and time. Thus no single Riemann normal frame can be used to track the evolution for an extended period of time. We will have no choice but to jump periodically to a new frame. This can be elegantly handled in the moving frame formalism. Thus our task reduces to finding a new set of evolution equations for the frame components R xyxy and R xzxz based on the equations given above for R xyxy and R xzxz .
Let e µ a , a = t, x, y, z be an orthonormal tetrad ¶ , tied to the worldline of the central vertex and aligned to the coordinate axes. Thus we have e µ t as the future pointing tangent vector to the worldline while e µ z points along the z-axis. Then Since our spacetime is spherically symmetric it is not hard to see that the tetrads of two consecutive cells (on the vertex worldline) are related by a boost in the t − z plane (arising from gradients in the lapse function). A simple calculation shows that In our frame we have chosen (g µν ) o = diag(−N 2 , 1, 1, 1), e µ a = δ µ a for a = x, y, z and e µ t = 1/N , thus we see that the last term in the previous equation is proportional to R txxz . But for the Schwarzschild spacetime we know that R txxz = 0 and thus we have     : The Schwarzschild areal coordinate r and the lattice L xx at the throat are related by r(t) = 2mL xx (t)/L xx (0). In geodesic slicing r(t) at the throat is described by the parametric equations r(t) = m(1 + cos η(t)), t(η) = m(η + sin η). These equations allow us to plot the exact evolution of r(t) (the smooth curve) against estimates from the lattice (solid points). The relative errors are seen to be very small and are dominated by the truncation errors in the Runge-Kutta scheme.   1.00e-03 Figure 6: As per figure (5) but this time we display the lapse profiles. Note how quickly the lapse collapse at the throat in the 1 + log slicing. This would likely cause serious underflow problems for t 300. Notice also the uniform spacing of the curves along the logarithmic axis for the maximal lapse. This show that the lapse collapses exponentially at the throat (a well known result for maximal slicing). The left pair of figures shows the evolution over the range t = 0 to 100 in steps of 10, while the right pair cover the range t = 0 to 1000 in steps of 100. This lapse does not appear (on this time scale) to exhibit an exponential collapse at the throat. The profiles for the algebraic slicings appear to propagate into the asymptotically flat regions far more rapidly than any of the differential slicings. This may be explain why an artificial viscosity was not needed for the algebraic slicings -they carry away any small numerical errors before they have chance to grow.  : According to equation (6.9) the quantity L 3 xx R xyxy should be conserved. In the top row we display the relative error, defined by 1 − C(t)/C(0) with C(t) = L 3 xx (t)R xyxy (t), for two choices of lapse. The errors for the 1 + log slicing are much larger than those for the algebraic slicing which we attribute to the action of the artificial viscosity terms. In the bottom row we display the momentum constraint for the same pair of slicings. This shows a slow growth in the momentum constraint over time (judging by the peaks in the plots, the growth appears to be linear in time).