# Well-balanced high order finite difference WENO schemes for a first-order Z4
formulation of the Einstein field equations

Dinshaw Balsara Department of Physics and Astronomy, University of Notre Dame,
Notre Dame, IN, United States  ACMS Department, University of Notre Dame,
Notre Dame, USA  Deepak Bhoriya Department of Physics and Astronomy,
University of Notre Dame, Notre Dame, IN, United States  Olindo Zanotti
Laboratory of Applied Mathematics, Department of Civil, Environmental and
Mechanical Engineering, University of Trento, Via Mesiano 77, I-38123 Trento,
Italy  Michael Dumbser Laboratory of Applied Mathematics, Department of Civil,
Environmental and Mechanical Engineering, University of Trento, Via Mesiano
77, I-38123 Trento, Italy

###### Abstract

In this work we aim at developing a new class of high order accurate well-
balanced finite difference (FD) Weighted Essentially Non-Oscillatory (WENO)
methods for numerical general relativity, which can be applied to any
first–order reduction of the Einstein field equations, even if
non–conservative terms are present. We choose the first–order non–conservative
Z4 formulation of the Einstein equations, which has a built–in cleaning
procedure that accounts for the Einstein constraints and that has already
shown its ability in keeping stationary solutions stable over long timescales.
Upon the introduction of auxiliary variables, the vacuum Einstein equations in
first order form constitute a PDE system of 54 equations that is naturally
non-conservative. While prior versions of FD WENO schemes were developed
exclusively for hyperbolic systems of conservation laws, we show how to design
FD WENO schemes that can handle non-conservative products. Different variants
of FD–WENO are discussed, with an eye to their suitability for higher order
accurate formulations for numerical general relativity. We successfully solve
a set of fundamental tests of numerical general relativity with an order of
accuracy that has been implemented up to ninth order in space. Due to their
intrinsic robustness, flexibility and ease of implementation, finite
difference WENO schemes can effectively replace traditional central finite
differencing in any first–order formulation of the Einstein field equations,
without any artificial viscosity. When used in combination with _well
balancing_ , the new numerical schemes allow to preserve stationary
equilibrium solutions of the Einstein equations exactly, such as isolated
rotating and non-rotating black holes. This is particularly relevant in view
of the numerical study of the quasi-normal modes of oscillations of relevant
astrophysical sources. In conclusion, general relativistic high energy
astrophysics could benefit from this new class of numerical schemes and the
ecosystem of other desirable capabilities that has been built around them.

first-order Z4 formulation of the Einstein field equations — Numerical general
relativity — High order finite difference schemes — WENO methods — well-
balancing

##  1 Introduction

In recent years, many of the equations that are needed in computational
astrophysics have seen the development of higher order methods for their
solution. This harkens to the fact that the corresponding high accuracy is
needed in response to more accurate observations; General Relativity (GR) is
an emblematic example of this fact. Indeed, the third generation of
gravitational wave detectors (Lück et al., 2022), such as the Einstein
Telescope in Europe and the Cosmic Explorer in US, with a planned sensitivity
$h$ in the range $10^{-25}\div 10^{-24}\,\rm{Hz}^{-1/2}$, will soon call for
more accurate numerical solution methods in Numerical Relativity (NR).

As is well known, the Einstein equations present inbuilt challenges that do
not show up in many of the other governing equations commonly considered in
computational astrophysics. For one thing, the system naturally arises as a
second order PDE system, whereas our best numerical tools have been developed
for first order PDE systems. Possibly due to the great results made possible
by the second–order implementation of the BSSNOK formulation of the Einstein
equations (Shibata & Nakamura, 1995; Baumgarte & Shapiro, 2010), most of the
first–order analysis have instead remained at an exploratory or theoretical
level (Fischer & Marsden, 1972; Anderson & York, 1999; Kidder et al., 2001;
Bona et al., 2003; Alvi, 2002; Alekseenko & Arnold, 2003; Buchman & Bardeen,
2003; Bona & Palenzuela-Luque, 2005; Olivares et al., 2022). Notable but rare
examples of practical implementations of first order formulations into
numerical codes were provided by Bona et al. (2004a), Brown et al. (2012),
Buchman & Bardeen (2005), Dumbser et al. (2018). The latter, in particular,
applied Discontinuous Galerkin (DG) schemes to the conformal and covariant Z4
formulation of the Einstein equations (FO-CCZ4) of Alic et al. (2012). More
recently, Dumbser et al. (2024) have obtained very promising results after
using DG methods in a well–balanced first–order implementation of the damped
version of the Z4 formulation proposed by Gundlach et al. (2005), henceforth
referred to as the FO-Z4 formulation. We recall that the Z4 formulation of the
Einstein equations was originally proposed by Bona et al. (2003, 2004a) to
automatically account for a proper treatment of the Einstein constraints
through the addition of a four vector $z^{\mu}$, in a rather similar way to
what is done in the divergence cleaning approach by Munz et al. (2000) and
Dedner et al. (2002) for the Maxwell and magnetohydrodynamics (MHD) equations.
Later on, this approach has been combined with the conformal decomposition of
the metric, which is missing in the original Z4 formulation, to obtain the
so–called CCZ4 and Z4c formulations (Alic et al., 2012; Hilditch & Richter,
2012; Hilditch et al., 2012; Alic et al., 2013; Dumbser et al., 2018; Peterson
et al., 2023).

On the computational side, one can expect that, for various reasons, finite
difference (FD) schemes will remain the preferred choice by the NR community
still for many years to come. When a second–order formulation of the Einstein
equations is adopted, central finite difference schemes are the most natural
choice, which is indeed the case for such popular codes as Einstein–Toolkit
(Löffler et al., 2012), LazEv (Zlochower et al., 2005; Lousto & Healy, 2023),
BAM (Brügmann et al., 2008; Thierfelder et al., 2011), McLachlan (McLachlan,
2008), GRChombo (Clough et al., 2015), AMReX (Peterson et al., 2023), SACRA
(Yamamoto et al., 2008; Kiuchi et al., 2017). On the other hand, when a
first–order formulation of the Einstein equations is available, a natural
temptation arises: namely to migrate from central finite difference schemes to
finite difference Weighted Essentially Non-Oscillatory (WENO) methods, due to
their superior robustness in the presence of shock waves and singularities,
maintaining at the same time an excellent performance in terms of accuracy.
However, original WENO methods were specifically devised for first–order
hyperbolic systems in conservation form, and they do therefore not fit
straightforwardly into the non–conservative form of the first–order Einstein
equations. In fact, the usage of FD–WENO methods in the relativistic context
has been so far limited to the solution of the (conservative) term
$\nabla_{\mu}T^{\mu\nu}=0$, either in special relativity or in general
relativistic but stationary spacetimes (Del Zanna et al., 2007; Wu & Tang,
2015; Inghirami et al., 2016). It is interesting to note that the Z4
formulation of the Einstein equations was actually originally proposed in
first–order conservative form (Bona et al., 2004a; Bona & Palenzuela-Luque,
2005), though it does not seem that it has ever been implemented by using WENO
finite difference methods.

As we have just discussed, the very large size of the FO-Z4 system, as well as
its non-conservative aspect, restricts our choice of numerical methods that
can be applied to it. As realized in Dumbser et al. (2018, 2020) and Dumbser
et al. (2024), path-conservative Discontinuous Galerkin (DG) schemes may be
one route to higher order for the Einstein equations, but DG schemes tend to
be prohibitively expensive, both in terms of their memory usage as well as in
terms of computational complexity. This is because a DG scheme retains all the
higher order modes of a hyperbolic system and for a gigantic hyperbolic
system, like the one being considered, that exacerbates the memory
consumption. All these modes have also to be evolved in time with the result
that the computational complexity of evaluating so many evolutionary equations
for each of the modes of a DG scheme can again become prohibitive. Finite
volume Weighted Essentially Non-Oscillatory (WENO) schemes are another
alternative route to higher order and recently some efficient options have
been offered, see Balsara et al. (2023b), and the supplement of that paper, as
well as Dumbser et al. (2013). However, the finite volume WENO schemes require
the reconstruction of the same number of modes as a DG scheme of the same
order. Consequently, finite volume WENO schemes also have a prohibitive memory
usage, though their computational complexity and timestep restriction are much
better than DG schemes of the same order. The only remaining path to high
order accuracy at low computational cost is therefore through finite
difference WENO schemes (Shu & Osher, 1988, 1989; Jiang & Shu, 1996; Balsara &
Shu, 2000; Balsara et al., 2016, 2023a, 2024a, 2024b). Compared to their
finite volume and DG cousins, finite difference schemes have a memory usage
and computational complexity that can be described as downright Spartan.
Unlike DG schemes, they offer a robust CFL even at higher orders. As shown in
Balsara et al. (2024b), the ninth order finite difference WENO schemes are not
much more expensive compared to their third order variants, with the result
that the path to progressively higher order does not incur much additional
expense in terms of computational complexity. This owes to the dimension-by-
dimension approach that is used in finite difference WENO schemes.

It is worth noting that there are two variants of finite difference WENO
schemes. The well-known one, which we can refer to as the classical finite
difference WENO scheme (here referred to as FD–WENO) stems from Shu & Osher
(1989). For a long time, it was only available in a form that was suitable for
conservation laws. In that form, it is quite useless for first–order
non–conservative formulations of the Einstein equations, such as those
considered by Dumbser et al. (2018), Dumbser et al. (2024) or Brown et al.
(2012). Recently, the FD–WENO schemes have also been extended to hyperbolic
PDE systems with non-conservative products (Balsara et al., 2023a). In a paper
that slightly predates Shu & Osher (1989), Shu & Osher (1988) had also
presented an alternative formulation of finite difference WENO schemes, but
for a long time these methods had remained inaccessible to the broader
community. This is partly owing to the mathematically recondite nature of the
original presentation and partly owing to the fact that the original
formulation in Shu & Osher (1988) was far from broadly usable. We will refer
to this alternative strain of finite difference WENO schemes as the AFD–WENO
schemes. The original AFD–WENO schemes were, therefore, eclipsed by the vastly
more popular FD–WENO schemes. A paper by Merriman (2003) tried to make the
underlying mathematics of the AFD–WENO schemes more accessible. These days, of
course, the mathematics that underlies AFD–WENO schemes has become very
accessible; please see Section 2 from Balsara et al. (2024a) as well as its
Appendix A, which provides a computer algebra system-based derivation of the
scheme. The impetus for AFD–WENO schemes came when Cai & Ladeinde (2008) and
Nonomura et al. (2010) showed that FD–WENO schemes could not preserve a free
stream condition on geometrically complex meshes. In Jiang et al. (2013, 2014)
it was realized that FD–WENO suffered from an inability to preserve a free
stream condition because it was based on flux reconstruction. AFD–WENO schemes
were found to be free of such limitations. However, even in Zheng et al.
(2021), controlling the Gibbs phenomenon that arises from the higher order
flux derivatives in AFD–WENO was found to be an elusive enterprise. A full
resolution of controlling the Gibbs phenomenon only emerged in Balsara et al.
(2024a) where a different type of WENO interpolation was invented for that
purpose. The paper by Balsara et al. (2024a) still presented AFD–WENO schemes
for conservation laws. The extension of these AFD–WENO methods to hyperbolic
systems with non-conservative products was finally accomplished in Balsara et
al. (2024b). The upshot of this paragraph is that we now have two very
proficient methods for treating hyperbolic PDEs with a large number of non-
conservative products, such as arises in NR. We have the FD–WENO methods from
Balsara et al. (2023a) and we have the AFD–WENO methods from Balsara et al.
(2024b).

Modern astrophysical codes rely on capabilities that go beyond the baseline
scheme. A good example would be well-balancing, which makes it possible to
approach steady state on moderately resolved meshes, see Käppeli (2022) for a
review and Gaburro et al. (2021) for a first successful application of well-
balancing to numerical general relativity. Because FD–WENO schemes in
conservation form have been around for a while, well-balancing has been
developed for these methods (Xing & Shu, 2011). In Balsara et al. (2023a) we
showed that FD–WENO schemes for non-conservative products also tend to be
well-balanced if contact-discontinuity preserving Riemann solvers are used.
For AFD–WENO schemes in conservation form, Xu & Shu (2024) have shown that
they can be formulated so as to be well-balanced and can also preserve moving
equilibria provided one is willing to go through the cumbersome process of
identifying those equilibria. In recent work, we have seen the development of
well-balanced methods that preserve moving equilibria (Xu & Shu, 2024).
Another capability that goes beyond the baseline scheme is positivity
preservation; which is also referred to these days as the physical condition
preserving (PCP) property. For FD–WENO schemes in conservation form, Hu et al.
(2013) have developed positivity preserving formulations. This author and his
collaborators have also developed PCP formulations for AFD–WENO schemes
(Bhoriya et al., 2024). Astrophysical codes also have geometrical constraints
that require divergence-preservation (Balsara & Spicer, 1999; Balsara, 2010,
2012, 2014) and curl-preservation (Balsara et al., 2021; Balsara & Käppeli,
2022) and the author and his collaborators are extending these capabilities to
include AFD–WENO methods. We see, therefore, that there is an ecosystem of
ancillary algorithms that have been built, and are being built, around finite
difference WENO methods that make them very useful for NR.

The plan of the paper is the following. In Section 2 we recall the essential
features of the first order Z4 formulation of the Einstein field equations;
Section 3 contains the core novelties of this paper, by presenting a new class
of finite difference WENO schemes, whose validation is reported in Section 4.
Finally, Section 5 contains the conclusions of our work.

We assume a signature $\\{-,+,+,+\\}$ for the space-time metric and we use
Greek letters $\mu,\nu,\lambda,\ldots$ (running from 0 to 3) for four-
dimensional space-time tensor components, while Latin letters $i,j,k,\ldots$
(running from 1 to 3) are employed for three-dimensional spatial tensor
components. Moreover, we set $G=c=1$.

##  2 The damped first–order Z4 system of the Einstein equations

As usual for 3+1 formulations of the Einstein equations, the spacetime is
foliated through $\Sigma_{t}=const$ hypersurfaces as (Alcubierre, 2008;
Baumgarte & Shapiro, 2010; Gourgoulhon, 2012; Rezzolla & Zanotti, 2013)

$ds^{2}=g_{\mu\nu}dx^{\mu}dx^{\nu}=-(\alpha^{2}-\beta_{i}\beta^{i})dt^{2}+2\beta_{i}dx^{i}dt+\gamma_{ij}dx^{i}dx^{j}\,,$
(1)

where $g_{\mu\nu}$ is the spacetime metric tensor, $\alpha$ is the lapse,
$\beta^{i}$ is the shift and $\gamma_{ij}$ is the metric of the three
dimensional space. Within this framework, we adopt the damped version of the
Z4 formulation of the vacuum Einstein equations, originally proposed by
Gundlach et al. (2005) and recently reformulated with minor modifications by
Dumbser et al. (2024), i.e.

$G_{\mu\nu}+\nabla_{\mu}z_{\nu}+\nabla_{\nu}z_{\mu}-\nabla_{\pi}z^{\pi}g_{\mu\nu}-\kappa_{1}(n_{\mu}z_{\nu}+n_{\nu}z_{\mu})-\kappa_{2}n_{\pi}z^{\pi}g_{\mu\nu}=0\,,$
(2)

where $G_{\mu\nu}$ is the Einstein tensor and $n^{\mu}$ is the four–velocity
of the Eulerian observer. The four vector $z^{\mu}=\Theta n^{\mu}+Z^{\mu}$ is
artificially introduced to _clean_ the violations of the Einstein constraints,
while the two constant coefficients $\kappa_{1}$ and $\kappa_{2}$ can act as
damping mechanisms over the four vector $z^{\mu}$. Not present in the original
formulation by Bona et al. (2003, 2004b), such coefficients were introduced by
Gundlach et al. (2005) and slightly modified by Dumbser et al. (2024) in such
a way to avoid their mutual multiplication in the resulting PDE system. The
_extrinsic curvature_ of the hypersurface $\Sigma_{t}$, a crucial quantity of
differential geometry, is one of the primary variables of the Z4 formulation
and it is defined as

$K_{\mu\nu}=-\gamma^{\alpha}_{\mu}\nabla_{\alpha}n_{\nu}=-\nabla_{\mu}n_{\nu}-n_{\mu}a_{\nu}\,,$
(3)

where $a_{\mu}=n^{\nu}\nabla_{\nu}n_{\mu}$ is the acceleration of the Eulerian
observer. Upon the introduction of 30 auxiliary variables involving first
derivatives of the metric terms, namely

$\displaystyle
A_{i}:=\partial_{i}\ln\alpha=\frac{\partial_{i}\alpha}{\alpha}\,,\qquad
B_{k}^{\,\,i}:=\partial_{k}\beta^{i}\,,\qquad
D_{kij}:=\frac{1}{2}\partial_{k}\gamma_{ij}\,,\qquad$ (4)

and after fixing the gauge conditions in a pretty standard way as (Baumgarte &
Shapiro, 2010; Faber et al., 2007)

$\displaystyle\partial_{t}\ln\alpha-\beta^{k}{\partial_{k}\ln\alpha}$
$\displaystyle=$ $\displaystyle-g(\alpha)\alpha(K-K_{0}-2c\Theta)\,,$ (5)
$\displaystyle\partial_{t}\beta^{i}$ $\displaystyle=$
$\displaystyle\frac{3}{4}b^{i},$ (6) $\displaystyle\partial_{t}b^{i}$
$\displaystyle=$ $\displaystyle\partial_{t}\hat{\Gamma}^{i}-\eta b^{i}\,,$ (7)

the vacuum Einstein equations (2) lead to the following first order system of
hyperbolic PDEs:

$\displaystyle\partial_{t}\gamma_{ij}-\beta^{k}\partial_{k}\gamma_{ij}=\gamma_{ik}B_{j}^{\,\,k}+\gamma_{kj}B_{i}^{\,\,k}-2\alpha
K_{ij}\,,$ (8)
$\displaystyle\partial_{t}K_{ij}-\beta^{k}\partial_{k}K_{ij}+\alpha\partial_{(i}{A}_{j)}-\alpha{\gamma}^{kl}\left(\partial_{(k}{D}_{i)jl}-\partial_{(k}{D}_{l)ij}\right)+\alpha{\gamma}^{kl}\left(\partial_{(j}{D}_{i)kl}-\partial_{(j}{D}_{l)ik}\right)-2\alpha\partial_{(i}{Z}_{j)}=K_{ki}B_{j}^{\,\,k}+K_{kj}B_{i}^{\,\,k}$
$\displaystyle-\alpha
A_{i}A_{j}+\alpha\Gamma^{k}_{ij}A_{k}+\alpha\bigg{[}-2{\gamma}^{kn}{\gamma}^{pl}D_{knp}\left(D_{ijl}+D_{jil}-D_{lij}\right)+2{\gamma}^{kn}{\gamma}^{pl}D_{jnp}\left(D_{ikl}+D_{kil}-D_{lik}\right)$
$\displaystyle+\Gamma^{m}_{lm}\Gamma^{l}_{ij}-\Gamma^{m}_{lj}\Gamma^{l}_{im}\bigg{]}-2\alpha\Gamma^{k}_{ij}Z_{k}-\alpha\Theta\gamma_{ij}(\kappa_{1}+\kappa_{2})-2\alpha
K_{il}\gamma^{lm}K_{mj}+\alpha K_{ij}(K-2\Theta)\,,$ (9)
$\displaystyle\partial_{t}\Theta-\beta^{k}\partial_{k}\Theta-\frac{1}{2}\alpha{e^{2}}\left[\gamma^{ij}{\gamma}^{kl}\left(\partial_{(k}{D}_{i)jl}-\partial_{(k}{D}_{l)ij}\right)-\gamma^{ij}{\gamma}^{kl}\left(\partial_{(j}{D}_{i)kl}-\partial_{(j}{D}_{l)ik}\right)+2\gamma^{ij}\partial_{i}Z_{j}\right]=$
$\displaystyle=\frac{\alpha}{2}\,e^{2}\,\bigg{[}-2\gamma^{ij}{\gamma}^{kn}{\gamma}^{pl}D_{knp}\left(D_{ijl}+D_{jil}-D_{lij}\right)+2\gamma^{ij}{\gamma}^{kn}{\gamma}^{pl}D_{jnp}\left(D_{ikl}+D_{kil}-D_{lik}\right)$
$\displaystyle+\gamma^{ij}\left(\Gamma^{m}_{lm}\Gamma^{l}_{ij}-\Gamma^{m}_{lj}\Gamma^{l}_{im}\right)+K^{2}-K_{ij}\,K^{ij}\bigg{]}+\alpha\left[-\gamma^{ij}\,\Gamma^{k}_{ij}Z_{k}-Z^{k}A_{k}\right]-\alpha\Theta
K-\alpha\Theta(2\kappa_{1}+\kappa_{2})\,,$ (10)
$\displaystyle\partial_{t}Z_{i}-\beta^{k}\partial_{k}Z_{i}-\alpha\partial_{i}\Theta-\alpha\left[\gamma^{jm}\partial_{j}K_{mi}-\gamma^{mn}\partial_{i}K_{mn}\right]=Z_{k}\,B_{i}^{\,\,k}+\alpha\,\bigg{[}-\gamma^{jm}(\Gamma^{n}_{jm}K_{ni}+\Gamma^{n}_{ji}K_{mn})$
$\displaystyle+\gamma^{mn}(\Gamma^{l}_{im}K_{ln}+\Gamma^{l}_{in}K_{ml})\bigg{]}+\alpha[-2\,{K_{i}}^{j}\,Z_{j}-\Theta\,A_{i}-\kappa_{1}Z_{i}]\,,$
(11) $\displaystyle\partial_{t}A_{i}-{\beta^{k}\partial_{k}A_{i}}+\alpha
g(\alpha)\left(\gamma^{mn}\partial_{i}K_{mn}-\partial_{i}K_{0}-2c\partial_{i}\Theta\right)=-\alpha
A_{i}\left(K-K_{0}-2\Theta c\right)\left(g(\alpha)+\alpha
g^{\prime}(\alpha)\right)$ (12) $\displaystyle+2\alpha
g(\alpha)K^{jk}D_{ijk}+B_{i}^{\,\,k}~{}A_{k}\,,$ (13)
$\displaystyle\partial_{t}B_{k}^{\,\,i}-s\left(\frac{3}{4}\partial_{k}b^{i}-\alpha^{2}\mu\,\gamma^{ij}\gamma^{nl}\left(\partial_{k}D_{ljn}-\partial_{l}D_{kjn}\right)\right)=0\,,$
(14)
$\displaystyle\partial_{t}D_{kij}-{\beta^{l}\partial_{l}D_{kij}}-\frac{1}{2}{\gamma}_{mi}\partial_{(k}{B}_{j)}^{m}-\frac{1}{2}{\gamma}_{mj}\partial_{(k}{B}_{i)}^{m}+\alpha\partial_{k}{K}_{ij}=B_{k}^{\,\,m}D_{mij}+B_{j}^{\,\,m}D_{kmi}+B_{i}^{\,\,m}D_{kmj}-\alpha
A_{k}{K}_{ij}\,.$ (15)

The equations (8)–(15), augmented by the gauge conditions (5)–(7), form a non-
conservative system which can be written as

$\frac{\partial\mathbf{U}}{\partial
t}+\mathbf{A}_{i}(\mathbf{U})\frac{\partial\mathbf{U}}{\partial
x_{i}}={\bf{S}}(\mathbf{U}),\qquad\textnormal{or,
equivalently,}\qquad\frac{\partial\mathbf{U}}{\partial
t}+\mathbf{A}(\mathbf{U})\cdot\nabla\mathbf{U}={\bf{S}}(\mathbf{U}),$ (16)

where $\mathbf{U}$ is the state vector, composed of 54 dynamical variables,
i.e. 10 for the lapse, the shift vector and the metric components, 6 for
$K_{ij}$, 4 for the $z^{\mu}$ four vector, 3 for $A_{i}$, 9 for
$B_{i}^{\,\,j}$, 18 for $D_{ijk}$, 1 for $K_{0}$ and 3 for $b^{i}$. The source
term $\bf S(\mathbf{U})$ contains algebraic terms only. The hyperbolic nature
of (16) has been proved by Dumbser et al. (2024) after a careful analysis of
the matrix $\mathbf{A}(\mathbf{U})$ via mathematical software packages, and
the interested reader is pointed to that work, and especially the Appendix in
that work, for further details. We just recall here that, as already noticed
by Dumbser et al. (2018), hyperbolicity is helped by moving any spatial
derivatives of the metric terms $\alpha$, $\beta^{i}$ and ${\gamma}_{ij}$ to
the right hand side via the auxiliary variables (4).

A few additional comments regarding the terms entering Eq. (8)–(15) are worth
giving

  * •

$\Gamma^{i}_{jk}={\gamma}^{kl}\left(D_{ijl}+D_{jil}-D_{lij}\right)$ are the
Christoffel symbols of the spatial metric $\gamma_{ij}$.

  * •

$K$ in Eq. (5) (and elsewhere) is the trace of the extrinsic curvature, but it
is not a primary variables, i.e. it does not belong to the vector
$\mathbf{U}$.

  * •

The function $g(\alpha)$ in Eq. (5) is set to $g(\alpha)=2/\alpha$ for the
_1+log gauge condition_ , and to $g(\alpha)=1$ for the _harmonic gauge
condition_.

  * •

The factor $c$ in Eq. (5) (and elsewhere) is always zero, except for the test
of the gauge wave discussed in Sect. 4.1.

  * •

The factor $s$ in Eq. (14), either 1 or 0, is used to switch the
_gamma–driver_ on or off.

  * •

The quantities $\hat{\Gamma}^{i}$ in Eq. (7) for the gamma driver are defined
as $\hat{\Gamma}^{i}=\gamma^{jk}\,\Gamma^{i}_{jk}+2\gamma^{ij}Z_{j}$, but they
are not primary variables. Hence their evolution is obtained from that of the
primary variables, i.e.

$\displaystyle\partial_{t}\Gamma^{i}_{jk}$ $\displaystyle=$
$\displaystyle\gamma^{im}\beta^{r}\left[\partial_{r}D_{jmk}+\partial_{r}D_{kjm}-\partial_{r}D_{mjk}\right]+\partial_{(j}B_{k)}^{\,\,i}-\alpha\gamma^{im}\left(\partial_{j}K_{mk}+\partial_{k}K_{jm}-\partial_{m}K_{jk}\right)+$
(17)
$\displaystyle+\gamma^{im}\big{[}D_{jmn}B_{k}^{\,\,n}+D_{nmk}B_{j}^{\,\,n}+D_{knm}B_{j}^{\,\,n}+D_{njm}B_{k}^{\,\,n}-D_{mjn}B_{k}^{\,\,n}-D_{mnk}B_{j}^{\,\,n}\big{]}$
$\displaystyle-\alpha\gamma^{im}\left(A_{j}K_{mk}+A_{k}K_{jm}-A_{m}K_{jk}\right)$
$\displaystyle+\big{[}-2\gamma^{ip}\gamma^{mq}\beta^{r}D_{rpq}-\gamma^{mr}B_{r}^{\,\,i}+2\alpha\gamma^{ip}\gamma^{mq}K_{pq}\big{]}\left(D_{jmk}+D_{kjm}-D_{mjk}\right).$

  * •

The factor $\eta$ in Eq. (7) is a damping parameter for the gamma–driver.

  * •

The factor $e$ in Eq. (2) is the cleaning speed of the Einstein energy
constraint.

We are now in a position to address the discretization of the system (16) via
FD–WENO schemes, to which Section 3 is entirely devoted.

##  3 A consideration of several finite difference WENO schemes as they
pertain to numerical relativity

We had briefly mentioned that finite difference WENO schemes offer the three-
fold advantages of very high accuracy, very low memory usage and low
computational complexity. To use them well, and also to make decisions about
their strengths and weaknesses, it is very important to understand the finite
difference WENO philosophy very briefly and why it offers these three-fold
benefits. Consider, therefore, the finite volume schemes that are quite common
in computational astrophysics. They all start with a volume-averaged
representation of the solution vector in each zone. By polling the neighboring
zones, a multidimensional high order finite volume scheme will use WENO
reconstruction to build all the higher order moments in multiple dimensions.
As the order of accuracy increases, the number of moments with cross-terms
that have to be built also increases. By contrast, consider a finite
difference approach to a simple one-dimensional conservation law, which we
write as:

${{\partial}_{t}}\mathbf{U}+{{\partial}_{x}}\mathbf{F}\left(\mathbf{U}\right)=0\
\text{ }\Leftrightarrow\text{
}{{\partial}_{t}}\mathbf{U}=-{{\partial}_{x}}\mathbf{F}\left(\mathbf{U}\right)\,.$
(18)

A finite difference scheme will start with a mesh function,
$\left\\{{{\mathbf{U}}_{i}}\right\\}$ that is provided in the form of point
values of the solution vector at each of the zone centers “$i$”. Here we
assume a uniform mesh with zones of size $\Delta x$, taking timesteps of size
$\Delta t$. Because we are starting with point values, the finite difference
scheme requires us to accurately evaluate the gradient
${{\partial}_{x}}\mathbf{F}\left(\mathbf{U}\right)$ at that same zone-centered
point. The transcription to multiple dimensions is, therefore, very easy
because we simply want all the flux gradients in all three directions to be
evaluated at the same point with sufficiently high accuracy. This shows us
that finite difference schemes operate on a dimension-by-dimension basis. This
is the reason why it was acceptable to only show eqn. (18) in one dimension;
because additional dimensions only contribute additively. Unlike high order
finite volume approaches where a 3D finite volume scheme can be substantially
more than three times costlier than a one-dimensional scheme, a 3D finite
difference scheme will only be three times costlier than a one-dimensional
one. We also see that the memory usage is very favorable for finite difference
schemes. This is because we do not need to reconstruct all the higher order
moments in all dimensions and store them in computer memory, as we would do
for a finite volume scheme. A finite difference scheme only needs to retain
the point values of the solution vector at the zone centers. Any higher
moments that are needed will only be needed in a dimension-by-dimension
fashion and can be discarded from computer memory once that dimension has been
processed. By the same token, it is also possible to show that the
computational complexity of a finite difference WENO scheme does not increase
by much as one progresses to higher orders of spatial accuracy; please see
Table V in subsection 5.4 of Balsara et al. (2024b).

Now that the advantages of finite difference WENO have been documented, we
provide very brief descriptions of these schemes here so that the reader has
sufficient background, concentrated all in one place, with which to understand
these schemes as they will be used for the FO-Z4 system from Sect. 2.

###  3.1 FD-WENO Schemes for Conservation Laws

Although our immediate goal is not to understand conservation laws, it helps
to quickly document the WENO philosophy that led to FD-WENO schemes for
conservation laws. (Without an understanding of conservation laws it is
impossible to understand systems in non-conservative form.) The primary goal
of FD-WENO for conservation laws is to try and write eqn. (18) for any zone
“$i$” on a one-dimensional mesh as

${{\partial}_{t}}\mathbf{U}_{i}=-\frac{{{{\mathbf{\hat{F}}}}_{i+1/2}}-{{{\mathbf{\hat{F}}}}_{i-1/2}}}{\Delta
x}.$ (19)

Here ${\mathbf{\hat{F}}}_{i+1/2}$ and ${{\mathbf{\hat{F}}}_{i-1/2}}$ are the
reconstructed numerical fluxes at the zone boundaries ${{x}_{i}}+{\Delta
x}/{2}\;$ and ${{x}_{i}}-{\Delta x}/{2}\;$. For this to work out as a viable
high order scheme, we need the finite difference approximation (FDA) given by
${\left({{{\mathbf{\hat{F}}}}_{i+1/2}}-{{{\mathbf{\hat{F}}}}_{i-1/2}}\right)}/{\Delta
x}\;$ to approximate ${{\left({{\partial}_{x}}\mathbf{F}\right)}_{i}}$ with
very high order of accuracy. Therefore, the problem can be stated as follows.
We have to start with the so–called _physical fluxes_ evaluated from the point
values on the mesh, i.e.
$\left\\{\mathbf{F}\left({{\mathbf{U}}_{i}}\right)\right\\}$, and obtain
_numerical fluxes_ such that
${\left({{{\mathbf{\hat{F}}}}_{i+1/2}}-{{{\mathbf{\hat{F}}}}_{i-1/2}}\right)}/{\Delta
x}\;={{\left({{\partial}_{x}}\mathbf{F}\right)}_{i}}+O\left(\Delta{{x}^{k}}\right)$
for a spatially $k^{th}$ order accurate scheme.

The denouement of the previously-posed problem comes from the fundamental
theorem of integral calculus, as was first realized by Shu & Osher (1989).
Consider a function $\mathbf{f}\left(x\right)$ that is defined in terms of
another functional $\mathcal{F}\left(\xi\right)$ as follows:-

$\mathbf{f}\left(x\right)\equiv\frac{1}{\Delta x}\int\limits_{x-\Delta
x/2}^{x+\Delta x/2}{\mathcal{F}\left(\xi\right)d\xi}.$ (20)

Then the fundamental theorem of integral calculus tells us that:-

$\left({{\partial}_{x}}\mathbf{f}\right){{|}_{x={{x}_{i}}}}=\frac{1}{\Delta
x}\left[\mathcal{F}\left({{x}_{i}}+\frac{\Delta
x}{2}\right)-\mathcal{F}\left({{x}_{i}}-\frac{\Delta x}{2}\right)\right].$
(21)

We see, therefore, that the left hand side of eqn. (21) is exactly of the same
form as the right hand side of the second equation in eqn. (18) whereas the
right hand side of eqn. (21) is exactly of the same form as the right hand
side of eqn. (19). We have, therefore, found a way of connecting the gradient
in eqn. (18) to the FDA in eqn. (19). Consequently, we seek a function
$\mathcal{F}\left(\xi\right)$ whose value
$\mathcal{F}\left({{x}_{i+1/2}}\right)$ at the zone boundary approximates
${{\mathbf{\hat{F}}}_{i+1/2}}$ with a sufficiently high order of accuracy.
Now, the very familiar problem of “reconstruction by primitive” that was
introduced in Colella & Woodward (1984) (and extended to ENO schemes with
accreditation in Harten et al. (1986)) comes to the rescue. It says that the
reconstruction polynomial for the zone boundary “$i+1/2$” that we seek should
be such that it matches the condition

$\int\limits_{\xi={{x}_{i-j}}-\Delta x/2}^{\xi={{x}_{i-j}}+\Delta
x/2}{\mathcal{F}\left(\xi\right)d\xi}=\Delta x\
\mathbf{F}\left({{\mathbf{U}}_{i-j}}\right)$ (22)

for some zones “$j$” that are adjacent to the zone boundary “$i+1/2$” under
consideration. If the polynomial is chosen to be of sufficient order of
accuracy, the task that was identified in the previous paragraph is
accomplished. We can make the polynomial sufficiently accurate by making the
stencil operations alluded to in eqn. (22) wide enough. We choose
$\mathcal{F}\left(\xi\right)$ to be a consistent polynomial that is of degree
“$k-1$” in order to obtain a $k^{th}$ order accurate scheme. This polynomial
will be pinned down more precisely in subsequent paragraphs.

While the process of identifying the polynomial associated with reconstructing
the fluxes was described in the previous paragraph, more is needed in order to
obtain a successful scheme. Such a scheme needs to be upwinded and the
reconstruction process that was briefly described in eqn. (22) needs to be
based on non-linear hybridization. To appreciate the concept of upwinding,
realize that any flux “$\mathbf{F}\left(\mathbf{U}\right)$” can be split as

$\mathbf{F}\left(\mathbf{U}\right)=\frac{1}{2}\left(\mathbf{F}\left(\mathbf{U}\right)+S\text{
}\mathbf{U}\right)+\frac{1}{2}\left(\mathbf{F}\left(\mathbf{U}\right)-S\text{
}\mathbf{U}\right)={{\mathbf{F}}^{\mathbf{+}}}\left(\mathbf{U}\right)+{{\mathbf{F}}^{-}}\left(\mathbf{U}\right)$
(23)

with the definitions

${{\mathbf{F}}^{\mathbf{+}}}\left(\mathbf{U}\right)\equiv\frac{1}{2}\left(\mathbf{F}\left(\mathbf{U}\right)+S\text{
}\mathbf{U}\right)\quad\text{ and
}\quad{{\mathbf{F}}^{-}}\left(\mathbf{U}\right)\equiv\frac{1}{2}\left(\mathbf{F}\left(\mathbf{U}\right)-S\text{
}\mathbf{U}\right)$ (24)

Here the wave speed “$S$” is the maximum of the absolute values of all the
signal speeds from the flux. This choice of “$S$” ensures that the Jacobian
${\partial{{\mathbf{F}}^{\mathbf{+}}}\left(\mathbf{U}\right)}/{\partial\mathbf{U}}\;$
will always have non-negative eigenvalues and the Jacobian
${\partial{{\mathbf{F}}^{-}}\left(\mathbf{U}\right)}/{\partial\mathbf{U}}\;$
will always have non-positive eigenvalues. (As a shorthand, we will often use
${{\mathbf{F}}^{\mathbf{+}}}\left(\mathbf{U}\right)\to{{\mathbf{F}}^{\mathbf{+}}}$
and ${{\mathbf{F}}^{-}}\left(\mathbf{U}\right)\to{{\mathbf{F}}^{-}}$.) The
flux splitting in the above two equations is known as the LLF (locally
Lax–Friedrichs) flux splitting. Now please focus on Fig. 1. We see a mesh with
zones, and zone boundaries, as well as a mesh function that is specified at
zone centers. At each zone center, we can use the point value of the mesh
function to also evaluate the zone-centered point value for the flux. These
zone-centered flux values are to be used to obtain the numerical flux
${{\mathbf{\hat{F}}}_{i+1/2}}$ at the purple zone boundary in Fig. 1. But
realize that the numerical flux ${{\mathbf{\hat{F}}}_{i+1/2}}$ is constituted
by the sum of $\mathbf{\hat{F}}_{i+1/2}^{+}$ (which carries the right-going
flux contributions) and $\mathbf{\hat{F}}_{i+1/2}^{-}$ (which carries the
left-going flux contributions). If $\mathbf{\hat{F}}_{i+1/2}^{+}$ and
$\mathbf{\hat{F}}_{i+1/2}^{-}$ are obtained via a high order, one-dimensional
finite-volume-style reconstruction then the entire scheme shown in eqn. (19)
will be high order accurate in space. This paragraph has, therefore, shown how
the flux splitting is to be used in FD-WENO schemes for conservation laws.
Notice that in light of the discussion surrounding eqns. (20), (21) and (22),
the reconstruction has to be a finite-volume style reconstruction. In such a
reconstruction, the zone averages are known, as shown in the right hand side
of eqn. (22), and the reconstruction polynomial has to be constructed so as to
match the zone averages.

Figure 1: The schematic structure of a fifth order FD-WENO scheme in
conservation form. Please focus on the purple zone boundary at $i+1/2$. The
numerical flux at that zone boundary is made up of two split flux
contributions $\mathbf{F}^{+}$ and $\mathbf{F}^{-}$. The stencils contributing
to the reconstruction of the right-going flux $\mathbf{F}^{+}$ at the purple
zone boundary are shown. The stencils contributing to the reconstruction of
the left-going flux $\mathbf{F}^{-}$ at the purple zone boundary are also
shown. .

The previous paragraph addressed the issue of flux splitting. But it did not
address the twin issues of upwinding and non-linear hybridization. We know
that these two concepts are essential for obtaining a stable scheme. Upwinding
requires that any stencils that contribute to the right-going flux
$\mathbf{\hat{F}}_{i+1/2}^{+}$ must include the zone that is left of the zone
boundary “$i+1/2$”. As a result, those stencils must include zone “$i$” in
Fig. 1. We see that this requirement is met by the stencils shown in Fig. 1.
Upwinding also requires that any stencils that contribute to the left-going
flux $\mathbf{\hat{F}}_{i+1/2}^{-}$ must include the zone that is right of the
zone boundary “$i+1/2$”. As a result, those stencils must include zone “$i+1$”
in Fig. 1; and we see that this requirement is met by the stencils shown in
Fig. 1.

Now that flux splitting and upwinding have been clarified, let us address the
issue of non-linear hybridization. Over the years, there have been many
interpretations of what it means to non-linearly hybridize the available
stencils in a WENO scheme (Jiang & Shu, 1996; Balsara & Shu, 2000; Levy et
al., 2000; Henrick et al., 2006; Borges et al., 2008; Zhu & Qiu, 2016; Balsara
et al., 2016; Cravero & Semplice, 2016; Zhu & Shu, 2018). Fig. 1 shows the
WENO–AO approach of Balsara et al. (2016); where “AO” stands for adaptive
order. Fig. 1 focuses on a fifth order WENO-AO reconstruction, but the methods
have been extended up to $11^{th}$ order. The three smaller left-biased, zone-
centered and right-biased stencils cover the zones of interest; i.e., zone
“$i$” and zone “$i+1$” respectively. Their non-linear hybridization would give
a stable stencil that is third order accurate. For example, consider the
stencils with the right arrow in Fig. 1 that contribute to
$\mathbf{\hat{F}}_{i+1/2}^{+}$. The smaller stencils include the left-biased
stencil which includes zones $\left\\{i-2,i-1,i\right\\}$, the centered
stencil which includes zones $\left\\{i-1,i,i+1\right\\}$ and the right-biased
stencil which includes zones $\left\\{i,i+1,i+2\right\\}$. Together, these
three smaller stencils, if they are non-linearly hybridized amongst
themselves, would yield a stable, spatially third order accurate
reconstruction over the zone “$i$” which must be included in the upwinding of
$\mathbf{\hat{F}}_{i+1/2}^{+}$. To get a fifth order accurate reconstruction
over the zone “$i$”, one must make a non-linear hybridization of these smaller
stencils with the larger stencil in Fig. 1 which includes zones
$\left\\{i-2,i-1,i,i+1,i+2\right\\}$. The smaller stencils guarantee
stability, when stability becomes an issue. The larger stencil provides higher
(fifth) order accuracy, when the solution is smooth enough to justify the
higher accuracy. The stencils for obtaining a fifth order in space
approximation for the right-going flux are shown by the stencils with a right-
pointing arrow in Fig. 1. Fig. 1 also displays the smaller and larger stencils
used for the spatially fifth order accurate, upwinded reconstruction of the
left-going flux $\mathbf{\hat{F}}_{i+1/2}^{-}$. The stencils for obtaining a
fifth order in space approximation for the left-going flux are shown by the
stencils with a left-pointing arrow in Fig. 1. We can see that all those
stencils in Fig. 1 cover the zone “$i+1$”. The final FD-WENO scheme for
conservation laws is then given by:-

${{\partial}_{t}}\mathbf{U}_{i}=-\frac{\left(\mathbf{\hat{F}}_{i+1/2}^{+}+\mathbf{\hat{F}}_{i+1/2}^{-}\right)-\left(\mathbf{\hat{F}}_{i-1/2}^{+}+\mathbf{\hat{F}}_{i-1/2}^{-}\right)}{\Delta
x}$ (25)

The reader can obtain further detail associated with WENO-AO reconstruction
shown in Fig. 1 from Sections 2 and 3 of Balsara et al. (2016). We should also
mention that the reconstruction should be done in characteristic variables in
order to obtain the best quality solution. This completes our very brief
description of FD-WENO for conservation laws. This FD-WENO scheme for
conservation laws will prove indispensable in the next section for
understanding FD-WENO schemes for systems in non-conservation form like the
FO-Z4 formulation of the Einstein equations.

###  3.2 FD-WENO Schemes for Hyperbolic Systems with Non-Conservative Products

For a very long time, FD-WENO schemes were only available for conservation
laws. However, first order formulations of the Einstein equations, as are
nowadays presented (Brown et al., 2012; Dumbser et al., 2018, 2024), are cast
in a non conservative form like (16). Here, to simplify the description of the
numerical scheme, we disregard the source terms on the right hand side and we
consider

$\displaystyle{{\partial}_{t}}\mathbf{U}+\mathbf{A}\left(\mathbf{U}\right){{\partial}_{x}}\mathbf{U}=0\
\text{ }\Leftrightarrow\text{
}{{\partial}_{t}}\mathbf{U}=-\mathbf{A}\left(\mathbf{U}\right){{\partial}_{x}}\mathbf{U}\,.$
(26)

Such situations are best handled in fluctuation form. The problem was that,
until the advent of Balsara et al. (2023a), fluctuation form-based FD-WENO
schemes were not available for eqn. (26). Eqn. (16) in Section 4 of Balsara et
al. (2023a) yields a scheme that can be applied to eqn. (26). In Section 5 of
Balsara et al. (2023a) we do provide a derivation of the scheme. However, that
derivation is somewhat harder to understand, so we present a much simpler
derivation here. We give this derivation in the limit of a linear flux,
$\mathbf{F}=\mathbf{AU}$, with “$\mathbf{A}$” is a constant matrix, but the
final form that we will obtain will be in fluctuation form so that it can be
extended to include any Riemann solver that accommodates the non-linearities
at the zone boundaries. As in Balsara et al. (2023a), we work in the limit of
the LLF flux.

Figure 2: Panel shows part of the mesh around zone “i”. The mesh functions are
collocated at the zone centers, as shown by the thick dots. The zone
boundaries are shown by the vertical lines. The figure also shows the stencils
associated with the zone “$i$” for the fifth order WENO–AO
reconstruction/interpolation. We have three smaller third order stencils and a
large fifth order stencil. The reconstructed/interpolated variables at the
zone boundaries are shown with a caret. The variables with a superscript star
are resolved states obtained by the pointwise application of a simple HLL or
LLF Riemann solver at the zone boundaries.

Consider Fig. 2 where the solution vector has itself been reconstructed.
Because we are only reconstructing the solution vector, only one
reconstruction is needed within each zone. (This is different from what is
shown in Fig. 1, which uses two reconstructions for each of the right-going
and left-going fluxes per zone.) As a result, at each zone boundary, say the
zone boundary “$i+1/2$” in Fig. 2, we have the left state
$\mathbf{\hat{U}}_{i+1/2}^{-}$ and the right state
$\mathbf{\hat{U}}_{i+1/2}^{+}$. Please try to understand how
$\mathbf{\hat{F}}_{i+1/2}^{-}$ in Fig. 1 relates to
$\mathbf{\hat{U}}_{i+1/2}^{+}$ in Fig. 2. Likewise, please try to understand
how $\mathbf{\hat{F}}_{i+1/2}^{+}$ in Fig. 1 relates to
$\mathbf{\hat{U}}_{i+1/2}^{-}$ in Fig. 2, With “$\mathbf{A}$” held constant,
we first write

$\displaystyle\mathbf{\hat{F}}_{i+1/2}^{-}$
$\displaystyle=\frac{1}{2}\left(\mathbf{A}-S\mathbf{I}\right)\mathbf{\hat{U}}_{i+1/2}^{+}\quad;\quad\mathbf{\hat{F}}_{i+1/2}^{+}=\frac{1}{2}\left(\mathbf{A}+S\mathbf{I}\right)\mathbf{\hat{U}}_{i+1/2}^{-}\quad;\quad$
(27) $\displaystyle\mathbf{\hat{F}}_{i-1/2}^{-}$
$\displaystyle=\frac{1}{2}\left(\mathbf{A}-S\mathbf{I}\right)\mathbf{\hat{U}}_{i-1/2}^{+}\quad;\quad\mathbf{\hat{F}}_{i-1/2}^{+}=\frac{1}{2}\left(\mathbf{A}+S\mathbf{I}\right)\mathbf{\hat{U}}_{i-1/2}^{-}$

It is then easy to show that:-

$\displaystyle-\left(\mathbf{\hat{F}}_{i+1/2}^{-}+\mathbf{\hat{F}}_{i+1/2}^{+}\right)$
$\displaystyle=-\frac{1}{2}\left(\mathbf{A}-S\mathbf{I}\right)\mathbf{\hat{U}}_{i+1/2}^{+}-\frac{1}{2}\left(\mathbf{A}+S\mathbf{I}\right)\mathbf{\hat{U}}_{i+1/2}^{-}$
(28)
$\displaystyle=-\frac{1}{2}\left(\mathbf{A}-S\mathbf{I}\right)\mathbf{\hat{U}}_{i+1/2}^{+}+\frac{1}{2}\left(\mathbf{A}-S\mathbf{I}\right)\mathbf{\hat{U}}_{i+1/2}^{-}-\mathbf{A\hat{U}}_{i+1/2}^{-}$
$\displaystyle=-\frac{1}{2}\left(\mathbf{A}-S\mathbf{I}\right)\left(\mathbf{\hat{U}}_{i+1/2}^{+}-\mathbf{\hat{U}}_{i+1/2}^{-}\right)-\mathbf{A\hat{U}}_{i+1/2}^{-}$
$\displaystyle=-\mathbf{D}_{i+1/2}^{-}-\mathbf{A\hat{U}}_{i+1/2}^{-}$

It is also easy to show via a similar derivation that:-

$\displaystyle\left(\mathbf{\hat{F}}_{i-1/2}^{-}+\mathbf{\hat{F}}_{i-1/2}^{+}\right)$
$\displaystyle=\frac{1}{2}\left(\mathbf{A}-S\mathbf{I}\right)\mathbf{\hat{U}}_{i-1/2}^{+}+\frac{1}{2}\left(\mathbf{A}+S\mathbf{I}\right)\mathbf{\hat{U}}_{i-1/2}^{-}$
(29)
$\displaystyle=-\frac{1}{2}\left(\mathbf{A}+S\mathbf{I}\right)\mathbf{\hat{U}}_{i-1/2}^{+}+\frac{1}{2}\left(\mathbf{A}+S\mathbf{I}\right)\mathbf{\hat{U}}_{i-1/2}^{-}+\mathbf{A\hat{U}}_{i-1/2}^{+}$
$\displaystyle=-\frac{1}{2}\left(\mathbf{A}+S\mathbf{I}\right)\left(\mathbf{\hat{U}}_{i-1/2}^{+}-\mathbf{\hat{U}}_{i-1/2}^{-}\right)+\mathbf{A\hat{U}}_{i-1/2}^{+}$
$\displaystyle=-\mathbf{D}_{i-1/2}^{+}+\mathbf{A\hat{U}}_{i-1/2}^{+}$

In the above two equations, and for a linear system, we can define the left-
going and right-going fluctuations as $\mathbf{D}_{i+1/2}^{-}$ $\equiv$
$\left[\left(\mathbf{A}-S\mathbf{I}\right)/2\right]$
$\left(\mathbf{\hat{U}}_{i+1/2}^{+}-\mathbf{\hat{U}}_{i+1/2}^{-}\right)$ and
$\mathbf{D}_{i-1/2}^{+}\equiv\left[\left(\mathbf{A}+S\mathbf{I}\right)/2\right]\left(\mathbf{\hat{U}}_{i-1/2}^{+}-\mathbf{\hat{U}}_{i-1/2}^{-}\right)$
respectively. We will soon point to references where their non-linear
extensions can be obtained. Putting the above two equations together, we can
now show that

$-\frac{\left(\mathbf{\hat{F}}_{i+1/2}^{+}+\mathbf{\hat{F}}_{i+1/2}^{-}\right)-\left(\mathbf{\hat{F}}_{i-1/2}^{+}+\mathbf{\hat{F}}_{i-1/2}^{-}\right)}{\Delta
x}=-\frac{1}{\Delta
x}\left(\mathbf{D}_{i+1/2}^{-}+\mathbf{D}_{i-1/2}^{+}\right)-\mathbf{A}\frac{\left(\mathbf{\hat{U}}_{i+1/2}^{-}-\mathbf{\hat{U}}_{i-1/2}^{+}\right)}{\Delta
x}.$ (30)

While the above derivation was done in a simpler context, it is now easy to
identify the full update equation:-

${{\partial}_{t}}{{\mathbf{U}}_{i}}=-\frac{1}{\Delta
x}\left(\mathbf{D}_{HLLI}^{-}\left(\mathbf{\hat{U}}_{i+1/2}^{-},\mathbf{\hat{U}}_{i+1/2}^{+}\right)+\mathbf{D}_{HLLI}^{+}\left(\mathbf{\hat{U}}_{i-1/2}^{-},\mathbf{\hat{U}}_{i-1/2}^{+}\right)\right)-\mathbf{A}\left({{\mathbf{U}}_{i}}\right)\frac{\left(\mathbf{\hat{U}}_{i+1/2}^{-}-\mathbf{\hat{U}}_{i-1/2}^{+}\right)}{\Delta
x}$ (31)

The fluctuation
$\mathbf{D}_{HLLI}^{-}\left(\mathbf{\hat{U}}_{i+1/2}^{-},\mathbf{\hat{U}}_{i+1/2}^{+}\right)$
captures the contribution from the left-going waves at zone boundary
“$i+1/2$”. The fluctuation
$\mathbf{D}_{HLLI}^{+}\left(\mathbf{\hat{U}}_{i-1/2}^{-},\mathbf{\hat{U}}_{i-1/2}^{+}\right)$
captures the contribution from the right-going waves at zone boundary
“$i-1/2$”. The subscript “HLLI” indicates that we are using the HLLI Riemann
solver from Dumbser & Balsara (2016), which can indeed handle non-linearities
in PDE systems of the form shown in eqn. (26). Appendix C of Dumbser & Balsara
(2016) provides further details on evaluating fluctuations for PDEs with non-
linearities.

In the previous sub-section we have mentioned that it is advantageous to carry
out the reconstruction in the characteristic variables. In a _tour de force_ ,
the 54 eigenvectors of FO-Z4 system for the full Einstein–Euler system have
been derived in Dumbser et al. (2024). However, it would be prohibitive to
project the entire system into that eigenspace for carrying out a
reconstruction of the solution vector. Fortunately, the eigenspace splits into
a 5 component subspace for the hydrodynamic system and an eigenspace for the
rest of the system. Since shocks form only in the hydrodynamic sector, while
in the spacetime part all the corresponding fields are linearly degenerate and
no shocks can form, we treat the reconstruction of the variables in the GR
sector without any characteristic projection. This is not a fundamental
limitation, it is only a time-saving expedient in light of the fact that the
first order hyperbolic formulation of GR results in such a large system.

It is also useful to observe that eqn. (31) is not in exact flux conservative
form. An exact flux conservative form is essential (because of the Lax-
Wendroff theorem) for accurately predicting shock locations when the system is
conservative, as is the case for hydrodynamics. The modern trend is to seek
out schemes that revert to a conservation form when such is present. However,
eqn. (31) does derive from eqn. (25) which is in conservation form, so for
most practical cases it does a rather good job in predicting shock locations.
Besides, since the solutions of the Einstein equations are smooth, the lack of
a conservation form is not a major impediment. As a result, eqn. (31) has an
important place in the solution of the Einstein equations. This completes our
description of FD-WENO schemes as they are extended to hyperbolic PDEs, like
the FO-Z4 system, that are not in conservation form.

###  3.3 AFD-WENO Schemes for Conservation Laws

The previous Sub-section has shown us how to obtain a scheme that is fully in
non-conservative form. It is very easy to implement and very suitable for the
Einstein equations in vacuum, i.e. when evolving the hydrodynamics is not a
priority. However, it is often the case that one has to simultaneously evolve
the equations for relativistic hydrodynamics in addition to the Einstein
equations. In that case, we will need a scheme that can handle non-
conservative products and is, nevertheless, versatile enough to revert to a
conservation form when such a conservation form is present in the problem. The
scheme in the previous Sub-section does not meet this requirement. We have
also seen that a study of conservation laws is a first step towards deriving
methods that can handle non-conservative products. Therefore, in this Sub-
section we will present an alternative finite difference WENO (AFD-WENO)
scheme for conservation laws. This will be an essential first step towards
deriving schemes that can seamlessly accommodate both a conservation law as
well as non-conservative products – that will be the task of the next section.

This and the next Sub-section rely on interpolation, which is different from
the reconstruction that was used in the previous two Sub-sections. It is,
therefore, worthwhile to make a distinction between reconstruction and
interpolation. Reconstruction is used in all finite volume astrophysical codes
and also in the previous FD-WENO schemes where the fluxes or states are
reconstructed. It consists of starting with the zone averages in a given
stencil and obtaining therefrom the high degree polynomial whose integration
over each of the zones of the stencil matches the original zone averages.
Interpolation is used less often in the numerical solution of conservation
laws, however it is the approach that will be needed in this and the next Sub-
section. It consists of starting with the point values at each of the zone
centers of a stencil and obtaining therefrom the high degree polynomial that
matches those point values. Therefore, the two words, reconstruction and
interpolation, carry different connotations. When applied to the same stencil,
reconstruction and interpolation produce polynomials with the same degree.
However, the underlying polynomial coefficients that are produced by invoking
reconstruction or interpolation on a given stencil can indeed be very
different. Standard WENO concepts like linear weights, smoothness indicators,
normalized non-linear weights etc. are often the same for reconstruction and
interpolation; so there is indeed a beneficial transference of knowledge
between them. WENO reconstruction, especially as it relates to FD-WENO
schemes, has been amply documented in the literature starting from Jiang & Shu
(1996), Balsara & Shu (2000) and continuing through Balsara et al. (2016),
where it was presented in its most polished form. The WENO interpolation, as
it relates to AFD-WENO, has also been recently documented in a very polished
form in Sections 3 and 4 of Balsara et al. (2024a).

Figure 3: Panel shows part of the mesh around zone boundary “$i+1/2$”. The
fluxes are evaluated pointwise at the zone centers, as shown by the thick
dots. The zone boundaries are shown by the vertical lines. The figure also
shows the stencils associated with the zone boundary “$i+1/2$” for the third
and fifth order AFD-WENO schemes. We have two smaller third order stencils and
a large sixth order stencil. For a third order AFD-WENO scheme, the two
smaller stencils can be non-linearly hybridized. In that case, the second
derivatives of the flux can be obtained at the zone boundary when the
smoothness in the solution warrants it. For fifth order AFD-WENO, the two
smaller stencils can be non-linearly hybridized along with the larger stencil.
In that case, the second and fourth derivatives of the flux can be obtained at
the zone boundary when the smoothness in the solution warrants it. The process
described here can be done for Adaptive Order and Multiresolution WENO
interpolation.

Let us say that we have a high order pointwise WENO interpolation strategy
that is applied to the solution vector of the mesh. At each zone boundary, say
“$i+1/2$” we will then have high order interpolants
$\mathbf{\hat{U}}_{i+1/2}^{-}$ and $\mathbf{\hat{U}}_{i+1/2}^{+}$. Fig. 2
shows us that the same stencils can be used for interpolation and
reconstruction. Say also that we invoke the Riemann solver (in pointwise
fashion) at each zone boundary to obtain
${{\mathbf{F}}^{*}}\left(\mathbf{\hat{U}}_{i+1/2}^{-},\mathbf{\hat{U}}_{i+1/2}^{+}\right)$.
The flux from the previous sentence will indeed be a suitably high order flux
at the zone boundary “$i+1/2$”. However, say that we naively assert a discrete
in space but continuous in time update in the zone “$i$” of the form

${{\partial}_{t}}{{\mathbf{U}}_{i}}=-\frac{1}{\Delta
x}\left({{\mathbf{F}}^{*}}\left(\mathbf{\hat{U}}_{i+1/2}^{-},\mathbf{\hat{U}}_{i+1/2}^{+}\right)-{{\mathbf{F}}^{*}}\left(\mathbf{\hat{U}}_{i-1/2}^{-},\mathbf{\hat{U}}_{i-1/2}^{+}\right)\right).$
(32)

We would find, to our chagrin, that eqn. (32) only has second order spatial
accuracy. Even if the solution had been very smooth, and even if the
interpolation had been carried out with very high order accuracy, the above
equation would only result in a spatially second order scheme. Understanding
why this is so will indeed show us the way out of this dilemma. We illustrate
this for the simplest case where we assume that we are trying to obtain a
third order accurate scheme. Because we have assumed a very smooth solution
and a very smooth flux, we can make the Taylor series expansion:-

$\mathbf{F}\left(x\right)={{f}_{0}}+x{{\left({{\partial}_{x}}f\right)}_{0}}+\frac{{{x}^{2}}}{2}{{\left(\partial_{x}^{2}f\right)}_{0}}+\frac{{{x}^{3}}}{6}{{\left(\partial_{x}^{3}f\right)}_{0}}+...$
(33)

All the terms of the Taylor series, ${{f}_{0}}$,
${{\left({{\partial}_{x}}f\right)}_{0}}$,
${{\left(\partial_{x}^{2}f\right)}_{0}}$,
${{\left(\partial_{x}^{3}f\right)}_{0}}$are all evaluated at $x=0$. We can
evaluate eqn. (33) and its higher derivatives at $x=\pm\Delta x/2$ to get

$\begin{aligned} &{{\left.\mathbf{F}\left(x\right)\right|}_{x=\pm\Delta
x/2}}={{f}_{0}}\pm\frac{\Delta
x}{2}{{\left({{\partial}_{x}}f\right)}_{0}}+\frac{\Delta{{x}^{2}}}{8}{{\left(\partial_{x}^{2}f\right)}_{0}}\pm\frac{\Delta{{x}^{3}}}{48}{{\left(\partial_{x}^{3}f\right)}_{0}}+...\quad;\quad\\\
&{{\left.{{\partial}_{x}}\mathbf{F}\left(x\right)\right|}_{x=\pm\Delta
x/2}}={{\left({{\partial}_{x}}f\right)}_{0}}\pm\frac{\Delta
x}{2}{{\left(\partial_{x}^{2}f\right)}_{0}}+\frac{\Delta{{x}^{2}}}{8}{{\left(\partial_{x}^{3}f\right)}_{0}}\pm...\quad;\quad\\\
&{{\left.\partial_{x}^{2}\mathbf{F}\left(x\right)\right|}_{x=\pm\Delta
x/2}}={{\left(\partial_{x}^{2}f\right)}_{0}}\pm\frac{\Delta
x}{2}{{\left(\partial_{x}^{3}f\right)}_{0}}+...\end{aligned}.$ (34)

Finite differencing the point values of the fluxes at $x=\pm\Delta x/2$
gives:-

$\frac{1}{\Delta x}\left[{{\left.\mathbf{F}\left(x\right)\right|}_{x=\Delta
x/2}}-{{\left.\mathbf{F}\left(x\right)\right|}_{x=-\Delta
x/2}}\right]={{\left({{\partial}_{x}}f\right)}_{0}}+\frac{\Delta{{x}^{2}}}{24}{{\left(\partial_{x}^{3}f\right)}_{0}}$
(35)

The error term, which is proportional to $\Delta{{x}^{2}}$ in eqn. (35), tells
us that the scheme is only second order accurate regardless of the accuracy of
the interpolation. Now realize that if this were a good $k^{th}$ order scheme,
it should have returned ${{\left({{\partial}_{x}}f\right)}_{0}}$ at $x=0$ with
an error term that is proportional to $\Delta{{x}^{k}}$. That is indeed the
true meaning of finite differencing the right hand side of eqn. (18). Instead,
the presence of
${\Delta{{x}^{2}}{{\left(\partial_{x}^{3}f\right)}_{0}}}/{24}\;$ prevents eqn.
(35) from even being a third order accurate expression. However, realize that
a third derivative at the zone center is equivalent to finite differencing two
second derivatives that are evaluated at the zone boundaries. Thus we can take
our numerical fluxes at $x=\pm\Delta x/2$ to be

$\mathbf{F}_{x=\pm\Delta
x/2}^{num}=\left[{{\left.\mathbf{F}\left(x\right)\right|}_{x=\pm\Delta
x/2}}\right]-\frac{\Delta{{x}^{2}}}{24}\left[{{\left.\partial_{x}^{2}\mathbf{F}\left(x\right)\right|}_{x=\pm\Delta
x/2}}\right]$ (36)

By finite differencing the numerical fluxes from eqn. (36) we can easily see
that

$\frac{1}{\Delta x}\left[\mathbf{F}_{x=\Delta x/2}^{num}-\mathbf{F}_{x=-\Delta
x/2}^{num}\right]={{\left({{\partial}_{x}}f\right)}_{0}}+O\left(\Delta{{x}^{4}}\right)$
(37)

In other words, the correction term in eqn. (36) with the higher order flux
derivative was essential for restoring accuracy.

Appendix A of Balsara et al. (2024a) shows us how to design an AFD-WENO scheme
for conservation laws up to any desired accuracy. For up to ninth order of
accuracy, the resulting AFD-WENO scheme is explicitly given by

$\displaystyle{{\partial}_{t}}{{\mathbf{U}}_{i}}=-\frac{1}{\Delta
x}\left\\{{{\mathbf{F}}^{*}}\left(\mathbf{\hat{U}}_{i+1/2}^{-},\mathbf{\hat{U}}_{i+1/2}^{+}\right)-{{\mathbf{F}}^{*}}\left(\mathbf{\hat{U}}_{i-1/2}^{-},\mathbf{\hat{U}}_{i-1/2}^{+}\right)\right\\}$
(38) $\displaystyle-\frac{1}{\Delta
x}\Bigg{\\{}\left[-\frac{1}{24}{{\left(\Delta
x\right)}^{2}}{{\left[\partial_{x}^{2}\mathbf{F}\right]}_{i+1/2}}+{\color[rgb]{1,0,0}\frac{7}{5760}{{\left(\Delta
x\right)}^{4}}{{\left[\partial_{x}^{4}\mathbf{F}\right]}_{i+1/2}}}-{\color[rgb]{0,0,1}\frac{31}{967680}{{\left(\Delta
x\right)}^{6}}{{\left[\partial_{x}^{6}\mathbf{F}\right]}_{i+1/2}}}+{\color[rgb]{1,0,1}\frac{127}{154828800}{{\left(\Delta
x\right)}^{8}}{{\left[\partial_{x}^{8}\mathbf{F}\right]}_{i+1/2}}}\right]$
$\displaystyle\quad\quad-\left[-\frac{1}{24}{{\left(\Delta
x\right)}^{2}}{{\left[\partial_{x}^{2}\mathbf{F}\right]}_{i-1/2}}+{\color[rgb]{1,0,0}\frac{7}{5760}{{\left(\Delta
x\right)}^{4}}{{\left[\partial_{x}^{4}\mathbf{F}\right]}_{i-1/2}}}-{\color[rgb]{0,0,1}\frac{31}{967680}{{\left(\Delta
x\right)}^{6}}{{\left[\partial_{x}^{6}\mathbf{F}\right]}_{i-1/2}}}+{\color[rgb]{1,0,1}\frac{127}{154828800}{{\left(\Delta
x\right)}^{8}}{{\left[\partial_{x}^{8}\mathbf{F}\right]}_{i-1/2}}}\right]\Bigg{\\}}$

Notice that eqn. (38) is still in conservation form, and therefore, it should
be able to capture shock locations accurately. The black terms in the above
equation yield a third order scheme. In that case, the second derivatives of
$\mathbf{F}$ have to be evaluated at the zone boundaries with a WENO
interpolation scheme that is at least third order accurate. If the red terms
are also included, in addition to the black terms, the scheme becomes fifth
order accurate. In that case, all the derivatives of $\mathbf{F}$ have to be
evaluated at the zone boundaries with a WENO interpolation scheme that is at
least fifth order accurate. If the blue terms are also included, in addition
to the black and red terms, we get a seventh order scheme. In that case, all
the derivatives of $\mathbf{F}$ have to be evaluated at the zone boundaries
with a WENO interpolation scheme that is at least seventh order accurate. If
the magenta terms are included, in addition to the black, red and blue terms,
we get a ninth order scheme. In that case, all the derivatives of $\mathbf{F}$
have to be evaluated at the zone boundaries with a WENO interpolation scheme
that is at least ninth order accurate. The WENO-based process by which these
higher order derivatives can be obtained without engendering any Gibbs
phenomenon is described in Section 4 of Balsara et al. (2024a). The stencils
that are used within a fifth order scheme to obtain the higher order flux
derivatives at the zone boundaries are also shown schematically in Fig. 3 of
this paper. This completes our discussion of AFD-WENO schemes for conservation
laws.

###  3.4 AFD-WENO Schemes for Hyperbolic Systems with Non-Conservative
Products

The key purpose of the previous section was to lead us to the AFD-WENO schemes
for hyperbolic systems with some flux terms and some non-conservative
products. Since first order treatments of GR produce hyperbolic systems that
are of this type, such AFD-WENO schemes are most useful for GR. Such schemes
have been obtained very recently by Balsara et al. (2024b) and we briefly show
the reader why they work in this Sub-section. In one dimension, such systems
can be formally written as:-

$\displaystyle{{\partial}_{t}}\mathbf{U}+{{\partial}_{x}}\mathbf{F}\left(\mathbf{U}\right)+\mathbf{C}\left(\mathbf{U}\right){{\partial}_{x}}\mathbf{U}=0.$
(39)

Here $\mathbf{C}\left(\mathbf{U}\right)$ is a solution-dependent matrix of
non-conservative products.

The trick for the derivation is to start with a conservation law of the same
form as eqn. (18) but to rewrite it differently as

$\displaystyle{{\partial}_{t}}\mathbf{U}+{{\partial}_{x}}\mathbf{F}\left(\mathbf{U}\right)=0\qquad\text{
with the flux splitting
}\qquad\mathbf{F}\left(\mathbf{U}\right)={{\mathbf{F}}_{C}}\left(\mathbf{U}\right)+{{\mathbf{F}}_{NC}}\left(\mathbf{U}\right).$
(40)

In other words, we will be effecting a trompe l’oeil where we will treat the
flux ${{\mathbf{F}}_{C}}\left(\mathbf{U}\right)$ in conservation form but we
will write ${{\mathbf{F}}_{NC}}\left(\mathbf{U}\right)$ as if it can only be
written as a non-conservative product. The Jacobians of the fluxes can now be
written as:-

$\displaystyle\mathbf{A}\left(\mathbf{U}\right)=\frac{\partial\left({{\mathbf{F}}_{C}}\left(\mathbf{U}\right)+{{\mathbf{F}}_{NC}}\left(\mathbf{U}\right)\right)}{\partial\mathbf{U}}=\mathbf{B}\left(\mathbf{U}\right)+\mathbf{C}\left(\mathbf{U}\right)\quad\text{with}\quad\mathbf{B}\left(\mathbf{U}\right)\equiv\frac{\partial{{\mathbf{F}}_{C}}\left(\mathbf{U}\right)}{\partial\mathbf{U}}\quad\text{and}\quad\mathbf{C}\left(\mathbf{U}\right)\equiv\frac{\partial{{\mathbf{F}}_{NC}}\left(\mathbf{U}\right)}{\partial\mathbf{U}}.$
(41)

The upshot is that eqn. (40) can now be written in the following equivalent
forms:-

$\displaystyle{{\partial}_{t}}\mathbf{U}+{{\partial}_{x}}\mathbf{F}\left(\mathbf{U}\right)=0\quad\Leftrightarrow\quad{{\partial}_{t}}\mathbf{U}+\mathbf{A}\left(\mathbf{U}\right){{\partial}_{x}}\mathbf{U}=0\quad\Leftrightarrow\quad{{\partial}_{t}}\mathbf{U}+{{\partial}_{x}}{{\mathbf{F}}_{C}}\left(\mathbf{U}\right)+\mathbf{C}\left(\mathbf{U}\right){{\partial}_{x}}\mathbf{U}=0.$
(42)

It is now easy to see that if we initially obtain a scheme for the rightmost
equation in eqn. (42), then we can eventually obtain a scheme for eqn. (39)
just by dropping the subscript “$C$” in
${{\mathbf{F}}_{C}}\left(\mathbf{U}\right)$.

We start our derivation by writing the resolved fluxes
${{\mathbf{F}}^{*}}\left(\mathbf{\hat{U}}_{i+1/2}^{-},\mathbf{\hat{U}}_{i+1/2}^{+}\right)$
and
${{\mathbf{F}}^{*}}\left(\mathbf{\hat{U}}_{i-1/2}^{-},\mathbf{\hat{U}}_{i-1/2}^{+}\right)$
in the form of fluctuations by using the following definitions for the left-
going and right-going fluctuations:-

$\begin{aligned}
{{\mathbf{F}}^{*}}\left(\mathbf{\hat{U}}_{i+1/2}^{-},\mathbf{\hat{U}}_{i+1/2}^{+}\right)&=\mathbf{D}^{*-}\left(\mathbf{\hat{U}}_{i+1/2}^{-},\mathbf{\hat{U}}_{i+1/2}^{+}\right)+{{\mathbf{F}}_{C}}\left(\mathbf{\hat{U}}_{i+1/2}^{-}\right)+{{\mathbf{F}}_{NC}}\left(\mathbf{\hat{U}}_{i+1/2}^{-}\right)\quad\text{
;}\\\
{{\mathbf{F}}^{*}}\left(\mathbf{\hat{U}}_{i-1/2}^{-},\mathbf{\hat{U}}_{i-1/2}^{+}\right)&=-\mathbf{D}^{*+}\left(\mathbf{\hat{U}}_{i-1/2}^{-},\mathbf{\hat{U}}_{i-1/2}^{+}\right)+{{\mathbf{F}}_{C}}\left(\mathbf{\hat{U}}_{i-1/2}^{+}\right)+{{\mathbf{F}}_{NC}}\left(\mathbf{\hat{U}}_{i-1/2}^{+}\right)\end{aligned}.$
(43)

Putting eqn. (43) in eqn. (38) allows us to write the intermediate equation:-

$\displaystyle{{\partial}_{t}}{{\mathbf{U}}_{i}}=-\frac{1}{\Delta
x}\left\\{\mathbf{D}^{*-}\left(\mathbf{\hat{U}}_{i+1/2}^{-},\mathbf{\hat{U}}_{i+1/2}^{+}\right)+\mathbf{D}^{*+}\left(\mathbf{\hat{U}}_{i-1/2}^{-},\mathbf{\hat{U}}_{i-1/2}^{+}\right)\right\\}$
(44) $\displaystyle\quad-\frac{1}{\Delta
x}\left\\{{{\mathbf{F}}_{C}}\left(\mathbf{\hat{U}}_{i+1/2}^{-}\right)-{{\mathbf{F}}_{C}}\left(\mathbf{\hat{U}}_{i-1/2}^{+}\right)\right\\}-\frac{1}{\Delta
x}\left\\{{{\mathbf{F}}_{NC}}\left(\mathbf{\hat{U}}_{i+1/2}^{-}\right)-{{\mathbf{F}}_{NC}}\left(\mathbf{\hat{U}}_{i-1/2}^{+}\right)\right\\}$
$\displaystyle\quad-\frac{1}{\Delta
x}\Bigg{\\{}\left[-{\color[rgb]{0,0,0}\frac{1}{24}{{\left(\Delta
x\right)}^{2}}{{\left[\partial_{x}^{2}\mathbf{F}\right]}_{i+1/2}}}+{\color[rgb]{1,0,0}\frac{7}{5760}{{\left(\Delta
x\right)}^{4}}{{\left[\partial_{x}^{4}\mathbf{F}\right]}_{i+1/2}}}-{\color[rgb]{0,0,1}\frac{31}{967680}{{\left(\Delta
x\right)}^{6}}{{\left[\partial_{x}^{6}\mathbf{F}\right]}_{i+1/2}}}+{\color[rgb]{1,0,1}\frac{127}{154828800}{{\left(\Delta
x\right)}^{8}}{{\left[\partial_{x}^{8}\mathbf{F}\right]}_{i+1/2}}}\right]$
$\displaystyle\quad\qquad\ -\left[-\frac{1}{24}{{\left(\Delta
x\right)}^{2}}{{\left[\partial_{x}^{2}\mathbf{F}\right]}_{i-1/2}}+{\color[rgb]{1,0,0}\frac{7}{5760}{{\left(\Delta
x\right)}^{4}}{{\left[\partial_{x}^{4}\mathbf{F}\right]}_{i-1/2}}}-{\color[rgb]{0,0,1}\frac{31}{967680}{{\left(\Delta
x\right)}^{6}}{{\left[\partial_{x}^{6}\mathbf{F}\right]}_{i-1/2}}}+{\color[rgb]{1,0,1}\frac{127}{154828800}{{\left(\Delta
x\right)}^{8}}{{\left[\partial_{x}^{8}\mathbf{F}\right]}_{i-1/2}}}\right]\Bigg{\\}}$

Now, we will write the term
$\left\\{{{\mathbf{F}}_{NC}}\left(\mathbf{\hat{U}}_{i+1/2}^{-}\right)-{{\mathbf{F}}_{NC}}\left(\mathbf{\hat{U}}_{i-1/2}^{+}\right)\right\\}$
as a term that looks like
$\mathbf{C}\left({{\mathbf{U}}_{i}}\right){{\left({{\partial}_{x}}\mathbf{\hat{U}}\right)}_{i}}$
along with some higher order terms which will eventually cancel off. We write

$\displaystyle-$ $\displaystyle\frac{1}{\Delta
x}\left\\{{{\mathbf{F}}_{NC}}\left(\mathbf{\hat{U}}_{i+1/2}^{-}\right)-{{\mathbf{F}}_{NC}}\left(\mathbf{\hat{U}}_{i-1/2}^{+}\right)\right\\}$
(45)
$\displaystyle\cong-{{\left({{\partial}_{x}}{{\mathbf{F}}_{NC}}\right)}_{i}}-\left\\{+\frac{\Delta{{x}^{2}}}{24}{{\left(\partial_{x}^{3}{{\mathbf{F}}_{NC}}\right)}_{i}}+{\color[rgb]{1,0,0}\frac{\Delta{{x}^{4}}}{1920}{{\left(\partial_{x}^{5}{{\mathbf{F}}_{NC}}\right)}_{i}}}+{\color[rgb]{0,0,1}\frac{\Delta{{x}^{6}}}{322560}{{\left(\partial_{x}^{7}{{\mathbf{F}}_{NC}}\right)}_{i}}}+{\color[rgb]{1,0,1}\frac{\Delta{{x}^{8}}}{92897280}{{\left(\partial_{x}^{9}{{\mathbf{F}}_{NC}}\right)}_{i}}}\right\\}$
$\displaystyle\cong-\mathbf{C}\left({{\mathbf{U}}_{i}}\right){{\left({{\partial}_{x}}\mathbf{\hat{U}}\right)}_{i}}-\frac{1}{\Delta
x}\Bigg{\\{}\left[\frac{1}{24}{{\left(\Delta
x\right)}^{2}}{{\left[\partial_{x}^{2}{{\mathbf{F}}_{NC}}\right]}_{i+1/2}}-{\color[rgb]{1,0,0}\frac{7}{5760}{{\left(\Delta
x\right)}^{4}}{{\left[\partial_{x}^{4}{{\mathbf{F}}_{NC}}\right]}_{i+1/2}}}\right.$
$\displaystyle\left.\qquad\qquad\qquad\qquad\qquad\quad+{\color[rgb]{0,0,1}\frac{31}{967680}{{\left(\Delta
x\right)}^{6}}{{\left[\partial_{x}^{6}{{\mathbf{F}}_{NC}}\right]}_{i+1/2}}}-{\color[rgb]{1,0,1}\frac{127}{154828800}{{\left(\Delta
x\right)}^{8}}{{\left[\partial_{x}^{8}{{\mathbf{F}}_{NC}}\right]}_{i+1/2}}}\right]$
$\displaystyle\qquad\qquad\qquad\qquad\qquad-\left[\frac{1}{24}{{\left(\Delta
x\right)}^{2}}{{\left[\partial_{x}^{2}{{\mathbf{F}}_{NC}}\right]}_{i-1/2}}-{\color[rgb]{1,0,0}\frac{7}{5760}{{\left(\Delta
x\right)}^{4}}{{\left[\partial_{x}^{4}{{\mathbf{F}}_{NC}}\right]}_{i-1/2}}}\right.$
$\displaystyle\left.\qquad\qquad\qquad\qquad\qquad\quad+{\color[rgb]{0,0,1}\frac{31}{967680}{{\left(\Delta
x\right)}^{6}}{{\left[\partial_{x}^{6}{{\mathbf{F}}_{NC}}\right]}_{i-1/2}}}-{\color[rgb]{1,0,1}\frac{127}{154828800}{{\left(\Delta
x\right)}^{8}}{{\left[\partial_{x}^{8}{{\mathbf{F}}_{NC}}\right]}_{i-1/2}}}\right]\Bigg{\\}}$

We insert eqn. (45) in eqn. (44) and use the fact that
${{\mathbf{F}}_{C}}\left(\mathbf{U}\right)=\mathbf{F}\left(\mathbf{U}\right)-{{\mathbf{F}}_{NC}}\left(\mathbf{U}\right)$
to get the penultimate form of our update equation:-

$\displaystyle{{\partial}_{t}}{{\mathbf{U}}_{i}}=-\frac{1}{\Delta
x}\left\\{\mathbf{D}^{*-}\left(\mathbf{\hat{U}}_{i+1/2}^{-},\mathbf{\hat{U}}_{i+1/2}^{+}\right)+\mathbf{D}^{*+}\left(\mathbf{\hat{U}}_{i-1/2}^{-},\mathbf{\hat{U}}_{i-1/2}^{+}\right)\right\\}-\frac{1}{\Delta
x}\left\\{{{\mathbf{F}}_{C}}\left(\mathbf{\hat{U}}_{i+1/2}^{-}\right)-{{\mathbf{F}}_{C}}\left(\mathbf{\hat{U}}_{i-1/2}^{+}\right)\right\\}-\mathbf{C}\left({{\mathbf{U}}_{i}}\right){{\left({{\partial}_{x}}\mathbf{\hat{U}}\right)}_{i}}$
(46) $\displaystyle-\frac{1}{\Delta
x}\Bigg{\\{}\left[-\frac{1}{24}{{\left(\Delta
x\right)}^{2}}{{\left[\partial_{x}^{2}{{\mathbf{F}_{C}}}\right]}_{i+1/2}}+{\color[rgb]{1,0,0}\frac{7}{5760}{{\left(\Delta
x\right)}^{4}}{{\left[\partial_{x}^{4}{{\mathbf{F}_{C}}}\right]}_{i+1/2}}}-{\color[rgb]{0,0,1}\frac{31}{967680}{{\left(\Delta
x\right)}^{6}}{{\left[\partial_{x}^{6}{{\mathbf{F}}_{C}}\right]}_{i+1/2}}}+{\color[rgb]{1,0,1}\frac{127}{154828800}{{\left(\Delta
x\right)}^{8}}{{\left[\partial_{x}^{8}{{\mathbf{F}}_{C}}\right]}_{i+1/2}}}\right]$
$\displaystyle\quad\ \ -\left[-\frac{1}{24}{{\left(\Delta
x\right)}^{2}}{{\left[\partial_{x}^{2}{{\mathbf{F}}_{C}}\right]}_{i-1/2}}+{\color[rgb]{1,0,0}\frac{7}{5760}{{\left(\Delta
x\right)}^{4}}{{\left[\partial_{x}^{4}{{\mathbf{F}}_{C}}\right]}_{i-1/2}}}-{\color[rgb]{0,0,1}\frac{31}{967680}{{\left(\Delta
x\right)}^{6}}{{\left[\partial_{x}^{6}{{\mathbf{F}}_{C}}\right]}_{i-1/2}}}+{\color[rgb]{1,0,1}\frac{127}{154828800}{{\left(\Delta
x\right)}^{8}}{{\left[\partial_{x}^{8}{{\mathbf{F}}_{C}}\right]}_{i-1/2}}}\right]\Bigg{\\}}$

Next, we will derive the final forms of the equations that we need.

The above equation would give us the update strategy for the last equation in
eqn. (42). But we want to transition from the last equation in eqn. (42) to
eqn. (40). This is easily accomplished by dropping the subscript “$C$” in
${{\mathbf{F}}_{C}}\left(\mathbf{U}\right)$ in eqn. (46) to get one of our
final equations as:-

$\displaystyle{{\partial}_{t}}{{\mathbf{U}}_{i}}=-\frac{1}{\Delta
x}\left\\{\mathbf{D}^{*-}\left(\mathbf{\hat{U}}_{i+1/2}^{-},\mathbf{\hat{U}}_{i+1/2}^{+}\right)+\mathbf{D}^{*+}\left(\mathbf{\hat{U}}_{i-1/2}^{-},\mathbf{\hat{U}}_{i-1/2}^{+}\right)\right\\}-\frac{1}{\Delta
x}\left\\{{{\mathbf{F}}}\left(\mathbf{\hat{U}}_{i+1/2}^{-}\right)-{{\mathbf{F}}}\left(\mathbf{\hat{U}}_{i-1/2}^{+}\right)\right\\}-\mathbf{C}\left({{\mathbf{U}}_{i}}\right){{\left({{\partial}_{x}}\mathbf{\hat{U}}\right)}_{i}}$
(47) $\displaystyle\quad-\frac{1}{\Delta
x}\Bigg{\\{}\left[-\frac{1}{24}{{\left(\Delta
x\right)}^{2}}{{\left[\partial_{x}^{2}{{\mathbf{F}}}\right]}_{i+1/2}}+{\color[rgb]{1,0,0}\frac{7}{5760}{{\left(\Delta
x\right)}^{4}}{{\left[\partial_{x}^{4}{{\mathbf{F}}}\right]}_{i+1/2}}}-{\color[rgb]{0,0,1}\frac{31}{967680}{{\left(\Delta
x\right)}^{6}}{{\left[\partial_{x}^{6}{{\mathbf{F}}}\right]}_{i+1/2}}}+{\color[rgb]{1,0,1}\frac{127}{154828800}{{\left(\Delta
x\right)}^{8}}{{\left[\partial_{x}^{8}{{\mathbf{F}}}\right]}_{i+1/2}}}\right]$
$\displaystyle\qquad\quad-\left[-\frac{1}{24}{{\left(\Delta
x\right)}^{2}}{{\left[\partial_{x}^{2}{{\mathbf{F}}}\right]}_{i-1/2}}+{\color[rgb]{1,0,0}\frac{7}{5760}{{\left(\Delta
x\right)}^{4}}{{\left[\partial_{x}^{4}{{\mathbf{F}}}\right]}_{i-1/2}}}-{\color[rgb]{0,0,1}\frac{31}{967680}{{\left(\Delta
x\right)}^{6}}{{\left[\partial_{x}^{6}{{\mathbf{F}}}\right]}_{i-1/2}}}+{\color[rgb]{1,0,1}\frac{127}{154828800}{{\left(\Delta
x\right)}^{8}}{{\left[\partial_{x}^{8}{{\mathbf{F}}}\right]}_{i-1/2}}}\right]\Bigg{\\}}$

The above equation is very illustrative. Notice that when the fluxes are
absent, i.e. when the entire system is represented by the non-conservative
products, the above equation reduces to

$\displaystyle{{\partial}_{t}}{{\mathbf{U}}_{i}}=-\frac{1}{\Delta
x}\left\\{\mathbf{D}^{*-}\left(\mathbf{\hat{U}}_{i+1/2}^{-},\mathbf{\hat{U}}_{i+1/2}^{+}\right)+\mathbf{D}^{*+}\left(\mathbf{\hat{U}}_{i-1/2}^{-},\mathbf{\hat{U}}_{i-1/2}^{+}\right)\right\\}-\mathbf{C}\left({{\mathbf{U}}_{i}}\right){{\left({{\partial}_{x}}\mathbf{\hat{U}}\right)}_{i}}$
(48)

which is almost identical to eqn. (31). (Recall that
$\mathbf{C}\left({{\mathbf{U}}_{i}}\right)\to\mathbf{A}\left({{\mathbf{U}}_{i}}\right)$
in that limit.) In other words, when solving just the Einstein equations,
without additional hydrodynamic equations, AFD-WENO and FD-WENO become
practically identical to one another. Of course, eqn. (31) relies on
reconstruction whereas eqn. (48) relies on interpolation. A form that is
completely equivalent to eqn. (47) at the analytical level, but has slightly
better stability properties, can be written as:-

$\displaystyle{{\partial}_{t}}{{\mathbf{U}}_{i}}=-\frac{1}{\Delta
x}\left\\{\mathbf{D}_{i+1/2}^{*-}\left(\mathbf{\hat{U}}_{i+1/2}^{-},\mathbf{\hat{U}}_{i+1/2}^{+}\right)+\mathbf{D}_{i-1/2}^{*+}\left(\mathbf{\hat{U}}_{i-1/2}^{-},\mathbf{\hat{U}}_{i-1/2}^{+}\right)\right\\}-\frac{1}{\Delta
x}\left\\{\mathbf{F}\left(\mathbf{\hat{U}}_{i+1/2}^{-}\right)-\mathbf{F}\left(\mathbf{\hat{U}}_{i-1/2}^{+}\right)\right\\}$
(49) $\displaystyle\quad-\frac{1}{\Delta
x}\mathbf{C}\left({{\mathbf{U}}_{i}}\right)\left(\mathbf{\hat{U}}_{i+1/2}^{-}-\mathbf{\hat{U}}_{i-1/2}^{+}\right)-\frac{1}{\Delta
x}\mathbf{C}\left({{\mathbf{U}}_{i}}\right)\Bigg{\\{}\left[-\frac{1}{24}{{\left(\Delta
x\right)}^{2}}{{\left[\partial_{x}^{2}{\mathbf{\hat{U}}}\right]}_{i+1/2}}+{\color[rgb]{1,0,0}\frac{7}{5760}{{\left(\Delta
x\right)}^{4}}{{\left[\partial_{x}^{4}\mathbf{\hat{U}}\right]}_{i+1/2}}}\right.$
$\displaystyle\left.\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad-{\color[rgb]{0,0,1}\frac{31}{967680}{{\left(\Delta
x\right)}^{6}}{{\left[\partial_{x}^{6}\mathbf{\hat{U}}\right]}_{i+1/2}}}+{\color[rgb]{1,0,1}\frac{127}{154828800}{{\left(\Delta
x\right)}^{8}}{{\left[\partial_{x}^{8}\mathbf{\hat{U}}\right]}_{i+1/2}}}\right]$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad-\left[-\frac{1}{24}{{\left(\Delta
x\right)}^{2}}{{\left[\partial_{x}^{2}{\mathbf{\hat{U}}}\right]}_{i-1/2}}+{\color[rgb]{1,0,0}\frac{7}{5760}{{\left(\Delta
x\right)}^{4}}{{\left[\partial_{x}^{4}\mathbf{\hat{U}}\right]}_{i-1/2}}}\right.$
$\displaystyle\left.\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad-{\color[rgb]{0,0,1}\frac{31}{967680}{{\left(\Delta
x\right)}^{6}}{{\left[\partial_{x}^{6}\mathbf{\hat{U}}\right]}_{i-1/2}}}+{\color[rgb]{1,0,1}\frac{127}{154828800}{{\left(\Delta
x\right)}^{8}}{{\left[\partial_{x}^{8}\mathbf{\hat{U}}\right]}_{i-1/2}}}\right]\Bigg{\\}}$
$\displaystyle\quad-\frac{1}{\Delta
x}\Bigg{\\{}\left[-\frac{1}{24}{{\left(\Delta
x\right)}^{2}}{{\left[\partial_{x}^{2}{{\mathbf{F}}}\right]}_{i+1/2}}+{\color[rgb]{1,0,0}\frac{7}{5760}{{\left(\Delta
x\right)}^{4}}{{\left[\partial_{x}^{4}{{\mathbf{F}}}\right]}_{i+1/2}}}-{\color[rgb]{0,0,1}\frac{31}{967680}{{\left(\Delta
x\right)}^{6}}{{\left[\partial_{x}^{6}{{\mathbf{F}}}\right]}_{i+1/2}}}+{\color[rgb]{1,0,1}\frac{127}{154828800}{{\left(\Delta
x\right)}^{8}}{{\left[\partial_{x}^{8}{{\mathbf{F}}}\right]}_{i+1/2}}}\right]$
$\displaystyle\qquad\quad-\left[-\frac{1}{24}{{\left(\Delta
x\right)}^{2}}{{\left[\partial_{x}^{2}{{\mathbf{F}}}\right]}_{i-1/2}}+{\color[rgb]{1,0,0}\frac{7}{5760}{{\left(\Delta
x\right)}^{4}}{{\left[\partial_{x}^{4}{{\mathbf{F}}}\right]}_{i-1/2}}}-{\color[rgb]{0,0,1}\frac{31}{967680}{{\left(\Delta
x\right)}^{6}}{{\left[\partial_{x}^{6}{{\mathbf{F}}}\right]}_{i-1/2}}}+{\color[rgb]{1,0,1}\frac{127}{154828800}{{\left(\Delta
x\right)}^{8}}{{\left[\partial_{x}^{8}{{\mathbf{F}}}\right]}_{i-1/2}}}\right]\Bigg{\\}}$

This completes our derivation of AFD-WENO schemes for systems with non-
conservative products as they arise in GR. Along with the derivations, we have
striven to provide insights on the strength of each scheme and its suitability
for NR. Notice that the evaluation of the flux derivatives in eqns. (47) or
(49) is the same as was sketched in Fig. 3. The evaluation of the higher
derivatives of the states at the zone boundaries in eqn. (49) also follows a
process that is practically identical to Fig. 3 and is described in detail in
Balsara et al. (2024b).

Some overarching comments that connect the multiple finite difference WENO
formulations are in order.

  1. 1.

Realize that the scheme in Sub-section 3.1 relies on a flux splitting into an
LLF flux and can, therefore, be somewhat dissipative. The schemes described in
Sub-sections 3.2, 3.3 and 3.4 can use improved Riemann solvers. Any good
Riemann solver can be used in the schemes described in those three Sub-
scections. The HLLI Riemann solver from Dumbser & Balsara (2016) might be very
beneficial because it can handle all different types of intermediate waves
including the contact discontinuities that are needed in well-balancing of
hydrodynamic problems.

  2. 2.

The AFD-WENO scheme in this sub-section has an advantage over the FD-WENO
scheme in Sub-section 3.2 in that it recovers the conservation form when such
a form is present. This can be very useful if the hydrodynamical equations are
to be evolved in conjunction with the Einstein equations.

  3. 3.

The term ${{\left({{\partial}_{x}}\mathbf{\hat{U}}\right)}_{i}}$ in eqn. (47)
has to be evaluated using the same WENO interpolation that was used for the
state vector; as a result, eqn. (47) usually shows one lower order of accuracy
compared to eqn. (49). Eqn. (47) may have some advantages for GR nevertheless
because when the solution is smooth, it gives very good quality solutions at a
lower computational cost.

  4. 4.

When strong shocks are present, eqn. (49) has slightly better stability
properties compared to eqn. (47) because the term
$\mathbf{C}\left({{\mathbf{U}}_{i}}\right)\left(\mathbf{\hat{U}}_{i+1/2}^{-}-\mathbf{\hat{U}}_{i-1/2}^{+}\right)$
is itself treated in a finite difference form. However, since strong shocks
are only present in the hydrodynamic sector in a NR calculation, the treatment
of flux terms is identical in both equations.

  5. 5.

We have only addressed spatial accuracy in this Section. Like all finite
difference WENO schemes, the schemes presented here are to be coupled with
strong stability preserving Runge-Kutta (SSP-RK) timestepping. Temporally
high-order SSP-RK timestepping has been described in Shu & Osher (1988) or
Spiteri & Ruuth (2002, 2003).

  6. 6.

In general, the treatment of the source terms in GR is such that the source
terms are always non-stiff. As a result, we shouldn’t need variants of SSP-RK
timestepping that include stiff source terms.

  7. 7.

It would be useful to keep an eye on the ecosystem of ideas that are emerging
around AFD-WENO methods. Such ideas include well-balancing, physical
realizability, and preservation of divergence and curl constraints. These
enhancements retain the low cost of the original AFD-WENO scheme.

This completes our discussion of FD-WENO and AFD-WENO schemes for systems with
non-conservative products as they arise in GR.

###  3.5 Well balancing of FD-WENO and AFD-WENO schemes

Preserving the stationarity of an equilibrium solution $\mathbf{U}_{e}$ is far
from trivial, and, unless suitable steps are undertaken, a general numerical
scheme will not be able to satisfy the condition
$\partial_{t}\mathbf{U}_{e}=0$ exactly over arbitrarily long timescales. This
is particularly frustrating, for instance, if one is interested in studying
numerically the oscillation modes of astrophysical objects in equilibrium. In
these circumstances, if the order of magnitude of the discretization errors
introduced by the numerical scheme are close to those of the physical initial
perturbation under investigation, it will be impossible to distinguish among a
physical effect and a numerical artifact. In the worse cases, the accumulation
of numerical errors can even spoil completely the equilibrium solution over
sufficiently long timescales. Hence, a whole field of research exists to
ensure the so-called _well–balancing_ of a numerical scheme. Initially
proposed in the context of the shallow water equations, well–balanced
numerical schemes have been proposed in a variety of different techniques,
sometimes tailored to specific PDE systems (Bermudez & Vázquez-Cendón, 1994;
LeVeque, 1998; Gosse, 2001; Audusse et al., 2004; Botta et al., 2004; Castro
et al., 2017; Castro & Parés, 2020; Gaburro et al., 2021; Xu & Shu, 2024;
Bhoriya et al., 2024).

Here we follow the pragmatic approach chosen by Dumbser et al. (2024), who,
motivated by the complexity of the Einstein equations, subtracted the
equilibrium solution alltogether from the governing PDE, see also Pareschi &
Rey (2017); Berberich et al. (2021), therefore solving the following augmented
system

$\displaystyle{{\partial}_{t}}\mathbf{U}+{{\partial}_{x}}\mathbf{F}\left(\mathbf{U}\right)-{{\partial}_{x}}\mathbf{F}\left(\mathbf{U}_{e}\right)+\mathbf{C}\left(\mathbf{U}\right){{\partial}_{x}}\mathbf{U}-\mathbf{C}\left(\mathbf{U}_{e}\right){{\partial}_{x}}\mathbf{U}_{e}=\mathbf{S}(\mathbf{U})-\mathbf{S}(\mathbf{U}_{e});\qquad{{\partial}_{t}}\mathbf{U}_{e}=0.$
(50)

In this way, the discretization errors introduced by the numerical scheme are
removed and exact equilibrium can be maintained. There is a price to pay, of
course, in so far the vector of evolved quantities is effectively doubled,
becoming $\tilde{\mathbf{U}}=[\mathbf{U},\mathbf{U}_{e}]$. The equilibrium
sector of the extended vector $\tilde{\mathbf{U}}$, i.e. $\mathbf{U}_{e}$, is
ultimately excluded from the time evolution, but it enters all the spatial
discretization procedures of the numerical scheme.

We have added this property to our numerical schemes as an extra feature that
can be activated, or not, depending on the problem considered. For more
details on the numerical approach, see Dumbser et al. (2024); Pareschi & Rey
(2017); Berberich et al. (2021).

##  4 Numerical tests

In this section, we explore two distinct computational schemes outlined in
Section 3. First, we consider the scheme described by Eqn. (31), explained in
Section 3.2. It is noteworthy that employing a $k^{th}$-order accurate WENO–AO
reconstruction in Eqn. (31) results in a $k^{th}$-order accurate FD–WENO
scheme. We label the resultant scheme as FD–WENO-k. Next, we explore the
AFD–WENO scheme outlined in Section 3.4 and defined by Eqn. (47). The scheme
given by Eqn. (47) usually shows one order of accuracy lower compared to eqn.
(49); nevertheless, it may have some advantages for GR because, when the
solution is smooth, it gives very good quality solutions at a lower
computational cost. Implementing a WENO–AO interpolation of order ’$k$’ in
Eqn. (47) yields a $(k-1)^{th}$ order scheme. We label the resultant scheme as
AFD–WENO-(k-1). For both schemes, we explore the scenarios where $k$ is set to
3, 5, 7, and 9, providing a comprehensive examination across different order
schemes.

In the following, we show the performances of the new schemes over a
representative sample of standard tests in NR. While doing so, we are also
going to monitor the energy constraint $H$ and the momentum constraints
$M_{i}$, also known as ADM constraints, which are defined as

$\displaystyle H$ $\displaystyle=$ $\displaystyle R-K_{ij}K^{ij}+K^{2}\,,$
(51) $\displaystyle M_{i}$ $\displaystyle=$
$\displaystyle\gamma^{jl}\left(\partial_{l}K_{ij}-\partial_{i}K_{jl}-\Gamma^{m}_{jl}K_{mi}+\Gamma^{m}_{ji}K_{ml}\right)\,.$
(52)

###  4.1 Gauge wave

We start our validation by considering the propagation of a gauge wave, which
is generated after a suitable change of coordinates in a flat Minkowski
spacetime (Alcubierre et al., 2004), producing the following metric

$ds^{2}=-H(x,t)\,dt^{2}+H(x,t)\,dx^{2}+dy^{2}+dz^{2},\quad\text{where}\quad
H(x,t)=1-A\,\sin\left(2\pi(x-t)\right)\,.$ (53)

In practice, it describes a sinusoidal wave of amplitude $A$ that propagates
along the $x$-axis, with the lapse $\alpha=\sqrt{H}$ and the extrinsic
curvature initialized through its definition
$K_{ij}=-\partial_{t}\gamma_{ij}/(2\alpha)$. As reported by several authors,
this test cannot be evolved stably and accurately over long timescales by the
the classical BSSNOK system, either in second or first-order formulations
(Babiuc et al., 2008; Brown et al., 2012), while it has been successfully
solved within the first–order CCZ4 and Z4 formulations (Dumbser et al., 2018,
2024), in both cases with no need to activate the damping mechanisms. We have
also solved this problem by selecting $\kappa_{1}=0$, $\kappa_{2}=0$, no
_gamma–driver_ ($s=0$), harmonic gauge conditions, i.e. $g(\alpha)=1$ in Eq
(5), and $c=1$. The numerical domain is the rectangle
$[-0.5,0.5]\times[-0.05,0.05]$, that is covered by a $N_{x}\times N_{y}$
uniform grid. Periodic boundary conditions are selected.

Scheme | Grid size |  $L_{1}-$errors | Order |  $L_{\infty}-$errors | Order
---|---|---|---|---|---
FD–WENO-3 | 32 $\times$ 2 | 1.32436e-04 | – | 2.68393e-04 | –
64 $\times$ 4 | 1.65963e-05 | 3.00 | 3.41976e-05 | 2.97
128 $\times$ 8 | 2.07707e-06 | 3.00 | 4.29149e-06 | 2.99
256 $\times$ 16 | 2.59697e-07 | 3.00 | 5.36977e-07 | 3.00
FD–WENO-5 | 32 $\times$ 2 | 1.41137e-06 | – | 4.11999e-06 | –
64 $\times$ 4 | 4.42906e-08 | 4.99 | 1.31668e-07 | 4.97
128 $\times$ 8 | 1.38883e-09 | 5.00 | 4.13672e-09 | 4.99
256 $\times$ 16 | 4.34472e-11 | 5.00 | 1.29485e-10 | 5.00
FD–WENO-7 | 16 $\times$ 1 | 9.83235e-06 | – | 2.75912E-05 | –
32 $\times$ 2 | 4.15493e-08 | 7.89 | 1.50863E-07 | 7.51
48 $\times$ 3 | 2.43860e-09 | 6.99 | 8.53961E-09 | 7.08
64 $\times$ 4 | 3.29132e-10 | 6.96 | 1.13484E-09 | 7.02
FD–WENO-9 | 8 $\times$ 1 | 1.70885e-04 | – | 3.55896e-04 | –
16 $\times$ 2 | 8.48349e-07 | 7.65 | 2.66925e-06 | 7.06
24 $\times$ 3 | 2.57895e-08 | 8.62 | 9.82247e-08 | 8.14
32 $\times$ 4 | 2.17566e-09 | 8.60 | 8.25465e-09 | 8.61
Table 1: Gauge wave: $L_{1}$ errors, $L_{\infty}$ errors, and the corresponding order of convergence for the third, fifth, seventh and ninth order accurate FD–WENO schemes. The lapse $\alpha$ variable has been shown. A final time of $t=1.0$ has been used for the study. Scheme | Grid size |  $L_{1}-$errors | Order |  $L_{\infty}-$errors | Order
---|---|---|---|---|---
AFD–WENO-2 | 32 $\times$ 2 | 3.45649E-04 | – | 5.68896E-04 | –
64 $\times$ 4 | 8.25823E-05 | 2.07 | 1.34221E-04 | 2.08
128 $\times$ 8 | 2.03870E-05 | 2.02 | 3.30125E-05 | 2.02
256 $\times$ 16 | 5.07952E-06 | 2.00 | 8.21743E-06 | 2.01
AFD–WENO-4 | 32 $\times$ 2 | 3.69051E-06 | – | 8.72790E-06 | –
64 $\times$ 4 | 1.98678E-07 | 4.22 | 4.35846E-07 | 4.32
128 $\times$ 8 | 1.22491E-08 | 4.02 | 2.63489E-08 | 4.05
256 $\times$ 16 | 7.63165E-10 | 4.00 | 1.63497E-09 | 4.01
AFD–WENO-6 | 16 $\times$ 1 | 9.80463e-06 | – | 2.80057e-05 | –
32 $\times$ 2 | 7.37080e-08 | 7.06 | 2.24641e-07 | 6.96
48 $\times$ 3 | 6.23277e-09 | 6.09 | 1.81759e-08 | 6.20
64 $\times$ 4 | 1.13379e-09 | 5.92 | 3.27951e-09 | 5.95
AFD–WENO-8 | 8 $\times$ 1 | 2.77564e-04 | – | 5.40091e-04 | –
16 $\times$ 2 | 8.74626e-07 | 8.31 | 2.78709e-06 | 7.60
24 $\times$ 3 | 3.62176e-08 | 7.85 | 1.31702e-07 | 7.53
32 $\times$ 4 | 4.60690e-09 | 7.17 | 1.68872e-08 | 7.14
Table 2: Gauge wave: $L_{1}$ errors, $L_{\infty}$ errors, and the
corresponding order of convergence for the second, fourth, sixth, and eighth
order accurate AFD–WENO schemes. The lapse $\alpha$ variable has been shown. A
final time of $t=1.0$ has been used for the study.

This test is also rather useful for validating the order of accuracy of the
numerical scheme. We have therefore considered a battery of tests with a short
final time, $t=1.0$, producing results that are reported in Tab. 1–2,
corresponding to our FD–WENO and AFD–WENO schemes, respectively. As it can be
appreciated, the nominal order of convergence is successfully reproduced.

Figure 4, on the contrary, reports the results obtained with FD–WENO schemes,
by showing, in the first column, the profiles of the lapse $\alpha$ in a long
run with a final time $t=1000$, compared to the exact solution, for various
orders, and using a wave amplitude $A=0.1$. In the right column we display the
evolution of the ADM constraints, which remain perfectly under control all
along the simulation. Figure 5 serves the same purpose, but this time for the
results obtained with AFD–WENO schemes. It is interesting to note that on such
a long simulation both FD–WENO-3 and AFD–WENO-2 schemes fail to reproduce the
correct solution (top–left panels of Fig. 4 and Fig. 5), with a slightly
better performance of FD–WENO-3 with respect to AFD–WENO-2. For the higher-
order FD–WENO and AFD–WENO schemes ($k=5,7,9$), the numerical solutions match
perfectly well the analytic ones.

Figure 4: Gauge wave: The left panel (Figs. 4a, 4c, 4e, 4g) shows profiles of
the lapse $\alpha$ at the final time $t=1000$ for various order accurate
FD–WENO schemes. The right panel (Figs. 4b, 4d, 4f, 4h) shows the
$L_{2}$-error evolution of the ADM constraints for the corresponding order
scheme. Figure 5: Gauge wave: The left panel (Figs. 5a, 5c, 5e, 5g) shows
profiles of the lapse $\alpha$ at the final time $t=1000$ for various order
accurate AFD–WENO schemes. The right panel (Figs. 5b, 5d, 5f, 5h) shows the
$L_{2}$-error evolution of the ADM constraints for the corresponding order
scheme.

###  4.2 Robust stability test

The so–called _robust stability test_ is used to provide an empiric indication
of the hyperbolicity of the PDE system, and it was introduced when there were
very few theoretical analysis of the second order BSSNOK formulation (Beyer &
Sarbach, 2004). In view of the huge progress witnessed by NR in the last
twenty years, and especially because strong hyperbolicity of the present
first–order Z4 formulation has been verified by Dumbser et al. (2024), this
test has lost part of its relevance. Nevertheless, it is still quite useful to
exclude possible exponential unstable modes which could be in principle
triggered by a given numerical algorithm (Calabrese et al., 2006). We should
also comment on the fact that, to the best of our knowledge, none of the
first–order formulations of the Einstein equations have shown any problem with
this test (Bona et al., 2004a; Brown et al., 2012; Dumbser et al., 2018,
2024). We have therefore prepared the usual setup, given by a flat Minkowski
spacetime in two space dimensions, in which all the metric terms are perturbed
with a random perturbation of amplitude $\pm 10^{-7}/\varrho^{2}$. The integer
$\rho$ is also used to control the resolution in a set of increasingly refined
meshes. The computational domain is given by the square
$[-0.5;0.5]\times[-0.5;0.5]$. As customary for this test, the _gamma-driver_
is activated, with $\eta=0.2$ in Eq. (7). The time evolution of the Einstein
constraints for the FD-WENO and for the AFD-WENO schemes are shown in Fig. 6
and Fig. 7, respectively, establishing quite a successful outcome for this
test.

Figure 6: Robust stability test: Time evolution of $L_{2}-$norms of Einstein
constraints for the Robust Stability Test with a random initial perturbation
of amplitude $10^{-7}/\rho^{2}$ in all quantities on a sequence of
successively refined meshes on the square domain in 2D. The final time is
$t=1000$, and we used a fifth-order accurate FD-WENO scheme (FD-WENO-5) to
obtain the results. Figure 7: Robust stability test: Time evolution of
$L_{2}-$norms of Einstein constraints for the Robust Stability Test with a
random initial perturbation of amplitude $10^{-7}/\rho^{2}$ in all quantities
on a sequence of successively refined meshes on the square domain in 2D. The
final time is $t=1000$, and we used a fourth-order accurate AFD-WENO scheme
(AFD-WENO-4) to obtain the results.

###  4.3 Gowdy wave

There is a relatively simple solution to test the numerical behavior of the
PDE system in the strong field regime, and this is given by the so–called
Gowdy solution (Gowdy, 1971). It describes a polarized gravitational wave (for
us along the x-direction) via a spacetime metric given by

$ds^{2}=t^{-1/2}e^{Q/2}(-dt^{2}+dx^{2})+t(e^{P}dy^{2}+e^{-P}dz^{2})\,,$ (54)

with

$\displaystyle P(x,t)$ $\displaystyle=$ $\displaystyle J_{0}(2\pi t)\cos(2\pi
x)\,,$ (55) $\displaystyle Q(x,t)$ $\displaystyle=$ $\displaystyle\pi
J_{0}(2\pi)J_{1}(2\pi)-2\pi\,tJ_{0}(2\pi t)J_{1}(2\pi t)\cos^{2}(2\pi
x)+2\pi^{2}t^{2}[J_{0}^{2}(2\pi t)+J_{1}^{2}(2\pi
t)-J_{0}^{2}(2\pi)-J_{1}^{2}(2\pi)]\,,$ (56)

where $J_{0}$ and $J_{1}$ are the Bessel functions (Bona et al., 2004a;
Alcubierre et al., 2004). This test becomes interesting when used to simulate
the formation of a singularity. This is better obtained after performing the
time change of coordinate $t\rightarrow\tau$, i.e.

$t=t_{0}e^{-\tau/\tau_{0}}\,,$ (57)

implying that the singularity is reached in the limit of
$\tau\rightarrow+\infty$, $\tau$ now being the new time coordinate of the
code. Since the roots of $J_{0}(2\pi\,t)$ determines the times at which the
lapse is momentarily constant (in space), we conform to the standard choice
and we select $t_{0}$ as the 20th root of $J_{0}(2\pi\,t)$, giving
$t_{0}=9.8753205829098$. The corresponding value of the lapse resulting from
the coordinate transformation (57) is $\alpha=t^{3/4}\,e^{Q/4}/\tau_{0}$, and,
after imposing $\alpha=1$ at $\tau=0$ ($t=t_{0}$), we get the value
$\tau_{0}=471.806749033034$. The time coordinate transformation does not
affect the metric terms $\gamma_{ij}$, which follow directly from (54), while
the extrinsic curvature can be easly computed from the definition
$K_{ij}=-1/(2\alpha)\,\partial_{\tau}\gamma_{ij}$. The simulation then starts
at $\tau=0$, reaching the singularity in the limit of $\tau\rightarrow+\infty$
($t\rightarrow 0)$.

We have solved this test as a one-dimensional problem using a mesh composed of
120 zones, up to the final time $\tau=1000$. We impose periodic boundary
conditions. In Fig. 8 and Fig. 9 we show the results obtained with our FD-WENO
and AFD-WENO schemes, respectively. The left panels report the variable
$K_{xx}$, and they show a perfect matching to the exact solution at the final
time. The right panels, on the other hand, show the evolution of the Einstein
constraints, which, in spite of a linear increase already documented in the
literature (Bona et al., 2004a), remain under control.

Figure 8: Gowdy wave: The left panel (Figs. 8a,8c,8e,8g) shows profiles of the
$K_{xx}$ variable at the final time $\tau=1000$ for various order accurate
FD–WENO schemes. The right panel (Figs. 8b,8d,8f,8h) shows the $L_{2}$-error
evolution of the ADM constraints at the corresponding order. Figure 9: Gowdy
wave: The left panel (Figs. 9a,9c,9e,9g) shows profiles of the $K_{xx}$
variable at the final time $\tau=1000$ for various order accurate AFD–WENO
schemes. The right panel (Figs. 9b,9d,9f,9h) shows the $L_{2}$-error evolution
of the ADM constraints at the corresponding order.

###  4.4 Stationary black holes

Contrary to naive expectations, keeping a simple Schwarzschild or Kerr black
hole stationary in three dimensions via a NR code is far from trivial. To the
best of our knowledge, this has never been obtained through the second–order
BSSNOK formulation, apart from a promising attempt performed by Alcubierre &
Brügmann (2001), who, however, performed the calculations on one octant only
(because of spherical symmetry) while they reported the appearance of an
unstable mode when they extended the simulation in the full grid space. In
2001, after using an alternative first–order formulation of the Einstein
equations, henceforth the KST formulation, Kidder et al. (2001) were able to
evolve a Schwarzschild black hole in three space dimensions, though the
momentum constraints did not seem to be fully under control. Better results
were obtained by Brügmann (2013), who kept a Schwarzschild black hole in
spherical Kerr–Schild coordinates stationary for very long timescales, using
spectral methods in combination with a first–order formulation of the Einstein
equations in the generalized harmonic gauge. Again in this formulation, Tichy
et al. (2023) obtained rather good results as well for a non-rotating black
hole, but using DG methods. Finally, after using a well–balanced version of
the first–order Z4 formulation of the Einstein equations, in combination with
an ADER–DG scheme, Dumbser et al. (2024) were able to preserve perfect
stationary equilibria of even extreme Kerr black holes over arbitrarily long
timescales.

We have therefore considered the stationary solution represented by the Kerr
spacetime in pseudo-Cartesian Kerr–Schild coordinates, as they can be found,
for instance, in de Felice & Clarke (1990) (Sect. 11.4) or in Wiltshire et al.
(2009). Here, we reconsider that problem by adopting a seventh-order FD-WENO
scheme and a sixth-order AFD-WENO scheme, applied in three-dimensional
(spatial) simulations. The three-dimensional computational domain is given by
$\Omega=[-5;5]^{3}$ covered by a $80\times 80\times 80$ uniform grid, while,
as boundary conditions, we impose the stationary equilibrium solution at the
outer border. We also recall that, in the coordinates adopted, the physical
singularity is mapped into the the so–called ring singularity in z = 0 plane,
with a radius given by the spin of the black hole. For this reason we have
excised the ring singularity from the numerical domain, which remains inside
an appropriate excision box.

Moreover, we have found crucial for a successful evolution of these
equilibrium tests over arbitrarily long times to activate the _well balancing_
property described in Sect. 3.5. Without well balancing, small but ever
growing deviations from equilibrium where detected on times $t\geq 100$, even
at the sixth order of convergence of our numerical schemes. Under these
settings, we have considered the following two representative cases:

  * •

A Schwarzschild black hole ($a=0$). Fig. 10 and Fig. 11 show the results of a
long run with $t_{f}=1000\,M$, obtained with the FD-WENO scheme and with the
AFD-WENO scheme, respectively. The central top panel reports a contour-color
representation of the lapse $\alpha$ through the $z=0$ plane, while the arrows
visualize the vector field of the shift $\beta^{i}$, which, in spite of the
black hole being non/rotating, is non zero due to the choice of the
coordinates. The remaining panels, with 1D profiles along the $z$ axis, show
the comparison with the equilibrium solution and prove the ability of the code
in preserving the stationarity of the solution.

  * •

A Kerr black hole ($a=0.9$). In a rather similar way, Fig. 12 and Fig. 13 show
the results of a long run with $t_{f}=1000\,M$, obtained with the FD-WENO
scheme and with the AFD-WENO scheme, respectively. The vector field of the
shift, represented in the central top panels of those figures, manifests the
clockwise rotation of the black hole.

Finally, the bottom right panels of all figures from Fig. 10 to Fig. 13
monitor the evolution of the Einstein constraints, which remain perfectly
constant all along the simulation. This is a remarkable result that can only
be obtained thanks to the well balancing of the numerical scheme.

Figure 10: Schwarzschild black hole: Figs. 10a-10d show $\alpha$, $K_{zz}$ and
$\gamma_{zz}$ profiles for the Schwarzschild black hole ($a=0$). Fig. 10e
shows the evolution of ADM-constraints in terms of the $L_{2}$-norm. The final
time is $t=1000$, and we have used a seventh-order accurate FD-WENO scheme.
Figure 11: Schwarzschild black hole: Figs. 11a-11d show $\alpha$, $K_{zz}$ and
$\gamma_{zz}$ profiles for the Schwarzschild black hole ($a=0$). Fig. 11e
shows the evolution of ADM-constraints in terms of the $L_{2}$-norm. The final
time is $t=1000$, and we have used a sixth-order accurate AFD-WENO scheme.
Figure 12: Kerr black hole: Figs. 12a-12d show $\alpha$, $K_{zz}$ and
$\gamma_{zz}$ profiles for the Kerr black hole ($a=0.9$). Fig. 12e shows the
evolution of ADM-constraints in terms of the $L_{2}$-norm. The final time is
$t=1000$, and we have used a seventh-order accurate FD-WENO scheme. Figure 13:
Kerr black hole: Figs. 13a-13d show $\alpha$, $K_{zz}$ and $\gamma_{zz}$
profiles for the Kerr black hole ($a=0.9$). Fig. 13e shows the evolution of
ADM-constraints in terms of the $L_{2}$-norm. The final time is $t=1000$, and
we have used a sixth-order accurate AFD-WENO scheme.

###  4.5 Head on collision of two black holes

The correct simulation of binary mergers (either black holes or neutron stars)
is of course the final goal of numerical relativity, which is proving to be an
invaluable scientific tool in the preset gravitational wave era (Abbott et
al., 2019). In this work, still at the level of an exploratory calculation,
before a dedicated analysis is performed, we have considered the head-on
collision of two nonrotating black holes, which are modelled as two moving
punctures. We compute the initial conditions through the TwoPunctures initial
data code of Ansorg et al. (2004), selecting two black holes with equals
masses, $M_{1}=M_{2}=1$, initial positions given by
$\boldsymbol{x}^{-}=(-1,0,0)$ and $\boldsymbol{x}^{+}=(+1,0,0)$, zero
individual spins and zero linear momenta.

The three-dimensional computational domain is given by $\Omega=[-15;15]^{3}$
covered by a $80\times 80\times 80$ uniform grid, while outgoing boundary
conditions are imposed at the border of the computational domain. As shown by
Dumbser et al. (2024), at time $t=0$ we find it useful to filter the lapse as

$\alpha=\frac{\alpha r^{6}+\epsilon\alpha_{min}}{r^{6}+\epsilon}\,,$ (58)

where $\alpha_{min}=0.01$, $\epsilon=10^{-4}$. Moreover, during the evolution,
all the metric terms are smoothed out so as to avoid metric spikes, but rather
reaching a smooth maximum value at $\gamma_{max}\sim 25$. With these caveats,
we have solved this test using $\kappa_{1}=0.2$, $\kappa_{2}=0.2$, $c=0$,
$\mu=0.1$ and with the _gamma–driver_ activated. The results of our
calculations are reported in Fig. 14 and Fig. 15, respectively, where the
contour plots of $\gamma_{xx}$, as well as the shift vector field, are shown
at four representative times. Fig. 14 refers to the calculations performed
with the FD-WENO scheme, while Fig. 15 refers to those obtained with the AFD-
WENO scheme.

By the time $t\sim 10$, the merger of the two black holes is almost complete,
and a single stationary black hole is hence produced. The vector field of the
shift $\beta^{i}$ manifests the expected behaviour, being oriented in opposite
direction with respect to the merger.

Figure 14: Head on collision of two black holes: Profiles for the metric term
$\gamma_{xx}$ shown at various time levels $(t=0,5,10,20)$. The computations
were performed using a ninth–order accurate FD–WENO scheme. Figure 15: Head on
collision of two black holes: Profiles for the metric term $\gamma_{xx}$ shown
at various time levels $(t=0,5,10,20)$. The computations were performed using
an eighth–order accurate AFD–WENO scheme.

##  5 Conclusions

Advances in the sensitivities of gravitational wave interferometers are such
that they will make it possible to make very accurate observations. As a
result, we will soon need very high order schemes for numerical relativity
that can match those accurate observations. While second order (in space
derivatives) formulations, such as BSSNOK, are routinely adopted, the
numerical technology for seamlessly formulating them at very high orders of
spatial and temporal accuracy is lacking. This forces us to look at suitable
high order accurate discretizations of first order formulations of the
Einstein field equations. Because first order formulations have a hyperbolic
PDE structure that is common to other very well-explored PDE systems, higher
order solution methods have been very well-developed for this class of
problem. Therefore, it is very attractive to look at high accuracy methods for
solving first order hyperbolic PDE systems. The well-developed higher order
metholologies include DG schemes, finite volume WENO schemes and finite
difference WENO schemes. All three of these classes of schemes present natural
pathways to high order of accuracy in space. However, the first order
formulations for numerical GR have the additional complication that they
result in extremely large non-conservative hyperbolic PDE systems. These
systems are indeed so large that the memory requirements for DG, and perhaps
even finite volume WENO schemes, are quite prohibitive. This forces us to take
a hard look at recent advances in finite difference WENO schemes.

The Z4 formulation of the Einstein equations, written as a first–order system,
naturally gives rise to a hyperbolic PDE system, that we have indicated with
the acronym FO-Z4 and that is briefly described in Section 2. While the
original Z4 system was proposed in first–order _conservative_ form Bona et al.
(2004a); Bona & Palenzuela-Luque (2005), a close analysis performed by Dumbser
et al. (2024) has shown that its hyperbolicity is better preserved when it is
cast as a hyperbolic PDE system in _non-conservation form_. This is of course
problematic from the numerical viewpoint, because finite difference WENO
methods for systems that are in non-conservation form are not well known to
the astrophysical community and general-purpose solution methods for PDEs with
non-conservative products have only been recently discovered. It is,
therefore, of great value to the astrophysical community that these methods
should all be concatenated in one place and shown to work for the FO-Z4 system
of general relativity. While we focus on the FO-Z4 system, the methods that we
present here are much more broadly useful for any first–order formulation of
GR. Indeed, the Introduction documents other first–order formulations for GR
that also result in hyperbolic PDEs with non-conservative products.

Section 3 shows us the different finite difference WENO technologies that are
at our disposal. These bifurcate broadly into FD-WENO schemes that are based
on higher order reconstruction and AFD-WENO schemes that are based on higher
order interpolation. FD-WENO schemes are described in Sub-sections 3.1 and
3.2. The FD-WENO scheme in Sub-section 3.2 is particularly easy to implement,
very low cost, and especially useful when the relativistic hydrodynamical
equations don’t need to be solved in conjunction with the Einstein field
equations. But we also realize that the relativistic hydrodynamics equations
are naturally written in conservation-law form. As a result, we will quite
frequently need the AFD-WENO formulation that is described in Sub-sections 3.3
and 3.4. While AFD-WENO schemes cost marginally more than their FD-WENO
cousins, they have the great advantage that they can retrieve a conservation
form when conservation needs to be respected in the physical problem. We also
present discussions on reconstruction, interpolation, non-linear hybridization
and upwinding which are all the algorithmic niceties that a robust high order
scheme must have if it is to operate stably. Visual descriptions of the
algorithms have also been provided in Figs. 1, 2 and 3 in order to make the
methods accessible to the greater community.

It is also worth pointing out that these days one desires not just a baseline
scheme but also an ecosystem of supporting ideas around the baseline scheme.
Examples of this ecosystem of ideas include well-balancing, physical
constraint preservation, and divergence- and curl-free evolution of vector
fields. Wherever possible, we have documented the additional availability of
this ecosystem of supporting capabilities for the FD-WENO and AFD-WENO schemes
described here.

We have therefore implemented the FD-WENO and AFD-WENO schemes in a three
dimensional code that uses Cartesian coordinates within the standard $3+1$
foliation of spacetime. We have successfully reproduced some of the most
relevant standard tests of NR in vacuum, as outlined in Alcubierre et al.
(2004), which include a gauge wave, the robust stability test, the Gowdy wave,
stationary isolated black holes (with excision) and a simple head-on collision
of two puncture black holes (Ansorg et al., 2004). Concerning the gauge wave,
we could evolve it on a long time scale with no deviations from the analytic
solution, provided the order of the scheme is at least 4. This points once
more to the conclusion that high order of accuracy might become absolutely
crucial in NR. Concerning single stationary black holes, we could evolve them
stably over long timescales, thus confirming the positive results obtained by
Dumbser et al. (2024) who also adopted the same FO-Z4 formulation, but with an
entirely different numerical scheme. In this test case, just like in Dumbser
et al. (2024), we had to augment our numerical scheme with a suitable _well-
balancing_ feature, in such a way to preserve perfect stationarity. We recall
that, on the contrary, the standard BSSNOK formulation has known problems in
performing such a test. Academic as it may seem at first sight, this test is
nevertheless very important since it opens the door to a rather accurate
numerical study of the normal modes of oscillations of black holes (Berti et
al., 2009; Mamani et al., 2022). Finally, concerning binary black holes
mergers, we confirm the ability of the FO-Z4 formulation in treating simple
set-ups with two moving punctures. However, the lack of the conformal factor,
that is typical of this formulation, must be compensated by appropriate
filtering of the metric spikes that might be produced during the merger.
Further investigations are required to make this approach more efficient while
at the same time physically neutral.

##  6 Acknowledgments

This work was financially supported by the Italian Ministry of Education,
University and Research (MIUR) in the framework of the PRIN 2022 project High
order structure-preserving semi-implicit schemes for hyperbolic equations and
via the Departments of Excellence Initiative 2018–2027 attributed to DICAM of
the University of Trento (grant L. 232/2016). DSB acknowledges support via NSF
grant NSF-AST-2009776, NASA grant NASA-2020-1241 and NASA grant 80NSSC22K0628.

We acknowledge the CINECA award under the ISCRA initiative, for the
availability of high-performance computing resources and support.

## References

  * Abbott et al. (2019) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019, Phys. Rev. X, 9, 031040, doi: 10.1103/PhysRevX.9.031040
  * Alcubierre (2008) Alcubierre, M. 2008, Introduction to 3+1 numerical relativity, Vol. 140 (Oxford University Press)
  * Alcubierre et al. (2004) Alcubierre, M., Allen, G., & Bona et al., C. 2004, Classical and Quantum Gravity, 21, 589, doi: 10.1088/0264-9381/21/2/019
  * Alcubierre & Brügmann (2001) Alcubierre, M., & Brügmann, B. 2001, Phys. Rev. D, 63, 104006, doi: 10.1103/PhysRevD.63.104006
  * Alekseenko & Arnold (2003) Alekseenko, A. M., & Arnold, D. N. 2003, Phys. Rev. D, 68, 064013, doi: 10.1103/PhysRevD.68.064013
  * Alic et al. (2012) Alic, D., Bona-Casas, C., Bona, C., Rezzolla, L., & Palenzuela, C. 2012, Physical Review D, 85, 064040
  * Alic et al. (2013) Alic, D., Kastaun, W., & Rezzolla, L. 2013, Physical Review D, 88, 064049
  * Alvi (2002) Alvi, K. 2002, Classical and Quantum Gravity, 19, 5153, doi: 10.1088/0264-9381/19/20/309
  * Anderson & York (1999) Anderson, A., & York, J. W. 1999, Physical Review Letters, 82, 4384, doi: 10.1103/PhysRevLett.82.4384
  * Ansorg et al. (2004) Ansorg, M., Brügmann, B., & Tichy, W. 2004, Phys. Rev. D, 70, 064011
  * Audusse et al. (2004) Audusse, E., Bouchut, F., Bristeau, M., Klein, R., & Perthame, B. 2004, SIAM Journal on Scientific Computing, 25, 2050
  * Babiuc et al. (2008) Babiuc, M. C., Husa, S., Alic, D., et al. 2008, Classical and Quantum Gravity, 25, 125012, doi: 10.1088/0264-9381/25/12/125012
  * Balsara (2010) Balsara, D. S. 2010, Journal of Computational Physics, 229, 1970
  * Balsara (2012) —. 2012, Journal of Computational Physics, 231, 7504
  * Balsara (2014) —. 2014, Journal of Computational Physics, 277, 163
  * Balsara et al. (2023a) Balsara, D. S., Bhoriya, D., Shu, C.-W., & Kumar, H. 2023a, Communications on Applied Mathematics and Computation, doi: 10.1007/s42967-023-00275-9
  * Balsara et al. (2024a) —. 2024a, Communications on Applied Mathematics and Computation (To be appear)
  * Balsara et al. (2024b) —. 2024b, Communications on Applied Mathematics and Computation (To be appear)
  * Balsara et al. (2016) Balsara, D. S., Garain, S., & Shu, C.-W. 2016, Journal of Computational Physics, 326, 780
  * Balsara & Käppeli (2022) Balsara, D. S., & Käppeli, R. 2022, Communications on Applied Mathematics and Computation, 4, 945
  * Balsara et al. (2021) Balsara, D. S., Käppeli, R., Boscheri, W., & Dumbser, M. 2021, Communications on Applied Mathematics and Computation, 1
  * Balsara et al. (2023b) Balsara, D. S., Samantaray, S., & Subramanian, S. 2023b, Communications on Applied Mathematics and Computation, 5, 428
  * Balsara & Shu (2000) Balsara, D. S., & Shu, C.-W. 2000, Journal of Computational Physics, 160, 405
  * Balsara & Spicer (1999) Balsara, D. S., & Spicer, D. S. 1999, Journal of Computational Physics, 149, 270
  * Baumgarte & Shapiro (2010) Baumgarte, T. W., & Shapiro, S. L. 2010, Numerical relativity: solving Einstein’s equations on the computer (Cambridge University Press)
  * Berberich et al. (2021) Berberich, J. P., Chandrashekar, P., & Klingenberg, C. 2021, Computers & Fluids, 104858
  * Bermudez & Vázquez-Cendón (1994) Bermudez, A., & Vázquez-Cendón, M. 1994, Computers & Fluids, 23, 1049
  * Berti et al. (2009) Berti, E., Cardoso, V., & Starinets, A. O. 2009, Classical and Quantum Gravity, 26, 163001, doi: 10.1088/0264-9381/26/16/163001
  * Beyer & Sarbach (2004) Beyer, H., & Sarbach, O. 2004, Phys. Rev. D, 70, 104004, doi: 10.1103/PhysRevD.70.104004
  * Bhoriya et al. (2024) Bhoriya, D., Balsara, D. S., & Shu, C.-W. 2024, Submitted to Journal of Scientific Computing
  * Bona et al. (2003) Bona, C., Ledvinka, T., Palenzuela, C., & Zácek, M. 2003, Phys. Rev. D, 67, 104005
  * Bona et al. (2004a) —. 2004a, Phys. Rev. D, 69, 64036
  * Bona et al. (2004b) —. 2004b, Phys. Rev. D, 69, 64036
  * Bona & Palenzuela-Luque (2005) Bona, C., & Palenzuela-Luque, C. 2005, Elements of Numerical Relativity (Berlin: Springer-Verlag), doi: 10.1007/b135928
  * Borges et al. (2008) Borges, R., Carmona, M., Costa, B., & Don, W. S. 2008, Journal of computational physics, 227, 3191
  * Botta et al. (2004) Botta, N., Klein, R., Langenberg, S., & Lützenkirchen, S. 2004, Journal of Computational Physics, 196, 539
  * Brown et al. (2012) Brown, J. D., Diener, P., Field, S. E., et al. 2012, Physical Review D - Particles, Fields, Gravitation and Cosmology, 85
  * Brügmann et al. (2008) Brügmann, B., González, J. A., Hannam, M., et al. 2008, Phys. Rev. D, 77, 024027, doi: 10.1103/PhysRevD.77.024027
  * Brügmann (2013) Brügmann, B. 2013, Journal of Computational Physics, 235, 216, doi: https://doi.org/10.1016/j.jcp.2012.11.007
  * Buchman & Bardeen (2003) Buchman, L. T., & Bardeen, J. M. 2003, Physical Review D, 67, 084017, doi: 10.1103/PhysRevD.67.084017
  * Buchman & Bardeen (2005) Buchman, L. T., & Bardeen, J. M. 2005, Phys. Rev. D, 72, 124014, doi: 10.1103/PhysRevD.72.124014
  * Cai & Ladeinde (2008) Cai, X., & Ladeinde, F. 2008, in 46th AIAA Aerospace Sciences Meeting and Exhibit, 36
  * Calabrese et al. (2006) Calabrese, G., Hinder, I., & Husa, S. 2006, Journal of Computational Physics, 218, 607, doi: https://doi.org/10.1016/j.jcp.2006.02.027
  * Castro et al. (2017) Castro, M. J., de Luna, T. M., & Parés, C. 2017, in Handbook of Numerical Analysis, Vol. 18 (Elsevier), 131–175
  * Castro & Parés (2020) Castro, M. J., & Parés, C. 2020, Journal of Scientific Computing, 82, 1
  * Clough et al. (2015) Clough, K., Figueras, P., Finkel, H., et al. 2015, Classical and Quantum Gravity, 32, 245011, doi: 10.1088/0264-9381/32/24/245011
  * Colella & Woodward (1984) Colella, P., & Woodward, P. R. 1984, Journal of computational physics, 54, 174
  * Cravero & Semplice (2016) Cravero, I., & Semplice, M. 2016, Journal of Scientific Computing, 67, 1219
  * de Felice & Clarke (1990) de Felice, F., & Clarke, C. J. S. 1990, Relativity on curved manifolds (Cambridge University Press)
  * Dedner et al. (2002) Dedner, A., Kemm, F., Kröner, D., et al. 2002, Journal of Computational Physics, 175, 645, doi: 10.1006/jcph.2001.6961
  * Del Zanna et al. (2007) Del Zanna, L., Zanotti, O., Bucciantini, N., & Londrillo, P. 2007, Astronomy & Astrophysics, 473, 11
  * Dumbser & Balsara (2016) Dumbser, M., & Balsara, D. S. 2016, Journal of Computational Physics, 304, 275
  * Dumbser et al. (2020) Dumbser, M., Fambri, F., Gaburro, E., & Reinarz, A. 2020, Journal of Computational Physics, 404, 109088
  * Dumbser et al. (2018) Dumbser, M., Guercilena, F., Köppel, S., Rezzolla, L., & Zanotti, O. 2018, Physical Review D, 97, 084053
  * Dumbser et al. (2024) Dumbser, M., Zanotti, O., Gaburro, E., & Peshkov, I. 2024, J. Comput. Phys., 504, 112875, doi: 10.1016/j.jcp.2024.112875
  * Dumbser et al. (2013) Dumbser, M., Zanotti, O., Hidalgo, A., & Balsara, D. 2013, J. Comput. Phys., 248, 257
  * Faber et al. (2007) Faber, J. A., Baumgarte, T. W., Etienne, Z. B., Shapiro, S. L., & Taniguchi, K. 2007, Phys. Rev. D, 76, 104021, doi: 10.1103/PhysRevD.76.104021
  * Fischer & Marsden (1972) Fischer, A. E., & Marsden, J. E. 1972, Communications in Mathematical Physics, 28, 1, doi: 10.1007/BF02099369
  * Gaburro et al. (2021) Gaburro, E., Castro, M. J., & Dumbser, M. 2021, SIAM Journal on Scientific Computing, 43, B1226, doi: 10.1137/21M1399154
  * Gosse (2001) Gosse, L. 2001, Mathematical Models and Methods in Applied Sciences, 11, 339
  * Gourgoulhon (2012) Gourgoulhon, E. 2012, 3+1 Formalism in General Relativity, Vol. 846, doi: 10.1007/978-3-642-24525-1
  * Gowdy (1971) Gowdy, R. H. 1971, Physical Review Letters, 27, 826, doi: 10.1103/PhysRevLett.27.826
  * Gundlach et al. (2005) Gundlach, C., Martin-Garcia, J. M., Calabrese, G., & Hinder, I. 2005, Class. Quantum Grav., 22, 3767, doi: 10.1088/0264-9381/22/17/025
  * Harten et al. (1986) Harten, A., Osher, S., Engquist, B., & Chakravarthy, S. R. 1986, Applied Numerical Mathematics, 2, 347
  * Henrick et al. (2006) Henrick, A. K., Aslam, T. D., & Powers, J. M. 2006, Journal of Computational Physics, 213, 311
  * Hilditch et al. (2012) Hilditch, D., Bernuzzi, S., Thierfelder, M., et al. 2012, arXiv:1212.2901.  https://arxiv.org/abs/1212.2901
  * Hilditch & Richter (2012) Hilditch, D., & Richter, R. 2012, Physical Review D, 86, 123017, doi: 10.1103/PhysRevD.86.123017
  * Hu et al. (2013) Hu, X. Y., Adams, N. A., & Shu, C.-W. 2013, Journal of Computational Physics, 242, 169
  * Inghirami et al. (2016) Inghirami, G., Del Zanna, L., Beraudo, A., et al. 2016, European Physical Journal C, 76, 659, doi: 10.1140/epjc/s10052-016-4516-8
  * Jiang & Shu (1996) Jiang, G.-S., & Shu, C.-W. 1996, Journal of computational physics, 126, 202
  * Jiang et al. (2013) Jiang, Y., Shu, C.-W., & Zhang, M. 2013, SIAM Journal on Scientific Computing, 35, A1137, doi: 10.1137/120889885
  * Jiang et al. (2014) —. 2014, Methods and applications of analysis, 21, 1
  * Käppeli (2022) Käppeli, R. 2022, Living Reviews in Computational Astrophysics, 8, 2
  * Kidder et al. (2001) Kidder, L. E., Scheel, M. A., & Teukolsky, S. A. 2001, Phys. Rev. D, 64, 064017, doi: 10.1103/PhysRevD.64.064017
  * Kiuchi et al. (2017) Kiuchi, K., Kawaguchi, K., Kyutoku, K., et al. 2017, Phys. Rev. D, 96, 084060, doi: 10.1103/PhysRevD.96.084060
  * LeVeque (1998) LeVeque, R. 1998, Journal of Computational Physics, 146, 346
  * Levy et al. (2000) Levy, D., Puppo, G., & Russo, G. 2000, Applied Numerical Mathematics, 33, 415 – 421, doi: 10.1016/S0168-9274(99)00108-7
  * Lousto & Healy (2023) Lousto, C. O., & Healy, J. 2023, Classical and Quantum Gravity, 40, 09LT01, doi: 10.1088/1361-6382/acc7ef
  * Lück et al. (2022) Lück, H., Smith, J., & Punturo, M. 2022, Third-Generation Gravitational-Wave Observatories (Singapore: Springer Nature Singapore), 283–300, doi: 10.1007/978-981-16-4306-4_7
  * Löffler et al. (2012) Löffler, F., Faber, J., Bentivegna, E., et al. 2012, Classical and Quantum Gravity, 29, 115001, doi: 10.1088/0264-9381/29/11/115001
  * Mamani et al. (2022) Mamani, L. A. H., Masa, A. D. D., Sanches, L. T., & Zanchin, V. T. 2022, European Physical Journal C, 82, 897, doi: 10.1140/epjc/s10052-022-10865-1
  * McLachlan (2008) McLachlan. 2008, McLachlan, a Public BSSN Code.  http://www.cct.lsu.edu/~eschnett/McLachlan/
  * Merriman (2003) Merriman, B. 2003, Journal of Scientific Computing, 19, 309
  * Munz et al. (2000) Munz, C., Omnes, P., Schneider, R., Sonnendrücker, E., & Voss, U. 2000, J. Comput. Phys., 161, 484
  * Nonomura et al. (2010) Nonomura, T., Iizuka, N., & Fujii, K. 2010, Computers & Fluids, 39, 197
  * Olivares et al. (2022) Olivares, H., Peshkov, I. M., Most, E. R., Guercilena, F. M., & Papenfort, L. J. 2022, Phys. Rev. D, 105, 124038, doi: 10.1103/PhysRevD.105.124038
  * Pareschi & Rey (2017) Pareschi, L., & Rey, T. 2017, Computers & Fluids, 156, 329, doi: https://doi.org/10.1016/j.compfluid.2017.07.013
  * Peterson et al. (2023) Peterson, A. J., Willcox, D., & Mösta, P. 2023, Classical and Quantum Gravity, 40, 245013, doi: 10.1088/1361-6382/ad0b37
  * Rezzolla & Zanotti (2013) Rezzolla, L., & Zanotti, O. 2013, Relativistic hydrodynamics (Oxford University Press)
  * Shibata & Nakamura (1995) Shibata, M., & Nakamura, T. 1995, Physical Review D, 52, 5428
  * Shu & Osher (1988) Shu, C.-W., & Osher, S. 1988, Journal of computational physics, 77, 439
  * Shu & Osher (1989) —. 1989, Journal of Computational Physics, 83, 32
  * Spiteri & Ruuth (2002) Spiteri, R. J., & Ruuth, S. J. 2002, SIAM Journal on Numerical Analysis, 40, 469
  * Spiteri & Ruuth (2003) —. 2003, Mathematics and Computers in Simulation, 62, 125
  * Thierfelder et al. (2011) Thierfelder, M., Bernuzzi, S., & Brügmann, B. 2011, Phys. Rev. D, 84, 044012, doi: 10.1103/PhysRevD.84.044012
  * Tichy et al. (2023) Tichy, W., Ji, L., Adhikari, A., Rashti, A., & Pirog, M. 2023, Classical and Quantum Gravity, 40, 025004, doi: 10.1088/1361-6382/acaae7
  * Wiltshire et al. (2009) Wiltshire, D. L., Visser, M., & Scott, S. M. 2009, The Kerr Spacetime
  * Wu & Tang (2015) Wu, K., & Tang, H. 2015, Journal of Computational Physics, 298, 539, doi: 10.1016/j.jcp.2015.06.012
  * Xing & Shu (2011) Xing, Y., & Shu, C.-W. 2011, Advances in Water Resources, 34, 1026
  * Xu & Shu (2024) Xu, Z., & Shu, C.-W. 2024
  * Yamamoto et al. (2008) Yamamoto, T., Shibata, M., & Taniguchi, K. 2008, Phys. Rev. D, 78, 064054, doi: 10.1103/PhysRevD.78.064054
  * Zheng et al. (2021) Zheng, F., Shu, C.-W., & Qiu, J. 2021, Journal of Computational Physics, 445, 110597
  * Zhu & Qiu (2016) Zhu, J., & Qiu, J. 2016, Journal of Computational Physics, 318, 110
  * Zhu & Shu (2018) Zhu, J., & Shu, C.-W. 2018, Journal of Computational Physics, 375, 659
  * Zlochower et al. (2005) Zlochower, Y., Baker, J. G., Campanelli, M., & Lousto, C. O. 2005, Phys. Rev. D, 72, 024021, doi: 10.1103/PhysRevD.72.024021