# Implementation of advanced Riemann solvers in a neutrino-radiation
magnetohydrodynamics code in numerical relativity and its application to a
binary neutron star merger

Kenta Kiuchi  Max Planck Institute for Gravitational Physics (Albert Einstein
Institute), Am Mühlenberg, Potsdam-Golm, 14476, Germany  Center for
Gravitational Physics and Quantum Information, Yukawa Institute for
Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan  Loren E. Held
Max Planck Institute for Gravitational Physics (Albert Einstein Institute), Am
Mühlenberg, Potsdam-Golm, 14476, Germany  Yuichiro Sekiguchi  Department of
Physics, Toho University, Funabashi, Chiba 274-8510, Japan  Center for
Gravitational Physics and Quantum Information, Yukawa Institute for
Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan  Masaru Shibata
Max Planck Institute for Gravitational Physics (Albert Einstein Institute), Am
Mühlenberg, Potsdam-Golm, 14476, Germany  Center for Gravitational Physics and
Quantum Information, Yukawa Institute for Theoretical Physics, Kyoto
University, Kyoto 606-8502, Japan

###### Abstract

We implement advanced Riemann solvers HLLC and HLLD Mignone and Bodo (2005);
Mignone  _et al._ (2009) together with an advanced constrained transport
scheme Gardiner and Stone (2008) in a numerical-relativity neutrino-radiation
magnetohydrodynamics code. We validate our implementation by performing a
series of one- and multi-dimensional test problems for relativistic
hydrodynamics and magnetohydrodynamics in both Minkowski spacetime and a
static black hole spacetime. We find that the numerical solutions with the
advanced Riemann solvers are more accurate than those with the HLLE solver Del
Zanna  _et al._ (2003), which was originally implemented in our code. As an
application to numerical relativity, we simulate an asymmetric binary neutron
star merger leading to a short-lived massive neutron star both with and
without magnetic fields. We find that the lifetime of the rotating massive
neutron star formed after the merger and also the amount of the tidally-driven
dynamical ejecta are overestimated when we employ the diffusive HLLE solver.
We also find that the magnetorotational instability is less resolved when we
employ the HLLE solver because of the solver’s large numerical diffusivity.
This causes a spurious enhancement both of magnetic winding resulting from
large scale poloidal magnetic fields, and also of the energy of the outflow
induced by magnetic pressure.

##  I Introduction

The first direct detection of gravitational waves from a binary neutron star
merger (GW170817) and its electromagnetic counterparts (AT 2017gfo/SGRB
170817A) heralded the beginning of multimessenger astronomy including
gravitational waves Abbott  _et al._ (2017a, b). In this event, the tidal
deformability of the neutron star binary was measured for the first time and
found to be in the interval $100\lesssim\tilde{\Lambda}\lesssim 800$, with an
accurate measurement of the total mass of the binary yielding
$2.73^{+0.04}_{-0.01}M_{\odot}$ Abbott  _et al._ (2017a, 2019, 2018); De  _et
al._ (2018) 111The precise value of $\tilde{\Lambda}$ depends on the details
of the analysis.. Any viable neutron star matter equations of state must
satisfy this observational constraint on tidal deformability. In this event it
was also shown that the binary neutron star merger drives a short gamma-ray
burst Goldstein  _et al._ (2017); Abbott  _et al._ (2017b); Savchenko  _et
al._ (2017); Mooley  _et al._ (2018), thus providing the first ‘smoking gun’
for supporting the hypothesis that binary mergers can be the central engine of
short gamma-ray bursts Paczynski (1986); Goodman (1986); Eichler  _et al._
(1989); Narayan  _et al._ (1992). Finally, this event indicated that neutron-
rich matter is likely to be ejected during the merger and heavy elements are
synthesized within these ejecta by means of the rapid neutron capture process
on nuclei (the $r$-process) Metzger  _et al._ (2010); Lattimer and Schramm
(1974); Eichler  _et al._ (1989); Wanajo  _et al._ (2014). It had been
predicted that the $r$-process nucleosynthesis subsequently causes so-called
kilonova emission via the radioactive decay of unstable $r$-process elements
Tanaka and Hotokezaka (2013); Barnes and Kasen (2013), and a kilonova was
indeed observed after the merger in the near infrared, optical, and
ultraviolet bands Arcavi  _et al._ (2017); Chornock  _et al._ (2017); Coulter
_et al._ (2017); Cowperthwaite  _et al._ (2017); Drout  _et al._ (2017); Kasen
_et al._ (2017); Kasliwal  _et al._ (2017); Kilpatrick  _et al._ (2017);
McCully  _et al._ (2017); Nicholl  _et al._ (2017); Shappee  _et al._ (2017);
Smartt  _et al._ (2017); Soares-Santos  _et al._ (2017); Tanaka  _et al._
(2017); Tanvir  _et al._ (2017).

References Abbott  _et al._ (2020, 2021) reported the detection of a second
binary neutron star merger event (GW190425) and measured a total binary mass
of $3.4^{+0.1}_{-0.1}M_{\odot}$, which is much larger than the total mass
measured in binary pulsars observed in our Galaxy Farrow  _et al._ (2019). The
merger dynamics, mass ejection process, and resultant electromagnetic emission
due to $r$-process nucleosynthesis could be different from those in GW170817
and AT 2017gfo Kyutoku  _et al._ (2020); Dudi  _et al._ (2021). Although an
electromagnetic counterpart was not observed in GW190425, either due to poor
sky localization or due to intrinsically dimmer emission Hosseinzadeh  _et
al._ (2019); Coughlin  _et al._ (2019, 2020), the existence of a massive
binary neutron star suggests that the binary neutron star merger and
associated mass ejection process could have a diversity of mechanisms. The new
observation run O4 is planned to commence at the end of 2022 O4- , and could
lead to the observation of binary neutron star mergers and associated
electromagnetic counterparts that are qualitatively different from those
observed in GW170817. This motivates building binary neutron star merger
models based on reliable numerical relativity simulations for predicting and
interpreting gravitational wave events in preparation for the upcoming
observational run Shibata  _et al._ (2017).

Recent axisymmetric neutrino-radiation viscous-hydrodynamics simulations of
binary neutron star merger remnants in numerical relativity suggest that the
amount of post-merger ejecta launched from the merger remnant due to viscous
effects (which were facilitated in these simulations by an effective ‘alpha’
viscosity parameter) can be larger than the dynamical ejecta launched during
the merger itself Fujibayashi  _et al._ (2018, 2020a, 2020b, 2020c); Shibata
_et al._ (2017). The timescale of the post-merger mass ejection is $O(1)$
second, and depends on the value of the viscosity parameter. Plausible values
of the viscosity parameter are inferred from three-dimensional
magnetohydrodynamics simulations of the binary neutron star merger remnant in
which angular momentum transport is facilitated in a self-consistent manner by
the magnetorotational instability Kiuchi  _et al._ (2018) (see also Refs.
Christie  _et al._ (2019); Fernández  _et al._ (2019) for magnetohydrodynamics
simulations of a massive torus in a stationary black hole spacetime). The
electron fraction of the post-merger ejecta and the resultant $r$-process
nucleosynthesis also depends on this viscosity parameter Fujibayashi  _et al._
(2020b, a, c), although the electron fraction of the post-merger ejecta is
appreciably larger than that of the dynamical ejecta.

Furthermore, very recently we performed neutrino-radiation
magnetohydrodynamics simulations of black hole-neutron star mergers in
numerical relativity Hayashi  _et al._ (2021). We found post-merger mass
ejection due to magnetorotational instability-driven turbulence and the launch
of a Poynting flux-dominated outflow. The post-merger mass ejection and the
Poynting-flux dominated outflow sets in at several $100$ ms after the merger
and lasts for $1$–$2$ seconds after the merger. These timescales are
determined by the strength of the effective viscosity associated with both
magnetorotational-instability turbulence and neutrino cooling Hayashi  _et
al._ (2021).

All these recent studies show that for modeling future gravitational wave
events it is necessary to perform self-consistent (i.e. in which turbulence is
sustained by the magnetorotational instability) three-dimensional neutrino-
radiation magnetohydrodynamics simulations of binary neutron star mergers in
general relativity for the durations of $O(1)$ second. In particular, it is
crucial to reproduce a magneto-turbulent state driven by the magnetorotational
instability inside the merger remnant because the resultant effective
turbulent viscosity transports angular momentum outwards and heats up the
matter via viscous heating Balbus and Hawley (1998).

Finite volume methods are a popular combination of numerical schemes for
simulations of astrophysical fluid dynamics due to their inherent conservation
properties and ability to capture sharp discontinuities in the flow such as
shocks Toro (2013). Central to these schemes is the solution of the so-called
Riemann problem in which one considers two constant states separated by a
discontinuity. The solution consists of three waves in hydrodynamics and seven
waves in magnetohydrodynamics. As exact Riemann solvers are computationally
expensive Giacomazzo and Rezzolla (2006), approximate Riemann solvers are
often used. One such family of approximate Riemann solvers is the HLL-based
Riemann solvers, in which only a subset of the full seven waves in the Riemann
fan are considered. The HLL(E) solver, for example, takes into account only
shocks/rarefactions and omits the contact discontinuity Harten  _et al._
(1983).

At present, the Riemann solver and constrained transport scheme implemented in
existing numerical relativity magnetohydrodynamics codes, e.g. Hayashi  _et
al._ (2021); Radice  _et al._ (2021); Bernuzzi and Dietrich (2016); Most  _et
al._ (2019); Mösta  _et al._ (2014); Etienne  _et al._ (2015); Viganò  _et
al._ (2020); Cipolletta  _et al._ (2021); Foucart  _et al._ (2021), are based
on the HLLE solver Kurganov and Tadmor (2000); Harten  _et al._ (1983); Del
Zanna  _et al._ (2003). (An exception is the SpECTRE Deppe  _et al._ (2021)
code, which is based on the discontinuous Galerkin method.) This Riemann
solver is known to be very diffusive Toro  _et al._ (1994); Mignone  _et al._
(2009); Held and Latter (2018). The numerical diffusion inherent in the
Riemann solver adversely affects the accuracy of the numerical solution, in
particular for long-term simulations of compact object mergers of $O(1)$
second. Although Refs. Bernuzzi and Dietrich (2016); Most  _et al._ (2019);
Radice  _et al._ (2014) reported the implementation of fourth-order accurate
Riemann solvers in their numerical relativity codes, these solvers are based
on the finite difference method. Therefore, it is unclear how accurate these
finite difference-based Riemann solvers are for the problem of astrophysical
turbulence.

This paper reports a new implementation of advanced Riemann solvers in our
neutrino-radiation magnetohydrodynamics numerical relativity code Kiuchi  _et
al._ (2012); Sekiguchi  _et al._ (2012) based on the finite volume method. We
implement the HLLC solver for relativistic hydrodynamics, which restores the
contact discontinuity Mignone and Bodo (2005), and the HLLD solver for
relativistic magnetohydrodynamics, which takes into account five of the seven
waves in the Riemann fan Mignone  _et al._ (2009). Both these Riemann solves
are known to be less diffusive than the HLLE solver Del Zanna  _et al._
(2003). In addition, the constrained transport scheme in Ref. Gardiner and
Stone (2008), which relies on the solution given by a Riemann solver,
significantly suppresses numerical diffusion compared to the HLLE-constrained
transport scheme proposed in Ref. Del Zanna  _et al._ (2003) (see Ref. Mignone
and Del Zanna (2021) for a detailed comparison of different implementations of
the constrained transport scheme). Thus, in addition to implementing the
advanced Riemann solvers, HLLC and HLLD, we also implement the novel
constrained transport scheme of Ref. Gardiner and Stone (2008) in our code.

This paper is organized as follows. Section II summarizes the equations of
motion for general relativistic neutrino-radiation magnetohydrodynamics.
Section III is devoted to the numerical algorithm for general relativistic
magnetohydrodynamics: the finite volume method, the constrained transport
method (for enforcing divergence-free condition of the magnetic field), the
tetrad transformation (which enables us to use Riemann solvers designed for
special relativistic flows in full general relativity), the implementation of
the HLLC solver Mignone and Bodo (2005), that of the HLLD solver Mignone  _et
al._ (2009), and the electric field evaluation (which is used by the
constrained transport algorithm) Gardiner and Stone (2008). In Sec. IV, we
validate our implementation of the new Riemann solvers by performing one- and
multi-dimensional test problems both in Minkowski spacetime and in curved, but
static, spacetime in both relativistic hydrodynamics and magnetohydrodynamics.
Finally, in Sec. V we apply our new solvers in general relativity to a
dynamical spacetime. We first present the results of binary neutron star
merger simulations in the absence of magnetic fields (which are run up to
$\approx 40$–$50$ ms after the formation of the black hole), and subsequently
the evolution of the merger remnant with a magnetic field. Section VI
summarizes our results. Throughout this paper, we use geometrical units in
which $c=G=1$. Greek and Latin indices without hats denote the spacetime and
purely spatial components, respectively. Those with hats indicate tetrad
components.

##  II Governing equations for general relativistic neutrino-radiation
magnetohydrodynamics

In this section, we briefly summarize the set of basic equations of general
relativistic neutrino-radiation magnetohydrodynamics using the 3+1 formalism.
The reader can find a more comprehensive derivation of these equations in,
e.g., Ref. Shibata (2016).

We begin by introducing a unit vector normal to a spatial hypersurface of
constant coordinate time, $t$,

$\displaystyle n^{\mu}$
$\displaystyle=\left(\frac{1}{\alpha},-\frac{\beta^{i}}{\alpha}\right),~{}n_{\mu}=\left(-\alpha,0\right),$
(1)

where $\alpha$ and $\beta^{i}$ are the lapse function and shift vector,
respectively. With this vector, the four dimensional metric can be decomposed
into

$\displaystyle g_{\mu\nu}$
$\displaystyle=\left(\begin{array}[]{cc}-\alpha^{2}+\beta_{i}\beta^{i}&\beta_{i}\\\
\beta_{i}&\gamma_{ij}\end{array}\right),$ (4)

where $\gamma_{ij}$ is the three-dimensional spatial metric.

The stress-energy-momentum tensor for ideal magnetohydrodynamics and for a
free-streaming neutrino-radiation field are, respectively, given by

$\displaystyle T^{\mu\nu}_{\text{(MHD)}}=\rho
hu^{\mu}u^{\nu}+Pg^{\mu\nu}+\left(u^{\mu}u^{\nu}+\frac{1}{2}g^{\mu\nu}\right)b^{2}-b^{\mu}b^{\nu},$
$\displaystyle
T^{\mu\nu}_{\text{(Rad,s,$\nu_{i}$)}}=E_{(\nu_{i})}n^{\mu}n^{\nu}+F^{\mu}_{(\nu_{i})}n^{\nu}+F^{\nu}_{(\nu_{i})}n^{\mu}+P^{\mu\nu}_{(\nu_{i})},$
(5)

where $\rho$, $P$, $u^{\mu}$, $b^{\mu}$, $E_{(\nu_{i})}$,
$F^{\mu}_{(\nu_{i})}$, and $P^{\mu\nu}_{(\nu_{i})}$ are, respectively, the
rest-mass density, pressure, four-velocity, magnetic field (measured in the
fluid rest frame), radiation energy density, radiation momentum, and radiation
stress-energy-momentum tensor of the neutrino species $\nu_{i}$ in the
Eulerian frame. $h=1+\varepsilon+P/\rho$ denotes the relativistic specific
enthalpy with $\varepsilon$ the specific internal energy. We consider the
electron neutrino $\nu_{e}$, electron antineutrino $\bar{\nu}_{e}$, and the
total of $\mu$ and $\tau$ neutrinos and antineutrinos collectively denoted by
$\nu_{x}$ Sekiguchi  _et al._ (2012); Shibata (2016). Note that we assume that
the stress-energy-momentum tensor of the neutrino-radiation field is split
into a trapped component and a free-streaming component. The stress-energy-
momentum tensor of the trapped neutrinos is then absorbed into that for the
ideal magnetohydrodynamics fluid because trapped neutrinos are strongly
coupled to the fluid Sekiguchi  _et al._ (2012); Shibata (2016).

The conserved mass density, total momentum density, and total energy density
of an electrically conducting fluid are defined by

$\displaystyle D$ $\displaystyle\equiv\rho w$ (6) $\displaystyle J_{i}$
$\displaystyle\equiv-{\gamma^{\mu}}_{i}n^{\nu}T_{\mu\nu}^{\text{(MHD)}}$
$\displaystyle=\rho whu_{i}+\frac{B^{2}u_{i}-(B^{j}u_{j})B_{i}}{w},$ (7)
$\displaystyle\rho_{\rm H}$ $\displaystyle\equiv
n^{\mu}n^{\nu}T_{\mu\nu}^{\text{(MHD)}}$ $\displaystyle=\rho
w\left(hw-\frac{P}{\rho
w}\right)+\left(w^{2}-\frac{1}{2}\right)b^{2}-(B^{i}u_{i})^{2},$ (8)

where $w\equiv-n_{\mu}u^{\mu}=\alpha u^{t}$ is the Lorentz factor measured by
an Eulerian observer and $B^{i}$ is the magnetic field measured in the
Eulerian frame and satisfies $B^{\mu}n_{\mu}=0$ (i.e., $B^{t}=0$). The
relation between $b^{\mu}$ and $B^{i}$ is given by

$\displaystyle
b^{t}=\frac{B^{i}u_{i}}{\alpha},~{}b_{i}=\frac{B_{i}+\left(B^{j}u_{j}\right)u_{i}}{w},$
(9)

and thus,

$\displaystyle b^{2}=\frac{B^{2}+(B^{i}u_{i})^{2}}{w^{2}},$ (10)

where $B^{2}=B_{i}B^{i}$.

The equations of motion of ideal magnetohydrodynamics and of the free-
streaming neutrino-radiation field are derived from the conservation of the
stress-energy-momentum tensor, the continuity equations for rest-mass density,
electron fraction, electron neutrino fraction, electron antineutrino fraction,
and heavy neutrino fraction, and the Maxwell equations. These conservation
laws are written as

$\displaystyle\nabla_{\mu}{\left(T_{\text{(MHD)}}\right)^{\mu}}_{\nu}=-\sum_{\nu_{i}=\nu_{e},\bar{\nu}_{e},\nu_{x}}G_{\nu}^{\text{($\nu_{i}$,leak)}},$
(11)
$\displaystyle\nabla_{\mu}{\left(T_{\text{(Rad,s,$\nu_{i}$)}}\right)^{\mu}}_{\nu}=G_{\nu}^{\text{($\nu_{i}$,leak)}},$
(12) $\displaystyle\nabla_{\mu}(\rho u^{\mu})=0,$ (13)
$\displaystyle\nabla_{\mu}(\rho u^{\mu}Y_{L})=\rho\gamma_{L},$ (14)
$\displaystyle\nabla_{\mu}{{}^{*}F}^{\mu\nu}=0,$ (15)

where $L=e,\nu_{e},\bar{\nu}_{e},$ and $\nu_{x}$ denotes electrons, electron
neutrinos, electron antineutrinos, and heavy neutrinos, respectively. $Y_{L}$
and $\gamma_{L}$ denote the fractions with respect to the baryon and the
source term for the number of the species $L$, respectively.
$G_{\nu}^{\text{($\nu_{i}$,leak)}}$ is an interaction term between the fluid
and free-streaming neutrino-radiation field of the neutrino species $\nu_{i}$
in the framework of a general relativistic neutrino leakage scheme Sekiguchi
_et al._ (2012); Sekiguchi (2010). Here ${}^{*}F^{\mu\nu}$ is the Hodge dual
of the Faraday tensor, which is given by
${}^{*}F^{\mu\nu}=b^{\mu}u^{\nu}-b^{\nu}u^{\mu}$ in ideal
magnetohydrodynamics.

Equations (11), (13), and (15) can be written in conservative form as

$\displaystyle\partial_{t}\left(\sqrt{\gamma}Q_{A}\right)+\partial_{j}\left(\sqrt{\gamma}F^{j}_{A}\right)=M_{A},$
(16) $\displaystyle\partial_{i}\left(\sqrt{\gamma}B^{i}\right)=0,$ (17)

where the flow quantities are given by the state vector
$Q_{A}=(D,J_{i},\rho_{\rm H},B^{k})^{T}$ with $A\in[0,7]$. The corresponding
fluxes are given by

$F^{j}_{A}=\\\ \left(\begin{array}[]{c}Dv^{j}\\\
\displaystyle{J_{i}v^{j}+\alpha\left(P+\frac{b^{2}}{2}\right){\delta^{j}}_{i}-\frac{\alpha}{w^{2}}B^{j}\left[B_{i}+(B^{k}u_{k})u_{i}\right]}\\\
\displaystyle{\rho_{\rm
H}v^{j}+\left(P+\frac{b^{2}}{2}\right)(v^{j}+\beta^{j})-\frac{\alpha}{w}(B^{k}u_{k})B^{j}}\\\
B^{k}v^{j}-B^{j}v^{k}\end{array}\right),$

and the source terms are

$\displaystyle M_{A}=\left(\begin{array}[]{c}0\\\ -\sqrt{\gamma}\rho_{\rm
H}\partial_{i}\alpha+\sqrt{\gamma}J_{k}\partial_{i}\beta^{k}+\frac{\alpha}{3}{S^{k}}_{k}\partial_{i}\sqrt{\gamma}-\frac{1}{2}\alpha\gamma^{1/6}S_{jk}\partial_{i}\tilde{\gamma}^{jk}-\alpha\sqrt{\gamma}G_{\text{(leak)}}^{\mu}{\gamma}_{\mu
i}\\\
\frac{\alpha}{3}\sqrt{\gamma}K{S_{k}}^{k}+\alpha\gamma^{1/6}\hat{S}_{ij}\tilde{A}^{ij}-\sqrt{\gamma}J_{k}D^{k}\alpha+\alpha\sqrt{\gamma}G_{\text{(leak)}}^{\mu}n_{\mu}\\\
0\end{array}\right),$ (22)

where $v^{j}\equiv u^{j}/u^{t}$,
$G^{\mu}_{\text{(leak)}}=\sum_{\nu_{i}=\nu_{e},\bar{\nu}_{e},\nu_{x}}G^{\mu}_{(\nu_{i},\text{leak})}$
, and the spatial components of the stress-energy-momentum tensor are given by

$\displaystyle
S_{ij}\equiv{\gamma^{\mu}}_{i}{\gamma^{\nu}}_{j}T^{\text{(MHD)}}_{\mu\nu}$
$\displaystyle=\left(\rho
h+b^{2}\right)u_{i}u_{j}+\left(P+\frac{b^{2}}{2}\right)\gamma_{ij}-b_{i}b_{j},$
(23)
$\displaystyle\hat{S}_{ij}=S_{ij}-\left(P+\frac{b^{2}}{2}\right)\gamma_{ij}.$
(24)

We also introduce the conformal metric
$\tilde{\gamma}_{ij}=\psi^{-4}\gamma_{ij}$ and the trace-free conformal
extrinsic curvature
$\tilde{A}_{ij}=\psi^{-4}\left(K_{ij}-\frac{1}{3}K\gamma_{ij}\right)$, where
$\psi$ and $K_{ij}$ are the conformal factor and the extrinsic curvature,
respectively. The explicit forms for $\gamma_{L}$ and
$G^{\text{($\nu_{i}$,leak)}}_{\mu}$ and for the equation of motion of the
free-streaming neutrino-radiation field can be found in Refs. Fujibayashi  _et
al._ (2017); Shibata (2016). The high resolution shock capturing scheme for
the neutrino-radiation field Eq. (12) is the same as that in Ref. Shibata  _et
al._ (2011).

##  III Numerical algorithm

In this section, we describe the numerical algorithms which we implemented in
our code. In Sec. III.1 we present the finite volume algorithm and
discretization scheme, and in Sec. III.2 we discuss the transformation to
Minkowski spacetime used to implement the HLLC and HLLD solvers in general
relativity. The implementation of the HLLC and HLLD solvers themselves is
presented in Sec. III.3 and III.4, respectively. Finally, the evaluation of
the electric field used by the constrained transport algorithm is discussed in
Sec. III.5.

###  III.1 Finite volume method

####  III.1.1 Fluid and magnetic field at cell center

Let $\Omega$ be a region of a given four-dimensional manifold $\cal M$,
bounded by a closed three-dimensional surface $\partial\Omega$, where
$\partial\Omega$ denotes the surface of a four-dimensional parallelepiped
composed of two spacelike surfaces $\\{\Sigma_{t},\Sigma_{t+\Delta t}\\}$ and
three sets of two timelike surfaces $\\{\Sigma_{x^{i}},\Sigma_{x^{i}+\Delta
x^{i}}\\}$ that connect the two temporal slices Font (2003). The timelike
surface, e.g., $\Sigma_{x}$, may also be regarded as a time series of
constant-$(t,x)$ surfaces, $S_{x}(t)$. We integrate Eq. (16) over the domain
of $\Omega$:

$\displaystyle\int_{\Omega}\frac{1}{\sqrt{-g}}\partial_{t}\left(\sqrt{\gamma}Q_{A}\right)d\Omega+\int_{\Omega}\frac{1}{\sqrt{-g}}\partial_{i}\left(\sqrt{\gamma}F^{i}_{A}\right)d\Omega$
$\displaystyle=\int_{\Omega}\frac{1}{\sqrt{-g}}M_{A}d\Omega,$ (25)

where $d\Omega=\sqrt{-g}dtdxdydz$.

Using Gauss’s theorem, this equation can be integrated to give

$\displaystyle(\bar{Q}_{A}\Delta V)_{t+\Delta t}-(\bar{Q}_{A}\Delta V)_{t}=$
$\displaystyle-\left(\int_{\Sigma_{x+\Delta
x}}\sqrt{\gamma}F^{x}_{A}dtdydz-\int_{\Sigma_{x}}\sqrt{\gamma}F^{x}_{A}dtdydz\right)$
$\displaystyle-\left(\int_{\Sigma_{y+\Delta
y}}\sqrt{\gamma}F^{y}_{A}dtdxdz-\int_{\Sigma_{y}}\sqrt{\gamma}F^{y}_{A}dtdxdz\right)$
$\displaystyle-\left(\int_{\Sigma_{z+\Delta
z}}\sqrt{\gamma}F^{z}_{A}dtdxdy-\int_{\Sigma_{z}}\sqrt{\gamma}F^{z}_{A}dtdxdy\right)$
$\displaystyle+\int_{\Omega}\frac{1}{\sqrt{-g}}M_{A}d\Omega,$ (26)

where

$\displaystyle\bar{Q}_{A}\equiv\frac{1}{\Delta
V}\int\sqrt{\gamma}Q_{A}dxdydz,$ (27) $\displaystyle\Delta
V\equiv\int\sqrt{\gamma}dxdydz,$ (28)

are, respectively, the three-dimensional proper volume-averaged conserved
quantities and the proper volume. Let us now define a cell consisting of
$[x_{j}-\Delta x/2:x_{j}+\Delta x/2]\times[y_{k}-\Delta y/2:y_{k}+\Delta
y/2]\times[z_{l}-\Delta z/2:z_{l}+\Delta z/2]$ (see Fig. 1). We next consider
a numerical flux, which approximates a time-averaged flux at the cell
interface and depends on the solution of the Riemann problem at the interface.
For example, in the $x$-direction the flux across the right-hand interface is
given by

$\displaystyle({\tilde{F}^{x}}_{A})_{j+\frac{1}{2},k,l}\approx\frac{1}{\Delta
t}\int^{t^{n+1}}_{t^{n}}F^{x}_{A}(Q_{A}(x_{j+1/2},y_{k},z_{l},t))dt,$ (29)

where $t^{n+1}=t^{n}+\Delta t$. With this numerical flux, Eq. (26) can be
discretized as

$\displaystyle(\bar{Q}_{A}\Delta V)^{n+1}_{j,k,l}-(\bar{Q}_{A}\Delta
V)^{n}_{j,k,l}=$ $\displaystyle-\Delta t\Big{[}\left(\Delta
A_{x}\right)_{j+\frac{1}{2},k,l}\left({\tilde{F}^{x}}_{A}\right)_{j+\frac{1}{2},k,l}$
$\displaystyle-\left(\Delta
A_{x}\right)_{j-\frac{1}{2},k,l}\left({\tilde{F}^{x}}_{A}\right)_{j-\frac{1}{2},k,l}\Big{]}$
$\displaystyle-\Delta t\Big{[}\left(\Delta
A_{y}\right)_{j,k+\frac{1}{2},l}\left({\tilde{F}^{y}}_{A}\right)_{j,k+\frac{1}{2},l}$
$\displaystyle-\left(\Delta
A_{y}\right)_{j,k-\frac{1}{2},l}\left({\tilde{F}^{y}}_{A}\right)_{j,k-\frac{1}{2},l}\Big{]}$
$\displaystyle-\Delta t\Big{[}\left(\Delta
A_{z}\right)_{j,k,l+\frac{1}{2}}\left({\tilde{F}^{z}}_{A}\right)_{j,k,l+\frac{1}{2}}$
$\displaystyle-\left(\Delta
A_{z}\right)_{j,k,l-\frac{1}{2}}\left({\tilde{F}^{z}}_{A}\right)_{j,k,l-\frac{1}{2}}\Big{]}$
$\displaystyle+\int M_{A}dtdxdydz,$ (30)

where

$\displaystyle\left(\Delta
A_{x}\right)_{j\pm\frac{1}{2},k,l}=\int\sqrt{\gamma(x_{j\pm\frac{1}{2}},y_{k},z_{l})}dydz,$
(31) $\displaystyle\left(\Delta
A_{y}\right)_{j,k\pm\frac{1}{2},l}=\int\sqrt{\gamma(x_{j},y_{k\pm\frac{1}{2}},z_{l})}dxdz,$
(32) $\displaystyle\left(\Delta
A_{z}\right)_{j,k,l\pm\frac{1}{2}}=\int\sqrt{\gamma(x_{j},y_{k},z_{l\pm\frac{1}{2}})}dxdy.$
(33)

We also assume that the determinant of the spatial metric does not change
significantly during the time step. If we introduce the volume- or surface
area-averaged determinant of the spatial metric, denoted by $\bar{\gamma}$,
this equation is reduced to

$\displaystyle(\sqrt{\bar{\gamma}}\bar{Q}_{A})^{n+1}_{j,k,l}-(\sqrt{\bar{\gamma}}\bar{Q}_{A})^{n}_{j,k,l}=$
$\displaystyle-\frac{\Delta t}{\Delta
x}\left[\left(\sqrt{\bar{\gamma}}\right)_{j+\frac{1}{2},k,l}\left({\tilde{F}^{x}}_{A}\right)_{j+\frac{1}{2},k,l}-\left(\sqrt{\bar{\gamma}}\right)_{j-\frac{1}{2},k,l}\left({\tilde{F}^{x}}_{A}\right)_{j-\frac{1}{2},k,l}\right]$
$\displaystyle-\frac{\Delta t}{\Delta
y}\left[\left(\sqrt{\bar{\gamma}}\right)_{j,k+\frac{1}{2},l}\left({\tilde{F}^{y}}_{A}\right)_{j,k+\frac{1}{2},l}-\left(\sqrt{\bar{\gamma}}\right)_{j,k-\frac{1}{2},l}\left({\tilde{F}^{y}}_{A}\right)_{j,k-\frac{1}{2},l}\right]$
$\displaystyle-\frac{\Delta t}{\Delta
z}\left[\left(\sqrt{\bar{\gamma}}\right)_{j,k,l+\frac{1}{2}}\left({\tilde{F}^{z}}_{A}\right)_{j,k,l+\frac{1}{2}}-\left(\sqrt{\bar{\gamma}}\right)_{j,k,l-\frac{1}{2}}\left({\tilde{F}^{z}}_{A}\right)_{j,k,l-\frac{1}{2}}\right]$
$\displaystyle+\left(\bar{M}_{A}\right)_{j,k,l},$ (34)

where

$\displaystyle\left(\sqrt{\bar{\gamma}}\right)_{j,k,l}\equiv\frac{1}{\Delta
x\Delta y\Delta z}\left(\Delta V\right)_{j,k,l},$ (35)
$\displaystyle\left(\sqrt{\bar{\gamma}}\right)_{j\pm\frac{1}{2},k,l}\equiv\frac{1}{\Delta
y\Delta z}\left(\Delta A_{x}\right)_{j\pm\frac{1}{2},k,l},$ (36)
$\displaystyle\left(\sqrt{\bar{\gamma}}\right)_{j,k\pm\frac{1}{2},l}\equiv\frac{1}{\Delta
x\Delta z}\left(\Delta A_{y}\right)_{j,k\pm\frac{1}{2},l},$ (37)
$\displaystyle\left(\sqrt{\bar{\gamma}}\right)_{j,k,l\pm\frac{1}{2}}\equiv\frac{1}{\Delta
x\Delta y}\left(\Delta A_{z}\right)_{j,k,l\pm\frac{1}{2}},$ (38)

and

$\displaystyle\left(\bar{M}_{A}\right)_{j,k,l}\equiv\frac{1}{\Delta x\Delta
y\Delta z}\int M_{A}dtdxdydz.$ (39)

####  III.1.2 Magnetic fields at cell surface

To ensure that the divergence-free condition (17) is maintained, we employ the
constrained transport method introduced by Evans and Hawley Evans and Hawley
(1988). In this method, the magnetic-field components are defined at the cell
surfaces, and the electric field components are defined at the cell edges (see
Fig. 1).

We then integrate Eq. (16) for $A\in[5,7]$ on $\Sigma_{x^{i}+\Delta x^{i}}$.
For example, through the surface $\Sigma_{z+\Delta z}$, we have

$\displaystyle\int_{\Sigma_{z+\Delta
z}}\frac{1}{\sqrt{-g}}\partial_{t}\left(\sqrt{\gamma}B^{z}\right)dS_{\Omega_{z}}$
$\displaystyle+\int_{\Sigma_{z+\Delta
z}}\frac{1}{\sqrt{-g}}\partial_{j}\left(\sqrt{\gamma}\epsilon^{zjk}E_{k}\right)dS_{\Omega_{z}}=0,$
(40)

where $E_{k}=-\epsilon_{kij}v^{i}B^{j}$, $\epsilon_{ijk}$ is the three-
dimensional Levi-Civita tensor, and $dS_{\Omega_{z}}=\sqrt{-g}dtdxdy$. Using
Stokes’ theorem, this equation is integrated to give

$\displaystyle\left(\bar{B}^{z}\Delta A_{z}\right)_{t+\Delta
t}-\left(\bar{B}^{z}\Delta A_{z}\right)_{t}$
$\displaystyle=-\int_{t}^{t+\Delta t}\oint_{\partial S_{z+\Delta
z}}\sqrt{\gamma}E_{i}dx^{i}dt,$ (41)

where

$\displaystyle\bar{B}^{z}\equiv\frac{1}{\Delta A_{z}}\int_{S_{z+\Delta
z}}\sqrt{\gamma}B^{z}dxdy,$

is the surface-averaged magnetic field. Similarly, through the surfaces
$\Sigma_{x+\Delta x}$ and $\Sigma_{y+\Delta y}$, respectively, we have

$\displaystyle\left(\bar{B}^{x}\Delta A_{x}\right)_{t+\Delta
t}-\left(\bar{B}^{x}\Delta A_{x}\right)_{t}$
$\displaystyle=-\int_{t}^{t+\Delta t}\oint_{\partial S_{x+\Delta
x}}\sqrt{\gamma}E_{i}dx^{i}dt,$ (42) $\displaystyle\left(\bar{B}^{y}\Delta
A_{y}\right)_{t+\Delta t}-\left(\bar{B}^{y}\Delta A_{y}\right)_{t}$
$\displaystyle=-\int_{t}^{t+\Delta t}\oint_{\partial S_{y+\Delta
y}}\sqrt{\gamma}E_{i}dx^{i}dt,$ (43)

where

$\displaystyle\bar{B}^{x}\equiv\frac{1}{\Delta A_{x}}\int_{S_{x+\Delta
x}}\sqrt{\gamma}B^{x}dydz,$ $\displaystyle\bar{B}^{y}\equiv\frac{1}{\Delta
A_{y}}\int_{S_{y+\Delta y}}\sqrt{\gamma}B^{y}dxdz.$

We next consider a cell surface consisting of $[x_{j}-\Delta x/2:x_{j}+\Delta
x/2]\times[y_{k}-\Delta y/2:y_{k}+\Delta y/2]$, $[y_{k}-\Delta
y/2:y_{k}+\Delta y/2]\times[z_{l}-\Delta z/2:z_{l}+\Delta z/2]$,
$[x_{j}-\Delta x/2:x_{j}+\Delta x/2]\times[z_{l}-\Delta z/2:z_{l}+\Delta z/2]$
and a numerical flux which approximates a time-averaged electric field at the
cell edge, given by

$\displaystyle\left(\tilde{E}_{x}\right)_{j,k+\frac{1}{2},l+\frac{1}{2}}\approx\frac{1}{\Delta
t}\int^{t^{n+1}}_{t^{n}}E_{x}(Q_{A}(x_{j},y_{k+1/2},z_{l+1/2}))dt,$ (44)
$\displaystyle\left(\tilde{E}_{y}\right)_{j+\frac{1}{2},k,l+\frac{1}{2}}\approx\frac{1}{\Delta
t}\int^{t^{n+1}}_{t^{n}}E_{y}(Q_{A}(x_{j+1/2},y_{k},z_{l+1/2}))dt,$ (45)
$\displaystyle\left(\tilde{E}_{z}\right)_{j+\frac{1}{2},k+\frac{1}{2},l}\approx\frac{1}{\Delta
t}\int^{t^{n+1}}_{t^{n}}E_{z}(Q_{A}(x_{j+1/2},y_{k+1/2},z_{l}))dt.$ (46)

With these averaged electric fields, Eqs. (III.1.2)–(43) are discretized as

$\displaystyle\left(\bar{B}^{x}\Delta
A_{x}\right)^{n+1}_{j+\frac{1}{2},k,l}-\left(\bar{B}^{x}\Delta
A_{x}\right)^{n}_{j+\frac{1}{2},k,l}$ $\displaystyle=\Delta
t\Big{[}\left(\Delta
l_{y}\right)_{j+\frac{1}{2},k,l+\frac{1}{2}}\left(\tilde{E}_{y}\right)_{j+\frac{1}{2},k,l+\frac{1}{2}}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}-\left(\Delta
l_{y}\right)_{j+\frac{1}{2},k,l-\frac{1}{2}}\left(\tilde{E}_{y}\right)_{j+\frac{1}{2},k,l-\frac{1}{2}}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}-\left(\Delta
l_{z}\right)_{j+\frac{1}{2},k+\frac{1}{2},l}\left(\tilde{E}_{z}\right)_{j+\frac{1}{2},k+\frac{1}{2},l}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}+\left(\Delta
l_{z}\right)_{j+\frac{1}{2},k-\frac{1}{2},l}\left(\tilde{E}_{z}\right)_{j+\frac{1}{2},k-\frac{1}{2},l}\Big{]},$
(47) $\displaystyle\left(\bar{B}^{y}\Delta
A_{y}\right)^{n+1}_{j,k+\frac{1}{2},l}-\left(\bar{B}^{y}\Delta
A_{y}\right)^{n}_{j,k+\frac{1}{2},l}$ $\displaystyle=\Delta
t\Big{[}\left(\Delta
l_{z}\right)_{j+\frac{1}{2},k+\frac{1}{2},l}\left(\tilde{E}_{z}\right)_{j+\frac{1}{2},k+\frac{1}{2},l}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}-\left(\Delta
l_{z}\right)_{j-\frac{1}{2},k+\frac{1}{2},l}\left(\tilde{E}_{z}\right)_{j-\frac{1}{2},k+\frac{1}{2},l}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}-\left(\Delta
l_{x}\right)_{j,k+\frac{1}{2},l+\frac{1}{2}}\left(\tilde{E}_{x}\right)_{j,k+\frac{1}{2},l+\frac{1}{2}}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}+\left(\Delta
l_{x}\right)_{j,k+\frac{1}{2},l-\frac{1}{2}}\left(\tilde{E}_{x}\right)_{j,k+\frac{1}{2},l-\frac{1}{2}}\Big{]},$
(48) $\displaystyle\left(\bar{B}^{z}\Delta
A_{z}\right)^{n+1}_{j,k,l+\frac{1}{2}}-\left(\bar{B}^{z}\Delta
A_{z}\right)^{n}_{j,k,l+\frac{1}{2}}$ $\displaystyle=\Delta
t\Big{[}\left(\Delta
l_{x}\right)_{j,k+\frac{1}{2},l+\frac{1}{2}}\left(\tilde{E}_{x}\right)_{j,k+\frac{1}{2},l+\frac{1}{2}}$
$\displaystyle~{}~{}~{}~{}~{}-\left(\Delta
l_{x}\right)_{j,k-\frac{1}{2},l+\frac{1}{2}}\left(\tilde{E}_{x}\right)_{j,k-\frac{1}{2},l+\frac{1}{2}}$
$\displaystyle~{}~{}~{}~{}~{}-\left(\Delta
l_{y}\right)_{j+\frac{1}{2},k,l+\frac{1}{2}}\left(\tilde{E}_{y}\right)_{j+\frac{1}{2},k,l+\frac{1}{2}}$
$\displaystyle~{}~{}~{}~{}~{}+\left(\Delta
l_{y}\right)_{j-\frac{1}{2},k,l+\frac{1}{2}}\left(\tilde{E}_{y}\right)_{j-\frac{1}{2},k,l+\frac{1}{2}}\Big{]},$
(49)

where

$\displaystyle\left(\Delta
l_{x}\right)_{j,k\pm\frac{1}{2},l\pm\frac{1}{2}}=\int\sqrt{\gamma(x_{j},y_{k\pm
1/2},z_{l\pm 1/2})}dx,$ (50) $\displaystyle\left(\Delta
l_{y}\right)_{j\pm\frac{1}{2},k,l\pm\frac{1}{2}}=\int\sqrt{\gamma(x_{j\pm
1/2},y_{k},z_{l\pm 1/2})}dy,$ (51) $\displaystyle\left(\Delta
l_{z}\right)_{j\pm\frac{1}{2},k\pm\frac{1}{2},l}=\int\sqrt{\gamma(x_{j\pm
1/2},y_{k\pm 1/2},z_{l})}dz.$ (52)

If we introduce a line-averaged determinant of the spatial metric by

$\displaystyle\left(\sqrt{\bar{\gamma}}\right)_{j,k\pm\frac{1}{2},l\pm\frac{1}{2}}\equiv\frac{\left(\Delta
l_{x}\right)_{j,k\pm\frac{1}{2},l\pm\frac{1}{2}}}{\Delta x},$
$\displaystyle\left(\sqrt{\bar{\gamma}}\right)_{j\pm\frac{1}{2},k,l\pm\frac{1}{2}}\equiv\frac{\left(\Delta
l_{y}\right)_{j\pm\frac{1}{2},k,l\pm\frac{1}{2}}}{\Delta y},$
$\displaystyle\left(\sqrt{\bar{\gamma}}\right)_{j\pm\frac{1}{2},k\pm\frac{1}{2},l}\equiv\frac{\left(\Delta
l_{z}\right)_{j\pm\frac{1}{2},k\pm\frac{1}{2},l}}{\Delta z},$ (53)

then, together with the surface area-averaged spatial metric given by Eqs.
(36)–(38), Eqs. (47)–(49) are reduced to

$\displaystyle\left(\sqrt{\bar{\gamma}}\bar{B}^{x}\right)^{n+1}_{j+\frac{1}{2},k,l}-\left(\sqrt{\bar{\gamma}}\bar{B}^{x}\right)^{n}_{j+\frac{1}{2},k,l}$
$\displaystyle=\frac{\Delta t}{\Delta
z}\Big{[}\left(\sqrt{\bar{\gamma}}\right)_{j+\frac{1}{2},k,l+\frac{1}{2}}\left(\tilde{E}_{y}\right)_{j+\frac{1}{2},k,l+\frac{1}{2}}$
$\displaystyle~{}~{}~{}~{}~{}~{}-\left(\sqrt{\bar{\gamma}}\right)_{j+\frac{1}{2},k,l-\frac{1}{2}}\left(\tilde{E}_{y}\right)_{j+\frac{1}{2},k,l-\frac{1}{2}}\Big{]}$
$\displaystyle-\frac{\Delta t}{\Delta
y}\Big{[}\left(\sqrt{\bar{\gamma}}\right)_{j+\frac{1}{2},k+\frac{1}{2},l}\left(\tilde{E}_{z}\right)_{j+\frac{1}{2},k+\frac{1}{2},l}$
$\displaystyle~{}~{}~{}~{}~{}~{}-\left(\sqrt{\bar{\gamma}}\right)_{j+\frac{1}{2},k-\frac{1}{2},l}\left(\tilde{E}_{z}\right)_{j+\frac{1}{2},k-\frac{1}{2},l}\Big{]},$
(54)
$\displaystyle\left(\sqrt{\bar{\gamma}}\bar{B}^{y}\right)^{n+1}_{j,k+\frac{1}{2},l}-\left(\sqrt{\bar{\gamma}}\bar{B}^{y}\right)^{n}_{j,k+\frac{1}{2},l}$
$\displaystyle=\frac{\Delta t}{\Delta
x}\Big{[}\left(\sqrt{\bar{\gamma}}\right)_{j+\frac{1}{2},k+\frac{1}{2},l}\left(\tilde{E}_{z}\right)_{j+\frac{1}{2},k+\frac{1}{2},l}$
$\displaystyle~{}~{}~{}~{}~{}~{}-\left(\sqrt{\bar{\gamma}}\right)_{j-\frac{1}{2},k+\frac{1}{2},l}\left(\tilde{E}_{z}\right)_{j-\frac{1}{2},k+\frac{1}{2},l}\Big{]},$
$\displaystyle-\frac{\Delta t}{\Delta
z}\Big{[}\left(\sqrt{\bar{\gamma}}\right)_{j,k+\frac{1}{2},l+\frac{1}{2}}\left(\tilde{E}_{x}\right)_{j,k+\frac{1}{2},l+\frac{1}{2}}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}-\left(\sqrt{\bar{\gamma}}\right)_{j,k+\frac{1}{2},l-\frac{1}{2}}\left(\tilde{E}_{x}\right)_{j,k+\frac{1}{2},l-\frac{1}{2}}\Big{]},$
(55)
$\displaystyle\left(\sqrt{\bar{\gamma}}\bar{B}^{z}\right)^{n+1}_{j,k,l+\frac{1}{2}}-\left(\sqrt{\bar{\gamma}}\bar{B}^{z}\right)^{n}_{j,k,l+\frac{1}{2}}$
$\displaystyle=\frac{\Delta t}{\Delta
y}\Big{[}\left(\sqrt{\bar{\gamma}}\right)_{j,k+\frac{1}{2},l+\frac{1}{2}}\left(\tilde{E}_{x}\right)_{j,k+\frac{1}{2},l+\frac{1}{2}}$
$\displaystyle~{}~{}~{}~{}~{}~{}-\left(\sqrt{\bar{\gamma}}\right)_{j,k-\frac{1}{2},l+\frac{1}{2}}\left(\tilde{E}_{x}\right)_{j,k-\frac{1}{2},l+\frac{1}{2}}\Big{]}$
$\displaystyle-\frac{\Delta t}{\Delta
x}\Big{[}\left(\sqrt{\bar{\gamma}}\right)_{j+\frac{1}{2},k,l+\frac{1}{2}}\left(\tilde{E}_{y}\right)_{j+\frac{1}{2},k,l+\frac{1}{2}}$
$\displaystyle~{}~{}~{}~{}~{}~{}-\left(\sqrt{\bar{\gamma}}\right)_{j-\frac{1}{2},k,l+\frac{1}{2}}\left(\tilde{E}_{y}\right)_{j-\frac{1}{2},k,l+\frac{1}{2}}\Big{]}.$
(56)

The magnetic-field distribution inside the cell is reconstructed from the
magnetic fields at the cell surface. Practically, we reconstruct the magnetic
field at the cell center in Eq. (34) by

$\displaystyle\left(\bar{B}^{x}\right)_{j,k,l}=\frac{1}{2}\left[\left(\bar{B}^{x}\right)_{j+\frac{1}{2},k,l}+\left(\bar{B}^{x}\right)_{j-\frac{1}{2},k,l}\right],$
(57)
$\displaystyle\left(\bar{B}^{y}\right)_{j,k,l}=\frac{1}{2}\left[\left(\bar{B}^{y}\right)_{j,k+\frac{1}{2},l}+\left(\bar{B}^{y}\right)_{j,k-\frac{1}{2},l}\right],$
(58)
$\displaystyle\left(\bar{B}^{z}\right)_{j,k,l}=\frac{1}{2}\left[\left(\bar{B}^{z}\right)_{j,k,l+\frac{1}{2}}+\left(\bar{B}^{z}\right)_{j,k,l-\frac{1}{2}}\right].$
(59)

###  III.2 Tetrad frame

To evaluate the numerical fluxes through cell interfaces (e.g. Eq. (29)), we
implement HLL-type Riemann solvers Mignone and Bodo (2005); Mignone  _et al._
(2009). Because these Riemann solvers are designed to solve a Riemann problem
in Minkowski spacetime (except for the HLLE solver, which we have implemented
directly in curved spacetime, see, e.g., Ref. Shibata and Sekiguchi (2005)),
it is necessary to transform all the equations into a tetrad frame in order to
apply these methods to a general relativistic framework.

Following Ref. White  _et al._ (2016), we define a tetrad basis in the
$x$-direction, for example, by

$\displaystyle{e_{(\hat{t})}}^{\mu}=n^{\mu},$ (60)
$\displaystyle{e_{(\hat{x})}}^{\mu}=\hat{B}\left(0,\gamma^{xi}\right),$ (61)
$\displaystyle{e_{(\hat{y})}}^{\mu}=\hat{D}\left(0,0,\gamma_{zz},-\gamma_{yz}\right),$
(62) $\displaystyle{e_{(\hat{z})}}^{\mu}=\hat{C}\left(0,0,0,1\right),$ (63)

where

$\displaystyle\hat{B}=\frac{1}{\sqrt{\gamma^{xx}}},$ (64)
$\displaystyle\hat{C}=\frac{1}{\sqrt{\gamma_{zz}}},$ (65)
$\displaystyle\hat{D}=\frac{1}{\sqrt{\gamma_{zz}\left(\gamma_{yy}\gamma_{zz}-\gamma_{yz}^{2}\right)}}.$
(66)

With this basis, we can perform a transformation from the Eulerian frame to
the tetrad frame by

$\displaystyle V_{(\hat{\mu})}={e_{(\hat{\mu})}}^{\mu}V_{\mu},$ (67)
$\displaystyle
Q_{(\hat{\mu})(\hat{\nu})}={e_{(\hat{\mu})}}^{\mu}{e_{(\hat{\nu})}}^{\nu}Q_{\mu\nu},$
(68)

where $V_{\mu}$ and $Q_{\mu\nu}$ denote a covariant vector and tensor,
respectively, in the Eulerian frame. The covariant components of the tetrad
basis are

$\displaystyle e_{(\hat{t})\mu}$ $\displaystyle=n_{\mu},$ (69) $\displaystyle
e_{(\hat{x})\mu}$
$\displaystyle=\hat{B}\left(\beta^{x},{\delta_{i}}^{x}\right),$ (70)
$\displaystyle e_{(\hat{y})\mu}$
$\displaystyle=\hat{D}\Big{(}\beta_{y}\gamma_{zz}-\beta_{z}\gamma_{yz},\gamma_{xy}\gamma_{zz}-\gamma_{xz}\gamma_{yz},$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\gamma_{yy}\gamma_{zz}-\gamma_{yz}^{2},0\Big{)},$
(71) $\displaystyle e_{(\hat{z})\mu}$
$\displaystyle=\hat{C}\left(\beta_{z},\gamma_{iz}\right).$ (72)

With this basis, we can then perform the transformation from the tetrad frame
to the Eulerian frame by

$\displaystyle V_{\mu}=e_{(\hat{\mu})\mu}V^{(\hat{\mu})},$ (73) $\displaystyle
Q_{\mu\nu}=e_{(\hat{\mu})\mu}e_{(\hat{\nu})\nu}Q^{(\hat{\mu})(\hat{\nu})}.$
(74)

With this tetrad basis the procedure to obtain the numerical flux
$\left({\tilde{F}^{x}}_{A}\right)_{j+\frac{1}{2},k,l}$ is as follows: first,
we calculate the tetrad component of $u_{(\hat{\imath})},v^{(\hat{\imath})}$,
and $B^{(\hat{\imath})}$ by

$\displaystyle
u_{(\hat{\imath})}=e_{(\hat{\imath})\mu}u^{\mu}=\frac{w}{\alpha}\left(e_{(\hat{\imath})t}+e_{(\hat{\imath})j}v^{j}\right),$
(75) $\displaystyle
v^{(\hat{\imath})}\equiv\frac{u^{(\hat{\imath})}}{u^{(\hat{t})}}=\frac{{e^{(\hat{\imath})}}_{\mu}u^{\mu}}{{e^{(\hat{t})}}_{\nu}u^{\nu}}=\frac{e_{(\hat{\imath})t}+e_{(\hat{\imath})j}v^{j}}{\alpha},$
(76) $\displaystyle
B^{(\hat{\imath})}={e^{(\hat{\imath})}}_{\mu}B^{\mu}={e^{(\hat{\imath})}}_{j}B^{j}.$
(77)

Second, we solve a Riemann problem in the locally Minkowski spacetime to
obtain the numerical flux
$\left({\tilde{f}^{(\hat{x})}}_{A}\right)_{j+\frac{1}{2},k,l}$ and the
conserved quantities $\left(q_{A}\right)_{j+\frac{1}{2},k,l}$ at the cell
interface (see the next section for more detail on the Riemann problem).
Finally, we transform back to the Eulerian frame from the tetrad frame by

$\displaystyle({\tilde{F}^{x}}_{0})_{j+\frac{1}{2},k,l}=\left(Dv^{x}\right)_{j+\frac{1}{2},k,l}$
$\displaystyle=\left(\alpha\left({e_{(\hat{t})}}^{x}D+{e_{(\hat{x})}}^{x}\tilde{f}^{(\hat{x})}_{0}\right)\right)_{j+\frac{1}{2},k,l},$
(78)
$\displaystyle({\tilde{F}^{x}}_{1})_{j+\frac{1}{2},k,l}=\left(\alpha{T^{x}}_{x}\right)_{j+\frac{1}{2},k,l}$
$\displaystyle=\Big{(}\alpha\Big{(}{e_{(\hat{t})}}^{x}e_{(\hat{\imath})x}J_{(\hat{\imath})}+{e_{(\hat{x})}}^{x}e_{(\hat{\imath})x}\tilde{f}^{(\hat{x})}_{i}\Big{)}\Big{)}_{j+\frac{1}{2},k,l},$
(79)
$\displaystyle({\tilde{F}^{x}}_{2})_{j+\frac{1}{2},k,l}=\left(\alpha{T^{x}}_{y}\right)_{j+\frac{1}{2},k,l}$
$\displaystyle=\Big{(}\alpha\Big{(}{e_{(\hat{t})}}^{x}e_{(\hat{\imath})y}J_{(\hat{\imath})}+{e_{(\hat{x})}}^{x}e_{(\hat{\imath})y}\tilde{f}^{(\hat{x})}_{i}\Big{)}\Big{)}_{j+\frac{1}{2},k,l},$
(80)
$\displaystyle({\tilde{F}^{x}}_{3})_{j+\frac{1}{2},k,l}=\left(\alpha{T^{x}}_{z}\right)_{j+\frac{1}{2},k,l}$
$\displaystyle=\Big{(}\alpha\Big{(}{e_{(\hat{t})}}^{x}e_{(\hat{z})z}J_{(\hat{z})}+{e_{(\hat{x})}}^{x}e_{(\hat{z})z}\tilde{f}^{(\hat{x})}_{3}\Big{)}\Big{)}_{j+\frac{1}{2},k,l},$
(81)
$\displaystyle({\tilde{F}^{x}}_{4})_{j+\frac{1}{2},k,l}=\left(-\alpha{T^{x}}_{\mu}n^{\mu}\right)_{j+\frac{1}{2},k,l}$
$\displaystyle=\left(\alpha\left({e_{(\hat{t})}}^{x}\rho_{\rm
H}+{e_{(\hat{x})}}^{x}\tilde{f}^{(\hat{x})}_{4}\right)\right)_{j+\frac{1}{2},k,l},$
(82) $\displaystyle(\tilde{F}^{x}_{5})_{j+\frac{1}{2},k,l}=0,$ (83)
$\displaystyle(\tilde{F}^{x}_{6})_{j+\frac{1}{2},k,l}=\left(-\tilde{E}_{z}\right)_{j+\frac{1}{2},k,l}=\left(\alpha{{}^{*}F}^{yx}\right)_{j+\frac{1}{2},k,l}$
$\displaystyle=\Big{(}\alpha\Big{(}{e_{(\hat{\imath})}}^{y}{e_{(\hat{t})}}^{x}\bar{B}^{(\hat{\imath})}-{e_{(\hat{t})}}^{y}{e_{(\hat{x})}}^{x}\bar{B}^{(\hat{x})}$
$\displaystyle~{}~{}~{}~{}~{}+{e_{(\hat{y})}}^{y}{e_{(\hat{x})}}^{x}\tilde{f}^{(\hat{x})}_{6}\Big{)}\Big{)}_{j+\frac{1}{2},k,l},$
(84)
$\displaystyle(\tilde{F}^{x}_{7})_{j+\frac{1}{2},k,l}=\left(\tilde{E}_{y}\right)_{j+\frac{1}{2},k,l}=\left(\alpha{{}^{*}F}^{zx}\right)_{j+\frac{1}{2},k,l}$
$\displaystyle=\Big{(}\alpha\Big{(}{e_{(\hat{\imath})}}^{z}{e_{(\hat{t})}}^{x}\bar{B}^{(\hat{\imath})}-{e_{(\hat{t})}}^{z}{e_{(\hat{x})}}^{x}\bar{B}^{(\hat{x})}$
$\displaystyle~{}~{}~{}~{}~{}+{e_{(\hat{y})}}^{z}{e_{(\hat{x})}}^{x}\tilde{f}^{(\hat{x})}_{6}+{e_{(\hat{z})}}^{z}{e_{(\hat{x})}}^{x}\tilde{f}^{(\hat{x})}_{7}\Big{)}\Big{)}_{j+\frac{1}{2},k,l},$
(85)

where $\hat{\imath}=\hat{x},\hat{y},\hat{z}$ are contracted with $i=1,2,3$,
respectively, in the second term of the right-hand side of Eqs. (79) and (80).
Note that, from now on, we do not distinguish the upper- and lower-spatial
tetrad components, e.g., $B^{(\hat{\imath})}=B_{(\hat{\imath})}$. These
numerical fluxes are used to update the conserved quantities in Eq. (34). An
interface velocity is calculated by White  _et al._ (2016)

$\displaystyle
v^{(\hat{x})}_{\text{interface}}=\frac{d\hat{x}}{d\hat{t}}=\frac{\beta^{x}}{\alpha\sqrt{\gamma^{xx}}}.$
(86)

This velocity is used to calculate a numerical flux at the interface (see Eqs.
(103) and (137) in the next section). The tetrad basis and numerical fluxes in
the $y$\- and $z$-directions are summarized in Appendix A.

Figure 1: Schematic of a cell, cell interface, and cell edge for the finite
volume method with the constrained transport method. Fluid quantities,
$(\bar{Q}_{A})_{j,k,l}$, are defined at the cell center. The magnetic field
components,
$(\bar{B}^{x})_{j+\frac{1}{2},k,l},(\bar{B}^{y})_{j,k+\frac{1}{2},l},(\bar{B}^{z})_{j,k,l+\frac{1}{2}}$,
are defined at the cell interfaces. The electric field components,
$(\tilde{E}^{x})_{j,k+\frac{1}{2},l+\frac{1}{2}},(\tilde{E}^{y})_{j+\frac{1}{2},k,l+\frac{1}{2}},(\tilde{E}^{z})_{j+\frac{1}{2},k+\frac{1}{2},l}$,
are defined at the cell edges.  Figure 2: Riemann fan structure for the HLLC
solver for relativistic hydrodynamics (left), and for the HLLD solver for
relativistic magnetohydrodynamics (right) in the tetrad frame. In the HLLC
solver (left panel), the left-going nonlinear wave with $\lambda_{L}$, the
contact discontinuity with $\lambda_{c}$, and the right-going nonlinear wave
with $\lambda_{R}$, propagate from the discontinuity located at
$\hat{x}_{j+\frac{1}{2}}$, where $\lambda_{L,c,R}$ denotes the characteristic
speed of each wave. Consequently, the $L$, $cL$, $cR$, and $R$ states appear.
In the HLLD solver (right panel), the left/right-propagating fast wave with
$\lambda_{L}$/$\lambda_{R}$, the left/right-propagating Alfvén wave with
characteristic speed $\lambda_{aL}$/$\lambda_{aR}$, and the contact
discontinuity with $\lambda_{c}$, are taken into account. Consequently, the
$L$, $aL$, $cL$, $cR$, $aR$, and $R$ states appear. In the general
relativistic case, the interface initially located at
$\hat{x}_{j+\frac{1}{2}}$ may move with an interface velocity
$v^{(\hat{x})}_{\text{interface}}$ which is proportional to the shift vector
$\beta^{x}$.

###  III.3 HLLC solver for relativistic hydrodynamics

In the absence of electromagnetic fields, Eq. (34) with $A\in[0,4]$ are
reduced to those of relativistic hydrodynamics. In this case, one choice for
the Riemann solver is the HLLC solver proposed in Ref. Mignone and Bodo
(2005). We calculate the HLLC flux
$\left(\tilde{f}^{(\hat{x})}_{A}\right)_{j+\frac{1}{2},k,l}$ in the tetrad
frame by solving the source-free one-dimensional conservation law:

$\displaystyle\partial_{(\hat{t})}q_{A}+\partial_{(\hat{x})}f^{(\hat{x})}_{A}=0,$
(87) $\displaystyle q_{A}=\left(\begin{array}[]{c}D\\\ J_{(\hat{\imath})}\\\
\rho_{\rm H}\end{array}\right),$ (91) $\displaystyle
f^{(\hat{x})}_{A}=\left(\begin{array}[]{c}Dv^{(\hat{x})}\\\
J_{(\hat{\imath})}v^{(\hat{x})}+P{\delta^{(\hat{x})}}_{(\hat{\imath})}\\\
\rho_{\rm H}v^{(\hat{x})}+Pv^{(\hat{x})}\end{array}\right),$ (95)

where $\partial_{(\hat{\mu})}\equiv{e_{(\hat{\mu})}}^{\mu}\partial_{\mu}$.
Given an initial condition described by

$\displaystyle
q_{A}(x,0)=\left\\{\begin{array}[]{l}\left(q_{A}\right)_{L}~{}\text{if
$x<x_{j+\frac{1}{2}}$},\\\ \left(q_{A}\right)_{R}~{}\text{if
$x>x_{j+\frac{1}{2}}$},\end{array}\right.$ (98)

for $x_{j}\leq x\leq x_{j+1}$, three characteristic speeds and therefore four
states will appear in the Riemann fan (see Fig. 2). In the HLLC solver, one
needs to find the pressure in the intermediate states (the $cL$ and $cR$
states) which satisfies a jump condition. Then, the numerical flux is
calculated by (see the left panel of Fig. 2)

$\displaystyle\left(\tilde{f}^{(\hat{x})}_{A}\right)_{j+\frac{1}{2}}$
$\displaystyle=\left\\{\begin{array}[]{cll}&(f^{(\hat{x})}_{A})_{L}&\text{if
}\lambda_{L}>v^{(\hat{x})}_{\text{interface}}\\\
&(f^{(\hat{x})}_{A})_{cL}&\text{if
}\lambda_{L}<v^{(\hat{x})}_{\text{interface}}<\lambda_{c}\\\
&(f^{(\hat{x})}_{A})_{cR}&\text{if
}\lambda_{c}<v^{(\hat{x})}_{\text{interface}}<\lambda_{R}\\\
&(f^{(\hat{x})}_{A})_{R}&\text{if
}\lambda_{R}<v^{(\hat{x})}_{\text{interface}},\\\ \end{array}\right.$ (103)

where

$\displaystyle\left(f^{(\hat{x})}_{A}\right)_{L/R}=f^{(\hat{x})}_{A}\left(q_{L/R}\right),$
(104)
$\displaystyle\left(f^{(\hat{x})}_{A}\right)_{cL/cR}=\left(f^{(\hat{x})}_{A}\right)_{L/R}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\lambda_{L/R}\left(\left(q_{A}\right)_{cL/cR}-\left(q_{A}\right)_{L/R}\right),$
(105)

and $\lambda_{L/R}$ is the characteristic speed of the left/right-going
nonlinear wave. Equation (105) is obtained from the jump condition and
$\lambda_{c}$ is the characteristic speed of the contact discontinuity. By
imposing continuity of the pressure across the contact discontinuity, one
finds a quadratic equation for $\lambda_{c}$ Mignone and Bodo (2005):

$\displaystyle F^{\text{HLL}}_{\rho_{\rm
H}}\lambda_{c}^{2}-\left(\rho^{\text{HLL}}_{\rm
H}+F^{\text{HLL}}_{J_{(\hat{x})}}\right)\lambda_{c}+J^{\text{HLL}}_{(\hat{x})}=0,$
(106)

where $\rho^{\text{HLL}}_{\rm H}$, $J^{\text{HLL}}_{(\hat{x})}$,
$F^{\text{HLL}}_{\rho_{\rm H}}$, and $F^{\text{HLL}}_{J_{(\hat{x})}}$ denote
conserved quantities and fluxes in the HLL state:

$\displaystyle
J^{\text{HLL}}_{(\hat{x})}=\frac{\lambda_{R}J_{(\hat{x})}^{R}-\lambda_{L}J_{(\hat{x})}^{L}+f^{(\hat{x})}_{1,L}-f^{(\hat{x})}_{1,R}}{\lambda_{R}-\lambda_{L}},$
(107) $\displaystyle\rho^{\text{HLL}}_{\rm H}=\frac{\lambda_{R}\rho_{\rm
H}^{R}-\lambda_{L}\rho_{\rm
H}^{L}+f^{(\hat{x})}_{4,L}-f^{(\hat{x})}_{4,R}}{\lambda_{R}-\lambda_{L}},$
(108) $\displaystyle
F^{\text{HLL}}_{J_{(\hat{x})}}=\frac{\lambda_{R}f^{(\hat{x})}_{1,L}-\lambda_{L}f^{(\hat{x})}_{1,R}+\lambda_{R}\lambda_{L}\left(J_{(\hat{x})}^{R}-J_{(\hat{x})}^{L}\right)}{\lambda_{R}-\lambda_{L}},$
(109) $\displaystyle F^{\text{HLL}}_{\rho_{\rm
H}}=\frac{\lambda_{R}f^{(\hat{x})}_{4,L}-\lambda_{L}f^{(\hat{x})}_{4,R}+\lambda_{R}\lambda_{L}\left(\rho_{\rm
H}^{R}-\rho_{\rm H}^{L}\right)}{\lambda_{R}-\lambda_{L}}.$ (110)

Once we obtain the speed of the contact discontinuity $\lambda_{c}$, the
pressure in the intermediate state is determined by

$\displaystyle P_{c}\equiv P_{cL}=P_{cR}=-\lambda_{c}F^{\text{HLL}}_{\rho_{\rm
H}}+F^{\text{HLL}}_{J_{(\hat{x})}}.$ (111)

Then the conserved quantities in the $cL$ and $cR$ states are given by

$\displaystyle
D_{cL/cR}=\frac{D_{L/R}\left(\lambda_{L/R}-v^{(\hat{x})}_{L/R}\right)}{\lambda_{L/R}-\lambda_{c}},$
(112)
$\displaystyle\left(J_{(\hat{\imath})}\right)_{cL/cR}=\frac{1}{\lambda_{L/R}-\lambda_{c}}$
$\displaystyle\times\left[\left(J_{(\hat{\imath})}\right)_{L/R}\left(\lambda_{L/R}-v^{(\hat{x})}_{L/R}\right)+(P_{c}-P_{L/R}){\delta^{(\hat{x})}}_{(\hat{\imath})}\right],$
(113) $\displaystyle\left(\rho_{\rm H}\right)_{cL/cR}=\frac{\left(\rho_{\rm
H}\right)_{L/R}\left(\lambda_{L/R}-v^{(\hat{x})}_{L/R}\right)+P_{c}\lambda_{c}-P_{L/R}v^{(\hat{x})}_{L/R}}{\lambda_{L/R}-\lambda_{c}},$
(114)

where the subscripts $cL$ and $cR$ on the left-hand side of the equations
correspond to $L$ and $R$ on the right-hand side, respectively. These
quantities can be used to evaluate the flux in the $cL/cR$ state (105) and the
flux in the Eulerian frame (see, e.g., Eq. (78)).

For the left and right characteristic speeds $\lambda_{L/R}$, we apply Davis’s
estimate Mignone and Bodo (2005):

$\displaystyle\lambda_{L}=\min(\lambda^{-}(q_{L}),\lambda^{-}(q_{R})),$ (115)
$\displaystyle\lambda_{R}=\max(\lambda^{+}(q_{L}),\lambda^{+}(q_{R})),$ (116)

and

$\displaystyle\lambda^{\pm}(q_{A})=\frac{1}{1-v^{2}c_{s}^{2}}\left[v^{(\hat{x})}\left(1-c_{s}^{2}\right)\right.$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\left.\pm
c_{s}\sqrt{(1-v^{2})(1-v^{2}c_{s}^{2}-(1-c_{s}^{2})(v^{(\hat{x})})^{2})}\right],$
(117)

where

$\displaystyle v^{2}=v^{(\hat{\imath})}v_{(\hat{\imath})},$ (118)
$\displaystyle c_{s}^{2}=\frac{1}{h}\left[\frac{\partial
P}{\partial\rho}\Big{|}_{\varepsilon}+\frac{P}{\rho^{2}}\frac{\partial
P}{\partial\varepsilon}\Big{|}_{\rho}\right].$ (119)

The equivalent expressions in the $y$\- and $z$-directions are given by
permutation of the indices $x,y$, and $z$.

###  III.4 HLLD solver for relativistic magnetohydrodynamics

In the presence of an electromagnetic field, one choice for the Riemann solver
is the HLLD solver proposed in Ref. Mignone  _et al._ (2009). For this case,
we calculate the HLLD flux,
$\left(\tilde{f}^{(\hat{x})}_{A}\right)_{j+\frac{1}{2},k,l}$ in the tetrad
frame by solving the one-dimensional conservation law:

$\displaystyle\partial_{(\hat{t})}q_{A}+\partial_{(\hat{x})}f^{(\hat{x})}_{A}=0,$
(120)

where

$\displaystyle q_{A}=\left(\begin{array}[]{c}D\\\ J_{(\hat{\imath})}\\\
\rho_{\rm H}\\\ B^{(\hat{\imath})}\end{array}\right),$ (125) $\displaystyle
f^{(\hat{x})}_{A}=$ $\displaystyle\left(\begin{array}[]{c}Dv^{(\hat{x})}\\\
J_{(\hat{\imath})}v^{(\hat{x})}+P_{\text{tot}}{\delta^{(\hat{x})}}_{(\hat{\imath})}-\frac{B^{(\hat{x})}}{w^{2}}\left[B_{(\hat{\imath})}+(B^{(\hat{k})}u_{(\hat{k})})u_{(\hat{\imath})}\right]\\\
\rho_{\rm
H}v^{(\hat{x})}+P_{\text{tot}}v^{(\hat{x})}-\frac{1}{w}(B^{(\hat{k})}u_{(\hat{k})})B^{(\hat{x})}\\\
v^{(\hat{x})}B^{(\hat{\imath})}-v^{(\hat{\imath})}B^{(\hat{x})}\end{array}\right).$
(130)

Here, $q_{A}$ has seven components ($A=0,1,2,3,4,6,7$), and
$P_{\text{tot}}\equiv P+b^{2}/2$ is the total pressure (gas plus magnetic).
Note that the equation for $B^{(\hat{x})}$ is simply
$\partial_{(\hat{t})}B^{(\hat{x})}=0$, and thus, $B^{(\hat{x})}$ is constant
for the Riemann problem of the $x$-direction. Together with the initial
condition given by Eq. (98) for the relevant components, the full
magnetohydrodynamics Riemann fan consists of seven waves separating eight
states Font (2003). In the HLLD solver two of these seven waves (the slow
magnetosonic waves) are neglected. As a result, the Riemann fan with the HLLD
solver consists of five waves separating six states (see Fig. 2). In the HLLD
solver, we need to find the total pressure $P_{\text{tot}}$ which satisfies a
jump condition across the five waves. The numerical flux is then given by (see
the right-hand panel of Fig. 2)

$\displaystyle\left(\tilde{f}^{(\hat{x})}_{A}\right)_{j+\frac{1}{2}}$
$\displaystyle=\left\\{\begin{array}[]{lll}&(f^{(\hat{x})}_{A})_{L}&\text{if
}\lambda_{L}>v^{(\hat{x})}_{\text{interface}}\\\
&(f^{(\hat{x})}_{A})_{aL}&\text{if
}\lambda_{L}<v^{(\hat{x})}_{\text{interface}}<\lambda_{aL}\\\
&(f^{(\hat{x})}_{A})_{cL}&\text{if
}\lambda_{aL}<v^{(\hat{x})}_{\text{interface}}<\lambda_{c}\\\
&(f^{(\hat{x})}_{A})_{cR}&\text{if
}\lambda_{c}<v^{(\hat{x})}_{\text{interface}}<\lambda_{aR}\\\
&(f^{(\hat{x})}_{A})_{aR}&\text{if
}\lambda_{aR}<v^{(\hat{x})}_{\text{interface}}<\lambda_{R}\\\
&(f^{(\hat{x})}_{A})_{R}&\text{if
}\lambda_{R}<v^{(\hat{x})}_{\text{interface}},\\\ \end{array}\right.$ (137)

where

$\displaystyle\left(f^{(\hat{x})}_{A}\right)_{L/R}$
$\displaystyle=f^{(\hat{x})}_{A}\left(q_{L/R}\right),$ (138)
$\displaystyle\left(f^{(\hat{x})}_{A}\right)_{aL/aR}$
$\displaystyle=\left(f^{(\hat{x})}_{A}\right)_{L/R}$
$\displaystyle~{}+\lambda_{L/R}\left(\left(q_{A}\right)_{aL/aR}-\left(q_{A}\right)_{L/R}\right),$
(139) $\displaystyle\left(f^{(\hat{x})}_{A}\right)_{cL/cR}$
$\displaystyle=\left(f^{(\hat{x})}_{A}\right)_{aL/aR}$
$\displaystyle~{}+\lambda_{aL/aR}\left((q_{A})_{cL/cR}-(q_{A})_{aL/aR}\right).$
(140)

The latter two fluxes are obtained from the jump condition.

In the following subsections, we present specific quantities employed by the
HLLD solver: the characteristic speeds of the five waves, and the six states.

####  III.4.1 Characteristic speeds

For the fast waves, an approximate characteristic speed proposed in Refs.
Gammie  _et al._ (2003); Shibata and Sekiguchi (2005) is given by

$\displaystyle\lambda_{\text{FW}}^{\pm}(q_{A})=\frac{1}{1-v^{2}\zeta}\left[v^{(\hat{x})}\left(1-\zeta\right)\right.$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}\left.\pm\sqrt{\zeta}\sqrt{(1-v^{2})(1-v^{2}\zeta-(1-\zeta)(v^{(\hat{x})})^{2})}\right],$
(141)

where

$\displaystyle v^{2}=v^{(\hat{\imath})}v_{(\hat{\imath})},$ (142)
$\displaystyle\zeta=v_{A}^{2}+c_{s}^{2}-v_{A}^{2}c_{s}^{2},$ (143)
$\displaystyle v_{A}^{2}=\frac{b^{2}}{\rho h+b^{2}}.$ (144)

For the Alfvén wave, the characteristic speed is given by

$\displaystyle\lambda_{\text{Alf}}^{\pm}(q_{A})=\frac{b^{(\hat{x})}\pm
u^{(\hat{x})}\sqrt{\rho h+b^{2}}}{b^{(\hat{t})}\pm u^{(\hat{t})}\sqrt{\rho
h+b^{2}}},$ (145)

and for the contact wave by

$\displaystyle\lambda_{c}(q_{A})=v^{(\hat{x})}.$ (146)

####  III.4.2 L/R state

Given left- and right-state quantities, we first calculate the following
quantities which should be preserved when one crosses the fast waves:

$\displaystyle\left(R_{D}\right)_{L/R}=\left(\lambda
D-f^{(\hat{x})}_{0}\right)_{L/R},$ (147)
$\displaystyle\left(R_{J_{(\hat{\imath})}}\right)_{L/R}=\left(\lambda
J_{(\hat{\imath})}-f^{(\hat{x})}_{i}\right)_{L/R},$ (148)
$\displaystyle\left(R_{\rho_{\rm H}}\right)_{L/R}=\left(\lambda\rho_{\rm
H}-f^{(\hat{x})}_{4}\right)_{L/R},$ (149)
$\displaystyle\left(R_{B^{(\hat{k})}}\right)_{L/R}=\left(\lambda
B^{(\hat{k})}-f^{(\hat{x})}_{k}\right)_{L/R},$ (150)

where $\hat{\imath}=\hat{x},\hat{y},\hat{z}$ for $i=1,2,3$, respectively, in
Eq. (148). Also $\hat{k}=\hat{y},\hat{z}$ for $k=6,7$, respectively, in Eq.
(150). For the above quantities, we employ the characteristic speed defined by

$\displaystyle\lambda_{L}=\min\left(\lambda^{-}_{\text{FW}}(q_{L}),\lambda^{-}_{\text{FW}}(q_{R})\right),$
$\displaystyle\lambda_{R}=\max\left(\lambda^{+}_{\text{FW}}(q_{L}),\lambda^{+}_{\text{FW}}(q_{R})\right).$
(151)

####  III.4.3 aL/aR state

Given an initial guess for the unknown total pressure $P_{\text{tot}}$ (which
should be constant inside the Riemann fan), the three velocities in the $aL$
and $aR$ states are given by

$\displaystyle\left(v^{(\hat{x})}\right)_{aL/aR}$
$\displaystyle=\left(\frac{B^{(\hat{x})}(AB^{(\hat{x})}+\lambda
C)-(A+G)(P_{\text{tot}}+R_{J_{(\hat{x})}})}{X}\right)_{L/R},$ (152)
$\displaystyle\left(v^{(\hat{y})}\right)_{aL/aR}$
$\displaystyle=\left(\frac{QR_{J_{(\hat{y})}}+R_{B^{(\hat{y})}}[C+B^{(\hat{x})}(\lambda
R_{J_{(\hat{x})}}-R_{\rho_{\rm H}})]}{X}\right)_{L/R},$ (153)
$\displaystyle\left(v^{(\hat{z})}\right)_{aL/aR}$
$\displaystyle=\left(\frac{QR_{J_{(\hat{z})}}+R_{B^{(\hat{z})}}[C+B^{(\hat{x})}(\lambda
R_{J_{(\hat{x})}}-R_{\rho_{\rm H}})]}{X}\right)_{L/R},$ (154)

where

$\displaystyle A=R_{J_{(\hat{x})}}-\lambda R_{\rho_{\rm
H}}+P_{\text{tot}}(1-\lambda^{2}),$ (155) $\displaystyle
G=R_{B^{(\hat{y})}}R_{B^{(\hat{y})}}+R_{B^{(\hat{z})}}R_{B^{(\hat{z})}},$
(156) $\displaystyle
C=R_{J_{(\hat{y})}}R_{B^{(\hat{y})}}+R_{J_{(\hat{z})}}R_{B^{(\hat{z})}},$
(157) $\displaystyle Q=-A-G+(B^{(\hat{x})})^{2}(1-\lambda^{2}),$ (158)
$\displaystyle X=B^{(\hat{x})}(A\lambda B^{(\hat{x})}+C)-(A+G)(\lambda
P_{\text{tot}}+R_{\rho_{\rm H}}).$ (159)

Note that $aL$ and $aR$ on the left-hand side of Eqs. (152)–(154) correspond
to $L$ and $R$ for $R_{J_{(\hat{\imath})}}$, $R_{\rho_{\rm H}}$,
$R_{B^{(\hat{k})}}$, and $\lambda$ on the right-hand side of the same
equations, respectively. With these velocity components, the magnetic field is
calculated from the jump condition by

$\displaystyle\left(B^{(\hat{k})}\right)_{aL/aR}=\frac{\left(R_{B^{(\hat{k})}}\right)_{L/R}-B^{(\hat{x})}\left(v^{(\hat{k})}\right)_{aL/aR}}{\lambda_{L/R}-\left(v^{(\hat{x})}\right)_{aL/aR}}$
(160)

for $\hat{k}=\hat{y},\hat{z}$. The total enthalpy density is calculated by

$\displaystyle(\rho h_{\text{tot}})_{aL/aR}\equiv(\rho h+b^{2})_{aL/aR}$
$\displaystyle=P_{\text{tot}}+\frac{\left(R_{\rho_{\rm
H}}\right)_{L/R}-\left(v^{(\hat{\imath})}\right)_{aL/aR}\left(R_{J_{(\hat{\imath})}}\right)_{L/R}}{\lambda_{L/R}-\left(v^{(\hat{x})}\right)_{aL/aR}}.$
(161)

The conserved quantities necessary for the numerical flux in Eq. (139) and in
the Eulerian frame (see Eqs. (78)–(85)) are calculated by

$\displaystyle
D_{aL/aR}=\frac{\left(R_{D}\right)_{L/R}}{\lambda_{L/R}-\left(v^{(\hat{x})}\right)_{aL/aR}},$
(162) $\displaystyle\left(\rho_{\rm H}\right)_{aL/aR}=$
$\displaystyle\frac{\left(R_{\rho_{\rm
H}}\right)_{L/R}+P_{\text{tot}}\left(v^{(\hat{x})}\right)_{aL/aR}-\left(v^{(\hat{k})}B^{(\hat{k})}\right)_{aL/aR}B^{(\hat{x})}}{\lambda_{L/R}-\left(v^{(\hat{x})}\right)_{aL/aR}},$
(163) $\displaystyle\left(J_{(\hat{\imath})}\right)_{aL/aR}=\left((\rho_{\rm
H}+P_{\text{tot}})v^{(\hat{\imath})}-(v^{(\hat{k})}B^{(\hat{k})})B^{(\hat{\imath})}\right)_{aL/aR}.$
(164)

####  III.4.4 cL/cR state

Following Ref. Mignone  _et al._ (2009), we first define

$\displaystyle\sigma^{(\hat{\mu})}=\eta u^{(\hat{\mu})}+b^{(\hat{\mu})},$
(165) $\displaystyle\eta=\pm\text{sgn}(B^{(\hat{x})})\sqrt{\rho
h_{\text{tot}}},$ (166)

where the plus (minus) sign corresponds to the right (left) state. We then
define $K^{(\hat{k})}$ by

$\displaystyle
K^{(\hat{k})}\equiv\frac{\sigma^{(\hat{k})}}{\sigma^{(\hat{0})}}=v^{(\hat{k})}+\frac{B^{(\hat{k})}}{w\sigma^{(\hat{0})}}.$
(167)

Here $K^{(\hat{x})}$ is nothing other than the Alfvén wave speed in the
$x$-direction. From the jump condition one can find that $K^{(\hat{\imath})}$,
$\rho h_{\text{tot}}$, $D/w\sigma^{(\hat{0})}$, and $\eta$ do not change
across the Alfvén waves. Therefore, $\eta$, $K^{(\hat{\imath})}$, and the
total enthalpy density are calculated by

$\displaystyle\eta_{cL/cR}=\eta_{aL/aR},$ (168)
$\displaystyle(K^{(\hat{\imath})})_{cL/cR}=(K^{(\hat{\imath})})_{aL/aR}$
$\displaystyle=\frac{\left(R_{J_{(\hat{\imath})}}\right)_{L/R}+P_{\text{tot}}\delta_{(\hat{\imath})(\hat{x})}+\left(R_{B^{(\hat{i)}}}\right)_{L/R}\eta_{aL/aR}}{\lambda_{L/R}P_{\text{tot}}+\left(R_{\rho_{\rm
H}}\right)_{L/R}+B^{(\hat{x})}\eta_{aL/aR}},$ (169) $\displaystyle(\rho
h_{\text{tot}})_{cL/cR}=(\rho h_{\text{tot}})_{aL/aR},$ (170)

where $cL$ and $cR$ on the left-hand side of the equations correspond to $aL$
and $aR$ on the right-hand side of the same equations, respectively.

The magnetic field and the three velocity in the $cL$ and $cR$ states are
calculated by

$\displaystyle(B^{(\hat{k})})_{cL}=(B^{(\hat{k})})_{cR}$
$\displaystyle=\Big{[}\left\\{B^{(\hat{k})}(\lambda-v^{(\hat{x})})+B^{(\hat{x})}v^{(\hat{k})}\right\\}_{aR}$
$\displaystyle-\left\\{B^{(\hat{k})}(\lambda-v^{(\hat{x})})+B^{(\hat{x})}v^{(\hat{k})}\right\\}_{aL}\Big{]}\frac{1}{\lambda_{aR}-\lambda_{aL}},$
(171)
$\displaystyle(v^{(\hat{\imath})})_{cL/cR}=\left(K^{(\hat{\imath})}-\frac{B^{(\hat{\imath})}(1-K^{(\hat{k})}K^{(\hat{k})})}{\eta-K^{(\hat{l})}B^{(\hat{l})}}\right)_{cL/cR},$
(172)

and the characteristic speed is

$\displaystyle\lambda_{aL/aR}=K^{(\hat{x})}_{aL/aR}.$ (173)

We impose the continuity condition on the normal velocity across the contact
discontinuity, i.e., $v^{(\hat{x})}_{cL}=v^{(\hat{x})}_{cR}$, by

$\displaystyle\Delta
K^{(\hat{x})}\left[1-B^{(\hat{x})}\left(Y_{R}-Y_{L}\right)\right]=0,$ (174)
$\displaystyle
Y_{L/R}=\left(\frac{1-K^{(\hat{\imath})}K^{(\hat{\imath})}}{\eta\Delta
K^{(\hat{x})}-\Delta K^{(\hat{x})}K^{(\hat{j})}B^{(\hat{j})}}\right)_{cL/cR},$
(175)

where $\Delta K^{(\hat{x})}=K^{(\hat{x})}_{aR}-K^{(\hat{x})}_{aL}$. This
equation gives an improved guess of the total pressure in the next iteration
step. Then we go back to Eq. (152) and repeat the same procedure until it
converges with sufficient accuracy. In practice, we employ the Newton-Raphson
method to solve Eq. (174).

The conserved quantities necessary for the numerical flux in Eq. (140) and in
the Eulerian frame (see Eqs. (78)–(85)) are

$\displaystyle
D_{cL/cR}=D_{aL/aR}\frac{\lambda_{aL/aR}-v^{(\hat{x})}_{aL/aR}}{\lambda_{aL/aR}-v^{(\hat{x})}_{cL/cR}},$
(176) $\displaystyle\left(\rho_{\rm
H}\right)_{cL/cR}=\frac{1}{\lambda_{aL/aR}-v^{(\hat{x})}_{cL/cR}}$
$\displaystyle\times\Big{[}\lambda_{aL/aR}\left(\rho_{{\rm
H}0}\right)_{aL/aR}-\left(J_{(\hat{x})}\right)_{aL/aR}+P_{\text{tot}}v^{(\hat{x})}_{cL/cR}$
$\displaystyle~{}~{}~{}~{}~{}-(v^{(\hat{\imath})}B^{(\hat{\imath})})_{cL/cR}B^{(\hat{x})}\Big{]},$
(177)
$\displaystyle\left(J_{(\hat{\imath})}\right)_{cL/cR}=\left(\left(\rho_{{\rm
H}}\right)_{cL/cR}+P_{\text{tot}}\right)v^{(\hat{\imath})}_{cL/cR}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-(v^{(\hat{k})}B^{(\hat{k})})_{cL/cR}B^{(\hat{\imath})}_{cL/cR}.$
(178)

The equivalent expressions in the $y$\- and $z$-directions are given by
permutation of the indices $x,y$, and $z$.

###  III.5 Electric-field evaluation

The constrained transport method used to enforce the divergence-free condition
on the magnetic field requires us to evaluate the electric field defined at
the cell edges. Gardiner and Stone Gardiner and Stone (2008) proposed a method
for evaluating the electric-field by utilizing the numerical fluxes which are
obtained by the Riemann solver. In their method, for example, the
$z$-component of the electric field is evaluated by

$\displaystyle\tilde{E}^{z}_{j+\frac{1}{2},k+\frac{1}{2},l}=\frac{1}{4}\Big{(}\tilde{E}^{z}_{j+\frac{1}{2},k,l}+\tilde{E}^{z}_{j+\frac{1}{2},k+1,l}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\tilde{E}^{z}_{j,k+\frac{1}{2},l}+\tilde{E}^{z}_{j+1,k+\frac{1}{2},l}\Big{)}$
$\displaystyle+\frac{\Delta y}{8}\left(\left(\frac{\partial E^{z}}{\partial
y}\right)_{j+\frac{1}{2},k+\frac{1}{4},l}-\left(\frac{\partial E^{z}}{\partial
y}\right)_{j+\frac{1}{2},k+\frac{3}{4},l}\right)$ $\displaystyle+\frac{\Delta
x}{8}\left(\left(\frac{\partial E^{z}}{\partial
x}\right)_{j+\frac{1}{4},k+\frac{1}{2},l}-\left(\frac{\partial E^{z}}{\partial
x}\right)_{j+\frac{3}{4},k+\frac{1}{2},l}\right)$ (179)

where

$\displaystyle\left(\frac{\partial E^{z}}{\partial
y}\right)_{j+\frac{1}{2},k+\frac{1}{4},l}$
$\displaystyle=\left\\{\begin{array}[]{l}\frac{2\left(\tilde{E}^{z}_{j,k+\frac{1}{2},l}-E^{z}_{j,k,l}\right)}{\Delta
y}\text{ for }\tilde{v}^{x}_{j+\frac{1}{2},k,l}>0,\\\
\frac{2\left(\tilde{E}^{z}_{j+1,k+\frac{1}{2},l}-E^{z}_{j+1,k,l}\right)}{\Delta
y}\text{ for }\tilde{v}^{x}_{j+\frac{1}{2},k,l}<0,\\\
\frac{\left(\tilde{E}^{z}_{j,k+\frac{1}{2},l}-E^{z}_{j,k,l}+\tilde{E}^{z}_{j+1,k+\frac{1}{2},l}-E^{z}_{j+1,k,l}\right)}{\Delta
y}\text{ otherwise,}\\\ \end{array}\right.$ $\displaystyle\left(\frac{\partial
E^{z}}{\partial y}\right)_{j+\frac{1}{2},k+\frac{3}{4},l}$
$\displaystyle=\left\\{\begin{array}[]{l}\frac{2\left(E^{z}_{j,k+1,l}-\tilde{E}^{z}_{j,k+\frac{1}{2},l}\right)}{\Delta
y}\text{ for }\tilde{v}^{x}_{j+\frac{1}{2},k+1,l}>0,\\\
\frac{2\left(E^{z}_{j+1,k+1,l}-\tilde{E}^{z}_{j+1,k+\frac{1}{2},l}\right)}{\Delta
y}\text{ for }\tilde{v}^{x}_{j+\frac{1}{2},k+1,l}<0,\\\
\frac{\left(E^{z}_{j,k+1,l}-\tilde{E}^{z}_{j,k+\frac{1}{2},l}+E^{z}_{j+1,k+1,l}-\tilde{E}^{z}_{j+1,k+\frac{1}{2},l}\right)}{\Delta
y}\text{ otherwise,}\\\ \end{array}\right.$ $\displaystyle\left(\frac{\partial
E^{z}}{\partial x}\right)_{j+\frac{1}{4},k+\frac{1}{2},l}$
$\displaystyle=\left\\{\begin{array}[]{l}\frac{2\left(\tilde{E}^{z}_{j+\frac{1}{2},k,l}-E^{z}_{j,k,l}\right)}{\Delta
x}\text{ for }\tilde{v}^{y}_{j,k+\frac{1}{2},l}>0,\\\
\frac{2\left(\tilde{E}^{z}_{j+\frac{1}{2},k+1,l}-E^{z}_{j,k+1,l}\right)}{\Delta
x}\text{ for }\tilde{v}^{y}_{j,k+\frac{1}{2},l}<0,\\\
\frac{\left(\tilde{E}^{z}_{j+\frac{1}{2},k,l}-E^{z}_{j,k,l}+\tilde{E}^{z}_{j+\frac{1}{2},k+1,l}-E^{z}_{j,k+1,l}\right)}{\Delta
x}\text{ otherwise,}\\\ \end{array}\right.$ $\displaystyle\left(\frac{\partial
E^{z}}{\partial x}\right)_{j+\frac{3}{4},k+\frac{1}{2},l}$
$\displaystyle=\left\\{\begin{array}[]{l}\frac{2\left(E^{z}_{j+1,k,l}-\tilde{E}^{z}_{j+\frac{1}{2},k,l}\right)}{\Delta
x}\text{ for }\tilde{v}^{y}_{j+1,k+\frac{1}{2},l}>0,\\\
\frac{2\left(E^{z}_{j+1,k+1,l}-\tilde{E}^{z}_{j+\frac{1}{2},k+1,l}\right)}{\Delta
x}\text{ for }\tilde{v}^{y}_{j+1,k+\frac{1}{2},l}<0,\\\
\frac{\left(E^{z}_{j+1,k,l}-\tilde{E}^{z}_{j+\frac{1}{2},k,l}+E^{z}_{j+1,k+1,l}-\tilde{E}^{z}_{j+\frac{1}{2},k+1,l}\right)}{\Delta
x}\text{ otherwise.}\\\ \end{array}\right..$

Here $\tilde{v}^{x}_{j+\frac{1}{2},k,l}$ and
$\tilde{E}^{z}_{j+\frac{1}{2},k,l}$ are identical to the fluxes
$\left(\tilde{F}^{x}_{0}\right)_{j+\frac{1}{2},k,l}$ and
$\left(-\tilde{F}^{x}_{6}\right)_{j+\frac{1}{2},k,l}$ in Eqs. (78) and
(III.2), which are given by the Riemann solver in the $x$-direction.
Similarly, $\tilde{v}^{y}_{j,k+\frac{1}{2},l}$ and
$\tilde{E}^{z}_{j,k+\frac{1}{2},l}$ are given by the Riemann solver in the
$y$-direction. $E^{z}_{j,k,l}$ is calculated from the quantities defined at
the cell center, i.e., Eqs. (57)–(59) and the three velocity. Therefore, the
accuracy of this constrained transport scheme depends on the accuracy of an
employed Riemann solver. Equivalent expressions for the $x$\- and
$y$-components of the electric field are given by permutation of the indices
$x,y,$ and $z$. These electric fields are used to update the magnetic field in
Eqs. (54)–(56).

Table 1: Initial conditions used for special relativistic one-dimensional test problems. The third column shows the $\Gamma$ index and the second-to-last column shows the final time of the simulations, $t$. Test problem | State | $\Gamma$ | $\rho$ | $v^{x}$ | $v^{y}$ | $v^{z}$ |  $P$ |  $B^{x}$ |  $B^{y}$ | $B^{z}$ |  $t$ |  CFL
---|---|---|---|---|---|---|---|---|---|---|---|---
Problem HD1 | L | 4/3 | 1 | 0.9 | 0 | 0 | 1 | – | – | – | 0.4 | 0.8
| R |  | 1 | 0 | 0 | 0 | 10 | – | – | – |  |
Problem HD2 | L | 5/3 | 1 | $-0.6$ | 0 | 0 | 10 | – | – | – | 0.4 | 0.8
| R |  | 10 | 0.5 | 0 | 0 | 20 | – | – | – |  |
Problem HD3 | L | 5/3 | 10 | 0 | 0 | 0 | 40 | – | – | – | 0.4 | 0.8
| R |  | 1 | 0 | 0 | 0 | 3 | – | – | – |  |
Problem HD4 | L | 5/3 | 1 | 0 | 0 | 0 | $10^{3}$ | – | – | – | 0.4 | 0.8
| T |  | 1 | 0 | 0 | 0 | $10^{-2}$ | – | – | – | 0.4 | 0.8
Problem MHD1 | L | 5/3 | 10 | 0 | 0.7 | 0.2 | 1 | 5 | 1 | 0.5 | 1 | 0.8
| R |  | 1 | 0 | 0.7 | 0.2 | 1 | 5 | 1 | 0.5 |  |
Problem MHD2 | L | 5/3 | 1 | 0.4 | $-0.3$ | 0.5 | 1 | 2.4 | 1 | $-1.6$ | 1 | 0.8
| R |  | 1 | 0.377347 | $-0.482389$ | 0.424190 | 1 | 2.4 | $-0.1$ | $-2.1728213$ |  |
Problem MHD3 | L | 2 | 1 | 0 | 0 | 0 | 1 | 0.5 | 1 | 0 | 0.4 | 0.8
| R |  | 0.125 | 0 | 0 | 0 | 0.1 | 0.5 | $-1$ | 0 |  |
Problem MHD4 | L | 5/3 | 1.08 | 0.4 | 0.3 | 0.2 | 0.95 | 2 | 0.3 | 0.3 | 0.55 | 0.8
| R |  | 1 | $-0.45$ | $-0.2$ | 0.2 | 1 | 2 | $-0.7$ | $-0.5$ |  |
Problem MHD5 | L | 5/3 | 1 | 0.999 | 0 | 0 | 0.1 | 10 | 7 | 7 | 0.4 | 0.8
| R |  | 1 | $-0.999$ | 0 | 0 | 0.1 | 10 | $-7$ | $-7$ |  |
Problem MHD6 | L | 5/3 | 1 | 0 | 0.3 | 0.4 | 5 | 1 | 6 | 2 | 0.5 | 0.8
| R |  | 0.9 | 0 | 0 | 0 | 5.3 | 1 | 5 | 2 |  |

In the rest of this paper, we refer to this particular algorithm for
evaluating the electric field as CT_GS. On the other hand, the electric-field
evaluation algorithm which was originally implemented in our code, and which
is based on HLLE Del Zanna  _et al._ (2003); Shibata and Sekiguchi (2005), is
referred to as CT_HLLE Shibata and Sekiguchi (2005). For the base Riemann
solver, we use either HLLC, HLLD, or HLLE. Here the last one is the base
Riemann solver which was originally implemented in our code Shibata and
Sekiguchi (2005). In the hydrodynamics test simulations shown in the next
section, we refer to the particular combination of numerical schemes used in a
particular test problem in terms of the base solver, only. In the
magnetohydrodynamics test simulations, we describe a simulation both in terms
of the base Riemann solver and in terms of the algorithm used for the
evaluation of the electric-field. For example, HLLD-CT_GS means that the
(base) Riemann solver is HLLD and the electric-field evaluation is CT_GS.

##  IV Validation of the HLLC and HLLD solvers

In this section, we introduce various problems designed to test the
implementation of the advanced Riemann solvers and constrained transport
algorithm discussed in the previous section. We start with a common suite of
one-dimensional special relativistic shock-tube problems in both hydrodynamics
and magnetohydrodynamics (see Sec. IV.1). Next, in Sec. IV.2, we turn our
attention to multi-dimensional hydrodynamics and magnetohydrodynamics test
problems in special relativity (specifically, we consider a two-dimensional
hydrodynamical shock, a cylindrical hydrodynamical blast wave, a
magnetohydrodynamical current sheet, and the Kelvin-Helmholtz instability in
magnetohydrodynamics). In Sec. IV.3 we consider Bondi flow onto a black hole
(in both hydrodynamics and magnetohydrodynamics) as a test problem in general
relativity with a static spacetime.

For all the test problems we assume a $\Gamma$-law equation of state given by

$\displaystyle P=\left(\Gamma-1\right)\rho\varepsilon.$ (180)

We also employ a cell-centered grid structure in which the $x$-coordinate 222
In Ref. Kiuchi  _et al._ (2012), we employed a vertex-centered grid structure.
We updated the interpolation scheme of the metric and fluid at the refinement
boundary for cell-centered grid structure for a simulation in a dynamical
spacetime. is given by

$\displaystyle x_{j}=\left(j+\frac{1}{2}\right)\Delta x,$ (181)

with $j\in[-N_{x}-1,N_{x}]$ and grid spacing $\Delta x$ (and likewise for the
$y$\- and $z$-components). As a time integrator, we employ the fourth-order
Runge-Kutta method (RK4) in all our test simulations. For reconstruction of
the solution at cell-interfaces, we employ either 1st-order reconstruction or
3rd-order piecewise parabolic method (PPM) (Shibata and Sekiguchi, 2005;
Colella and Woodward, 1984). For the PPM reconstruction, we employ the min-mod
limiter function with a compression parameter which is generally set to $b=2$
Shibata and Sekiguchi (2005), though in some cases we employ different values
of $b$.

Figure 3: Left: Rest-mass density profile for Problem HD1 at $t=0.4$ (contact
wave located between left- and right-propagating shock waves). The blue and
green curves show the results with the HLLC and HLLE solvers, respectively.
The dashed and solid curves show the results with 1st-order reconstruction and
3rd-order (PPM) reconstruction, respectively. Right: Rest-mass density profile
for Problem HD2 at $t=0.4$ (contact wave located between left- and right-
propagating rarefaction waves). We employ RK4 for the time integration in all
the simulations. The blue and green solid curves are indistinguishable on the
scale of this plot.

###  IV.1 Special relativistic one-dimensional problems

First, we consider special relativistic problems in one spatial dimension.
With this setup, the tetrad basis in Sec. III.2 is reduced to a coordinate
vector in Minkowski spacetime. Thus, the setup is suitable for validating the
Riemann solvers described in Sec. III.3 and III.4. We assume Minkowski metric,
and thus turn off the solver for Einstein’s equations in the code. The initial
conditions for all the one-dimensional test problems are summarized in Table
1. We note that the test suite employed in this paper is the same as that
presented in Refs. Mignone  _et al._ (2009); Mattia and Mignone (2021).

Figure 4: Same as Fig. 3, but for Problem HD3 (left) and Problem HD4 (right).

####  IV.1.1 Hydrodynamics: one-dimensional shock tubes

The first special relativistic hydrodynamics test (HD1) is the computation of
a contact discontinuity. For this we prepare a simulation domain of
$x\in[-1,1]$ with $\Delta x=0.01$ and $N_{x}=100$. We integrate the numerical
solution up to $t=0.4$. In the left panel of Fig. 3 we plot the rest-mass
density profile at the end of the simulation. In this problem, left- and
right-propagating shock waves appear from the initial discontinuity, with a
contact discontinuity sandwiched between them. The blue and green curves
denote the numerical solution with the HLLC and HLLE solvers, respectively.
The solid and dashed curves denote the simulation results with 3rd-order PPM
reconstruction and 1st-order reconstruction, respectively. First, we consider
the results obtained with 1st-order reconstruction (dashed curves). With the
HLLC solver, the contact discontinuity located at $x\approx 0.2$ is more
sharply captured than with the HLLE solver. This behavior is expected since
the HLLC solver explicitly restores the contact wave inside the Riemann fan.
When we employ 3rd-order reconstruction, however, we find that there is no
qualitative difference between the two solvers. This suggests that the weak
point of a particular solver may be alleviated by using a high enough
resolution.

Figure 5: Left: Rest-mass density profile in Problem MHD1 (a problem with a
stationary contact discontinuity) at $t=1$. The blue, green, and cyan curves
present the result with the HLLD-CT_GS, HLLD-CT_HLLE, and HLLE-CT_HLLE
solvers, respectively. The solid and dashed curves show the results with 3rd-
order PPM cell reconstruction and 1st-order cell reconstruction, respectively.
Right: Same as the left panel, but for the $B^{y}$ profile in Problem MHD2 (a
problem with a stationary rotational discontinuity). The insets are a close-up
of the discontinuity with the HLLD-CT_GS solver and 1st-order cell
reconstruction.

The right panel of Fig. 3 shows the rest-mass density profile for the second
test problem (Problem HD2) listed in Table 1. The simulation domain and the
grid spacing are the same as those in Problem HD1. In this problem, left- and
right-going rarefaction waves propagate away from the initial discontinuity,
and a contact discontinuity appears between the two and is located at
$x\approx-0.1$. As in our first test problem, we find that the contact
discontinuity is more sharply captured with the HLLC solver than that with the
HLLE solver when 1st-order reconstruction is used, while we find no
qualitative difference between the numerical solutions obtained with the two
solvers when we employ 3rd-order-accurate reconstruction.

The third hydrodynamics test problem (Problem HD3 in Table 1) is the often-
employed shock-tube problem. Here, the simulation domain spans
$x\in[-0.5,0.5]$ with $\Delta x=0.005$ and $N_{x}=100$. In this problem, the
initial discontinuity decays into a left-propagating rarefaction wave and a
right-propagating shock wave. The contact discontinuity adjacent to the shock
wave also propagates to the right. The left panel of Fig. 4 shows the rest-
mass density profile at the end of the simulation for which the contact
discontinuity is located at $x\approx 0.25$. In this problem, we find that
there is no qualitative difference between the numerical solutions with the
two solvers irrespective of the cell reconstruction accuracy. This behaviour
is also reported in Ref. Mignone and Bodo (2005). For obtaining an accurate
result for this particular shock-tube problem it is necessary to employ an
accurate reconstruction method. This suggests that employing an accurate
reconstruction method is as important as employing an accurate solver in
numerical hydrodynamics at least in the one-dimensional problems.

For the fourth (final) hydrodynamics test problem (Problem HD4), we employ a
simulation domain of $x\in[-0.5,0.5]$ with a grid spacing of $\Delta
x=0.0025$, i.e., $N_{x}=200$. The solution consists of a left-propagating
rarefaction wave and a right-propagating shock wave. Note that the result
differs from that in Problem HD3 as the shock is much stronger compared to the
one in Problem HD3 because of the initial large pressure jump (see Table I). A
right-propagating contact discontinuity appears adjacent to the shock wave. We
plot the rest-mass density profile at $t=0.4$ in the right panel of Fig. 4. We
find that the contact discontinuity (located at $x\approx 0.35$) is more
sharply resolved with the HLLC solver than with the HLLE solver when we employ
3rd-order reconstruction. We find that the compression parameter $b$ for the
min-mod function in the PPM cell reconstruction should be reduced to be 1 in
this problem (i.e., a steep limiter function does not work; see, e.g., Ref.
Shibata and Sekiguchi (2005)). Otherwise, spurious waves appear irrespective
of which solver is used (not shown).

Figure 6: Profile of the rest mass density (top-left), the $x$-component of
the three velocity (top-right), the $y$-component of the three velocity
(bottom-left), and the $y$-component of the magnetic field (bottom-right) at
$t=0.4$ in Problem MHD3. The blue, green, and cyan curves denote the numerical
solution with the HLLD-CT_GS, HLLD-CT_HLLE, and HLLE-CT_HLLE solvers,
respectively. We employ RK4 with 3rd-order PPM cell reconstruction (solid
curves), and also with 1st-order cell reconstruction (dashed curves).

Figure 7: Profile of the rest mass density (top-left), the $x$-component of
the three velocity (top-right), the $y$-component of the three velocity
(middle-left), the $z$-component of the three velocity (middle-right), the
$y$-component of the magnetic field (bottom-left), and the $z$-component of
the magnetic field (bottom-right) at $t=0.55$ in Problem MHD4. The blue,
green, and cyan curves denote the numerical solution with the HLLD-CT_GS,
HLLD-CT_HLLE, and HLLE-CT_HLLE solvers, respectively. We employ RK4 with 3rd-
order PPM cell reconstruction (solid curves), and also with 1st-order
reconstruction (dashed curves). The insets show a close-up of the
discontinuity.

Figure 8: Same as Fig. 6, but for Problem MHD5.

Figure 9: Same as Fig. 7, but for Problem MHD6.

####  IV.1.2 Magnetohydrodynamics: one-dimensional shock tubes

In this section we consider six special relativistic magnetohydrodynamics test
problems in one spatial dimension. All the test problems except for Problem
MHD6 are carried out in a domain of size $x\in[-0.5,0.5]$ with grid spacing
$\Delta x=0.005$ (i.e., $N_{x}=100$). For Problem MHD6, the domain is
identical but we employ a higher resolution with $\Delta x=0.0025$ and
$N_{x}=200$.

In the first problem (Problem MHD1 in Table 1) the solution consists of a
stationary contact discontinuity. The left panel of Fig. 5 plots the rest-mass
density profile at the end of the simulation. Because the HLLD-CT_GS and HLLD-
CT_HLLE solvers exactly capture the contact discontinuity, the numerical
solutions remain stationary even when we employ 1st-order reconstruction (see
the inset in the left panel of Fig. 5). On the other hand, with the HLLE-
CT_HLLE solver, the initial contact discontinuity is broadened because this
solver neglects the contact discontinuity inside the Riemann fan. However when
we employ 3rd-order (PPM) reconstruction, this spurious broadening of the
contact discontinuity is suppressed, although the contact wave is still not
resolved as sharply as it is with the HLLD solver.

In the second magnetohydrodynamics test problem (Problem MHD2 in Table 1) we
model a stationary rotational discontinuity (i.e., an Alfvén wave). The right
panel of Fig. 5 presents the profile of the $y$-component of the magnetic
field at the end of the simulation. When we employ 1st-order cell
reconstruction, the HLLD-CT_GS solver reproduces the stationary solution (see
the inset in the right panel of Fig. 5). This is because the rotational
discontinuity is captured exactly by the HLLD solver, and the electric field
at the cell edge is evaluated with the numerical flux, i.e., the electric
field at the cell interface, given by the HLLD solver with the CT_GS scheme
(see, e.g., Eq. (179)). With the HLLD-CT_HLLE solver, on the other hand, the
rotational discontinuity is broadened because the electric field at the cell
interface given by the HLLD solver is not used to evaluate the electric field
at the cell edge in the CT_HLLE scheme. With the HLLE-CT_HLLE solver, the
rotational discontinuity inside the Riemann fan is not captured. As a result,
the initial rotational discontinuity is spuriously broadened. This drawback is
improved by employing 3rd-order PPM cell reconstruction in the HLLD-CT_HLLE
and HLLE-CT_HLLE solvers. Note that for this problem, we employ the
compression parameter $b=3$ in the min-mod function for the PPM reconstruction
in the HLLD-CT_GS run. Otherwise, we find the over- and under-shoot in the
vicinity of the initial rotational discontinuity (not shown) because the
default value of $b=2$ is not sufficient to capture the initial steep profile.

The third magnetohydrodynamics test problem (Problem MHD3) is the relativistic
extension of the Brio-Wu shock tube Brio and Wu (1988). In this problem, the
solution consists of a left-propagating rarefaction wave, a right-propagating
slow shock wave (located at $x\approx 0.18$ in Fig. 6), and a right-
propagating rarefaction wave. In addition there is a right-ward propagating
contact discontinuity located at $x\approx 0.15$ in Fig. 6 adjacent to the
shock wave. Finally at $x\approx 0$, a compound wave appears. When we use 1st-
order reconstruction, the contact discontinuity is captured more sharply with
the HLLD-CT_GS solver than with the HLLD-CT_HLLE solver, while the HLLE-
CT_HLLE solver cannot capture the contact discontinuity at all if the 1st-
order reconstruction is used. The slow shock is also captured more sharply
with the HLLD-CT_GS solver than with the HLLD-CT_HLLE solver, while with the
HLLE-CT_HLLE solver the slow shock wave is significantly broadened. This
feature is also found for the compound wave. While the various waves are
better captured in 3rd-order PPM reconstruction irrespective of the chosen
solvers, we find that the higher order reconstruction method induces
artificial oscillatory behavior behind the compound wave (in the regions of
$-0.2\lesssim x\lesssim 0$) in both the rest-mass density and in the
$x$-component of the velocity (see also in the regions of $0.2\lesssim
x\lesssim 0.4$ in the $x$-component of the velocity). The amplitude of these
oscillations is reduced when we employ the diffusive compression parameter of
the PPM reconstruction $b=1$.

The fourth magnetohydrodynamics test problem (Problem MHD4) consists of a left
(right)-propagating fast wave located at $x\approx-0.4\,(+0.4)$, a left-
propagating rarefaction wave ($x\approx-0.3$), a contact discontinuity
($x\approx-0.04$), a right-propagating slow wave ($x\approx+0.2$), and,
finally, a left (right)-propagating Alfvén wave (located at
$x\approx-0.33\,(+0.22)$). See Fig. 7 for the solutions (the inset in the
$B^{y}$ ($B^{z}$) panel shows a close-up region of the right (left) Alfvén
waves). When we employ 1st-order reconstruction, both solvers are able to
capture the fast waves, but the contact discontinuity is captured more sharply
with the HLLD-CT_GS or HLLD-CT_HLLE solver than with the HLLE-CT_HLLE solver
(see the inset in the rest-mass density profile in Fig. 7 which shows a close-
up of the contact discontinuity). Irrespective of the solvers, it is hard to
distinguish the slow wave and the right-propagating Alfvén wave, and also
between the the rarefaction wave and the left-propagating Alfvén wave. When we
employ 3rd-order PPM reconstruction, on the other hand, we find no qualitative
difference in the numerical solutions between the different solvers.

For our fifth magnetohydrodynamics test problem (Problem MHD5) we consider the
relativistic collision of two streams. Figure 8 shows the result at $t=0.4$.
In this problem the solution consists of left (right)-propagating fast waves
located at $x\approx-0.3\,(+0.3)$, and left (right)-propagating slow waves
located at $x\approx-0.06\,(+0.06)$. When we employ 1st-order cell
reconstruction, the slow waves are captured more sharply with the HLLD
solvers, i.e., the HLLD-CT_GS or HLLD-CT_HLLE solvers than with the HLLE-
CT_HLLE solver. On the other hand, the resolution across the outermost fast
waves is essentially the same for all solvers. Irrespective of the solver or
reconstruction method used, a spurious undershoot in the rest-mass density
appears at $x\approx 0$. This is known as the wall-heating problem Noh (1987):
it is well-known that Godunov-type schemes cannot avoid this pathological
behavior. As reported in Ref. Mignone  _et al._ (2009), the undershoot is
shallower with the HLLE-CT_HLLE solver due to the solver’s larger numerical
diffusion. When we employ 3rd-order PPM reconstruction, both the HLLD and HLLE
solvers are equally capable of capturing the slow waves as well as the fast
waves.

In the final problem (Problem MHD6) in our one-dimensional suite, the solution
consists of all seven waves Giacomazzo and Rezzolla (2006). The numerical
results are shown in Fig. 9. In this problem, a contact discontinuity appears
at $x\approx 0.05$, a rarefaction wave propagates to the left of the contact
discontinuity, which can be seen at $x\approx-0.4$, and the rotational
discontinuity at $x\approx-0.06$ and the slow shock at $x\approx-0.04$ follow
the rarefaction wave (see the inset in the panel for $B^{y}$ in Fig. 9). To
the right of the contact discontinuity, a fast shock propagates up to
$x\approx 0.4$. The rotational discontinuity at $x\approx 0.08$ and the slow
shock at $x\approx 0.06$ follow the fast shock (again, this is most easily
seen in the inset in the panel for $B^{y}$ in Fig. 9).

When we use 1st-order cell reconstruction (dashed curves), the contact
discontinuity is resolved only with the HLLD solvers, i.e. HLLD-CT_GS or HLLD-
CT_HLLE, (see the panel for $\rho$ in Fig. 9). With 1st-order reconstruction,
however, it is difficult to disentangle the left/right-propagating rotational
discontinuities and left/right-propagating slow shocks, even with the HLLD
solver (see the dashed curves in the inset in the panel for $B^{y}$ in Fig.
9). When we employ 3rd-order PPM reconstruction (solid curves), on the other
hand, the difference between the various Riemann solvers is striking. With the
HLLD-CT_GS solver, the left/right-propagating rotational discontinuities and
slow shocks are captured as plotted in the inset in the panel for $B^{y}$ in
Fig. 9. With the HLLD-CT_HLLE or HLLE-CT_HLLE solvers, the right-propagating
rotational discontinuity and the left/right-propagating slow shock are
captured, but the left-propagating rotational discontinuity is not. This
demonstrates the ability of the HLLD-CT_GS solver to properly capture all
seven of the emergent waves.

Figure 10: Logarithmic contour of the rest-mass density in the two-dimensional
Riemann problem at $t=0.9$ with the HLLC solver (left) and HLLE solver
(right). We employ RK4 and 3rd-order PPM cell reconstruction.

Figure 11: Rest-mass density profile in the cylindrical blast wave problem at
$t=18$ with the HLLC solver (left) and HLLE solver (right). We employ RK4 and
3rd-order PPM cell reconstruction.

Figure 12: Thermal pressure profile for the magnetized current sheet problem
at three different times: $t=10.03t_{A}$ (left), $t=30.08t_{A}$ (middle), and
$t=50.06t_{A}$ (right). Numerical results were obtained with the HLLD-CT_GS
(top row), the HLLD-CT_HLLE (middle row), and the HLLE-CT_HLLE solvers (bottom
row), respectively. The white curves denote the magnetic-field lines. We
employ RK4 with 3rd-order PPM reconstruction and a resolution of
$(N_{x},N_{y})=(512,256)$ in all simulations.

###  IV.2 Special relativistic multi-dimensional problems

####  IV.2.1 Hydrodynamics: two-dimensional shock tube

For our first multi-dimensional special relativistic test problem, we consider
the two-dimensional Riemann problem first proposed in Ref. Zanna and
Bucciantini (2002). The simulation domain spans $\in[-1,1]$ in both the $x$\-
and $y$-directions. We set $\Delta x=\Delta y=0.01$ and $N_{x}=N_{y}=100$. We
impose outflow boundary conditions in both directions. We use an adiabatic
index of $\Gamma=5/3$. We employ 3rd-order PPM reconstruction and set the CFL
number to be 0.45. Finally, the initial condition is given by

$\displaystyle\left(\rho,v^{x},v^{y},P\right)=\left\\{\begin{array}[]{l}(0.1,0,0,0.01)~{}\text{for}~{}x,y>0,\\\
(0.1,0.99,0,1)~{}\text{for}~{}x<0,~{}y>0,\\\
(0.5,0,0,1)~{}\text{for}~{}x,y<0,\\\
(0.1,0,0.99,1)~{}\text{for}~{}x>0,~{}y<0.\end{array}\right.$ (186)

Figure 10 show the logarithmic contour of the rest-mass density at $t=0.9$
with the HLLC solver (left panel) and with the HLLE solver (right panel). The
most notable difference in the solutions between the two solvers appears
around the two tangential discontinuities in the lower-left portion of the
simulation domain. With the HLLC solver (left panel), the initial tangential
discontinuities remain sharp. With the HLLE solver, on the other hand,
spurious waves propagate along each axis from the initial tangential
discontinuities due to numerical diffusion. Unlike the one-dimensional
problems, the spurious diffusion out of the initial tangential discontinuities
that occurs with the HLLE solver cannot be mitigated even when we employ 3rd-
order PPM reconstruction. Thus in this multi-dimensional test problem we
observe a qualitative difference in the solutions between the HLLC and HLLE
solvers that cannot be removed by resorting to higher-order reconstruction.

####  IV.2.2 Hydrodynamics: two-dimensional cylindrical explosion

For the second special relativistic multi-dimensional test problem, we
consider a cylindrical blast wave in two dimensions. For this problem, we
choose the simulation domain to span $x\in[-2,2]$ and $y\in[-3,3]$, and set
$\Delta x=\Delta y=0.02$, i.e., $(N_{x},N_{y})=(100,150)$ in the $x$\- and
$y$-directions, respectively. Periodic boundary conditions are imposed at the
$x$\- and $y$-boundaries. We set the adiabatic index to $\Gamma=5/3$, employ
3rd-order PPM reconstruction, and set the CFL number to 0.45. The initial
condition is given by

$\displaystyle\rho=1,~{}P=\left\\{\begin{array}[]{l}2.5~{}\text{for}~{}\sqrt{x^{2}+y^{2}}<0.5,\\\
0.1~{}\text{for}~{}\sqrt{x^{2}+y^{2}}>0.5.\\\ \end{array}\right.$ (189)

Figure 11 shows the rest-mass density profile at $t=18$ with the HLLC solver
(left panel) and the HLLE solver (right panel). By this time, the blast wave
has intersected itself many times, and consequently a Rayleigh-Taylor-like
instability (known in this context as the Richtmyer-Meshkov instability) has
developed White  _et al._ (2016). With the HLLC solver, the Richtmyer-Meshkov
instability is well resolved, and as a result the density irregularity around
the elliptical figure is sharply captured. By contrast, with the HLLE solver,
the fine structure around the elliptical figure is not captured well due to
the large numerical diffusivity. This demonstrates an effective improvement in
spatial resolution with the HLLC solver compared to that with the HLLE solver.

####  IV.2.3 Magnetohydrodynamics: two-dimensional magnetized current sheet

Next we consider a two-dimensional problem in relativistic
magnetohydrodynamics: that of a magnetized current sheet, studied recently by
Refs. Mignone and Del Zanna (2021); White  _et al._ (2016). The initial
profile for the magnetic field is given by

$\displaystyle B^{x}=B_{0}\tanh\left(\frac{y}{a}\right),$ (190)

where $B_{0}=1$ and $a=0.04$. The density is uniform with $\rho=1$ and the
fluid is at rest with $v^{i}=0$. The thermal pressure is determined from the
force balance with the magnetic pressure and its profile is given by

$\displaystyle P=\frac{B_{0}^{2}}{2}\left(\beta+1\right)-\frac{B_{x}^{2}}{2},$
(191)

where $\beta$ is the initial plasma-beta parameter, which we set to unity. The
equilibrium magnetic field is initially perturbed and the perturbation is
given by the $z$-component of the vector potential as

$\displaystyle\delta A_{z}=\epsilon
B_{0}\cos\left(\frac{k_{y}y}{2}\right)\cos\left(k_{x}x\right),$ (192)

where $k_{x}=2\pi/L_{x}$, $k_{y}=2\pi/L_{y}$, $\epsilon=10^{-3}$, and $L_{x}$
and $L_{y}$ denote the domain size in the $x$\- and $y$-directions,
respectively. We employ a simulation domain consisting of $x\in[-0.5,0.5]$ and
$y\in[-0.25,0.25]$. To check convergence, we carry out simulations at three
different resolutions: $(N_{x},N_{y})=(512,256)$, $(256,128)$, and $(128,64)$
. We set the CFL number to 0.8 in all simulations. We impose a periodic
boundary condition in the $x$-direction, and a reflective boundary condition
in the $y$-direction. With this setup, the maximum Alfvén wave speed is
$\approx 0.557$ and the Alfvén timescale is $t_{\rm A}\approx 1.78$.

Figure 13: Magnetic-field energy dissipation in the magnetized current sheet
problem as a function of time. Blue, green, and cyan curves denote numerical
solutions with the HLLD-CT_GS, HLLD-CT_HLLE, and HLLE-CT_HLLE solvers,
respectively. The solid, dashed, and dotted curves denote a resolution of
$(N_{x},N_{y})=(512,256)$, $(256,128)$, and $(128,64)$, respectively. We
employ RK4 and 3rd-order PPM cell reconstruction in all runs.

Figure 12 displays colorplots of the thermal pressure together with the
magnetic-field lines at three different times: $t=10.03t_{A}$ (left panel),
$t=30.08t_{A}$ (center), and $t=50.06t_{A}$ (right). The top, middle, and
bottle panels show the numerical solutions with the HLLD-CT_GS, HLLD-CT_HLLE,
and HLLE-CT_HLLE solvers, respectively. The snapshots are all taken from our
highest resolution runs with $(N_{x},N_{y})=(512,256)$. Magnetic field lines
reconnect at $y\approx 0$ due to the numerical resistivity inherent both in
the Riemann solvers as well as in the constrained transport scheme. Once
reconnection starts, the profile of the magnetic-field lines changes, and as a
result, the thermal pressure profile is modified, leading to the formation of
island-like structures.

The timescale of the reconnection depends on how large the numerical
resistivity is. Figure 12 indicates that the HLLD-CT_GS solver is accompanied
with the smallest numerical resistivity because the formation of the islands
is delayed. It is found that HLLD-CT_HLLE solver has the largest numerical
resistivity, leading to rapid formation of the islands. This does not agree
with one’s naive expectation, because the HLLE-CT_HLLE solver is actually less
dissipative than the HLLD-CT_HLLE solver. In other words, we observe an
unexpected hierarchy between the HLLD-CT_HLLE and HLLE-CT_HLLE solvers. This
stems from the algorithm of the CT_HLLE solver. In this constrained transport
scheme, dissipation terms which are proportional to the maximum absolute value
of the characteristic speed appear in the electric-field evaluation (see,
e.g., Eq. (44) in Ref. Del Zanna  _et al._ (2003)). This characteristic speed
is then obtained from the (global) Riemann solver. We find that the HLLD
solver returns a larger characteristic speed than the HLLE solver. As a
result, the HLLD-CT_HLLE solver ends up being more diffusive than the HLLE-
CT_HLLE solver, as can be seen in this test problem.

Figure 13 shows the fraction of the initial magnetic-field energy that is
dissipated as a function of time. With the HLLD-CT_GS solver (blue curves),
the magnetic-field energy dissipates only gradually. Also, the dissipation
rate is suppressed when we employ higher resolution: the energy increases by
an order of magnitude only over $50$ Alfvén timescales. This feature is also
found for the HLLE-CT_HLLE solver (cyan curves), although the dissipation rate
steeply rises at a later time, $t\approx 40t_{A}$, even in our highest
resolution run. With the HLLD-CT_HLLE solver (green curves), magnetic
reconnection commences immediately after the simulation starts. We conclude
that for problems involving strong magnetic field gradients (current sheets)
accurate evolution can be modeled only when the HLLD solver is paired with
CT_GS for the constrained transport.

Figure 14: Perturbed velocity difference $\Delta
v^{y}\equiv(v_{\mathrm{max}}^{y}-v_{\mathrm{min}}^{y})/2$ in the special-
relativistic magnetohydrodynamical Kelvin-Helmholtz instability as a function
of time. The blue curves denote results obtained with the HLLD-CT_GS solver
and the cyan curves denote results HLLE-CT_HLLE solver. The solid, dashed, and
dotted curves correspond to resolutions of $(N_{x},N_{y})=(200,400)$,
$(100,200)$, and $(64,128)$, respectively.  Figure 15: Density snapshots taken
near nonlinear saturation of the special relativistic magnetohydrodynamical
Kelvin-Helmholtz instability. Top row: the results with the HLLD-CT_GS solver
at three different resolutions. Bottom row: same as the top panel but with the
HLLE-CT_HLLE solver. The white lines indicate velocity field streamlines. The
left, middle, and right panels show the low, middle, and high resolution runs,
respectively.

####  IV.2.4 Magnetohydrodynamics: two-dimensional Kelvin-Helmoltz instability

The second two-dimensional problem in special relativistic
magnetohydrodynamics is the Kelvin-Helmholtz instability, as proposed in Refs.
Mignone  _et al._ (2009); Bucciantini and Del Zanna (2006). For this, we
prepare a simulation domain which spans $x\in[-0.5,0.5]$ and $y\in[-1,1]$. To
check the convergence, we perform the simulations with three different
resolutions: $(N_{x},N_{y})=(200,400)$ (‘high’ resolution),
$(N_{x},N_{y})=(100,200)$ (‘medium’ resolution), and $(N_{x},N_{y})=(64,128)$
(‘low’ resolution). The simulations are carried out with either the HLLD-CT_GS
or HLLE-CT_HLLE solvers, and we employ 3rd-order PPM reconstruction for all
the simulations. We impose a periodic boundary condition in the $x$-direction,
and an outflow boundary condition in the $y$-direction. The CFL number is set
to 0.4 in all the simulations.

As the initial condition, we give a tanh-shaped shear velocity profile for the
$x$-component,

$\displaystyle v^{x}=-v_{\text{sh}}\tanh(y/a),$ (193)

where $v_{\text{sh}}=0.25$ and $a=0.02$. The thickness $a$ of the shear layer
is covered by around $2,4,$ and $8$ grid cells at the low, medium, and high
resolutions, respectively. We employ a uniform density of $\rho=1$, and a
uniform gas pressure with $P=20$. The adiabatic index is taken to be $4/3$.
Note that our setup is different from that employed in the recent test
simulation for the Kelvin-Helmholtz instability in special relativistic
magnetohydrodynamics of Ref. Mattia and Mignone (2021), in which the authors
employ a non-uniform density field, a smaller shear-layer thickness of
$a=0.01$, and an amplitude of the $x$-component of the velocity
($v_{\text{sh}}=0.5$) which is twice that used in our runs.

The magnetic field at $t=0$ is given by

$\displaystyle(B^{x},B^{y},B^{z})=\left(\sqrt{2\sigma_{\text{pol}}P},0,0\right),$
(194)

i.e. the magnetic field is initially uniform and parallel to the velocity in
the lower-half of the $xy$-plane. We set $\sigma_{\text{pol}}=0.01$. The shear
layer is perturbed by the motion in the $y$-direction as

$\displaystyle v^{y}=\frac{1}{40000}\sin(2\pi x)\exp\left(-100y^{2}\right),$
(195)

while $v^{z}=0$.

Figure 14 shows the perturbed velocity difference $\Delta
v^{y}\equiv(v_{\mathrm{max}}^{y}-v_{\mathrm{min}}^{y})/2$ as a function of
time taken from six simulations at three different resolutions and employing
either the HLLD-CT_GS or HLLE-CT_HLLE solver. All the simulations start from
perturbations of size $\sim 10^{-5}$. We find exponential growth followed by
nonlinear saturation at the end of the linear phase at $t\sim 10$. The
behavior during the linear phase depends strongly on the solver, particularly
at low resolutions. Nonlinear saturation occurs more quickly in the
simulations with the HLLD-CT_GS solver than in those with the HLLE-CT_HLLE
solver, but the saturation amplitude depends only weakly on the solver and
resolution. The growth rate is higher with the less diffusive HLLD-CT_GS
solver than with the HLLE-CT_HLLE solver, but the results converge between the
two solvers as the resolution is improved. This result is consistent with that
in Ref. Mattia and Mignone (2021) (see their Fig. 14). The evolution after the
nonlinear saturation is not sensitive to the solver or resolution, although at
late times (not shown) the velocity difference decays more quickly in the
simulations with the (more diffusive) HLLE-CT_HLLE solver than with the HLLD-
CT_GS.

In Fig. 15 we show snapshots of the density at nonlinear saturation $t\sim 10$
from the six simulations. The top row shows results from the low, medium, and
high resolution runs using the HLLD-CT_GS solver, while the bottom row shows
the corresponding snapshots from runs that employ the HLLE-CT_HLLE solver.
Using the HLLD-CT_GS solvers, we observe the formation of a single vortex
together with two neighboring, stretched secondary vortices that are well-
resolved at all resolutions, whereas with the HLLE-CT_HLLE solver we see the
formation of only a single large vortex at the shear interface, mirroring the
behaviour of the Kelvin-Helmholtz instability in the simulations of Ref.
Bucciantini and Del Zanna (2006) which employed the HLLE-CT_HLLE solver. Our
results show that, at least at low resolutions, the HLLE solver is not
appropriate for studying phenomena in which the Kelvin-Helmholtz instability
plays an important role.

###  IV.3 General relativistic problems in a fixed background spacetime

Figure 16: Top: Radial rest-mass density profile for (hydrodynamic) Bondi flow
in a non-rotating black hole spacetime at $t=22.5M$ (the end of the
simulation). The mass accretion rate is fixed at $\dot{M}_{\text{acc}}=0.797$.
The solid black curve indicates the analytical solution, while the blue and
green dots denote the numerical solution obtained with the HLLC and HLLE
solvers, respectively. The inset show the solution near the inner boundary. In
the inset, the red and yellow dots denote the location of the horizon and of
the inner boundary, respectively. Middle: Same as the top panel, but showing
the radial velocity profile. Bottom: $L_{1}$ norm of the error in the rest-
mass density as a function of the spatial grid spacing. The blue and green
dots denote the error of the numerical solutions with the HLLC and HLLE
solvers, respectively. The dotted line denotes 2nd-order convergence.

####  IV.3.1 Hydrodynamics: Bondi flow

As a test problem in a curved (but static) spacetime, we consider spherical
accretion (ingoing Bondi flow) onto a non-rotating black hole Hawley  _et al._
(1984). The Bondi flow in Schwarzschild coordinates has been extensively
discussed in the literature (see, e.g., Ref. White  _et al._ (2016)).
Following previous work Gammie  _et al._ (2003); Shibata and Sekiguchi (2005),
we adopt the parameters for this problem as follows: an adiabatic index of
$\Gamma=4/3$, an adiabat of $K=1$, and a critical radius of $r_{\text{c}}=8M$,
where $M$ denotes the black hole mass. With this setup, the mass accretion
rate $\dot{M}_{\text{acc}}$ is 0.797. We perform simulations both with the
HLLE and HLLC solvers, and employ 3rd-order PPM reconstruction.

Our numerical-relativity code employs the so-called puncture formalism, and
hence, in the presence of black holes, the black-hole spacetime is foliated in
most cases by the so-called limiting hypersurface Estabrook  _et al._ (1973).
Thus, for preparing a practical setup in this test problem, a non-rotating
black hole should be described in the so-called maximal trumpet geometry
rather than in Schwarzschild coordinates or in isotropic coordinates on slices
of constant Schwarzschild time Baumgarte and Naculich (2007); Miller and
Baumgarte (2017). Note that in both of these latter two coordinate systems,
the fluid four-velocity exhibits pathological behavior near the horizon Miller
and Baumgarte (2017) 333The pathological behavior at the horizon is avoidable
if one employs Kerr-Schild coordinates. However, we employ the maximal trumpet
geometry in this test problem because of its high compatibility with our
numerical relativity code.. In Appendix B, we describe the explicit coordinate
transformation from the Schwarzschild coordinates to the maximal trumpet
geometry. With this geometry, the radial component of the shift vector is non-
zero. Therefore, the tetrad basis (see, e.g., Eq. (60)) does not agree any
longer with a coordinate basis in the Minkowski spacetime, and the cell
interface may be dragged by the shift vector as discussed in Sec. III.2.

We employ a simulation domain in Cartesian coordinates spanning
$x,y,z\in[0,L]$ with $L=10M$. The grid spacing of the simulation is given by
$\Delta=\Delta x=\Delta y=\Delta z=0.1M$ with $N=N_{x}=N_{y}=N_{z}=100$ as the
number of grid cells in each direction. We also check convergence by
increasing the resolution to $N=200$ and $N=400$, which correspond to grid
spacings of $\Delta=0.05M$ and $0.025M$, respectively. We set the CFL number
to 0.45 and integrate the numerical solution up to $t=22.5M$. We impose a
stationary boundary condition at the outer and inner boundaries, with the
latter located at $r_{\text{in}}=0.4M$. Note that the horizon in this geometry
is located at $r_{\text{BH}}\approx 0.78M$. We also impose octant symmetry at
the $x,y$, and $z=0$ planes.

Figure 16 shows radial profiles of the rest-mass density and the radial
velocity calculated by the HLLC solver with the blue dots and by the HLLE
solver with the green dots on top of the analytic solution Hawley  _et al._
(1984). The profiles are along the diagonal direction, i.e., $x=y=z$ in the
simulation domain. This figure demonstrates that our implementation of the
HLLC solver in curved spacetime works well. It also shows that, for this
particular problem, the HLLE solver works as well as the HLLC solver because
of the smoothness of the accretion flow, as many other previous
implementations have shown; e.g., Refs. Gammie  _et al._ (2003); Shibata and
Sekiguchi (2005); Mösta  _et al._ (2014); Etienne  _et al._ (2010).

In the lower panel of Fig. 16 we plot the $L_{1}$ norm of the error in the
rest-mass density as a function of the spatial grid spacing. The convergence
order of the $L_{1}$ norm of the error is $\approx 2$ both for the HLLC and
HLLE solvers, because our Riemann solver is 2nd-order accurate. One likely
reason for the slight deviation from the expected accuracy is that spherical
symmetry of the accretion flow is not perfectly preserved during the evolution
because we simulate it in the Cartesian geometry. This plot also shows that
the numerical solution with the HLLC solver is more accurate than that with
the HLLE solver. Our interpretation of this is that with the tetrad
transformation (see Sec. III.2) the frame-dragging effect of the cell
interface is taken into account with a better accuracy (see also Fig. 2)
444Note that our HLLE solver (which was the only Riemann solver present in our
original formulation), does not employ the tetrad transformation, but is
instead formulated directly in a curved spacetime. See Ref. Shibata and
Sekiguchi (2005) for details..

Figure 17: Same as Fig. 16, but for a magnetized Bondi flow in a non-rotating
black hole spacetime at $t=22.5M$. The blue and cyan dots denote numerical
solutions with the HLLD-CT_GS and HLLE-CT_HLLE solvers, respectively. The
bottom-left panel shows the radial magnetic field profile.

####  IV.3.2 Magnetohydrodynamics: Magnetized Bondi flow

The next test problem in a curved spacetime is magnetized Bondi flow onto a
non-rotating black hole. It is known that a purely radial magnetic field does
not alter the flow profile of non-magnetized Bondi flow Shibata (2016).
Therefore, we employ the same flow profile used in the previous section. From
the divergence-free condition (17), the radial magnetic field should be
$B^{R}\propto f/R^{2}$ in Schwarzschild coordinates (see Appendix B for the
definition of $f$ and the transformation to the maximal trumpet geometry). The
amplitude of the magnetic field is chosen to be such that $b^{2}/\rho=1$ at
$R=3M$. We perform two simulations, one with the HLLD-CT_GS solver and the
other with the HLLE-CT_HLLE solver. We employ RK4 and 3rd-order PPM
reconstruction in both cases.

We employ a simulation domain in Cartesian coordinates spanning
$x,y,z\in[0,L]$ with $L=12.8M$. The grid spacing of the simulation is
$\Delta=\Delta x=\Delta y=\Delta z=0.1M$ with $N=N_{x}=N_{y}=N_{z}=128$ being
the number of the grid cells in each direction. To check convergence, we
perform better-resolved simulations with $N=256$ and $N=512$, i.e.,
$\Delta=0.05M$ and $0.025M$, respectively. We impose octant symmetry at the
$x$, $y$, and $z=0$ planes, and a stationary condition at the outer and inner
boundaries, with the latter located at $r_{\text{in}}=0.4M$. Numerical
simulations are performed up to $t=22.5M$.

Figure 17 shows the radial profiles of the rest-mass density (top-left), the
radial velocity (top-right), and the radial magnetic field (bottom-left).
Numerical solutions with the HLLD-CT_GS solver are indicated by the blue dots,
while those with the HLLE-CT_HLLE solver are indicated by the cyan dots. As in
the non-magnetized cases, the flow profiles agree with the analytic solution
Hawley  _et al._ (1984) (see also the insets in Fig. 17 which show the
solution close to the inner boundary). The rest mass density inside the
horizon slightly deviates from the analytic solution. However, the deviation
decreases as the spatial resolution is increased. This demonstrates that our
HLLD solver works just as well as our HLLC solver in a curved spacetime. As in
the hydrodynamic case, we find no qualitative difference in the numerical
solutions between the HLLD-CT_GS and HLLE-CT_HLLE solvers because of the
smoothness of the accretion flow. The bottom-right panel in Fig. 17 plots the
$L_{1}$ norm of the error in the rest-mass density as a function of the
spatial grid spacing. It shows that (i) the numerical solution with the HLLD-
CT_GS solver is more accurate than that with the HLLE-CT_HLLE solver, and (ii)
the order of the convergence is $\approx 2$. These results are essentially the
same as those in the previous subsection. Again, the deviation from the formal
accuracy of the Riemann solver is likely to be an artifact of the Cartesian
geometry which we employ.

##  V Application to a dynamical spacetime

Finally, we apply our new Riemann solvers in general relativity to a dynamical
spacetime. We simulate a binary neutron star merger, both with and without
magnetic fields. We turn on the solver for Einstein’s equations and the
neutrino-radiation hydrodynamics solver in the simulations shown in this
section (see Eqs. (11)–(15)).

###  V.1 Hydrodynamics: binary neutron star merger

####  V.1.1 Setup

First, we consider non-magnetized asymmetric binary neutron stars with masses
of $1.2$ and $1.5M_{\odot}$. We utilize the spectral method library LORENE
LORENE ; Gourgoulhon  _et al._ (2001); Taniguchi and Shibata (2010); Taniguchi
and Gourgoulhon (2002, 2003) to generate a quasi-equilibrium configuration of
the irrotational binary neutron star. We also employ an eccentricity reduction
prescription to generate an initial condition that has low orbital
eccentricity Kyutoku  _et al._ (2014). The initial orbital angular velocity is
set to be $m_{0}\Omega_{0}=0.028$ where $m_{0}=2.7M_{\odot}$ is the total mass
of the binary.

Our solver for Einstein’s equations implements the BSSN-puncture formulation
Shibata and Takashi (1995); Baumgarte and Shapiro (1998); Campanelli  _et al._
(2006); Baker  _et al._ (2006), locally incorporating the Z4c prescription for
constraint propagation Hilditch  _et al._ (2013). We employ 4th-order centered
finite differencing for the spatial derivative of the metric, a lop-sided
finite difference for the advection term associated with the shift vector, and
4th-order Runge-Kutta for the time integrator. For the relativistic
hydrodynamics solver, we employ either the HLLC or HLLE solver, together with
3rd-order PPM cell reconstruction.

We employ the SFHo equation of state for relatively high-density nuclear
matter Steiner  _et al._ (2013), and the Timmes (Helmholtz) equation of state
for the low-density part Timmes and Swesty (2000). Because high-resolution
shock-capturing schemes cannot treat the vacuum state, we need to implement an
artificial atmosphere outside the neutron stars. In this simulation, we set a
constant atmospheric density of $\rho_{\text{atm}}=10^{3}~{}{\rm g/cm^{3}}$
for the inner part of the finest fixed mesh refinement (FMR) domain, for which
the refinement boundary along each axis is typically located at
$L_{\text{fin}}=38.7\,{\rm km}$ (see below for the FMR setup in detail). We
also set a power-law profile of the atmospheric density of
$\rho_{\text{atm}}=10^{3}(L_{\text{fin}}/r)^{3}\,{\rm g/cm^{3}}$ for
$r>L_{\text{fin}}$ and as far as the atmospheric density is larger than the
floor value which is determined by the employed equation of state. In our
present table for the equation of state, this floor is $\approx 0.17~{}{\rm
g/cm^{3}}$ and if $\rho_{\text{atm}}$ becomes smaller than this value, we set
the the atmospheric density to the floor value. The atmospheric temperature is
set to be $10^{-3}~{}{\rm MeV}$.

We also explicitly solve the radiation-hydrodynamics equations for neutrinos
in time using an approximate neutrino-transfer scheme based on a leakage
scheme Sekiguchi (2010) and the truncated moment formalism Thorne (1981);
Shibata  _et al._ (2011). The cooling source terms are computed using a
general-relativistic leakage scheme Sekiguchi  _et al._ (2012), and heating
source terms due to neutrino capture processes are computed by the method
presented in Ref. Fujibayashi  _et al._ (2017).

The computational region consists of $13$ levels of FMR half-cubic domains.
The size of each FMR domain is $\in[-L/2^{l-1},L/2^{l-1}]$ for $x$ and $y$,
and $z\in[0,L/2^{l-1}]$ with $l=1,2,\cdots,13$. Note that in the $z$-direction
we impose reflection symmetry with respect to the equatorial plane, $z=0$. We
set the overall domain size to $L\approx 158,000~{}{\rm km}$ and
$N=N_{x}=N_{y}=N_{z}=258$. Thus the grid spacing of the finest FMR domain is
$\Delta x_{13}=\Delta y_{13}=\Delta z_{13}=150$ m. To check convergence, we
also perform simulations with lower resolutions of $N=196$ and $N=158$, for
which the grid spacing of the finest FMR domain is $\Delta x_{13}=200$ m and
$\Delta x_{13}=250$ m, respectively. For the HLLC solver, we perform an
additional simulation with $N=377$ and $\Delta x_{13}=100$ m. By virtue of the
cell-centered grid structure, the cell interface of the parent FMR domain
coincides with that of the child FMR domain. We employ the reflux prescription
during time marching of the Berger-Oliger type mesh refinement algorithm to
ensure the conservation of baryonic mass.

####  V.1.2 Inspiral phase

Figure 18: Maximum rest-mass density (left) and density-weighted Hamiltonian
constraint violation (right) as functions of time in non-magnetized
simulations of inspiraling binary neutron stars. The blue and green curves
denote the results with the HLLC and HLLE solvers, respectively. The solid,
dashed, dotted, and dot-dashed curves denote the results with (finest-level)
grid spacings of $\Delta x_{13}=100$ m, $150$ m, $200$ m, and $250$ m,
respectively.

Figure 19: (Top) Orbital separation as a function of time in hydrodynamic
simulations of a binary neutron star merger. The color code and line styles
have the same meaning as in Fig. 18. The cross symbols denote the final moment
at which two density maxima can be identified. (Bottom) Colorplots of the
rest-mass density in the orbital plane at $t=14.64$ ms with the HLLC solver
(left panel) and with the HLLE solver (right panel). Both simulations have
been run with a (finest-level) grid-spacing of $\Delta x_{13}=150$ m.

Figure 20: Same as Fig. 18, but for $t\geq 15$ ms. The inset in the left panel
shows a close-up over a short duration of the solution obtained with the HLLC
and HLLE solvers with $\Delta x_{13}=150$ m.

The left and right panels of Fig. 18, respectively, show the time-evolution
(during the inspiral phase) of the maximum rest-mass density and the density-
weighted Hamiltonian constraint violation (see Eqs. (29) and (30) in Ref.
Shibata and Sekiguchi (2005) for definitions). The blue and green curves
denote the results with the HLLC and HLLE solvers, respectively, and the
solid, dashed, dotted, and dot-dashed curves indicate the resolution (i.e.
$\Delta x_{13}=100$ m, $150$ m, $200$ m, and $250$ m, respectively). During
inspiral, the maximum rest-mass density oscillates due to numerical error
regardless of which solver is used. It also decreases partly due to the
numerical error. However, the degree of the decrease is much more prominent
with the HLLE solver than it is with the HLLC solver, especially at the
coarsest resolution. This is due to the large numerical diffusion inherent in
the HLLE solver. Specifically, this solver is more subject to spurious
broadening of the density profile near the stellar surface (not shown),
leading to a higher degree of spurious neutron-star expansion and to a
resultant decrease in the maximum rest-mass density. However, this artifact is
mitigated with the HLLC solver, because of its stronger capability of
capturing irregular surfaces, i.e., the stellar surface.

The right panel of Fig. 18 shows that, for a given grid resolution, the time-
averaged value of the constraint violation during the inspiral phase is
smaller with the HLLC solver than with the HLLE solver. This demonstrates that
the numerical result with the HLLC solver is more accurate than that with the
HLLE solver. We find that the order of convergence of the density-weighted
Hamiltonian constraint violation is $1.7$–$1.8$, irrespective of which Riemann
solver is used. Note that this convergence is slow compared to that achieved
using a higher-order finite difference scheme Bernuzzi and Dietrich (2016);
Radice  _et al._ (2014); Most  _et al._ (2019), but could be improved if we
were to employ a more accurate reconstruction scheme such as MP5 Suresh and
Huynh (1997). However, the implementation of such a scheme is beyond the scope
of this paper.

The top panel of Fig. 19 shows the orbital separation of the binary as a
function of time. Here we define ‘orbital separation’ as the coordinate
distance in the orbital plane between the two rest-mass density maxima. The
cross symbols denote the final moment at which we can unambiguously identify
the two rest-mass density maxima. At this point the less massive neutron star
has been significantly tidally elongated, and we define this time as being the
time of onset of the merger. This plot shows that the merger time found in the
simulation with the HLLC solver is later than that with the HLLE solver (the
reason for this will be described shortly). The bottom panels display contour
plots of the rest-mass density in the orbital plane at the moment of merger
for the runs with the HLLC solver (left panel) and the HLLE solver (right
panel) (for both cases, $\Delta x_{13}=150$ m). The orbital phase with the
HLLE solver slightly larger compared to that with the HLLC solver. This
implies that the neutron star simulated with the HLLE solver is more subject
to artificial tidal deformation than the neutron star with the HLLC solver,
because the HLLE solver (since it cannot accurately resolve the irregularities
at the stellar surface) results in a larger spurious expansion of the neutron
star. Note that the tidal elongation of the low-density part of the less
massive neutron star is more enhanced with the HLLE solver than with the HLLC
solver, as found from the comparison of the two contour plots. It is this
enhanced (but artificial) tidal elongation with the HLLE solver that
ultimately results in the earlier merger time observed when when we employ
that solver.

Figure 21: (Left) Dimensionless spin of the black hole as a function of time
in non-magnetized binary neutron star merger simulations. (Right)
Gravitationally bound baryonic mass outside the apparent horizon as a function
of time in non-magnetized binary neutron star merger simulations. The color
code and the line style of the legend are the same as in Fig. 18.

Figure 22: Luminosity of electron neutrinos (left) and electron antineutrinos
(right) as functions of time in non-magnetized binary neutron star mergers.
The color code and the line style of the legend are the same as in Fig. 18.
Figure 23: Dynamical ejecta mass as a function of time. In the inset, the
vertical axis is shown using a logarithmic scale and the shaded region
represents the violation of the baryonic mass conservation in a run using the
HLLC solver and at a resolution of $\Delta x_{13}=100$ m. The color code and
the line style of the legend are the same as in Fig. 18.

####  V.1.3 Post-merger phase

Having presented various diagnostics from the inspiral phase, we now turn our
attention to the post-merger phase. Figure 20 shows the maximum rest-mass
density and the density-weighted Hamiltonian constraint violation as functions
of time during the post-merger phase. The existence of oscillations in the
density after the merger indicates the formation of a massive neutron star
remnant rather than a direct collapse to a black hole. The remnant massive
neutron star gradually contracts due to angular momentum transport by the
gravitational torque, neutrino cooling, and gravitational-wave emission.
Eventually the remnant collapses to a black hole, indicated by the blow-up of
the maximum rest-mass density (see the figure at $30$–$40$ ms). Because the
remnant massive neutron star is a meta-stable object, its stability is
influenced significantly by the numerical truncation error, by the randomness
with which the collapse time does not become a monotonic function of the
spatial grid spacing. (The non-monotonic behavior of the black-hole formation
time with respect to the grid spacing is also reported in the literature (see,
e.g, Ref. Baiotti and Rezzolla (2017)).) Specifically, with the HLLC solver,
the collapse time of the remnant coincides for both the $\Delta x_{13}=150$ m
and $\Delta x_{13}=250$ m runs, is earliest for the $\Delta x_{13}=100$ m run,
and is latest for the $\Delta x_{13}=200$ m run. With the HLLE solver, the
collapse time is earliest for the $\Delta x_{13}=200$ m run and the latest for
$\Delta x_{13}=250$ m run. The inset in the left panel of Fig. 20 shows a
close-up of the results with the HLLC and HLLE solvers for $\Delta x_{13}=150$
m.

However, the collapse time of the remnant is systematically earlier for runs
with the HLLC solver. This is related to the evolution of the oscillation
amplitude of the remnant neutron star. For $t\lesssim 20$ ms, the oscillation
amplitude of the maximum rest-mass density is approximately identical for the
two solvers (see the left-hand panel of Fig. 20). After that, however, the
oscillations are noticeably damped when we use the HLLE solver. This implies
that the oscillation energy is dissipated by the numerical diffusion inherent
in the HLLE solver. Thus, the lifetime of the remnant massive neutron star is
significantly overestimated when the more diffusive HLLE solver is used.

The right panel of Fig. 20 shows that the density-weighted Hamiltonian
constraint violation is of order $10^{-4}$ during the remnant massive neutron
star phase and of order $10^{-1}$ after the black hole formation. The
constraint violation only slowly decreases with increased resolution in the
post-merger phase, and does so regardless of which solver is used. The reason
for this is that during the merger phase, shocks are formed inside a large
portion of the neutron stars. Because shocks are always computed with first-
order accuracy in numerical hydrodynamics, the overall accuracy of the
solution deteriorates and the convergence becomes slow.

The left panel of Fig. 21 shows the evolution of the dimensionless spin of the
remnant black hole 555The methods we use to estimate the non-dimensional spin
and the mass of the black hole are the same as in Refs. Kiuchi  _et al._
(2009); Shibata and Taniguchi (2008).. We find spurious spin-down of the black
hole due to numerical diffusion, in particular when the simulations are
performed at low resolutions with $\Delta x_{13}=200$ m or $250$ m. We measure
the spin-down rate in the HLLC run and estimate that the dimensionless spin
decreases by $\gtrsim 0.1$ in 1 s if $r_{\rm AH}/\Delta x_{13}\lesssim 15$
where $r_{\rm AH}$ denotes the minimum radius of the apparent horizon.
However, the spurious spin-down rate decreases approximately at the 4th order,
reflecting the order of the accuracy in the solver for Einstein’s equations.
This implies that the spurious decrease of the dimensionless spin will be
suppressed to the required level if we perform a simulation with a
sufficiently high resolution. In low-resolution runs, however, the spurious
spin down will influence the evolution of the disk because the specific
angular momentum at the inner stable circular orbit will increase as a result