# Quantifying uncertainties in general relativistic magnetohydrodynamic codes

Pedro L. Espino  Department of Physics, The Pennsylvania State University,
University Park, PA 16802, USA  Department of Physics, University of
California, Berkeley, CA 94720, USA  Gabriele Bozzola  Department of
Astronomy, University of Arizona, Tucson, AZ 85721,USA  Vasileios Paschalidis
Department of Astronomy, University of Arizona, Tucson, AZ 85721,USA
Department of Physics, University of Arizona, Tucson, AZ 85721, USA

###### Abstract

In this paper, we show that similar open-source codes for general relativistic
magnetohydrodynamic (GRMHD) produce different results for key features of
binary neutron star mergers. First, we present a new open-source version of
the publicly available IllinoisGRMHD code that provides support for realistic,
finite temperature equations of state. After stringent tests of our upgraded
code, we perform a code comparison between GRHydro, IllinoisGRMHD, Spritz, and
WhiskyTHC, which implement the same physics, but slightly different
computational methods. The benefit of the comparison is that all codes are
embedded in the EinsteinToolkit suite, hence their only difference is
algorithmic. We find similar convergence properties, fluid dynamics, and
gravitational waves, but different merger times, remnant lifetimes, and
gravitational wave phases. Such differences must be resolved before the post-
merger dynamics modeled with such simulations can be reliably used to infer
the properties of nuclear matter especially in the era of precision
gravitational wave astronomy.

##  I Introduction

The detection of gravitational waves (GWs) from the merger of a binary neutron
star (BNS) Abbott  _et al._ (2017) has established the need for a clearer
understanding of BNS systems. To gain insight into the physics of BNS mergers,
observations associated with such events – be they in the GW or
electromagnetic (EM) spectrum – must be analyzed with either simplified
analytical models which parameterize the dominant physical phenomena in the
system, e.g., Raaijmakers  _et al._ (2021); Easter  _et al._ (2020); Barbieri
_et al._ (2020); Breschi  _et al._ (2021); Nicholl  _et al._ (2021) or with
the use of accurate numerical relativity (NR) simulations.
Parametric/phenomenological models of BNS mergers cannot reliably capture the
physics of the most extreme stages of the merger, including the merger itself
and the environment immediately following, where the matter is in a highly
dynamical state and the spacetime curvature is the strongest. For a reliable,
first-principles understanding of these stages of BNS mergers, the use of NR
is the only recourse. Numerical simulations allow for the systematic isolation
of different physical phenomena/physics, which provides a powerful tool for
deducing the role of magnetic fields, e.g., Paschalidis  _et al._ (2015); Ruiz
_et al._ (2018, 2016a); Kiuchi  _et al._ (2018); Shibata  _et al._ (2021);
Ciolfi and Kalinani (2020); Mosta  _et al._ (2020); Ciolfi (2020a), equation
of state (EOS) effects, e.g., Abdikamalov  _et al._ (2009); Bauswein  _et al._
(2013); Palenzuela  _et al._ (2015); Sekiguchi  _et al._ (2015); Lehner  _et
al._ (2016); Radice  _et al._ (2018a); Most  _et al._ (2019a); Dietrich  _et
al._ (2018); Vincent  _et al._ (2020); Raithel  _et al._ (2021a), and neutrino
transport, e.g., Sekiguchi  _et al._ (2015); Radice  _et al._ (2018a); Nedora
_et al._ (2019); Foucart  _et al._ (2020).

There are many important physical effects to consider in simulations of BNS
mergers (see Paschalidis (2017); Baiotti and Rezzolla (2017); Duez and
Zlochower (2019); Radice  _et al._ (2020); Bernuzzi (2020); Ciolfi (2020b) for
reviews); such as finite temperature EOSs and magnetic fields. Regarding the
use of realistic EOSs in BNS merger simulations, an important aspect is the
consistent evolution of the electron fraction, as it provides crucial
information for the description of nucleosynthesis processes that take part
during and after merger, and is necessary to fully describe kilonova (KN)
afterglows Kasen  _et al._ (2013); Barnes and Kasen (2013). Additionally, the
thermodynamic properties of the post-merger remnant and ejecta are important
and not well approximated by cold EOSs (see, e.g., Rezzolla  _et al._ (2018);
Raithel  _et al._ (2021b)). Magnetic fields are also a crucial ingredient to
consider in BNS merger simulations. In particular, an area which remains
poorly understood is the effect of magnetic fields on ejecta properties. Large
scale magnetic fields are important in the evolution of the post-merger
remnant and in the kilonova signal associated with BNS mergers, as they are
expected to power relativistic jets Blandford and Znajek (1977); Van Putten
(1999); Bucciantini  _et al._ (2012); Ruiz  _et al._ (2016b); Mosta  _et al._
(2020) and drive relativistic outflow Siegel  _et al._ (2014); Fernandez  _et
al._ (2019); Metzger  _et al._ (2018). The interplay between finite
temperature EOSs and strong magnetic fields remains insufficiently
characterized in BNS merger simulations. For example, realistic EOSs have been
shown to produce lower velocity ejecta when compared to simple analytic EOSs
Wanajo  _et al._ (2014); Sekiguchi  _et al._ (2015); Perego  _et al._ (2017);
Vincent  _et al._ (2020); Bernuzzi (2020) while, on the other hand,
simulations which account for magnetic effects suggest that, if magnetic
fields are strong enough, shocked dynamical ejecta during merger could be
boosted to higher velocities Kiuchi  _et al._ (2018); Radice  _et al._
(2018b); Shibata  _et al._ (2021); Ciolfi and Kalinani (2020); Mosta  _et al._
(2020); Ciolfi (2020a). Most modern simulations of BNS mergers consider only
some of the aforementioned physical effects, and there have only been a
handful of simulations which consider all effects at once Kiuchi  _et al._
(2018); Shibata  _et al._ (2021); Ciolfi and Kalinani (2020); Mosta  _et al._
(2020); Ciolfi (2020a); Sun  _et al._ (2022). Even in the case of more
complete simulations, all treatments of neutrino transport to-date must make
some level of approximation since solving the full 6+1 dimensional neutrino
transport equation at is not possible with current schemes and computational
resources.

There are currently available several codes with support to different levels
of physics. For instance, many state-of-the-art codes with microphysical EOS
compatibility do not implement a constrained-transport treatment of magnetic
fields. While variety in the numerical methods and codes used in the study of
BNS mergers is a strong-point from the perspective of understanding numerical
systematic errors, key differences between codes can lead to a range of
predictions for the relevant observables, depending on the physics included in
simulations. For example, although there is general agreement in the ejecta
properties predicted by different codes Palenzuela  _et al._ (2015); Sekiguchi
_et al._ (2015); Lehner  _et al._ (2016); Shibata  _et al._ (2017); Perego
_et al._ (2017); Bovard  _et al._ (2017); Radice  _et al._ (2018); Nedora  _et
al._ (2019); Vincent  _et al._ (2020), dynamical ejecta masses range over
$${10}^{-4}\text{\,}\mathrm{M}_{\odot}$\lesssim M_{\rm
ej}\lesssim${10}^{-2}\text{\,}\mathrm{M}_{\odot}$$ while speeds range over
$0.1c\lesssim v_{\rm ej}\lesssim 0.3c$, depending on the total binary mass,
EOS, and numerical methods considered. Simulations with even the highest
resolutions can have numerical uncertainties of over $40\%$ Foucart  _et al._
(2016); Radice  _et al._ (2018); Shibata  _et al._ (2017); Bernuzzi  _et al._
(2020). The need for accurate simulations and detailed estimates of the error
budgets is becoming more and more urgent as future GW detectors are expected
to be more sensitive than current systematic errors Pürrer and Haster (2020).
Future investigations of BNS mergers would ideally require a set of catch-all,
accurate, open source GRMHD codes which can reliably simulate BNS systems
while considering as many microphysical phenomena as possible. At the very
least, reliable future predictions using NR codes require detailed
comparisons, including cross-checks of the results obtained using different
codes and of the varied numerical methods used within the codes themselves.

In order to begin addressing each of these needs, we first present a new
implementation of realistic, finite-temperature EOS support in
IllinoisGRMHD111Recently, Werneck  _et al._ (2022) implemented a similar
extension to IllinoisGRMHD. (as extensively discussed below), and then perform
a systematic comparison between different open-source GRMHD codes:
IllinoisGRMHD Etienne  _et al._ (2015), GRHydro Mosta  _et al._ (2014), Spritz
Giacomazzo  _et al._ (2020), and WhiskyTHC Radice and Rezzolla (2012); Radice
_et al._ (2013, 2014). These codes implement the same physics and have largely
overlapping numerical methods, e.g., they are all housed within the
EinsteinToolkit Loffler  _et al._ (2012); Etienne  _et al._ (2021), thereby we
can isolate and compare the impact on the evolution and uncover the existence
of systematic errors arising only from the GRMHD codes. This allows us to test
the robustness of certain results with respect to different computational
choices. In most cases, such choices are about failsafes for when the main
algorithms fail, as in the case of the conservative-to-primitive scheme used
in the artificial atmosphere. In our comparison, we focused on BNS mergers and
found that important quantities disagree across the codes. For example, codes
do not agree on whether the remnant is stable or undergoes collapse very
shortly after merger, which has profound implications for EM observations.
This calls for more detailed studies and comparisons across different codes.

This paper is split in two parts. First, we discuss in detail the formalism
employed by IllinoisGRMHD, the extensions required to reach feature-parity
with the other codes, and the tests we performed. In particular, in Sec. II we
review the basic equations that are relevant to BNS systems, and provide
detail on their numerical treatment within IllinoisGRMHD. In Sec. III we
describe our methods for implementing realistic, finite temperature EOS
capability within IllinoisGRMHD. In Sec. IV we discuss stringent dynamical
tests of the extended version of IllinoisGRMHD which cover a broad range of
scenarios. Readers that are familiar with GRMHD simulations may skip this part
and focus on the results. In the second part, starting from Section IV.5, we
discuss the results from our simulations of BNS mergers.

Our extension to IllinoisGRMHD to allow for the use of realistic, finite-
temperature EOSs is public. The extensions we have made to IllinoisGRMHD are
crucial for understanding the interplay between magnetic fields and the EOS in
the post-merger environment, and is a first step to including neutrino
transport schemes in the code. The current state of IllinoisGRMHD (including
the extensions we describe here) opens up many avenues of investigation
surrounding compact object mergers with strong magnetic fields. For instance,
the new code capabilities will make possible the calculation of
nucleosynthesis rates and help elucidate the role of microphysical, finite
temperature EOSs in BNS mergers with strong magnetic fields, among many other
interesting phenomena.

Throughout the work we use geometrized units, where $G=c=1$, unless otherwise
stated. In addition, in cases where we use logarithmic scales, we assume that
$\log\equiv\log_{10}$, unless otherwise noted. All visualizations and post-
process analyses in this work were carried out with the kuibit software
package Bozzola (2021).

##  II Basic equations

Throughout this work we will predominantly work with numerical codes that
solve the Einstein field equations,

$G_{\mu\nu}=8\pi T_{\mu\nu},$ (1)

(where $G^{\mu\nu}$ and $T^{\mu\nu}$ are the Einstein and stress-energy
tensors, respectively), coupled to the equations of ideal relativistic
(magneto)hydrodynamics. In particular, we focus on the use of IllinoisGRMHD to
solve the equations of ideal relativistic magnetohydrodynamics. IllinoisGRMHD
takes advantage of the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation
Nakamura  _et al._ (1987); Shibata and Nakamura (1995); Baumgarte and Shapiro
(1999) of the 3+1 Arnowitt-Deser-Misner (ADM) formalism, which recasts the
Einstein equations in the form of an initial-value problem (for more details,
see textbooks on the subject, e.g. Alcubierre (2008); Bona and Palenzuela-
Luque (2005); Baumgarte and Shapiro (2010); Shibata (2016)), in which the
spacetime line element is given by

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

where $\alpha$ is the lapse, $\beta^{i}$ is the shift vector,
$\gamma_{\mu\nu}=g_{\mu\nu}+n_{\mu}n_{\nu}$ represents the induced metric on
spacelike hypersurfaces, and $n^{\mu}=(1/\alpha,-\beta^{i}/\alpha)$ is the
future-pointing unit vector orthogonal to each spatial slice. The solution to
Eq. (1) involves the evolution of the magnetohydrodynamic variables that
appear on the right hand sides of the ADM equations. An approach which is well
suited to the 3+1 formalism is the Eulerian (or Valencia) formulation of
relativistic hydrodynamics Font (2008); Rezzolla and Zanotti (2013). In this
formulation, the evolution equations for the relevant fluid variables arise
from several conservation laws, including the continuity equation

$\nabla_{\mu}(\rho_{\rm b}u^{\mu})=0,$ (3)

where $\rho_{\rm b}$ is the rest mass density and $u^{\mu}$ is the fluid four-
velocity, the lepton number conservation (when neutrino effects are ignored)

$u^{\mu}\nabla_{\mu}(Y_{\rm e})=0,$ (4)

which can be rewritten by use of Eq. (3) as

$\nabla_{\mu}(\rho_{b}Y_{\rm e}u^{\mu})=0,$ (5)

where $Y_{\rm e}\equiv n_{\rm e}/n_{\rm b}$ is the electron fraction and
$n_{\rm b}$ ($n_{\rm e}$) is the baryon (electron) number density, the
conservation of stress-energy

$\nabla_{\mu}T^{\mu\nu}=0,$ (6)

and the homogeneous Maxwell’s equations

$\nabla_{\nu}F^{*\mu\nu}=0,$ (7)

where $F^{*\mu\nu}=\frac{1}{2}\varepsilon^{\mu\nu\alpha\beta}F_{\alpha\beta}$
is the dual to the electromagnetic field strength tensor $F_{\alpha\beta}$ and
$\varepsilon^{\mu\nu\alpha\beta}$ is the rank-4 Levi-Civita symbol. The matter
variables are evolved using Eqs. (3)-(7) once they have been cast in flux-
conservative form,

$\partial_{t}\mathbf{C}+\mathbf{\nabla}\cdot\mathbf{F}=\mathbf{S},$ (8)

where $\mathbf{C},\mathbf{F}$, and $\mathbf{S}$ are vectors which contain the
conservative, flux, and source terms, respectively. The vector of primitive
variables $\mathbf{P}$

$\displaystyle\mathbf{P}$ $\displaystyle=\begin{bmatrix}\rho_{\rm b}\\\ P\\\
v^{i}\\\ B^{i}\\\ Y_{\rm e}\end{bmatrix},$ (9)

contains information on the physical state of the fluid, where $P$ is the
fluid pressure, $v^{i}=u^{i}/u^{0}$ is the fluid three-velocity and $B^{i}$
are the spatial components of the magnetic field as measured by a normal
observer. The conservative variables $\mathbf{C}$ are determined in terms of
the primitive variables, the lapse $\alpha$, and the metric as

$\displaystyle\mathbf{C}$ $\displaystyle=\begin{bmatrix}\rho_{*}\\\
\tilde{\tau}\\\ \tilde{S}_{i}\\\ \tilde{B}^{i}\\\ \tilde{Y}_{e}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\alpha\sqrt{\gamma}\rho_{\rm b}u^{0}\\\
\alpha^{2}\sqrt{\gamma}T^{00}-\rho_{*}\\\ (\rho_{*}h+\alpha
u^{0}b^{2})u_{i}-\alpha\sqrt{\gamma}b^{0}b_{i}\\\ \sqrt{\gamma}B^{i}\\\
\rho_{*}Y_{\rm e}\end{bmatrix},$ (10)

where $b^{\mu}=B^{\mu}_{(u)}/\sqrt{4\pi}$ (where $B^{\mu}_{(u)}$ is the
magnetic field as measured by an observer in the fluid rest-frame), $\gamma$
is the determinant of the 3-metric, and $h=1+\epsilon+P/\rho_{0}$ is the
specific enthalpy (where $\epsilon$ is the specific internal energy). We note
that IllinoisGRMHD uses the coordinate three-velocity, $v^{i}\equiv
u^{i}/u^{0}$, unlike many other evolution codes, which also adopt the Valencia
formalism, but use the Eulerian three-velocity $v^{i}_{(n)}=u^{i}/(\alpha
u^{0})+\beta^{i}/\alpha$ Banyuls  _et al._ (1997); Anton  _et al._ (2006),
i.e., the velocity measured by normal observers. Additionally, we note that
IllinoisGRMHD works with an ideal fluid stress-energy tensor of the form

$T_{\mu\nu}=(\rho_{\rm
b}h+b^{2})u_{\mu}u_{\nu}+(P+b^{2}/2)g_{\mu\nu}-b_{\mu}b_{\nu}.$ (11)

The evolution equations for the relevant fluid variables are determined by
using the aforementioned conservation laws and casting them in the form of Eq.
(8). It is useful to pay special attention to the evolution of the magnetic
field, due to how it is treated in IllinoisGRMHD. Specifically, the evolution
equation for the magnetic field is

$\partial_{t}\tilde{B}^{i}+\partial_{j}(v^{j}\tilde{B}^{i}-v^{i}\tilde{B}^{j})=0,$
(12)

where $\tilde{B}^{i}$ is the conservative variable corresponding to the
magnetic field $B^{i}$, as defined in Eq. (10). The use of Eq. (12) to evolve
the magnetic field results in terms that violate the no-monopole constraint
($\mathbf{\nabla}\cdot\mathbf{B}=0)$, which is addressed in IllinoisGRMHD by
instead considering the evolution of the vector potential
$\mathcal{A}_{\mu}=\Phi n_{\mu}+A_{\mu}$ (where $A_{i}$ is the magnetic vector
potential and $\Phi$ is electric scalar potential) and recovering the magnetic
field as $\tilde{B}^{i}=\varepsilon^{ijk}\partial_{j}A_{k}$. The evolution
equation for $A_{i}$ is

$\partial_{t}A_{i}=\varepsilon_{ijk}v^{j}\tilde{B}^{k}-\partial_{i}(\alpha\Phi-\beta^{j}A_{j}).$
(13)

IllinoisGRMHD works in the generalized Lorenz gauge
$\nabla_{\mu}\mathcal{A}^{\mu}=\xi n_{\mu}\mathcal{A}^{\mu}$ Farris  _et al._
(2012), where $\xi$ is chosen such that the Courant-Friedrich-Lewy (CFL)
condition corresponding to it is satisfied at all times given the grid
choices. In the following we provide additional descriptions of the algorithms
employed within IllinoisGRMHD for the implementation of finite temperature EOS
compatibility. We direct the reader to Etienne  _et al._ (2015) for a detailed
description of all of the additional algorithms employed within IllinoisGRMHD.

##  III Implementation of realistic equation of state compatibility within
IllinoisGRMHD

The current open-source version of IllinoisGRMHD solves the equations of GRMHD
by assuming simple, analytic EOSs, such as $\Gamma$-law or piecewise
polytropic EOSs with hybrid thermal effects. These EOSs can only provide a
qualitative understanding of the state of matter during a BNS merger Rezzolla
_et al._ (2018); Raithel  _et al._ (2021a). However, efforts to model
parametrically both the thermal and the cold component of the nuclear EOS in
BNS mergers are under way, see, e.g., Raithel  _et al._ (2021a); Most and
Raithel (2021); Raithel and Paschalidis (2022); Raithel  _et al._ (2022)

The implementation of realistic EOSs within IllinoisGRMHD allows us to
understand, in a more detailed manner, the interplay between the EOS, thermal
effects, and magnetic fields in these systems. Moreover, it is a crucial first
step toward implementing additional important microphysics within
IllinoisGRMHD, such as neutrino transport. The inclusion of realistic EOS
capability within IllinoisGRMHD required two steps: the implementation of the
evolution equation for the electron fraction and the implementation of
algorithms which can perform the non-trivial inversion of the evolved
conservative variables $\mathbf{C}$ to the physical primitive variables
$\mathbf{P}$, which we refer to as conservative-to-primitive inversion. In the
following we discuss the implementation of the evolution equation for the
electron fraction $Y_{\rm e}$ within our extended version of IllinoisGRMHD. We
also discuss the implementation of state-of-the-art routines within
IllinoisGRMHD for conservative-to-primitive inversion and present relevant
tests. For the remainder of this work, we refer to the currently available
version of IllinoisGRMHD as OriginalIllinoisGRMHD (abbreviated as OIL). Our
extended version of IllinoisGRMHD will be referred to as
MicrophysicalIllinoisGRMHD (abbreviated as MIL). In cases where we discuss
features which are common between the two codes, we will refer to them jointly
as the IllinoisGRMHD code. Algorithmically, OIL and MIL are identical, except
for the changes highlighted in this work which are required for realistic EOS
compatibility.

###  III.1 Evolution of the electron fraction

In flux-conservative form, the electron fraction evolution equation is

$\partial_{t}(\tilde{Y}_{\rm e})+\partial_{j}(v^{j}\tilde{Y}_{\rm e})=0.$ (14)

With the inclusion of Eq. (14), the full set of equations solved within MIL is

$\displaystyle\partial_{t}\begin{bmatrix}\rho_{*}\\\ \tilde{Y}_{\rm e}\\\
\tilde{\tau}\\\ \tilde{S}_{i}\\\ A_{i}\end{bmatrix}$
$\displaystyle+\partial_{j}\begin{bmatrix}\rho_{*}v^{j}\\\ \tilde{Y}_{\rm
e}v^{j}\\\ \alpha^{2}\sqrt{\gamma}T^{0j}-\rho_{*}v^{j}\\\
\alpha\sqrt{\gamma}T^{j}_{i}\\\ \alpha\Phi-\beta^{j}A_{j}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}0\\\ 0\\\ s\\\
\frac{1}{2}\alpha\sqrt{\gamma}T^{\alpha\beta}g_{\alpha\beta,i}\\\
\epsilon_{ijk}v^{j}\tilde{B}^{k}\end{bmatrix},$ (15)

where

$\displaystyle
s=\alpha\sqrt{\gamma}[(T^{00}\beta^{i}\beta^{j}+2T^{0i}\beta^{j}+T^{ij})K_{ij}$
(16) $\displaystyle-(T^{00}\beta^{i}+T^{0i})\partial_{i}\alpha],$

and $K_{ij}$ is the extrinsic curvature. The evolution of $\tilde{Y}_{\rm e}$
follows that of the other conservative variables, which begins with the
determination of initial conditions. Presently, we allow for $Y_{\rm e}$ to be
initialized in two possible ways:

  1. 1.

Linear profile: we set $Y_{\rm e}=\Upsilon\rho_{\rm b}$ (where $\Upsilon$ is a
constant that ensures proper dimensionality), such that the electron fraction
profile is linear with respect to $\rho_{\rm b}$, in order to consider a
simple profile where gradients are non-zero throughout the solution grid. This
initialization is unphysical and only useful for testing the advection of
$Y_{\rm e}$ in situations without realistic EOSs, i.e., passive advection of
the electron fraction.

  2. 2.

$\beta$-equilibrium profile: we set $Y_{\rm e}(\rho_{\rm b})$ according to the
conditions for $\beta$-equilibrium in cold neutron star (NS) matter. This
initial profile is suitable for realistic descriptions of isolated stars as
well as for binaries that are initially separated at large enough distances
such that the components are cold. All of the BNS initial data considered in
our tests are built for quasi-equilibrium systems, in which the assumption
that the components are cold, $\beta$-equilibrated stars is well justified. We
use the same assumption of cold $\beta$-equilibrium to initialize all other
hydrodynamical variables at $t=0$.

Once the initial data are specified for the primitive variables at all grid
points, the conservative variables are obtained through the simple algebraic
relations provided in Eq. (10), which provides $\mathbf{C}$ at all grid
points. For the evolution of $\mathbf{C}$, we employ 3 ghost-zones at the
outer boundary of each adaptive mesh refinement (AMR) grid. We fill all buffer
zones at the refinement level boundary with data interpolated from neighboring
rougher or finer levels of refinement using standard prolongation or
restriction operators, respectively. Filling in the buffer zones in this
manner results in $\mathbf{C}$ being either prolongated or restricted after
calculation Etienne  _et al._ (2015). To ensure consistency between
$\mathbf{C}$ and $\mathbf{P}$, the primitives are recovered from the
conservatives using a root-finding algorithm at $t=0$, which we discuss in
Sec. III.2. Next, we evaluate the flux term $\mathbf{F}$ in preparation for
the next time step. To this end, we must reconstruct the primitives between
grid points (i.e, at cell interfaces). IllinoisGRMHD employs the piecewise
parabolic method (PPM) Colella and Woodward (1984) for primitives
reconstruction. Reconstruction is used to evaluate the primitives on the left
and right interfaces of all grid points, $\mathbf{P}_{L,R}$, in all
directions. These interface values are then used to calculate the
corresponding conservative variables at cell interfaces, which are in turn
used to calculate the flux term in Eq. (14) using a second-order, finite-
volume, high-resolution shock-capturing (HRSC) scheme. The handling of fluxes
at grid interfaces $\mathbf{F}_{L,R}$ requires a solution to a Riemann
problem. IllinoisGRMHD employs the standard Harten-Lax-vanLeer (HLL) Harten
_et al._ (1983) approximate Riemann solver, where for a given direction the
electron fraction flux is given by

$F^{\rm HLL}(Y_{\rm
e})=\dfrac{c^{-}F_{R}+c^{+}F_{L}-c^{+}c^{-}(\tilde{Y}_{{\rm
e},R}-\tilde{Y}_{{\rm e},L})}{c^{+}+c^{-}},$ (17)

where $c^{\pm}=\pm\text{max}\left(0,c^{R}_{\pm},c^{L}_{\pm}\right)$ and
$c^{L,R}_{\pm}$ are the maximum ($+$) and minimum ($-$) characteristic speeds
at the left (L) and right (R) cell interfaces (see Etienne  _et al._ (2010,
2015) for further algorithmic details). The derivatives of the fluxes are then
determined and summed independently for each direction. For instance, the flux
along the $x$-direction takes the form

${(\partial_{x}F^{x})}_{ijk}=\dfrac{F^{HLL,x}_{i+\frac{1}{2}jk}(Y_{\rm
e})-F^{HLL,x}_{i-\frac{1}{2}jk}(Y_{\rm e})}{\Delta x}.$ (18)

The fluxes along the $y$\- and $z$-directions take a similar form, but we
instead consider finite differencing along the $j$ and $k$ indices,
respectively.

The evolution equation for $Y_{\rm e}$ does not include source terms in the
absence of neutrinos, so the right-hand-side of Eq. (18) is then passed to the
method of lines (MoL) thorn, which integrates the conservative variable
$\tilde{Y}_{\rm e}$ forward in time. At this point, the updated conservative
variables would be known at all grid points except the outer boundary. The
next step is to recover the primitive variables given these evolved
conservative variables (see Sec. III.2). After the primitives have been
recovered, they are checked for physicality and marginally modified if they
are outside of their physical ranges Etienne  _et al._ (2012). For example, in
the case of the electron fraction, we check that

$Y_{\rm e,\text{lower}}\leq Y_{\rm e}\leq Y_{\rm e,upper},$ (19)

where $Y_{\rm e,\text{lower}}$ ($Y_{\rm e,\text{upper}}$) corresponds to the
lowest (highest) value for $Y_{\rm e}$ available in an EOS table. Next, outer
boundary conditions are placed on the recovered primitives to fill the
necessary 3 ghost-zones in each direction. We apply zero-derivative outflow
outer boundary conditions as described in Etienne  _et al._ (2015). Up to this
point, the primitives $\mathbf{P}$ are known at the new time step on all grid
points. The final step we take is to recompute the conservatives on all grid
points using Eq. (10) for consistency between $\mathbf{P}$ and $\mathbf{C}$,
and the evolution algorithm is allowed to proceed.

###  III.2 Conservative-to-primitive solvers

Figure 1: Top panel: Relative error in the recovery of primitives for a
selection of the solvers originally implemented in Siegel  _et al._ (2018) and
ported to theEinsteinToolkit within the ConservativeToPrimitive thorn (labeled
C2P). We test the LS220 EOS O’Connor and Ott (2010); Lattimer and Swesty
(1991) with three of the available solvers. For all tests, we fix a subset of
the primitives and provide random initial data for the remainder.
Specifically, we fix $P_{\rm mag}/P=10^{-3}$, $W=2$, and $Y_{\rm e}=0.1$ while
scanning the physically allowed range for $T$ and $\rho_{\rm b}$. Bottom
panel: Same as the top panel, but using the original implementation of each
solver in the code packaged with Siegel  _et al._ (2018) (labeled S).

At every step of the evolution an inversion from the evolved conservative
variables $\mathbf{C}$ to the physical primitive variables $\mathbf{P}$ is
required to know the state of the fluid. Eq. (10) presents a system of non-
linear algebraic expressions which can be solved for 9 relevant fluid
variables ($\rho_{\rm b}$, $v^{i}$, $B^{i}$, $Y_{\rm e}$, and either $h$, or
$\epsilon$). These 9 variables, along with an EOS, in turn provide all of the
information required to determine $\mathbf{P}$, along with other fluid
variables of interest. For example, a solution to Eq. (10) can provide the 5
main variables ($\rho_{0}$, $v_{i}$, $\epsilon$) and, trivially, $B_{i}$ and
$Y_{\rm e}$. We can then determine the remaining variable, $P$ with the use of
an EOS and incidentally obtain information on other physical variables such as
the temperature $T$ and specific entropy $s_{\rm b}$. As we require a solution
for 5 _main_ variables in order to determine $\mathbf{P}$, the primitives
inversion problem is fundamentally a non-trivial 5D problem that cannot be
solved analytically. 5D schemes which solve Eq. (10) were originally
implemented in early MHD codes such as HARM Gammie  _et al._ (2003). However,
these schemes were eventually found to be inefficient and inaccurate, which
led to the development of methods which solve for only two auxiliary variables
Noble  _et al._ (2006) and thereby reduce the dimensionality of the problem to
2D. The dimensionality of the problem can be further reduced to 1D; modern 1D
algorithms which provide reliable and efficient solutions have been developed
Newman and Hamlin (2014); Palenzuela  _et al._ (2015) and are widely used in
GRMHD codes.

In order to consider strongly magnetized systems which include realistic
descriptions of the dense matter EOS, we implement state-of-the-art
conservative-to-primitive solvers within MIL. Our implementation includes
porting the solvers discussed in Siegel  _et al._ (2018) for use in the
EinsteinToolkit, packaged within a new thorn ConservativeToPrimitive which can
interface with MIL and its associated thorns, but also works as a standalone
thorn which can be used with other GRMHD codes that operate within the Cactus
infrastructure. We focus on a subset of the solvers implemented in Siegel  _et
al._ (2018), due to their reliability, speed, and algorithmic similarity to
the original solvers used in OIL (see Kastaun  _et al._ (2021) for another
possible approach). In particular we focus on the 2D solver of Noble  _et al._
(2006) (which we label Noble), the 1D solver of Palenzuela  _et al._ (2015)
(which we label Palenzuela), and the 1D solver of Newman and Hamlin (2014)
(which we label Newman). We note that the Noble solver uses the same algorithm
as the solvers in the original version of IllinoisGRMHD and that other GRMHD
codes with realistic EOS capability rely on the Palenzuela algorithm
Giacomazzo  _et al._ (2020); Most  _et al._ (2019b); Kalinani  _et al._
(2022). We refer the reader to Siegel  _et al._ (2018) for a review of the
algorithms used in each solver within ConservativeToPrimitive111We point out
that Siegel  _et al._ (2018) appears to have typographical errors in the
algorithm description for the Newman solver, when compared to the original
paper of Newman and Hamlin (2014). In particular, there is a difference in the
calculation of the auxiliary variable $\mathcal{M}^{2}$ (see Eq. (47) of
Siegel  _et al._ (2018) as compared to Eq. (4.7) of Newman and Hamlin (2014)),
where Newman and Hamlin (2014) correctly calculates it as
$\mathcal{M}^{2}=m_{i}m^{i}=\tilde{S}_{i}\tilde{S}^{i}$, where $\tilde{S}_{i}$
is the conservative variable associated with the momentum density which
appears in the left-hand-side of Eq. (15). Eq. (47) of Siegel  _et al._ (2018)
inconsistently calculates this variable as
$\mathcal{M}^{2}={(B^{i}v_{i})}^{2}/\sqrt{\gamma}$, where $\gamma$ is the
determinant of the 3-metric $\gamma_{ij}$, $B^{i}$ is the magnetic field, and
$v_{i}$ is the fluid 3-velocity. We also note that the first term of Eq. (55)
in Siegel  _et al._ (2018) misses a factor of the auxiliary variable $a$, when
compared to the analogous Eq. (5.11) in Newman and Hamlin (2014). We point out
that despite these typographical errors, the numerical implementation of the
Newman solver within the code of Siegel  _et al._ (2018) is consistent with
the correct algorithmic steps presented in Newman and Hamlin (2014)..

We employ a set of preliminary tests to confirm that our port of these solvers
to the EinsteinToolkit behaves as intended. In particular, we test the
aforementioned solvers using the same tests as Siegel  _et al._ (2018) where a
set of primitives $\mathbf{P}$ are initialized, randomly perturbed, used to
calculate a set of conservatives $\mathbf{C}$, and finally recovered into a
new set $\mathbf{P}$’. The primitives recovery for each solver is then
assessed by considering the relative error between the original set
$\mathbf{P}$ and the recovered set $\mathbf{P}^{\prime}$ (along with other
diagnostics including the number of interpolation calls to the EOS table and
number of algorithm iterations). In these tests, a subset of the primitives is
varied over the physically allowed range while holding others constant. We
focus on the case of the LS220 realistic, finite temperature, tabulated EOS
O’Connor and Ott (2010); Lattimer and Swesty (1991) and employ tests where we
prescribe the ratio of magnetic to fluid pressure $P_{\rm
mag}/P=b^{2}/(2P)=0.001$, the Lorentz factor $W=2$, and the electron fraction
$Y_{\rm e}=0.1$, while scanning over the allowed range of rest mass density
$\rho_{\rm b}$ and temperature $T$ for this EOS. For each value of $\rho_{\rm
b}$ and $T$ that we consider, we determine the relative error between the
original and recovered set of primitives. We present the results for this set
of tests in the top panel of Fig. 1. As a comparison, we also show the results
of these tests using the implementation of the solvers within the code of
Siegel  _et al._ (2018) in the bottom panel of Fig. 1.

For all tests, we find very good agreement between the original solvers and
our implementation of them within the EinsteinToolkit. Points on Fig. 1 which
appear in white indicate the failure of a given solver to recover a consistent
set of primitives. Points at which the solvers fail which are surrounded by
successful recovery are usually due to the random perturbations we induce on
$\mathbf{P}$ before recovery, which occurs for both the original code of
Siegel  _et al._ (2018) and our implementation within ConservativeToPrimitive
(see the top left and lower right panels of Fig. 1). However, there are
regions of variable space where solvers consistently fail regardless of the
initial perturbations (in both our implementation of the solvers within
ConservativeToPrimitive and the original implementation in Siegel  _et al._
(2018)). The most reliable solver across variable space is the Newman solver
although, as indicated by the rightmost panels of Fig. 1, it may fail for cold
matter ($T\lesssim$10\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$) in the high
density regime $\rho_{\rm
b}\approx${10}^{14}\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$$, which is
a highly relevant area of variable space for BNS mergers. It is important to
note that regions where a given solver may fail are often well covered by at
least one of the other solvers. Therefore, in practice the best use of these
solvers should employ a hybrid approach, where a solver which fails in a given
region of variable space is substituted by a solver which provides reliable
solutions in that same region. In the following we discuss the choice of
solvers we employ moving forward.

###  III.3 Solution algorithm

The preliminary tests presented in Fig. 1 suggest that some solvers are better
suited for systems involving neutron stars with realistic EOSs over others.
Nevertheless, we note that the use of _any_ solver is still likely to produce
robust solutions. Therefore, the chief property that we use in determining the
primary solver to utilize is the performance speed and agreement with the
results of the OIL code. For realistic EOSs, the bottleneck in limiting solver
speed is typically the number of times the solver must interpolate the EOS.
Among the solvers in ConservativeToPrimitive, the Noble solver requires the
fewest such EOS calls, resulting in the most efficient solutions Siegel  _et
al._ (2018). Nevertheless, the Palenzuela and Newman solvers require a similar
number of EOS calls to Noble, and the Palenzuela solver has the added benefit
of strong agreement with OIL and of reliable coverage across densities and
temperatures, as suggested by the leftmost panels of Fig. 1. Due to the
reliability and agreement with OIL, and because the performance is not much
worse than other solvers, we consider the Palenzuela solver as the optimal
choice for our primary solver. In the remainder of this work, we adopt the
following solution algorithm for conservative-to-primitive inversion in MIL

  1. 1.

We employ the Palenzuela solver in our first attempt at a solution. We find
that the Palenzuela solver is robust for both analytic and microphysical EOSs
and comparably efficient to other solvers.

  2. 2.

If the Palenzuela solver fails, as it is disposed to for highly relativistic
matter, we try to recover a solution with either the Newman or Noble solvers.
The Noble solver is robust for highly magnetized flows and highly efficient,
but begins to fail at high Lorentz factors. The Newman solver is reliable for
both highly magnetized and highly relativistic flows.

  3. 3.

If all of the aforementioned solvers fail we impose a fix of the velocities,
as done in the OIL code Etienne  _et al._ (2015), which assumes a cold EOS
($T=$0.1\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$); we refer to this fix as
the ‘Font’ fix, which has been generalized from its original form in OIL for
generic EOSs.

  4. 4.

Finally, if all other recovery attempts have failed, it is likely that the
fluid element is at low densities, so we set the recovered primitives to
values which are sensible for atmospheric conditions (cold temperature $T_{\rm
atm}=$0.1\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$, low density $\rho_{\rm
b,atm}\leq 10^{-7}\rho_{\rm b,max}$, fixed electron fraction $Y_{\rm
e,atm}=0.25$, and a pressure which is consistent with the EOS $P_{\rm
atm}=P(\rho_{\rm b,atm},T_{\rm atm},Y_{\rm e,atm})$). We note that the other
open-source codes we consider, with the exception of GRHydro which uses the
same treatment as MIL, only fix atmospheric values of $\rho_{\rm b}$ and $T$,
but not $Y_{\rm e}$. We emphasize, however, that the values of the primitive
variables in the atmosphere are generally dynamically unimportant, and are set
for numerical stability.

We find that the above solution algorithm leads to a robust and efficient
evolution across the physically allowed range of variables for a finite
temperature, realistic EOS. We note that the order of solvers used is
completely customizable at run-time and can be changed dynamically for any
simulation, and that any ordering of the Palenzuela, Noble, and Newman solvers
results in a reliable solution. The order of solvers chosen here is optimal
when considering efficiency and reliability for finite temperature, tabulated
EOSs. In the remainder of this work, we utilize the MIL code with the above
solution algorithm, unless otherwise noted.

##  IV Dynamical tests of MIL and code comparisons

###  IV.1 Summary of tests, initial data, and grid hierarchies

Table 1: Grid hierarchy and simulation settings for the tests presented in Sec. IV, in the case of the lowest-resolution grids employed. For each case we list the test name, the finest-level grid resolution $dx_{\rm fin}$, the number of refinement levels (RL), the number of grid points per NS radius (GP), and the half-side length of each nested grid (HSL) listed from finest to coarsest. We also list the evolution codes (Evol. code), Riemann solver (Riemann sol.), reconstruction scheme (Recon. scheme), and conservative-to-primitive error tolerance (C-to-P tol.) used in each simulation. In terms of the Riemann solver, we use either the Harten-Lax-vanLeer (HLL) Harten  _et al._ (1983) or Harten-Lax-vanLeer-Einfeldt (HLLE) Einfeldt (1988), depending on the algorithm available in each code.  Test | $dx_{\rm fin}$ | RL | GP | HSL ($M_{\odot}$) | Evol. code | Riemann sol. | Recon. scheme | C-to-P tol.
---|---|---|---|---|---|---|---|---
G2TOV | 0.25 | 1 | 32 | $(10.0)$ | MIL | HLL | PPM | $10^{-8}$
|  |  |  |  | GRHydro | HLLE | PPM | $10^{-8}$
$\rm LSTOV$ | 0.111 | 3 | 64 | $(8.9,17.8,35.6)$ | MIL | HLL | PPM | $10^{-8}$
|  |  |  |  | Spritz | HLLE | PPM | $10^{-8}$
G2BNS | 0.155 | 7 | 40 | $(13.0,17.9,26.1,60.0,120.0,240.0,396.0)$ | OIL | HLL | PPM | $10^{-8}$
|  |  |  |  | MIL | HLL | PPM | $10^{-8}$
|  |  |  |  | Spritz | HLLE | PPM | $10^{-8}$
|  |  |  |  | GRHydro | HLLE | PPM | $10^{-8}$
|  |  |  |  | WhiskyTHC | HLLE | MP5 | $10^{-8}$
LSBNS | 0.14 | 7 | 50 | $(8.7,17.4,34.8,69.6,139.2,278.4,556.8)$ | MIL | HLL | PPM | $10^{-8}$
|  |  |  |  | GRHydro | HLLE | PPM | $10^{-8}$
|  |  |  |  | WhiskyTHC | HLLE | MP5 | $10^{-8}$

In this section we focus on code comparisons between MIL and other open-source
codes GR(M)HD codes. As tests of MIL, we employ a number of evolutions for
which the system behavior is well understood. Where relevant, we compare the
results of the MIL code to other publicly available codes with similar
capabilities. In particular we compare with the GRHydro Mosta  _et al._
(2014), Spritz Giacomazzo  _et al._ (2020), WhiskyTHC Radice and Rezzolla
(2012); Radice  _et al._ (2013, 2014), and OIL Etienne  _et al._ (2015) codes.

We present the results of these tests as follows:

  1. 1.

G2TOV: We first discuss tests which employ a TOV star built with a $\Gamma=2$
polytrope, which we label test G2TOV. For test G2TOV we compare the solution
to results obtained with the GRHydro code. This test was designed to mirror
those considered in Etienne  _et al._ (2015).

  2. 2.

LSTOV: Next, we discuss tests where we consider a magnetized TOV star
constructed using the finite temperature, tabulated LS220 EOS O’Connor and Ott
(2010); Lattimer and Swesty (1991), which we label test LSTOV. For this test
we compare our solutions to those obtained with the publicly available Spritz
code, which uses an implementation of the Palenzuela solver for finite
temperature EOSs.

  3. 3.

G2BNS: Next, we test our code in the case of a BNS merger, assuming a
$\Gamma=2$ polytropic EOS. We compare the solution against the results of the
GRHydro, Spritz, WhiskyTHC, and OIL codes. We consider both magnetized and
unmagnetized systems, and refer to this test as G2BNS.

  4. 4.

LSBNS: Finally, we consider a BNS built using the LS220 EOS. Along with the
MIL code, we simulate this system with the GRHydro, and WhiskyTHC codes. The
binary configuration is built such that there is a transient hypermassive
neutron star (HMNS) post-merger remnant, which allows us to compare properties
of the merger remnant between codes. We also consider evolution with the
Spritz code, but we do not include those results, because the public version
of the code is still under development for cases with tabulated EOS support in
the context of binary neutron star mergers. We label this test as LSBNS.

The use of MIL and ConservativeToPrimitive requires the initialization of
three new required variables when evolving finite temperature, realistic EOSs.
Namely, we must initialize the specific internal energy $\epsilon$,
temperature $T$, and electron fraction $Y_{\rm e}$ such that they are
consistent with the initial data. For tests which employ a cold polytropic EOS
we initialize the specific internal energy as

$\epsilon=\dfrac{P}{(\Gamma-1)\rho_{0}}.$ (20)

For polytropic EOSs, the initialization of $T$ and $Y_{\rm e}$ are independent
to that of $\epsilon$ and do not play a role in the fluid evolution, so we
simply initialize $Y_{\rm e}$ based on the first type of profile listed in
Sec. III.1 and fix $T=$0.1\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$, although
the temperature is a passive variable in such tests. For tests which employ
realistic, finite temperature EOSs, we are interested in equilibrium initial
data, which requires cold, uniform temperatures and $\beta$-equilibrium. To
construct cold, $\beta$-equilibrium initial data, we build barotropic tables
that provide $P(\rho_{\rm b})$. Due to the fact that ConservativeToPrimitive
makes use of EOS tables in the format of the StellarCollapse repository
O’Connor and Ott (2010) via the EOS driver thorn EOS_Omni, we construct cold,
$\beta$-equilibrium tables required for initial data using the StellarCollapse
tables. The StellarCollapse tables provide several fluid variables (including
the pressure $P$, specific internal energy $\epsilon$, and constituent
chemical potentials $\mu_{i}$, where $i=p,n,e$ is an index over the particle
species for protons, neutrons, and electrons, respectively) at discrete values
of the rest mass density $\rho_{\rm b}$, electron fraction $Y_{\rm e}$, and
temperature $T$. Each of the fluid variables may be interpolated to arbitrary
values of the triplet $(\rho_{\rm b},Y_{\rm e},T)$. To extract the cold,
$\beta$-equilibrated, barotropic functions required for initial data
($P(\rho_{\rm b})$, and $\epsilon(\rho_{\rm b})$), we take the following
steps:

  1. 1.

We fix the value of $T=T_{\rm
cold}=$0.1\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$. The lowest possible
temperature values in these tables are typically
$T=0.01-$0.05\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$, but entries near the
table limits are not finely sampled, which often leads to interpolation errors
in the primitives recovery process. As such, we employ a finite temperature
which is still cold compared to the Fermi energy but avoids the table
boundaries.

  2. 2.

For a given value of $\rho_{\rm b}=\rho_{\rm b,curr}$, we scan the available
range of $Y_{\rm e}$ and evaluate the chemical potentials relevant for beta
equilibrium ($\mu_{i}(Y_{\rm e},\rho_{\rm b,curr},T_{\rm cold})$).

  3. 3.

We employ root-finding to locate the value of $Y_{\rm e}$ which corresponds to
$\beta$-equilibrium, such that

$\mu_{n}-\mu_{p}-\mu_{e}=0.$ (21)

At this value of $Y_{\rm e}=Y_{\rm e,\beta}$, we record the pressure and
specific internal energy. We continue the algorithm from step 2 above,
scanning the available range of $\rho_{\rm b}$ and thereby constructing
tabulated functions for the barotropic pressure $P$ and specific internal
energy $\epsilon$ as functions of the rest mass density $\rho_{\rm b}$ for
$\beta$-equilibrium matter.

In Tab. 1 we present the grid hierarchy corresponding to each dynamical test
we consider. For each case we list the test name, the finest-level grid
resolution $dx_{\rm fin}$, the number of refinement levels (RL), the number of
grid-points per NS radius (GP), the half-side length of each nested grid (HSL)
listed from finest to coarsest, and the evolution codes used. In all cases we
evolve the spacetime using the McLachlan spacetime evolution code within the
EinsteinToolkit (via the ML_BSSN thorn) Brown  _et al._ (2009); Reisswig  _et
al._ (2011), unless otherwise noted. ML_BSSN solves the Einstein equations
within the BSSN formulation of the ADM 3+1 formalism. Our gauge conditions
consist of the “1+log” slicing condition for the lapse Bona  _et al._ (1995)
and the “Gamma-driver” condition for the shift Alcubierre  _et al._ (2003).
Integration in time is carried out using a fourth-order accurate Runge-Kutta
(RK4) scheme, using the MoL thorn, with a Courant factor of 0.5, unless
otherwise stated.

###  IV.2 Diagnostics

We employ several diagnostics to assess the quality of our simulations. In
cases where several resolution simulations are considered and some solution is
known, we monitor convergence of the solution by calculating the convergence
order $n_{\rm con}$ as

$n_{\rm
con}=\dfrac{\log\left(\dfrac{L_{i}-L_{0}}{L_{j}-L_{0}}\right)}{\log\left(\Delta
x_{i}/\Delta x_{j}\right)},$ (22)

where $L_{0}$ is the value that the solution should approach in the continuum
limit, $i$ is a label for a grid of a given resolution, $L_{i}$ corresponds to
the solution on grid $i$, and $\Delta x_{i}$ corresponds to the grid
resolution for grid $i$. For all cases involving fluid variables we typically
consider the reference solution $L_{0}$ in Eq. (22) to be the initial data
itself. The use of such a reference solution can prove troublesome in cases
with unreliable or faulty initial data. In several of the cases we consider,
we find a lack of convergence toward the continuum solution even at the start
of the simulation, suggesting non-convergent initial data. In such cases we
consider the self-convergence of a given solution, which is demonstrated when

$\displaystyle(L_{\rm CR}-L_{\rm HR})\left[{\left(\dfrac{\Delta x_{\rm
MR}}{\Delta x_{\rm HR}}\right)}^{n}-1\right]=$ (23) $\displaystyle(L_{\rm
MR}-L_{\rm HR})\left[{\left(\dfrac{\Delta x_{\rm CR}}{\Delta x_{\rm
HR}}\right)}^{n}-1\right],$

where $n$ is the expected convergence order, $L$ is the solution on a given
grid, $\Delta x$ is the grid resolution for a given grid, and the labels CR,
MR, and HR correspond to the canonical, medium, and high resolution grids,
respectively. We refer to quantities scaled according to Eq. (23) using the
symbol $c[L]$, for quantity $L$.

We monitor the maximum rest mass density $\rho_{\rm b,max}$ and minimum of the
lapse function $\alpha_{\rm min}$ to qualitatively assess stability and
collapse. We also consider spacetime quantities, such as the Hamiltonian
constraint violations, in order to assess convergence of the solution. We
consider 2D snapshots of several fluid quantities, including the rest mass
density $\rho_{\rm b}$, temperature $T$, and different components of the
magnetic field strength $B_{i}$, where relevant. When considering rotational
properties of the fluid, we compute the angular velocity in the equatorial
plane as

$\Omega=\dfrac{Xv^{y}-Yv^{x}}{\varpi^{2}},$ (24)

where $\varpi\equiv\sqrt{X^{2}+Y^{2}}$, $X\equiv(x-x_{\rm com})$,
$Y\equiv(y-y_{\rm com})$, $x$ and $y$ are the Cartesian coordinates, and
$x_{\rm com}$ and $y_{\rm com}$ correspond to the coordinates of the Newtonian
center-of-mass.

Figure 2: Left panel: Results of test ${\rm G2TOV}$ in the case of evolution
with the MIL code. In the top panel we show the relative change in central
rest mass density compared to the value at $t=0$, where $\Delta\rho_{\rm
b,c}\equiv\rho_{\rm b,c}(t)-\rho_{\rm b,c}(0)$. We depict results for grid
resolutions that employ 32, 40, and 50 grid points per NS radius using solid
magenta, dashed blue, and dotted green lines, respectively. In the bottom
panel we show the convergence rate, calculated using Eq. (22), for the
quantity depicted in the top panel in cases where we compare the
canonical/medium and medium/high resolution results using solid red and dashed
blue lines, respectively. Right panel: same as left panels but for the GRHydro
code.

In cases where we consider gravitational radiation, we extract GWs using the
Newman-Penrose formalism Newman and Penrose (1962); Penrose (1963), with focus
on the $s=-2$ spin-weighted spherical harmonic decompositions of the Newman-
Penrose scalar $\Psi_{4}$. The coefficients of the spin-weighted decomposition
are labeled $\Psi_{4}^{l,m}$, where $l$ and $m$ are the degree and order for
the spherical harmonics. We extract $\Psi_{4}^{l,m}$ from the numerical
solution at fixed concentric spheres with increasing coordinate radii $r_{\rm
ex}=\eta~{}$\text{\,}\mathrm{k}\mathrm{m}$$, where $\eta$ takes on several
discrete values $70\lesssim\eta\lesssim 450$. When extracting GWs, we ensure
to use the values of $\Psi_{4}^{l,m}$ which are in the wave zone. Where
relevant, we compute the GW strain $h$ from

$\Psi_{4}=\ddot{h}_{+}-i\ddot{h}_{\times}.$ (25)

using the fixed-frequency integration (FFI) method Reisswig  _et al._ (2011).
We identify the time of merger $t_{\rm mer}$ as the time when the GW amplitude
reaches an absolute maximum during the merger.

###  IV.3 G2TOV: Test using a TOV star with polytropic EOS

We first consider the case of TOV initial data constructed using a $\Gamma=2$
polytropic EOS. We construct initial data using the TOVSolver thorn within the
EinsteinToolkit with a central rest mass density of $\rho_{\rm b,c}=0.129$ and
polytropic constant of $\kappa=1$. We consider a static spacetime and do not
employ mesh refinement (i.e., the test is carried out in the Cowling
approximation). This is the only test we consider under the Cowling
approximation. The half-side length of the solution grid extends to
$1.25R_{\rm NS}$, where $R_{\rm NS}\approx$8\text{\,}\mathrm{M}_{\odot}$$ is
the NS radius. Our canonical, medium, and high resolution grids employ grid
spacings such that there are 32, 40, and 50 grid points covering the NS
radius, respectively. We evolve with the MIL code and with the GRHydro code as
a comparison. In the case of the GRHydro evolution we employ the Harten-Lax-
vanLeer-Einfeldt (HLLE) Einfeldt (1988) Riemann solver and PPM reconstruction
for maximum algorithmic overlap with MIL.

In the top panel of Fig. 2 we show the evolution of the relative change in the
central rest mass density, where $\Delta\rho_{\rm b,c}\equiv\rho_{\rm
b,c}(t)-\rho_{\rm b,c}(0)$ over a timescale of $20\,t_{\rm dyn}$, where

$t_{\rm dyn}=1/\sqrt{\rho_{\rm b,c}(0)}$ (26)

is the dynamical time. We show results for the low, medium, and high
resolution grids, using solid magenta, dashed blue, and dotted green lines,
respectively. We find that the evolution of a polytropic TOV star proceeds as
expected, with oscillations in the rest mass density not exceeding $\sim
0.2\%$ and decreasing over time as the model settles. In the lower panel of
Fig. 2 we also depict the convergence order of $\Delta\rho_{\rm b,c}/\rho_{\rm
b,c}(0)$, with comparisons of the canonical/medium and medium/high resolution
results depicted with the solid red and dashed blue lines, respectively. As we
increase the grid resolution we find that the quantity $\Delta\rho_{\rm b,c}$
converges to zero at the expected second-order rate. In the right panel of
Fig. 2 we show the results of the ${\rm G2TOV}$ test in the case of evolution
with the GRHydro code. We find behavior in the GRHydro code which is
consistent with that of the MIL code. When evolving with GRHydro we find that
oscillations in the rest mass density are of larger amplitude when compared to
the MIL code, but the solution nonetheless converges at the expected rate.

###  IV.4 LSTOV: Test using a cold, magnetized TOV star with realistic EOS

Figure 3: Left panel: Results of test ${\rm LSTOV}$ in the case of evolution
with the MIL code. We depict the central rest mass density $\rho_{\rm b,c}$
scaled by the value at $t=0$ (top panel), central temperature $T_{\rm c}$
(lower panel, solid lines), and maximum temperature $T_{\rm max}$ (lower
panel, dashed lines) as functions of time (scaled by the dynamical time) for a
magnetized TOV star built with the LS220 EOS using the Cook code Cook  _et
al._ (1994a). We show the results corresponding to evolution with the MIL
(magenta lines) and the Spritz (blue lines) codes. Right panel: Same as the
left panel, but in the case where initial data was built using the Lorene code
Gourgoulhon  _et al._ (2016).

In this set of tests we consider a cold, magnetized TOV star built using a
realistic EOS. We construct initial data with several codes as a way to
understand the systematic error introduced at the level of the initial data.
In particular, we use the Cook code Cook  _et al._ (1994a), the RNSID code
Stergioulas and Friedman (1995); Nozawa  _et al._ (1998) packaged within the
EinsteinToolkit, and the Lorene code Gourgoulhon  _et al._ (2016). In the
following we focus on initial data constructed with the Cook and Lorene codes,
as the RNSID code produces similar results. We construct initial data for a
cold, barotropic, $\beta$-equilibrated EOS obtained using the LS220 EOS. The
initial configuration has a central energy (rest mass) density of
$0.8\text{\times}{10}^{15}\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$
$($0.74\text{\times}{10}^{15}\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$)$,
ADM (rest) mass of
$$1.4\text{\,}\mathrm{M}_{\odot}$($1.55\text{\,}\mathrm{M}_{\odot}$)$, and
circumferential radius of $R_{\rm NS}=$12.69\text{\,}\mathrm{k}\mathrm{m}$$.
We employ 3 levels of fixed mesh refinement, with the innermost, finest
refinement level extending to $1.5R_{\rm NS}$, and every subsequent level
extending to twice the distance of the adjacent finer level. For all tests we
employ a finest-level grid resolution of $dx_{\rm fin}=R_{\rm NS}/64$ (which
we refer to as canonical resolution), such that there are 64 grid points
covering the NS radius. We superimpose a purely poloidal magnetic on non-
magnetized initial data which is confined to the interior of the NS, using the
Seed_Magnetic_Fields_BNS code within the EinsteinToolkit. The field structure
approximately corresponds to that generated by a current loop Paschalidis  _et
al._ (2013), with the toroidal component of the vector potential taking the
form

$A_{\phi}=\dfrac{\pi
r_{0}^{2}I_{0}\varpi^{2}}{{(r_{0}^{2}+r^{2})}^{3/2}}\left(1+\dfrac{15r_{0}^{2}(r_{0}^{2}+\varpi^{2})}{8{(r_{0}^{2}+r^{2})}^{2}}\right),$
(27)

where $r_{0}$ is the current loop radius, $I_{0}$ is the loop current,
$r^{2}=\varpi+z^{2}$, $\varpi={(x-x_{\rm NS})}^{2}+{(y-y_{\rm NS})}^{2}$, and
$x_{\rm NS}$ and $y_{\rm NS}$ are the initial coordinates of the NS center-of-
mass. This choice of magnetic field exhibits a $1/r^{3}$ fall-off outside the
NS, as expected for a dipole field, but as mentioned above we do not extend
the field to the exterior; the Seed_Magnetic_Fields_BNS code employs a cut-off
pressure $P_{\rm cut}$ below which the vector potential is set to zero, which
we set to $P_{\rm cut}\approx 0.15P_{\rm max}(0)$ (where $P_{\rm max}(0)$ is
the maximum fluid pressure at the start of the simulation). The initial field
strength corresponds to a ratio of magnetic to fluid pressure of $P_{\rm
mag}/P\approx 0.001$.

As a comparison, we also show the evolution of this model with the Spritz code
Giacomazzo  _et al._ (2020), which uses an implementation of the Palenzuela
solver for finite temperature, tabulated EOSs. We employ the Spritz code with
PPM reconstruction and the HLLE Riemann solver for maximum algorithmic overlap
with IllinoisGRMHD. In Fig. 3 we show the evolution of the central rest mass
density $\rho_{\rm b,c}$ (top panels), maximum temperature $T_{\rm max}$
(bottom panels using dashed lines), and central temperature $T_{c}$ (bottom
panels using solid lines) as functions of time. In the left panel of Fig. 3 we
focus on cases where initial data was constructed with the Cook code. We
depict results corresponding to MIL and Spritz using magenta and blue lines,
respectively. We find that the evolution of the rest mass density is
qualitatively the same between the MIL the Spritz codes, exhibiting
oscillations of at most $\sim 1\%$. The configuration shows an initial
relatively large oscillation in the rest mass density of $\gtrsim 1\%$ which
peaks at $t\approx 1-2\,t_{\rm dyn}$. This initial large oscillation is likely
due to the creation of a cold barotropic EOS from the finite temperature LS220
EOS table as previously discussed, which requires interpolation and root-
finding routines (we note that similar oscillations happens in test cases with
unmagnetized initial configurations, which suggests that the superimposed
magnetic field is likely not the culprit behind the oscillations). It is
possible that the use of a cold EOS created from a finite temperature EOS
results in initial data with percent-level numerical errors, which is
consistent with the oscillations of amplitude $\sim 1.01\rho_{\rm b,c}(0)$
shown in the top left panel of Fig. 3. We note that initial data produced with
the Cook code uses a compactified radial grid (true for both RNSID and LORENE)
in spherical polar coordinates, and as such interpolation of the initial data
onto the Cartesian grids of the Cactus infrastructure results in additional
numerical errors in the initial configuration. As the initial configuration
settles to the nearest equilibrium, it begins to exhibit smaller oscillations
in the rest mass density of $\sim 0.5\%$. These oscillations decay over a
timescale of $t\approx 20\,t_{\rm dyn}$. In the bottom left panel of Fig. 3,
we show the evolution of the maximum and central temperatures using dashed and
solid lines, respectively, for initial data constructed with the Cook code.
The initial configuration is cold for all solvers tested, with a central
(maximum) temperature of $T_{\rm
c}\approx$0.1\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$ ($T_{\rm
max}\approx$0.5\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$). We note that
initial data is constructed for an isothermal star at
$T=$0.1\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$. However, once the initial
data for $\rho_{b}$, $\epsilon$, $Y_{\rm e}$, and $P$ are interpolated onto
the solution grid, we recover the temperature for self-consistency among the
fluid variables, which leads to errors in the temperature profile at $t=0$,
near the stellar surface. Early in the simulations, regardless of the
evolution code used, the model develops a warm atmosphere, which becomes
increasingly warm as the star undergoes small oscillations. The central
configuration remains cold as the temperature of the atmosphere increases and
saturates to $T_{\rm max}\approx$8\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$.

Figure 4: Left panel: Results of test ${\rm LSTOV}$ in the case of evolution
with the MIL code. We depict the convergence of the L2 norm of the Hamiltonian
constraints $||H||$ in the case of a magnetized TOV star built using the LS220
EOS with the Lorene code, evolved with the MIL code and Palenzuela solver. The
top panel depicts $||H||$ for the low, medium, and high resolution grids using
solid red, dashed green, and dotted blue lines, respectively. The bottom panel
depicts $||H||$ scaled according to second-order convergence while comparing
the canonical and medium (CR,MR) and medium and high (MR,HR) resolution
results using red and blue lines, respectively. Right panel: Same as the left
panel but in the case of evolution with the Spritz code.

To test the role of initial data on the evolution of the model, we consider
the same initial configuration but instead built with the Lorene and RNSID
codes. We consider the same central rest mass density and grid hierarchy as in
the case considering initial data constructed with the Cook code. In the right
panel of Fig. 3 we present results of the evolution with the MIL and Spritz
codes using magenta and blue lines, respectively, with focus on the case of
initial data built with the Lorene code. In all cases considered, we observe
initial oscillations of the rest mass density which do not exceed $2\%$ and
subsequent oscillations which are at the level of $\sim 1\%$. We generally
find that lower temperatures when using initial data constructed with Lorene,
in both the bulk of the star and in the atmosphere. We find that evolution
with the Spritz code results in a central region which remains cold
($T\leq$0.1\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$, shown using the blue
solid line in the lower right panel of Fig. 3) throughout the simulation,
whereas the MIL code produces slightly warmer central regions of up to $T_{\rm
c}=$0.3\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$. We note that in the MIL code
we impose a lower bound on the temperature of
$0.1\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$ to avoid EOS interpolations near
the table limits. Such a low-temperature floor is not imposed in the Spritz
code, which may be the reason for the slight difference in $T_{\rm c}$ between
the two codes. A central temperature of $T_{\rm
c}=$0.3\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$ is still significantly cold
from a nuclear perspective and tests of TOV star evolution without imposing
the lower bound on $T$ results in a temperature profile which is consistent
with the Spritz evolutions. The evolution of the rest mass density and
temperature presented in Fig. 3 is typical of stable TOV stars. The fact that
a realistic TOV star remains stable and cold over a timescale of $t\approx
20\,t_{\rm dyn}$ and the qualitatively similar behavior between the MIL and
Spritz codes indicates that our implementation of the solvers within
ConservativeToPrimitive is reliable. Tests with RNSID initial data produce
very similar results to those presented in Fig. 4.

We consider the models constructed using the Cook code for a resolution study.
We evolve the Cook initial data in two additional higher resolution
simulations with both the MIL and Spritz codes. For medium (high) resolution
tests, we employ grid resolutions of $dx_{\rm med}=dx_{\rm can}/1.25$
($dx_{\rm hi}=dx_{\rm can}/1.25^{2}$) such that there are 80 (100) grid points
covering the NS radius, where $dx_{\rm can}=dx_{\rm fin}$ as listed in Tab. 1
for the LSTOV test. In the top panel of Fig. 4 we show the L2 norm of the
Hamiltonian constraints $||H||$ for the canonical (solid red lines), medium
(dashed green lines) and high (dotted blue lines) resolution grids, in cases
where we evolve with the MIL (left panel) and Spritz (right panel) codes. We
find that $||H||$ is generally small and converges toward zero (the expected
value in the continuum limit) for higher resolution grids. We also note that
the evolution of $||H||$ is in agreement between the MIL and Spritz codes.
However, we find that $||H||$ converges toward zero at a rate which is slower
than expected, for both evolution codes considered. Specifically, we find a
linear convergence rate of $||H||$ toward zero. We note that the expected
convergence rate of the evolutions we consider is second-order, and that
second-order convergence is demonstrated for polytropic initial data (see Fig.
2). It is likely that the lower order convergence observed for tests employing
realistic EOSs is due to errors in the initial data. In particular, the
initial data codes we consider make use of compactified spherical polar
coordinate grids, which are not the same as the evolution grids. Numerical
errors are therefore inevitable when interpolating solutions from their
original compactified coordinates to the Cartesian grids employed within the
EinsteinToolkit. We point out that the evolution of magnetized TOV stars
presented in this section is consistent with stable evolution, and coincides
with the results presented in Giacomazzo  _et al._ (2020) for Spritz. As there
are likely errors at the level of the initial data when using realistic EOSs
we instead consider the self-convergence of the solution. In the lower panels
of Fig. 4 we show $||H||$ scaled according to second order convergence (i.e.,
$n=2$ in Eq. (23)) for test LSTOV while comparing the canonical and high
resolution (solid red lines) and medium and high resolution solutions (dashed
blue lines). We again show results for evolution with the MIL and Spritz codes
using the left and right panels, respectively. The lower panels of Fig. 4
demonstrate that the solution exhibits self-convergence at the expected order.
We also consider the same test of convergence in cases with initial
configurations built with the Lorene and RNSID codes and find similar results
to those presented in Fig. 4. We point out that the construction of TOV
initial data is algorithmically very similar between the Cook, RNSID, and
Lorene codes, as they all make use of compactified spherical coordinate grids.
As such, interpolation errors that results from a change in the coordinate
system will arise in all tested cases.

###  IV.5 G2BNS: Test of a BNS with a Gamma-law EOS

Figure 5: Left panel: Maximum rest mass density $\rho_{\rm b}$, scaled by the
value at $t=0$, for an unmagnetized BNS assuming a $\Gamma=2$ polytropic EOS.
We depict the results corresponding to the GRHydro (solid purple line), OIL
(dashed magenta line), MIL with and without the use of the solvers within
ConservativeToPrimitive (dash-dotted golden and dotted orange lines,
respectively), Spritz (dash-dash-dotted line), and WhiskyTHC (blue dash-dash-
dot-dotted line) codes. Right panel: Same format as the left panel, but
instead depicting the physical time produced by each simulation per hour of
run time, across 224 Intel Haswell V3 cores. We use identical initial data,
grid structures, conservative-to-primitive error tolerance, and computational
resources in all cases, as well as similar flux treatments and reconstruction
schemes. The tests presented here are as close as possible to a one-to-one
comparison between codes that we could achieve.

As an additional test of the solvers implemented within
ConservativeToPrimitive, we consider the evolution of a BNS system built
assuming a polytropic EOS. We consider this test because the OIL code can
reliably simulate such a system and because it provides a stringent test of
the MIL code new algorithms in a highly relativistic scenario. We assume a
polytropic constant of $\kappa=123.61$ and adiabatic index of $\Gamma=2$. We
construct initial data for an equal mass BNS with an orbital separation of
$45\text{\,}\mathrm{k}\mathrm{m}$ using the Lorene code; the ADM mass and
radius of each binary component are $M_{\rm
NS}=$1.4\text{\,}\mathrm{M}_{\odot}$$ and $R_{\rm
NS}=$10.36\text{\,}\mathrm{k}\mathrm{m}$$, respectively. During the evolution,
we assume a $\Gamma$-law EOS with adiabatic index $\Gamma=2$. In this test we
employ AMR, with two sets of nested boxes centered on each NS and each
employing 7 levels of refinement. We employ an additional set of nested boxes
with 7 levels of refinement centered at the origin of the solution to reliably
resolve the BNS merger remnant. The outer grid boundary is located at $R_{\rm
out}=$582.12\text{\,}\mathrm{k}\mathrm{m}$\approx 38.4R_{\rm NS}$ and the
half-side length of the finest level around each star extends to $\sim
1.25R_{\rm NS}$ with each subsequent coarser level extending to $\eta R_{\rm
NS}$ where $\eta\in(1.7,2.5,5.8,11.6,23.2)$. The finest-level grid resolution
is set to $dx_{\rm fin}\approx R_{\rm NS}/40$, such that we initially resolve
each binary component with at least 40 grid points along the radius. For test
G2BNS, we consider both an unmagnetized, and magnetized scenario, because not
all of the evolution codes considered can handle a constrained-transport
evolution of the magnetic field. In cases where a magnetic field is also
considered, we seed purely poloidal fields which are confined to the interior
of each NS, with vector potential of the form Eq. (27), such that the initial
maximum field strength is $B_{\rm
max}\approx$2.3\text{\times}{10}^{15}\text{\,}\mathrm{G}$$.

We simulate the system with the use of MIL and track the evolution of the rest
mass density $\rho_{\rm b}$, fluid 3-velocity $v^{i}$, and the magnetic field
$B^{i}$, when relevant. We approximate the merger time $t_{\rm mer}$ by
considering the time when the rest mass density reaches a maximum during the
evolution. We compare to several other open-source codes in the literature,
including the OIL, GRHydro, WhiskyTHC, and Spritz codes. Similar to other
tests, we employ all evolution codes considered with the HLL/HLLE Riemann
solver. For the IllinoisGRMHD, GRHydro, and Spritz cases we use PPM
reconstruction. For the WhiskyTHC case we use MP5 reconstruction SURESH199783.
We also note that WhiskyTHC uses an internal rescaling of the rest mass
density, and as such the initial conditions between WhiskyTHC and other cases
is not identical. However, this internal rescaling of the rest mass density
results in a maximal relative difference of less than $0.4\%$ between cases.

####  IV.5.1 Unmagnetized BNS

In cases where we do not consider the evolution of the magnetic field, we are
free to compare the results provided by all of the open-source GRHD codes
considered, as they are all able to handle an unmagnetized BNS system built
with a polytropic EOS. As such, we compare the results of simulations using
the MIL, OIL, GRHydro, WhiskyTHC, and Spritz codes. We emphasize that the
numerical grids, initial data, reconstruction scheme, conservative-to-
primitive schemes, and flux treatment for all of these simulations are chosen
for maximum overlap, with the latter three likely having differences in
implementation between codes. The comparison of results presented in this
section was the closest scenario to a one-to-one comparison between open-
source numerical codes that we could obtain.

Figure 6: Top panel: Equatorial snapshots of the rest mass density $\rho_{\rm
b}$ (scaled by the maximum value at $t=0$) for test G2BNS at
$t\approx$7.3\text{\,}\mathrm{m}\mathrm{s}$$. From left to right we show
results for the OIL code, MIL code with Palenzuela solver, and Spritz code.
The rightmost panel depicts a snapshot for an evolution using the MIL code
wherein we use an alternative calculation of the specific enthalpy $h$ (see
Eqs. (28) and (29)). Bottom panel: Equatorial snapshots of the $z-$component
of the magnetic field $B_{z}$ (scaled by the maximum value at $t=0$) for test
${\rm G2BNS_{mag}}$ at $t\approx$3.54\text{\,}\mathrm{m}\mathrm{s}$$ for the
same set of codes as the top panel. We point out common features that develop
in the magnetic field structure using black arrows.

In the left panel of Fig. 5 we show the maximum rest mass density in the
simulations for the unmagnetized case of test G2BNS. We show results for the
GRHydro (solid purple line), OIL (dashed magenta line), MIL with and
without222We consider an evolution with the MIL code that does not employ the
solvers within ConservativeToPrimitive, and instead uses the original solvers
within the OIL code. We label this test as ${\tt MIL_{\rm noC2P}}$ in Fig. 5.
The fact that the ${\tt MIL_{\rm noC2P}}$ evolution proceeds identically to
the OIL case ensures that other changes made to the code (e.g., the changes
that allow for passive $Y_{\rm e}$ advection) do not affect the dynamics. the
use of the solvers within ConservativeToPrimitive (dash-dotted golden and
dotted orange lines, respectively), Spritz (dash-dash-dotted line), and
WhiskyTHC (blue dash-dash-dot-dotted line) codes. We find that during the
inspiral, all codes behave qualitatively the same, with the maximum rest mass
density exhibiting oscillations of $\lesssim 2\%$ for the first $t\approx
5\text{ ms}$ of the inspiral. We find that the merger time for these
simulations all fall within $9.81\pm 0.1\text{ ms}$, resulting in a difference
of at most $2\%$. On the other hand, we find significant deviations of the
simulation results after the merger. Each code considered results in different
size oscillations of the rest mass density for the first $10\text{ ms}$ after
merger. The GRHydro code results in the largest post-merger oscillations in
$\rho_{\rm b,max}$, and is the only case where $\rho_{\rm b,max}$ shows a
tendency to grow. In all cases besides GRHydro, the evolution of $\rho_{\rm
b,max}$ is qualitatively similar, and the post-merger value of $\rho_{\rm
b,max}$ oscillates between $0.9\rho_{\rm b,max}(0)-\rho_{\rm b,max}(0)$. We
find qualitatively similar GW emission in all cases, but the frequency
corresponding to peak GW emission shows a spread of at most $\sim 15\%$
between cases. Given that all codes use similar flux treatments (we note that
in all other tested codes, the closest Riemann solver to the HLL solver
adopted in IllinoisGRMHD is the positivity-preserving HLLE solver Einfeldt
(1988), which is what we adopt in evolutions with Spritz, GRHydro, and
WhiskyTHC) and reconstruction schemes (in the case of WhiskyTHC, we use MP5
while in all other cases PPM is used), along with identical numerical grid,
initial data, conservative-to-primitive error tolerance, and spacetime
evolution scheme, it is unclear why the GRHydro code deviates so strikingly
from all other codes, and why there is a general lack of quantitative
agreement in the post-merger evolution between most codes. Key differences
between codes may arise in the implementation of conservative-to-primitive
algorithms, which may be the chief reason for discrepancies between codes. We
leave the investigation of algorithmic implementations within open-source GRHD
codes, and a more detailed code comparison, to future work.

In the right panel of Fig. 5 we show a measure of the performance for each
code. Specifically, we show the amount of physical time produced per hour for
each simulation, using the same amount of computational resources
(specifically, 224 Intel Haswell V3 cores). We find that, for this case study,
the WhiskyTHC code performs the fastest (producing around 0.07 ms/hour
throughout the inspiral and in the post-merger phase, respectively), while the
Spritz code performs the slowest (producing around 0.03 ms/hour during the
inspiral and 0.035 ms/hour in the post-merger phase, respectively). All other
codes perform similarly, with 0.05 ms/hour during the inspiral and 0.07
ms/hour in the post-merger. We note that the hierarchy suggested by the right
panel of Fig. 5 is not a final verdict on the performance of each code, and
only demonstrates the performance for the particular case considered. We find
that changes to the specific settings governing the conservative-to-primitive
inversion and atmosphere control can significantly impact the performance of
the codes considered. For example, we find that the conservative-to-primitive
routines within Spritz often perform as reliably as in other codes with the
use of fewer algorithmic iterations, which serves to speed up the code.

####  IV.5.2 Magnetized BNS

In cases where we _do_ consider a magnetic field for test G2BNS, we compare
the evolution provided by the MIL code to that provided by the OIL and Spritz
codes (note that the Spritz code uses an implementation of the Noble solver
for polytropic EOSs). In the top panel of Fig. 6 we show equatorial snapshots
of the rest mass density for test G2BNS at time
$t\approx$7.3\text{\,}\mathrm{m}\mathrm{s}$$, corresponding to first contact
between the binary components. From left to right in Fig. 6 we show snapshots
for evolution with the OIL, MIL, and Spritz codes. In the rightmost panel of
Fig. 6 we also show a snapshot corresponding to the MIL code in the case where
we use an alternative set of reconstructed variables to calculate the specific
enthalpy $h$ (labeled $\texttt{MIL}_{\tt b}$), as will be discussed in further
detail below. We find that all codes considered produce very similar results
as far as the merger time and features during inspiral are concerned. We find
that the merger times agree to within 2%, specifically $t_{\rm
mer}=47.85t_{\rm dyn}$, $t_{\rm mer}=47.85t_{\rm dyn}$, and $t_{\rm
mer}=49.04t_{\rm dyn}$ for the OIL, MIL and Spritz cases, respectively. At a
time close to and preceding merger, pictured in the top panel of Fig. 6, the
main difference between the results of each code can be observed in the region
between the binary components, likely due to how low density, near-atmosphere
regions are treated in the conservative-to-primitive solvers implemented in
each code.

In the lower panel of Fig. 6 we show snapshots of the $z-$component of the
magnetic field during the inspiral at time
$t\approx$3.54\text{\,}\mathrm{m}\mathrm{s}$$. Moving from left to right, we
again show results for the OIL, MIL, and Spritz codes; the rightmost panel
corresponds to the MIL code with primitive reconstruction that agrees with the
OIL code (see Eq. (28)), which we label $\texttt{MIL}_{\tt b}$. Specifically,
in OIL the calculation of specific enthalpy in the flux computation uses an
analytic expression for $\epsilon$, based on a $\Gamma-$law EOS. In the case
where the adiabatic index for the thermal component of the EOS is $\Gamma_{\rm
th}=2$, the specific enthalpy becomes

$h=1+\epsilon(P_{\rm r})+\dfrac{P_{\rm r}}{\rho_{\rm b,r}}=1+\dfrac{2P_{\rm
r}}{\rho_{\rm b,r}},$ (28)

where $P_{\rm r}$ and $\rho_{\rm b,r}$ are the reconstructed pressure and rest
mass density, respectively. On the other hand, the solution within MIL and the
Spritz codes calculate these specific enthalpies based on the reconstructed
values of $\epsilon$

$h=1+\epsilon_{\rm r}+\dfrac{P(\epsilon_{\rm r},Y_{\rm e,r},\rho_{\rm
b,r})}{\rho_{\rm b,r}},$ (29)

where $P(\epsilon_{\rm r},Y_{\rm e,r},\rho_{\rm b,r})$ is the pressure for a
generic EOS based on the reconstructed values of the specific internal energy
$\epsilon_{\rm r}$, electron fraction $Y_{\rm e,r}$, and rest mass density
$\rho_{\rm b,r}$. In other words, the OIL code uses the reconstructed pressure
to calculate $h$, whereas the MIL and Spritz codes use the reconstructed
specific internal energy. The calculation of $h$ based on Eq. (29) is required
for self-consistency among the fluid variables in the case of generic EOSs,
but in the case of analytic EOSs, it is preferable to proceed with calculation
of $h$ based on Eq. (28) for accuracy. The differences in use of reconstructed
variables plays a role in the preservation of the magnetic field structure. We
find that the OIL code tends to advect the magnetic field structure accurately
for longer times when compared to the other tested codes. The MIL and Spritz
evolutions all produce similar features in the magnetic field near the low
density parts of the star. In particular, solutions employing generic solvers,
which are compatible with realistic EOSs, produce ‘kinks’ in the magnetic
field structure, as highlighted by the arrows in the lower panel of Fig. 6. We
note that, despite the artifacts seen in the magnetic field structure for
solvers compatible with realistic EOSs, the magnetic field is generally
properly advected during inspiral in all cases. Moreover, the sizes of the
artifacts decrease (and the time at which they appear increases) with the use
of higher resolutions grids. The agreement of MIL with the well-tested Spritz
and OIL codes suggests that the fluid and magnetic field solutions these codes
provide are reliable for highly dynamical systems. For the sake of consistency
with the OIL code, we adopt the same methods within MIL for calculating
specific enthalpies in cases where we consider analytic EOSs. We find that
this results in a field structure that is conserved throughout more of the
inspiral, as illustrated by the disappearance of the ‘kinks’ in the field
structure in the rightmost snapshot in Fig. 6, wherein we use the same method
for calculating specific enthalpies as in the OIL code. In other words, moving
forward, within the MIL code we calculate $h$ based on Eqs. (28) and Eq. (29)
for analytic and finite temperature EOSs, respectively.

###  IV.6 LSBNS: Tests using a BNS with finite temperature equations of state

Figure 7: Maximum rest mass density $\rho_{\rm b,max}$ as a function of time
for the inspiral and merger of a BNS using the LS220 EOS (test LSBNS). We
depict the evolution of $\rho_{\rm b,max}$ in the case of the MIL, WhiskyTHC,
and GRHydro codes using solid purple, dashed orange, and dotted green lines,
respectively. Figure 8: GW signals corresponding to test LSBNS. We depict the
characteristic strain at a distance of
$10\text{\,}\mathrm{M}\mathrm{p}\mathrm{c}$ from the source in the case of the
MIL, WhiskyTHC, and GRHydro codes using solid purple, dashed orange, and
dotted green lines, respectively. For a better comparison, we use the same
duration signals in all cases, corresponding to times ranging from
$t\approx$10\text{\,}\mathrm{m}\mathrm{s}$$ to
$t\approx$20\text{\,}\mathrm{m}\mathrm{s}$$ (with signal duration $\tau_{\rm
signal}\approx$10\text{\,}\mathrm{m}\mathrm{s}$$), which allows us to focus on
the post-merger signals. The vertical lines of the same color and pattern mark
the peak frequency for each case.

MIL

WhiskyTHC

GRHydro

Figure 9: Top row: Equatorial snapshots of the rest mass density (left half)
and temperature (right half) for test LSBNS, in the case of evolution with the
MIL code. From left to right we depict snapshots at times
$t=$5.32\text{\,}\mathrm{m}\mathrm{s}$$,
$t=$8.37\text{\,}\mathrm{m}\mathrm{s}$$,
$t=$12.70\text{\,}\mathrm{m}\mathrm{s}$$, and
$t=$18.88\text{\,}\mathrm{m}\mathrm{s}$$, respectively. The top (bottom)
colorbar corresponds to the rest mass density (temperature). Middle row: same
as top tow but in the case of evolution with the WhiskyTHC code. We outline
regions above which the rest mass density takes on values of $\rho_{\rm
b}=${10}^{13}\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$$ and $\rho_{\rm
b}=${10}^{14}\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$$ using dashed and
solid black contours, respectively. Bottom row: same as top row but in the
case of evolution with the GRHydro code.

In this section we consider the evolution of BNS systems built using a finite-
temperature EOS. Specifically, we consider the LS220 EOS and construct initial
data for an unmagnetized, equal-mass BNS with ADM mass $M_{\rm
ADM}=$2.7\text{\,}\mathrm{M}_{\odot}$$, and initial separation of
$45\text{\,}\mathrm{k}\mathrm{m}$. Our evolution grid extends to approximately
$830\text{\,}\mathrm{k}\mathrm{m}$ at the outer boundary and consists of 3
additional sets of 6 nested cubes (corresponding to 7 levels of refinement);
we center one set of cubes at the origin of the solution grid, and the other
two sets of cubes are used to track the center of each neutron star during the
inspiral. The grid resolution of each nested cube is half that of the adjacent
larger cube, and the finest-level grid resolution is set such that each NS is
resolved with at least $50$ grid points per radius during the inspiral;
generally, the grid resolution of a given refinement level is
$dx_{l}=2^{(7-l)}R_{\rm NS}/50$, where $l=7(1)$ labels the finest (coarsest)
resolution grid and $R_{\rm NS}\approx$10.31\text{\,}\mathrm{k}\mathrm{m}$$.
We note that although these grid resolutions are lower than the standard used
for production-level runs, such grid resolutions are routinely considered in
the context of BNS mergers as low-resolution simulations which are useful in
convergence studies Sekiguchi  _et al._ (2012); Hayashi  _et al._ (2022). We
employ these grids because we were predominantly interested in comparisons of
different code features. In the context of BNS mergers, we leave comparisons
using high-resolution simulations and a full convergence study to future work,
since each of the codes we consider has been tested extensively for
convergence. The half-side length of the finest-resolution grid is set to
$r_{\rm HSL,7}=1.25R_{\rm NS}\approx$12.95\text{\,}\mathrm{k}\mathrm{m}$$;
generally, the half-side length corresponding to a given refinement level is
$r_{{\rm HSL},l}=1.25R_{\rm NS}2^{(7-l)}$. In the case of the MIL code, we use
the HLL Riemann solver. In all other cases we use the HLLE Riemann solver. We
use PPM for primitive reconstruction for all evolution codes considered except
WhiskyTHC, in which case we use MP5.

As a way to compare different open-source GRHD codes currently available, we
consider the evolution of this system with use of the MIL, WhiskyTHC, and
GRHydro codes. We present the maximum rest mass density for each evolution in
Fig. 7, using a purple solid line, orange dashed line, and green dotted line
for the case of MIL, WhiskyTHC, and GRHydro, respectively.

We identify the merger time $t_{\rm mer}$ as the time where the $l=2$, $m=2$
mode of the dominant polarization of the GW strain peaks; for reference, we
show these GWs in Fig. 8. We find merger times of $t_{\rm
mer}=$13.57\text{\,}\mathrm{m}\mathrm{s}$$, $t_{\rm
mer}=$11.72\text{\,}\mathrm{m}\mathrm{s}$$, and $t_{\rm
mer}=$12.14\text{\,}\mathrm{m}\mathrm{s}$$ in the case of the MIL, WhiskyTHC,
and GRHydro evolutions, respectively, which results in an approximately
$13.6\%$ relative difference, at most, in merger times between cases. In Fig.
8 we show the characteristic strain Moore  _et al._ (2015) extracted from the
simulations in test LSBNS, assuming a source at a distance of
$10\text{\,}\mathrm{M}\mathrm{p}\mathrm{c}$. We find very similar inspiral GW
signals from all codes considered, with relative differences in the phase of
the gravitational wave strains of at most $\sim 2\%$. We find that the largest
discrepancy between GW signals arises in the immediate post-merger
environment. The post-merger environment is when the matter is in a highly
turbulent and dynamical state; in this stage of the merger that the relative
difference in the phase of the GW strain can reach as high $\sim 20\%$ (when
comparing the MIL and WhiskyTHC cases), but is typically closer to $\sim
10\%$. After this highly dynamical state, as the remnant begins to settle, all
codes again produce very similar GW signals. The frequency corresponding to
peak post-merger GW emission $f_{\rm peak}^{2,2}$ is similar between all
codes, with the MIL, WhiskyTHC, and GRHydro evolutions resulting in $f_{\rm
peak}^{2,2}=$3.11\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$$, $f_{\rm
peak}^{2,2}=$3.06\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$$, and $f_{\rm
peak}^{2,2}=$3.13\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$$.

In the top, middle, and lower panels of Fig. 9 we show snapshots of the rest
mass density and temperature in the equatorial plane at key points during the
merger for the case of MIL, WhiskyTHC, and GRHydro, respectively. For each
evolution code considered we depict 4 snapshots, ordered by time from left to
right at different stages of the merger corresponding to times
$t\in(4.25,8.37,12.70,18.88)\text{ ms}$. In the left-most panels in Fig. 9 we
depict the system at a time during the inspiral. We generally find that the
binary components remain cold during the inspiral, with typical temperatures
inside the stars of $T\lesssim 1-$2\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$.
Temperatures of $\mathcal{O}($1\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$)$
within the stars are expected during the inspiral Raithel  _et al._ (2021b),
and may be due to the growth of numerical error seeded at the level of initial
data or interpolation errors from the EOS tables. We find that the stellar
atmospheres become significantly warm as the stars inspiral, reaching
temperatures of $T_{\rm atmo}\approx
5-$10\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$ before even the first orbit is
complete (note that we employ a density cutoff to define the atmosphere, such
that all matter with below $\rho_{\rm
b,atmo}\approx$6.17\text{\times}{10}^{07}\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$\approx
10^{-7}\rho_{\rm b,max}(0)$ is set to the atmosphere as detailed in Sec.
III.2). This temperature profile (cold temperatures of
$T\approx$1\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$ within the stars and
higher temperatures $T_{\rm
atmo}\approx$10\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$ in the atmospheres)
is maintained during the entire inspiral up to merger, as depicted by the
leftmost panels in Fig. 9. We find that, during the inspiral, the warmest
temperatures develop in the case of the GRHydro and WhiskyTHC evolutions,
reaching as high as $T\approx$5\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$ in
the bulk of the star and $T_{\rm
atmo}\approx$15\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$ in the atmosphere.
These cases also exhibit warm, low-density clouds which develop around the
binary components throughout the inspiral. However, the atmosphere is
dynamically unimportant.

Figure 10: Dominant coefficients of the spin-weighted decomposition of the
Newman-Penrose scalar $\Psi_{4}^{2,m}$ (scaled with the extraction radius
$r_{\rm ex}$ and ADM mass $M$) as a function of retarded time $t_{\rm ret}$
(shifted by the time of merger $t_{\rm mer}$) for the $l=2$, $m=1$ (dashed
lines) and $l=2$, $m=2$ (solid lines) GW modes for test LSBNS. In the left,
center, and right panels we show results for evolution with MIL (purple
lines), WhiskyTHC (orange lines), and GRHydro (green lines), respectively

In the second-from-left column in Fig. 9 we depict the systems at a time near
merger. At this point in the evolution the first significant shock heating
happens as the cores of the two stars collide. In all cases, temperatures in
the shear layer between the two stellar cores can climb to $T\approx
30-$40\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$. As the two cores continue to
orbit and merge, the warm shocked material is redistributed toward the outside
of the merger remnant, resulting in a warm envelope of
$T\approx$10\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$ surrounding the merged
core which can reach temperatures of
$T\approx$40\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$. We depict such states
in the third column of Fig. 9, which roughly correspond to a few milliseconds
after merger. Finally, in the rightmost panels of Fig. 9 we show the state of
the merger remnants at a time corresponding to approximately
$7-$9\text{\,}\mathrm{m}\mathrm{s}$$ after merger. As the merger remnant
settles, the final configurations approach a warm central object with
temperatures $T\approx$10\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$. We find
that in all cases this central object is surrounded by a distribution of
matter which resembles a ring of hot material (with temperatures ranging
$T\approx 40-$50\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$), which becomes
cooler with increased radial distance from the central configuration (dropping
to temperatures $T\approx 10-$20\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$).

Figure 11: Time-averaged, azimuthally-averaged angular velocity as a function
of radial distance on the equatorial plane for the case where we use the MIL
(purple solid line) and GRHydro (dashed green line) codes in test LSBNS. We
time-average over a period of approximately $10\text{\,}\mathrm{m}\mathrm{s}$,
ranging from simulation times of $t\approx$20\text{\,}\mathrm{m}\mathrm{s}$$
to $t\approx$30\text{\,}\mathrm{m}\mathrm{s}$$. The purple and green contours
depicts the range of values for the instantaneous angular velocity in the time
window we consider for the MIL and GRHydro evolutions, respectively. The
profiles indicate that the remnant is differentially rotating by the end of
simulation.

For a clear view of remnant properties, we consider rest mass density cutoffs
which allow us to identify regions of the remnant. We identify the remnant
itself and its core as regions in the post-merger environment above $\rho_{\rm
b}\geq${10}^{13}\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$$ and
$\rho_{\rm b}\geq${10}^{14}\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$$,
respectively Kiuchi  _et al._ (2018). We highlight the remnant and its core
using dashed and solid black lines in Fig. 9, respectively. We find that, in
all cases, immediately following the merger a dense core is formed near the
origin which typically extends to coordinate distance of
$\sim$15\text{\,}\mathrm{k}\mathrm{m}$$ from the origin and is surrounded by
lower density material which extends to coordinate distance
$\sim$20\text{\,}\mathrm{k}\mathrm{m}$$. This HMNS remnant is surrounded by a
low-density, disk-like structure.

The two key factors that determine the final fate of a BNS merger remnant are
the remnant rest mass and the maximum rest mass of a uniformly rotating NS
(i.e., the ‘supramassive’ limit mass) $M_{\rm supra}$ for the given EOS. We
note that, for hadronic EOSs like the ones considered in this work, the value
of $M_{\rm supra}$ is closely related to the value of the maximum rest mass of
a non-rotating star $M_{\rm TOV}$ and is expected to be about $20\%$ larger
$M_{\rm supra}\approx 1.2M_{\rm TOV}$ Cook  _et al._ (1994a, b); Lasota  _et
al._ (1996); Morrison  _et al._ (2004); Breu and Rezzolla (2016) (although
generally this mass increase can reach up to about $35\%$ Bozzola  _et al._
(2019)); in the case of the LS220 EOS, $M_{\rm
supra}=$2.83\text{\,}\mathrm{M}_{\odot}$$, and the total system rest mass in
test LSBNS is $M_{\rm b}=$2.89\text{\,}\mathrm{M}_{\odot}$$. One of the main
sources of support against gravitational collapse in BNS merger remnants comes
from the large amounts of differential rotation in the system (thermal support
can also contribute significantly against gravitational collapse Raithel  _et
al._ (2021b)). Differential rotation can greatly increase the maximum mass a
star can support Ansorg  _et al._ (2009); Bozzola  _et al._ (2018); Espino and
Paschalidis (2019), but, as differential rotation is removed from the system
(for example, via braking by magnetic fields Shapiro (2000); Shibata  _et al._
(2021)), the remnant will approach a uniformly rotating NS. If the total
remnant mass is below $M_{\rm supra}$ for a given EOS, we may expect a long-
lived NS remnant Margalit and Metzger (2017). On the other hand, for systems
with $M_{\rm tot}>M_{\rm supra}$, we may expect collapse. In all cases, we
observe a transient HMNS remnant, as expected for this total binary mass and
EOS. A key difference in the post-merger evolution for the simulations we
consider is in the survival time of the HMNS remnant. In the MIL case, we find
that the remnant does not collapse to a black hole (BH) by the end of
evolution ($t_{\rm fin}\sim 40\text{ ms}$). On the contrary, the remnant
collapses after $t_{\rm coll}\approx$21.6\text{\,}\mathrm{m}\mathrm{s}$$ and
$t_{\rm coll}\approx$25.2\text{\,}\mathrm{m}\mathrm{s}$$ in the case of the
WhiskyTHC and GRHydro evolutions, respectively (demonstrated by the sharp rise
in the rest mass density in Fig. 7 and corroborated by the collapse of the
lapse function $\alpha$ in both cases).

The remnant survival time $\tau_{\rm remnant}$ is an important quantity which
significantly affects several observables associated with BNS mergers. For
instance, neutrino irradiation from a metastable HMNS may affect ejecta
properties, which would in turn affect the KN signal Metzger and Fernandez
(2014); Margalit and Metzger (2017); Holmbeck  _et al._ (2021). Additionally,
$\tau_{\rm remnant}$ is expected to affect the nucleosynthesis in accretion
disk outflows, with shorter HMNS lifetimes resulting in relatively low
lanthanide and actinide abundances Lippuner  _et al._ (2017) due to the
shorter neutrino irradiation times Metzger and Fernandez (2014); Martin  _et
al._ (2015a, b); Miller  _et al._ (2020). The significant differences in
remnant survival times $\tau_{\rm remnant}$ exhibited by each code suggests
that additional sutdies must be carried out which rely on the value of
$\tau_{\rm remant}$ extracted from NR simulations. The remnant survival time
may sensitively depend on the EOS Hotokezaka  _et al._ (2013); Lucca and
Sagunski (2020), neutrino treatment (relevant for cooling and emergent viscous
effects) Duez  _et al._ (2004); Shibata and Kiuchi (2017); Palenzuela  _et
al._ (2022), magnetic field treatment (relevant for angular momentum
transport) Shapiro (2000); Radice (2017); Shibata  _et al._ (2021); Palenzuela
_et al._ (2022), and (as suggested by the tests presented here) differences in
the implementation of algorithmically equivalent numerical methods.
Additional, in-depth studies remain to be carried to systematically understand
all of the aforementioned effects on $\tau_{\rm remnant}$.

The survival time of the remnant can also play a role in the development of
different fluid instabilities. Fluid instabilities that uniquely develop in
the post-merger environment - such as the one-arm spiral instability and
magneto-rotational instability (MRI) - would be arrested upon gravitational
collapse to a black hole. For instance, a potential observable signature which
arises from the one-arm spiral instability Centrella  _et al._ (2001); Saijo
_et al._ (2003) is the long-lived, sustained powering of the $l=2$, $m=1$ GW
mode Paschalidis  _et al._ (2015), which typically has half the characteristic
frequency of the initially dominant but decaying $l=2$, $m=2$ GW mode
Paschalidis  _et al._ (2015); Radice  _et al._ (2016). Over relatively long
timescales, and as long as the post-merger massive neutron star remnant has
not collapsed to a black hole, the one-arm spiral mode can continuously power
the emission of the $l=2$, $m=1$ GW mode Paschalidis  _et al._ (2015);
simulations which produce shorter remant survival times may not allow the one-
arm spiral instability to develop. In Fig. 10 we show the $l=2$, $m=2$ (using
solid lines) and $l=2$, $m=1$ (using dashed lines) GW modes for test LSBNS. In
all simulations we see increased emission of the $m=1$ mode at a time close to
merger, suggesting the potential early development of the one-arm spiral
instability. In the WhiskyTHC and GRHydro evolutions, the amplitude of the
$m=1$ mode remains roughly constant until collapse, at which point the
emission in this mode ceases. However, in the case of the MIL evolution, which
has not collapsed by the end of simulation, GW emission in the $m=1$ mode
continues and is further amplified, while emission in the $m=2$ mode decays.
Due to the different remnant collapse times, the $m=1$ mode is allowed to
continue growing in the MIL evolution, but not in the WhiskyTHC and GRHydro
evolutions.

Another set of instabilities that is relelvant in the post-merger environment,
and are thereby affected by the remnant collapse time, are magnetic
instabilities. In particular, the MRI can develop in the strongly
differentially rotating environments following a BNS merger Siegel  _et al._
(2013); Kiuchi  _et al._ (2014); Kaplan  _et al._ (2014); Guilet  _et al._
(2016); Ruiz  _et al._ (2016a); Paschalidis (2017). As long as differential
rotation with an appropriate profile is sustained in a massive neutron star
post-merger remnant, the MRI may develop and produce exponential magnetic
field amplification (see Paschalidis and Stergioulas (2017); Ciolfi (2020b)
and references therein). To consider the amount of differential rotation in
the post-merger remnants depicted in Fig. 9 we calculate the approximate
angular velocity in the orbital plane $\Omega=v^{\phi}$, using Eq. (24). In
Fig. 11 we show the radial profile of the angular velocity $\Omega(\varpi)$,
averaged along the azimuthal direction and in time, for test LSBNS in the case
of evolution with the MIL code (depicted using purple lines and contours) and
the GRHydro code (depicted using green lines and contours). We time-average
over a time window of approximately $10\text{\,}\mathrm{m}\mathrm{s}$ ranging
from $t\approx$20\text{\,}\mathrm{m}\mathrm{s}$$ to
$t\approx$30\text{\,}\mathrm{m}\mathrm{s}$$ for the MIL case and from
$t\approx$14\text{\,}\mathrm{m}\mathrm{s}$$ to
$t\approx$24\text{\,}\mathrm{m}\mathrm{s}$$ in the GRHydro case. We find that
by the end of each simulation the remnants still exhibits significant amounts
of differential rotation. The supramassive limit for EOS LS220 falls below the
total system mass. Therefore, we may reasonably expect a delayed collapse once
the configuration is no longer supported by differential rotation. The
timescale on which differential rotation is removed in BNS merger remnants is
typically $\sim\mathcal{O}($100\text{\,}\mathrm{m}\mathrm{s}$)$ Ciolfi  _et
al._ (2019); Ciolfi and Kalinani (2020), which is significantly longer than
the timescales probed by our dynamical simulations. However, such arguments
assume that typically a large fraction of the angular momentum must be lost or
transport for collapse to proceed, which is not generally the case. Therefore,
although we may expect the same final remnants – a black hole – in all of our
simulations for test LSBNS, the timescales over which this happens varies
significantly across codes. We note that similar amounts of differential
rotation are observed in the GRHydro case as well, suggesting that it was not
the lack of differential rotation support which led to the collapse of the
post-merger remnants in that case. As the massive NS remnants for both the MIL
and GRHydro cases exhibit a significant amount of differential rotation by the
end of simulation (or until collapse in the GRHydro case), we may expect the
MRI to develop, if sufficiently resolved, in analogous simulations which
include magnetic fields. However, our results suggest that the MRI could have
a significantly longer time to develop in evolutions with the MIL code, due to
a potentially significantly longer remnant lifetime than we can probe. On the
other hand, the magnetic field can also reduce or prolong the lifetime of the
remnant. The comparison cases presented in this section suggest that
understanding the error introduced by using different numerical methods in NR
simulations of BNS mergers is crucially important to the questions of magnetic
field amplification mechanisms, fluid instabilities, and observables
associated with such systems.

##  V Conclusion and Outlook

In this work we have demonstrated the capabilities of the MIL code, a code
which builds upon the current open source version of IllinoisGRMHD by
including necessary changes for microphysical EOS compatibility. We outlined
the algorithmic structure of IllinoisGRMHD and expanded on the additions that
MIL provides. Specifically, we discussed the implementation of the electron
fraction advection equation and the implementation of new state-of-the-art
conservatives-to-primitives solvers within the code. The main products of
these aforementioned additions are the MIL code itself and a new standalone
thorn from primitives recovery, ConservativeToPrimitive, which is in principle
able to interface with any GRMHD code that uses the Cactus infrastructure.
Over a barrage of tests we established that the MIL code works as well as
other publicly available GRMHD codes with similar capabilities. Crucially, for
the systems considered, the MIL code agrees to a large extent with the Spritz
code GRHydro, and WhiskyTHC codes.

Throughout the work we considered the numerical solutions provided by MIL for
a wide range of astrophysical relevance. In particular, we consider several
TOV and BNS configurations, all of which provide stringent tests of the MIL
code. Where relevant, we compare with other open source GRMHD codes. We find
general agreement and similar convergence properties to the GRHydro, Spritz,
and WhiskyTHC codes, as well as to the current open source version of
IllinoisGRMHD (referred to as OIL throughout the work). In some cases we find
improved performance for MIL when compared to GRHydro (a property which is
carried over from the OIL code as explored in Etienne  _et al._ (2015)). We
find general quantitative agreement between all codes considered in scenarios
where the matter remains at low temperatures throughout the evolution. The
largest differences observed between the results of each code appear during
highly dynamical scenarios such as during BNS mergers and in the immediate
post-merger evolution. Although we find agreement in the post-merger peak
frequency of the characteristic strain, other quantitative differences arise
in these highly dynamical scenarios, likely due to the differences in the
implementation of conservative-to-primitive algorithms. For example, some
differences we highlight include: (1) the GRHydro and WhiskyTHC codes leading
to higher temperatures than the MIL code in a BNS merger scenario, typically
by up to $\Delta T\approx$10\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$; (2)
merger times which differ by approximately $10\%$ between codes; (3)
significantly different survival times of the post-merger remnant; (4) large
oscillations in the rest mass density for GRHydro when comparing to other
codes, and general lack of quantitative agreement in the post-merger
evolution; (5) better preservation of the initial magnetic field structure
during the inspiral in the OIL and MIL codes when compared to the Spritz code.

Given that all codes considered allow some freedom in the flux treatment and
reconstruction scheme, we ensure in all our tests that these two aspects of
the codes have maximal algorithmic overlap. However, even if two codes employ
the same algorithm, the implementation methods will be different in each case.
These differences in algorithmic implementation may lead to the quantitative
differences noted throughout the work. Despite the quantitative differences
highlighted, we emphasize that all codes nevertheless produce qualitatively
consistent (magneto)hydrodynamics evolutions and result in quantitatively
similar gravitational waves. We note that the MIL code provides additional
freedom over GRHydro, Spritz, and WhiskyTHC in the primitives inversion
algorithm used during any given simulation through the thorn
ConservativeToPrimitive, allowing for suitable comparisons with different
codes and the ability to use the optimal conservative-to-primitive solver for
the problem at hand. Additionally, the use of several solvers within
ConservativeToPrimitive for a given simulation allows for a robust and
efficient primitives recovery solution in the case of finite temperature EOSs,
regardless of the variable space probed.

The updates included in MIL make possible the investigation of important, open
questions surrounding BNS mergers. In particular, MIL will allow us to
investigate the interplay between the effects of magnetic fields and
microphysical and finite temperature EOSs in BNS mergers, the calculation of
nucleosynthesis yields associated with BNS mergers, and provides the
groundwork for the inclusion of neutrino effects within IllinoisGRMHD. Future
work utilizing the full capabilities of MIL will investigate all of the
aforementioned phenomena. For instance, the tests presented in Sec. IV.6 will
be expanded upon for a full investigation of the interplay between magnetic
field and realistic EOS effects. Additionally, because the MIL code provides a
dynamical description of several fluid and thermodynamic quantities –
including the rest mass density, temperature, specific entropy, and electron
fraction – its output may be used with nuclear reaction network codes such as
SkyNet Lippuner and Roberts (2017) which only require such quantities as
input. In turn, nuclear reaction networks provide insight into the
nucleosynthesis rates associated with different stages of BNS mergers. Full
GRMHD simulations of BNS mergers with magnetic field constrained transport
have not received much attention in the context of nucleosynthesis, and this
is now possible with MIL. Finally, because the electron fraction is known at
all times during evolutions with MIL, the code may be expanded further to
include neutrino effects. For instance, one may include neutrino leakage by
interfacing MIL with the ZelmaniLeak thorn Ott  _et al._ (2009), or include
more accurate descriptions of neutrino transport that move beyond leakage
Foucart  _et al._ (2020); Nedora  _et al._ (2021, 2020); Weih  _et al._
(2020); Radice  _et al._ (2022). Such an extension to MIL will help elucidate
the combined effect of magnetic and neutrino driven outflows after a BNS
merger. Specifically, neutrino effects are expected to have a significant role
on the electron fraction of the outflow Radice  _et al._ (2022). The
constrained transport magnetic field evolution and finite temperature
description of MIL coupled with accurate neutrino transport could provide some
of the most detailed insight into the physics following a BNS merger. We note
that another recent extension to IllinoisGRMHD similar to ours (albeit using
alternative methods and implementations) was considered in Werneck  _et al._
(2022). A detailed comparison of the codes considered in this work along with
the results presented in Werneck  _et al._ (2022) is a crucial step toward
understanding systematic errors in NR simulations of BNS mergers. We hope to
consider these and other investigations in future work with the MIL code. We
conclude the present work with a note about the availability of the MIL code.
The source code for MIL, along with the supporting thorns, is publicly
available at https://github.com/pilambdaepsilon/UAThorns.We also provide
examples of parameter files, initial data, and equation of state files for all
of the examples considered in this work.

##  VI Acknowledgments

This work was supported by NSF Grants PHY-1912619 and PHY-2145421 to the
University of Arizona. PE also acknowledges support from NSF Grant PHY-2020275
(Network for Neutrinos, Nuclear Astrophysics, and Symmetries (N3AS)). During
completion of this work, PE was in part supported by the Marshall Foundation
Dissertation Scholarship. GB is supported by NASA Grant 80NSSC20K1542 to the
University of Arizona. We wish to thank Bruno Giacomazzo, Federico Cipolletta,
and David Radice for discussions. Simulations were in part performed at the
Comet and Expanse clusters at SDSC, and the Stampede2 cluster at TACC through
XSEDE grant TG-PHY190020. Simulations were also performed on the ElGato,
Ocelote, and Puma clusters at the University of Arizona. kuibit Bozzola (2021)
is based on NumPy Harris  _et al._ (2020), SciPy Virtanen  _et al._ (2020),
and h5py Collette (2013).

## References

  * Abbott  _et al._ (2017) B. P. Abbott _et al._ (Virgo, LIGO Scientific), Phys. Rev. Lett. 119, 161101 (2017), arXiv:1710.05832 [gr-qc] .
  * Raaijmakers  _et al._ (2021) G. Raaijmakers _et al._ , (2021), arXiv:2102.11569 [astro-ph.HE] .
  * Easter  _et al._ (2020) P. J. Easter, S. Ghonge, P. D. Lasky, A. R. Casey, J. A. Clark, F. H. Vivanco, and K. Chatziioannou, Phys. Rev. D 102, 043011 (2020), arXiv:2006.04396 [astro-ph.HE] .
  * Barbieri  _et al._ (2020) C. Barbieri, O. S. Salafia, M. Colpi, G. Ghirlanda, and A. Perego, (2020), arXiv:2002.09395 [astro-ph.HE] .
  * Breschi  _et al._ (2021) M. Breschi, R. Gamba, and S. Bernuzzi, (2021), arXiv:2102.00017 [gr-qc] .
  * Nicholl  _et al._ (2021) M. Nicholl, B. Margalit, P. Schmidt, G. P. Smith, E. J. Ridley, and J. Nuttall, (2021), arXiv:2102.02229 [astro-ph.HE] .
  * Paschalidis  _et al._ (2015) V. Paschalidis, M. Ruiz, and S. L. Shapiro, Astrophys. J. Lett. 806, L14 (2015).
  * Ruiz  _et al._ (2018) M. Ruiz, S. L. Shapiro, and A. Tsokaros, Phys. Rev. D97, 021501 (2018), arXiv:1711.00473 [astro-ph.HE] .
  * Ruiz  _et al._ (2016a) M. Ruiz, R. N. Lang, V. Paschalidis, and S. L. Shapiro, The Astrophysical Journal 824, L6 (2016a).
  * Kiuchi  _et al._ (2018) K. Kiuchi, K. Kyutoku, Y. Sekiguchi, and M. Shibata, Phys. Rev. D 97, 124039 (2018), arXiv:1710.01311 [astro-ph.HE] .
  * Shibata  _et al._ (2021) M. Shibata, S. Fujibayashi, and Y. Sekiguchi, Phys. Rev. D 103, 043022 (2021), arXiv:2102.01346 [astro-ph.HE] .
  * Ciolfi and Kalinani (2020) R. Ciolfi and J. V. Kalinani, Astrophys. J. Lett. 900, L35 (2020), arXiv:2004.11298 [astro-ph.HE] .
  * Mosta  _et al._ (2020) P. Mosta, D. Radice, R. Haas, E. Schnetter, and S. Bernuzzi, Astrophys. J. Lett. 901, L37 (2020), arXiv:2003.06043 [astro-ph.HE] .
  * Ciolfi (2020a) R. Ciolfi, Mon. Not. Roy. Astron. Soc. 495, L66 (2020a), arXiv:2001.10241 [astro-ph.HE] .
  * Abdikamalov  _et al._ (2009) E. B. Abdikamalov, H. Dimmelmeier, L. Rezzolla, and J. C. Miller, Mon. Not. Roy. Astron. Soc. 394, 52 (2009), arXiv:0806.1700 [astro-ph] .
  * Bauswein  _et al._ (2013) A. Bauswein, T. W. Baumgarte, and H. T. Janka, Phys. Rev. Lett. 111, 131101 (2013), arXiv:1307.5191 [astro-ph.SR] .
  * Palenzuela  _et al._ (2015) C. Palenzuela, S. L. Liebling, D. Neilsen, L. Lehner, O. Caballero, E. O?Connor, and M. Anderson, Physical Review D 92 (2015), 10.1103/physrevd.92.044045.
  * Sekiguchi  _et al._ (2015) Y. Sekiguchi, K. Kiuchi, K. Kyutoku, and M. Shibata, Phys. Rev. D 91, 064059 (2015), arXiv:1502.06660 [astro-ph.HE] .
  * Lehner  _et al._ (2016) L. Lehner, S. L. Liebling, C. Palenzuela, O. L. Caballero, E. O’Connor, M. Anderson, and D. Neilsen, Class. Quant. Grav. 33, 184002 (2016), arXiv:1603.00501 [gr-qc] .
  * Radice  _et al._ (2018a) D. Radice, A. Perego, K. Hotokezaka, S. A. Fromm, S. Bernuzzi, and L. F. Roberts, Astrophys. J. 869, 130 (2018a), arXiv:1809.11161 [astro-ph.HE] .
  * Most  _et al._ (2019a) E. R. Most, L. J. Papenfort, V. Dexheimer, M. Hanauske, S. Schramm, H. Stoecker, and L. Rezzolla, Phys. Rev. Lett. 122, 061101 (2019a), arXiv:1807.03684 [astro-ph.HE] .
  * Dietrich  _et al._ (2018) T. Dietrich, D. Radice, S. Bernuzzi, F. Zappa, A. Perego, B. Brügmann, S. V. Chaurasia, R. Dudi, W. Tichy, and M. Ujevic, Class. Quant. Grav. 35, 24LT01 (2018), arXiv:1806.01625 [gr-qc] .
  * Vincent  _et al._ (2020) T. Vincent, F. Foucart, M. D. Duez, R. Haas, L. E. Kidder, H. P. Pfeiffer, and M. A. Scheel, Phys. Rev. D 101, 044053 (2020), arXiv:1908.00655 [gr-qc] .
  * Raithel  _et al._ (2021a) C. Raithel, V. Paschalidis, and F. Özel, Phys. Rev. D 104, 063016 (2021a), arXiv:2104.07226 [astro-ph.HE] .
  * Nedora  _et al._ (2019) V. Nedora, S. Bernuzzi, D. Radice, A. Perego, A. Endrizzi, and N. Ortiz, Astrophys. J. Lett. 886, L30 (2019), arXiv:1907.04872 [astro-ph.HE] .
  * Foucart  _et al._ (2020) F. Foucart, M. D. Duez, F. Hebert, L. E. Kidder, H. P. Pfeiffer, and M. A. Scheel, The Astrophysical Journal 902, L27 (2020).
  * Paschalidis (2017) V. Paschalidis, Class. Quant. Grav. 34, 084002 (2017), arXiv:1611.01519 [astro-ph.HE] .
  * Baiotti and Rezzolla (2017) L. Baiotti and L. Rezzolla, Reports on Progress in Physics 80, 096901 (2017), arXiv:1607.03540 [gr-qc] .
  * Duez and Zlochower (2019) M. D. Duez and Y. Zlochower, Reports on Progress in Physics 82, 016902 (2019), arXiv:1808.06011 [gr-qc] .
  * Radice  _et al._ (2020) D. Radice, S. Bernuzzi, and A. Perego, Annual Review of Nuclear and Particle Science 70, 95 (2020).
  * Bernuzzi (2020) S. Bernuzzi, Gen. Rel. Grav. 52, 108 (2020), arXiv:2004.06419 [astro-ph.HE] .
  * Ciolfi (2020b) R. Ciolfi, Gen. Rel. Grav. 52, 59 (2020b), arXiv:2003.07572 [astro-ph.HE] .
  * Kasen  _et al._ (2013) D. Kasen, N. R. Badnell, and J. Barnes, The Astrophysical Journal 774, 25 (2013).
  * Barnes and Kasen (2013) J. Barnes and D. Kasen, Astrophys. J. 775, 18 (2013), arXiv:1303.5787 [astro-ph.HE] .
  * Rezzolla  _et al._ (2018) L. Rezzolla, E. R. Most, and L. R. Weih, Astrophys. J. 852, L25 (2018), arXiv:1711.00314 [astro-ph.HE] .
  * Raithel  _et al._ (2021b) C. Raithel, V. Paschalidis, and F. Ozel, “Realistic finite-temperature effects in neutron star merger simulations,” (2021b), arXiv:2104.07226 [astro-ph.HE] .
  * Blandford and Znajek (1977) R. D. Blandford and R. L. Znajek, Mon. Not. Roy. Astron. Soc. 179, 433 (1977).
  * Van Putten (1999) M. H. P. M. Van Putten, Science 284, 115 (1999).
  * Bucciantini  _et al._ (2012) N. Bucciantini, B. D. Metzger, T. A. Thompson, and E. Quataert, Mon. Not. Roy. Astron. Soc. 419, 1537 (2012), arXiv:1106.4668 [astro-ph.HE] .
  * Ruiz  _et al._ (2016b) M. Ruiz, R. N. Lang, V. Paschalidis, and S. L. Shapiro, Astrophys. J. 824, L6 (2016b).
  * Siegel  _et al._ (2014) D. M. Siegel, R. Ciolfi, and L. Rezzolla, Astrophys. J. Lett. 785, L6 (2014), arXiv:1401.4544 [astro-ph.HE] .
  * Fernandez  _et al._ (2019) R. Fernandez, A. Tchekhovskoy, E. Quataert, F. Foucart, and D. Kasen, Mon. Not. Roy. Astron. Soc. 482, 3373 (2019), arXiv:1808.00461 [astro-ph.HE] .
  * Metzger  _et al._ (2018) B. D. Metzger, T. A. Thompson, and E. Quataert, Astrophys. J. 856, 101 (2018), arXiv:1801.04286 [astro-ph.HE] .
  * Wanajo  _et al._ (2014) S. Wanajo, Y. Sekiguchi, N. Nishimura, K. Kiuchi, K. Kyutoku, and M. Shibata, Astrophys. J. Lett. 789, L39 (2014), arXiv:1402.7317 [astro-ph.SR] .
  * Perego  _et al._ (2017) A. Perego, D. Radice, and S. Bernuzzi, Astrophys. J. Lett. 850, L37 (2017), arXiv:1711.03982 [astro-ph.HE] .
  * Radice  _et al._ (2018b) D. Radice, A. Perego, K. Hotokezaka, S. Bernuzzi, S. A. Fromm, and L. F. Roberts, Astrophys. J. Lett. 869, L35 (2018b), arXiv:1809.11163 [astro-ph.HE] .
  * Sun  _et al._ (2022) L. Sun, M. Ruiz, S. L. Shapiro, and A. Tsokaros, Phys. Rev. D 105, 104028 (2022), arXiv:2202.12901 [astro-ph.HE] .
  * Shibata  _et al._ (2017) M. Shibata, S. Fujibayashi, K. Hotokezaka, K. Kiuchi, K. Kyutoku, Y. Sekiguchi, and M. Tanaka, Phys. Rev. D96, 123012 (2017), arXiv:1710.07579 [astro-ph.HE] .
  * Bovard  _et al._ (2017) L. Bovard, D. Martin, F. Guercilena, A. Arcones, L. Rezzolla, and O. Korobkin, Phys. Rev. D 96, 124005 (2017), arXiv:1709.09630 [gr-qc] .
  * Radice  _et al._ (2018) D. Radice, A. Perego, S. Bernuzzi, and B. Zhang, ArXiv e-prints (2018), arXiv:1803.10865 [astro-ph.HE] .
  * Foucart  _et al._ (2016) F. Foucart, E. O’Connor, L. Roberts, L. E. Kidder, H. P. Pfeiffer, and M. A. Scheel, Phys. Rev. D94, 123016 (2016), arXiv:1607.07450 [astro-ph.HE] .
  * Radice  _et al._ (2018) D. Radice, A. Perego, K. Hotokezaka, S. Bernuzzi, S. A. Fromm, and L. F. Roberts, The Astrophysical Journal 869, L35 (2018).
  * Bernuzzi  _et al._ (2020) S. Bernuzzi _et al._ , Mon. Not. Roy. Astron. Soc. 497, 1488 (2020), arXiv:2003.06015 [astro-ph.HE] .
  * Pürrer and Haster (2020) M. Pürrer and C.-J. Haster, Phys. Rev. Res. 2, 023151 (2020), arXiv:1912.10055 [gr-qc] .
  * Werneck  _et al._ (2022) L. R. Werneck _et al._ , (2022), arXiv:2208.14487 [gr-qc] .
  * Etienne  _et al._ (2015) Z. B. Etienne, V. Paschalidis, R. Haas, P. Moesta, and S. L. Shapiro, Class. Quant. Grav. 32, 175009 (2015), arXiv:1501.07276 [astro-ph.HE] .
  * Mosta  _et al._ (2014) P. Mosta, B. C. Mundim, J. A. Faber, R. Haas, S. C. Noble, T. Bode, F. Loffler, C. D. Ott, C. Reisswig, and E. Schnetter, Class. Quant. Grav. 31, 015005 (2014), arXiv:1304.5544 [gr-qc] .
  * Giacomazzo  _et al._ (2020) B. Giacomazzo, F. Cipolletta, J. Kalinani, R. Ciolfi, L. Sala, B. Giudici, and E. Giangrandi, “The Spritz Code,”  (2020).
  * Radice and Rezzolla (2012) D. Radice and L. Rezzolla, Astronomy and Astrophysics 547, A26 (2012).
  * Radice  _et al._ (2013) D. Radice, L. Rezzolla, and F. Galeazzi, Monthly Notices of the Royal Astronomical Society: Letters 437, L46 (2013).
  * Radice  _et al._ (2014) D. Radice, L. Rezzolla, and F. Galeazzi, Classical and Quantum Gravity 31, 075012 (2014).
  * Loffler  _et al._ (2012) F. Loffler, J. Faber, E. Bentivegna, T. Bode, P. Diener, R. Haas, I. Hinder, B. C. Mundim, C. D. Ott, E. Schnetter, G. Allen, M. Campanelli, and P. Laguna, Classical and Quantum Gravity 29, 115001 (2012).
  * Etienne  _et al._ (2021) Z. Etienne, S. R. Brandt, P. Diener, W. E. Gabella, M. Gracia-Linares, R. Haas, A. Kedia, M. Alcubierre, D. Alic, G. Allen, M. Ansorg, M. Babiuc-Hamilton, L. Baiotti, W. Benger, E. Bentivegna, S. Bernuzzi, T. Bode, G. Bozzola, B. Brendal, B. Bruegmann, M. Campanelli, F. Cipolletta, G. Corvino, S. Cupp, R. D. Pietri, H. Dimmelmeier, R. Dooley, N. Dorband, M. Elley, Y. E. Khamra, J. Faber, T. Font, J. Frieben, B. Giacomazzo, T. Goodale, C. Gundlach, I. Hawke, S. Hawley, I. Hinder, E. A. Huerta, S. Husa, S. Iyer, D. Johnson, T. Kellermann, A. Knapp, M. Koppitz, P. Laguna, G. Lanferman, F. Loffler, J. Masso, L. Menger, A. Merzky, J. M. Miller, M. Miller, P. Moesta, P. Montero, B. Mundim, A. Nerozzi, S. C. Noble, C. Ott, R. Paruchuri, D. Pollney, D. Radice, T. Radke, C. Reisswig, L. Rezzolla, D. Rideout, M. Ripeanu, L. Sala, J. A. Schewtschenko, E. Schnetter, B. Schutz, E. Seidel, E. Seidel, J. Shalf, K. Sible, U. Sperhake, N. Stergioulas, W.-M. Suen, B. Szilagyi, R. Takahashi, M. Thomas, J. Thornburg, M. Tobias, A. Tonita, P. Walker, M.-B. Wan, B. Wardell, L. Werneck, H. Witek, M. Zilhao, B. Zink, and Y. Zlochower, “The einstein toolkit,”  (2021), to find out more, visit http://einsteintoolkit.org.
  * Bozzola (2021) G. Bozzola, Journal of Open Source Software 6, 3099 (2021).
  * 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 [gr-qc] .
  * Alcubierre (2008) M. Alcubierre,  _Introduction to 3+1 Numerical Relativity_ (Oxford University Press, UK, 2008).
  * Bona and Palenzuela-Luque (2005) C. Bona and C. Palenzuela-Luque,  _Elements of Numerical Relativity_, Vol. 673 (2005).
  * Baumgarte and Shapiro (2010) T. W. Baumgarte and S. L. Shapiro,  _Numerical Relativity: Solving Einstein’s Equations on the Computer_ (Cambridge University Press, 2010).
  * Shibata (2016) M. Shibata,  _Numerical Relativity_ (World Scientific Publishing Co, 2016).
  * Font (2008) J. A. Font, Living Rev. Rel. 11, 7 (2008).
  * Rezzolla and Zanotti (2013) L. Rezzolla and O. Zanotti,  _Relativistic Hydrodynamics_ (2013).
  * 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).
  * Anton  _et al._ (2006) L. Anton, O. Zanotti, J. A. Miralles, J. M. Marti, J. M. Ibanez, J. A. Font, and J. A. Pons, Astrophys. J. 637, 296 (2006), arXiv:astro-ph/0506063 .
  * Farris  _et al._ (2012) B. D. Farris, R. Gold, V. Paschalidis, Z. B. Etienne, and S. L. Shapiro, Physical Review Letters 109, 221102 (2012), arXiv:1207.3354 [astro-ph.HE] .
  * Most and Raithel (2021) E. R. Most and C. A. Raithel, Phys. Rev. D 104, 124012 (2021), arXiv:2107.06804 [astro-ph.HE] .
  * Raithel and Paschalidis (2022) C. A. Raithel and V. Paschalidis, Phys. Rev. D 106, 023015 (2022), arXiv:2204.00698 [gr-qc] .
  * Raithel  _et al._ (2022) C. A. Raithel, P. Espino, and V. Paschalidis,  (2022), 10.1093/mnras/stac2450, arXiv:2206.14838 [astro-ph.HE] .
  * Colella and Woodward (1984) P. Colella and P. R. Woodward, J. Comput. Phys. 54, 174 (1984).
  * Harten  _et al._ (1983) A. Harten, P. D. Lax, and B. V. Leer, SIAM Review 25, 35 (1983).
  * Etienne  _et al._ (2010) Z. B. Etienne, Y. T. Liu, and S. L. Shapiro, Phys. Rev. D 82, 084031 (2010).
  * Etienne  _et al._ (2012) Z. B. Etienne, Y. T. Liu, V. Paschalidis, and S. L. Shapiro, Phys. Rev. D 85, 064029 (2012), arXiv:1112.0568 [astro-ph.HE] .
  * Siegel  _et al._ (2018) D. M. Siegel, P. Mosta, D. Desai, and S. Wu, Astrophys. J. 859, 71 (2018), arXiv:1712.07538 [astro-ph.HE] .
  * O’Connor and Ott (2010) E. O’Connor and C. D. Ott, Class. Quant. Grav. 27, 114103 (2010), arXiv:0912.2393 [astro-ph.HE] .
  * Lattimer and Swesty (1991) J. M. Lattimer and D. F. Swesty, Nucl. Phys. A 535, 331 (1991).
  * Gammie  _et al._ (2003) C. F. Gammie, J. C. McKinney, and G. Toth, Astrophys. J. 589, 444 (2003), arXiv:astro-ph/0301509 [astro-ph] .
  * 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 [astro-ph] .
  * Newman and Hamlin (2014) W. I. Newman and N. D. Hamlin, SIAM Journal on Scientific Computing 36, B661 (2014), https://doi.org/10.1137/140956749 .
  * Kastaun  _et al._ (2021) W. Kastaun, J. V. Kalinani, and R. Ciolfi, Phys. Rev. D 103, 023018 (2021), arXiv:2005.01821 [gr-qc] .
  * Most  _et al._ (2019b) E. R. Most, L. J. Papenfort, and L. Rezzolla, Mon. Not. Roy. Astron. Soc. 490, 3588 (2019b), arXiv:1907.10328 [astro-ph.HE] .
  * Kalinani  _et al._ (2022) J. V. Kalinani, R. Ciolfi, W. Kastaun, B. Giacomazzo, F. Cipolletta, and L. Ennoggi, Phys. Rev. D 105, 103031 (2022), arXiv:2107.10620 [astro-ph.HE] .
  * Einfeldt (1988) B. Einfeldt, SIAM Journal on Numerical Analysis 25, 294 (1988), https://doi.org/10.1137/0725021 .
  * Brown  _et al._ (2009) D. Brown, P. Diener, O. Sarbach, E. Schnetter, and M. Tiglio, Phys. Rev. D 79, 044023 (2009), arXiv:0809.3533 [gr-qc] .
  * Reisswig  _et al._ (2011) C. Reisswig, C. D. Ott, U. Sperhake, and E. Schnetter, Physical Review D 83 (2011), 10.1103/physrevd.83.064008.
  * Bona  _et al._ (1995) C. Bona, J. Masso, E. Seidel, and J. Stela, Phys. Rev. Lett. 75, 600 (1995), arXiv:gr-qc/9412071 [gr-qc] .
  * Alcubierre  _et al._ (2003) M. Alcubierre, B. Bruegmann, P. Diener, M. Koppitz, D. Pollney, E. Seidel, and R. Takahashi, Phys. Rev. D67, 084023 (2003), arXiv:gr-qc/0206072 [gr-qc] .
  * Newman and Penrose (1962) E. Newman and R. Penrose, J. Math. Phys. 3, 566 (1962).
  * Penrose (1963) R. Penrose, Phys. Rev. Lett. 10, 66 (1963).