{'fraction_non_alphanumeric': 0.030222544166767274, 'fraction_numerical': 0.013535085248232659, 'mean_word_length': 4.6578627808136, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 6, '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}}
Compact binary evolutions with the Z4c formulation

Numerical relativity simulations of compact binaries with the Z4c and BSSNOK formulations are compared. The Z4c formulation is advantageous in every case considered. In simulations of non-vacuum spacetimes the constraint violations due to truncation errors are between one and three orders of magnitude lower in the Z4c evolutions. Improvements are also found in the accuracy of the computed gravitational radiation. For equal-mass irrotational binary neutron star evolutions we find that the absolute errors in phase and amplitude of the waveforms can be up to a factor of four smaller. The quality of the Z4c numerical data is also demonstrated by a remarkably accurate computation of the ADM mass from surface integrals. For equal-mass non-spinning binary puncture black hole evolutions we find that the absolute errors in phase and amplitude of the waveforms can be up to a factor of two smaller. In the same evolutions we find that away from the punctures the Hamiltonian constraint violation is reduced by between one and two orders of magnitude. Furthermore, the utility of gravitational radiation controlling, constraint preserving boundary conditions for the Z4c formulation is demonstrated. The evolution of spacetimes containing a single compact object confirm earlier results in spherical symmetry. The boundary conditions avoid spurious and non-convergent effects present in high resolution runs with either formulation with a more naive boundary treatment. We conclude that Z4c is preferable to BSSNOK for the numerical solution of the 3+1 Einstein equations with the puncture gauge.


I. INTRODUCTION
This is the conclusion of a series of papers [1][2][3][4][5] about the development of a formulation of general relativity (GR), called Z4c, that attempts to combine the strengths of two popular evolution systems for applications in freeevolution numerical relativity. Here we summarize the logical development of the formulation.
In the BSSNOK formulation [6][7][8] there is a zerospeed characteristic variable in the constraint subsystem, which can result in large Hamiltonian constraint violations in numerical applications; the removal of this mode is one of the key advantages of the generalized harmonic formulation [9][10][11][12] over BSSNOK. The generalized harmonic formulation also possesses a constraint damping scheme [13], which exponentially damps away small, high-frequency constraint violations at the continuum level. Furthermore the trivial wave-like nature of the generalized harmonic subsystem allows for the convenient construction of constraint preserving boundary conditions [14][15][16]. On the other hand, the key advantages in the BSSNOK formulation are the choice of conformal variables, and the fact that the formulation does not come tied to a particular gauge, which allows for the selection of the moving puncture gauge [17][18][19][20][21][22]. The combination of these two strengths allows for the evolution of black holes represented by coordinate singularities on the grid without severe numerical difficulties. For BSSNOK radiation controlling constraint preserving boundary conditions have been proposed [23], but to our knowledge have not been successfully used in numerical applications. With these considerations in mind a conformal decompo-sition of the Z4 formulation [24][25][26][27][28][29][30] was proposed in [1]. Because of the close relationship between the Z4 and generalized harmonic formulations, this conformal decomposition inherits all of the strengths outlined in the previous discussion. In [1] a set of spherically symmetric (1D) tests involving single black hole and neutron star spacetimes demonstrated that, indeed, the Z4c formulation guarantees the robustness of, and better constraint preservation than BSSNOK, especially for non-vacuum spacetimes. It was also found that the main advantage of Z4c in the bulk of the computational domain, namely the propagating constraint violations, presents problems at the outer boundary if the constraints are not absorbed but rather reflected.
A first attempt to tackle this problem was presented in [2], in which high derivative order constraint preserving boundary conditions based on those of [16] were proposed for Z4c. Well-posedness of the Z4c constraint subsystem initial boundary value problem was demonstrated, and numerical results in explicit spherical symmetry demonstrated the efficacy of the constraint preserving boundary conditions. In this work we use constraint preserving boundary conditions, motivated by a forthcoming study [5], in which well-posedness of high order boundary conditions that are constraint preserving and control incoming gravitational radiation [31,32], is analyzed.
In [3], the performance of the Z4 constraint damping scheme of [13] applied to Z4c in black hole and neutron star spacetimes was studied in detail in spherical symmetry. The constraint damping scheme is effective, as expected, in the non-linear system provided that the con-straint violation is sufficiently small and resolved on the numerical grid; in the case of grid noise the combination of artificial dissipation and damping helps to suppress constraint violations. But it was found that care should be taken in the choice of damping parameters. Success in spherical tests does not necessarily guarantee that of 3D simulations. Preliminary simulations with Z4c and Sommerfeld BCs in 3D with pure box in box mesh refinement showed poor behavior at the outer boundaries in long evolutions. Therefore to try and bridge the gap between the evolutions in spherical symmetry and full applications in astrophysical spacetimes, numerical stability of Z4c in 3D evolutions was studied in [4], where numerical stability of the linearized Z4c system coupled to the puncture gauge with a novel discretization was shown, and numerical evidence from the Apples-with-Apples [33][34][35] tests was presented with both the standard and the novel discretization. The importance of algebraic constraint projection was highlighted, and limitations of the puncture gauge for applications in cosmology were observed.
The first application of Z4c was the study of the endstate of a collapsing neutron star with the puncture gauge [36] in 1D and 3D simulations of spherical configurations. The puncture gauge handles collapsing matter without the need for (matter or metric) excision because, during the process, the shift condition pushes the spatial coordinates off of the matter region. This study was followed up with a similar discussion of dust collapse in [37].
A variation of the conformal decomposition, CCZ4, was proposed in [38]. The difference between Z4c and CCZ4 is that CCZ4 includes parametrized constraint addition that, for some choice, correspond to the original four-covariant Z4 formulation. In the generalized harmonic formulation it is known that the inclusion of these non-principal terms can be problematic in some test cases [39], but constraint growth can be mitigated by the use of constraint damping, which is the approach of [38]. For a particular choice of the constraint addition parameters not corresponding to the four-covariant Z4 system, evolutions of binary black hole spacetimes were presented and shown to reduce Hamiltonian constraint violation, by a factor of around four or five (see Fig. 4 of [38]), relative to such simulations for the BSSNOK formulation. On the other hand, in the Z4c system nonprincipal constraint addition that makes the equations of motion as close as possible to those of the BSSNOK formulation whilst still obtaining the desired PDE properties in the constraint subsystem is chosen. It would be interesting to know in more generality how the addition of non-principal constraints affects their evolution.
In this paper we present the first long-term 3D evolutions of black hole and neutron star binaries with the Z4c formulation. In section II the Z4c equations of motion and boundary conditions are presented. We then discuss, in section III, changes to the BAM numerical code since [40]. In particular we describe the implementation of spherical patches for the wave zone [41] and the radiation controlling constraint preserving conditions of [5].
In sections IV and V we present our simulations of single and binary compact objects, respectively. Appendices A and B contain, respectively, descriptions of the spherical patch implementation and evolution of Teukolsky waves, the latter of which we use for code validation. We conclude in section VI.
We use units G = c = 1 throughout, unless otherwise stated.

II. THE Z4C EQUATIONS OF MOTION
In this section we summarize for completeness the Z4c equations of motion, constraints, and boundary conditions. The evolved quantities of the formulation are the conformal spatial metricγ ij , the lapse α, the shift vector β i , the conformal tracefree part of the extrinsic curva-tureÃ ij , the constraint Θ and, up to constraint addition, the trace of the extrinsic curvatureK = K − 2Θ. Finally we evolve the the conformal contracted Christoffel sym-bolsΓ i , which are initially set according toΓ i = −∂ jγ ij . Evolution equations. The equations of motion for the Z4c formulation are for the metric components, for the extrinsic curvature components and for the remaining variables. Here the intrinsic curvature associated with the ADM metric γ ij = χ −1γ ij is written as and we employ the shorthand The derivative operator D i is that compatible with the ADM metric. Numerical evolutions are performed with a particular flavor of the puncture gauge [17,18,21,22] Constraints. The system is subject to constraints of which the latter two, the algebraic constraints, are explicitly imposed in numerical integration. Boundary conditions. When the spacetime manifold has a smooth boundary with spacelike unit (with respect to the ADM metric) normal s i we choose for the trace of the extrinsic curvature the boundary condition where= denotes equality only in the boundary, we use the shorthand ∂ s ≡ s i ∂ i and we use upper case Latin indices A, B, C to denote those that have been projected with the operator q i j = δ i j − s i s j . We can alternatively write this as a second order derivative boundary condition on the lapse. Note that in the expressions for the boundary conditions we never commute the unit normal, projection operator q i j , or physical projection operator q (P ) kl ij = q k i q l j − 1 2 q ij q kl through any derivative operator. For example, ∂ tΓ s = s i ∂ tΓ i . So we have for vectors v i and symmetric tensors S ij , respectively. We refer to the various components of the time derivative under this 2 + 1 decomposition as the scalar, vector and tensor sectors in the obvious way. If we are given conformally flat initial data with β i s i = 0 and α 2 µ S is constant, as is always the case for the data evolved in this work, we choose the boundary conditions for the longitudinal part of the shift. We take for constraint preservation in the scalar sector, where for exampleD i =γ ij ∂ j when acting on scalars, so that the tilde denotes that the conformal metric was used in the contraction. In the vector sector we have for the gauge, and for the constraints, where we denote projected indices by upper case characters starting from the beginning of the alphabet. Here we denote R qq = q ij R ij . Finally we take for the tensor sector. The boundary conditions in the scalar and vector sectors are designed to absorb outgoing constraint violations. The last of the conditions (25) are equivalent to the requirement that Ψ 0= 0 (see [31,32,42] for more details). This condition could also be used with the BSSNOK formulation, but for constraint preservation more work may be needed to adapt the other conditions. In the development of this work we have tried alternative conditions onK andΓ i with only small differences in the outcome. We do not claim that these gauge boundary conditions are optimal. In our numerical experiments we sometimes also employ the more naively constructed Sommerfeld conditions Note that since the trace constraint onÃ ij is constantly imposed, the last of these conditions (30) constitute only five boundary conditions. Since the Z4c formulation coupled to the puncture gauge has ten incoming characteristics in the weak field region, these conditions are not overdetermined, in contrast to the standard conditions used with BSSNOK with box-in-box mesh refinement. Regardless of the formulation, or whether Sommerfeld conditions are taken for every evolved field or just the subsetK,Γ i , Θ,Ã ij , they are not constraint preserving and do not control incoming gravitational radiation.

III. NUMERICAL METHOD AND PARAMETERS
In this section we describe the numerical technique employed in this work. We use the BAM code [40,[43][44][45], a Cartesian-based adaptive mesh refinement (AMR) infrastructure optimized for the evolution of BBH and BNS spacetimes in 3+1 GR. Vacuum spacetime evolutions have also been performed with the AMSS-NCKU code [46], which employs the same methods but with an independent implementation.
AMSS-NCKU and BAM basics. Before discussing the upgrades to the codes used in this work, we summarize the main points of the numerical methods used by AMSS-NCKU and BAM. The evolution algorithm is based on the method-of-lines with explicit Runge-Kutta (RK) time integrators (in this work we employed fourth order RK for vacuum spacetimes and third order RK for non-vacuum spacetimes) and finite differences approximation of the spatial derivatives. The numerical domain is made of a hierarchy of cell-centered nested Cartesian grids (nested boxes centered on the punctures [45,47]). The hierarchy consists of L levels of refinement labeled by l = 0, ..., L−1. A refinement level consists of one or more Cartesian grids with constant grid spacing h l on level l. A refinement factor of two is used such that h l = h 0 /2 l . The grids are properly nested in that the coordinate extent of any grid at level l, l > 0, is completely covered by the grids at level l − 1. Some of the mesh refinement levels can be dynamically moved and adapted during the time evolution according to the technique of "moving boxes". The Berger-Oliger algorithm is employed for the time stepping [48], though only on the inner levels [44]. Interpolation in Berger-Oliger time stepping is performed at second order. A Courant-Friedrich-Lewy factor of 0.25 is employed in all the runs. We refer the reader to [40,46] for more details.
Numerical treatment of the field equations. The BSS-NOK and Z4c equations of motion are implemented numerically in the same way. Fields derivatives are approximated by centered finite difference expressions (fourth order in this work), except for the shift advection terms which are instead computed with lop-sided expressions [40,[49][50][51]. Algebraic constraints are enforced after every time step in BAM and after every Runge-Kutta substep in AMSS-NCKU. The gauge parameters in our numerical simulations are fixed to µ L = 2/α, µ S = 1/α 2 and η = 2/M ADM , unless otherwise stated. Note that this is not a choice of parameters for which the calculations of [5] are expected to guarantee well-posedness of the initial boundary value problem because we have not carefully taken µ S in such a way as to avoid either some generically distinct speeds in the system clashing, or sets of measure zero on which the evolution equations may be only weakly hyperbolic. In earlier studies [38] the choice c > 0 with µ S ≃ c/α 2 was not found to greatly affect the behavior of the scheme in applications. As highlighted in Appendix A of [4] it is challenging to identify problems caused by weak hyperbolicity in applications, even if the system is weakly hyperbolic everywhere in space. If the degeneracy happens on sets of measure zero we therefore expect that in practical applications it will be nearly impossible to identify as the cause of any concrete numerical problem, although in principle such degeneracy should of course be avoided. In the Z4c simulations presented in this work the constraint damping scheme with the values κ 1 = 0.02 and κ 2 = 0 is used. These values have been suggested in the detailed 1D numerical analysis of [3]. Preliminary exploratory runs in 3D indicated the combined use of artificial dissipation and constraint damping terms is important (in some cases essential) to avoid instabilities at the interfaces between boxes and spherical patches, thus confirming some of the expectations from the 1D runs (see also the discussion in Appendix A). OpenMP support in terms of a hybrid OpenMP/MPI implementation has been added inside BAM to key functions with high computational cost, such as evolution equations, interpolation and wave extraction, which improves the efficiency of memory management of the code.
Hydrodynamics treatment (BAM only). The algorithm implemented for the general relativistic hydro-dynamics (GRHD, referred to hereafter as "matter" for brevity) is a robust high-resolution-shock-capturing (HRSC) method [43] based on primitive reconstruction and the Local-Lax-Friedrichs (LLF) central scheme for the numerical fluxes, see e.g. [52]. Metric variables are interpolated in space by means of sixth order Lagrangian polynomials, while matter ones by a fourth order weighted-essentially-non-oscillatory scheme. Primitive reconstruction is performed here with the 5th order WENO scheme of [53], which has been found important for long term accuracy [54,55].
Spherical patches for the wave zone. Both the AMSS-NCKU and BAM codes have been upgraded for this work.
The Cartesian box-in-box mesh refinement has been extended with spherical patches ("cubed spheres") [41,56,57] for the wave zone. They provide us with adapted coordinates for waves, and, as demonstrated in appendix A, improved accuracy in GW extraction. The presence of a spherical outer boundary furthermore allows a straightforward implementation of BCs (either Sommerfeld or CP) due to the absence of corners. When the spherical patches are being used the Cartesian moving boxes, as previously implemented, are employed only in the strong-field region for the simulation of the binary orbital motion, or a collapsing star, while the propagation of gravitational and possibly also electromagnetic waves distant from the source is simulated on cubed spheres [41,57] (see also [58][59][60]). The technique is based on the covering of the sphere by six patches, each patch having local coordinates that are then mapped to Cartesian ones in such a way to avoid pathologies associated with standard spherical polar coordinates. As opposed to the box-in-box setup, spherical patches allow constant radial resolution with linear scaling in the number of grid points, while the boxes result in effectively constant angular resolution as well. In practice, the l = 0 Cartesian box is substituted with six patches overlapping with the Cartesian box at level l = 1 and among themselves. The resolution of h 1 is also the radial resolution employed in the patches. A grid configuration is specified by the number of: i) levels L, ii) grid points in each non-moving box per direction n, iii) grid points in each moving box per direction n mv , iv) the coarsest box per direction h 1 = h, v) grid points in each patch, n r and n θ, φ , which are typically chosen as n θ, φ = n/2. More details on the implementation are given in Appendix A). Finally, GWs are extracted using the Newman-Penrose formalism, in particular by computing the ψ 4 scalar on coordinate spheres in the wave zone (see Sec. III of [40]). Mode decomposition of ψ 4 is performed by projections onto spin weighted spherical harmonics and integration on the spheres with a Simpson algorithm.
Boundary conditions. The radiation controlling, constraint preserving boundary conditions (17)(18)(19)(20)(21)(22)(23)(24)(25) are implemented according to the following simple recipe. Inside every Runge-Kutta substep Lagrange extrapolation, of sixth order in our experiments, is used to populate enough ghostzones in a neighborhood of the boundary, so that the same finite difference and dissipation operators used in the bulk may be evaluated at the boundary. The metric components α, β i , χ,γ ij are updated at the boundary with their standard equations of motion from the bulk, whereas the boundary conditions (17)(18)(19)(20)(21)(22)(23)(24)(25) are used in place of the standard equations of motion forK,Γ i , Θ,Ã ij . Since the evolution system is not symmetric hyperbolic, at least inside a large class of symmetrizers [4], with the standard puncture gauge, we can not rely on a discrete energy method to guarantee numerical stability, even in the linear constant coefficient approximation. We will see however that this implementation of the boundary conditions is numerically well-behaved in the experiments we perform. Sommerfeld boundary conditions are implemented in a similar way; instead of replacing the equations of motion forK,Γ i , Θ,Ã ij by (17)(18)(19)(20)(21)(22)(23)(24)(25) we choose (26)(27)(28)(29)(30) and likewise for BSSNOK, but without the Θ boundary condition.

IV. LONG-TERM EVOLUTION OF SINGLE COMPACT OBJECTS
In this section we present 3D evolutions of single isolated compact objects and compare systematically BSS-NOK and Z4c runs. The main focus is on the behavior of the Hamiltonian constraint violation, due to truncation and artificial BCs errors. We assess long term stability of the Z4c formulation in 3D evolutions of puncture black holes (both non-spinning and rapidly spinning) and compact stars, and demonstrate an overall improvement in constraint preservation and in some instances accuracy of physical quantities (absolute numerical errors at finite resolution) when the Z4c formulation is employed. The spurious effect of Sommerfeld BCs and the improvement obtained with the new BCs are discussed in detail. Our results for spherically symmetric spacetimes are interpreted in view of previous results obtained in [1,2] by means of 1D simulations. The use of Z4c does not significantly improve the computation of the gravitational radiation emitted by a rapidly spinning Bowen-York puncture [61]. In particular we find that in both cases we do not obtain a clear pointwise convergence of the waves at the resolutions used in these tests. This result is possibly related to a lack of resolution, or to the use of punctures in the black hole description, rather than to deficiencies in either formulation.

A. Non-spinning puncture
Evolutions of a non-spinning puncture of mass M = 1 with BSSNOK and Z4c and Sommerfeld and constraint preserving BCs are compared. We confirm qualitatively the results obtained with 1D simulations [1,2]. In particular the new BCs reduce significantly the spurious constraint violation incoming from the boundary in the case of Z4c. The largest single constraint violation however occurs at the puncture, where both formulations give similar results. Setup. We employ five levels of box-in-box mesh refinement, and attach the shells at r ∼ 21.5 M . Each box has n = 40, and the resolution of level l = 5 is h 5 = 0.0625 M per direction. We choose n θ, φ = 40 angular and n r = 40 radial points in each spherical patch so that the outer boundary is located at r = 58.5 M . No symmetries are imposed. The puncture is evolved with a precollapsed lapse as discussed in [40]. In real applications the outer boundary is typically placed further out, perhaps at 500 M or 1000 M . Although this does not solve issues caused by the boundary in principle, and is not an efficient treatment, in practice it reduces some of the features we will encounter here. By design however, in these tests, we aim to see the behavior of the constraint preserving conditions.
Constraint violation in the strong-field region. In each test we find that the Hamiltonian constraint violation in a neighborhood of the puncture is large, taking a maximum value, at this resolution of ∼ 10 3 at around t = 5 M , before being rapidly suppressed to ∼ 3. These values should not be taken particularly seriously because of the finite regularity of the solution at the puncture. Since much of the physics is concentrated around the puncture large numerical error and therefore large constraint violation are to be expected in this region. These violations converge away more slowly than others in the simulations. The Z4c evolutions do not reduce this violation. As observed in spherical symmetry (see Fig. 10 of [3]) we see small amounts of the puncture constraint violation propagating out of the horizon with either formulation. Unsurprisingly, the rectangular mesh-refinement boundaries in the 3D evolutions obscure the feature. The Z4c evolutions do not qualitatively affect this behavior, which seems to improve with resolution.
Constraint violation at large radii. In Fig. 1 the Hamiltonian constraint violation on the spherical shells in space at times t = 75 M and t = 1000 M is plotted. We find that the incoming constraint violation of the Z4c Sommerfeld evolution is comparable to that of the BSS-NOK evolution, although at late times most of the constraint violation caused by the outer boundary has been absorbed by the boundary; the problem is that the constraint violation induced by the Sommerfeld boundary condition, with either BSSNOK or Z4c, does not converge away with resolution. Note that comparing the BSSNOK Sommerfeld evolution with the Z4c CPBC data is not entirely fair, because there is every possibility that constraint preserving boundary conditions for BSSNOK, see for example [23], could also be implemented. In any case, it is evident that the Z4c constraint preserving boundary conditions are helpful in reducing this violation. It is possible to reduce further the constraint violation by using a large constraint damping term, as for example in [38], but in spherical symmetry we found that values of κ 1 > 0.1 M can adversely affect the dynamics of the evolutions with Z4c at late times. We are therefore wary of using large damping parameters, but make no claim about the effect of such damping terms on other conformal decompositions of Z4.

B. Non-rotating star
Our comparison of the new three dimensional numerics with earlier findings in spherical symmetry is completed by the evolution of a single stable neutron star. We find Z4c advantageous in reducing constraint violation that accumulates in grid regions occupied by the matter. In particular the L 2 norm of the Hamiltonian constraint is found to be three orders of magnitude smaller for standard resolutions. The spurious ringing effect due to Sommerfeld BCs pointed out in spherical symmetry is visible also in 3D simulations, but does not dominate the error budget at typical resolutions. In our earlier study [1] we also evolved an unstable single star in a so-called mi-gration test. We suppress such an experiment here; our tests of strong field dynamics are instead performed with binary spacetimes in section V.
Setup. We use three levels of box mesh refinement, and attach the spherical grids at r ∼ 18 M . The lowest resolution runs use n = 48, the star is completely covered by the innermost grid level l = 2 with a resolution of h 2 = 0.36 M per direction. We choose n θ, φ = 48 angular and n r = 48 radial points in each spherical patch so that the outer boundary is located at r = 50 M . Additional runs at twice the resolution (n = 96) are performed. No symmetries are imposed. The star is described by a Γ = 2 polytropic EOS with gravitational mass M = 1.4 M ⊙ and radius 5.7 M, which is exactly the same model evolved in [1] in bespoke spherical symmetry. The evolutions are characterized by oscillations triggered by truncation errors.
Constraint violation and truncation error dynamics. The key findings of the earlier numerical experiments in spherical symmetry were: (i) in the BSSNOK region containing matter, the Hamiltonian constraint inside the star grows throughout the evolution, whereas in the Z4c evolution it does not; (ii) the BSSNOK simulations have larger oscillations (larger truncation errors) than Z4c; (iii) if Sommerfeld BCs and sufficient resolution are employed, incoming modes from the boundary perturb the star amplifying unphysically the oscillations (see Fig. 1 [2]). In 3D simulations we find some similar features. The growth of the Hamiltonian constraint is the dominant one, because of lack of resolution the oscillations are not significantly affected, and the effect of Sommerfeld BCs is less pronounced but visible. The norm of the Hamiltonian constraint is plotted as a function of time in the left panel of Fig. 2. By the end of the simulation the constraint violation in the Z4c evolution is about two orders of magnitude smaller than the BSS-NOK violation and interestingly smaller even than the initial violation. The behavior is similar to that shown in Fig. 5 of [1]. We will demonstrate the growth of the Hamiltonian constraint also pointwise on the grid in section V, in the case of binary spacetimes. Considering momentum constraint violations one finds, as in spherical symmetry, that the differences are far less dramatic. The main difference in this case is that the momentum constraint violations with BSSNOK are more dynamical, and slightly larger, although this is probably just because of reflections from the Sommerfeld outer boundary condition, which propagate many times over the computational domain. With either BSSNOK or Z4c, the largest persistent violation in the momentum constraint occurs at the surface of the star. In the right panel of Fig. 2 we show the oscillations of the central density during the simulation time. It is not possible to distinguish significant differences, probably because of the low resolution employed. The large effect seen in [1] was with a resolution ten times higher than here.
The Sommerfeld boundary kick. In the left panel of Fig. 2 one sees the effect of the Sommerfeld BCs with Z4c on the norm of the Hamiltonian constraint. As in the single puncture evolution, roughly when the outer boundary becomes causally connected to the central body, there is a large incoming pulse of constraint violation. This pulse perturbs the central object, but unlike in the spherical case (see Fig. 1 of [2]), the violation is not the dominating effect on the dynamics. In the right panel of Fig. 2 we do not see a significant jump in the central density. A blow-up of the central density plot is however shown in Fig. 3, together with data from a shorter run at twice the resolution. The figure demonstrates that at approximately 60M the star is slightly perturbed by the Sommerfeld BCs. More importantly, the figure shows that the size of the perturbation is not converging away with resolution, whereas the amplitude of the oscillations does, so we expect this error to be dominant at higher resolutions. By contrast CPBCs are characterized by smaller reflections. The Hamiltonian constraint violations propagating out from the star appear to converge at approximately second order, in line with our expectation for our hydrodynamics scheme. This rate of convergence is maintained by the constraint preserving boundary conditions (see also Fig. 5 of [2]).

C. Rapidly spinning puncture
Here we compare the evolutions of a single Bowen-York spinning puncture [61][62][63] with spin a = S/M 2 ≃ 0.92. We choose such a comparatively large value of the spin parameter (see also [64]) to test performance for large values of the conformal factor. For a ≃ 0.92 the puncture contributes only 24% of the mass, while the remainder is in the Brill wave contribution of the conformal factor. Both formulations give comparable results in terms of stability, the norms of the constraint violation, the majority of which occurs at a few points near the puncture, and gravitational waves. With either system we observe that pointwise fourth order convergence is not achieved in the GWs at the resolutions at which we performed the tests, although global errors scale between third and fourth order. We expect that this behavior is related to inaccuracies caused either by lack of resolution, or intrinsic to the puncture description of the black hole.
Setup. We use five levels of box mesh refinement, and attach the spherical grids at r ∼ 50 M . In the lowest resolution runs each box has n = 48 points per direction, and the resolution of level l = 5 is h = 0.065 M . We choose n θ, φ = 48 angular and n r = 48 radial points in each spherical patch so that the outer boundary is located at r = 150 M . Runs at resolutions n = 64, 96, 128 (with the grid spacing scaled in order to maintain the same grid setup) are performed. Bitant symmetry is imposed. The initial data for this test is constructed in with an ad-hoc method, used elsewhere in the literature, in which the BAM spectral binary puncture initial data solver is applied to a single puncture with large spin and another, unboosted and unspinning, located very close   [1,2]. This effect does not converge away with resolution, but at these resolutions is not the dominant effect.
by with a relative uncorrected mass of 10 −12 . There is no sign of the second puncture on the grid. Problems with convergence do however make this construction a point to address. Basic features of the dynamics. The (2, 0) multipolar mode of r ψ 4 emitted during the evolution of this initial data are shown in Fig. 4, for the lowest resolution runs. The waves lie on top of each other; as resolution is increased they converge to each other. At late times, t = 350 M , a boundary effect is visible in the BSS-NOK data. The BCs of Z4c improves this behavior significantly. Note that t = 350 M corresponds roughly to the time needed by a wave initially near the puncture to propagate to the outer boundary, be reflected from the strong field region and travel out once again to the extraction sphere at 90 M . No feature is visible as the wave passes from the outer boundary at around t = 210 M , perhaps because ψ 4 is by construction insensitive to incoming gravitational radiation. We computed the ADM mass integral (see Eq. (32) later) and find that in the BSSNOK runs a drift begins exactly when the outer boundary becomes causally connected with each observer, although this drift is then swamped by the effect of the outgoing gravitational waves. The effect of the Sommerfeld boundary condition on physical quantities is discussed in more detail for binary neutron star simulations later. times, up to t ∼ 100, only second order pointwise convergence is observed with either BSSNOK or Z4c. Later, the errors scale at third-to-fourth order rate in magnitude (norm) but pointwise convergence is lost. The use of either more resolution, properly constructed initial data, or simply more refinement levels, i.e. high resolution at the puncture, might improve this behavior. In [40] similar convergence tests were performed but with a lower spin a = 0.194, with two more, a total of seven, levels of mesh refinement, but with a smaller range in resolutions, and a lower maximum resolution in a neighborhood of the puncture. Although the plots presented in [40] (see Fig. 6) are scaled for fourth order convergence, the difference in resolutions makes it hard to distinguish between different convergence factors with confidence. The earlier study also found pointwise convergence is not maintained in the gravitational waves, thus our findings are consistent. A detailed discussion of spin and higher order finite differencing can be found in [64].

V. EVOLUTION OF COMPACT BINARY INITIAL DATA
In this section we assess the performance of the Z4c formulation in the simulation of the merger of two compact objects. Two standard initial configurations are studied, with a set of approximately 60 simulations, in detail: an equal-mass, non spinning, black hole binary (BBH), and an equal-mass, irrotational neutron star binary (BNS).
We discuss evolutions of about three orbits and compare systematically the results from several convergence series with the corresponding ones obtained with BSSNOK. Presentation rubric. In the presentation of the results, we describe the main features of the dynamics, then the Hamiltonian constraint violations are compared and finally the accuracy of the GWs and other physical quantities are discussed in relation with the Hamiltonian constraint violations. For brevity we discuss only the main emission channel, the (2, 2) multipole of the radiation. As is standard, the (complex) ψ 4 projection r ψ 4 22 , extracted on a coordinate sphere S r of radius r, is decomposed into amplitude and phase according to On the same coordinate sphere we compute the integral the spatial metric γ ij = χ −1γ ij which approximately represents the ADM energy. The ADM energy (or mass) of the system is given by and it is a conserved quantity. On a given sphere in the wave zone, however, E ADM (r) is expected to deviate from M ADM due to the gravitational energy radiated away from the sphere, Note that the outer integral in Eq. (34) is performed by a simple Riemann sum. The angular momentum is computed with similar ADM-like integrals, see e.g. [65] and references therein, but note that this quantity remains ambiguously defined (gauge dependent) also in the asymptotic limit r → ∞ for a generic asymptotically flat spacetime. We stress that the observed differences, up to outer boundary effects, seem to converge to the same continuum solution. They are however significant at fairly high resolutions in most of the cases we studied. Particularly relevant are the differences in Hamiltonian constraint violation in the evolution of non-vacuum spacetimes. We show evidence that these violations are correlated with the quality of the numerical waveforms. Convergence tests. Standard three level selfconvergence tests are performed using simulations at different grid resolutions. These tests can be biased by the choice of the triplet. As discussed in [55] we consider only triplets in which i) the ratios between the low and medium and medium and high resolutions are h L /h M ≃ h M /h H > 1 (ideally > 2); and ii) the scaling factor is at least of order two, where p is the convergence rate. We find that if these criteria are not met then the measured convergence order from different triplets is not consistent. Additionally, in our experience we found it important to verify that the use of different triplets gives consistent results. Even when the criteria above are satisfied, different triplets can give differing convergence factors. The result is acceptable if the rates consistently improve as higher resolution triplets are taken.
Initial data. Before discussing the results we summarize the initial data and grids employed. Initial data for BBHs are conformally flat puncture initial data constructed using the Brandt-Brügmann puncture method [66] with the BAM implementation of the spectral puncture initial data solver [67]. The holes have equal mass, an initial separation of d = 7 M and are placed in a quasi-circular configuration, on which the Padé resummed eccentricity reduction algorithm of [68] was applied. The initial data are interpolated onto the AMSS-NCKU and BAM grids by eighth order Lagrangian barycentric interpolation. Initial data for BNS assume a conformally flat metric and irrotational flow.  I: Summary of the grid configurations and of the runs. Columns: name of the configuration, maximum refinement level, moving levels are those with l > L mv , number of points per direction in the moving levels, resolution in the level l = L−1, number of points per direction in the non-moving levels, resolution in the levels l = 1 (radial resolution in the shells), number of radial points in the shells, number of angular points in the shells, outer boundary. Note that the resolution is given in units of M for BBH runs, but in units of M⊙ for BNS runs. Runs marked with "*" were reproduced with the AMSS-NCKU code.
Name L L mv n mv hL−1 n h1 nr n θ, φ rout BBH0* 9 2 48 0. 0182  brary and have been already evolved in several places [43,55,70]. The data are interpolated onto the BAM grid by spectral interpolation. Grid setup. All of the runs performed at different resolutions are listed in Tab. I, together with details of the grid setup. Note that several of the BBH simulations, those marked with a " * " in Tab. I, have been performed with both the BAM and AMSS-NCKU codes, indicating, at least for the chosen grid settings, the robustness of our findings. Note that our lowest resolution vacuum setup BBH0 has maximum resolution comparable to that of the highest resolution in the earlier BAM calibration paper [40], although the boxes in [40] were larger which also affects accuracy. The resolutions of the setups BBH2 and BBH3 are comparable to the ones of recent BBH simulations [71,72] that required only moderate accu-racy, while the maximum BBH9 resolution is similar to what was used to obtain some of the BAM waveforms for the NINJA-2 catalog [73]. The highest resolution in the matter simulations BNS6, BNS6r approaches those in the accurate runs of [54,55]. Other details about gauge conditions, damping parameters, etc. have been already given in Sec. III.
A. Equal-mass, non-spinning BBH Basic features of the dynamics. Let us discuss the evolution of BBH initial data. The black holes evolve for about 2.5 orbits before merging, radiating energy and angular momentum in gravitational waves. In Fig. 6 the puncture tracks are plotted for the grid BBH2 in Tab. I. The gravitational waves have exactly the same basic profile, but unlike in the evolution of a single spinning puncture (see Fig. 4) they are visually distinguishable, with the absolute maximum of the (2, 2) multipole mode of the real part of r ψ 4 occurring, approximately 2.2 M , earlier in the BSSNOK data, in accordance with the expectation from Fig. 6. The differences that accumulate over the evolution converge away with resolution. In the BBH4 data the delay is only 1.1 M . Also in contrast to the single spinning puncture case, presumably because the outer boundary is placed so far out at r ≃ 2090 M , no boundary feature is visible in the BSSNOK waves, at least up to a radius of r = 400 M within the run-time of the simulation.
Constraint violation. In Fig. 7 the Hamiltonian constraint log 10 |H| is plotted in space for the BBH1 runs on the orbital plane at a simulation time t = 146 M (roughly 1.5 orbits) on refinement levels l = 3, 4, 5. As in case of the single puncture evolution the Hamiltonian constraint violation is dominated by the punctures and, for this data, at this time, has a maximum absolute value of ∼ 10 −2 regardless of the formulation. However, the Hamiltonian constraint violation differs in the strong field region depending on the formulation. In the BSS-NOK case a significant Hamiltonian violation extends to level l = 3 with an almost spherical pattern. In the Z4c case the violation is mainly restricted to the highest level around the puncture. Note the effect of Cartesian mesh refinement in the plot. For the same data the momentum constraint away from the puncture is also roughly an order of magnitude smaller in the Z4c data. The highest resolution runs, BBH9, have the smallest constraint violation even in the strong field region, but the violation there is dominated by that of the puncture. The difference between the Hamiltonian constraint violation of the BBH1 and the BBH9 data in this region is at most a factor of three for either formulation despite the difference in resolution.
Gravitational wave accuracy. Quantitative differences are observed in the gravitational radiation computed in BSSNOK and Z4c simulations. The differences can be appreciated in self-convergence tests for phase  Tab. I. In the continuum limit, upto outer boundary effects, the tracks should agree perfectly because we are evolving the same data with the same gauge choice regardless of the formulation. At finite resolution however different formulations and discretizations may give different results. Initially the tracks agree, but a difference accumulates; the BSSNOK punctures merge slightly sooner than the Z4c ones. This plot does not indicate either set of numerical data is better than the other, only that there is qualitative agreement between the results. and amplitude. Typical results are presented in Fig. 8. The phase errors accumulated to merger (t ∼ 270 M ) are shown in bottom panels. When the Z4c simulation is employed we observe a factor three of improvement. This behavior is common to all of the simulations performed. Similarly, the amplitude error improves by a factor twoto-three, as can be appreciated from the top panel of the figure. The results obtained with the Z4c are not only characterized by smaller errors, but also by the fact that convergence is maintained for longer time, in particular in phase errors. All the errors are observed to converge with increasing resolution. The Z4c and BSSNOK data sets also appear to converge to each other, suggesting the same continuum limit is approached. We conclude that at finite resolution the use of Z4c results in the computation of more accurate waveforms.
Convergence issues. For the converge tests presented in Fig. 8 third order convergence is obtained, despite the use of fourth order operators for the bulk derivatives (see below). In our experiments with the described grid settings we find it difficult to achieve fourth order convergence in GWs from orbiting puncture runs, regardless of the formulation employed. On the other hand clearer (order approximately 3.5) convergence is demonstrated in Appendix B for the Teukolsky wave, a simpler problem which however retains some of the important features, namely, nonlinearity, in the sense that we evolve with the full BAM Z4c infrastructure, constraint violation, because the wave only satisfies the constraints to linear order in perturbation theory, present in the BBH calculations. On the other hand the Teukolsky wave test has neither moving boxes or punctures. We were unable to identify the precise reason for the behavior in presence of punctures, however we point out the following well-known facts: i) punctures have finite regularity, which can obviously affect the formal order of convergence of high-order finite-differencing stencils; ii) moving box simulations have several source of errors (second order time interpolation, regridding, the number of mesh refinement buffer zones, 3D interpolation, data between boxes and shells, spheres and mode projection), a clearly defined error budget is very complicated to construct (for a discussion see e.g. [40,41,74]); iii) simulations have several freely specifiable parameters (grid parameters and/or gauge parameters), and systematically tuning all of them is beyond the aim of this work; iv) the convergence order observed actually depends on the triplets chosen in the self-converge test. The experimentally measured convergence factors lie between two and four for both the formulations during the inspiral. Around the time of the merger the BSSNOK convergence rate drifts up to five or six. The drift is on average smaller for higher resolution triplets. This does not mean that the BSSNOK data are converging at a high order, but rather that the errors in the simulation do not allow us to judge the rate in a meaningful way with a simple error model.
We do not claim that waveforms from either formulation are unreliable, but we are reluctant to estimate the absolute errors by Richardson extrapolation, which assumes a certain rate of convergence.
AMSS-NCKU-BAM comparison. The qualitative features of these results have also been observed in simulations with the AMSS-NCKU code.
Starting with the puncture tracks as in Fig. 6, we find that the BSSNOK also merge slightly earlier at a given resolution, roughly 3.2 M for BBH0 and 1.0 M for BBH2, in the same ball-park as the values obtained with BAM. Note that we do not expect to obtain identical values from the two codes, because although they share many ingredients, some specifics, for example the shells implementation and grid placement, differ. As in Fig. 7 the Hamiltonian constraint appears to be smaller in the Z4c data when using AMSS-NCKU. To test the robustness of this finding we have made a number of experiments with different precollapsed initial profiles for the lapse. We also tried different constraint damping factors for κ 1 , up to κ 1 = 0.1. Neither change seems to make a significant difference; the Z4c Hamiltonian constraint is always smaller. Note that in the earlier stability work [4], the Hamiltonian constraint was computed directly from the conformal variables. Now, as in BAM, we compute it by transforming first to the ADM variables. The two versions of the calculation differ by additions of Θ andZ i , and obviously the finite difference approximation.
The AMSS-NCKU code appears to produce slightly larger errors at mesh refinement boundaries than those of BAM. The reason for this is currently unclear, but may be the cause of the larger errors that we find in the AMSS-NCKU wave-forms. Nevertheless we still find, as in Fig. 8 that the Z4c wave forms are roughly twice as accurate in phase and amplitude.

B. Equal-mass, irrotational BNS
Basic features of the dynamics. Let us discuss the evolution of BNS. The binary evolves approximately two orbits before contact, then merges forming a hypermassive neutron star (HMNS) which finally collapses on dynamical timescales. Dynamics and related gravitational wave emission were described in detail in our previous work [43] so we do not repeat them here. We only mention that the gravitational emission is characterized by approximately six cycles during which the GW fre-quency increases monotonically (even after contact), the merger time is defined conventionally at peak of the amplitude's (2, 2) mode, nonlinear (quasi-radial and nonaxisymmetric) oscillations of the HMNS generate the post-merger signal which decays exponentially after collapse. Figure 9 shows snapshots of the rest-mass density on the orbital plane around merger for a typical simulation obtained with the two formulation (and same setup). There are visible differences in the "dynamics", and in the Z4c run the contact and merger happens at a later simulation time. As in the case of the BBH puncture tracks, these plots are gauge dependent but they refer to the same gauge choice up to truncation errors and so they can be compared. In this respect, note that the centroids of the stars are offset, this is a coordinate effect due to the choice η = 2 in the shift condition and was studied  Fig. 6 apply; namely, it seems that the compact objects merge earlier in the BSSNOK data. In this case however there is already an expectation that the Z4c data will be more accurate because in earlier work [43] it was found that the objects merge later and later as resolution is increased with the BSSNOK formulation.
Constraint violation. The L 2 norm of the Hamiltonian constraint violation during the evolution is reported in Fig. 10. On the refinement level l = 2 one can observe an improvement of a factor ∼ 100 during the whole evolution if Z4c is employed. By contrast the norms of each component of the momentum constraint agree for both evolutions, an almost constant violation around 10 −7 is observed. The Hamiltonian violation on the orbital plane around merger time is shown in Fig. 11, which includes now refinement levels l = 2 and l = 3. Even using the highest resolution run (BNS6) we register a two-to-three order magnitude difference in the absolute value of the Hamiltonian constraint. It is noteworthy that the violation has an almost spherical pattern on and around the strong field region of the binary, which suggests that the violation is not dominated by rectangular mesh refinement boundaries.
Gravitational wave accuracy. The differences at finite resolution between BSSNOK and Z4c described so far have an impact on the computation of physical quantities: gravitational waves and ADM mass. Figure 14 shows the result of a standard three-level self-convergence test for phase and amplitude of the (2, 2) mode of the gravitational radiation emitted. For the particular triplet shown BSSNOK data cease converging before contact. An analogue result was already found at similar resolu-tions in [43] (see also the detailed discussion in [54,55]). On the other hand the Z4c data are found to converge at approximately second order rate beyond contact and up to merger (t ∼ 650 M ). Similarly to the BBH case, Z4c waveforms are found to be more accurate by a factor three-to-four in accumulated phase and amplitude to merger time. We relate this behavior with the improvement obtained in Hamiltonian constraint preservation. Because truncation errors and the Hamiltonian violation (especially in the matter region) in BNS simulations are generically larger, the use of Z4c makes a significant difference for the accuracy of these simulations.
Mass and angular momentum conservation. As in the case of a single spinning puncture we consider the ADM mass integral, Eq. (32), and the energy radiated in GWs, Eq. (34). In Fig. 12 we show that Z4c permits a reliable computation of E ADM (r). When corrected for the GW energy, the conservation of the ADM mass is of order of 0.1%. By contrast BSSNOK data do not even allow for a reliable estimate of M ADM . Note that the ADM mass can be also estimated by means of a volume integral rather than a surface integral, see e.g. [75,76]. The more expensive volume integral computation can be more accurate and is found to give better results for BSS-NOK in the case of black hole binaries [76]. Figure 13 shows the ADM angular momentum integral, both with and without a correction by radiated angular momen-  Fig. 9. At this resolution the Z4c violation is between one and two orders of magnitude smaller throughout the evolution.
tum, at an extraction radius of 300 M , for both BSS-NOK and Z4c with CPBCs. The outer boundary for this run (BNS2a in Tab. I) is at approximately 360 M . Evidently at either radius the BSSNOK data is very poorly behaved, and has large error. In contrast in the Z4c evolution the corrected angular momentum are well conserved until the merger signal, which we expect to be much less accurate, reaches the extraction sphere. Early in the simulation the BSSNOK data are also conserved, but there is a jump at exactly the time t = 60 M when the outer boundary becomes causally connected to the extraction radius. We checked that this feature holds on every extraction sphere. This demonstrates clearly that the Sommerfeld boundary condition has an effect on the behavior of physical quantities inside BSSNOK simulations. One way to avoid the feature would be to place the outer boundary further away, which, depending on the physics of interest, may not be prohibitively expensive with the spherical shells for the wave zone, but it is not desirable to discard every extraction sphere as soon as it becomes causally connected to the outer boundary, because a larger domain would ideally be used for more reliable wave extraction, rather than as a buffer for poor boundary conditions. For comparison it would be very interesting to see data from BSSNOK simulations, where there is still a zero-speed mode in the constraints, with an implementation of the constraint preserving conditions of [23]. It is natural to compare the early times in Figures 12 and 13 to Fig. 11 of [77], a similar plot, computed with the spectral Einstein code, for inspiralling black hole neutron star binary data. The Z4c results are competitive, although it remains to be seen if they will continue to be so over many orbits.

VI. CONCLUSIONS
This paper is the conclusion of a series [1][2][3][4][5], the aim of which was to bring the advantages of the generalized harmonic formulation to the moving puncture method. We presented here for the first time 3D numerical relativity simulations of compact binaries performed with Z4c, a conformal decomposition of the Z4 formulation.
We started with evolutions of single compact objects and found that the expectations obtained by earlier work in spherical symmetry are largely borne out in the 3D numerics. The most striking feature in these tests is that, in the evolution of a single stable star by t = 1000 M at the resolutions used in our tests, is that the norm of the Hamiltonian constraint is approximately three orders of magnitude smaller in the Z4c data. For the first time we have presented results which combine radiation controlling, constraint preserving outer boundary conditions with the moving puncture method. These boundary conditions removed a perturbation to the central rest mass density of the star which is present when using Sommerfeld boundary conditions. At the resolutions of our tests however the outer boundary is typically not the leading order contribution to numerical error. In evolutions of a single spinning puncture we found that the Z4c with the new outer boundary conditions remove certain constraint violating features present in the BSSNOK data, but at least at early times the qualitative physical picture is unaltered.
We then compared evolutions of compact binary spacetimes. In these tests we placed the outer boundary much farther away from the central body so as to simplify the discussion. Throughout the evolution of binary neutron star initial data, we find that at the same resolutions the Z4c formulation has between one and two orders of magnitude less Hamiltonian constraint violation in the norm. Interestingly the Hamiltonian constraint violation in the Z4c tests in this case stays at or below the level in the initial data, at least until the stars merge, and the simulations may therefore even be competitive with those of a constrained formulation. We find similar, albeit much less pronounced effects in the constraint violation in binary black hole simulations, but the change of formulation does nothing to cure the dominant constraint violation at the punctures themselves. The higher quality of the Z4c data is also apparent in physically meaningful quantities. In terms of gravitational wave accuracy we find that with any triplet, satisfying certain criteria, for either binary neutron star or binary black hole data, the absolute error in either amplitude or phase of the extracted gravitational waves is between two and four times smaller in the Z4c evolutions. The difference, in the evolution of compact binaries, between conservation of the ADM mass integral with the two formulations is remarkable. In the BSSNOK simulations one can not reliably correct the integral with the radiated gravitational wave energy to arrive at a constant. In the Z4c simulations near perfect conservation is achieved. Furthermore,  Fig. 7 is to the BBH data. In this case there is no sign of the Cartesian grid structure in the violation. Most of violation in the BSSNOK simulation appears inside the stars. The Z4c violation seems to be on average one or two orders of magnitude smaller than the BSSNOK violation, as seen more clearly in Fig. 9.  [40]. The left panel shows results for BSSNOK, the right panel for Z4c; in comparison the latter demonstrates remarkable conservation.
despite placing the outer boundary at a large coordinate radius, we find that the BSSNOK data are corrupted, for example in the angular momentum, by the Sommerfeld boundary condition, whereas the Z4c data are free of this problem. If nothing else this motivates the use of constraint preserving boundary conditions with the BSS-NOK formulation.
In summary we have presented a large suite of numerical experiments in which the Z4c formulation was shown to give more accurate results than BSSNOK, both in terms of constraint violation and extracted physical quantities. We therefore expect that the Z4c formulation will become a standard tool for numerical relativity. The radiated angular momentum JGW is computed according to standard formulas, see for example [65]. The left panel shows results for BSSNOK, the right panel for Z4c. The jumps in the BSSNOK data happen exactly as the extraction spheres become causally connected to the outer boundary at 360 M . Other extraction radii show the same feature. ful discussions. We are especially grateful to Denis Pollney for his patient explanation of the Llama multi-patch communication algorithm. We are grateful to Tim Dietrich for providing the initial data for the spinning puncture evolutions. This work was supported in part by DFG grant SFB/Transregio 7 "Gravitational Wave Astronomy" as well as by NSF grants PHY-0855315 and PHY-1204334. Computations were performed primarily on JUROPA (JSC, Jülich), LRZ (Munich), and LOUHI (CSC) cluster. Computer time on LOUHI is granted under PRACE-Tier-1 2011.

Appendix A: Spherical patches
In this appendix the implementation of the spherical patches in the BAM and AMSS-NCKU codes is described. We follow closely [41,57], to whom we refer for further details. We describe first the BAM implementation and then highlight the differences in the AMSS-NCKU code.
The grid structure. The BAM grid is made of a hierarchy of nested Cartesian grid boxes, each level is labeled by the integer l = 0, 1, .... The level l = 0 (outermost box) is replaced with six patches with local coordinates a j = {R, φ, θ}. The grid of each patch is uniformly spaced in local coordinates. The maps between the local coordinates and the Cartesian coordinates are ±x patches: φ = arctan(y/x), θ = arctan(z/x) (A1) ±y patches: φ = arctan(x/y), θ = arctan(z/y) (A2) ±z patches: φ = arctan(x/z), θ = arctan(y/z) (A3) where (φ, θ) ∈ (−π/4 : π/4) × (−π/4 : π/4). The local radial coordinate R is a fisheye coordinate [78][79][80] as function of the radius r = x 2 + y 2 + z 2 . The particular function function f (r) implemented is where the parameters A and B determine the stretching of the coordinate and ǫ the transition region. The condition f (r) = r at interface between spherical patches and the boxes is important to avoid a step behavior and minimize interpolation errors. In this work we did not employ the fisheye coordinate, i.e. we set A = 0, B = 1. Spatial derivatives. Field derivatives at a grid point inside the patches are calculated by finite differences on the uniformly spaced local grid. The derivatives in Cartesian coordinates are obtained using the chain rule and the Jacobian of the transformation x i (a j ), note that the last term in these equations fixes a typo in equation (5b) in [41]. Grid schematic. The grid structure is sketched in Fig. 15. The dark solid lines represent the physical grid, the lighter lines denote the ghost points which are needed for the finite differences. The green shaded regions denote the ghost zones populated by inter-patch interpolation, and overlap with the neighbor patches. The other colored areas overlap with the l = 0 level, and the resolution of the Cartesian box is the same as the one of the radial one in the patch. The number of grid points in all the ghost zones (cyan, green, yellow) is equal. The distance between the points r 0 and r 1 is equal to that between r 1 and r 2 .
Data communication. Spherical patches overlap on ghosts zones. Two neighboring patches share the radial coordinate and the angular coordinate perpendicular to the mutual boundary. Therefore only a 1D interpolation parallel to the boundary has to be performed in the green regions of Fig. 15. A Lagrange interpolation (sixth order in this work) which uses the most centered possible stencil is employed. Interpolation between patches and the l = 0 level is performed with Lagrangian polynomial interpolation in 3D (colored regions of Fig. 15). The various interpolations are done in the following order: i) interpolate from box to shell (yellow region), ii) interpolate between shells (green region), iii) interpolate from shell to box (cyan and red regions), iv) set symmetries in box. For simplicity, grid symmetries are not applied in the shells during evolution. Each patch is evolved entirely and only afterward values at symmetry points are overwritten by copying them. We note that because the interpolation of the red region in Fig. 15 depends on points in the cyan region, the ordering of the different interpolation can give, in principle, different results. On the other hand, not interpolating the red region results in a "double evolution", which we find leads to high frequency oscillations.
Dissipation. Artificial dissipation is necessary for numerical stability during the evolution. In particular to maintain stable box-shell interface regions. Experimentally we found it important to apply different amounts of dissipation in box and spherical patches. In particular we use a lower dissipation on the spherical patches than in the box. For this work we tested only sixth order Kreiss-Oliger dissipation operators, using artificial dissipation coefficients σ = 0.5 in the bulk and σ = 0.1 in the shells (see [40] for our terminology). The number of angular points has been chosen according to the size of the l = 0 box, n θ, φ ∼ n/2. Placing the box-shell interface close to the strong field region may also cause large growth of the error, even causing the code to crash, or affect the choice of the dissipation parameters required for stability.
Parallelization. In order to reduce global MPI communication during shells synchronization, we use an opti-mized parallelization. The computations on the patches are distributed only in the radial direction. Every processor has six patch parts with the same radial extension, hence the synchronization in the angular directions is performed locally by each MPI job. Note that this method implies a minimum number of radial points is necessary for a given number of processors. In some early tests we found that avoiding interpolation in this way can result in computations that are an order of magnitude faster in the shells.
Variations of the numerical method in the AMSS-NCKU code. There are several differences in the implementation of shells in AMSS-NCKU with respect to BAM. In AMSS-NCKU there is no option to use a radial fish-eye coordinate. Relating to the interpolation between box and shell, the scheme is somewhat different to BAM. We set the length between the point r 0 and r 1 as the length of six points and we take these points as buffer points. We set the distance between r 1 and r 2 larger than the distance between r 0 and r 1 . We take this part of region as the double cover region of both box and shell. We set six points for the box outside of r 2 and take them as buffer points. Similar to the treatment of the mesh refinement interface, we fill the buffer points for box and shell only at the end of a full RK4 step. As opposed to BAM we do not interpolate the points between r 1 and r 2 . Due to this double cover region, we can interpolate between box and shell in parallel fashion. But before that we have to synchronize the data between different shell patches. In the MPI communication, we divide the data in both radial direction and angular directions but try to make resulting blocks as cubic as possible for minimization of inter processor data change.
Box-Shell comparison. Let us briefly discuss code performance and waveform quality by comparing runs which employ spherical patches ("shell runs" hereafter) in the wave zone to those that do not ("box runs" hereafter). For shell runs, radial and angular resolution can be adjusted separately, at linear scaling in the number of grid points in radius and quadratic scaling for both angular directions. Nested boxes imply effectively constant angular resolution for increasing radius and decreasing radial resolution. Increasing radial resolution leads to a cubic scaling in the number of grid points, because angular resolution increases simultaneously. This can be partially compensated by choosing larger boxes in the wave zone, which are comparatively cheap in run-time because of Berger-Oliger time adaptivity, and typical runs are limited by run-time rather than memory. Constant radial resolution is helpful for tracking waves traveling to infinity, on the other hand decreasing radial resolution is a recipe sometimes employed deliberately for filtering features going to and coming from the outer boundary. An entirely different type of filter is effectively in use when extracting only (2, 2) modes. Features due to rectangular outer boundaries and rectangular refinement boundaries (cmp. Fig. 2) are not visible at l = 2, but start to show at l = 4. A spherical outer boundary has the advan- tage that a unique normal vector to the boundary exists, in fact our early (and incomplete) experiments with the new Z4c boundary conditions failed for a cubical outer boundary. However, a spherical boundary can also lead to a focused inward reflection of outgoing waves, while a rectangular boundary scatters spherical waves leading to a diffuse reflection. For example, the back-reflection feature in Fig. 4 is much smaller for cubical outer boundaries, as are oscillations of a neutron star induced by boundary reflections. A clean treatment of a spherical boundary is nevertheless preferable since the reflections can be minimized, and even though a constraint violation from a cubical outer boundary may not be as visible, it is still present in the inner domain. The choice between box or shell runs actually depends on the waveform accuracy goal, and requires a balance between accuracy and computational cost. We discuss two examples employing the BSSNOK formulation. shell runs when the four coarsest levels are substituted by spherical patches that approximately extend to the same outer radius. The computational cost of the shell runs is approximately twice the box runs. Using this triplet, third order self-convergence is achieved in the shell runs at all the extraction radii, while in the box runs self-convergence is lower than three at small radii and progressively degrades when the extraction is performed on coarser levels. The waveform in box runs is more noisy. Phase and amplitude of the waveform can be extrapolated according to where f (u, R) is the quantity to extrapolate extracted at radius R, u is a retarded time, and R the Schwarzschild radius, see e.g. [55,81]. A measure of the extrapolation error is the difference between the extrapolated function and the last radius values. This error is reported in the upper row of Fig. 16 for box runs and the lower row for shell runs, using the radii R = {150, 140, 130, 120, 110, 100, 90, 80, 70, 60, 50} and different values of K, i.e. different polynomials. As evident from the figures the extrapolation for shell runs is more accurate, and results for different K are consistent with each other (no oscillations and overshooting). The "zig-zag" behavior in the left panel of the lower row of Fig. 16 is due to the limited number of digits used in the output. For these data the choice K = 1 seems to be optimal.

Appendix B: Evolution of Teukolsky waves
In this appendix we evolve, for code validation, Teukolsky's wave solution with the BAM code.
Motivation. The codes used in this work have sources of error of different polynomial order in the grid spacing, but the operations performed most often, finite differencing and time integration, are performed at fourth order. It is perhaps not surprising that we never find clear fourth order convergence in our evolutions of compact binaries.
For the hydrodynamics simulations this behavior can be excused because the HRSC scheme we employ is only second order accurate. But for our vacuum simulations the situation is less clear. We therefore also studied the evolution of Teukolsky's solution [82] to the linearized Einstein equations. This solution represents a weak gravitational wave propagating on a flat background. It satisfies the constraint equations of General Relativity at linear order. We use this solution as initial data and then evolve forward in time with the BAM code. Note that extensive convergence testing demonstrating fourth order convergence, in test cases, was shown for the AMSS-NCKU code in [4].
Teukolsky wave initial data. The metric has the form ds 2 = −dt 2 + (1 + Af rr )dr 2 + 2Bf rθ rdrdθ + 2Bf rφ r sin θdrdφ + (1 + Cf (1) θθ + Af (2) θθ )r 2 dθ 2 + 2(A − 2C)f θφ r 2 sin θdθdφ + (1 + Cf (1) θθ + Af (2) θθ )r 2 sin 2 θdφ 2 (B1) where we have f rr = sin 2 θ(cos 2 φ − sin 2 φ) , f rθ = sin θ cos θ(cos 2 φ − sin 2 φ) , (B3) f rφ = −2 sin θ sin φ cos φ , and also This corresponds to an outgoing wave that is a pure m = 2 mode. The generating function F is given by Derivatives of order n of F are denoted by F (n) . Apart from choosing a different generating function our wave is very similar to what was used in [83]. For our runs we have chosen a = 10 −4 , N = 10 and λ = 10. Setup. We have evolved such a wave with an initial lapse of one and a shift of zero. For the evolution, as in the rest of the work, we have chosen the standard 1+log and Gamma driver evolution equation for lapse and shift (11)(12) with µ L = 2/α , µ S = 1/α 2 . We use four levels of box mesh refinement, and attach the spherical grids at r ∼ 30. In the lowest resolution runs each box has n = 48 points per direction, and the resolution of level l = 4 is h 4 = 0.167. We choose n θ, φ = n/2 angular and n r = n radial points in each spherical patch so that the outer boundary is located at r = 100. Runs at resolutions n = 48, 64, 80 (with the grid spacing scaled in order to maintain the same grid setup) are performed. Results. We find that both the BSSNOK and the Z4c system can successfully be used to evolve these waves. For Z4c we have set κ 1 = 0.02 and κ 2 = 0. Since there are no strong gravitational fields we can extract the waves at any radius. Figure 17 shows a convergence plot for waves extracted at a radius of r = 80. The solid lines shows the difference between ψ 4 at low and medium resolution (L − M ), while the dash dotted lines show the difference between medium and high resolution (M − H). When scaled with the proper factor of 3.66 for fourth order convergence, we see that M − H coincides approximately with L − M . For this data set the observed order of convergence is around 3.5. The left panel shows the results for BSSNOK, the right panel shows those for Z4c.