Conformal and covariant Z4 formulation of the Einstein equations: strongly hyperbolic first-order reduction and solution with discontinuous Galerkin schemes

We present a strongly hyperbolic first-order formulation of the Einstein equations based on the conformal and covariant Z4 system (CCZ4) with constraint-violation damping, which we refer to as FO-CCZ4. As CCZ4, this formulation combines the advantages of a conformal and traceless formulation, with the suppression of constraint violations given by the damping terms, but being first order in time and space, it is particularly suited for a discontinuous Galerkin (DG) implementation. The strongly hyperbolic first-order formulation has been obtained by making careful use of first and second-order ordering constraints. A proof of strong hyperbolicity is given for a selected choice of standard gauges via an analytical computation of the entire eigenstructure of the FO-CCZ4 system. The resulting governing partial differential equations system is written in non-conservative form and requires the evolution of 58 unknowns. A key feature of our formulation is that the first-order CCZ4 system decouples into a set of pure ordinary differential equations and a reduced hyperbolic system of partial differential equations that contains only linearly degenerate fields. We implement FO-CCZ4 in a high-order path-conservative arbitrary-high-order-method-using-derivatives (ADER)-DG scheme with adaptive mesh refinement and local time-stepping, supplemented with a third-order ADER-WENO subcell finite-volume limiter in order to deal with singularities arising with black holes. We validate the correctness of the formulation through a series of standard tests in vacuum, performed in one, two and three spatial dimensions, and also present preliminary results on the evolution of binary black-hole systems. To the best of our knowledge, these are the first successful three-dimensional simulations of moving punctures carried out with high-order DG schemes using a first-order formulation of the Einstein equations.


I. INTRODUCTION
Large scale, fully general-relativistic numerical simulations have emerged in the last decade as a very powerful tool for the study of astrophysical systems, following the breakthrough calculations of the inspiral and merger of binary black holes [1][2][3].
The interest for such numerical techniques and the results they can produce has been only strengthened by the recent direct detection of gravitational waves [4], which paves the way for the forthcoming era of gravitational-wave astronomy.
General-relativistic simulations require (among other issues) stable and accurate methods for evolving the spacetime, i.e., for solving the Einstein field equations. The development of hyperbolic formulations of the Einstein equations that allow for longterm simulations of generic spacetimes, including the ones encompassing the physical singularities arising in the presence of black holes, has been therefore of great importance in numerical relativity. The first step in this direction has been the derivation of the Arnowitt-Deser-Misner (ADM) formulation (originally introduced in [5], but see [6][7][8][9][10][11] for a more modern perspective). While this formulation splits time and space and naturally presents general relativity as an initial boundary-value problem, suitable for numerical implementation, it is known to be not hyperbolic -at least when usual gauge choices are considered (see [12] for a discussion) -and therefore unstable in numerical applications.
Subsequently, a lot of effort has been devoted to find hyperbolic formulations of the Einstein equations. These efforts have lead to the derivation of the Baumgarte-Shapiro-Shibata-Nakamura-Oohara-Kojima (BSSNOK) formulation [13][14][15][16], which achieves hyperbolicity via a conformal transformation of the 3-metric and the promotion of some contractions of the Christoffel symbols to independently evolved variables and, most importantly, by inserting the momentum and Hamiltonian constraint expressions in the evolution system. A general-covariant alternative is the Z4 formulation of [17][18][19], which has been presented both in firstand second-order form in the spatial derivatives. More successful have been formulations based on the Z4 one that include a conformal transformation of the metric. These are the Z4c formulation, that removes some source terms in the Einstein equations in order to bring the evolution equations into a form which is closer to the BSSNOK system [20], and the conformal and covariant CCZ4 formulation [21,22] (see also [23,24] for some recent and slight variants). The Z4 family of formulations also admits mechanisms to damp constraint violations as they arise during the evolution [21,22,25,26].
For completeness, it should be mentioned that the 3+1 formalism is not the only way to develop a formulation of the Einstein equations suitable for numerical implementation: alternatives are the generalized-harmonic formalism e.g., [1,[27][28][29]; the characteristic-evolution formalism [30] the conformal approach [31,32] and fully-constrained formulations [33]. These approaches, however, are not the subject of the present work.
Parallel to the quest for better formulations of the equations, the development and implementation of better numerical methods has been a main priority of ongoing research. While most general-relativistic codes use finite-differences (e.g., [34]) or spectral methods (e.g., [35]) for the spacetime evolution, increasing interest is being focused towards DG methods (see, e.g., [36] for an introduction and review). DG methods are very attractive due to their excellent scalability and wave-propagation properties. The latter allow the propagation of smooth linear and nonlinear waves over long distances with little dissipation and dispersion errors, and should thus be in principle particularly well suited for the solution of the Einstein equations, where (apart from physical singularities in black holes) the fields are smooth and high accuracy can be achieved. So far, however, only a rather limited number of attempts have been made to solve the Einstein equations with DG methods. Field et al. [37] tested a second-order BSSNOK formulation, while Brown et al. [38] developed a first-order formulation of BSSNOK, however both works were limited to spherical symmetry and vacuum spacetimes. The first DG implementation in non-vacuum spacetimes was published by Radice & Rezzolla [39], but was still restricted to spherical symmetry. The first threedimensional (3D) implementation, albeit in a fixed spacetime and focused on hydrodynamics was developed by Bugner et al. [40]. More recently, Miller and Schnetter [41] proposed an operator-based DG method suitable also for second-order systems and applied it to the BSSNOK system, while Kidder et al. [42] developed a task based relativistic magnetohydrodynamics code.
In this work we propose a novel first-order form of the CCZ4 system, which we refer to as FO-CCZ4. We thoroughly study its eigenstructure and in particular show that it is strongly hyperbolic for two typical choices of gauges, namely zero shift with harmonic lapse and the Gamma-driver with 1+log slicing. We then implement this formulation in a fully three-dimensional code, using an ADER-DG algorithm with adaptive mesh refinement (AMR) and local time-stepping (LTS), supplemented with a high order ADER-WENO [43][44][45] finite-volume subcell limiter [46][47][48] to deal with singularities in black-hole spacetimes. This family of schemes has already been successfully applied also to the classical and special relativistic MHD equations (see [49,50]).
We test the stability and accuracy of the ADER-DG discretization applied to our novel FO-CCZ4 formulation in a series of standard tests for general-relativistic codes [51,52]. We also verify that our scheme converges at the expected order of accuracy and we provide evidence of long-time robustness and stability. Finally we apply the method to the long-term evolution of single black-hole spacetimes, showing that we are able to stably evolve a puncture black-hole spacetime for a time scale of ∼ 1000 M (M being the mass of the black hole). We also present preliminary results for the head-on collision of two black holes. To the best of our knowledge, these are the first simulations of black-hole spacetimes performed in three spatial dimensions with a high-order DG code.
The structure of the paper is as follows. In Section II we derive the full set of first-order evolution equations and prove the strong hyperbolicity for common gauge choices by computing the full eigenstructure. In Section III we introduce the numerical scheme intended to solve the partial differential equations (PDE) system. In Section IV we show a number of benchmark results to demonstrate correctness of both the formulation and the numerical solver. Finally, the conclusions are summarised in Section V.
We work in a geometrized set of units, in which the speed of light and the gravitational constant are set to unity, i.e., c = G = 1. Greek indices run from 0 to 3, Latin indices run from 1 to 3 and we use the Einstein summation convention of repeated indices.

II. A FIRST-ORDER STRONGLY HYPERBOLIC CCZ4 SYSTEM: FO-CCZ4
A. The original second-order CCZ4 system The second-order CCZ4 system can be derived from the Z4 Lagrangian L = g µν (R µν +2∇ µ Z ν ), which adds terms dependent on the Z µ vector to the classical Einstein-Hilbert Lagrangian (see [7] for a complete derivation). The Einstein field equations are recovered by minimizing the corresponding action and the algebraic constraints Z µ = 0. Additional constraint-damping terms can be introduced [25], so that the Einstein equations of the constraint-damped Z4 system in vacuum read where R µν is the Ricci tensor and n is the unit vector normal to the spatial hypersurfaces. Here, κ 1 , κ 2 are tuning constants related to the characteristic time of the exponential damping the of constraint violations. In order to formulate a well-posed Cauchy problem, we apply the 3+1 decomposition of space time (see, e.g., [7,8]), so that the line element reads with lapse α, shift β i and 3-metric γ ij . The 3+1 split leads to evolution equations for γ ij as well as the extrinsic curvature K ij = − 1 2 L n γ ij , L n being the Lie derivative along n µ ; because of the gauge freedom of general relativity, the functions α and β can be in principle freely specified. The four constraint equations of the ADM system (generally formulated as an elliptic system, but see, e.g., [53] for an alternative formulation) become four evolution equations for the Z µ vector.
The CCZ4 formulation, as presented in [21], introduces the conformal factor φ := (det(γ ij )) −1/6 to define the conformal 3-metricγ ij := φ 2 γ ij , with unit determinant. As in the BSSNOK system, the extrinsic curvature is decomposed into its trace K = K ij γ ij and a trace-free partÃ ij , which are promoted to primary evolution variables i.e., The second-order version of the vacuum CCZ4 equations, including the evolution equations for the 1 + log slicing [Eq. (4g)] and Gamma-driver shift condition [Eqs. (4h)-(4i)], is reported here for clarity, using essentially the same notation as in [21] with the contracted Christoffel symbolsΓ i :=γ jkΓi jk =γ ijγkl ∂ lγjk , the shorthandΓ i :=Γ i + 2γ ij Z j , and the use of the upper index TF to indicate a quantity whose trace has been removed.
We recall that the four-vector Z µ is an extra dynamical field specifically introduced to account for the energy and momentum constraints of the Einstein equations [17,18,54]. Its temporal component is Z 0 = Θ/α and the indices of its spatial part may be raised and lowered with the spatial physical metric γ ij . Following [21], the Hamiltonian constraint H and the momentum constraint M i of the CCZ4 system read as usual, namely where of course H = 0 = M i in the continuum limit.

B. Introduction of the auxiliary variables and resulting ordering constraints
We introduce the following 33 auxiliary variables, which involve first spatial derivatives of the metric terms, An immediate consequence of (6) and the Schwarz theorem on the symmetry of second-order derivatives are the following second order ordering constraints [55], which read: SinceÃ ij is by construction trace-free, the following additional constraint holds:γ ijÃ ij = 0, and thus These relations will be important later in order to derive a strongly hyperbolic system in first-order form. Furthermore, from the constraint det(γ ij ) = 1 and via the Jacobi formula ∂ k det(A) = tr(det(A)A −1 ∂ k A) on the derivatives of the determinant of a matrix, we obtain the following additional algebraic constraints on the auxiliary variables D kij (see also [38]) From Eq. (9), another differential constraint follows, namely, In practical implementations, however, we have not found particular benefits from making use of this additional constraint in the FO-CCZ4 formulation. The evolution equations for the auxiliary quantities are obtained by applying the temporal derivative operator ∂ t to equations (6), by subsequently exchanging the spatial and temporal derivatives on the right-hand side of the resulting equations and by making use of the PDEs (4a), (4c), (4g) and (4h).
Many different first-order formulations of the CCZ4 system are possible, since any non-purely algebraic term in the original second-order system can be written as a combination of conservative terms and non-conservative products (see [55,56] for a parametric study of such families of systems). In this work, we considered the two extreme cases: a first one, where as many terms as possible are written in a conservative flux-divergence form (see, e.g., [19], as an example for the first-order Z4 system) and a second formulation, similar to the ideas outlined in [6], making maximum use of the first-order ordering constraints, so that the variables defining the 4-metric (α, β i , φ andγ ij ) are only evolved by a nonlinear system of ordinary differential equations (ODEs) and where the rest of the dynamics is written in terms of non-conservative products. The coefficients of these nonconservative products are only functions of α, β i , φ andγ ij and no differential terms in these variables appear. The dynamical variables of the FO-CCZ4 system with Gamma-driver shift condition are then:Ã ij , K, Θ,Γ i , b i (the b i vector is an auxiliary field used to write the Gamma-driver gauge condition [6,21]) and the auxiliary variables A k , B i k , P k and D kij . In this paper we will follow the second approach, i.e., the final system of 58 evolution equations will consist of 11 ODEs and 47 PDEs and will have a very special structure discussed later in Section II D.
C. Strongly hyperbolic first-order form of the CCZ4 system The most natural first-order formulation of the CCZ4 system is non-conservative and appears in the following form discussed later in more detail where one has the state vector Q, the system matrices A i and the purely algebraic source terms S(Q). To obtain a strongly hyperbolic first-order system from the second-order CCZ4 formulation of Alic et al. [21] given by (4a)-(4i) we systematically use the constraints (7) and (8) and make maximum possible use of the auxiliary variables Eq. (6). In other words, our first-order CCZ4 system does not contain any spatial derivatives of α, β i ,γ ij and φ any more, but all these terms have been moved to the purely algebraic source term S(Q) by using (6). This has the immediate consequence that the evolution equations (12a), (12b), (12c) and (12d) reduce to ordinary differential equations instead of partial differential equations. Our final non-conservative first-order CCZ4 system reads as follows: with the PDEs for the auxiliary variables given by: Indicated in red in the equations above are those terms that have been added to the PDE to obtain an approximate symmetrization of the sparsity pattern of the system matrices (see discussion in Sec. II D and Fig. 1). A few remarks should be made now. First, the function g(α) in the PDE for the lapse α controls the slicing condition, where g(α) = 1 leads to harmonic slicing and g(α) = 2/α leads to the so-called 1 + log slicing condition, see [57]. Second, in order to obtain the advective terms along the shift vector in the evolution equations of the auxiliary variables, we have used the identities (7). We stress that it is important to use the second-order ordering constraints (7) in an appropriate way to guarantee strong hyperbolicity, since a naive first-order formulation of the second-order CCZ4 system that just uses the auxiliary variables in order to remove the second-order spatial derivatives will only lead to a weakly hyperbolic system (see [55] for a detailed discussion on the use of second-order ordering constraints in second order in space first order in time hyperbolic systems). Third, we have found that the use of first and second-order ordering constraints alone is not enough, but that one must also literally derive the PDE (12l) for D kij from (4a) by explicitly exploiting the fact thatÃ ij is trace-free via the use of the constraint T k by adding Eq. (8) to Eq. (12l). Without the use of T k in Eq. (12l), the system immediately loses its strong hyperbolicity (see also [58] for a similar observation in the Z4c system). Once again, these important additional terms in the FO-CCZ4 system related to the constraints (7) and (8) have been highlighted in red in Eqs. (12a)-(12m).
We also have introduced several additional constants compared to the original second-order CCZ4 system. In particular: • the constant τ is a relaxation time to enforce the algebraic constraints on the determinant ofγ ij and on the trace ofÃ ij "weakly" (see the discussion in [21]).
• the constant e is a cleaning speed for the Hamiltonian constraint, following the ideas of the generalized Lagrangian multiplier (GLM) approach of Dedner et al. [59]. As the cleaning is a non-physical process, e > 1 is in principle allowed; this leads to faster constraint transport and thus can be used to obtain a better satisfaction of the constraints for purely numerical purposes, but e = 1 breaks the covariance of the FO-CCZ4 system.
• the constant µ > 0 appears in Eq. (12k) and allows one to adjust the contribution of second-order ordering constraints.
• the constant s contributes to the evolution equations for b i , β i and B i k and allows to turn on or off the evolution of the shift. For s = 0 we have the simple gauge condition ∂ t β i = 0, while for s = 1 the usual Gamma-driver gauge condition is obtained. discrete level. Note also that when treating black holes as punctures, the lapse would vanish at the puncture location and its logarithm diverge. We therefore impose a positive lower limit in our numerical implementation. Since we employ a DG scheme where the solution in every element is represented by an interpolating polynomial (see section III), in an element surrounding the puncture the polynomial might actually reach values lower than the limit due to the Runge phenomenon; even in this case, however, the logarithm would not diverge. Furthermore, we have the following expressions and identities for various terms appearing in the evolution equations: Here, we have again made use of the second-order ordering constraints (7) by symmetrizing the spatial derivatives of the auxiliary variables as follows: Many of the above terms will contribute to the purely algebraic source term, as well as to the non-conservative product. For example, in the spatial derivatives of the Christoffel symbols of the conformal metric (16), the first bracket contributes only to the purely algebraic source term, while the second bracket is a non-conservative product. In a practical implementation, it is therefore necessary to carefully separate each contribution. We also stress that in our FO-CCZ4 formulation, the Ricci tensor R ij is directly calculated from the Riemann tensor R m ikj and the Christoffel symbols and their derivatives ab definitionem, without making use of the typical splitting of the Ricci tensor as e.g., used in [21]. We also compute the contracted Christoffel symbolsΓ i directly from their definition, without making use of the fact that the determinant ofγ ij is unity, since in general this cannot be guaranteed to hold exactly at the discrete level, unless the algebraic constraints are rigorously enforced.
From a more formal and mathematical point of view, the additional use of the second-order ordering constraints (7) and the constraint T k (the terms colored in red) can be motivated by looking at the structure of the sparsity pattern of the system matrix A · n = A 1 n 1 + A 2 n 2 + A 3 n 3 with and without the use of these constraints. In Fig. 1 we report the sparsity pattern of the system matrix in the normal direction n = 1/ √ 3(1, 1, 1) for the Gamma-driver shift condition and the 1 + log slicing condition for a randomly perturbed flat Minkowski spacetime, neglecting all matrix entries whose absolute value is below a threshold of 10 −7 . The blue dots represent the original sparsity structure without the use of the second-order ordering constraints (7) and without using the constraint (8), while the combination of the blue and the red dots shows the sparsity pattern after the terms colored in red have been added to the PDE system. Our approach for finding a suitable form of the ordering constraints to be added is based on approximate symmetrization of the sparsity pattern of the system matrix, in order to avoid Jordan blocks, which cannot be diagonalized. Such Jordan blocks are evident in the sparsity pattern given by the blue dots alone in Fig. 1.
We are not aware of works in which the constraint T k has been used in conformal first-order hyperbolic formulations of the 3+1 Einstein equations, but its effect becomes rather clear from Fig. 1. It is also directly evident from Fig. 1 that the first 11 quantitiesγ ij , α, β i and φ are only evolved by ODEs and that the entire system does not depend on spatial derivatives of these variables, since all entries in the first 11 rows and columns of the system matrix are zero. However, we explicitly stress here that our FO-CCZ4 system is not symmetric hyperbolic in the sense of Friedrichs [60], like for example the family of symmetric hyperbolic and thermodynamically compatible systems of Godunov and Romenski [61][62][63]. Further work in this direction will be necessary to try and achieve a symmetric hyperbolic form of FO-CCZ4 with a convex extension.
In summary and as an aid to the reader, we list below the key ideas that have been used in order to obtain the strongly hyperbolic FO-CCZ4 system proposed in this paper: 1. maximum use of the first-order ordering constraints (6) in order to split the complete system into 11 pure ODEs (29) for the evolution of the quantities defining the 4-metric (α, β i ,γ ij and φ), and with no spatial derivatives of these quantities appearing in the remaining PDE system (30). However, if we want to keep this very particular split structure of the PDE system, it is not possible to add damping terms proportional to the first-order ordering constraints (6) to the system, since this would make spatial derivatives of α, β i ,γ ij , φ appear again and may eventually lead to Jordan blocks which cannot be diagonalized. We therefore explicitly refrain from adding these terms, in contrast to what has been done in [38]. Following the philosophy above, also writing the system in a flux-conservative form like in [19,64] is not possible, since the fluxes will in general depend on the 4-metric and thus, after application of the chain rule, spatial derivatives of α, β i ,γ ij and φ would appear again in the quasi-linear form. We note that not adding any damping terms proportional to the first-order ordering constraints (6) may lead to a rapid growth of these constraints on the discrete level (see also [28]). This effect, however, may be reduced by a periodic reinitialization of the auxiliary variables with appropriate discrete versions of Eq. (6), either after a certain number of timesteps, or if a large growth of the first-order constraint violations is detected. However, in this paper this option has not been further investigated.
2. approximate symmetrization of the sparsity pattern of the system matrix A · n by appropriate use of the second-order ordering constraints (7) and the constraint (8), i.e., by adding the terms highlighted in red in PDEs (12a)-(12m). Symmetrization of the first derivatives of the auxiliary variables by using (28), apart from the advective terms along the shift vector.
3. introduction of an independent constraint propagation speed e for the Hamiltonian constraint H in the PDE (12g) for the variable Θ, following the ideas of the GLM approach of Dedner et al. [59].
4. use of the logarithms of α and φ as evolution variables, in order to guarantee positivity for α and φ in a simple and natural way. These evolution quantities are consistent with the definitions of the auxiliary variables A k and P k .

D. Strong hyperbolicity
As already shown briefly above, the FO-CCZ4 system (12a)-(12m) can be written in compact matrix-vector form (11), where the complete state vector is given by k , D kij , P k , containing a total of 58 variables that have to be evolved in time. For clarity, we show the full sequential form of all 58 variables in vector Q in Appendix A. The vector Q can be split as Q T = (V T , U T ) into a vector V of the 11 quantities that define the 4-metric, V T := (γ ij , ln α, β i , ln φ), and a vector U of the remaining 47 dynamic variables Fig. 1 it is obvious that the vector V is evolved in time only via ODEs of the type where S (Q) contains the first 11 elements of the vector of purely algebraic source terms S(Q). Therefore, the eigenvalues associated with the ODE subsystem for V are trivially zero. Since in our formulation of the FO-CCZ4 system we have made maximum use of the first-order ordering constraints, Eqs. (12a)-(12m) do not contain any spatial derivative of the quantities in V , so that the columns in the matrices of the related eigenvectors are trivially the unit vectors. The remaining reduced system that needs to be analyzed contains the vector U of the dynamic quantities and has the very particular structure where the source term S (Q) contains the remaining elements of the source vector S(Q) and where the system matrices B i depend only on the vector V defining the 4-metric and do not depend on the vector U . The non-trivial eigenvectors of the FIG. 1: Sparsity pattern of the system matrix A · n with n = (1, 1, 1)/ √ 3 for randomly perturbed flat Minkowski spacetime using the Gamma-driver shift condition (s = 1) and 1 + log slicing (g(α) = 2/α), without the use of the constraints (7) and (8) (blue dots) and with the use of these constraints (blue & red dots). The achieved approximate symmetrization of the sparsity pattern is evident. Note also the complete absence of non-zero entries in the first 11 lines and columns corresponding to the variablesγij,α, β i and φ, which clearly highlights the special structure of our FO-CCZ4 system that can be split into a set of pure ODEs and a reduced PDE system, as discussed in Section II D.
complete system (11) can thus be obtained from those of the reduced system (30) by simply adding zeros corresponding to the quantities contained in V .
An immediate consequence of the very particular splitting of (11) into the ODEs (29) and the reduced PDEs (30) is that all waves appearing in the system (30) and thus in (11) are linearly degenerate (see [65] for a detailed discussion), since the eigenvalues λ i depend only on V and not on U and hence ∂λ i /∂Q · r i = 0, ∀ λ i . This also means that the FO-CCZ4 system cannot generate shock waves, since the formation of classical shock waves requires the compression of characteristics and thus the presence of genuinely nonlinear fields [10,65].
In order to prove strong hyperbolicity of the FO-CCZ4 system proposed in this paper, we compute the entire eigenstructure of the system matrix B 1 in the x 1 direction for two standard gauge choices: i) zero shift β i = 0 (hence s = 0) with harmonic slicing, i.e., g(α) = 1 and ii) the gamma driver shift condition (s = 1) with 1+log slicing, i.e. g(α) = 2/α. Note that, in principle, the eigenstructure of the principal symbol of the system should be computed for every normal direction n = 0 in space. However, this is not necessary in this case, since the Einstein equations are isotropic (see [12]).
For the first shift condition, there is no need to evolve the quantities b i and B i k , whose corresponding PDEs can therefore be neglected in the following analysis (the associated eigenvalues are simply zero and the eigenvectors are the unit vectors). For zero shift the vector U can thus be furthermore reduced to only 35 remaining dynamic quantities U T = (Ã ij , K, Θ,Γ i , A k , D kij , P k ). In this case the 35 eigenvalues of matrix B 1 in x 1 direction are The associated complete set of 35 right eigenvectors defining the right eigenvector matrix R as well as the inverse right eigenvector matrix (L = R −1 ) that defines the so-called left eigenvectors are given in the first section of the Appendix A.
The fact that the FO-CCZ4 system has only real eigenvalues and a complete set of linearly independent eigenvectors (where the matrix of eigenvectors is uniformly bounded) is a necessary and sufficient condition for strong hyperbolicity. Note that for harmonic lapse the eigenvectors r 22,23 are only linearly independent of r 24,···35 if c = 1, ∀e > 0 or for e = 1, ∀c ≥ 0. The choice c = 1 and e = 1 corresponds to the standard setting typically used for second order Z4 and CCZ4 systems, and the importance of using c = 1 has already been shown in the hyperbolicity analysis for the first and second order Z4 system carried out in [18,54], i.e. our results on the FO-CCZ4 system confirm previous findings made in the literature. For the gamma driver shift condition, the hyperbolicity analysis is much more complex and requires the computation of all 47 eigenvectors of the reduced dynamical system (30), this time including also the quantities b i and B i k . After tedious calculations it was possible to obtain analytical expressions for the eigenvalues and all 47 eigenvectors also in this case. The results are reported in the second section of the Appendix A. To the best of our knowledge, this is the first time that a hyperbolicity analysis of a first-order reduction of the CCZ4 system including the gamma driver shift condition has been carried out. An analysis of the FO-CCZ4 system with other shift conditions, such as the generalized harmonic shift [54,66], is left to future work.
At this point, we would like to add the following clarifying remark. The hyperbolicity analysis has been carried out for the FO-CCZ4 evolution system (12a)-(12m), which in principle admits violations of the algebraic constraints det(γ ij ) = 1, γ ijÃ ij = 0 andγ ij D kij = 0. Hence, compared to the original Z4 system [17][18][19], it has an augmented solution space. Since our hyperbolicity analysis has been made without enforcing the algebraic constraints, it is valid for the FO-CCZ4 system with the augmented solution space, but should not be regarded as an analysis of the original Z4 system. However, if the initial data satisfies the algebraic constraints, a direct consequence of the system (12a)-(12m) is that the constraints will remain satisfied for all times, so that our hyperbolicity analysis also covers solutions that satisfy the algebraic constraints.

Unlimited ADER-DG scheme and Riemann solvers
As mentioned above, the FO-CCZ4 system (12a)-(12m) above can be written as a non-conservative first-order hyperbolic system of the symbolic form given by Eq. (11) (see also [46][47][48]), where the matrices A i are the system matrices in the coordinate direction x i and their eigenstructure has been analyzed in the previous section.
When solving numerically the system (11), the computational domain Ω is covered by a set of non-overlapping Cartesian tensor-product elements indicates the barycenter of cell Ω i and ∆x i = (∆x i , ∆y i , ∆z i ) defines the size of Ω i in each spatial coordinate direction. Furthermore, the domain Ω is the union of all elements, i.e., Ω = Ω i . Adaptive mesh refinement (AMR) has been implemented in a cell-by-cell framework based on a tree structure [67], together with time-accurate local time-stepping (LTS; see Refs. [49,50,[68][69][70] for details). In the DG finite-element framework, the discrete solution of the PDE system (11) is denoted by u h in the following and is defined in the space of tensor products of piecewise polynomials of degree N in each spatial direction, denoted U h in the following. At time t n , in each element Ω i the discrete solution is written in terms of some spatial basis functions Φ l (x) and some unknown degrees of freedomû n i,l as follows, where l := (l 1 , l 2 , l 3 ) is a multi-index and the spatial basis functions Φ l (x) = ϕ l1 (ξ)ϕ l2 (η)ϕ l3 (ζ) are generated via the tensor product of one-dimensional basis functions ϕ k (ξ) on the reference element [0, 1]. The mapping from physical coordinates For the one-dimensional basis functions ϕ k (ξ) we use the Lagrange interpolation polynomials passing through the Gauss-Legendre quadrature nodes ξ j of an N + 1 point Gauss quadrature formula (see Fig. 2). Hence, the basis polynomials satisfy the interpolation property ϕ k (ξ j ) = δ kj , where δ kj is the usual Kronecker symbol. Due to this particular choice of a nodal tensor-product basis, the entire scheme can be written in a dimension-by-dimension fashion, where all integral operators can be decomposed into a sequence of one-dimensional operators acting only on the N + 1 degrees of freedom in the respective dimension.
To derive the ADER-DG method, we first multiply the governing equations (11) by a test function Φ k ∈ U h , identical to the spatial basis functions of Eq. (32). After that, we integrate over the spacetime control volume Ω i × [t n ; t n+1 ] and obtain t n+1 with dx = dx dy dz, i.e., we integrate over coordinate volumes rather than over physical volumes. Since the solution is discontinuous across element interfaces, the resulting jump terms have to be taken properly into account. This is done in our numerical scheme with the aid of the path-conservative approach, first developed by Castro and Parés in the finite-volume framework [71,72] and later extended also to the DG finite-element framework in [46,47,73]. In the ADER-DG framework, higher order in time is achieved with the use of an element-local spacetime predictor, denoted by q h (x, t), and which will be discussed later.
Using (32), integrating the first term by parts in time, taking into account the jumps between elements and making use of the local predictor solution q h instead of Q, the weak formulation (33) can be rewritten as (34) In (34), the first integral leads to the so-called "element mass matrix", which is diagonal for our choice of the basis and test functions, the second integral accounts for the smooth part of the discrete solution in the interior Ω • i = Ω i \∂Ω i of the element Ω i , the boundary integral accounts for the jumps across the element interfaces and the term on the right-hand side accounts for the presence of the purely algebraic source terms S. Following the path-conservative approach [48,71,72], the jump terms are defined via a path-integral in phase space between the boundary extrapolated states at the left q − h and at the right q + h of the interface as follows: with A · n = A 1 n 1 + A 2 n 2 + A 3 n 3 the system matrix in normal direction and where we have used the simple segment path In Eq. (35) Θ > 0 denotes an appropriate numerical viscosity matrix. According to [46][47][48], the path integral appearing in (35) is simply computed numerically via a Gauss-Legendre quadrature formula of sufficient order of accuracy. In this paper, we use one to three Gaussian quadrature points to approximate the path integral above. For a simple path-conservative Rusanov-type method [46,74], the viscosity matrix reads and where s max denotes the maximum wave speed found at the interface. Furthermore, In order to reduce numerical dissipation for the quantities evolved via ODEs, i.e., for α, β i ,γ ij and φ, in alternative to the Rusanov scheme we also employ the recently-proposed path-conservative Harten-Lax-van Leer-Einfeldt-Munz (HLLEM) method [75], which is a generalization of the original HLLEM method [76,77], and for which the jump terms on the element boundary read with Here, R * and L * are the rectangular matrices containing only the right and left eigenvectors of the internal waves associated with the eigenvalues Λ * that one wants to resolve exactly in the HLLEM Riemann solver. In our case, these internal waves are exactly the 11 stationary contact waves associated with the 11 ODEs for α, β i , φ andγ ij , hence their wave speed is zero and thus Λ * = 0 and δ * = I 11×11 . The associated 11 eigenvectors are the unit vectors, hence R * = I 58×11 and L * = I 11×58 .
With the left and right signal speeds simply chosen as s L = −s max and s R = +s max and with s max computed as in Eq. (37), the HLLEM scheme takes the same form of Eq. (35) with the viscosity matrix given by The choice of the approximate Riemann solver closes the description of the numerical scheme (34). Next, we will briefly describe the computation of the local spacetime predictor solution q h needed in Eq. (34) and (35).

Spacetime predictor
The element-local spacetime predictor solution q h (x, t) is computed from the known discrete solution u h (x, t n ) at time t n using a solution of the Cauchy problem "in the small", i.e., without considering the interaction with the neighbors, according to the terminology introduced by Harten et al. in [43]. In the ENO scheme of Harten et al. [43], and in the original ADER approach of Toro and Titarev [78][79][80], the so-called Cauchy-Kovalewski procedure was used. This procedure is very cumbersome and is based on local Taylor series expansions in space and time and where time derivatives are replaced by the known space derivatives at time t n by successively differentiating the governing PDE system with respect to space and time and inserting the resulting terms into the Taylor series. For an explicit example of the Cauchy-Kovalewski procedure applied to the three-dimensional Euler equations of compressible gas dynamics, see [81]. However, it is obvious that a highly complex PDE system as the FO-CCZ4 model (12a)-(12m) is not amenable to such an approach, which heavily relies on analytical manipulations of the PDE system. Therefore, we use an alternative local spacetime DG predictor [46,82,83], which only requires the pointwise computation of source terms and non-conservative products. The solution q h is expanded into a local spacetime basis with the multi-index l = (l 0 , l 1 , l 2 , l 3 ) and where the spacetime basis functions θ l (x, t) = ϕ l0 (τ )ϕ l1 (ξ)ϕ l2 (η)ϕ l3 (ζ) are again generated from the same one-dimensional nodal basis functions ϕ k (ξ) as before, i.e., the Lagrange interpolation polynomials of degree N passing through N + 1 Gauss-Legendre quadrature nodes. The spatial mapping x = x(ξ) is also the same as before and the coordinate time is mapped to the reference time τ ∈ [0, 1] via t = t n + τ ∆t. Multiplication of the PDE system (11) with a test function θ k and integration over Ω i × [t n , t n+1 ] yields the following weak form in space and time, which is different from (33), since now the test and basis functions are also time dependent: Using the local spacetime ansatz (41), Eq. (43) becomes an element-local nonlinear system for the unknown degrees of freedom q i,l of the spacetime polynomials q h . The solution of (43) can be easily found via a simple and fast converging, fixed-point iteration detailed e.g., in [46,69,84]. This completes the description of the unlimited ADER-DG scheme.

ADER-WENO finite-volume subcell limiter
In general the spatial metric is smooth, which justifies the use of the high-order unlimited ADER-DG scheme discussed in the previous sections. However, in the presence of black holes physical singularities appear, which can lead to numerical instabilities or can even lead to a failure of the computation. Following the ideas outlined in Refs. [49,50,69], we therefore supplement our ADER-DG scheme with a high-order ADER-WENO subcell finite-volume limiter, which is much more robust than the unlimited DG scheme, but which is at the same time still high-order accurate in space and time.
While in Refs. [49,50,69] a sophisticated a posteriori limiting strategy has been proposed, in this paper for simplicity we fix the limited cells a priori for the entire duration of the simulation, since we assume to know the location of the black holes. Future simulations with moving black holes will require a dynamic adjustment of the limiter, as suggested in [49,50,69]. In practice, each computational cell Ω i that has been marked for limiting is split into (2N + 1) 3 finite-volume subcells, which are denoted by Ω i,s and that satisfy Ω i = s Ω i,s (see Fig. 2). Note that this very fine division of a DG element into finite-volume subcells does not reduce the timestep of the overall ADER-DG scheme, since the Courant-Friedrichs-Lewy (CFL) coefficient of explicit DG schemes scales with 1/ (2N + 1), while the CFL of finite-volume schemes (used on the subgrid) is of the order of unity. The discrete solution in the subcells Ω i,s is represented at time t n in terms of piecewise constant subcell averagesū n i,s , i.e.,ū n i,s := These subcell averages are now evolved in time with a third-order ADER-WENO finite-volume scheme that looks very similar to the ADER-DG scheme (34), namely Here, we use again a spacetime predictor solution q h , which is now computed from an initial condition given by a WENO [45] reconstruction polynomial w h (x, t n ) computed from the cell averagesū n i,s via a multi-dimensional WENO reconstruction operator detailed in [68,85]. The values at the cell interfaces q − h and q + h are again computed as the boundary extrapolated values from the left and the right subcell adjacent to the interface.
To summarize our nonlinear WENO reconstruction briefly: for each subcell Ω i,s we compute several reconstruction polynomials w k h (x, t n ) requiring integral conservation of w k h on a set of different reconstruction stencils S k i,s , i.e., This system is solved via a constrained least-squares algorithm requiring at least exact conservation in the cell Ω i,s itself (see [85] for details). From the set of reconstruction polynomials w k h , the final WENO reconstruction polynomial w h is obtained by using a classical nonlinear weighted combination of the polynomials w k h (see [45,85]) as follows: where the oscillation indicators σ k are computed as usual from The small parameter in (47), which is only needed to avoid division by zero, is set to = 10 −14 and the exponent in the denominator is chosen as r = 8. The linear weights are λ 1 = 10 5 for the central stencil (i.e., k = 1), while all other stencils (i.e., k > 1) have linear weight λ k = 1. This choice corresponds also to the one made in [85].
In a practical implementation it is very convenient to write also the WENO reconstruction polynomials in terms of some reconstruction basis functions ψ l (x) as w h (x, t n ) = Ψ l (x)ŵ n l . In this paper, following [68], the basis functions Ψ l are defined exactly as the Φ l , i.e., as tensor products of Lagrange interpolation polynomials through the Gauss-Legendre quadrature nodes. For the limiter, we only use a piecewise quadratic reconstruction, leading to a nominally third-order accurate scheme. As already mentioned before, the predictor is computed according to (43), where the initial data u h (x, t n ) is replaced by w h (x, t n ) and the spatial control volumes Ω i are replaced by the subcells Ω i,s .
Once all subcell averagesū n+1 i,s inside a cell Ω i have been computed according to (45), the limited DG polynomial u h (x, t n+1 ) at the next time level is obtained again via a classical constrained least squares reconstruction procedure requiring and Gauss-Legendre basis at order N (N + 1 nodes) N S = 2N + 1 Finite Volume subcells N = 2: where (50) is a constraint and imposes conservation at the level of the control volume Ω i . This completes the brief description of the ADER-WENO subcell finite-volume limiter used in this paper. We should emphasize that all our attempts to simulate puncture black holes with the unlimited ADER-DG scheme described in Sec. IV E have failed after only a few timesteps, and only with the aid of the limiter described in this section we were able to carry out stable long-time evolutions of puncture black holes.

Summary of the path conservative ADER-DG method with subcell ADER-WENO finite volume limiter
For the sake of clarity, in this section we briefly summarize the description of the ADER-DG method with subcell finitevolume limiter outlined in the previous sections. In particular, to illustrate the algorithm flow and its practical implementation, we list below the various stages of the algorithm over a full timestep.
For each element Ω i , the algorithm to obtain the solution at time t n+1 from the data at time t n proceeds as follows: 1. at all Gauss-Legendre quadrature nodes in each element, the data at time t n for all evolved variables are present, as computed from the previous timestep of the scheme (or from the initial data); if the cell Ω i is flagged as troubled (the subcell limiter is active), also all the subcell finite volume averages are present; 2. the discrete data are modified to strictly enforce the algebraic constraints of the CCZ4 system (see section IV); 3. the degrees of freedomq i,l of the spacetime predictor are computed from (43) using a fixed-point iteration method; the spatial derivatives ∇q h of the solution needed in the predictor are computed by differentiating the DG polynomial (a so-called pseudo-spectral derivative). In unlimited cells, the initial data in (43) are given by the DG polynomials u h , in limited cells the initial data are given by the reconstruction polynomials w h obtained via a nonlinear WENO reconstruction operator acting on the subcell averages; inside the spacetime predictor, no information from the neighbouring elements is necessary; 4. the predictor is used to compute the volume and boundary integrals of the explicit fully discrete corrector stage, i.e. Eq. (34) for the ADER-DG scheme and Eq. (45) for the ADER finite volume limiter; in both cases the jump terms at the element interfaces are evaluated using the path-conservative approach of Parés and Castro [71,72], see Eq. (35); 5. the volume and boundary integrals are used to compute the solution at time t n+1 in a fully discrete one-step time-update, reminiscent of the forward Euler method; 6. finally, if a cell is flagged for limiting, the finite volume subcell averages are used to reconstruct the limited DG polynomial u h (x, t n+1 ).
The algorithm is simplified by the assumption that the cells to be limited and evolved via the ADER-WENO subcell finitevolume scheme rather than the ADER-DG scheme, are fixed a-priori throughout the entire simulation. Otherwise, additional steps would be present in order to determine if a particular cell is developing pathologies in the numerical solution and the limiting procedure should be activated, or not. Furthermore, in this work we do not employ a dynamic adaptive mesh-refinement (AMR) approach; instead, for simplicity the refined grid structure is fixed at the beginning of the simulation. Note finally that the spacetime predictor is computed "in the small", disregarding contributions from the neighbouring cells. This means in particular that also boundary conditions are not supplied to the spacetime predictor, but only to the corrector scheme which carries out the fully-discrete time update of each cell. Combined with the compact stencil of the ADER-DG method at any order of accuracy, which involves only the cell and its direct neighbors, this potentially allows for a very efficient parallel implementation of the algorithm.

IV. NUMERICAL TESTS
In the following we present a battery of standard tests that explore the ability of our formulation to carry out long-term stable evolutions of a number of different spacetimes with increasing degree of curvature. If not stated otherwise, in all of the tests we set initially Θ = 0,Γ i =Γ i and b i = 0 and the HLLEM method is used, i.e., Eq. (35) with the viscosity matrix (40). In all tests the algebraic constraints on the unit determinant ofγ ij , the zero trace ofÃ ij as well as the constraintγ ij D kij = 0 (which is a consequence of |γ ij | = 1) have all been rigorously enforced in the discrete solution u h (x, t n ) at the beginning of each timestep, but they have not been enforced during the computation of the spacetime predictor q h . Note that the predictor q h is only an auxiliary quantity that is overwritten after each timestep and which has a role similar to the evolution stage to the half timelevel in second-order MUSCL-Hancock type TVD finite-volume schemes. We therefore set τ → ∞ and thus neglect the corresponding source terms. In tests involving black holes, the lower limit on the lapse is set to be ln(α) ≥ −20. We will use the notation P N to indicate an ADER-DG scheme using piecewise polynomials of degree N to represent u h .

A. Linearized gravitational-wave test
The first test problem is a simple one-dimensional wave-propagation test problem in the linearized regime. The computational setup follows the one suggested by Alcubierre et al. in Ref. [86]. The computational domain is Ω = [−0.5, 0.5] with periodic boundary conditions in the x direction and two simulations are run until a final time of t = 1000: (i) a first one using 4 ADER-DG P 5 elements (i.e., a total number of 24 degrees of freedom) and (ii) a second one using only 2 ADER-DG P 9 elements (i.e., only 20 degrees of freedom). This test is run with the unlimited version of the ADER-DG scheme. The exact solution of the metric of the problem is given by and the wave amplitude = 10 −8 is chosen small enough in order to stay in the linear regime, so that terms O( 2 ) can be neglected. Since the shift is zero in the metric (51) (β i = 0), we set s = 0 in our FO-CCZ4 system and furthermore harmonic slicing is used, i.e., g(α) = 1. We also set K 0 = 0, c = 0, e = 2 and use the undamped version of the system, setting κ 1 = κ 2 = κ 3 = η = 0. Using the metric (51), the definition of the extrinsic curvature reduces to K ij = − 1 2 ∂ t γ ij /(α), so that the various components are given by K xx = K xy = K xz = K yz = 0, K yy = − 1 2 ∂ t h and K zz = + 1 2 ∂ t h. From this information, the conformal factor φ, the conformal spatial metricγ ij , the traceless conformal extrinsic curvatureÃ ij and all auxiliary variables can be computed by a direct calculation according to their definitions.
In Fig. 3 we report the temporal evolution of all ADM constraints (Hamiltonian and momentum constraints) as well as the errors of the algebraic constraints on the determinant of the conformal metric and the error in the trace ofÃ ij in both cases, i.e., using the ADER-DG P 5 and P 9 scheme. A comparison of the extrinsic-curvature componentÃ 22 with the exact solution is also provided at the final time t = 1000, showing overall an excellent agreement between numerical and exact solution. The quality of the results obtained with the ADER-DG schemes used in this paper, which are uniformly high-order accurate in both space and time, is significantly superior to the results shown in Ref. [86] for the same test problem using a finite difference scheme with much more grid points (between 50 and 200) compared to the very coarse mesh containing only 20 to 24 degrees of freedom used in our simulations. Note that a fair comparison between high order finite-difference and DG schemes must be made in terms of points per wavelength for finite-difference methods and in degrees of freedom per wavelength for DG schemes.

B. Gauge-wave test
Also this classical test problem has been taken from the collection of standard tests of Ref. [86]. The metric in this case is given by The metric (52) implies zero shift (β i = 0), hence we use once more s = 0 together with harmonic slicing g(α) = 1. Also for this test we employ the undamped version of the FO-CCZ4 system, setting κ 1 = κ 2 = κ 3 = η = 0. The computational domain in this case is two-dimensional, with Ω = [−0.5, 0.5] × [−0.05, 0.05] with periodic boundary conditions in all directions. Since β i = 0, the extrinsic curvature is again given by K ij = −∂ t γ ij /(2α), i.e., K yy = K zz = K xy = K xz = K yz = 0 and the remaining primary variables are We furthermore set K 0 = 0. The auxiliary variables can be obtained from their definition via a straightforward calculation. We first simulate this test problem with a perturbation amplitude of A = 0.1 until t = 1000 with an unlimited ADER-DG P 3 scheme and using 100 × 10 elements to cover the domain Ω. We run this physical setup twice, once with the default parameters e = c = 1, according to the original second order CCZ4 system [21] and a modified setting with e = 2 and c = 0 to obtain an improved cleaning of the Hamiltonian constraint. In both cases the system is strongly hyperbolic. The time evolution of the ADM constraints is reported in Fig. 4, showing only a very moderate growth of the constraint M 2 that is sublinear in time and close to machine precision. The other constraints H and M 1 remain essentially constant during the entire simulation. We emphasize that we have used the undamped version of the FO-CCZ4 system, and nevertheless obtain stable results, while the original second-order CCZ4 formulation was reported to fail for this test problem in the undamped version, and only the damped CCZ4 system was stable (see [21] for details). It is also worth recalling that both the first-and the second-order formulation of the BSSNOK system fail for this test case after a rather short time (see [21,38]). In Fig. 4 we also provide a direct comparison of the solution after 1000 crossing times for the conformal factor φ as well as for the trace of the extrinsic curvature K. Note the overall very good agreement between the numerical solution and the exact one. For the sake of clarity, in the plots of the waveforms we also report the numerical error computed as the difference between the numerical solution and the exact solution at the final time t = 1000. It can be clearly noticed from the computational results shown in Fig. 4 that the constraints and the phase errors in the waveforms are significantly smaller for the modified setting e = 2, which may justify the use of a faster cleaning speed of the Hamiltonian constraint e > 1 for purely numerical purposes. In any case, our FO-CCZ4 system behaves well also with the default setting e = c = 1, which is typically used in the standard second order CCZ4 system [21].
Since the gauge-wave test has a smooth nontrivial exact analytical solution and is also valid in the nonlinear regime of the equations, we can use it in order to perform a numerical convergence study. For this purpose, we run the test again with different   unlimited ADER-DG P N schemes on a sequence of successively refined meshes. To make the test more difficult, we choose a very large perturbation amplitude of A = 0.9, which takes the system in the highly nonlinear regime, although in the end the test consists only in a nonlinear re-parametrization of the flat Minkowski spacetime. For thise case we use c = 0 and e = 2. We set the final simulation time to t = 10 and continue using the undamped version of the FO-CCZ4 system. The L 2 error norms of the conformal factor φ, the lapse α and the trace of the extrinsic curvature K, together with the observed order of accuracy of the different ADER-DG schemes are reported in Table I. We observe essentially the expected order of accuracy of the scheme for N = 3 and N = 4, while a superconvergence is observed for N = 5 and N = 7. We think that this is due to the strong nonlinearities of the PDE system appearing in the regime in which we run this test case with A = 0.9 and that some leading errors may be dominated by quadratic terms in the metric and the conformal factor, which can lead to a faster error decay than N + 1 for coarse meshes. However, we expect that this superconvergence will disappear on sufficiently refined meshes; but since the absolute errors are already getting close to machine accuracy on the meshes used here, it is not possible to refine the mesh much more with double-precision arithmetics, at least in the N = 7 case. For the ADER-DG P 5 scheme using 100 × 10 elements a comparison between numerical and exact solution of the nonlinear waveforms for φ, α, K and D xxx is provided in Fig. 5 at t = 10, where we can note again an excellent agreement between exact and numerical solution.

C. Robust stability test
The so-called robust stability test is the last standard test problem that we take from Ref. [86]. While in the previous test problems we have used a simple frozen shift condition ∂ t β i = 0 by setting s = 0 in the FO-CCZ4 system, here we employ the classical Gamma-driver shift condition. Furthermore, we employ the 1 + log slicing condition, setting the slicing function to g(α) = 2/α and the parameter f of the Gamma driver to f = 0.75, which is also the typical value used for the BSSNOK system and for the classical second-order CCZ4 system (see [21] for details). We further set e = 2, κ 1 = κ 2 = κ 3 = 0, K 0 = 0, c = 1 and η = 0.
As customary in this test, we start from the flat Minkowski metric We then add uniformly distributed random perturbations to all quantities of the FO-CCZ4 system, i.e., to all primary and auxiliary variables and also to Θ andΓ i . The two-dimensional computational domain is Ω = [−0.5, 0.5] 2 and we run different simulations with an unlimited ADER-DG P 3 scheme on four successively refined meshes composed of 10ρ × 10ρ elements, corresponding to 40ρ × 40ρ degrees of freedom, where ρ ∈ {1, 2, 4, 8} is the refinement factor. The perturbation amplitude is = 10 −7 /ρ 2 , which corresponds to perturbation amplitudes that are three orders of magnitude larger that those suggested in Ref. [86].
The time evolution of the ADM constraints is reported in Fig. 6 for all four simulations. One can observe that after an initial decay the constraints remain essentially constant in time for all different grid resolutions, indicating that our FO-CCZ4 system indeed passes the robust stability test with the standard Gamma driver and 1 + log gauge conditions (see [58] for similar tests with the Z4c system).

D. Convergence tests on three-dimensional black-hole spacetimes
In this test we consider the evolution of isolated Schwarzschild and Kerr black holes in 3D Cartesian Kerr-Schild coordinates, with M = 1 the mass of the black hole and a the dimensionless spin. The metric in these coordinates is known analytically and thus the primary variables of our evolution system are given by x 2Hl x l y 2Hl x l z 2Hl x l y 1 + 2Hl 2 y 2Hl y l z 2Hl x l z 2Hl y l z 1 + 2Hl 2 with H := M r 3 r 4 + a 2 z 2 , S := 1 + 2H , l x := rx + ay r 2 + a 2 , l y := ry − ax r 2 + a 2 , l z := z r , and r := (x 2 + y 2 + z 2 − a 2 )/2 + ((x 2 + y 2 + z 2 − a 2 )/2) 2 + z 2 a 2 .
We furthermore use the fact that the solution is stationary, i.e., ∂ t γ ij = 0, hence the extrinsic curvature K ij is computed as follows (see [10]) The function K 0 is chosen as K 0 = K − β k ∂ k α / α 2 g(α) , so that ∂ t α = 0 and in this test the Gamma-driver shift condition is simplified to , with the consequence that the above exact solution corresponds to a stationary solution of the FO-CCZ4 system. In other words, we remove the advection terms from the evolution equations of the shift β i and the variable b i (see also [6]). The conformal factor φ and the auxiliary variables can be computed according

Schwarzschild black hole (a = 0)
Kerr black hole (a = 0.9) to their definition. The computational domain is chosen as Ω = [1, 5] 3 M 3 , and the exact solution given by the initial condition is imposed on all boundaries in all variables at all times. Note that this choice of boundary conditions is appropriate to study convergence since the exact solution is also a stationary solution of our PDE system. Note also that the black hole is centered at x = y = z = 0, so that we evolve only a section of the domain offset from the singularity, but encompassing regions both inside and outside of the event horizon; this effectively amounts to employing an excision of the black-hole interior. We furthermore set e = 2, c = 1, η = 0, and consider the undamped CCZ4 system with the 1 + log slicing, i.e., we set κ 1 = κ 2 = κ 3 = 0, f = 0.75 and g(α) = 2/α.
The simulations were performed with different ADER-DG schemes on a sequence of successively refined meshes until a final time of t = 10 M . The Rusanov method is used as approximate Riemann solver at the element interfaces. In the case of the Schwarzschild black hole we use a = 0, while for the Kerr black hole we set a = 0.9. The corresponding numerical convergence rates are reported for both cases in Table II, where we observe that the designed order of accuracy N + 1 of our high-order fully-discrete one-step ADER-DG schemes has been properly reached.

E. Evolution of a single puncture black hole
We next have applied the FO-CCZ4 formulation to a single puncture black hole [87] with mass M = 1 and dimensionless spin a = 0 located at the origin of a three-dimensional computational domain Ω = [−150, 150] 3 M 3 with periodic boundary conditions everywhere. The domain is discretized with an AMR mesh with grid spacing ∆x = ∆y = ∆z = 2.5 M within the inner box Ω b = [−15, 15] 3 M 3 , while ∆x = ∆y = ∆z = 7.5 M is used in the outer part of the domain. In the innermost zone Ω l = [−3, 3] 3 M 3 the third-order subcell ADER-WENO finite-volume limiter is activated throughout the entire simulation. For details on the AMR framework and the subcell finite-volume limiter we refer the interested reader again to [49,50,69]. We also stress that this simulation can be run only after activating the finite-volume subcell limiter, since a robust scheme is needed in order to deal with the puncture singularity. Without such a limiter, i.e., with a pure DG scheme, the code crashes after a few timesteps since the high-order unlimited DG scheme is not robust enough to deal with the puncture metric. In our simulation we use an ADER-DG P 3 scheme (N = 3), which leads to 2N + 1 = 7 finite-volume subcells per DG element, i.e., the effective mesh spacing in terms of points (cell averages) inside the domain Ω l is ∆x = ∆y = ∆z = 0.357 M . Note that we set up the mesh so that the puncture is located at the boundary of the DG elements; given the location of the degrees of freedom in the subcell grid (see Fig. 2), no grid point coincides with the puncture. We set the CCZ4 parameters to κ 1 = 0.1, κ 2 = 0, κ 3 = 0.5 and η = 0. The constant µ accounting for the second-order ordering constraints in the evolution of B i k is set to µ = 1/5, while for this test we use c = 1, f = 0.75 and e = 1 to be as close as possible to a standard second-order CCZ4 formulation, where the cleaning of the Hamiltonian constraint is done at the speed of light. The initial metric and lapse are provided by the TwoPunctures initial data code [88] (part of the Einstein Toolkit software [89]). Explicitly, the lapse is set initially to where r * := (r 4 + 10 −24 ) 1 4 and r is the coordinate distance of a grid point from the puncture. The auxiliary quantities (which are spatial derivatives of the primary quantities) are obtained via a simple fourth order central finite difference applied to the primary variables α and γ ij . Initially the shift and the extrinsic curvature are set to zero, i.e., β i = 0 and K ij = 0.
The evolution was carried out until a final time of t = 1000 M and Fig. 7 reports the evolution of the average L 2 error of the ADM constraints, which we define as where denotes the local error of each of the ADM quantities, i.e., Hamiltonian H and momentum constraints M i . In Fig. 7 also a view of the 3D grid setup is shown together with a zoom into the center region with the contour colors of the lapse function at a time of t = 200 M . It is probably worth recalling that, to the best of our knowledge, these are the first results obtained for a puncture black-hole spacetime using a fully three-dimensional DG finite-element method with AMR and LTS. Previous results obtained with highorder DG schemes for black-hole spacetimes were essentially limited to the one-dimensional case (see, e.g., Refs. [37,38,41]).

F. Preliminary results for moving punctures
The last test considered is a preliminary application of the FO-CCZ4 system to a binary system of two moving punctures. In particular, we consider a head-on collision of two nonrotating black holes of equal mass M = 1 with zero linear momentum initially located at x − = (−1, 0, 0) and x + = (+1, 0, 0). The three-dimensional computational domain is given by Ω = [−25, 25] 3 M 3 and flat Minkowski spacetime is imposed as boundary condition everywhere. The CCZ4 parameters are set to κ 1 = 0.1, κ 2 = 0, κ 3 = 0.5, η = 0 and furthermore we choose c = 1, e = 1, f = 1 and µ = 1/5. Again, the initial metric and the lapse are provided by the TwoPunctures initial data code [88], with the lapse set initially to where r * − and r * + are the coordinate distances of a grid point from either puncture (defined analogously to the previous section) and m − and m + are the so-called bare masses of the two black holes (see [88]) and in this case they are equal. The auxiliary quantities are computed from the primary variables via a fourth-order central finite-difference method. We use the simple and robust Rusanov method as approximate Riemann solver on the element boundaries. The shift and extrinsic curvature are initially set to β i = 0 and K ij = 0.
The domain is discretized with an AMR mesh of mesh spacing ∆x = ∆y = ∆z = 5/12 M within the inner box Ω b = [−2.5, 2.5] 3 M 3 , while ∆x = ∆y = ∆z = 1.25 M is used in the outer part of the domain. In the innermost zone Ω l = [−5/3, 5/3] 3 M 3 the third-order subcell ADER-WENO finite-volume limiter is activated throughout the entire simulation. As for a single puncture, we use an ADER-DG P 3 scheme (N = 3), whose 2N + 1 = 7 finite-volume subcells lead to an effective mesh spacing inside the domain Ω l of ∆x = ∆y = ∆z = 0.0595. Once again we remark that the use of the finite-volume subcell limiter is essential in order to obtain a stable evolution.
The simulation is run until a final time of t = 60 M and the evolution of the contour surfaces of the lapse and the shift vector are reported in Fig. 8. The contour surfaces of the conformal factor at the final time as well as the evolution of the ADM constraints are depicted in Fig. 9. Clearly, no sign of growth in the violation of the constraints appears after the two punctures have merged at t 10 M .
Although these results are meant mostly as a proof-of-concept rather than as a realistic modelling of the inspiral and merger on binary black-hole systems, they provide convincing evidence that binary systems of puncture black holes can be evolved stably with our path-conservative ADER-DG scheme with ADER-WENO subcell finite-volume limiter on AMR grids based on the FO-CCZ4 formulation proposed here. A more detailed and systematic investigation, which includes the emission of gravitational waves from binary systems of rotating black holes in quasi-circular orbits (see, e.g., [21]), will be the subject of future work.

V. CONCLUSIONS
We have proposed a possible FO-CCZ4 formulation of the Einstein equations, that is, a strongly hyperbolic first-order formulation of the conformal covariant Z4 (CCZ4) system of Alic et al. [21]. The system comprises 58 evolution equations for the complete state vector given by To the best of our knowledge, this is the first time that a first-order hyperbolic formulation of the CCZ4 system has been proposed and has been employed in a systematic series of numerical tests of increasing complexity in one, two and three spatial dimensions.
The key idea for obtaining strong hyperbolicity in the new formulation is the approximate symmetrization of the sparsity pattern of the system matrix, that is, the appropriate use of ordering constraints and by using the fact that the trace ofÃ ij is zero, in order to avoid the appearance of Jordan blocks that cannot be diagonalized. Another important idea employed to obtain the FO-CCZ4 formulation is the use of first-order ordering constraints in a way that reduces the evolution equations of the lapse α, the shift β i , the conformal metricγ ij and the conformal factor φ to a pure system of ordinary differential equations [6]. In other words, whenever differential terms with respect to α, β i , φ andγ ij appear, they are replaced by the corresponding auxiliary variables A k , B i k , P k and D kij and thus become algebraic source terms. This leads to a very particular split structure of the system that also greatly simplifies the analysis of the resulting FO-CCZ4 system, since the eigenvalues and eigenvectors associated with α, β i , φ andγ ij become trivial, being zero and unit, respectively. This has also the advantage that for the rest of the analysis a reduced system of partial differential equations relative to only 47 dynamic variables U T = (Ã ij , K, Θ,Γ i , b i , A k , B i k , D kij , P k ) can be considered. Furthermore, the matrix of the reduced system in the dynamic variables is only a function of α, β i , φ andγ ij , and not of the dynamic variables themselves, which not only substantially simplifies the hyperbolicity analysis but which also leads to the important result that all fields of our FO-CCZ4 system are linearly degenerate.
When compared to the first-order Z4 system proposed in Refs. [19,64], our entire FO-CCZ4 system is written in a fully non-conservative form, which is another key idea of our FO-CCZ4 formulation. We stress that the previously mentioned simplifications are not possible if a conservative formulation of the system based on the divergence of fluxes is used, e.g., like the one proposed in [19,64], since the Jacobian ∂F /∂Q of the flux F (Q) will also depend on the dynamical variables and the quasi-linear form of the system will also contain differential terms in α, β i , φ andγ ij , while in our genuinely non-conservative formulation, no differential terms of the latter quantities appear.
We have also provided a proof of strong hyperbolicity of our FO-CCZ4 system for two standard gauge choices, namely i) harmonic lapse and zero shift and ii) 1+log slicing combined with the Gamma-driver. In both cases we have computed the analytical expressions of all eigenvalues and all right eigenvectors of the system. For zero shift and harmonic lapse it was also possible to provide the inverse of the right eigenvector matrix, i.e. the so-called left eigenvectors of the system. To the best knowledge of the authors, this is the first time that the hyperbolicity of a first-order reduction of the CCZ4 system is analyzed, in particular including the Gamma-driver shift condition.
We have numerically solved the FO-CCZ4 system after discretizing it with the aid of a family of high-order fully-discrete one-step ADER discontinuous Galerkin (DG) schemes, supplemented with an ADER-WENO finite-volume limiter in order to deal with the physical singularities arising with black holes. The non-conservative formulation of the system is naturally treated within the framework of path-conservative schemes, first proposed by Castro and Parés in the finite-volume context [71,72] and later extended also to ADER-DG schemes in [46,47]. Furthermore, in order to ensure positivity of the numerical solution in terms of α and φ, we have decided to evolve the logarithms ln α and ln φ of these quantities in time, rather than the quantities themselves.
As customary for novel formulations of the Einstein equations, we have applied the strongly hyperbolic FO-CCZ4 system to a series of standard test cases suggested in Ref. [86], such as the gauge-wave test, the robust stability test and the linear-wave test bed. Besides providing evidence that the new system is able to reproduce the analytic solution accurately, we have carried out numerical convergence studies of the method on the gauge-wave test in the highly nonlinear regime of the equations, as well as on the Schwarzschild and a Kerr black hole using 3D Cartesian Kerr-Schild coordinates. We have also provided numerical evidence that our ADER-DG scheme with ADER-WENO finite-volume subcell limiter is able to perform a long time integration of a single puncture black hole with the usual Gamma driver and 1 + log gauge conditions that are typically used in simulations carried out with the BSSNOK evolution system. Finally, we have also shown some first preliminary results for two moving puncture black holes. To the best of our knowledge, the numerical results shown in this paper represent the first simulations of the 3+1 Einstein equations ever done with high-order DG and WENO finite-volume schemes on three-dimensional adaptive grids. All previous simulations of black-hole spacetimes with high-order DG schemes, in fact, were limited to the one-dimensional case only.