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

We implement advanced Riemann solvers HLLC and HLLD \cite{Mignone:2005ft,MUB:2009} together with an advanced constrained transport scheme \cite{Gardiner:2007nc} 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 \cite{DelZanna:2002rv}, 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 [5,6]. 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 Λ 800, with an accurate measurement of the total mass of the binary yielding 2.73 +0.04 −0.01 M [5, 7-9] [10]. 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 [6,[11][12][13], thus providing the first 'smoking gun' for supporting the hypothesis that binary mergers can be the central engine of short gamma-ray bursts [14][15][16][17]. 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) [16,[18][19][20]. It had been predicted that the r-process nucleosynthesis subsequently causes so-called kilonova emission via the radioactive decay of unstable r-process elements [21,22], and a kilonova was indeed observed after the merger in the near infrared, optical, and ultraviolet bands [23][24][25][26][27][28][29][30][31][32][33][34][35][36][37].
References [38,39] 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 , which is much larger than the total mass measured in binary pulsars observed in our Galaxy [40]. The merger dynam-ics, mass ejection process, and resultant electromagnetic emission due to r-process nucleosynthesis could be different from those in GW170817 and AT 2017gfo [41,42]. Although an electromagnetic counterpart was not observed in GW190425, either due to poor sky localization or due to intrinsically dimmer emission [43][44][45], 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 [46], 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 [47].
Recent axisymmetric neutrino-radiation viscoushydrodynamics 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 [47][48][49][50][51]. 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 threedimensional 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 [52] (see also Refs. [53,54] 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 [49][50][51], although the electron fraction of the post-merger ejecta is appreciably larger than that of the dynamical ejecta.
Furthermore, very recently we performed neutrinoradiation magnetohydrodynamics simulations of black hole-neutron star mergers in numerical relativity [55]. 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 [55].
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 [56].
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 [57]. 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 [58], 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 [59].
At present, the Riemann solver and constrained transport scheme implemented in existing numerical relativity magnetohydrodynamics codes, e.g. [55,[60][61][62][63][64][65][66][67], are based on the HLLE solver [4,59,68]. (An exception is the SpECTRE [69] code, which is based on the discontinuous Galerkin method.) This Riemann solver is known to be very diffusive [2,70,71]. 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. [61,62,72] 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 [73,74] based on the finite volume method. We implement the HLLC solver for relativistic hydrodynamics, which restores the contact discontinuity [1], and the HLLD solver for relativistic magnetohydrodynamics, which takes into account five of the seven waves in the Riemann fan [2]. Both these Riemann solves are known to be less diffusive than the HLLE solver [4]. In addition, the constrained transport scheme in Ref. [3], which relies on the solution given by a Riemann solver, significantly suppresses numerical diffusion compared to the HLLE-constrained transport scheme proposed in Ref. [4] (see Ref. [75] 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. [3] 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 [1], that of the HLLD solver [2], and the electric field evaluation (which is used by the constrained transport algorithm) [3]. In Sec. IV, we validate our implementation of the new Riemann solvers by performing one-and multidimensional 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 ≈ 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. [76].
We begin by introducing a unit vector normal to a spatial hypersurface of constant coordinate time, t, where α and β i are the lapse function and shift vector, respectively. With this vector, the four dimensional metric can be decomposed into where γ ij is the three-dimensional spatial metric.
The stress-energy-momentum tensor for ideal magnetohydrodynamics and for a free-streaming neutrinoradiation field are, respectively, given by T µν (MHD) = ρhu µ u ν + P g µν + u µ u ν + T µν (Rad,s,νi) = E (νi) n µ n ν + F µ (νi) n ν + F ν (νi) n µ + P µν (νi) , (2.3) where ρ, P , u µ , b µ , E (νi) , F µ (νi) , and P µν (ν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 ν i in the Eulerian frame. h = 1 + ε + P/ρ denotes the relativistic specific enthalpy with ε the specific internal energy. We consider the electron neutrino ν e , electron antineutrinoν e , and the total of µ and τ neutrinos and antineutrinos collectively denoted by ν x [74,76]. 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-energymomentum tensor of the trapped neutrinos is then absorbed into that for the ideal magnetohydrodynamics fluid because trapped neutrinos are strongly coupled to the fluid [74,76].
The conserved mass density, total momentum density, and total energy density of an electrically conducting fluid are defined by where w ≡ −n µ u µ = α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 µ n µ = 0 (i.e., B t = 0). The relation between b µ and B i is given by and thus, The equations of motion of ideal magnetohydrodynamics and of the free-streaming neutrino-radiation field are derived from the conservation of the stress-energymomentum 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 where L = e, ν e ,ν e , and ν x denotes electrons, electron neutrinos, electron antineutrinos, and heavy neutrinos, respectively. Y L and γ L denote the fractions with respect to the baryon and the source term for the number of the species L, respectively. G (νi,leak) ν is an interaction term between the fluid and free-streaming neutrino-radiation field of the neutrino species ν i in the framework of a general relativistic neutrino leakage scheme [74,77]. Here * F µν is the Hodge dual of the Faraday tensor, which is given by Equations (2.9), (2.11), and (2.13) can be written in conservative form as 14) where the flow quantities are given by the state vector 7]. The corresponding fluxes are given by , and the source terms are , and the spatial components of the stress-energy-momentum tensor are given by We also introduce the conformal metricγ ij = ψ −4 γ ij and the trace-free conformal extrinsic curvatureÃ ij = ψ −4 K ij − 1 3 Kγ ij , where ψ and K ij are the conformal factor and the extrinsic curvature, respectively. The explicit forms for γ L and G (νi,leak) µ and for the equation of motion of the free-streaming neutrino-radiation field can be found in Refs. [76,78]. The high resolution shock capturing scheme for the neutrino-radiation field Eq. (2.10) is the same as that in Ref. [79].

III. NUMERICAL ALGORITHM
In this section, we describe the numerical algorithms which we implemented in our code. In Sec. III A we present the finite volume algorithm and discretization scheme, and in Sec. III B 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 C and III D, respectively. Finally, the evaluation of the electric field used by the constrained transport algorithm is discussed in Sec. III E.
A. Finite volume method

Fluid and magnetic field at cell center
Let Ω be a region of a given four-dimensional manifold M, bounded by a closed three-dimensional surface ∂Ω, where ∂Ω denotes the surface of a four-dimensional parallelepiped composed of two spacelike surfaces {Σ t , Σ t+∆t } and three sets of two timelike surfaces {Σ x i , Σ x i +∆x i } that connect the two temporal slices [80]. The timelike surface, e.g., Σ x , may also be regarded as a time series of constant-(t, x) surfaces, S x (t). We integrate Eq. (2.14) over the domain of Ω: where dΩ = √ −gdtdxdydz.
Using Gauss's theorem, this equation can be integrated to give ∆V ≡ √ γdxdydz, (3.4) are, respectively, the three-dimensional proper volumeaveraged conserved quantities and the proper volume.
Let us now define a cell consisting of [x j − ∆x/2 : x j + ∆x/2] × [y k − ∆y/2 : y k + ∆y/2] × [z l − ∆z/2 : z l + ∆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 where t n+1 = t n +∆t. With this numerical flux, Eq. (3.2) can be discretized as where 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γ, this equa-tion is reduced to where and

Magnetic fields at cell surface
To ensure that the divergence-free condition (2.15) is maintained, we employ the constrained transport method introduced by Evans and Hawley [81]. 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. (2.14) for A ∈ [5,7] on Σ x i +∆x i . For example, through the surface Σ z+∆z , we have Levi-Civita tensor, and dS Ωz = √ −gdtdxdy. Using Stokes' theorem, this equation is integrated to give is the surface-averaged magnetic field. Similarly, through the surfaces Σ x+∆x and Σ y+∆y , respectively, we have We next consider a cell surface consisting of [x j − ∆x/2 : x j +∆x/2]×[y k −∆y/2 : y k +∆y/2], [y k −∆y/2 : y k + ∆y/2] × [z l − ∆z/2 : z l + ∆z/2], [x j − ∆x/2 : x j + ∆x/2] × [z l − ∆z/2 : z l + ∆z/2] and a numerical flux which approximates a time-averaged electric field at the cell edge, given by . (3.32) 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. (3.10) by

B. Tetrad frame
To evaluate the numerical fluxes through cell interfaces (e.g. Eq. (3.5)), we implement HLL-type Riemann solvers [1,2]. 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. [82]), 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. [83], we define a tetrad basis in the xdirection, for example, by With this basis, we can perform a transformation from the Eulerian frame to the tetrad frame by where V µ and Q µν denote a covariant vector and tensor, respectively, in the Eulerian frame. The covariant components of the tetrad basis are With this basis, we can then perform the transformation from the tetrad frame to the Eulerian frame by With this tetrad basis the procedure to obtain the numerical flux F x A j+ 1 2 ,k,l is as follows: first, we calculate the tetrad component of u (î) , v (î) , and B (î) by Second, we solve a Riemann problem in the locally Minkowski spacetime to obtain the numerical flux f (x) A j+ 1 2 ,k,l and the conserved quantities (q A ) j+ 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 xf (x) 7 j+ 1 2 ,k,l , (3.61) whereî =x,ŷ,ẑ are contracted with i = 1, 2, 3, respectively, in the second term of the right-hand side of Eqs. (3.55) and (3.56). Note that, from now on, we do not distinguish the upper-and lower-spatial tetrad components, e.g., B (î) = B (î) . These numerical fluxes are used to update the conserved quantities in Eq. (3.10). An interface velocity is calculated by [83] v This velocity is used to calculate a numerical flux at the interface (see Eqs. (3.67) and (3.87) in the next section). The tetrad basis and numerical fluxes in the y-and zdirections are summarized in Appendix A.

C. HLLC solver for relativistic hydrodynamics
In the absence of electromagnetic fields, Eq. (3.10) with A ∈ [0, 4] are reduced to those of relativistic hydrodynamics. In this case, one choice for the Riemann solver x z y (B y ) j,k+1/2,l (E z ) j+1/2,k+1/2,l (B z ) j,k,l+1/2 Schematic of a cell, cell interface, and cell edge for the finite volume method with the constrained transport method.
Fluid quantities, (QA) j,k,l , are defined at the cell center.
The electric field components, is the HLLC solver proposed in Ref. [1]. We calculate the in the tetrad frame by solving the source-free one-dimensional conservation law: where ∂ (μ) ≡ e (μ) µ ∂ µ . Given an initial condition described by for x j ≤ x ≤ 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) 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 λL, the contact discontinuity with λc, and the right-going nonlinear wave with λR, propagate from the discontinuity located at x j+ 1 2 , where λ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 λL/λR, the left/right-propagating Alfvén wave with characteristic speed λaL/λaR, and the contact discontinuity with λ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 atx j+ 1 2 may move with an interface (3.69) and λ L/R is the characteristic speed of the left/rightgoing nonlinear wave. Equation (3.69) is obtained from the jump condition and λ 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 λ c [1]: where ρ HLL H , J HLL (x) , F HLL ρH , and F HLL J (x) denote conserved quantities and fluxes in the HLL state: Once we obtain the speed of the contact discontinuity λ c , the pressure in the intermediate state is determined by Then the conserved quantities in the cL and cR states are given by where the subscripts cL and cR on the left-hand side of the equations correspond to L and R on the righthand side, respectively. These quantities can be used to evaluate the flux in the cL/cR state (3.69) and the flux in the Eulerian frame (see, e.g., Eq. (3.54)).
For the left and right characteristic speeds λ L/R , we apply Davis's estimate [1]: λ R = max(λ + (q L ), λ + (q R )), (3.80) and The equivalent expressions in the y-and z-directions are given by permutation of the indices x, y, and z.

D. 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. [2]. For this case, we calculate the HLLD flux, f (x) A j+ 1 2 ,k,l in the tetrad frame by solving the onedimensional conservation law: Here, q A has seven components (A = 0, 1, 2, 3, 4, 6, 7), and P tot ≡ P + b 2 /2 is the total pressure (gas plus magnetic). Note that the equation for B (x) is simply ∂ (t) B (x) = 0, and thus, B (x) is constant for the Riemann problem of the x-direction. Together with the initial condition given by Eq. (3.66) for the relevant components, the full magnetohydrodynamics Riemann fan consists of seven waves separating eight states [80]. 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 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) f (x) 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.

Characteristic speeds
For the fast waves, an approximate characteristic speed proposed in Refs. [82,84] is given by For the Alfvén wave, the characteristic speed is given by and for the contact wave by Given left-and right-state quantities, we first calculate the following quantities which should be preserved when one crosses the fast waves: whereî =x,ŷ,ẑ for i = 1, 2, 3, respectively, in Eq. (3.98). Alsok =ŷ,ẑ for k = 6, 7, respectively, in Eq. (3.100). For the above quantities, we employ the characteristic speed defined by

aL/aR state
Given an initial guess for the unknown total pressure P tot (which should be constant inside the Riemann fan), the three velocities in the aL and aR states are given by Note that aL and aR on the left-hand side of Eqs. (3.102)-(3.104) correspond to L and R for R J (î) , R ρH , R B (k) , and λ on the right-hand side of the same equations, respectively. With these velocity components, the magnetic field is calculated from the jump condition by fork =ŷ,ẑ. The total enthalpy density is calculated by The conserved quantities necessary for the numerical flux in Eq. (3.89) and in the Eulerian frame (see Eqs. (3.54)-(3.61)) are calculated by where the plus (minus) sign corresponds to the right (left) state. We then define K (k) by Here K (x) is nothing other than the Alfvén wave speed in the x-direction. From the jump condition one can find that K (î) , ρh tot , D/wσ (0) , and η do not change across the Alfvén waves. Therefore, η, K (î) , and the total enthalpy density are calculated by 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 122) and the characteristic speed is We impose the continuity condition on the normal velocity across the contact discontinuity, i.e., v aL . This equation gives an improved guess of the total pressure in the next iteration step. Then we go back to Eq. (3.102) and repeat the same procedure until it converges with sufficient accuracy. In practice, we employ the Newton-Raphson method to solve Eq. (3.124).
The conserved quantities necessary for the numerical flux in Eq. (3.90) and in the Eulerian frame (see Eqs. (3.54)-(3.61)) are The equivalent expressions in the y-and z-directions are given by permutation of the indices x, y, and z.

E. 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 [3] 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 bỹ . Similarly,ṽ y j,k+ 1 2 ,l andẼ z j,k+ 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. (3.33)-(3.35) 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. (3.30)-(3.32).
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 [4,82], is referred to as CT HLLE [82]. 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 [82]. 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 A). Next, in Sec. IV B, 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 C 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 Γ-law equation of state given by We also employ a cell-centered grid structure in which the x-coordinate [85] is given by and grid spacing ∆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) [82,86]. For the PPM reconstruction, we employ the min-mod limiter function with a compression parameter which is generally set to b = 2 [82], though in some cases we employ different values of b.
A. Special relativistic one-dimensional problems First, we consider special relativistic problems in one spatial dimension. With this setup, the tetrad basis in Sec. III B is reduced to a coordinate vector in Minkowski spacetime. Thus, the setup is suitable for validating the TABLE I. Initial conditions used for special relativistic one-dimensional test problems. The third column shows the Γ index and the second-to-last column shows the final time of the simulations, t. Riemann solvers described in Sec. III C and III D. 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 I. We note that the test suite employed in this paper is the same as that presented in Refs. [2,87]. 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 3rdorder 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 ≈ 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.
The right panel of Fig. 3 shows the rest-mass density profile for the second test problem (Problem HD2) listed in Table I. 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 ≈ −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 I) is the often-employed shock-tube problem. Here, the simulation domain spans x ∈ [−0.5, 0.5] with ∆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 ≈ 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. [1]. 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 onedimensional problems.
For the fourth (final) hydrodynamics test problem (Problem HD4), we employ a simulation domain of x ∈ [−0.5, 0.5] with a grid spacing of ∆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 ≈ 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. [82]). Otherwise, spurious waves appear irrespective of which solver is used (not shown).

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 ∈ [−0.5, 0.5] with grid spacing ∆x = 0.005 (i.e., N x = 100). For Problem MHD6, the domain is identical but we employ a higher resolution with ∆x = 0.0025 and N x = 200.
In the first problem (Problem MHD1 in Table I)  lution 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 I) 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. (3.129)). 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 [88]. In this problem, the solution consists of a left-propagating rarefaction wave, a right-propagating slow shock wave (located at x ≈ 0.18 in Fig. 6), and a right-propagating rarefaction wave. In addition there is a right-ward propagating contact discontinuity located at x ≈ 0.15 in Fig. 6 adjacent to the shock wave. Finally at x ≈ 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 1storder 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  Fig. 7 for the solutions (the inset in the B y (B z ) panel shows a closeup 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 recon- struction, 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 ≈ −0.3 (+0.3), and left (right)propagating slow waves located at x ≈ −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 ≈ 0. This is known as the wallheating problem [89]: it is well-known that Godunov-type schemes cannot avoid this pathological behavior. As reported in Ref. [2], 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 onedimensional suite, the solution consists of all seven waves [58]. The numerical results are shown in Fig. 9. In this problem, a contact discontinuity appears at x ≈ 0.05, a rarefaction wave propagates to the left of the contact discontinuity, which can be seen at x ≈ −0.4, and the rotational discontinuity at x ≈ −0.06 and the slow shock at x ≈ −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 ≈ 0.4. The rotational discontinuity at x ≈ 0.08 and the slow shock at x ≈ 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 ρ in Fig. 9). With 1st-order reconstruction, however, it is difficult to disentangle the left/rightpropagating rotational discontinuities and left/rightpropagating 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 rightpropagating rotational discontinuity and the left/rightpropagating slow shock are captured, but the leftpropagating rotational discontinuity is not. This demonstrates the ability of the HLLD-CT GS solver to properly capture all seven of the emergent waves.
(4.3) Figure 10 show the logarithmic contour of the restmass 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 higherorder reconstruction.

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 ∈ [−2, 2] and y ∈ [−3, 3], and set ∆x = ∆y = 0.02, i.e., (N x , N y ) = (100,150) in the xand y-directions, respectively. Periodic boundary conditions are imposed at the x-and y-boundaries. We set the adiabatic index to Γ = 5/3, employ 3rd-order PPM reconstruction, and set the CFL number to 0.45. The initial condition is given by ρ = 1, P = 2.5 for x 2 + y 2 < 0.5, 0.1 for x 2 + y 2 > 0.5. (4.4) 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-Taylorlike instability (known in this context as the Richtmyer-Meshkov instability) has developed [83]. 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.

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. [75,83]. The initial profile for the magnetic field is given by where B 0 = 1 and a = 0.04. The density is uniform with ρ = 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 where β 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 zcomponent of the vector potential as where k x = 2π/L x , k y = 2π/L y , = 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 ∈ [−0.5, 0.5] and y ∈ [−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 ≈ 0.557 and the Alfvén timescale is t A ≈ 1.78. 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 ≈ 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. [4]). 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 magneticfield 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 ≈ 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.

Magnetohydrodynamics: two-dimensional Kelvin-Helmoltz instability
The second two-dimensional problem in special relativistic magnetohydrodynamics is the Kelvin-Helmholtz instability, as proposed in Refs. [2,91]. For this, we prepare a simulation domain which spans x ∈ [−0.5, 0.5] and y ∈ [−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, v x = −v sh tanh(y/a), (4.8) where v 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 ρ = 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. [87], 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 sh = 0.5) which is twice that used in our runs.
The magnetic field at t = 0 is given by i.e. the magnetic field is initially uniform and parallel to the velocity in the lower-half of the xy-plane. We set σ pol = 0.01. The shear layer is perturbed by the motion in the y-direction as v y = 1 40000 sin(2πx) exp −100y 2 , (4.10) while v z = 0. Figure 14 shows the perturbed velocity difference ∆v y ≡ (v y max − v y min )/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 ∼ 10 −5 . We find exponential growth followed by nonlinear saturation at the end of the linear phase at t ∼ 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. [87] (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 ∼ 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. [91] 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.
C. General relativistic problems in a fixed background spacetime

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 [92]. The Bondi flow in Schwarzschild coordinates has been extensively discussed in the literature (see, e.g., Ref. [83]). Following previous work [82,84], we adopt the parameters for this problem as follows: an adiabatic index of Γ = 4/3, an adiabat of K = 1, and a critical radius of r c = 8M , where M denotes the black hole mass. With this setup, the mass accretion rateṀ 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 [93]. Thus, for preparing a practical setup in this test problem, a nonrotating 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 [94,95]. Note that in both of these latter two coordinate systems, the fluid fourvelocity exhibits pathological behavior near the horizon [95] [96]. 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. (3.36)) 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 B.
We employ a simulation domain in Cartesian coordinates spanning x, y, z ∈ [0, L] with L = 10M . The grid spacing of the simulation is given by ∆ = ∆x = ∆y = ∆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 ∆ = 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 in = 0.4M . Note that the horizon in this geometry is located at r BH ≈ 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 [92]. 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. [63,82,84,97].
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 ≈ 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 B) the frame-dragging effect of the cell interface is taken into account with a better accuracy (see also Fig. 2) [98].

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 [76]. Therefore, we employ the same flow profile used in the previous section. From the divergence-free condition (2.15), the radial magnetic field should be B R ∝ 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 /ρ = 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 ∈ [0, L] with L = 12.8M . The grid spacing of the simulation is ∆ = ∆x = ∆y = ∆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., ∆ = 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 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 [92] (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 bottomright 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 ≈ 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. (2.9)-(2.13)).
A. Hydrodynamics: binary neutron star merger

Setup
First, we consider non-magnetized asymmetric binary neutron stars with masses of 1.2 and 1.5M . We utilize the spectral method library LORENE [99-103] 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 [104]. The initial orbital angular velocity is set to be m 0 Ω 0 = 0.028 where m 0 = 2.7M is the total mass of the binary.
Our solver for Einstein's equations implements the BSSN-puncture formulation [105][106][107][108], locally incorporating the Z4c prescription for constraint propagation [109]. 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 3rdorder PPM cell reconstruction.
We employ the SFHo equation of state for relatively high-density nuclear matter [110], and the Timmes (Helmholtz) equation of state for the low-density part [111]. 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 ρ atm = 10 3 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 fin = 38.7 km (see below for the FMR setup in detail). We also set a power-law profile of the atmospheric density of ρ atm = 10 3 (L fin /r) 3 g/cm 3 for r > L 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 ≈ 0.17 g/cm 3 and if ρ 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 MeV.
We also explicitly solve the radiation-hydrodynamics equations for neutrinos in time using an approximate neutrino-transfer scheme based on a leakage scheme [77] and the truncated moment formalism [79,112]. The cooling source terms are computed using a general-relativistic leakage scheme [74], and heating source terms due to neutrino capture processes are computed by the method presented in Ref. [78].
The computational region consists of 13 levels of FMR half-cubic domains. The size of each FMR domain is ∈ [−L/2 l−1 , L/2 l−1 ] for x and y, and z ∈ [0, L/2 l−1 ] with l = 1, 2, · · · , 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 ≈ 158, 000 km and N = N x = N y = N z = 258. Thus the grid spacing of the finest FMR domain is ∆x 13 = ∆y 13 = ∆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 ∆x 13 = 200 m and ∆x 13 = 250 m, respectively. For the HLLC solver, we perform an additional simulation with N = 377 and ∆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.

Inspiral phase
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. [82] 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. ∆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 [61,62,72], but could be improved if we were to employ a more accurate reconstruction scheme such as MP5 [113]. 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, ∆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.

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 viola-tion 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   Fig. 20 shows a close-up of the results with the HLLC and HLLE solvers for ∆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 20 ms, the oscillation ampli-tude 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 densityweighted 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 [115]. We find spurious spin-down of the black hole due to numeri-  We measure the spin-down rate in the HLLC run and estimate that the dimensionless spin decreases by 0.1 in 1 s if r AH /∆x 13 15 where r 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 of the spin-down, which in turn will result in spurious mass accretion. Thus the grid resolution must be chosen carefully when the main aim is to quantitatively explore the evolution of the disk and subsequent mass ejection.
The right panel of Fig. 21 shows the gravitationally bound baryonic mass outside the apparent horizon [116]. Irrespective of which Riemann solver we employ, the bounded baryonic mass is not a monotonic function of the grid spacing. Before the formation of the black hole, the non-axisymmetric density structure of the remnant massive neutron star exerts a gravitational torque on the fluid elements. As a result, angular momentum is transported outwards. Thus the longer lifetime of the remnant massive neutron star results in the formation of a more massive torus after the neutron star remnant collapses to the black hole. Because the lifetime of the remnant massive neutron star is not a monotonic function of the grid spacing it is a natural consequence that we find that the baryonic mass of bound material does not converge as the resolution is increased. Nevertheless, in the simulations with the HLLC solver, the gravitationally bound baryonic mass is found to lie in a narrow range, between 0.055M and 0.075M , at the time of formation of the black hole (at t ∼ 30 ms). When we employ the HLLE solver, the bound baryonic mass is systematically larger than that with the HLLC solver (between 0.100M and 0.125M ). This is because the lifetime of the remnant massive neutron star is systematically longer in the simulations with the HLLE solver than with the HLLC solver, as already mentioned. Therefore, when one employs the HLLE solver, one should keep in mind that the bound baryonic mass could be overestimated with a systematic error of O(10 −2 M ). Figure 22 shows the time-evolution of the luminosity of electron neutrinos (left panel) and of electron antineutrinos (right panel). These plots show that the luminosity increases quickly after merger, reaching a peak value of ≈ 1.2 × 10 53 erg/s for the electron neutrinos and ≈ 1.9 × 10 53 erg/s for the electron antineutrinos at t ≈ 20 ms. These values agree broadly with our previous results [117]. After the formation of the black hole, the luminosity quickly decreases because the high density and temperature regions of the remnant massive neutron star are swallowed into the black hole [117,118]. Note that the overall evolution of the neutrino luminosity in the remnant massive neutron star phase does not significantly depend either on the Riemann solver, nor on the spatial grid spacing.
Finally, Fig. 23 shows the time-evolution of the gravitationally unbound baryonic mass, i.e. the ejecta mass. In this model (i.e., the model with appreciable mass asymmetry in the binary), mass ejection is driven primarily by the tidal force from the heavier component to the lighter one. The blue and green curves denote results with the HLLC and HLLE solvers, respectively. The solid, dashed, dotted, and dot-dashed curves denote the results with grid spacings of ∆x 13 = 100 m, 150 m, 200 m, and 250 m, respectively. The inset depicts the ejecta-mass evolution on a logarithmic scale along the vertical axis, and the shaded region denotes the violation of baryonic mass conservation. We find that the spurious mass ejection during the inspiral phase is O(10 −7 M ), and it decreases as the resolution is enhanced. We also find that the error in baryonic mass conservation is below 10 −7 M even after the merger. This figure shows that the ejecta mass decreases as the grid spacing is improved from 250 m to 150 m. This is likely to be related to the spurious expansion of the less massive neutron star during the inspiral phase, which we discussed above. This spurious expansion is enhanced in the lower resolution runs. When we employ ∆x 13 = 100 m for the HLLC solver, the ejecta mass is approximately identical to that with ∆x 13 = 150 m. Therefore, the convergence for the ejecta mass is approximately achieved in this model. Figure 23 also shows that the amount of ejecta mass in the simulation with the HLLC solver is smaller than that with the HLLE solver for a given grid spacing. Quantitatively, the ejecta mass difference due to the Riemann solver is ≈ 10 −3 M for ∆x 13 = 150 m in this model. This difference arises from how accurately the employed Riemann solver can capture the neutron-star shape during the late inspiral phase. As we have already emphasized, with the HLLE solver the neutron star spuriously expands during the inspiral phase. As a result, the less massive neutron star is more subject to (partly artificial) tidal deformation, thereby ultimately increasing the tidallydriven ejecta mass. When we employ the HLLC solver together with a high grid resolution this artifact is mitigated. This is one of the advantages of using a more sophisticated Riemann solver in this problem.
We conclude that the HLLC solver is superior to the HLLE solver both during the inspiral and post-merger phases of the binary neutron star merger. In particular, we note that for the purpose of obtaining accurate and high-precision gravitational waveforms during the late inspiral phase over more than 10 orbits, the HLLE solver is likely not an appropriate choice.
B. Magnetohydrodynamics: binary neutron star merger (evolution of remnant)

Setup
As an application of the new Riemann solvers (paired with our new implementation of the constrained transport scheme) to relativistic magnetohydrodynamics in a dynamical spacetime, we consider the evolution of the magnetized torus surrounding a black hole formed after a binary neutron star merger. The initial condition is taken from the final moment of the hydrodynamics simulation for a binary neutron star merger presented in the previous section. Specifically, our initial condition is taken from the result of the simulation run with the HLLC solver at a resolution of ∆x 13 = 150 m and at t ≈ 76 ms. The grid setup is exactly the same as in the hydrodynamics simulation for the binary neutron star merger.
We initialize the magnetic field inside the torus with a vector potential of the form where i = x or y, x BH and y BH denote the x-and ycoordinates of the central black hole, P is the gas pressure, and P max is its maximum. We choose the amplitude A b such that the initial maximum magnetic field strength is 10 15 G. We employ the HLLD-CT GS and HLLE-CT HLLE solvers and compare the results. We also employ Balsara's method to ensure the divergence-free condition and magnetic flux conservation in the refinement boundary [73,119,120]. In our implementation, not only is the divergence-free condition preserved to machine-precision, but the magnetic flux is also preserved across the refinement boundary. We note that the vector potential method [121], which has been widely implemented in numerical relativity codes, does not ensure the latter property [62][63][64]66].

Post-merger evolution
The left and right panels of Fig. 24 show the time evolution of the electromagnetic energy and the time evolution of the magnetorotational-instability quality factor, respectively. The electromagnetic energy is defined by [122] The origin of the time-axis is the same as in Fig. 18. The blue and cyan curves denote results with the HLLD CT-GS and HLLE CT-HLLE solvers, respectively. The magnetorotational-instability quality factor is defined by where is the wavelength of the fastest growing mode of the axisymmetric magnetorotational instability [123,124]. Note that we introduce a cut-off density in the quality factor to determine in which part of the torus the magnetorotational instability is resolved. These panels show that the electromagnetic energy is amplified during the initial stage of t 84-85 ms primarily due to magnetic winding rather than the magnetorotational instability, because the fastest growing mode of the magnetorotational instability in the high-density regions of the torus is not well resolved at these early times (see the solid curves in the right panel with a cut-off density of ρ cut = 10 10 g cm −3 ). During this stage, the electromagnetic energy with the HLLD-CT GS solver is larger than that with the HLLE-CT HLLE solver because the large numerical diffusion inherent in the HLLE-CT HLLE solver results in the diffusion of magnetic field lines. In addition, in the orbital plane magnetic fields are forced to reconnect because we impose plane symmetry with respect to the equatorial plane. For the HLLE-CT HLLE solver, reconnection in this plane is also enhanced due to numerical diffusion, and thus reduces the electromagnetic energy even further (see also the magnetized current sheet problem in Fig. 12).
After t ≈ 84-85 ms by which the poloidal magneticfield strength has been enhanced nearly to saturation level due to winding and subsequent outgoing motion resulting from the enhanced magnetic-field pressure, magnetorotational instability-driven turbulence begins to develop in the high-density region of the torus, because the fastest growing mode is now resolved by more than ten grid points (see the right panel of Fig. 24). This then establishes a turbulent state. Figure 25 displays the magnetic-field structure in the x-z plane. This figure shows that by the time the magnetic-field strength has saturated, turbulence has developed and an outflow associated with the turbulent activity is driven from the disk. The middle panels of Fig. 25 show the magnetic-field structure at t ≈ 90 ms. With the HLLD-CT GS solver, the inside of the torus exhibits smaller-scale turbulence than that with the HLLE-CT HLLE solver (see, e.g., the region of x ∈ [20,50] km and z ∈ [0, 20] km). The larger structures seen in the colormap also suggest that magnetic-field lines are more coherent with the HLLE-CT HLLE solver than they are with the HLLD-CT GS solver. Our explanation for this is that with the HLLE-CT HLLE solver, the magnetorotational instability is less resolved, and thus, the small-scale turbulent structure is less developed. As a result, largescale magnetic fields appear to be spuriously enhanced with HLLE-CT HLLE compared to HLLD-CT GS.
As evidence for this explanation, we calculate the power spectrum density of the electromagnetic energy defined by andb * (k i ) is its complex conjugate. Here, k i is the wave vector with i = x, y, z and k 2 = i k 2 i . dΩ k is a solid angle in k-space. We employ the Python package fiNUFFT [125,126] Figure 26 plots the power spectrum density of the magnetic-field energy at t ≈ 90.2 ms. The blue and cyan curves denote solutions with the HLLD-CT GS and HLLE-CT HLLE solvers, respectively. With the help of the non-uniform Fast Fourier Transformation, we obtain a power spectrum density that spans three orders of magnitude. It clearly shows that the power spectrum amplitude around k/(2π) = 10 −6 cm −1 is larger in the HLLE-CT HLLE than in the HLLD-CT GS run. This implies that a relatively large-scale magnetic field with a scale of ≈ 10 6 cm is generated in the HLLE-CT HLLE run compared to the HLLD-CT GS run. On the other hand, at small scales (i.e. with k/(2π) 10 −5 cm −1 ), the power spectrum density is higher in the HLLD-CT GS run than in the HLLE-CT HLLE run. This shows the HLLD-CT GS solver is able to sustain smaller-scale magnetorotational instability-driven turbulence than the HLLE-CT HLLE solver. Figure 24 indicates that the electromagnetic energy is still increasing for t 90 ms. We find that (i) the growth is not exponential, and (ii) the growth rate with the HLLE-CT HLLE solver is higher than with the HLLD-CT GS solver. This indicates that magnetic winding of a coherent poloidal magnetic field proceeds more efficiently (though spuriously) in the simulation with the HLLE-CT HLLE solver than with the HLLD-CT GS solver. This in turn enhances the launch of a magnetic tower outflow in the polar direction, as shown in the bottom panels of Fig. 25. While this outflow is observed regardless of which solver is used, we observe a more powerful magnetic tower outflow with the HLLE-CT HLLE solver, which reflects the greater (but spurious) coherency of the magnetic-field lines when we use of this solver.
To quantify how powerful the magnetic tower outflow is, we plot the angular distribution of the Poynting flux − √ −g(T r t ) (EM) = − √ −g(b 2 u r u t − b r b t ) on a sphere of r ≈ 50 km in Fig. 27. The snapshot is taken at t ≈ 103 ms. In the polar region, the Poynting flux with the HLLE-CT HLLE solver is much stronger than with the HLLD-CT GS solver. This plot suggests that the power of the magnetic tower outflow is overestimated when we employ the HLLE-CT HLLE solver.

VI. SUMMARY AND CONCLUSION
We implemented the advanced Riemann solvers HLLC [1] and HLLD [2] in our numerical relativity neutrino-radiation magnetohydrodynamics code. We validated our implementation by performing one-and multidimensional test problems in both Minkowski spacetime and in a fixed background spacetime, both in relativistic hydrodynamics and relativistic (ideal) magnetohydrodynamics. In the relativistic hydrodynamics test problems, we found that the HLLC solver is always superior to the HLLE solver, in particular, for the multi-dimensional case: the spurious waves associated with the HLLE solver disappear, and the grid resolution is effectively improved, when we employ the HLLC solver. For relativistic magnetohydrodynamics test problems, we also found that the performance of the HLLD solver together with the constrained transport method proposed by Gardiner and Stone [3], which relies on the accuracy of a Riemann solver, is the best for both one-dimensional as well as multi-dimensional test problems.
We also performed simulations of a non-magnetized asymmetric binary neutron star merger in a dynamical spacetime with the HLLC and HLLE solvers. We found that spurious broadening of the neutron star surface during the inspiral phase can be mitigated by employing the HLLC solver. As a result, the less massive companion of the binary is less subject to tidal elongation during the late inspiral phase than when the HLLE solver is used. This point is particularly important for deriving a highprecision gravitational waveform during the late inspiral phase, because one has to compute the orbital evolution precisely, i.e. excluding spurious numerical effects for this problem. The solution with the HLLC solver also differs from that with the HLLE solver in the subsequent postmerger evolution. For example, the amount of dynamical ejecta driven by the tidal interaction of the two stars and the lifetime of the remnant massive neutron star are overestimated when we employ the HLLE solver.
The neutron-rich dynamical ejecta and post-merger ejecta, the latter of which is launched from the merger remnant by an effective turbulent viscosity due to the magnetorotational instability [49-51, 53, 54, 78], will shine by means of radioactive decay of r-process elements which have been freshly synthesized in the ejecta (see, e.g., [18,20,127]). One of the most important aims in the observation of binary neutron star mergers is to observe this signal and to infer the binary parameters by comparing the observational results with the theoretical prediction from numerical relativity simulations. Therefore, we conclude that employing a better solver (i.e., the HLLC solver rather than the HLLE solver) is crucial for reliable modeling of electromagnetic counterparts from binary neutron star mergers.
We also performed simulations of the binary neutron star merger remnant, i.e. a black hole surrounded by a massive torus, in the framework of neutrino-radiation magnetohydrodynamics. We embedded a purely poloidal magnetic-field loop inside the torus and performed simulations with the HLLD-CT GS and HLLE-CT HLLE solvers. We found that (i) artificial magnetic-field dissipation is suppressed, and (ii) a well-resolved magneto-turbulent state is reproduced, when we employed the HLLD-CT GS solver. On the other hand, when we employed the dissipative HLLE-CT HLLE solver, the coherency of the magnetic-field lines is artificially enhanced, resulting in the launch of a powerful magnetic tower outflow due to magnetic winding of this coherent poloidal field. The emergence of a Poynting flux-dominated outflow from the black hole-torus system could be a key ingredient for driving a short gamma-ray burst from the compact binary merger remnant [55]. Therefore, we conclude that employing the HLLD solver paired with the constrained transport method proposed by Gardiner and Stone [3] is crucial for reliable modeling of the central engine of short gamma-ray bursts.
As a future project, we plan to perform long-term simulations of binary neutron star mergers and black holeneutron star binary mergers, employing the advanced Riemann solvers which we have implemented in our code.

ACKNOWLEDGMENTS
We thank Tsz Lok Lam, Sho Fujibayashi, Shinya Wanajo, and the members of the Computational Relativistic Astrophysics division at the AEI for useful discussions. Kenta Kiuchi thanks Koutarou Kyutoku for providing the initial data for the binary neutron star simulations, and for checking the manuscript. Kenta Kiuchi also thanks Kota Hayashi for providing the script used to generate Figure 27. Numerical simulations were performed on the Sakura, Cobra, and Raven clusters at the Max Planck Computing and Data Facility and on the Cray XC50 at CfCA of the National Astronomical Observatory of Japan. This work was in part supported by the Grantin-Aid for Scientific Research (grant Nos. 18H01213, 19K14720, and 20H00158) of Japan MEXT/JSPS. For convenience, we explicitly show the tetrad basis for the Riemann problem in the y-and z-directions.
e (ŷ)µ =B (β y , δ i y ) , e (ẑ)µ =D β z γ xx − β x γ xz , 0, γ xx γ yz − γ xy γ xz , γ xx γ zz − γ 2 xz . (A11) Bondi flow is usually described in Schwarzschild coordinates. However, our numerical relativity code has a high affinity with puncture coordinates because the solver for Einstein's equations handles a black hole with the moving puncture gauge. In numerical relativity with this gauge condition, black holes relax to a stationary solution in the so-called limit hypersurface. This implies that the code test should be done employing this special stationary hypersurface. To do this, one needs to seek a coordinate transformation from the Schwarzschild coordinates to the puncture coordinates (i.e., the coordinates of the limit hypersurface). One simple way of doing this is to describe a black hole as the maximal trumpet black hole puncture solution described in Ref. [94]. In these coordinates, the fluid quantities are well-behaved on the horizon.

Maximal trumpet black hole puncture
The stationary solution of the Schwarzschild spacetime in the limiting hypersurface can be written by where Here, C is the integration constant and R is the circumferential radius. A number of numerical relativity simulations of a single black hole spacetime using the moving puncture gauge showed that the numerical solution settles down to a member of the family with C = 3 √ 3M 2 4 , which has a limiting surface at R = 3M/2 [93]. If we consider a transformation of this solution into the isotropic coordinates by identifying the spatial metric in both coordinates as f −2 dR 2 + R 2 dΩ 2 = ψ 4 (dr 2 + r 2 dΩ 2 ), one may find a solution for r and ψ as [94] r = 2R + M + (4R 2 + 4M R + 3M 2 ) 1/2 4 × (4 + 3 √ 2)(2R − 3M ) 8R + 6M + 3(8R 2 + 8M R + 6M 2 ) 1/2 where we assumed C = 3 √ 3M 2 /4. The lapse function, shift vector, and non-zero components of the extrinsic curvature are given by

Velocity field and magnetic field of Bondi flow
The velocity field of Bondi flow in Schwarzschild coordinates should be transformed into the isotropic coordinates described in the previous section. The radial component is obtained by where u R is the radial velocity of Bondi flow in Schwarzschild coordinates (see, e.g., Ref. [83]). The time component is obtained by the normalization of the four velocity: Note that the four velocity in these coordinates does not exhibit pathological behavior on the horizon, which can be confirmed by a Taylor expansion of Eq. (B13) near the horizon [95]. Note also that the lower components of the four velocity, u t and u r , are well-behaved at the horizon because the metric has a regular form in the maximal trumpet geometry (B5). For magnetized Bondi flow, the radial component of the magnetic field in the maximal trumpet geometry is given by where B R is the radial component of the magnetic field in Schwarzschild coordinates. In the case of a purely radial magnetic field, the divergence-free condition (2.15) requires the radial component of the magnetic field in Schwarzschild coordinates be