# Numerical relativity simulations of compact binaries: comparison of cell-
and vertex-centered adaptive meshes

Boris Daszuta1 William Cook1 Peter Hammond2 Jacob Fields2,3 Eduardo M.
Gutiérrez2,3 Sebastiano Bernuzzi1 David Radice2,3,4 1Theoretisch-
Physikalisches Institut, Friedrich-Schiller-Universität Jena, 07743, Jena,
Germany  2Institute for Gravitation and the Cosmos, The Pennsylvania State
University, University Park, PA 16802, USA  3Department of Physics, The
Pennsylvania State University, University Park, PA 16802, USA  4Department of
Astronomy & Astrophysics, The Pennsylvania State University, University Park,
PA 16802, USA

###### Abstract

Given the compact binary evolution problem of numerical relativity, in the
finite-difference, block-based, adaptive mesh refinement context, choices must
be made on how evolved fields are to be discretized. In GR-Athena++, the
space-time solver was previously fixed to be vertex-centered. Here, our recent
extensions to a cell-centered treatment, are described. Simplifications in the
handling of variables during the treatment of general relativistic magneto-
hydrodynamical (GRMHD) evolution are found. A novelty is that performance
comparison for the two choices of grid sampling is made within a single code-
base. In the case of a binary black hole inspiral-merger problem, by evolving
geometric fields on vertex-centers, an average $\sim 20\%$ speed increase is
observed, when compared against cell-centered sampling. The opposite occurs in
the GRMHD setting. A binary neutron star inspiral-merger-collapse problem,
representative of typical production simulations is considered. We find that
cell-centered sampling for the space-time solver improves performance, by a
similar factor.

###### pacs:

04.25.D-, 04.30.Db, 95.30.Sf, 95.30.Lz, 97.60.Jd 97.60.Lf

##  I Introduction

Fashioning accurate description of astrophysically sourced, binary constituent
encounters, under extreme conditions, is a pressing challenge for the
numerical relativity (NR) community. Multimessenger detection efforts Abbott
_et al._ (2017a); Goldstein  _et al._ (2017); Savchenko  _et al._ (2017)
benefit when supplemented by such data. Observational efforts of next-
generation gravitational-wave detectors Akutsu  _et al._ (2020); Amaro-Seoane
(2017); Punturo  _et al._ (2010); Abbott  _et al._ (2017b) will also greatly
benefit from higher fidelity reference numerical data, that spans a larger
range of parameter space (or e.g. models Nagar  _et al._ (2018) suitably
informed by). To confront simulation challenges, sophisticated general
relativistic, magneto-hydrodynamical (GRMHD) techniques (see e.g. Font
(2007)), and microphysics descriptions Shibata  _et al._ (2011); Radice  _et
al._ (2022) have been pursued. In tandem, the development of code
infrastructure, and algorithms, with the aim of improved performance scaling
on modern high-performance computing (HPC) resources, has also received
recent, particular attention. Some examples of modern NR codes under active
development with these goals include: bamps Hilditch  _et al._ (2016); Bugner
_et al._ (2016), Dendro-GR Fernando  _et al._ (2018), Einstein Toolkit Löffler
_et al._ (2012) (GRaM-X Shankar  _et al._ (2023)), GRChombo Clough  _et al._
(2015), NMesh Tichy  _et al._ (2023), SpECTRE Kidder  _et al._ (2017), and
SPHINCS_BSSN Rosswog and Diener (2021). A variety of numerical methods
Grandclement and Novak (2009); Doulis  _et al._ (2022); Alfieri  _et al._
(2018), algorithmic approaches Stout  _et al._ (1997); Berger and Oliger
(1984); Burstedde  _et al._ (2019); Morton (1966), and NR formulations
Nakamura  _et al._ (1987); Shibata and Nakamura (1995); Baumgarte and Shapiro
(1999); Bernuzzi and Hilditch (2010); Hilditch  _et al._ (2013); Friedrich
(1985); Pretorius (2005); Lindblom  _et al._ (2006) are employed, to list but
a few prominent examples. In spite of this, the actual time taken (HPC
resource consumed) in the calculation of numerical evolutions, particularly
for extreme scenarios that require highly resolved physical scales, can
require weeks, or months Lousto and Healy (2020). Clearly, even minor
performance optimizations, of a few percent, can have a significant impact on
such HPC wait times, together with broader, sustainability considerations
Lannelongue  _et al._ (2023).

In this work, we focus on finite-difference based evolution of the
$\mathrm{Z}4$c formulation Bernuzzi and Hilditch (2010); Hilditch  _et al._
(2013) of NR, as coupled to the general relativistic, hydrodynamical (GRHD)
formulation111For simplicity we do not investigate magnetic fields in this
work. of Banyuls  _et al._ (1997). This is precisely the setting treated by
our code GR-Athena++ Daszuta  _et al._ (2021); Cook  _et al._ (2023) (see also
Daszuta and Cook (2024)). In brief, this is our extension to the astrophysical
(radiation), GRMHD code Athena++ White  _et al._ (2016); Felker and Stone
(2018); Stone  _et al._ (2020), we have geared toward evolution with dynamical
space-time. Building upon the underlying octree, block-based adaptive-mesh-
refinement (AMR) infrastructure, with addition of our $\mathrm{Z}4$c based
space-time solver both in vacuum Daszuta  _et al._ (2021), and for GRMHD with
dynamical space-time Cook  _et al._ (2023), has led to demonstrable scaling
efficiencies in excess of $80\%$ on (pre-)exascale HPC infrastructure.

Specifically, here, we aim to address a seemingly innocuous question: should
we select the geometric fields, that is, the dynamical fields of
$\mathrm{Z}4$c, to be discretized on vertex-centers (VC) or cell-centers? Our
initial choice was based on considerations for binary black hole evolution,
where transfer of data between differing levels of refinement, is less
computationally expensive for VC in contrast to cell-centers. On the other
hand, in consideration of GR(M)HD, typically hydrodynamical sampling is on
cell-centers (see e.g. BAM Brügmann  _et al._ (2008), GRHydro Mösta  _et al._
(2014), and WhiskyTHC Radice  _et al._ (2014)). This choice is in part due to
greater amenability of the underlying hydrodynamical treatment as based on
conservative, high-resolution-shock-capture (HRSC) methods (see e.g.
Thierfelder  _et al._ (2011a)). Athena++ also handles hydrodynamical
description on cell-centers. Having previously fixed geometric fields to VC,
discretized fields must consequently be transferred between grids, which
incurs an overhead. A priori it is not necessarily clear as to how large this
is. Our aim here is to investigate switching to a cell-centered treatment of
the space-time solver in GR-Athena++, thus providing an answer. Such a code-
internal, investigation, does not appear to have been previously performed.
This gap we fill here. This paper proceeds as follows. In §II we briefly
explain the computational domain structure and the variable discretization
strategy. Subsequently, we present a simplified refinement strategy, tailored
to the binary merger problem, and independent of geometric field sampling
selection. In §III we demonstrate the robustness of the technique through
numerically solving BBH and BNS inspiral-merger problems, representative of
typical production simulations. Section §IV concludes.

##  II Method

GR-Athena++ Daszuta  _et al._ (2021); Cook  _et al._ (2023) builds upon
Athena++ White  _et al._ (2016); Felker and Stone (2018); Stone  _et al._
(2020) and thus inherits overall features of the latter framework. In brief,
problems are formulated over a target computational domain $\Omega$ (the Mesh)
which is partitioned into a collection of sub-domains $\Omega_{i}$ (MeshBlock
objects). We have $\Omega=\sqcup_{i\in Z}\Omega_{i}$ where $Z$ may be thought
of as an ordered parametrization of multi-dimensional coordinates to a space-
filling curve Morton (1966); Burstedde  _et al._ (2019). For a $d$-dimensional
problem the overall extent of $\Omega$ together with suitable boundary data on
$\partial\Omega$ must be provided. Following this, the underlying
discretization is fixed through specification of Mesh and MeshBlock sampling
$N_{M}=(N_{M_{1}},\,\dots,\,N_{M_{d}})$ and
$N_{B}=(N_{B_{1}},\,\dots,\,N_{B_{d}})$, respectively, where $N_{B}$ must
divide $N_{M}$ component-wise222For a uniform number of samples in each
dimension we typically use a single scalar to represent all component values
of the tuples.. In the case that (local) refinement is desired, MeshBlock
partitioning into $2^{d}$ new MeshBlock objects is recursively performed. This
involves halving size along each axis at a fixed number of samples, which
doubles the resolution at each step. The aforementioned is done while
maintaining an (at-most) $2:1$ refinement ratio between neighboring
$\\{\Omega_{i}\\}_{i\in Z}$. While the logical-rectilinear structure here
described may then be further mapped to suitable curvilinear coordinates, for
simplicity we focus on the Cartesian context.

###  II.1 Cell-centered extension

In GR-Athena++ Daszuta  _et al._ (2021); Cook  _et al._ (2023) we extended
sampling support to enable vertex-centered (VC) discretization of field
variables together with infrastructure enabling high-order approximations to
be carried out. In this work we focus on analogous extensions in description
of geometric quantities on cell-centers (CX). A variable $\mathcal{V}(x)$
where $x\in[a,\,b]$ $(a<b)$ is CX discretized333The grid sampling will also be
referred to as CC. On CX fields are to be thought of as pointwise sampled,
whereas on CC they carry a cell-averaged interpretation. to
$\mathcal{V}_{i+1/2}$ over $x_{i+1/2}:=a+(i+1/2)\delta x$ where $\delta
x{}:=(b-a)/N_{B}$ and $i=0,\,\ldots,\,N_{B}-1$. This yields $N_{B}$ total
samples. In practice the interval is extended by an additional
$\mathcal{N}_{\mathrm{g}}{}$ “ghost” nodes either-side444In the case that
$d=\dim(\Omega)>1$ we utilize a tensor product grid of such extended one-
dimensional discretizations.. This extension facilitates evaluation of e.g.
derivative approximants, transfer of field data between neighboring MeshBlock
objects that may differ in the level of refinement, and/or imposition of
boundary conditions.

In order to populate ghost nodes, communication is required. In the case that
$\mathcal{V}$ is discretized over a Mesh decomposed at a fixed, common level,
under CX sampling, this entails that nearest-neighbor MeshBlock objects can be
utilized to populate ghost-layers through a straight-forward copy. Suppose
instead that sub-domains differ in level $l$. In this situation, analogous to
the VC case, each sub-domain ${}^{l}\Omega_{i}$ (where prefix $l$ indicates
level) as described by a MeshBlock, features a fundamental and an additional
coarse representation of the sampled variable. This latter is one level above
at $l-1$. Thus if $N_{B}$ is the number of physical points on a MeshBlock
extended by $2\mathcal{N}_{\mathrm{g}}{}$ ghost nodes, then the complementary
coarse representation features $N_{B}/2$ physical points extended by
$2\mathcal{N}_{\mathrm{cg}}{}$ coarse ghosts. The situation is illustrated
(locally) in Fig.1.

Figure 1: One-dimensional illustration of inter-MeshBlock communication given
a refined local Mesh structure. Each MeshBlock features a fundamental
representation of a sampled field depicted through square nodes which is
extended in each direction through ghost nodes (circles) that facilitate
communication. Additionally a complementary coarse sampling is present which
is explicitly depicted for data on level $l$ and suppressed on level $l-1$.
Here $N_{B}=12$, $\mathcal{N}_{\mathrm{g}}{}=2$, and
$\mathcal{N}_{\mathrm{cg}}{}=3$. Prior to all communication data in filled
square nodes is available; the task is thus to fill data at locations
indicated by shaded circles and diamonds. Nodes are assigned MeshBlock-local
labels for ease of discussion.

The general communication strategy and ordering of operations closely follows
and extends that of Athena++. This minimizes the amount of infrastructure
changes required. Synchronization of data across levels may be understood and
implemented by considering the transfer operations: fine to coarse
(restriction $\mathcal{R}$), and coarse to fine (prolongation $\mathcal{P}$).
As an example consider synchronization of ghost-layer data between
${}^{l-1}\Omega_{A}$ and ${}^{l}\Omega_{B}$ as depicted in Fig.1. Initially
data on the fundamental grid of ${}^{l}\Omega_{B}$ over the nodes
$(\texttt{is},\,\texttt{is}+1,\,\ldots,\,\texttt{ie})$ is known. This is first
restricted to a subset of the complementary coarse grid of ${}^{l}\Omega_{B}$
over the nodes
$(\texttt{cime}+1,\,\ldots,\,\texttt{cime}+\mathcal{N}_{\mathrm{g}}{})$; see
green diamonds in the figure. Subsequently data at the nodes
$(\texttt{ie}+1,\,\ldots,\,\texttt{ipe})$ of ${}^{l-1}\Omega_{A}$ may be
filled through communication or a copy (green circles in figure). Conversely
using data from ${}^{l-1}\Omega_{A}$ at the nodes
$(\texttt{ie}-\mathcal{N}_{\mathrm{cg}}{}+1,\,\ldots,\,\texttt{ie})$ allows
for filling of $(\texttt{cims},\,\ldots,\,\texttt{cime})$ on
${}^{l}\Omega_{B}$ (purple circles) through communication or a copy. Finally,
ghost layer nodes $(\texttt{ims},\,\ldots,\,\texttt{is}-1)$ of the fundamental
grid over ${}^{l}\Omega_{B}$ may be filled (blue circles) through
$\mathcal{P}$ based on data local to the MeshBlock. The procedure is analogous
for the pair ${}^{l}\Omega_{C}$ and ${}^{l-1}\Omega_{D}$.

Comparing node positions in Fig.1 across differing levels of refinement
clearly shows that nodes do not align but are fully staggered555One could
instead consider $3:1$ refinement ratios to obviate this.. This is one aspect
where CX differs from VC sampling Daszuta  _et al._ (2021). The coarsening and
refining of sampled field data via $\mathcal{R}$ and $\mathcal{P}$
respectively, is based on polynomial interpolation under CX and VC. Due to the
inter-level node staggering occuring in CX, a non-trivial evaluation of
interpolants for every target node is required. This is in contrast to the
case of VC where a subset of operations reduces to injection (copy) Daszuta
_et al._ (2021). The distinction heralds a performance implication as will be
demonstrated in §III.1 and §III.2.

We close this subsection with a collection of technical points concerning
implementation aspects.

  * •

The formal order of accuracy for the $\mathcal{R}$ and $\mathcal{P}$
operations on CX is controlled by the number of nodes selected. We fix this
choice based on $\mathcal{N}_{g}$ selection so to match the approximation
order with that of the finite-differencing.

  * •

Given a target node, source data entering an interpolant describing CX
$\mathcal{P}$ has a directional bias due to grid staggering. This can be
understood by noting that the source data nodes are not symmetrically placed
about the location of the target node, but offset. Consider at a desired
formal order of accuracy, the initial coarse, complementary restriction needed
when communicating data to a finer neighboring MeshBlock. There may be
insufficient nodes available to evaluate with symmetrically spaced source
data. In this situation, we again bias the stencil selection.

  * •

As is well-known, floating point (FP) arithmetic is not necessarily
associative: $(a+b)+c\neq a+(b+c)$ for FP numbers $\\{a,\,b,\,c\\}$. If care
is not taken in how numerical expressions are ordered, errors can e.g. lead to
catastrophic cancellations, or result in spurious symmetry breaking during
simulations Fleischmann  _et al._ (2019). A combination of Kahan-Babuska-
Neumaier compensated summation Neumaier (1974), and explicit symmetrization of
arithmetic operations involving averaging can partially mitigate this
Fleischmann  _et al._ (2019). This can be optionally activated for interpolant
evaluation. We have not however found it necessary for the simulations
presented here.

  * •

There are now two styles of cell-centered sampling where discretized field
variables may be registered and stored: that inherited from Athena++, which we
denote CC, and the newly introduced CX. The principal distinction is that for
a Mesh featuring MeshBlock objects on distinct levels, the $\mathcal{R}$ and
$\mathcal{P}$ operators are distinct: for CC there is minmod based slope-
limiting Stone  _et al._ (2020), whereas for CX we utilize the Barycentric
Lagrange approach Berrut and Trefethen (2004). Special care has been taken to
ensure that distinct fields may be suitably registered as either CC or CX with
their respective, logic available concurrently.

  * •

The VC double-restriction operation described in Daszuta  _et al._ (2021) we
have found to not be strictly required. Its removal means that when utilizing
vertex-centered sampling we now have $N_{B}\geq\max(4,\,\mathcal{N}_{g})$.
This allows for smaller MeshBlock size, which offers greater flexibility in
tailoring resolution adapted regions.

###  II.2 Refinement strategy

In order to resolve dynamical features over widely disparate spatial scales we
utilize adaptive mesh refinement (AMR). Our purpose here is to fix notation,
and to point to two simple, minimal criteria that will facilitate consistency
tests of the new infrastructure. The aforementioned feature evolution of
initial data describing a binary black hole (BBH) (§III.1), and a binary
neutron star (§III.2). For a comprehensive investigation of AMR criteria in
GR-Athena++, and effect on waveform quality for BBH evolution, see Rashti  _et
al._ (2024).

The AMR conditions of GR-Athena++ are quite flexible. A user-specified
function is defined for each MeshBlock; this function returns a flag that
controls whether to (de)-refine. Consider puncture-based evolution Brügmann
_et al._ (2008) of a BBH system. Here a black hole center is described through
a puncture. The position of a puncture $\mathbf{x}_{p}(t)$ can be computed by
augmenting the dynamical system evolved through an additional ODE based on
gauge information Campanelli  _et al._ (2006). In short, this is one variety
of a so-called tracker. This allows us to construct a dynamically adjusted,
octree-based, pseudo-box-in-box (pBIB) refined Mesh as detailed in Daszuta
_et al._ (2021) (or $L_{\infty}$-norm refinement in the language of Rashti
_et al._ (2024)).

Let us suppose the Mesh is the Cartesian coordinatized, computational domain
$\Omega:=[-D,\,D]^{3}$. A parameter controlling the maximal level of
refinement $N_{L}$ is specified. A hierarchy of MeshBlock objects is then
constructed, which increases in resolution as a puncture is approached (see
e.g. Fig.6). On the finest resolution level the grid spacing will be $\delta
x=2D/(N_{M}2^{N_{L}-1})$ Daszuta  _et al._ (2021). As an alternative to pBIB,
we can exploit information from the trackers and directly impose a desired
region to be at fixed resolution (equivalently level). Denote a spherical
region of radius $R$, centered at $p$ as $\mathbb{S}_{R}(p)$. We center such a
region on a puncture and impose a desired local level of refinement $l$. This
we write as $\mathbb{S}_{l,R}(\mathbf{x}_{p}(t))$.

In a similar vein, consider a dynamical, user selected, control field $\chi$
featuring a local extrema $\chi^{*}$. Another type of tracker is possible to
construct based on following the time development of $\chi^{*}(t)$ through
$\mathbb{S}_{l,R}(\chi^{*})$. This allows for the region surrounding the
extrema to be captured as its position changes in time. If multiple such
regions are registered, occur at differing $l$, and happen to intersect, then
we take the finest level as the target level. This is the principal condition
that we employ for AMR during BNS evolution.

##  III Validation

Our first goal in this section is validating the extended cell-centered
implementation in the context of vacuum binary black hole evolution. The
space-time solver implemented in GR-Athena++ is based on the $\mathrm{Z}4$c
formulation Bernuzzi and Hilditch (2010); Hilditch  _et al._ (2013). Here the
space-time is constructed based on the evolution of the dynamical, geometric
variables
$\mathcal{Z}:=\big{(}\chi,\,\tilde{\gamma}{}_{ij},\,\hat{K},\tilde{A}{}_{ij},\,\Theta,\,\tilde{\Gamma}{}^{i}\big{)}$.
The definition and meaning of these variables is standard and may be found,
together with implementation details, in Daszuta  _et al._ (2021). Previously,
motivated by AMR efficiency considerations in the vacuum context, we have
chosen these to evolve over VC grid sampling. We compactly denote this through
$\mathcal{Z}_{\mathrm{VC}}$. This work details a relaxation of this
restriction, extending functionality to allow for CX sampling, i.e. sampling
and evolving directly as $\mathcal{Z}_{\mathrm{CX}}$. We detail tests of this
in §III.1.

Secondly, we investigate the binary neutron star (BNS) inspiral-merger
problem. Our recently implemented GRMHD treatment in GR-Athena++ is based on
the formulation of Banyuls  _et al._ (1997). For full details, and
conventions, see Cook  _et al._ (2023). Here we recall key ingredients of
constructing the matter side of the evolution as based on selecting
$\mathcal{Z}_{\mathrm{VC}}$. For simplicity, we restrict our exposition to the
GRHD sector. For the hydrodynamics, the (conserved) field variables are
$\mathcal{Q}_{A}:=\big{(}D,\,S{}_{j},\,\tau\big{)}$, and satisfy a balance
law. Recall that this is evolved based on conservative, finite-volume, HRSC
techniques Thierfelder  _et al._ (2011a). The natural choice of sampling, that
also leverages extant Athena++ infrastructure is cell-centered, i.e. we take
$\mathcal{Q}_{A,\mathrm{CC}}$. We also work with a set of complementary,
primitive variables
$\mathcal{R}_{A}:=\big{(}\rho,\,\tilde{u}{}^{i},\,p\big{)}$. The relation
$\mathcal{Q}_{A}(\mathcal{R})$ is non-linearly implicit, and is evaluated
through numerical inversion Noble  _et al._ (2006). In this work the
conservative to primitive variable mapping is handled through the use of
RePrimAnd Kastaun  _et al._ (2021). This inversion requires geometric fields
to be available on the same grid sampling as the matter. Hence an intergrid
interpolation $\mathcal{Z}_{\mathrm{VC}}\rightarrow\mathcal{Z}_{\mathrm{CC}}$
is required.

In order to evolve the quantities $\mathcal{Q}_{A,\mathrm{CC}}$ via HRSC we
require assembly of fluxes at cell-interfaces (that is, face-centers (FC)). In
this work, this is achieved based on the reconstruction of
$\mathcal{R}_{A,\mathrm{CC}}$ to FC with the fifth order method of Borges  _et
al._ (2008) denoted WENO5Z. Flux assembly also requires a subset of
$\mathcal{Z}_{\mathrm{FC}}$ and hence another interpolation. The interface
fluxes are combined according to the Local-Lax-Friedrichs (LLF) scheme (see
e.g. Zanna and Bucciantini (2002)). Finally the sources appearing in the
hydrodynamical balance law are assembled, which again requires
$\mathcal{Z}_{\mathrm{CC}}$. On the other hand, the gravitational matter
density, momenta, and stress quantities that recouple the system to the space-
time solver, must be provided on VC. This entails the interpolation of
hydrodynamical quantities to VC. It is clear from the above summary, that by
instead taking the evolved geometric variables as $\mathcal{Z}_{\mathrm{CX}}$
from the start, fewer transfers of geometric data between grids would be
required. We further validate our code functionality, for BNS, in §III.2.

In addition to ensuring consistency of the implementation, and that problems
are solved correctly, a practical performance assessment is needed. We
demonstrate that execution speed, and suitable choice of
$\mathcal{Z}_{\mathrm{VC}}$ ($\mathcal{Z}_{\mathrm{CX}}$) is contingent on the
class of physical problem simulated.

###  III.1 Binary Black Holes

We assess the performance of our code by focusing on the well-known
calibration binary black hole (BBH) problem of Brügmann  _et al._ (2008). The
$\mathrm{Z}4$c system is coupled to the moving puncture gauge using the
$1+\log$-slicing condition Bona  _et al._ (1996) and the Gamma driver shift
Alcubierre  _et al._ (2003); van Meter  _et al._ (2006), with parameters as
described in Daszuta  _et al._ (2021). Initial data is generated through the
use of the TwoPunctures 666See https://bitbucket.org/bernuzzi/twopuncturesc/.
library based on Ansorg  _et al._ (2004). We model two initial non-spinning
black holes of bare-mass $m_{\pm}=0.483~{}M$, located on the $x-$axis, with
$x^{1}_{\mathrm{p,\pm}}(t=0)=\pm 3.257~{}M$, and initial momenta directed
along the $y-$axis, $p^{2}_{\mathrm{\pm}}(t=0)=\mp 0.133~{}M$. These parameter
choices lead to an evolution featuring $\sim 3$ orbits to merger, which occurs
at an evolution time $t\sim 170~{}M$.

In these tests the Courant-Friedrich-Lewy (CFL) is fixed at
$\mathcal{C}_{1}=0.2$. This choice determines the global time-step taken on
each MeshBlock based on the finest grid-spacing (occuring in the vicinity of a
puncture and denoted $\delta x_{p}$). Kreiss-Oliger dissipation is selected as
$\sigma_{\mathrm{D}}=0.02$ throughout. The $\mathrm{Z}4$c constraint damping
parameters are taken as $(\kappa_{1},\,\kappa_{2})=(0.02,\,0)$ (see Daszuta
_et al._ (2021)). Unless otherwise stated, in this section we take
$\mathcal{N}_{g}=\mathcal{N}_{cg}$, that is, the number of ghosts on the
fundamental grid is set equal to the number on the complementary coarse
representation (see Fig.1). We select $\mathcal{N}_{g}=4$, which induces a
formal $6^{\mathrm{th}}$ order of spatial approximation. Time-evolution is
performed using the explicit $4^{\mathrm{th}}$ order RK$4()4[2S]$ Ketcheson
(2010).

The overall grid configurations based on pBIB refinement (see §II.2) are
summarized in Tab.(1).

$N_{M}$ | $\delta x_{p}\times 10^{-2}[M]$ | $\\#\mathrm{MB}$
---|---|---
$96$ | $3.125$ | $5144$
$128$ | $2.344$ | $6336$
$160$ | $1.875$ | $18640$
$192$ | $1.563$ | $22392$
Table 1: Properties of grid configurations utilized for vacuum tests. The
total number of MeshBlock objects ($\\#\mathrm{MB}$) at the initial time is
indicated. In all cases the Mesh is the Cartesian coordinatized domain
$\Omega:=[-D,\,D]^{3}$, where $D=1536M$ is selected. This entails that grid
spacing is uniform along each spatial axis. In each case we select the maximal
refinement through $N_{L}=11$. The MeshBlock sampling is selected as
$N_{B}=16$. Resultant puncture resolutions $\delta x_{p}$ are also provided.

In principle the calibration BBH problem possesses a underlying bitant
symmetry (i.e. reflection across the $z=0$ plane to which the co-orbit and
merger is confined). This allows for an approximate halving of the total
computational resources required for a given problem. When imposing bitant
symmetry, $\\#\mathrm{MB}$ as listed in Tab.1 is roughly halved, and in order
to keep the grid spacing fixed, we halve the Mesh sampling parameter along the
$z$-axis to $N_{M}/2$. As an initial consistency check, we select the
geometric sampling as $\mathcal{Z}_{\mathrm{VC}}$ and the $N_{M}=96$ grid of
Tab.1. We perform a comparable run with $\mathcal{Z}_{\mathrm{CX}}$, for which
we inspect the puncture tracker evolution in Fig.2. Excellent agreement can be
seen.

Figure 2: Tracker evolution during calibration BBH inspiral and merger. The
one-parameter families $\mathbf{x}_{\pm}(t)$ describe where punctures are
located. Consistent behaviour for cell-centered extended (CX) and vertex-
centered (VC) choices of sampling is observed. Parameters: $N_{M}=96$, and
$N_{B}=16$.

The BBH evolution considered here should result in an evolution with a strict
symmetry about the $z=0$ plane. During numerical evolution, however, this
condition is not necessarily preserved. The $z$ component of a puncture can
potentially slowly drift out of the plane. The behaviour is demonstrated in
Fig.3 for full $3\mathrm{D}$ runs.

Figure 3: Evolution of the $z$-component of the tracker $\mathbf{x}_{-}(t)$
during calibration BBH inspiral and merger. The CX and VC sampling choices are
compared at fixed resolution. Parenthesis
$(\mathcal{N}_{g},\,\mathcal{N}_{cg})$ indicate ghost node selection. The
“bit.ID” set evolved with the choice $\mathcal{Z}_{\mathrm{CX}}$ has ID
explicitly symmetrized. All runs of this figure are performed over full
$3\mathrm{D}$ grids, regardless of ID symmetrization. Parameters: $N_{M}=96$,
and $N_{B}=16$. See text for further discussion.

In the case of $\mathcal{Z}_{\mathrm{CX}}$, based on the structure of the
level-to-level transfer operators, we do not observe, nor expect, strict FP
symmetry preservation. Instead, during the full $3\mathrm{D}$ run we find a
very minor drift developing, where the $z$-component of the puncture reaches
$\sim\mathcal{O}(10^{-10})$. We also verify in Fig.3 that $\mathcal{N}_{cg}$
may indeed be independently selected, so as to increase the formal order of
accuracy in the prolongation operation for $\mathcal{P}$. We have similarly
verified this for $\mathcal{Z}_{\mathrm{CX}}$. As this does not appear to
substantively change simulation quality for the short evolution here, we
proceed with runs that fix $\mathcal{N}_{g}=\mathcal{N}_{cg}=4$.

One reason for this behaviour, is that the TwoPunctures library does not allow
for explicit enforcement of an underlying symmetry about $z=0$. Thus, when
interpolating an initial data (ID) set to a Cartesian grid, the result
violates the bitant property at the FP level. We have taken special care to
allow for symmetrization of the prepared ID by rewriting the interpolation
process. The purpose here is to provide a test on preservation of FP symmetry
across various elements of the infrastructure, such as, e.g. operator
approximants (finite-differencing kernels, interpolation, etc.).

In the case of $\mathcal{Z}_{\mathrm{VC}}$, explicit bitant symmetrization of
the ID does lead to an evolution with exact confinement of the punctures to
the $z=0$ plane. This is in contrast to runs in Fig.3.

Having passed a tracker-based consistency check, we now check gravitational
waveforms (GW). To do this, we construct the Weyl scalar $\Psi_{4}$ Brügmann
_et al._ (2008), and then interpolate onto geodesic spheres Wang and Lee
(2011). Subsequently, we project through numerical quadrature evaluation, onto
spin-weighted spherical harmonics Goldberg  _et al._ (1967) ${}_{s}Y_{lm}$ of
spin-weight $s=-2$. Our convention differs from the reference by a Condon-
Shortley phase of $(-1)^{m}$. The above allows us to extract $\psi{}_{lm}$.
The $(\ell,\,m)$ modes of the gravitational wave strain $h$ satisfy
$\psi_{\ell m}=\ddot{h}_{\ell m}$. The strain is then given by the mode-sum:

$R\left(h_{+}-ih_{\times}\right)=\sum_{\ell=2}^{\infty}\sum_{m=-\ell}^{\ell}h_{\ell
m}(t)\;{}_{-2}Y_{\ell m}(\vartheta,\varphi)\,.$ (1)

Following the convention of the LIGO algorithms library LIGO Scientific
Collaboration (2018) we set:

$Rh{}_{\ell m}=A{}_{\ell m}\exp(-i\phi{}_{\ell m}),$ (2)

and the gravitational-wave instantaneous frequency is defined through:

$\omega{}_{\ell m}=\frac{d}{dt}\phi{}_{\ell m}.$ (3)

Note that the peak amplitude of the $(2,\,2)$ mode of $h$ is conventionally
taken as the time of merger. We denote this $u_{\mathrm{merger}}$.

The quantity $\psi{}_{22}$ is illustrated for the Calibration BBH in Fig.4.

Figure 4: Consistent behaviour between $\mathcal{Z}_{\mathrm{VC}}$ and
$\mathcal{Z}_{\mathrm{CX}}$ for calibration BBH waveforms with and without
bitant symmetry reduction. Black profiles depict the amplitude of the given
mode based on calculations involving bitant symmetry reduction. Parameters:
$N_{M}=96$, and $N_{B}=16$; extraction radius $R=85M$.

Rather than just imposing ID bitant symmetrization, symmetry reduction on the
$3\mathrm{D}$ grid is also imposed for a subset of runs. The waveforms
continue to overlap consistently. This provides another important consistency
test on implementation details.

Beyond verification of qualitative behaviour, we next focus on the
gravitational strain, impose bitant symmetry reduction, and assess convergence
properties based on Appendix A. The result, for constructing $h_{lm}$, based
on time-domain integration (see e.g Albanesi  _et al._ (2024)), is shown in
Fig.5.

Figure 5: The $(2,\,2)$ mode of the gravitational strain (top-panel) together
with convergence properties in the associated phase for CX and VC sampling in
the middle and lower panels respectively. All quantities are taken as
functions of the Schwarzschild tortoise coodinate $u$. We assess
$6^{\mathrm{th}}{}$ order compatible convergence. Indicated in a common
vertical black line across panels is $u_{\mathrm{merger}}\simeq 158\,M$. We
select calculations involving bitant symmetry reduction based on the coarse
(c), medium (m), and fine (f) resolution triplet induced through selecting
$N_{M}\in\\{128,\,160,\,192\\}$. Extraction performed at $R=85M$.

A $6^{\mathrm{th}}{}$ order compatible convergence trend is observed for the
phase $\phi_{22}$. This is present for both $\mathcal{Z}_{\mathrm{VC}}$, and
$\mathcal{Z}_{\mathrm{CX}}$. We have also verified that behaviour persists,
consistently between grid sampling choices, for e.g. the $(4,\,4)$ mode.

As described in §II.2 AMR conditions may be freely selected. In the approach
to a puncture location, the conformal factor $\chi$ approaches a local minimum
(see e.g. Thierfelder  _et al._ (2011b)). Indeed, we have numerically verified
that the trackers shown in Fig.2 can also be constructed based on such local
minimum. On account of the underlying dynamical (approximate) symmetries of
the problem, together with their impact on conserved quantities Mewes  _et
al._ (2020), and GW extraction taking place on spheres, an approximately
matched AMR condition is motivated. To this end, spherically matched Mesh
refinement is induced by fixing levels based on the collection
$\\{\mathbb{S}{}_{11,3}(\chi^{*}_{\pm}),\,\mathbb{S}{}_{6,100}(\mathbf{0})\\}$.
That is, two spherical regions of radius $3\,M$, with $N_{L}=11$ follow both
punctures. A central fixed sphere of radius $100\,M$, and $N_{L}=6$ is set to
improve resolution in the wave-zone. This is compared against pBIB in Fig.6
and Fig.7.

Figure 6: Mesh structure and values of $\tilde{\gamma}_{xx}$ under
$\mathcal{Z}_{\mathrm{CX}}$ at $t\simeq 110M$ in the orbital plane. Each
quadrant depicts distinct resolutions (AMR choices) for the same physical
calibration problem configuration. Counter-clockwise, starting from the upper
right quadrant we have pBIB AMR with $N_{M}\in\\{96,\,128,\,192\\}$. The lower
right quadrant depicts evolution with the dynamical spherical AMR criterion,
and $N_{M}=96$ fixed. Each grey box (MeshBlock) is sampled with $N_{B}=16$
points along each axis. For the cases involving the pBIB AMR criterion observe
that reflection of outward propagating dynamical spherical features in
$\tilde{\gamma}_{xx}$ is present in the vicinity of a change in refinement
level. See text for further discussion.  Figure 7: Zoom of the strong-field
region of Fig.6. Trackers are overlaid where lighter color indicates
increasing $t$. Observe that the finest resolution in both $N_{M}=96$ runs is
the same, however the region on the finest level is larger in the bottom right
(dynamical spherical region) run. Note: color scalings are distinct between
this figure and Fig.6.

Let us consider the behaviour of $\tilde{\gamma}_{xx}$ in the vicinity of a
change in refinement level as shown in Fig.6. It is observed that small,
spurious reflection occurs, from the underlying rectilinear grid structure. In
contrast, by comparing the two $N_{M}=96$ runs, we see that approximate
spherically matched AMR partially mitigates these spurious reflections. The
magnitude of the effect can also be seen to converge away when comparing
$N_{M}=192$. The change in refinement strategy also has an effect on GW
phasing, which we show in Fig.8. At early times the amplitude and phasing
error of the spherically matched AMR with $N_{M}=96$ is comparable to that of
a pBIB run at $N_{M}=160$. This can be attributed to a combination of
increased wave-zone resolution and better symmetry matching of the AMR (See
also the $L_{2}$ versus $L_{\infty}$ AMR discussion of Rashti  _et al._
(2024)).

Another potential advantage of running the spherically matched AMR is a
reduction in execution time. Running the spherically matched AMR at $N_{M}=96$
in bitant results in $\\#\mathrm{MB}=4532$ initially, where $\delta x_{p}$ is
as in Tab.1. Compared to pBIB AMR under bitant symmetry reduction at
$N_{M}=160$, it results in a computational resource reduction of $\sim
3.4\times$ when coarser time-step $\delta t$ at fixed CFL is taken into
account.

Figure 8: Comparison of differences in the amplitude $A$ of $(2,\,2)$ mode of
the gravitational strain (top-panel) and phase $\phi$ (bottom-panel), for runs
at different resolutions and refinement strategy, under
$\mathcal{Z}_{\mathrm{CX}}$. We consider pBIB AMR with $(c,\,m,\,f)$
corresponding to $N_{M}\in\\{96,\,160,\,192\\}$. A single run with spherically
matched AMR at $N_{M}=96$ is denoted with $c(\dagger)$. Vertical black line
across panels indicates $u_{\mathrm{merger}}$. Extraction performed at
$R=100M$.

Finally, we provide a rough characterization of performance differences
between $\mathcal{Z}_{\mathrm{CX}}$ and $\mathcal{Z}_{\mathrm{VC}}$ through a
short timing test. To this aim we evolve the problem for the calibration BBH
on the pBIB grid configurations of Tab.1. The initial checkpoint is taken at a
time of $t_{0}=5\,M$, where the AMR condition has refined the Mesh. Timings
are as summarized in Tab.2.

$N_{M}$ | $t_{F}(\mbox{{\tt cx}})$ | $t_{F}(\mbox{{\tt vc}})$ |  $r_{t}$ ($2\,\mathrm{dp}$) | $\\#\mathrm{MB}$ | $\\#\mathrm{MB}/\mathrm{Core}$
---|---|---|---|---|---
$96$ | $26.48$ | $35.17$ | $0.71$ | $2404$ | $\sim 3$
$128$ | $18.01$ | $21.98$ | $0.77$ | $3280$ | $\sim 4$
$160$ | $10.41$ | $10.80$ | $0.93$ | $8774$ | $\sim 11$
$192$ | $8.49$ | $9.03$ | $0.87$ | $11490$ | $\sim 15$
Table 2: One hour timing test on SuperMUC-NG of LRZ utilizing $16$ nodes. The
timing ratio $r$ is defined as $r_{t}:=(t_{F}(\mbox{{\tt
cx}})-t_{0})/(t_{F}(\mbox{{\tt vc}})-t_{0})$. The number of MeshBlock objects
is provided at $t_{0}$. Bitant symmetry reduction has been performed. The AMR
is based on the pBIB criterion.

We see that for pure vacuum, the space-time solver has a distinct advantage
when run with the choice $\mathcal{Z}_{\mathrm{VC}}$ in terms of speed.
Indeed, we observe an approximate average $\sim 20\%$ later final evolution
time, across $N_{M}$, when compared to the $\mathcal{Z}_{\mathrm{CX}}$ choice.
This we attribute principally to a difference in complexity of the
$\mathcal{R}/\mathcal{P}$ operations. A cautionary remark is however in order:
the MeshBlock/CPU saturation has not been taken into account, which requires a
certain threshold for efficiency Daszuta  _et al._ (2021); Cook  _et al._
(2023). With regard to performance, these timing results nonetheless strongly
suggest that $\mathcal{Z}_{\mathrm{VC}}$ is the better option for BBH.

###  III.2 Binary Neutron Stars

In the case of binary neutron star (BNS) evolution, in addition to the space-
time evolution, we need to address the evolution of the hydrodynamical
variables. These variables satisfy a balance law of the form
$\partial_{t}[\mathcal{Q}_{A}]+\partial_{i}\mathcal{F}{}_{A}^{i}[\mathcal{Q}]=\mathcal{S}_{A}[\mathcal{Q}]$
Banyuls  _et al._ (1997). The source term $\mathcal{S}_{A}[\mathcal{Q}]$ does
not involve derivatives of $\mathcal{Q}_{A}$, but does contain geometric
fields and their derivatives. Here the finite-volume, HRSC methods we employ
lead us to take the evolved $\mathcal{Q}_{A}$ on cell-centers. In the
$\partial_{i}\mathcal{F}{}_{A}^{i}[\mathcal{Q}]$ approximation step, primitive
variables $\mathcal{R}_{A}[Q]$ are assembled and reconstructed to FC Zanna and
Bucciantini (2002). This process requires geometric data on FC. Information
about the evolved fluid variable state, must be available on the grid used for
the space-time solver. If we select $\mathcal{Z}_{\mathrm{VC}}$ the required
intergrid transfer operations can be summarized as:
$\mathcal{Z}_{\mathrm{VC}\rightarrow\mathrm{CC}}$,
$\mathcal{Z}_{\mathrm{VC}\rightarrow\mathrm{FC}}$, and finally
$\mathcal{R}_{A,\mathrm{CC}\rightarrow\mathrm{VC}}$. On the other hand,
selecting $\mathcal{Z}_{\mathrm{CX}}$ means we only need
$\mathcal{Z}_{\mathrm{CX}\rightarrow\mathrm{FC}}$.

The intergrid transfer operations summarized above are handled through the use
of symmetric Lagrange interpolation. Specifically, during evaluation, $N_{I}$
nearest neighbor nodes either side of a target interpolation point are used
along a salient axis. We have found it beneficial, for stability reasons, to
fix $N_{I}=1$ in a subset of operations. This we refer to as a hybrid mode. In
this mode $N_{I}$ is only allowed to vary freely for $\mathcal{Z}$ derivative
assembly entering $\mathcal{S}_{A}[\mathcal{Q}]$, and primitives
$\mathcal{R}_{A,\mathrm{CC}\rightarrow\mathrm{VC}}$ that recouple to the
space-time solver.

As initial consistency checks for geometry-matter field coupling and the above
operations we consider static star tests in Appendix A. A more representative
setting of binary neutron star (BNS) evolution problems is required. For this,
we consider a quasi-circular, equal-mass BNS inspiral-merger problem. Here
irrotational, constraint-satisfying initial data is provided by Lorene
Gourgoulhon  _et al._ (2001). Specifically, as a validation problem, we
consider the G2_I14vs14_D4R33_45km binary777Dataset available at
https://lorene.obspm.fr/.  which features baryon mass $M_{b}=3.2500{\rm
M_{\odot}}$ and gravitational mass $M=3.0297{\rm M_{\odot}}$ at an initial
separation of $45$ km. The Arnowitt-Deser-Misner (ADM) mass of the binary is
$M_{\rm ADM}=2.9984{\rm M_{\odot}}$, the angular momentum $J_{\rm
ADM}=8.8354{\rm M_{\odot}}^{2}$G/c, and the initial orbital frequency is
$f_{0}\simeq 294$ Hz. For these tests we utilize an ideal gas equation-of-
state (see Appendix A) with $K=123.64$. The evolution of this BNS system
proceeds for three orbits prior to the formation of a massive remnant,
thereafter leading to gravitational collapse.

The computational domain is again selected as $\Omega:=[-D,\,D]^{3}$, where
$D=1536\,{\rm M_{\odot}}{}\simeq 2268\,\mathrm{km}$. The AMR is induced
through tracking the local minimum of the lapse $\alpha$, as positioned on
each constituent of the BNS system. Resolution levels are fixed according to
$\\{\mathbb{S}_{8,15}(\alpha^{*}_{\pm}),\,\mathbb{S}_{7,30}(\alpha^{*}_{\pm}),\,\mathbb{S}_{6,110}(\mathbf{0})\\}$.
To perform a resolution study, we take a sequence of
$N_{M}\in\\{64,\,96,\,128\\}$. This translates into a sequence of finest
resolution levels of $(554,\,369,\,277)\,\mathrm{m}$. The MeshBlock sampling
is selected as $N_{B}=16$. For the time-evolution we make use of the
SSPRK$(3,3)$ method of Gottlieb  _et al._ (2009), with a CFL of
$\mathcal{C}_{1}=0.2$, and $\sigma_{\mathrm{D}}=0.5$. A tenuous, artificial
atmosphere is added, where $\rho_{\mathrm{atm}}=10^{-16}$. We do not make use
of thresholding Cook  _et al._ (2023). For the RePrimAnd-based primitive
recovery we set the tolerance to $10^{-10}$. In this section no symmetry
reductions are performed.

We select $\mathcal{Z}_{\mathrm{CX}}$ and evolve the initial data with the
aforementioned parameters. An additional run featuring
$\mathcal{Z}_{\mathrm{VC}}$ with $N_{M}=64$ and all else fixed is also
performed for internal comparison. In each case we find approximately three
orbits prior to the formation of a massive remnant. A snapshot of the fluid
rest mass-density at this point ($t\simeq 8.8\mathrm{ms}$) is shown in Fig.9
comparing the four runs.

Figure 9: Mesh structure and values of fluid rest mass-density
$\rho\,[\mathrm{g}\,\mathrm{cm}^{-3}]$ shortly post-merger in the orbital
plane for the BNS inspiral problem described in the text. Counter-clockwise,
starting from the upper right quadrant we have evolution with CX sampling with
$N_{M}\in\\{64,\,96,\,128\\}$. The lower right quadrant depicts evolution
under VC sampling for geometric fields with $N_{M}=64$. Each grey box
(MeshBlock) is sampled with $N_{B}=16$. The AMR criterion for each run has
driven the distribution of refined region resolution towards an origin-
centered approximate sphere. Trackers that follow local minima of $\alpha(t)$
are overlaid. Darker tracker color indicates increasing $t$.

We observe that the spiral arm structure is qualitatively similar for the two
runs at $N_{M}=64$, which provides an initial consistency check. The structure
becomes sharper and more compact at higher $N_{M}$. The extrema tracking based
AMR can also be seen to guide increased resolution towards regions that
suitably follow the binary consituents.

Figure 10: Relative change of (global) maximum density $\rho_{c}$ (upper
panel), baryon mass conservation $M_{b}$ (middle), and minimum of the lapse
$\alpha$ (lower panel), for the BNS model problem. A resolution triplet
induced by the choices $N_{M}\in\\{64,\,96,\,128\\}$ and selection of
$\mathcal{Z}_{\mathrm{CX}}$ is depicted, with a single run at $N_{M}=64$ and
$\mathcal{Z}_{\mathrm{VC}}$ shown for reference. Collapse times, for each run,
are signalled by: post-merger, global peak in the $\rho_{c}$ monitor, blow up
of the mass conservation monitor, and confirmed by reduction of lapse below
$\sim 0.2$. Vertical black line depicts merger time ($t\simeq
8.4\,\mathrm{ms}$), as based on the peak amplitude $h{}_{22}$ associated with
the $N_{M}=128$ run (C.f. Fig.11). See text for further discussion. Note:
legend common to all panels.

Quantitative scalar monitors are shown in Fig.10. We find that the time of
collapse falls in the range $t\in[16.5,\,20.2]\,\mathrm{ms}$ for
$\mathcal{Z}_{\mathrm{CX}}$ selection. This depends on the selection of
$N_{M}$. Importantly, we observe consistency to within $1\%$ when comparing a
verification run with $N_{M}=64$ based on $\mathcal{Z}_{\mathrm{VC}}$.
Violation in baryon mass conservation grows to a plateau of
$\sim\mathcal{O}(10^{-6})$ at merger for each run. This remains approximately
constant until the point of gravitational collapse. Again we see that both
choices of grid sampling for the geometric fields lead to consistent results.

Gravitational waveforms are the next item to verify. To construct the
$(2,\,2)$ mode of the strain, for the BNS system, we use fixed frequency
integration (FFI) Reisswig and Pollney (2011). The FFI proceeds by
transforming the signal to the frequency representation, applying a low-
frequency cutoff $f_{0}$, and integrating. We select $f_{0}=0.002M^{-1}$. The
result is shown in Fig.11.

Figure 11: Gravitational strain $(2,\,2)$ mode extracted at $R=100M$ for
binary neutron star run of §III.2. (Top): real part of the strain with
$\mathcal{Z}_{\mathrm{CX}}$ for $N_{M}\in\\{64,\,96,\,128\\}$ is depicted. For
$N_{M}=64$ the amplitude is also shown. In the case of $N_{M}=128$ we find
$u_{\mathrm{merger}}\simeq 1583.2\,M$ ($t_{\mathrm{merger}}\simeq 1699.8\,M$).
We observe consistent behaviour at $N_{M}=64$ for a comparable run based on VC
sampling. (Middle): Instantaneous frequency is shown. (Top, middle): vertical
dashed lines indicate time of merger for the different data sets. (Bottom):
Phase differences for the CX resolution triplet suggest convergence (to
merger) as compatible with a second order trend.

Based on the $N_{M}=128$ run, we find that the time of merger (retarded
coordinate) is $u_{\mathrm{merger}}\simeq 1583.2\,M$. This translates into
$t_{\mathrm{merger}}\simeq 1699.8\,M\simeq 7.9\,\mathrm{ms}$. Indeed, all CX
runs fall within a band of $\sim 0.2\mathrm{ms}$. This indicates strong
robustness of the simulations performed. The merger time for the
$\mathcal{Z}_{\mathrm{VC}}$ waveform is within $3\%$ that of the CX at
$N_{M}=96$. Similarly, the morphology of the instantaneous frequency for CX
when cross-compared against VC is seen to be alike. With the method described
in Appendix A we also assess convergence. As can be seen in Fig.11 we have a
strong indication for a second order trend in the convergence of the
$h{}_{22}$ associated phase, as $N_{M}$ is increased.

Overall, based on this discussion, when running with
$\mathcal{Z}_{\mathrm{CX}}$, good agreement is found with
$\mathcal{Z}_{\mathrm{VC}}$ here, and our prior results in Cook  _et al._
(2023). This provides confidence in the correctness of our implementation, and
the validity of the strategy outlined in §II.2.

To close this section, as in the BBH case, we provide a rough characterization
of performance differences between $\mathcal{Z}_{\mathrm{CX}}$ and
$\mathcal{Z}_{\mathrm{VC}}$ through a short timing test. The initial
checkpoint is taken at a time of $t_{0}=5\,M$ with timings summarized in
Tab.3.

$N_{M}$ | $t_{F}(\mbox{{\tt cx}})$ | $t_{F}(\mbox{{\tt vc}})$ |  $r_{t}$ ($2\,\mathrm{dp}$) | $\\#\mathrm{MB}$ | $\\#\mathrm{MB}/\mathrm{Core}$
---|---|---|---|---|---
$64$ | $170.45$ | $151.02$ | $1.13$ | $3676$ | $\sim 5$
$96$ | $70.84$ | $61.11$ | $1.17$ | $6768$ | $\sim 9$
$128$ | $34.78$ | $28.36$ | $1.27$ | $12384$ | $\sim 16$
Table 3: One hour timing test on SuperMUC-NG of LRZ utilizing $16$ nodes. The
timing ratio $r$ is defined as $r_{t}:=(t_{F}(\mbox{{\tt
cx}})-t_{0})/(t_{F}(\mbox{{\tt vc}})-t_{0})$. The number of MeshBlock objects
is provided at $t_{0}$. No symmetry reduction has been performed. The AMR is
based on spherical regions criteria.

We see that for BNS, the space-time solver has a distinct advantage when run
with the choice $\mathcal{Z}_{\mathrm{CX}}$ in terms of speed. Indeed, we
observe an approximate average $\sim 17\%$ later evolution time, across
$N_{M}$, when compared to the $\mathcal{Z}_{\mathrm{VC}}$ choice. This is in
direct contrast to the BBH case (Tab.2). While it is the case that the level-
to-level operators become more involved for CX, when contrasted against VC,
the reduction in the number of intergrid transfer operations appears to
compensate. Indeed, it appears that the gap between CX and VC ($r_{t}$ of
Tab.3) increases with increasing $N_{M}$. This is also in contrast to the BBH
case, where the gap appeared to decrease. Overall, we see that CX achieves
later final evolved times in the BNS problem. Hence we conclude that
$\mathcal{Z}_{\mathrm{CX}}$ is the favoured choice, for this class of problem.

##  IV Conclusion

In this work we have proposed a new, cell-centered (CX) variable treatment
within the context of our space-time solver embedded within GR-Athena++. This
has allowed us to solve the $\mathrm{Z}4$c system in
$\mathcal{Z}_{\mathrm{CX}}$ mode, in a fashion consistent with our earlier
work on vertex-centers, i.e. $\mathcal{Z}_{\mathrm{VC}}$ mode. Tests for both
vacuum and GRHD problems demonstrated robust functionality for either grid
sampling. The choice of $\mathcal{Z}_{\mathrm{CX}}$ for GR(M)HD yields a
substantial simplification in treatment of discretized variables.

To validate the new treatment, together with implementation details within the
code, first consistency tests based on BBH evolution were performed. A step-
by-step strategy was taken which featured, inter alia, a resolution study.
This allowed us to assess a trend compatible with $6^{\mathrm{th}}$ order
convergence in the error of the phase of the $(2,\,2)$ mode of the
gravitational wave strain, when working with evolution of
$\mathcal{Z}_{\mathrm{CX}}$. This is consistent with our prior work Daszuta
_et al._ (2021). In order to provide insight as to which grid sampling
strategy performs better, for BBH evolution, we considered a one hour timing
test on SuperMUC-NG of the Leibniz-Rechenzentrum Munich. It was found that
$\mathcal{Z}_{\mathrm{VC}}$, on average, is $\sim 20\%$ faster. This we
attribute to the differences in how the level-to-level refinement operators
are evaluated for the two grid samplings. Motivated by the “symmetry-seeking”
property of the $L_{2}$-AMR previously investigated in Rashti  _et al._
(2024), we proposed, a simplified, spherical-region based AMR strategy. This
was demonstrated in the BBH validation problem to potentially partially
mitigate some rectilinear grid imprinting.

Following this, a quasi-circular, equal-mass, BNS inspiral, merger, and
gravitational collapse problem was illustrated in the context of GRHD. We
again utilized the $(2,\,2)$ mode of the gravitational wave strain as a
consistency diagnostic. A trend compatible with $2^{\mathrm{nd}}$ order
convergence was assessed in a resolution study when working with
$\mathcal{Z}_{\mathrm{CX}}$. This trend is consistent with that observed in
Cook  _et al._ (2023) when working with $\mathcal{Z}_{\mathrm{VC}}$. With
regard to performance, in contrast to the BBH case, it was found that
$\mathcal{Z}_{\mathrm{CX}}$ is, on average, $\sim 17\%$ faster for the runs
considered. This we attribute to a reduction in the total number of intergrid
interpolation operations required. For these tests, we also leveraged the
simplified, spherical-region based AMR strategy. While we do not explicitly
demonstrate GRMHD runs in this work, we have found another potential advantage
of working with $\mathcal{Z}_{\mathrm{CX}}$. First tests, involving the
initial data of the BNS simulation presented here, allowed for evolution with
magnetic fields through merger and collapse, without apparent horizon-based
excision. This contrasts against the situation with
$\mathcal{Z}_{\mathrm{VC}}$ where we have found such excision techniques to be
crucial.

###### Acknowledgements.

BD and SB acknowledges support by the EU Horizon under ERC Consolidator Grant,
no. InspiReM-101043372. PH acknowledges funding from the National Science
Foundation under Grant No. PHY-2116686. EG, DR acknowledge funding from the
National Science Foundation under grant No. AST-2108467. EG acknowledges
funding from an Institute for Gravitation and Cosmology fellowship. JF, DR
acknowledge U.S. Department of Energy, Office of Science, Division of Nuclear
Physics under Award Number(s) DE-SC0021177. DR acknowledge support from NASA
under award No. 80NSSC21K1720 and from the Sloan Foundation. The authors are
indebted to A. Celentano’s PRİSENCÓLİNENSİNÁİNCIÚSOL. The numerical
simulations were performed on the national HPE Apollo Hawk at the High
Performance Computing Center Stuttgart (HLRS). The authors acknowledge HLRS
for funding this project by providing access to the supercomputer HPE Apollo
Hawk under the grant numbers INTRHYGUE/44215 and MAGNETIST/44288. Simulations
were also performed on SuperMUC_NG at the Leibniz-Rechenzentrum (LRZ) Munich.
The authors acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-
centre.eu) for funding this project by providing computing time on the GCS
Supercomputer SuperMUC-NG at LRZ (allocations pn68wi and pn36jo). Simulations
were also performed on TACC’s Frontera (NSF LRAC allocation PHY23001) and on
Perlmutter. This research used resources of the National Energy Research
Scientific Computing Center, a DOE Office of Science User Facility supported
by the Office of Science of the U.S. Department of Energy under Contract No.
DE-AC02-05CH11231. Postprocessing and development run were performed on the
ARA cluster at Friedrich Schiller University Jena. The ARA cluster is funded
in part by DFG grants INST 275/334-1 FUGG and INST 275/363-1 FUGG, and ERC
Starting Grant, grant agreement no. BinGraSp-714626.

## References

  * Abbott  _et al._ (2017a) B. P. Abbott _et al._ (Virgo, LIGO Scientific), Phys. Rev. Lett. 119, 161101 (2017a), arXiv:1710.05832 [gr-qc] .
  * Goldstein  _et al._ (2017) A. Goldstein _et al._ , Astrophys. J. 848, L14 (2017), arXiv:1710.05446 [astro-ph.HE] .
  * Savchenko  _et al._ (2017) V. Savchenko _et al._ , Astrophys. J. 848, L15 (2017), arXiv:1710.05449 [astro-ph.HE] .
  * Akutsu  _et al._ (2020) T. Akutsu _et al._ , arXiv:2009.09305 [astro-ph, physics:gr-qc] (2020), arXiv:2009.09305 [astro-ph, physics:gr-qc] .
  * Amaro-Seoane (2017) e. a. Amaro-Seoane, arXiv:1702.00786 [astro-ph] (2017), arXiv:1702.00786 [astro-ph] .
  * Punturo  _et al._ (2010) M. Punturo, M. Abernathy, F. Acernese, B. Allen, N. Andersson, _et al._ , Class.Quant.Grav. 27, 194002 (2010).
  * Abbott  _et al._ (2017b) B. P. Abbott _et al._ (LIGO Scientific), Class. Quant. Grav. 34, 044001 (2017b), arXiv:1607.08697 [astro-ph.IM] .
  * Nagar  _et al._ (2018) A. Nagar _et al._ , Phys. Rev. D98, 104052 (2018), arXiv:1806.01772 [gr-qc] .
  * Font (2007) J. A. Font, Living Rev. Rel. 11, 7 (2007).
  * Shibata  _et al._ (2011) M. Shibata, K. Kiuchi, Y.-i. Sekiguchi, and Y. Suwa, Prog.Theor.Phys. 125, 1255 (2011), arXiv:1104.3937 [astro-ph.HE] .
  * Radice  _et al._ (2022) D. Radice, S. Bernuzzi, A. Perego, and R. Haas, Mon. Not. Roy. Astron. Soc. 512, 1499 (2022), arXiv:2111.14858 [astro-ph.HE] .
  * Hilditch  _et al._ (2016) D. Hilditch, A. Weyhausen, and B. Brügmann, Phys. Rev. D93, 063006 (2016), arXiv:1504.04732 [gr-qc] .
  * Bugner  _et al._ (2016) M. Bugner, T. Dietrich, S. Bernuzzi, A. Weyhausen, and B. Brügmann, Phys. Rev. D94, 084004 (2016), arXiv:1508.07147 [gr-qc] .
  * Fernando  _et al._ (2018) M. Fernando, D. Neilsen, H. Lim, E. Hirschmann, and H. Sundar,  (2018), 10.1137/18M1196972, arXiv:1807.06128 [gr-qc] .
  * Löffler  _et al._ (2012) F. Löffler _et al._ , Class. Quant. Grav. 29, 115001 (2012), arXiv:1111.3344 [gr-qc] .
  * Shankar  _et al._ (2023) S. Shankar, P. Mösta, S. R. Brandt, R. Haas, E. Schnetter, and Y. de Graaf, Class. Quant. Grav. 40, 205009 (2023), arXiv:2210.17509 [astro-ph.IM] .
  * Clough  _et al._ (2015) K. Clough, P. Figueras, H. Finkel, M. Kunesch, E. A. Lim, and S. Tunyasuvunakool, (2015), arXiv:1503.03436 [gr-qc] .
  * Tichy  _et al._ (2023) W. Tichy, L. Ji, A. Adhikari, A. Rashti, and M. Pirog, Classical and Quantum Gravity 40, 025004 (2023), arxiv:2212.06340 [gr-qc] .
  * Kidder  _et al._ (2017) L. E. Kidder _et al._ , J. Comput. Phys. 335, 84 (2017), arXiv:1609.00098 [astro-ph.HE] .
  * Rosswog and Diener (2021) S. Rosswog and P. Diener, Class. Quant. Grav. 38, 115002 (2021), arXiv:2012.13954 [gr-qc] .
  * Grandclement and Novak (2009) P. Grandclement and J. Novak, Living Rev. Rel. 12, 1 (2009), arXiv:0706.2286 [gr-qc] .
  * Doulis  _et al._ (2022) G. Doulis, F. Atteneder, S. Bernuzzi, and B. Brügmann, Phys. Rev. D 106, 024001 (2022), arXiv:2202.08839 [gr-qc] .
  * Alfieri  _et al._ (2018) R. Alfieri, S. Bernuzzi, A. Perego, and D. Radice, Journal of Low Power Electronics and Applications Special Issue ”Energy Aware Scientific Computing on Low Power and Heterogeneous Architectures”, 8 (2018), 10.3390/jlpea8020015.
  * Stout  _et al._ (1997) Q. F. Stout, D. L. De Zeeuw, T. I. Gombosi, C. P. T. Groth, H. G. Marshall, and K. G. Powell, in  _Proceedings of the 1997 ACM/IEEE Conference on Supercomputing_, SC ’97 (Association for Computing Machinery, New York, NY, USA, 1997) pp. 1–10.
  * Berger and Oliger (1984) M. J. Berger and J. Oliger, J.Comput.Phys. 53, 484 (1984).
  * Burstedde  _et al._ (2019) C. Burstedde, J. Holke, and T. Isaac, Foundations of Computational Mathematics 19, 843 (2019).
  * Morton (1966) G. M. Morton,  _A computer oriented geodetic data base and a new technique in file sequencing_ , Tech. Rep. (1966).
  * Nakamura  _et al._ (1987) T. Nakamura, K. Oohara, and Y. Kojima, Prog. Theor. Phys. Suppl. 90, 1 (1987).
  * Shibata and Nakamura (1995) M. Shibata and T. Nakamura, Phys. Rev. D52, 5428 (1995).
  * Baumgarte and Shapiro (1999) T. W. Baumgarte and S. L. Shapiro, Phys. Rev. D59, 024007 (1999), arXiv:gr-qc/9810065 .
  * Bernuzzi and Hilditch (2010) S. Bernuzzi and D. Hilditch, Phys. Rev. D81, 084003 (2010), arXiv:0912.2920 [gr-qc] .
  * Hilditch  _et al._ (2013) D. Hilditch, S. Bernuzzi, M. Thierfelder, Z. Cao, W. Tichy, and B. Bruegmann, Phys. Rev. D88, 084057 (2013), arXiv:1212.2901 [gr-qc] .
  * Friedrich (1985) H. Friedrich, Communications in Mathematical Physics 100, 525 (1985).
  * Pretorius (2005) F. Pretorius, Phys. Rev. Lett. 95, 121101 (2005), arXiv:gr-qc/0507014 .
  * Lindblom  _et al._ (2006) L. Lindblom, M. A. Scheel, L. E. Kidder, R. Owen, and O. Rinne, Class.Quant.Grav. 23, S447 (2006), arXiv:gr-qc/0512093 [gr-qc] .
  * Lousto and Healy (2020) C. O. Lousto and J. Healy, Phys. Rev. Lett. 125, 191102 (2020), arXiv:2006.04818 [gr-qc] .
  * Lannelongue  _et al._ (2023) L. Lannelongue, H.-E. G. Aronson, A. Bateman, E. Birney, T. Caplan, M. Juckes, J. McEntyre, A. D. Morris, G. Reilly, and M. Inouye, Nature Computational Science 3, 514 (2023).
  * Banyuls  _et al._ (1997) F. Banyuls, J. A. Font, J. M. A. Ibanez, J. M. A. Marti, and J. A. Miralles, Astrophys. J. 476, 221 (1997).
  * Daszuta  _et al._ (2021) B. Daszuta, F. Zappa, W. Cook, D. Radice, S. Bernuzzi, and V. Morozova, Astrophys. J. Supp. 257, 25 (2021), arXiv:2101.08289 [gr-qc] .
  * Cook  _et al._ (2023) W. Cook, B. Daszuta, J. Fields, P. Hammond, S. Albanesi, F. Zappa, S. Bernuzzi, and D. Radice, (2023), arXiv:2311.04989 [gr-qc] .
  * Daszuta and Cook (2024) B. Daszuta and W. Cook, “GR-Athena++: Magnetohydrodynamical evolution with dynamical space-time,”  (2024), arxiv:2406.05126 [astro-ph, physics:gr-qc] .
  * White  _et al._ (2016) C. J. White, J. M. Stone, and C. F. Gammie, The Astrophysical Journal Supplement Series 225, 22 (2016), arXiv:1511.00943 .
  * Felker and Stone (2018) K. G. Felker and J. M. Stone, Journal of Computational Physics 375, 1365 (2018).
  * Stone  _et al._ (2020) J. M. Stone, K. Tomida, C. J. White, and K. G. Felker, The Astrophysical Journal Supplement Series 249, 4 (2020).
  * Brügmann  _et al._ (2008) B. Brügmann, J. A. Gonzalez, M. Hannam, S. Husa, U. Sperhake, _et al._ , Phys.Rev. D77, 024027 (2008), arXiv:gr-qc/0610128 [gr-qc] .
  * Mösta  _et al._ (2014) P. Mösta, B. C. Mundim, J. A. Faber, R. Haas, S. C. Noble, _et al._ , Class.Quant.Grav. 31, 015005 (2014), arXiv:1304.5544 [gr-qc] .
  * Radice  _et al._ (2014) D. Radice, L. Rezzolla, and F. Galeazzi, Class.Quant.Grav. 31, 075012 (2014), arXiv:1312.5004 [gr-qc] .
  * Thierfelder  _et al._ (2011a) M. Thierfelder, S. Bernuzzi, and B. Brügmann, Phys.Rev. D84, 044012 (2011a), arXiv:1104.4751 [gr-qc] .
  * Fleischmann  _et al._ (2019) N. Fleischmann, S. Adami, and N. A. Adams, Computers & Fluids 189, 94 (2019).
  * Neumaier (1974) A. Neumaier, ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik 54, 39 (1974).
  * Stone  _et al._ (2020) J. M. Stone, K. Tomida, C. J. White, and K. G. Felker, Astrophys. J. Suppl. 249, 4 (2020), arXiv:2005.06651 [astro-ph.IM] .
  * Berrut and Trefethen (2004) J.-P. Berrut and L. N. Trefethen, SIAM Review 46, 501 (2004).
  * Rashti  _et al._ (2024) A. Rashti, M. Bhattacharyya, D. Radice, B. Daszuta, W. Cook, and S. Bernuzzi, Class. Quant. Grav. 41, 095001 (2024), arXiv:2312.05438 [gr-qc] .
  * Campanelli  _et al._ (2006) M. Campanelli, C. O. Lousto, P. Marronetti, and Y. Zlochower, Phys. Rev. Lett. 96, 111101 (2006), arXiv:gr-qc/0511048 .
  * Noble  _et al._ (2006) S. C. Noble, C. F. Gammie, J. C. McKinney, and L. Del Zanna, Astrophys. J. 641, 626 (2006), arXiv:astro-ph/0512420 .
  * Kastaun  _et al._ (2021) W. Kastaun, J. V. Kalinani, and R. Ciolfi, Phys. Rev. D 103, 023018 (2021), arXiv:2005.01821 [gr-qc] .
  * Borges  _et al._ (2008) R. Borges, M. Carmona, B. Costa, and W. S. Don, Journal of Computational Physics 227, 3191 (2008).
  * Zanna and Bucciantini (2002) L. D. Zanna and N. Bucciantini, Astronomy & Astrophysics 390, 1177 (2002).
  * Bona  _et al._ (1996) C. Bona, J. Massó, J. Stela, and E. Seidel, in  _The Seventh Marcel Grossmann Meeting: On Recent Developments in Theoretical and Experimental General Relativity, Gravitation, and Relativistic Field Theories_ , edited by R. T. Jantzen, G. M. Keiser, and R. Ruffini (World Scientific, Singapore, 1996).
  * Alcubierre  _et al._ (2003) M. Alcubierre, B. Brügmann, P. Diener, M. Koppitz, D. Pollney, _et al._ , Phys.Rev. D67, 084023 (2003), arXiv:gr-qc/0206072 [gr-qc] .
  * van Meter  _et al._ (2006) J. R. van Meter, J. G. Baker, M. Koppitz, and D.-I. Choi, Phys. Rev. D73, 124011 (2006), arXiv:gr-qc/0605030 .
  * Ansorg  _et al._ (2004) M. Ansorg, B. Brügmann, and W. Tichy, Phys. Rev. D70, 064011 (2004), arXiv:gr-qc/0404056 .
  * Ketcheson (2010) D. I. Ketcheson, Journal of Computational Physics 229, 1763 (2010).
  * Wang and Lee (2011) N. Wang and J.-L. Lee, SIAM Journal of Scientific Computing 33, 2536 (2011).
  * Goldberg  _et al._ (1967) J. N. Goldberg, A. J. MacFarlane, E. T. Newman, F. Rohrlich, and E. C. G. Sudarshan, J. Math. Phys. 8, 2155 (1967).
  * LIGO Scientific Collaboration (2018) LIGO Scientific Collaboration, “LIGO Algorithm Library - LALSuite,” free software (GPL) (2018).
  * Albanesi  _et al._ (2024) S. Albanesi, A. Rashti, F. Zappa, R. Gamba, W. Cook, B. Daszuta, S. Bernuzzi, A. Nagar, and D. Radice, “Scattering and dynamical capture of two black holes: Synergies between numerical and analytical methods,”  (2024), arxiv:2405.20398 [gr-qc] .
  * Thierfelder  _et al._ (2011b) M. Thierfelder, S. Bernuzzi, D. Hilditch, B. Brügmann, and L. Rezzolla, Phys.Rev. D83, 064022 (2011b), arXiv:1012.3703 [gr-qc] .
  * Mewes  _et al._ (2020) V. Mewes, Y. Zlochower, M. Campanelli, T. W. Baumgarte, Z. B. Etienne, F. G. Lopez Armengol, and F. Cipolletta, Phys. Rev. D 101, 104007 (2020), arXiv:2002.06225 [gr-qc] .
  * Gourgoulhon  _et al._ (2001) E. Gourgoulhon, P. Grandclement, K. Taniguchi, J.-A. Marck, and S. Bonazzola, Phys.Rev. D63, 064029 (2001), arXiv:gr-qc/0007028 [gr-qc] .
  * Gottlieb  _et al._ (2009) S. Gottlieb, D. I. Ketcheson, and C.-W. Shu, Journal of Scientific Computing 38, 251 (2009).
  * Reisswig and Pollney (2011) C. Reisswig and D. Pollney, Class.Quant.Grav. 28, 195015 (2011), arXiv:1006.1632 [gr-qc] .
  * Font  _et al._ (2000) J. A. Font, M. A. Miller, W.-M. Suen, and M. Tobias, Phys. Rev. D61, 044011 (2000), gr-qc/9811015 .

##  Appendix A Assessing convergence with a static star

Evolving a stable, static star in the general relativistic setting provides
for a simple, sanity check on code correctness. To this end initial data (ID)
is prepared based on the assumption of a spherically symmetric, self-
gravitating matter distribution at equilibrium satisfying the Tolman-
Oppenheimer-Volkoff (TOV) equations (see e.g. Font  _et al._ (2000)). For ID
preparation a polytropic EOS $p=K\rho^{\Gamma}$ is assumed with adiabatic
index $\Gamma=2$, and polytropic constant $K=100$ as in Bernuzzi and Hilditch
(2010). The stable star model is taken to have a central density
$\rho_{c}=1.280\times 10^{-3}$, gravitational mass $M=1.400$, and
circumferential radius $R=9.586$. During evolution we make use of the more
general EOS $p=(\Gamma-1)\rho\epsilon$ with
$\epsilon=K\rho^{\Gamma-1}/(\Gamma-1)$. We perform an evolution of this data
to a final time of $T=1500$. This is done for both $\mathcal{Z}_{\mathrm{VC}}$
and $\mathcal{Z}_{\mathrm{CX}}$ choices, for a variety of intergrid
interpolation orders (see §III.2). For these runs we have taken a CFL of
$\mathcal{C}_{1}=0.25$, Kreiss-Oliger dissipation $\sigma_{\mathrm{D}}=0.5$,
and set artificial atmosphere to $\rho_{\mathrm{atm}}=10^{-18}$. The
$\mathrm{Z}4$c damping parameters are taken as
$(\kappa_{1},\,\kappa_{2})=(0.02,\,0)$. Scalar diagnostic quantities are shown
in Fig.12.

Figure 12: Comparison of scalar quantities in stable star evolution (TOV test)
on a fixed grid configuration involving CX and VC sampling where various
varieties of intergrid interpolation of §III.2 are utilized. The Mesh and
MeshBlock samplings are $N_{M}=64$ and $N_{B}=16$ respectively. The
computational domain is taken as $\Omega:=[-1024,\,1024]^{3}$ with a central,
cubic, region $\Omega_{r}:=[-10,\,10]^{3}$ refined with $N_{L}=7$. The region
$\Omega_{r}$ completely contains the star. As intergrid interpolation order is
increased it can be seen (upper panel) that central density is better
preserved, with oscillation magnitude of $\rho_{c}(t)$ reduced. Hamiltonian
constraint violation (lower panel) is also reduced as interpolation order is
increased. We find common mass conservation (middle panel) properties after
$\sim t=200$ regardless of interpolation order, or geometric sampling.
Vertical grey line indicates first time of causal contact of $\partial\Omega$
with the origin.

To assess solution convergence we consider numerical evolution repeated at a
triplet of coarse, medium, and fine resolutions $(\delta x_{c},\,\delta
x_{m},\,\delta x_{f})$ that satisfy $\delta x_{c}>\delta x_{m}>\delta x_{f}$.
The convergence rate of an approximation to the corresponding field data
$\mathcal{F}$ may be investigated by comparing differences of solutions
$\delta\mathcal{F}_{ab}:=\mathcal{F}_{a}-\mathcal{F}_{b}$ at distinct
resolutions. For $\delta x^{n}$ one finds based on Taylor expansion that
$\delta\mathcal{F}_{cm}\simeq\delta\mathcal{F}_{mf}\times Q_{n}(c,m,m,f)$
where we have introduced the so-called convergence factor:

$\displaystyle Q_{n}(a,\,b,\,c,\,d):=\frac{\delta x_{a}^{n}-\delta
x_{b}^{n}}{\delta x_{c}^{n}-\delta x_{d}^{n}}.$ (4)

The $N_{I}=2$ simulations featuring CX and VC of Fig.12 show the smallest
deviation from the initial value $\rho_{c}$ during evolution, together with
the smallest constraint violation. We construct a triplet of resolutions for
convergence assessment by repeating the runs with $N_{L}\in\\{5,\,6\\}$. Based
on the refinement structure this entails $(\delta x_{c},\,\delta
x_{m},\,\delta x_{f})=(\delta x_{c},\,\delta x_{c}/2,\,\delta x_{c}/4)$ on the
finest level, covering the star, for each run. Behaviour for the central
density, together with the Hamiltonian constraint for our convergence
assessment is shown in Fig.13. As can be seen we observe $\sim
2^{\mathrm{nd}}$ order, and $\sim 5^{\mathrm{th}}$ order compatible trends for
the aforementioned quantities respectively.

Figure 13: Convergence of scalar quantities in the TOV test for a fixed grid
configuration involving CX and VC sampling. Overall grid structure is as in
Fig.12 where resolution triplets are constructed based selecting the level of
$\Omega_{r}$ refinement through $N_{L}\in\\{5,\,6,\,7\\}$. In dashed lines we
depict the rescaled coarse-medium difference according to $Q_{n}$ of Eq.(4).
As additionally verified through least-squares fitting on $t\in[128,\,950]$ we
find: $n\sim 2$ i.e. $\sim 2^{\mathrm{nd}}$ order convergence for $\rho_{c}$
(upper panel); $n\sim 5$ for the Hamiltonian constraint (lower panel).

Based on the consistency between the qualitative features of the CX and VC
scalar diagnostics depicted in Fig.12, together with consistent order of
convergence demonstrated in Fig.13 we consider the TOV test to pass.