{'fraction_non_alphanumeric': 0.035256698772766826, 'fraction_numerical': 0.013292525579860173, 'mean_word_length': 4.146386658431131, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, '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}}
Numerical stability of the AA evolution system compared to the ADM and BSSN systems

We explore the numerical stability properties of an evolution system suggested by Alekseenko and Arnold. We examine its behavior on a set of standardized testbeds, and we evolve a single black hole with different gauges. Based on a comparison with two other evolution systems with well-known properties, we discuss some of the strengths and limitations of such simple tests in predicting numerical stability in general.


I. INTRODUCTION
In recent years the quest for finding a numerically stable formulation of the Einstein evolution equations has become more and more intense, see e.g. [1,2] for reviews. Effort has been put into finding first order symmetric hyperbolic formulations of the evolution equations, since the properties of such systems can be analyzed mathematically. However, even if an evolution system is symmetric hyperbolic there is no guarantee that its numerical implementation is stable when evolving a highly dynamic black hole spacetime. The Baumgarte-Shapiro-Shibata-Nakamura (BSSN) [3,4] system is an example of a system that is not first order symmetric hyperbolic but has nice stability properties. Thus, current mathematical analysis is not sufficient to explore the properties of an evolution system. The system must be implemented and tested numerically before we can draw definite conclusions about its viability for a particular physical application.
In [5], Alekseenko and Arnold (AA) suggest a first order formulation of the evolution equations that is symmetric hyperbolic when considering only a subset of the variables. In particular, the metric itself and some of its first spatial derivatives are not considered part of the evolution system and are treated as given functions when showing hyperbolicity. It is argued that this is sensible since the metric is derivable from an ordinary differential equation. A distinguishing feature of the AA system is that a minimal number of first derivatives of the metric are introduced as independent variables. The system has only 20 unknowns and no parameters that have to be fixed for hyperbolicity, so it is relatively simple for a symmetric hyperbolic system.
We have chosen to implement the AA system numerically and to compare its numerical properties with those of the Arnowitt-Deser-Misner (ADM) [6,7] and BSSN systems. The ADM system is known to be unstable in many situations, but we include it here since for finite time intervals of evolutions it has been used rather successfully in practice, e.g. in 3D black hole simulations [8,9,10,11,12,13,14,15]. The BSSN system can be obtained from a trace-conformal decomposition of the ADM equations, and with suitable techniques it is very stable for black hole evolutions. For example, the first stable evolution for more than 100000M of a single Schwarzschild black hole in a (3+1)-dimensional code without adapted coordinates was obtained with a modified BSSN system [16]. Both the ADM and BSSN systems are not first order. However, there exist first order versions of the BSSN system which are symmetric hyperbolic [2,17,18,19], if the densitized lapse and shift are considered given functions. Straightforward first order forms of ADM are only weakly hyperbolic, which implies certain numerical instabilities (see e.g. [20]).
The second order version of the BSSN system shares with the AA system the property that a subsystem of it is indeed symmetric hyperbolic. A crucial step in the construction of BSSN is the introduction of the contracted Christoffel symbol of the conformal metric,Γ i , as an independent variable. The evolution equation for the extrinsic curvature then contains derivatives of the metric only in the form of a Laplace operator, so that the metric obeys a wave equation ifΓ i is considered a prescribed variable, ignoring its own evolution equation. This partial hyperbolicity of the BSSN system may well be a crucial ingredient in its success as evolution system. Hence the question arises whether symmetric hyperbolicity in a subsystem implies a numerical advantage for the AA system.
There is a variety of numerical tests that one can perform on an evolution system. A complex and important issue is that seemingly minor changes in the test conditions and implementation of the system can lead to very different conclusions about stability. In [21], an important step is taken towards creating a set of tests that can serve as a standardized benchmark for stability. Our discussion of the AA system in comparison to ADM and BSSN can also be viewed as a contribution to the ongoing development of such benchmarks. On the one hand, we report on the performance of our particular implementation of ADM and BSSN, where a large body of prior work allows us to judge how representative the current benchmark is. On the other hand, we apply these tests to a new evolution system about which nothing is known so far from numerical experiments, and the question is, for example, whether one can predict the usefulness of the AA system for black hole evolutions. We therefore also discuss results for the evolution of a single black hole that go beyond the current set of benchmarks in [21].
The paper is organized as follows. In Sec. II, we write out the AA system explicitly, in notation that is more familiar to numerical relativists, and we show that the complete system is not symmetric hyperbolic. We introduce a time-independent conformal rescaling of the AA system that we will need later for black hole evolutions. In Sec. III, we perform three tests from the recently proposed test suite [21] for numerical evolutions, the robust stability test, the gauge stability test, and the linear wave stability test. Finally, we evolve a single black hole space time in Sec. IV, since this is the physics example that we are most interested in. We conclude with a discussion in Sec. V.

II. THE AA SYSTEM
In [5], the evolution system is given in compact notation. Here we write it out in the form in which we have implemented it: g ij,s g pq,r g pq g rs + g pq,r g si,j g pq g rs + g pq,r g sj,i g pq g rs ) , ji − g pj,q g rs ,s g ri g pq + g pq,j g rs ,s g ri g pq − g pi,q g rs ,s g rj g pq + g pq,i g rs ,s g rj g pq + (5) g pq,r g rs ,j g si g pq + g pq,r g rs ,i g sj g pq − g pq,r g qs ,j g si g pr − g pq,r g qs ,i g sj g pr − g pj,q g qr ,s g ri g ps − g pq,j g qr ,s g ri g ps − g pi,q g qr ,s g rj g ps − g pq,i g qr ,s g rj g ps + 2g pj,q g pr ,s g ri g qs + 2g pi,q g pr ,s g rj g qs − 4g ij,s g pq,r g pr g qs + 4g ij,s g pq,r g pq g rs ) , ij − 2α ,i,j − α ,p g ij,q g pq + α ,p g qi,j g pq + α ,p g qj,i g pq + 2αKK ij + (6) 2β p ,j K pi + 2β p ,i K pj − 4αg pq K pi K qj , c 5) = − 1 4 g pq,r 2g qr ,s g ps − g pq ,s g rs + g st,u g pq g ru g st (7) ij − α g ij,s g pq,r (− (g ps g qr ) + g pq g rs ) + g ij c 5) ,q (g rj,p g ik − g ri,p g jk ) g qr + β p ,j (g ik,p − g qr,p g ik g qr ) + β p ,i (−g jk,p + g qr,p g jk g qr ) , α (g pj,q g pq K ik + g pq,j g pq K ik + g pi,q g pq K jk − g pq,i g pq K jk + g pj,q g ik K pq − g pq,j g ik K pq − g pi,q g jk K pq + g pq,i g jk K pq ) + g pj,q g pq w ik − g pq,j g pq w ik − g pi,q g pq w jk + g pq,i g pq w jk − g pj,q g ik w pq + g pq,j g ik w pq + g pi,q g jk w pq − g pq,i g jk w pq )) , αg pq ,i g jk K pq + 2α ,p g jk g pq K qi − 2α ,p g ik g pq K qj − αg pq,r g jk g pq g rs K si + αg pq,r g ik g pq g rs K sj , ,q (f pij + f pji ) + f pij,q g pq + f pji,q g pq − g pq ,r f psj g qi g rs + g pq ,r f sjp g qi g rs − g pq ,r f psi g qj g rs + g pq ,r f sip g qj g rs , The operator ∂ 0 = ∂ t − β p ∂ p is defined in terms of ordinary partial derivatives for all variables.
The new variable f ijk is anti-symmetric in i and j and also satisfies the cyclic property f ijk + f jki + f kij = 0. This means that f ijk represents 8 independent variables. The total number of 20 variables is thus smaller than in the case of most first order formulations (which usually involve 30 or more variables), but larger than for the BSSN system, which has 15 independent variables.
In addition, the equations in the AA system are more complex than in the BSSN system, so that the AA evolution system takes roughly twice as long to run the same number of iterations, even though we have made an effort to implement the AA system in a numerically optimal way. The simulations were carried out with the BAM code, which is a rewritten version of the code used in [9]. Part of BAM is a Mathematica script to convert tensor equations to C code, and the AA equations were generated this way using Mathematica and MathTensor.
A. Hyperbolicity of the AA system Looking at Eqs. (1)- (14) and assuming that α and β i are prescribed given functions, we see that the system has the form is the solution vector with spatial indices suppressed, and is a matrix which contains a symmetric submatrix In addition, Eq. (15) contains the two vector valued functions v and w, where the argument of w depends on and its transpose u T ,i , and also on another matrix Q(u), which is of the simple form so that the argument of w in fact only depends on squares of first derivatives of g ij . From Eq. (15) it is immediately apparent that the system is of first order form, as no second derivatives of u appear. Nevertheless the system differs from the standard first order form by the term w. This means that the standard theorems about well-posedness (see e.g. [1]) are not applicable, and that we cannot easily compute characteristic speeds and modes of the system. However, if we linearize the system around any background u B , all terms of the form w drop out. This can be seen as follows. Assume that where u L is a small perturbation to the background u B . Then and hence Eq. (15) becomes which is now of standard first order form, but with a modified matrix Looking at Eq. (23) and (24), it is immediately clear that the system is not symmetric hyperbolic. In fact, since the matrix S i (u B ) has several zero eigenvalues, the full system in general is unlikely to have a complete set of eigenvectors and thus to be strongly hyperbolic. The exception is linearizion around flat space wherer i ands i in Eq. (24) vanish, so that the system (23) becomes symmetric hyperbolic in this special case. Also, since S i (u) is symmetric, the subsystem of K ij and f ijk with g ij considered as a prescribed variable is symmetric hyperbolic, which holds true even for the non-linearized system (15). Alekseenko and Arnold [5] emphasize this last point and stress that the evolution equation for g ij is (as usual) just an ordinary differential equation. Yet in numerical evolutions the latter property may be unimportant, since g ij has to be evolved along with the other variables and cannot be prescribed. In addition, it is not clear if the fact that the evolution system (15) contains a symmetric hyperbolic subsystem has any bearing on the stability properties of the full system.

B. Conformal version of the AA system
In Sec. IV, we report on numerical evolutions of black hole puncture data [9,22,23]. In order to numerically evolve puncture data without excision the AA system has to be modified. The reason is that e.g. for two punctures the metric has the form where r 1 and r 2 are distances from the punctures, m 1 and m 2 are the bare puncture masses, and u is finite. From Eq. (25) it is apparent that the physical metric diverges like at each puncture, so that finite differencing calculations across or near any puncture will fail. The same problem also occurs in the variable For this reason we do not evolve the variables g ij and f ijk directly, rather we rescale g ij and f ijk by the time independent conformal factor such that the rescaled quantities become finite at the puncture. Furthermore, we also introduce the rescaled extrinsic curvaturē in order to remove divergences in K ij , which in the case of punctures with spin behaves as Since ψ is an a priori prescribed time independent function, the principal part of the system of evolution equations (1)- (14) remains unchanged if we use the rescalings in Eqs. (30), (31), and (32) to express the system in terms of the new variablesḡ ij ,K ij , andf ijk . In addition, all terms involvingḡ ij ,K ij ,f ijk , and their derivatives are finite in the new system so that finite differencing can be used without trouble. There are, however, additional terms containing spatial derivatives of ψ, which cannot be computed using finite differencing. Yet, since ψ is a known function we can use analytic expressions for its derivatives. Furthermore, all such derivative terms also contain appropriate powers of ψ which make them finite.

III. STABILITY TESTS
We have used some of the stability testbeds suggested in [21], also known as the "apples with apples" tests, to explore the properties of the AA system. We also show test results for the ADM and BSSN systems for comparison.

A. Robust stability test
The purpose of the robust stability test is to determine how an evolution system will handle random errors that are bound to occur at machine-precision level. Random constraint violating initial data in the linearized regime is used to simulate this machine error.
All apples with apples test are run on a full-3D grid, but in this case with periodic boundary conditions in all directions. The 3D domain is a 3-torus, and here we only use 3 grid-points in the y-and z-directions. The parameters are: • Number of grid points in each direction: n x = 50ρ, n y = n z = 3, with ρ = 1, 2, 4 • Courant factor = 0.5 • Iterations: 100000ρ, output every 100ρ iterations (corresponding to 1 crossing time) The initial data is given by where ε ij is a random number with a probability distribution, which is uniform in the interval (−10 −10 /ρ 2 , +10 −10 /ρ 2 ). After a small number of timesteps, the random noise in g ij will have propagated into all other evolved quantities, except for β i , which will remain identically zero for all time. Hence our initial data differ slightly form the ones in [21], who add random noise to all quantities which need initialization.
There is one obvious problem with the robust stability test. The initial data are random, and we must carefully check that the random number generator we use does not introduce errors because the numbers are not totally random. We use C's "random" function on Linux Redhat 7.3, a pseudo-random number generator based on a non-linear additive feedback algorithm that avoids the short-comings of some of the older implementations of the "rand" function. As a seed for the random number generator we use the time where the subroutine is called plus the clock cycle of the processor on which it is called. Since our code is parallelized and each processor uses its own seed, the actual random numbers on the grid depend on the number of processors. This could be avoided by additional coding, however, the qualitative result of the robust stability test is not expected to depend on the actual random numbers. We have run the test several times on different numbers of processors. We have also tried increasing the size of the domain, both in the x-direction and the y-and z-directions (which results in a different number of random numbers being generated). We find that although this does change the results, many features are robust and do not change with these variations. We will only comment on these robust features of the test.
ADM and BSSN results are similar to those reported in [21]. ADM grows exponentially (Fig. 1), while BSSN is stable (Fig. 3). Note that ADM clearly shows the signature of a grid mode, since when plotting versus the number of iterations almost identical exponential growth is obtained for the three different resolutions.
In [24], the ADM system is run with a Courant factor of only 0.03 with the result that it is stable for much longer, up to 200 crossing times. We have lowered the Courant factor in the ADM run to 0.25 and found that we do not see exponential growth in the Hamiltonian constraint. When plotting trK we see some growth, similar to Fig. 2 in [24], but we do not encounter a blowup (we ran the coarsest resolution to 10000 crossing times, and did not encounter exponential growth). The results are shown in Fig. 2. We also tried with Courant factor 0.03 as in [24], and this did not show any exponential blowup (the run was stopped at 600 crossing times). Preliminary experiments indicate that the transition between stable and unstable Courant factors for ADM lies between 0.4 and 0.5. This deserves further investigation.
The results for AA are shown in Fig. 4. AA is qualitatively stable for the 1000 crossing times suggested for this test, but oscillations in the maximum and minimum of the Hamiltonian occur. There are two types of oscillations, one with a long wavelength and one with a short wavelength. We can see that the short wavelength oscillations dampen out. The long wavelength of the medium and high resolution runs seems to be the same, while the wavelength in the low resolution run is about half of the wavelength in the other runs. As pointed out earlier, these features are robust when changing the domain-size, but the amplitude of the oscillations are affected by this change.
Given the presence of these slow oscillations, we ran another set of tests where we evolve for 10000 crossing times. Fig. 5 shows the results. We can see that at late times, the constraint violation for all 3 resolutions are growing, and this suggests that the AA system is stable for random noise initial data for a long time, but it will eventually become unstable. We also ran the test of the BSSN system to 10000 crossing times but found no instabilities or any indication of longterm growing oscillations in this case.
We tried lowering the Courant factor for the AA runs. We ran the coarsest resolution to 10000 crossing times with a Courant factor of 0.25 and found that the growth in the Hamiltonian constraint still happens, but it happens later than in the run with the Courant factor of 0.5. Lowering the Courant factor for the AA system does not change the general features of the oscillations.

B. Gauge wave stability test
In this test we look at the ability of the evolution systems to handle gauge dynamics. This is done by considering flat Minkowski space in a slicing where the 3-metric g ij is time dependent. The gauge wave is then given by g yy = g zz = 1, (39) g xy = g xz = g yz = 0, Here d is the size of the domain in the x-direction and A is the amplitude of the wave. Since this wave propagates along the x-axis and all derivatives are zero in the y-and z-directions, the problem is essentially one dimensional.
These are the parameters we use in our gauge stability tests: • Simulation domain: x ∈ [−0.5, +0.5] • Points: n x = 50ρ, n y = n z = 3, ρ = 1, 2, 4 • Courant factor = 0.25 Since we know the analytical solution at all times, we can compare our numerical results to the analytical results, see also [20]. As well as testing if the system has exponentially growing modes, we can check the convergence of the numerical result to the analytical solution with increasing resolution. Fig. 6 shows that for a gauge wave amplitude A = 0.1 the AA system converges as expected for a finite time interval, but there is exponential growth and the run crashes after about 100 crossing times. For A = 0.01 (Fig. 7) it takes longer for this non-convergence to show, about 1000 crossing times, but it is still there. If the AA system together with the given gauge equations was symmetric hyperbolic, and if one had a stable discretization scheme, then the result would be convergent. Conversely, assuming that ICN is an appropriate discretization scheme, we would conclude that the AA system in harmonic gauge is not symmetric hyperbolic, which agrees with our result in Sec. II A.
The ADM system shows no growth of the constraints at all, and the constraint violation remains at machine precision. This is probably because for a 1-dimensional gauge wave the ADM system simplifies such that no constraints violating modes are possible.
The results for the BSSN system are shown in Fig. 8. The BSSN system becomes unstable and non-convergent rather quickly (on the order of 100 crossing times). This is somewhat surprising since the BSSN system has been very successful in single black hole evolutions and neutron star evolutions [23,25], so we expect it to be able to handle gauge dynamics. However, part of the robustness of BSSN can be attributed to the fact that constraint violating modes have a finite speed [26,27] and can propagate out of the grid for example for radiative boundary conditions. However, on our grid with periodic boundaries, when a constraint violating mode appears it will stay on the grid, which is probably the reason that the BSSN system fails this test. Thus, this test by itself cannot determine whether a system can handle gauge dynamics, but must be accompanied by other tests without periodic boundaries. In particular, there is no contradiction to the observed stability of BSSN for single black holes with radiative boundaries.

C. Linear wave stability test
In this section we investigate the ability of the evolution systems to propagate the amplitude and phase of a traveling gravitational wave. The amplitude of this wave is sufficiently small so that we are in the linear regime. This test reveals effects of numerical dissipation and other sources of inaccuracy in the evolution algorithm.
Figs. 9,10, show the results for the ADM, AA, and BSSN systems, respectively. Fig. 9 is a comparison of the numerical wave to the analytical wave at 100 crossing times. We see that the AA numerical wave travels slightly faster than the ADM and BSSN numerical waves. Fig. 10 shows the L 2 -norm of the difference between the numerical and analytical linear waves as a function of time, and again the wave in the AA system travels faster than in the other systems. The Hamiltonian constraint has a value of about 10 −8 at the end of the run for the AA system, so we are still well within the linear regime. It is also worth noting that in the ADM system the constraint violation is constant throughout the evolution, while the constraints grow slightly for the other two systems (but the maximum violation is of the order 10 −8 for the entire run, for all resolutions).

IV. NUMERICAL TESTS INVOLVING A SINGLE SCHWARZSCHILD BLACK HOLE
In this Section we present numerical results for the evolution of a single Schwarzschild black hole in two different gauges.

A. Geodesic slicing
The initial data for this test is a Schwarzschild black hole in isotropic coordinates. We use geodesic slicing to evolve the initial data. In geodesic slicing, the coordinate lines correspond to freely falling observers, and the resulting slicing of Schwarzschild can be expressed analytically in terms of Novikov coordinates, e.g. [28], Chapter 2.7.2, and [29], which we transform to isotropic coordinates for direct comparison with the numerical results. To this end, we numerically solve the following implicit equation for the Schwarzschild area radial coordinate R = R(τ, R max ), where R = R max is the position at τ = 0 of a freely falling observer starting at rest, and M is the mass of the black hole. To transform from the R max radial coordinate to the isotropic radial coordinate r we use Then the rr-component of the metric in isotropic coordinates is given by where and R is computed from Eq. (50). Here ψ(r) = 1 + M 2r . Note that there is a typo in [29], Eq. (16).
Analytically, the metric becomes infinite at time which is the time for an observer that starts from rest at the horizon to reach the Schwarzschild singularity. This results in a crash in the numerical computations. We run this test on a 3D grid in the so-called cartoon mode [30], because the spherical symmetry of the problem means that we can do a computation using information only on a line passing trough the center of the black hole. The computational domain is z ∈ [0, 80M ], x = 0, y = 0. We use second order accurate finite differencing together with an iterative Crank-Nicholson scheme with a Courant factor of 0.25 for evolution. Since we never run longer than t crash = πM , the outer boundary at 80M has no effect on the black hole located at the origin. In order to compute the order of convergence, we run the test at the resolutions h = 0.01M , 2h = 0.02M , and 4h = 0.04M . We compute σ on each gridpoint present in the coarsest run, and f h is the value of the quantity under consideration (here the metric component g zz ) for resolution h. Our grid-points are not in the same places for all three resolutions, because we always stagger the puncture between two grid-points. This means that we must interpolate to get the values at all the coarse grid-points, for which we use 4th order Lagrange interpolation. Fig. 11 shows σ for the AA system at time t = π 2 M and time t = 3.0M in the region close to the black hole where the metric deviates the most from the flat conformal metric. We see that we have second order convergence close to the black hole in both cases. The spikes in the curve at later times are due to the fact that the curves of g zz for the 3 resolutions cross.

B. 1+log lapse, gamma driver shift
We have implemented 1+log lapse and gamma driver shift in our code. This gauge choice makes the BSSN system stable for more than 1000M for certain single black hole runs [23]. While the geometric motivation of this gauge choice, namely singularity avoidance and reduction of slice-stretching, should be valid for any evolution system, it is unclear whether the AA system will be as stable as BSSN or unstable in this test case. In particular, note that the gamma driver shift is designed to control the evolution of one of the BSSN variables,Γ i . We have not implemented proper outer boundary conditions for the AA system, so in our simulations we have waves coming in from the outer boundaries. We use the same parameters as for the geodesic runs, i.e. the domain in cartoon mode is z ∈ [0, 80M ], x = 0, y = 0, the resolution is h = 0.01M , and we use a Courant factor of 0.25.
Our AA run crashes at around 20M . Very steep gradients appear near the black hole and they eventually kill the run. The reason for this seems to be that the 1+log lapse depends on the trace of the physical extrinsic curvature. Since we evolve the conformal extrinsic curvature, finite difference errors are enlarged by a factor of ψ 4 when computing the physical extrinsic curvature, and these account for the differences in the lapse evolution between BSSN and AA runs. When using the AA system, the lapse drops very fast, leaving a sharp gradient between the frozen region and the region where the fields can evolve. The gradients in the physical variables that kill the run appear where the frozen region meets the dynamic region. In the BSSN system this gradient is more shallow and we do not see this effect.
This demonstrates that, as expected, a gauge choice that leads to numerically stable evolutions with one evolution system may well be unstable for others.

V. DISCUSSION
We have investigated the numerical properties of the previously untested AA evolution system by a variety of numerical experiments, which for comparison we also performed for the ADM and BSSN systems. Since one aspect of the test suite put forward in [21] is that different numerical implementations of one and the same evolution system can be compared, let us point out that our results for ADM and BSSN agree with [21] for the specific results shown there (see also [31]).
However, if the Courant factor in the robust stability test is lowered from 0.5 (as suggested in [21]) to 0.25 we find that the stability of ADM dramatically improves (compare Figs. 1 and 2), while the stability properties of AA remain unchanged. This implies that the Courant factor used in [21] is too large for the ADM system causing immediate exponential growth, instead of the initially linear growth expected for a weakly hyperbolic system. This confirms a similar observation already described in [24]. We use iterative Crank-Nicholson in all our simulations, but [24] also point out that dissipation can mask the linear growth expected for ADM. So one should repeat these experiments with a less dissipative scheme like third order Runge-Kutta to look for linear growth in ADM, but also in the AA systems, see Sec. II A.
One property of BSSN that has not been explicitly stated in [21] is its drastic failure on periodic domains. This observation is consistent with, and actually strengthens, the notion that BSSN performs well because it is able to propagate modes off the grid [26], in particular for isolated systems with radiative boundary conditions. It will be interesting to see how AA, ADM, and BSSN perform on a gauge wave test with outer boundary conditions that let the gauge wave propagate away from the domain.
Our findings for the AA system can be summarized as follows. In the robust stability test, the AA system is stable for a long time, but eventually does go unstable. In this case the AA system does not do as well as the BSSN system. For gauge waves in periodic geometry, the AA system is unstable, but the runs last much longer than the corresponding BSSN runs. The linear wave test shows that the AA system produces a larger drift in the phase than the ADM and BSSN systems, causing the linear wave to propagate faster compared to the analytical solution in the AA evolution than in the other two evolutions.
In the black hole runs, AA does as well as ADM and BSSN for geodesic slicing runs, but fails when trying to use gauges that makes BSSN runs long term stable. There are two distinct aspects of this gauge choice. While the gauge is appropriate geometrically (singularity avoiding and countering slice-stretching), there is no reason to expect numerical stability for the AA system. Since it was non-trivial to find a numerically stable gauge choice for BSSN, it is not too surprising that additional work is required to find a gauge choice for the AA system that allows long run times for a single black hole, if indeed such a choice exists. We have found that the tests published in [21] are helpful in exploring the properties of the AA system, but also that, not unexpectedly, these tests cannot by themselves determine whether an evolution system is worth exploring further for a particular physical system, say black hole evolutions. Note that one should expect that for different physical initial data sets different evolution systems are optimal, see e.g. [32]. Since the relationship between choice of evolution system, gauge choice, outer boundary conditions, and the physical properties of the problem we are trying to simulate is complicated, it is not clear that sufficient understanding of the numerics needed to evolve binary black holes can be gained by singling out for example the evolution system and ignoring the other issues involved. Ultimately, if one wants to evolve black holes, one should evolve black holes. While the AA system has had some success in our numerical experiments, tests with black holes, proper outer boundaries, and appropriate gauge choices are needed to determine if the AA system has a future in numerical relativity. FIG. 2: Robust stability test for the ADM system, with Courant factor 0.25. The legend is: solid: ρ = 1, dotted: ρ = 2, dash-dot: ρ = 4. We show L∞(ham) on the left and L∞(trK) on the right. The exponentially growing mode seen in Fig. 1 is not present.
FIG. 5: Robust stability for the AA system, run until 10000 crossing times. The legend is: solid: ρ = 1, dotted: ρ = 2, dash-dot: ρ = 4. We see that at late times the constraint violations are growing, which will probably cause a crash if we evolve long enough.
FIG. 6: Gauge stability test for the AA system with a wave amplitude A = 0.1. The legend is: solid: ρ = 1, dotted: ρ = 2, dash-dot: ρ = 4. As expected, the metric converges at early times, but there is exponential growth and the run crashes after about 100 crossing times.
FIG. 9: Linear wave stability test for (from left to right) the AA system, the BSSN system, and the ADM system. We compare the numerical wave at 100 crossing times to the analytical solution. The legend is: solid: analytical solution, dotted: ρ = 1, dash-dot: ρ = 2, dash: ρ = 4.
FIG. 10: Linear wave stability test for (from left to right) the AA system, the BSSN system, and the ADM system. We show the L2 norm of the difference between the analytical and numerical values at different crossing times. The legend is: solid: ρ = 1, dotted: ρ = 2, dash-dot: ρ = 4.
FIG. 11: Geodesic slicing of Schwarzschild. Shown is the order of convergence σ of the conformal metric componentḡzz for the AA system in the computational domain near the black hole for t = π 2 M and t = 3.0M . We observe second order convergence near the black hole.  Fig. 11). The analytical solution is shown at 3 different times: t = 1.55 (solid), t = 3.0 (dotted) and t = 3.125 (dash-dot). The other 3 panels show the absolute value of the relative difference of the numerical and analytical solutions at the 3 different times. In these plots the AA solution is solid, ADM is dotted, and BSSN is dash-dot.