Entropy Symmetrization and High-Order Accurate Entropy Stable Numerical Schemes for Relativistic MHD Equations

TL;DR

This paper develops high-order accurate, entropy stable numerical schemes for the relativistic MHD equations by symmetrizing the system and carefully designing fluxes and dissipation to ensure stability and accuracy.

Contribution

It introduces a symmetrizable RMHD system with a convex entropy pair and constructs high-order entropy stable schemes based on this formulation.

Findings

Schemes demonstrate high accuracy in numerical tests.

The methods are robust for complex RMHD simulations.

Entropy stability is successfully achieved with high-order discretizations.

Abstract

This paper presents entropy symmetrization and high-order accurate entropy stable schemes for the relativistic magnetohydrodynamic (RMHD) equations. It is shown that the conservative RMHD equations are not symmetrizable and do not admit a thermodynamic entropy pair. To address this issue, a symmetrizable RMHD system, equipped with a convex thermodynamic entropy pair, is proposed by adding a source term into the equations, providing an analogue to the non-relativistic Godunov--Powell system. Arbitrarily high-order accurate entropy stable finite difference schemes are developed on Cartesian meshes based on the symmetrizable RMHD system. The crucial ingredients of these schemes include (i) affordable explicit entropy conservative fluxes which are technically derived through carefully selected parameter variables, (ii) a special high-order discretization of the source term in the symmetrizable RMHD system, and (iii) suitable high-order dissipative operators based on essentially non-oscillatory reconstruction to ensure the entropy stability. Several numerical tests demonstrate the accuracy and robustness of the proposed entropy stable schemes.

Figures

Tables

Table 1. Table 1: Example 6 : l 1 superscript 𝑙 1 l^{1} -errors and l 2 superscript 𝑙 2 l^{2} -errors in v 2 subscript 𝑣 2 v_{2} at t = 0.5 𝑡 0.5 t=0.5 , and corresponding convergence rates for the EC4 and ES4 schemes at different grid resolutions.

$N$	EC4				ES4
	$l^1$-error	order	$l^2$-error	order	$l^1$-error	order	$l^2$-error	order
8	3.14e-3	–	3.53e-3	–	4.23e-3	–	4.62e-3	–
16	2.10e-4	3.91	2.33e-4	3.92	2.34e-4	4.17	2.60e-4	4.15
32	1.33e-5	3.98	1.48e-5	3.98	1.38e-5	4.09	1.54e-5	4.08
64	8.34e-7	3.99	9.27e-7	3.99	8.46e-7	4.03	9.40e-7	4.03
128	5.22e-8	4.00	5.80e-8	4.00	5.25e-8	4.01	5.83e-8	4.01
256	3.26e-9	4.00	3.62e-9	4.00	3.27e-9	4.00	3.63e-9	4.00
512	2.04e-10	4.00	2.26e-10	4.00	2.04e-10	4.00	2.27e-10	4.00

Table 2. Table 2: Same as Table 1 except for the EC6 and ES5 schemes.

$N$	EC6				ES5
	$l^1$-error	order	$l^2$-error	order	$l^1$-error	order	$l^2$-error	order
8	3.85e-4	–	4.36e-4	–	5.64e-3	–	6.16e-3	–
16	7.11e-6	5.76	7.96e-6	5.78	2.62e-4	4.43	2.94e-4	4.39
32	1.18e-7	5.91	1.31e-7	5.93	8.59e-6	4.93	9.53e-6	4.95
64	1.86e-9	5.98	2.07e-9	5.98	2.68e-7	5.00	2.97e-7	5.00
128	2.92e-11	6.00	3.24e-11	5.99	8.35e-9	5.00	9.27e-9	5.00
256	4.81e-13	5.92	5.34e-13	5.92	2.61e-10	5.00	2.89e-10	5.00

Table 3. Table 3: Initial data of the four 1D RPs in Example 1 .

		$\rho$	$v_1$	$v_2$	$v_3$	$B_1$	$B_2$	$B_3$	$p$
RP I	left state	1.08	0.4	0.3	0.2	2	0.3	0.3	0.95
	right state	1	-0.45	-0.2	0.2	2.0	-0.7	0.5	1
RP II	left state	1	0	0	0	5	6	6	30
	right state	1	0	0	0	5	0.7	0.7	1
RP III	left state	1	0	0.3	0.4	1	6	2	5
	right state	0.9	0	0	0	1	5	2	5.3
RP IV	left state	1	0	0	0	0	2	0	10
	right state	0.1	0	0	0	0	0	0	1

Equations

$$\frac{\partial \mathbf{U}}{\partial t}+\sum _{i=1}^d\frac{\partial {\mathbf{F}}_i(\mathbf{U})}{\partial x_i}=0,$$

$$\nabla \cdot \mathbf{B}:=\sum _{i=1}^d\frac{\partial B_i}{\partial x_i}=0,$$

$${\bf F}_{i}({\bf U})=\left(Dv_{i},v_{i}{\bf m}^{\top}-B_{i}\big{(}\gamma^{-2}{\bf B}^{\top}+({\bf v}\cdot{\bf B}){\bf v}^{\top}\big{)}+p_{tot}{\bf e}_{i}^{\top},v_{i}{\bf B}^{\top}-B_{i}{\bf v}^{\top},m_{i}\right)^{\top},$$

$${\mathcal{Q}}_i^{\prime }(\mathbf{U})={\mathcal{E}}^{\prime }(\mathbf{U}){\mathbf{F}}_i^{\prime }(\mathbf{U}),i=1,\dots ,d,$$

$$0={\mathcal{E}}^{\prime }(\mathbf{U})\left(\frac{\partial \mathbf{U}}{\partial t}+\sum _{i=1}^d\frac{\partial {\mathbf{F}}_i(\mathbf{U})}{\partial x_i}\right)={\mathcal{E}}^{\prime }(\mathbf{U})\frac{\partial \mathbf{U}}{\partial t}+\sum _{i=1}^d{\mathcal{Q}}_i^{\prime }(\mathbf{U})\frac{\partial \mathbf{U}}{\partial x_i}=\frac{\partial \mathcal{E}}{\partial t}+\sum _{i=1}^d\frac{\partial {\mathcal{Q}}_i}{\partial x_i}.$$

$$\frac{\partial \mathcal{E}}{\partial t}+\sum _{i=1}^d\frac{\partial {\mathcal{Q}}_i}{\partial x_i}\le 0,$$

$$\frac{\partial \mathbf{U}}{\partial \mathbf{W}}\frac{\partial \mathbf{W}}{\partial t}+\sum _{i=1}^d\frac{\partial {\mathbf{F}}_i}{\partial \mathbf{W}}\frac{\partial \mathbf{W}}{\partial x_i}=0,$$

$$\frac{\partial (\rho \gamma S)}{\partial t}+\nabla \cdot (\rho \gamma S\mathbf{v})+(\Gamma -1)\frac{\rho \gamma (\mathbf{v}\cdot \mathbf{B})}{p}\nabla \cdot \mathbf{B}=0.$$

$$\mathcal{E}(\mathbf{U})=-\frac{\rho \gamma S}{\Gamma -1},{\mathcal{Q}}_i(\mathbf{U})=-\frac{\rho \gamma S{v}_i}{\Gamma -1}$$

$${\bf W}={\mathcal{E}}^{\prime}({\bf U})^{\top}=\left(\frac{\Gamma-S}{\Gamma-1}+\frac{\rho}{p},~{}\frac{\rho\gamma}{p}{\bf v}^{\top},~{}\frac{\rho\gamma}{p}\Big{(}(1-|{\bf v}|^{2}){\bf B}^{\top}+({\bf v}\cdot{\bf B}){\bf v}^{\top}\Big{)},~{}-\frac{\rho\gamma}{p}\right)^{\top}.$$

$$\mathbf{V}={\left(\rho ,{\mathbf{v}}^{\top },{\mathbf{B}}^{\top },p\right)}^{\top }.$$

$${\mathcal{E}}^{\prime }(\mathbf{V})=\frac{1}{\Gamma -1}\left(\gamma (\Gamma -S),-\rho S{\gamma }^3{\mathbf{v}}^{\top },0,0,0,-\frac{\rho \gamma }{p}\right),$$

$$\frac{\partial \mathbf{U}}{\partial \mathbf{V}}=\left(\begin{array}{cccc}\gamma & \rho {\gamma }^3{\mathbf{v}}^{\top }& 0_3^{\top }& 0\\ {\gamma }^2\mathbf{v}& (\rho h{\gamma }^2+|\mathbf{B}{|}^2){\mathbf{I}}_3-\mathbf{B}{\mathbf{B}}^{\top }+2\rho h{\gamma }^4\mathbf{v}{\mathbf{v}}^{\top }& -(\mathbf{v}\cdot \mathbf{B}){\mathbf{I}}_3+2\mathbf{v}{\mathbf{B}}^{\top }-\mathbf{B}{\mathbf{v}}^{\top }& \frac{\Gamma {\gamma }^2}{\Gamma -1}\mathbf{v}\\ 0_3& {\mathbf{O}}_3& {\mathbf{I}}_3& 0_3\\ {\gamma }^2& (2\rho h{\gamma }^4+|\mathbf{B}{|}^2){\mathbf{v}}^{\top }-(\mathbf{v}\cdot \mathbf{B}){\mathbf{B}}^{\top }& (1+|\mathbf{v}{|}^2){\mathbf{B}}^{\top }-(\mathbf{v}\cdot \mathbf{B}){\mathbf{v}}^{\top }& \frac{\Gamma {\gamma }^2}{\Gamma -1}-1\end{array}\right),$$

$${\mathcal{Q}}_i^{\prime }(\mathbf{U})={\mathcal{E}}^{\prime }(\mathbf{U}){\mathbf{F}}_i^{\prime }(\mathbf{U})+\frac{\rho \gamma }{p}(\mathbf{v}\cdot \mathbf{B})B_i^{\prime }(\mathbf{U}),i=1,\dots ,d,$$

$$\frac{\partial {\mathbf{F}}_1}{\partial \mathbf{V}}=v_1\frac{\partial \mathbf{U}}{\partial \mathbf{V}}+\left(\begin{array}{cccc}0& \rho \gamma {\mathbf{e}}_1^{\top }& 0_3^{\top }& 0\\ 0_3& (\rho h{\gamma }^2+|\mathbf{B}{|}^2)\mathbf{v}{\mathbf{e}}_1^{\top }+{\mathbf{M}}_1& {\mathbf{M}}_2& {\mathbf{e}}_1\\ 0_3& \mathbf{B}{\mathbf{e}}_1^{\top }-B_1{\mathbf{I}}_3& -\mathbf{v}{\mathbf{e}}_1^{\top }& 0_3\\ 0& (\rho h{\gamma }^2+|\mathbf{B}{|}^2){\mathbf{e}}_1^{\top }-B_1{\mathbf{B}}^{\top }+v_1{\left((\mathbf{v}\cdot \mathbf{B})\mathbf{B}-|\mathbf{B}{|}^2\mathbf{v}\right)}^{\top }& {\mathbf{\beta }}_1^{\top }-(\mathbf{v}\cdot \mathbf{B}){\mathbf{e}}_1^{\top }& v_1\end{array}\right),$$

$$\displaystyle{\bf M}_{1}={\bf e}_{1}\big{(}({\bf v}\cdot{\bf B}){\bf B}^{\top}-|{\bf B}|^{2}{\bf v}^{\top}\big{)}+B_{1}(2{\bf B}{\bf v}^{\top}-{\bf v}{\bf B}^{\top})-({\bf v}\cdot{\bf B})(B_{1}{\bf I}_{3}+{\bf B}{\bf e}_{1}^{\top})$$

$${\mathbf{M}}_2=(1-|\mathbf{v}{|}^2)({\mathbf{e}}_1{\mathbf{B}}^{\top }-\mathbf{B}{\mathbf{e}}_1^{\top })+(\mathbf{v}\cdot \mathbf{B})({\mathbf{e}}_1{\mathbf{v}}^{\top }-\mathbf{v}{\mathbf{e}}_1^{\top })-B_1\mathbf{v}{\mathbf{v}}^{\top }-B_1(1-|\mathbf{v}{|}^2){\mathbf{I}}_3.$$

$$\frac{\partial \mathbf{W}}{\partial \mathbf{V}}=\left(\begin{array}{cccc}h/p& 0_3^{\top }& 0_3^{\top }& -\frac{\rho }{p^2}-\frac{1}{p(\Gamma -1)}\\ \frac{\gamma }{p}\mathbf{v}& \frac{\rho {\gamma }^3}{p}{\mathbf{M}}_3& {\mathbf{O}}_3& -\frac{\rho \gamma }{p^2}\mathbf{v}\\ \frac{\gamma }{p}((1-|\mathbf{v}{|}^2)\mathbf{B}+(\mathbf{v}\cdot \mathbf{B})\mathbf{v})& \frac{\rho {\gamma }^3}{p}{\mathbf{M}}_4& \frac{\rho \gamma }{p}{\mathbf{M}}_3& -\frac{\rho \gamma }{p^2}((1-|\mathbf{v}{|}^2)\mathbf{B}+(\mathbf{v}\cdot \mathbf{B})\mathbf{v})\\ -\frac{\gamma }{p}& -\frac{\rho {\gamma }^3}{p}{\mathbf{v}}^{\top }& 0_3^{\top }& \frac{\rho \gamma }{p^2}\end{array}\right),$$

$$\frac{\partial \mathbf{V}}{\partial \mathbf{W}}=\left(\begin{array}{cccc}\rho & \left(\rho +\frac{p}{\Gamma -1}\right)\gamma {\mathbf{v}}^{\top }& 0_3^{\top }& \left(\rho +\frac{p}{\Gamma -1}\right)\gamma \\ 0_3& \frac{p}{\rho \gamma }{\mathbf{I}}_3& {\mathbf{O}}_3& \frac{p}{\rho \gamma }\mathbf{v}\\ 0_3& \frac{p\gamma }{\rho }{\mathbf{M}}_5& \frac{p\gamma }{\rho }\left({\mathbf{I}}_3-\mathbf{v}{\mathbf{v}}^{\top }\right)& \frac{p\gamma }{\rho }\left((1+|\mathbf{v}{|}^2)\mathbf{B}-2(\mathbf{v}\cdot \mathbf{B})\mathbf{v}\right)\\ p& ph\gamma {\mathbf{v}}^{\top }& 0_3^{\top }& ph\gamma \end{array}\right),$$

$${\mathcal{Q}}_i^{\prime }(\mathbf{V})=\left(\frac{(\Gamma -S)\gamma v_1}{\Gamma -1},-\frac{\rho {\gamma }^{3}S}{\Gamma -1}\left((1-|\mathbf{v}{|}^2){\mathbf{e}}_1^{\top }+v_1{\mathbf{v}}^{\top }\right),0,0,0,-\frac{\rho \gamma v_1}{p(\Gamma -1)}\right).$$

$${\mathcal{E}}^{\prime }(\mathbf{U})\frac{\partial {\mathbf{F}}_1}{\partial \mathbf{V}}=\left(\frac{(\Gamma -S)\gamma v_1}{\Gamma -1},-\frac{\rho {\gamma }^{3}S}{\Gamma -1}\left((1-|\mathbf{v}{|}^2){\mathbf{e}}_1^{\top }+v_1{\mathbf{v}}^{\top }\right),-\frac{\rho \gamma }{p}(\mathbf{v}\cdot \mathbf{B}),0,0,-\frac{\rho \gamma v_1}{p(\Gamma -1)}\right).$$

$${\mathcal{Q}}^{\prime}_{1}({\bf U})-{\mathcal{E}}^{\prime}({\bf U}){\bf F}^{\prime}_{1}({\bf U})-\frac{\rho\gamma}{p}({\bf v}\cdot{\bf B})B_{1}^{\prime}({\bf U})=\big{(}{\mathcal{Q}}^{\prime}_{1}({\bf V})-{\mathcal{E}}^{\prime}({\bf U}){\bf F}^{\prime}_{1}({\bf V})-\frac{\rho\gamma}{p}({\bf v}\cdot{\bf B})B_{1}^{\prime}({\bf V})\big{)}\frac{\partial{\bf V}}{\partial{\bf U}}={\bf 0}.$$

$$\frac{\partial \mathbf{U}}{\partial t}+\sum _{i=1}^d\frac{\partial {\mathbf{F}}_i(\mathbf{U})}{\partial x_i}=-\mathbf{S}(\mathbf{U})\nabla \cdot \mathbf{B},$$

$$\mathbf{S}(\mathbf{U}):={\left(0,(1-|\mathbf{v}{|}^2){\mathbf{B}}^{\top }+(\mathbf{v}\cdot \mathbf{B}){\mathbf{v}}^{\top },{\mathbf{v}}^{\top },\mathbf{v}\cdot \mathbf{B}\right)}^{\top }.$$

$${\varphi }^{\prime }(\mathbf{W})\cdot \mathbf{W}=\varphi (\mathbf{W}).$$

$${\mathbf{S}}^{\top }={\varphi }^{\prime }(\mathbf{W}).$$

$${\varphi }^{\prime }(\mathbf{V})=\left(\frac{\gamma }{p}(\mathbf{v}\cdot \mathbf{B}),\frac{\rho {\gamma }^3}{p}\left((1-|\mathbf{v}{|}^2){\mathbf{B}}^{\top }+(\mathbf{v}\cdot \mathbf{B}){\mathbf{v}}^{\top }\right),\frac{\rho \gamma }{p}{\mathbf{v}}^{\top },-\frac{\rho \gamma }{p^2}(\mathbf{v}\cdot \mathbf{B})\right),$$

$${\varphi }^{\prime }(\mathbf{W})={\varphi }^{\prime }(\mathbf{V})\frac{\partial \mathbf{V}}{\partial \mathbf{W}}=\left(0,(1-|\mathbf{v}{|}^2){\mathbf{B}}^{\top }+(\mathbf{v}\cdot \mathbf{B}){\mathbf{v}}^{\top },{\mathbf{v}}^{\top },\mathbf{v}\cdot \mathbf{B}\right)={\mathbf{S}}^{\top }.$$

$$\displaystyle\varphi:={\bf W}^{\top}{\bf U}-{\mathcal{E}}=\rho\gamma+\frac{\rho\gamma}{2p}\Big{(}(1-|{\bf v}|^{2})|{\bf B}|^{2}+({\bf v}\cdot{\bf B})^{2}\Big{)}$$

$${\psi }_i:={\mathbf{W}}^{\top }{\mathbf{F}}_i-{\mathcal{Q}}_i+\varphi B_i,i=1,\dots ,d$$

$$\mathbf{U}={\phi }^{\prime }(\mathbf{W}{)}^{\top },{\mathbf{F}}_i={\psi }_i^{\prime }(\mathbf{W}{)}^{\top }-B_i{\varphi }^{\prime }(\mathbf{W}),i=1,\dots ,d.$$

$$\varphi^{\prime\prime}({\bf W})\frac{\partial{\bf W}}{\partial t}+\sum_{i=1}^{d}\Big{(}\psi_{i}^{\prime\prime}({\bf W})-B_{i}\phi^{\prime\prime}({\bf W})\Big{)}\frac{\partial{\bf W}}{\partial x_{i}}={\bf 0},$$

Full text

{AMS}

35L65, 65M12, 65M06, 76W05, 76Y05

1 Introduction

Magnetohydrodynamics (MHDs) describe the dynamics of electrically-conducting fluids in the presence of magnetic field and play an important role in many fields including astrophysics, plasma physics and space physics. In many cases, astrophysics and high energy physics often involve fluid flow at nearly speed of light so that the relativistic effect should be taken into account. Relativistic MHDs (RMHDs) have applications in investigating astrophysical scenarios from stellar to galactic scales, e.g., gamma-ray bursts, astrophysical jets, core collapse super-novae, formation of black holes, and merging of compact binaries.

In the $d$-dimensional space, the governing equations of special RMHDs can be written as a system of hyperbolic conservation laws

[TABLE]

along with an additional divergence-free condition on the magnetic field

[TABLE]

where $d=1,$ $2$ or $3$. Here we employ the geometrized unit system so that the speed of light $c$ in vacuum is equal to one. In (1), the conservative vector ${\bf U}=(D,{\bf m}^{\top},{\bf B}^{\top},E)^{\top}$, and the flux in the $x_{i}$-direction is defined by

[TABLE]

with the mass density $D=\rho\gamma$, the momentum density vector ${\bf m}=(\rho h{\gamma^{2}}+|{\bf B}|^{2}){\bf v}-({\bf v}\cdot{\bf B}){\bf B}$, the energy density $E=\rho h\gamma^{2}-p_{tot}+|{\bf B}|^{2}$, and the vector ${\bf e}_{i}$ denoting the $i$-th column of the unit matrix of size $3$. Additionally, $\rho$ is the rest-mass density, ${\bf v}=(v_{1},v_{2},v_{3})^{\top}$ denotes the fluid velocity vector, $\gamma=1/\sqrt{1-|{\bf v}|^{2}}$ is the Lorentz factor, $p_{tot}$ is the total pressure containing the gas pressure $p$ and magnetic pressure $p_{m}:=\frac{1}{2}\left(\gamma^{-2}|{\bf B}|^{2}+({\bf v}\cdot{\bf B})^{2}\right)$, $h=1+e+\frac{p}{\rho}$ represents the specific enthalpy, and $e$ is the specific internal energy. We consider the ideal equation of state $p=(\Gamma-1)\rho e$ to close the system (1), where $\Gamma$ is constant and denotes the adiabatic index.

The system (1) involves strong nonlinearity, making its analytic treatment quite difficult. Numerical simulation is a primary approach to explore physical laws in RMHDs. The numerical study of RMHDs may date back to the 1990s, to the best of our knowledge. For instance, the time-dependent morphological evolution of a RMHD jet was computed in [61], based on a divergence-free formation of the Maxwell equations coupled to fluids [59] with the numerical method proposed in [60]. During the past few decades, various numerical schemes have been developed for the RMHD equations, including but not limited to: the total variation diminishing scheme [3], adaptive mesh methods [58, 35], discontinuous Galerkin methods [67, 68], physical-constraints-preserving schemes [66], entropy limited approach [32], etc. Systematic review of numerical RMHD schemes is beyond the scope of the present paper, and we refer interested readers to the review articles [26, 44]. Besides the standard difficulty in solving the nonlinear hyperbolic systems, an additional numerical challenge for the RMHD system (1) comes from the divergence-free condition (2), which is also involved in the non-relativistic ideal MHD system. Numerical preservation of (2) is highly non-trivial (for $d\geq 2$) but crucial for the robustness of numerical computations. Numerical experiments and analysis in the non-relativistic MHD case indicated that violating the divergence-free condition (2) may lead to numerical instability and nonphysical solutions [9, 4, 63]. Various numerical techniques were proposed to reduce or control the effect of divergence error; see, e.g., [20, 50, 56, 16, 4, 41, 42, 64, 65].

Due to the nonlinear hyperbolic nature of the RMHD equations (1), solutions of (1) can be discontinuous with the presence of shocks or contact discontinuities. This leads to the consideration of weak solutions. However, weak solutions may not be unique. To select the “physically relevant” solution among all weak solutions, entropy conditions are usually imposed as the admissibility criterion. In the case of RMHD equations (1), there is a natural entropy condition arising from the second law of thermodynamics which should be respected. It is natural to seek entropy stable numerical schemes which satisfy a discrete version of that entropy condition. Entropy stable numerical methods ensure that the entropy is conserved in smooth regions and dissipated across discontinuities. Thus, the numerics precisely follow the physics of the second law of thermodynamics and can be more robust. Moreover, entropy stable schemes also allow one to limit the amount of dissipation added to the schemes to guarantee the entropy stability. For the above reasons, developing entropy stable schemes for RMHD equations (1) is meaningful and highly desirable.

Entropy stability analysis was well studied for first-order accurate schemes and scalar conservation laws [14, 34, 47, 48]. For systems of hyperbolic conservation laws, much attention was paid to explore entropy stable schemes with entropy stability focused on single given entropy function. Tadmor [53, 54] established the framework of entropy conservative fluxes, which conserves entropy locally, and entropy stable fluxes for second-order schemes. Lefloch, Mercier and Rohde [40] proposed a procedure to construct higher-order accurate entropy conservative fluxes. Fjordholm, Mishra and Tadmor [23] developed a methodology for constructing high-order accurate entropy stable schemes, which combine high-order entropy conservative fluxes and suitable numerical dissipation operators based on essentially non-oscillatory (ENO) reconstruction that satisfies the sign property [24]. On the other hand, high-order entropy stable schemes have also been constructed via the summation-by-parts (SBP) procedure [21, 10, 28]. Entropy stable space–time discontinuous Galerkin (DG) schemes were studied in [5, 6, 37], where the exact integration is required for the proof of entropy stability. More recently, a framework for designing entropy stable high-order DG methods through suitable numerical quadrature was proposed in [13], where the SBP operators established in [21, 10, 28] were used and also generalized to triangles. There are other studies that address various aspects of entropy stability, including but not limited to [8, 25, 7, 36]. As a key ingredient in designing entropy stable schemes, the construction of affordable entropy conservative fluxes has received much attention. Although there is a general way to construct entropy conservative fluxes based on path integration [53, 54], the resulting fluxes may not have an explicit formula and can be computationally expensive. Explicit entropy conservative fluxes were derived for the Euler equations [38, 11, 51], shallow water equations [29], special relativistic hydrodynamics without magnetic field [19], and ideal non-relativistic MHD equations [12, 62], etc. Different from the Euler equations, the conservative non-relativistic MHD equations are not symmetrizable and do not admit an entropy [31, 5, 12]. Entropy symmetrization can be achieved by a modified system with an additional source [31, 5]. Based on the modified formulation, several entropy stable schemes were developed for non-relativistic MHDs; see, e.g., [12, 17, 62, 43, 18].

This paper aims to present entropy analysis and high-order accurate entropy stable schemes for the RMHD equations on Cartesian meshes. The difficulties in this work mainly come from the highly nonlinear coupling between RMHD equations (1), which leads to no explicit expression of the primitive variables $(\rho,{\bf v},p)$, entropy variables and fluxes ${\bf F}_{i}$ in terms of ${\bf U}$. The effort and findings in this work include the following:

We show that the conservative form (1) of RMHD equations is not symmetrizable and thus does not admit a thermodynamic entropy pair. A modified RMHD system is proposed by building the divergence-free condition (2) into the equations through adding a source term, which is proportional to $\nabla\cdot{\bf B}$. The modified RMHD system is symmetrizable, possesses a convex thermodynamic entropy pair, and is an analogue to the non-relativistic symmetrizable MHD system proposed by Godunov [31] and Powell [50]. The Godunov-type entropy symmetrization procedure was also briefly described in [57] and was first performed in [52] for RMHDs.111We thank Professor Yuri Trakhinin of Sobolev Institute of Mathematics for drawing our attention to references [52, 57] and for his explanations of these works which we were not aware of. 2. 2.

We derive a consistent two-point entropy conservative numerical flux for the symmetrizable RMHD equations. The flux has an explicit analytical formula, is easy to compute and thus is affordable. The key is to carefully select a set of parameter variables which can explicitly express the entropy variables and potential fluxes in simple form. Due to the presence of the source term in the symmetrizable RMHD system, the standard framework [54] of the entropy conservative flux is not applicable here and should be modified to take the effect of the source term into account. 3. 3.

We construct semi-discrete high-order accurate entropy conservative schemes and entropy stable schemes for symmetrizable RMHD equations on Cartesian meshes. The second-order entropy conservative schemes are built on the proposed two-point entropy conservative flux and a central difference discretization to the source term. Higher-order entropy conservative schemes are constructed by suitable linear combinations [40] of the proposed two-point entropy conservative flux. To ensure the entropy conservative property and high-order accuracy, a suitable discretization of the symmetrization source term should be employed, which is a key ingredient in these high-order schemes. Arbitrarily high-order accurate entropy stable schemes are obtained by adding suitable dissipation terms into the entropy conservative schemes.

It is noticed that high-order accurate entropy stable schemes were developed in a very recent work [19] for the relativistic hydrodynamics (RHDs) without the magnetic field, whose governing equations are sometimes also called the relativistic Euler equations. To the best of our knowledge, the present paper is the first to study entropy stable schemes for the RMHD equations. Note that, in the RMHD system, the magnetic field, evolved by the induction equations, is also nonlinearly involved in the momentum and energy equations. Due to the presence of magnetic field, mathematical structures of the RMHD system are more complicated than the relativistic Euler system, making the present derivations of entropy stable schemes more difficult. Different from the relativistic Euler equations whose conservative form admits a thermodynamic entropy pair [19], the conservative form (1) of the RMHD system does not admit the thermodynamic entropy pair, so that entropy symmetrization of the RMHD system must be derived to accommodate an entropy condition on the PDE level. Moreover, the symmetrization source term added in the modified RMHD system also renders additional technical challenges on the numerical level, since a suitable high-order discretization of that source term must be devised such that the resulting numerical schemes satisfy entropy stability. These make the exploration in this work significantly different from those in [19] and highly nontrivial. When the magnetic field is zero, the proposed entropy stable schemes for the RMHDs reduce to a class of entropy stable schemes for the RHDs. It is also worth mentioning that our entropy conservative fluxes are derived through a set of carefully selected parameter variables, which are different from those used in [19] for the RHDs, rendering the expression of our resulting fluxes simpler; see Remark 3.4 for details.

The paper is organized as follows. After giving the entropy analysis in Sect. 2, we derive the explicit two-point entropy conservative fluxes in Sect. 3. One-dimensional (1D) entropy conservative schemes and entropy stable schemes are constructed in Sect. 4, and their extensions to two dimensions (2D) are presented in Sect. 5. We conduct numerical tests in Sect. 6 to verify the performance of the proposed high-order accurate entropy stable schemes, before concluding the paper in Sect. 7.

2 Entropy Analysis

First, let us recall the definition of an entropy function.

Definition 2.1**.**

A convex function ${\mathcal{E}}({\bf U})$ is called an entropy for the system (1) if there exist entropy fluxes ${\mathcal{Q}}_{i}({\bf U})$ such that

[TABLE]

*where the gradients ${\mathcal{E}}^{\prime}({\bf U})$ and ${\mathcal{Q}}^{\prime}_{i}({\bf U})$ are written as row vectors, and ${\bf F}^{\prime}_{i}({\bf U})$ is the Jacobian matrix. The functions $({\mathcal{E}},{\mathcal{Q}}_{i})$ form an entropy pair. *

If a hyperbolic system of conservation laws admits an entropy pair, then the smooth solutions of the system should satisfy

[TABLE]

Solutions of a nonlinear hyperbolic system can be discontinuous. This leads to the consideration of weak solutions, which, however, may not be unique, and for (non-smooth) weak solutions the above identity does not hold in general. The following inequality (see, e.g., [15, page 83]) is usually imposed as an admissibility criterion to select the “physically relevant” solution among all weak solutions:

[TABLE]

which is interpreted in the sense of distribution and known as an entropy condition.

The existence of an entropy pair is closely related to the symmetrization of a hyperbolic system of conservation laws [30].

Definition 2.2**.**

The system (1) is said to be symmetrizable if there exists a change of variables ${\bf U}\to{\bf W}$ which symmetrizes it, that is, the equations (1) become

[TABLE]

*where the matrix $\frac{\partial{\bf U}}{\partial{\bf W}}$ is symmetric positive definite and $\frac{\partial{\bf F}_{i}}{\partial{\bf W}}$ is symmetric for all $i$. *

Lemma 2.3**.**

*A necessary and sufficient condition for the system (1) to possess a strictly convex entropy ${\mathcal{E}}({\bf U})$ is that there exists a change of dependent variables ${\bf U}\to{\bf W}$ which symmetrizes (1). *

The proof of Lemma 2.3 can be found in [30].

2.1 Entropy Function for RMHD Equations

It is natural to ask whether the RMHD equations (1) admit an entropy pair or not.

Let us consider the thermodynamic entropy $S=\ln(p\rho^{-\Gamma}).$ For smooth solutions, the RMHD equations (1) can be used to derive an equation for $\rho\gamma S$:

[TABLE]

Then under the divergence-free condition (2), the following quantities

[TABLE]

satisfy an additional conservation law, thus ${\mathcal{E}}$ may be an entropy function.

Theorem 2.4**.**

The entropy variables corresponding to the entropy function ${\mathcal{E}}$ are

[TABLE]

Proof 2.5**.**

Since ${\mathcal{E}}$ cannot be explicitly expressed by ${\bf U}$, direct derivation of ${\mathcal{E}}^{\prime}({\bf U})$ can be quite difficult. Here we consider the following primitive variables

[TABLE]

As ${\mathcal{E}}$ and $\bf U$ can be explicitly formulated in terms of ${\bf V}$, it is easy to derive that

[TABLE]

*where ${\bf 0}_{3}=(0,0,0)^{\top}$, and ${\bf O}_{3}$ and ${\bf I}_{3}$ denote the zero square matrix and the identity matrix of size $3$, respectively. One can verify that the vector $\bf W$ satisfies ${\bf W}^{\top}\frac{\partial{\bf U}}{\partial{\bf V}}={\mathcal{E}}^{\prime}({\bf V}),$ which implies ${\bf W}^{\top}={\mathcal{E}}^{\prime}({\bf V})\left(\frac{\partial{\bf U}}{\partial{\bf V}}\right)^{-1}={\mathcal{E}}^{\prime}({\bf V})\frac{\partial{\bf V}}{\partial{\bf U}}={\mathcal{E}}^{\prime}({\bf U}).$ *

However, the change of variable ${\bf U}\to{\bf W}$ fails to symmetrize the RMHD equations (1), and the functions $({\mathcal{E}},{\mathcal{Q}}_{i})$ defined in (7) do not form an entropy pair for the RMHD equations (1); see the following theorem.

Theorem 2.6**.**

For RMHD equations (1) and the entropy variables $\bf W$, one has

the change of variable ${\bf U}\to{\bf W}$ fails to symmetrize the RMHD equations (1), i.e., the matrix $\frac{\partial{\bf F}_{i}}{\partial{\bf W}}$ is not symmetric in general. 2. 2.

the functions $({\mathcal{E}},{\mathcal{Q}}_{i})$ defined in (7) satisfy

[TABLE]

which implies that $({\mathcal{E}},{\mathcal{Q}}_{i})$ do not satisfy the condition (3).

Proof 2.7**.**

We only prove the conclusions for $i=1$, as the proofs for $2\leq i\leq d$ are similar. Let us first show that $\frac{\partial{\bf F}_{i}}{\partial{\bf W}}$ is not symmetric in general. Because ${\bf F}_{1}$ cannot be formulated explicitly in terms of $\bf W$, we calculate the Jacobian matrix ${\bf F}_{1}({\bf W})$ with the aid of primitive variables ${\bf V}$. The Jacobian matrix $\frac{\partial{\bf F}_{1}}{\partial{\bf V}}$ is computed as

[TABLE]

with ${\bf e}_{1}=(1,0,0)^{\top}$, ${\bm{\beta}}_{1}=v_{1}(1-|{\bf v}|^{2}){\bf B}^{\top}+v_{1}({\bf v}\cdot{\bf B}){\bf v}^{\top}-B_{1}{\bf v}^{\top}$, and

[TABLE]

Since ${\bf W}$ can be explicitly expressed in terms of $\bf V$, one can derive that

[TABLE]

with ${\bf M}_{3}=(1-|{\bf v}|^{2}){\bf I}_{3}+{\bf v}{\bf v}^{\top},{\bf M}_{4}=({\bf v}\cdot{\bf B}){\bf v}{\bf v}^{\top}+\gamma^{-2}\left[({\bf v}\cdot{\bf B}){\bf I}_{3}+{\bf v}{\bf B}^{\top}-{\bf B}{\bf v}^{\top}\right].$ Then, we obtain $\frac{\partial{\bf V}}{\partial{\bf W}}$ by the inverse of the matrix $\frac{\partial{\bf W}}{\partial{\bf V}}$, i.e.

[TABLE]

with ${\bf M}_{5}=2{\bf B}{\bf v}^{\top}-({\bf v}\cdot{\bf B}){\bf I}_{3}-{\bf v}\left((1-|{\bf v}|^{2}){\bf B}^{\top}+({\bf v}\cdot{\bf B}){\bf v}^{\top}\right).$ By the chain rule $\frac{\partial{\bf F}_{1}}{\partial{\bf W}}=\frac{\partial{\bf F}_{1}}{\partial{\bf V}}\frac{\partial{\bf V}}{\partial{\bf W}}$, we get the expression of $\frac{\partial{\bf F}_{1}}{\partial{\bf W}}$ and find it is not symmetric in general. For example, the $(2,6)$ element of the Jacobian matrix $\frac{\partial{\bf F}_{1}}{\partial{\bf W}}$ is $\frac{p\gamma B_{2}}{\rho}(1+v_{1}^{2}-v_{2}^{2}-v_{3}^{2}),$ while the $(6,2)$ element is $\frac{p\gamma}{\rho}\left(B_{2}(1+v_{1}^{2}-v_{3}^{2})-B_{1}v_{1}v_{2}+B_{3}v_{2}v_{3}\right).$

Next, let us prove (11). Note that

[TABLE]

Using (8), (10) and (12), one can derive that

[TABLE]

It follows that

[TABLE]

*The proof for $i=1$ is complete. Similarly one can prove the conclusions for $2\leq i\leq d$. *

2.2 Modified RMHD Equations and Entropy Symmetrization

To address the above issue, we propose a modified RMHD system

[TABLE]

where

[TABLE]

The system (15) can be considered as the relativistic extension of the Godunov–Powell system [31, 50] in the ideal non-relativistic MHDs. The right-hand side term of (15) is proportional to $\nabla\cdot{\bf B}$. This means, at the continuous level, the modified form (15) and conservative form (1) are equivalent under the condition (2). However, the “source term” ${\bf S}({\bf U})\nabla\cdot{\bf B}$ modifies the character of the RMHD equations, making the system (15) symmetrizable and admit the convex entropy pair $({\mathcal{E}},{\mathcal{Q}}_{i})$, as shown below.

Lemma 2.8**.**

Let $\phi:=\frac{\rho\gamma}{p}({\bf v}\cdot{\bf B}).$ In terms of the entropy variables ${\bf W}$ in (8), $\phi({\bf W})$ is a homogeneous function of degree one, i.e.,

[TABLE]

In addition, the gradient of $\phi({\bf W})$ with respect to $\bf W$ equals ${\bf S}^{\top}$, i.e.,

[TABLE]

Proof 2.9**.**

Taking the gradient of $\phi$ with respect to the primitive variables ${\bf V}$ gives

[TABLE]

which together with (14) imply

[TABLE]

*Thus the identity (18) holds. Based on (8), (16) and (18), one has $\phi^{\prime}({\bf W})\cdot{\bf W}={\bf S}^{\top}\cdot{\bf W}=\frac{\rho\gamma}{p}({\bf v}\cdot{\bf B})=\phi,$ which gives (17). *

Theorem 2.10**.**

*For the entropy variables $\bf W$ in (8), the change of variable ${\bf U}\to{\bf W}$ symmetrizes the modified RMHD equations (15), and the functions $({\mathcal{E}},{\mathcal{Q}}_{i})$ defined in (7) form an entropy pair of the modified RMHD equations (15). *

Proof 2.11**.**

Define

[TABLE]

which satisfy $\psi_{i}=\varphi v_{i},~{}i=1,\dots,d.$ In terms of the variables ${\bf W}$, the gradients of $\varphi$ and $\psi_{i}$ satisfy the following identities

[TABLE]

Substituting (21) into (15), we can rewrite the modified RMHD equations as

[TABLE]

where the Hessian matrices $\varphi^{\prime\prime}({\bf W})$, $\psi_{i}^{\prime\prime}({\bf W})$ and $\phi^{\prime\prime}({\bf W})$ are all symmetric. Moreover, ${\mathcal{E}}({\bf U})$ is a convex function on ${\bf U}\in{\mathcal{G}}$ and the matrix $\varphi^{\prime\prime}({\bf W})$ is positive definite. Hence, the change of variables ${\bf U}\to{\bf W}$ symmetrizes the modified equations (15). According to Lemma 2.3, the system (15) possesses a strictly convex entropy ${\mathcal{E}}({\bf U})$.

We now show that the functions $({\mathcal{E}},{\mathcal{Q}}_{i})$ defined in (7) form an entropy pair of the modified RMHD system (15). The Jacobian matrix of the system (15) in $x_{i}$-direction is given by

[TABLE]

Thanks to (11), (17) and (18), we have

[TABLE]

*Thus, ${\mathcal{Q}}^{\prime}_{i}({\bf U})={\mathcal{E}}^{\prime}({\bf U}){\bf A}_{i}({\bf U}),$ and the functions $({\mathcal{E}},{\mathcal{Q}}_{i})$ form an entropy pair of (15). *

It is worth noting that Ruggeri and Strumia [52] (perhaps for the first time) also found the entropy variables and derived the entropy symmetrization for RMHD [52]. This symmetrization was clearly described by Trakhinin in [57, Section 3] and used for investigating the stability of shock waves in RMHD. We remark that the quasi-linear form (22) was also obtained by Anile and Pennisi [1] for studying the mathematical structure of test RMHDs. The main goal of this paper is to construct entropy-stable shock-capturing schemes, which are based on the modified RMHD system (15) instead of (22). Another important different form of symmetrization for the RMHD system, in term of variables $(p,\gamma{\bf v}^{\top},{\bf B}^{\top},S)^{\top}$ instead of entropy variables, was derived in [27] for investigating the stability of relativistic current-vortex sheets.

Remark 2.12**.**

*We note that, in the modified RMHD system (15), the induction equation is given by $\frac{\partial{\bf B}}{\partial t}+\nabla\cdot({\bf v}{\bf B}^{\top}-{\bf B}{\bf v}^{\top})+{\bf v}\nabla\cdot{\bf B}={\bf 0}.$ Taking the divergence of this equation gives $\nabla\cdot\left(\frac{\partial{\bf B}}{\partial t}+\nabla\cdot({\bf v}{\bf B}^{\top}-{\bf B}{\bf v}^{\top})+{\bf v}\nabla\cdot{\bf B}\right)=\frac{\partial}{\partial t}(\nabla\cdot{\bf B})+\nabla\cdot({\bf v}\nabla\cdot{\bf B})={\bf 0}.$ Combining the continuity equation of (15), it yields $\frac{\partial}{\partial t}\left(\frac{\nabla\cdot{\bf B}}{\rho\gamma}\right)+{\bf v}\cdot\nabla\left(\frac{\nabla\cdot{\bf B}}{\rho\gamma}\right)=0.$ This implies that the quantity $\frac{\nabla\cdot{\bf B}}{\rho\gamma}$ is constant along streamlines. A similar property also holds for the symmetrizable non-relativistic ideal MHD system proposed by Godunov [31] and Powell [50]. As first demonstrated numerically by Powell [50] in his eight-wave method, such a property implies that the error in divergence may be advected away by the flow, and suitable discretization of the symmetrizable MHD system may give robust numerical methods in the sense of controlled divergence error [50], entropy stability [12, 62, 43] and positivity preservation [64, 65], etc. *

Similar to the Powell source term (cf. [50, 12, 62, 43]) for the ideal MHD system, the source term added in the modified RMHD system (15) is also non-conservative but necessary to obtain a symmetrizable formulation and to accommodate an entropy condition on the PDE level. Hence, in order to achieve entropy stability, our schemes presented later are designed based on the modified RMHD system (15), which leads to additional technical challenges in discretizing the source term suitably to ensure its compatibility with the entropy stability. As mentioned in [12, 43] for the non-relativistic MHD system, there is still a conflict between the entropy stability which requires the non-conservative source term, and the conservation property which is lost due to the source term. The loss of conservation property leaves the possibility that it may lead to incorrect resolutions for some discontinuous problems, as observed in the ideal MHD case [56]. It will be interesting to explore whether entropy stable schemes can be constructed via the conservative formulation (1).

3 Derivation of Two-point Entropy Conservative Fluxes

In this section, we derive explicit two-point entropy conservative numerical fluxes, which will play an important role in constructing entropy conservative schemes and entropy stable schemes, for the RMHD equations (15). Similar to the non-relativistic MHD case [12], the standard definition [54] of the entropy conservative flux is not applicable here and should be slightly modified, due to the presence of the term ${\bf S}({\bf U})\nabla\cdot{\bf B}$ in the symmetrizable RMHD equations (15). Here we adopt a definition similar to the one proposed in [12] for the non-relativistic MHD equations.

Definition 3.1**.**

For $i=1,2,3,$ a consistent two-point numerical flux ${\bf F}_{i}^{\star}({\bf U}_{L},{\bf U}_{R})$ is entropy conservative if

[TABLE]

*where ${\bf W}$, $\phi$, $\psi_{i}$ and $B_{i}$ are the entropy variables defined in (8), the function defined in Lemma 2.8, the potential fluxes defined in (20), and the $x_{i}$ magnetic field component $B_{i}$, respectively. The subscripts $L$ and $R$ indicate that those quantities are corresponding to the “left” state ${\bf U}_{L}$ and the “right” state ${\bf U}_{R}$, respectively. *

Now, we would like to construct explicit entropy conservative fluxes ${\bf F}_{i}^{*}({\bf U}_{L},{\bf U}_{R})$ satisfying the condition (23). For notational convenience, we employ

[TABLE]

to denote, respectively, the jump and the arithmetic mean of a quantity. In addition, we also need the logarithmic mean

[TABLE]

which was first introduced in [38]. Then, one has following identities

[TABLE]

which will be frequently used in the following derivation.

Let us introduce the following set of variables

[TABLE]

with ${\bf u}:=\gamma{\bf v},$ ${\bf H}=\gamma^{-1}{\bf B}$ and $\beta=\rho/p$. Define

[TABLE]

An explicit two-point entropy conservative flux for $i=1$ is given below.

Theorem 3.2**.**

Let $\widehat{e}:=1+(\{\kern-3.125pt\{\beta\}\kern-3.125pt\}^{\ln})^{-1}/(\Gamma-1)$, $\vartheta:=\{\kern-3.125pt\{{\bf u}\cdot{\bf H}\}\kern-3.125pt\}$, $\widehat{\Theta}:=\{\kern-3.125pt\{\beta\}\kern-3.125pt\}\left(\{\kern-3.125pt\{\beta\}\kern-3.125pt\}+\widehat{\bf u}\cdot{\bm{\mu}}\right)>0$, $\Theta:=\widehat{\Theta}\left(|\widehat{\bf u}|^{2}-\{\kern-3.125pt\{\gamma\}\kern-3.125pt\}^{2}\right)<0$, $\widehat{p}_{tot}:={\{\kern-3.125pt\{\rho\}\kern-3.125pt\}}/{\{\kern-3.125pt\{\beta\}\kern-3.125pt\}}+\{\kern-3.125pt\{p_{m}\}\kern-3.125pt\}$, and

[TABLE]

Then the numerical flux given by

[TABLE]

*satisfies (23) for $i=1$. *

Proof 3.3**.**

Let ${\bf F}_{1}^{\star}({\bf U}_{L},{\bf U}_{R})=:(f_{1},f_{2},\dots,f_{8})^{\top}$ and ${\bf W}=(W_{1},\cdots,W_{8})^{\top}$. Then the condition (23) for $i=1$ can be rewritten as

[TABLE]

To determine the unknown components of ${\bf F}_{1}^{\star}$, we would like to expand each jump term in (29) into linear combination of the jumps of certain parameter variables. This will give us a linear algebraic system of eight equations for the unknown components $(f_{1},f_{2},\dots,f_{8})$. There are many options to choose different sets of parameter variables, which may result in different fluxes. Here we take ${\bf z}=(\rho,{\bf u}^{\top},{\bf H}^{\top},\beta)^{\top}$ as the parameter variables, to make the resulting formulation of ${\bf F}_{1}^{\star}$ simple.

In terms of the parameter variables ${\bf z}=(\rho,{\bf u}^{\top},{\bf H}^{\top},\beta)^{\top}$, the entropy variables $\bf W$ in (8) can be explicitly expressed as

[TABLE]

and $\psi_{1}$ and $\phi$ can be expressed as

[TABLE]

Then, using the identities (25)–(26), we rewrite the jump terms involved in (29) as

[TABLE]

with

[TABLE]

Substituting the above expressions of jumps into (29) gives

[TABLE]

where

[TABLE]

One can verify that the numerical flux ${\bf F}_{1}^{\star}({\bf U}_{L},{\bf U}_{R})$ given by (28) solves the linear system ${\bf M}{\bf F}_{1}^{\star}=\bm{\varsigma}$. Thus, ${\bf F}_{1}^{\star}$ satisfies (32), which is equivalent to (23) for $i=1$. Hence the numerical flux ${\bf F}_{1}^{\star}({\bf U}_{L},{\bf U}_{R})$ is entropy conservative in the sense of (23).

It is worth noting that

[TABLE]

which implies ${\bf F}_{1}^{\star}$ given by (28) is the unique solution to the linear system ${\bf M}{\bf F}_{1}^{\star}=\bm{\varsigma}$.

Finally, let us verify that the numerical flux ${\bf F}_{1}^{\star}({\bf U}_{L},{\bf U}_{R})$ given by (28) is consistent with the flux ${\bf F}_{1}$. If letting ${\bf U}_{L}={\bf U}_{R}={\bf U}$, then one has

[TABLE]

which imply ${\bf F}_{1}^{\star}({\bf U},{\bf U})={\bf F}_{1}({\bf U})$. Thus, the numerical flux ${\bf F}_{1}^{\star}({\bf U}_{L},{\bf U}_{R})$ given by (28) is consistent with the flux ${\bf F}_{1}$.

*Therefore, the numerical flux ${\bf F}_{1}^{\star}({\bf U}_{L},{\bf U}_{R})$ is consistent and is entropy conservative in the sense of (23). The proof is complete. *

Explicit entropy conservative fluxes ${\bf F}_{i}^{\star}({\bf U}_{L},{\bf U}_{R})$ for $i=2$ and $i=3$ can be constructed similarly or obtained by simply using a symmetric transformation based on the rotational invariance of the system (15). For example, ${\bf F}_{2}^{\star}({\bf U}_{L},{\bf U}_{R})$ is given by

[TABLE]

where

[TABLE]

Remark 3.4**.**

Taking ${\bf B}_{L}={\bf B}_{R}={\bf 0}$, we obtain a set of explicit entropy conservative fluxes for the relativistic hydrodynamic equations with zero magnetic field:

[TABLE]

*where $i=1,2,3$, and $\widehat{\rho h}:=\frac{\mathopen{\{\kern-3.49998pt\{}\rho\mathclose{\}\kern-3.49998pt\}}/\mathopen{\{\kern-3.49998pt\{}\beta\mathclose{\}\kern-3.49998pt\}}+\mathopen{\{\kern-3.49998pt\{}\rho\mathclose{\}\kern-3.49998pt\}}^{\ln}\widehat{e}}{\mathopen{\{\kern-3.49998pt\{}\gamma\mathclose{\}\kern-3.49998pt\}}^{2}-|\widehat{\bf u}|^{2}}.$ It is noticed that the expressions of the entropy conservative fluxes (35) are simpler than those derived in [19] via a different set of parameter variables $(\rho,\beta,{\bf v})$. In fact, the set of parameter variables (27) we employed are carefully selected from many possible candidate sets, so as to render the resulting fluxes in a simple form. *

4 Entropy Conservative Schemes and Entropy Stable Schemes in One Dimension

In this section, we construct entropy conservative schemes and entropy stable schemes for the 1D symmetrizable RMHD equations (15). To avoid confusing subscripts, we will use $x$ to denote the 1D spatial coordinate, $\bf F$ to represent the flux vector ${\bf F}_{1}$, and $\mathcal{Q}$ to represent the entropy flux ${\mathcal{Q}}_{1}$ in $x_{1}$-direction.

For simplicity, we consider a uniform mesh $x_{1}<x_{2}<\cdots<x_{N}$ with mesh size $x_{i+1}-x_{i}=\Delta x$. The midpoint values are defined as $x_{i+1/2}:=(x_{i}+x_{i+1})/2$ and the spatial domain is partitioned into cells $I_{i}=(x_{i-1/2},x_{i+1/2})$. A semi-discrete finite difference scheme of the 1D symmetrizable RMHD equations (15) can be written as

[TABLE]

where ${\bf U}_{i}(t)\approx{\bf U}(x_{i},t)$, the numerical flux $\widehat{\bf F}_{i+\frac{1}{2}}$ is consistent with the flux ${\bf F}({\bf U})$, and

[TABLE]

For notational convenience, the $t$ dependence of all quantities is suppressed below.

4.1 Entropy Conservative Schemes

The semi-discrete scheme (36) is said to be entropy conservative if its computed solutions satisfy a discrete entropy equality

[TABLE]

for some numerical entropy flux $\widetilde{\mathcal{Q}}_{i+\frac{1}{2}}$ consistent with the entropy flux ${\mathcal{Q}}$.

We introduce the following notations

[TABLE]

to denote the jump and the arithmetic mean of a quantity at the interface $x_{i+1/2}$.

4.1.1 Second-order entropy conservative scheme

Similar to the non-relativistic case [12], a second-order accurate entropy conservative scheme is obtained by taking $\widehat{\bf F}_{i+\frac{1}{2}}$ as the two-point entropy conservative flux and $\widehat{B}_{1,i+\frac{1}{2}}=\{\kern-3.125pt\{B_{1}\}\kern-3.125pt\}_{i+\frac{1}{2}}$.

Theorem 4.1**.**

If taking $\widehat{\bf F}_{i+\frac{1}{2}}$ as an entropy conservative numerical flux ${\bf F}_{1}^{\star}({\bf U}_{i},{\bf U}_{i+1})$ satisfying (23) and $\widehat{B}_{1,i+\frac{1}{2}}=\{\kern-3.125pt\{B_{1}\}\kern-3.125pt\}_{i+\frac{1}{2}}$, then the scheme (36), which becomes

[TABLE]

is entropy conservative, and the corresponding numerical entropy flux is given by

[TABLE]

*where ${\bf W}$, $\phi$ and $\psi_{1}$ are the entropy variables defined in (8), the function defined in Lemma 2.8, the potential flux defined in (20), respectively. *

Proof 4.2**.**

Using (20), one can easily verify that the above numerical entropy flux is consistent with the entropy flux $\mathcal{Q}$. Note that the numerical flux ${\bf F}_{1}^{\star}({\bf U}_{i},{\bf U}_{i+1})$ satisfies

[TABLE]

Using (8), (38), (18), (17) and (40) sequentially, we have

[TABLE]

*which implies the discrete entropy equality (37) for the numerical entropy flux (39). The proof is complete. *

We obtain an entropy conservative scheme (38) if the numerical flux ${\bf F}_{1}^{\star}({\bf U}_{i},{\bf U}_{i+1})$ given by (28) is used. Other entropy conservative fluxes satisfying (23) can also be used in (38) to obtain different entropy conservative schemes. Note that the second-order accuracy of the scheme (38) is only achieved on uniform mesh grids.

4.1.2 Higher-order entropy conservative schemes

The semi-discrete entropy conservative scheme (38) is only second-order accurate. By using the proposed entropy conservative flux (28) as building blocks, one can construct $2k$th-order accurate entropy conservative fluxes for any $k\in\mathbb{N}_{+}$; see [40]. These consist of linear combinations of second-order entropy conservative flux (28). Specifically, a $2k$th-order accurate entropy conservative flux is defined as

[TABLE]

where the constants $\alpha_{k,r}$ satisfy

[TABLE]

The symmetrizable RMHD equations (15) have a special source term, which should be treated carefully in constructing high-order accurate entropy conservative schemes. We find the key point is to accordingly approximate the spatial derivative $\partial_{x}B_{1}$ as $\frac{1}{\Delta x}\left({\widetilde{B}_{1,i+\frac{1}{2}}^{2k,\star}-\widetilde{B}_{1,i-\frac{1}{2}}^{2k,\star}}\right)\approx\left.\partial_{x}B_{1}\right|_{x=x_{i}},$ with $\widetilde{B}_{1,i+\frac{1}{2}}^{2k,\star}$ defined as a linear combination of $\frac{1}{2}(B_{1,i-s}+B_{1,i-s+r})$ similar to (41). Specifically, we set

[TABLE]

As an example, the fourth-order ($k=2$) version of $\widetilde{\bf F}_{i+\frac{1}{2}}^{2k,\star}$ and $\widetilde{B}_{1,i+\frac{1}{2}}^{2k,\star}$ is given by

[TABLE]

and the sixth-order ($k=3$) version is

[TABLE]

If taking $\widehat{\bf F}_{i+\frac{1}{2}}=\widetilde{\bf F}_{i+\frac{1}{2}}^{2k,\star}$ and $\widehat{B}_{1,i+\frac{1}{2}}=\widetilde{B}_{1,i+\frac{1}{2}}^{2k,\star}$, then the scheme (36), which becomes

[TABLE]

is entropy conservative and $2k$th-order accurate.

Theorem 4.3**.**

The scheme (46) is entropy conservative, and the corresponding numerical entropy flux is given by

[TABLE]

where the constants $\alpha_{k,r}$ are defined in (42), and the function $\widetilde{\mathcal{Q}}$ is defined as

[TABLE]

Proof 4.4**.**

First, one can use (42) to verify that the numerical entropy flux $\widetilde{\mathcal{Q}}_{i+\frac{1}{2}}^{2k,\star}$ is consistent with the entropy flux $\mathcal{Q}$. Using (8), (46), (18) and (17), we obtain

[TABLE]

It is observed that

[TABLE]

Therefore,

[TABLE]

with \Pi_{i,i+r}={\bf W}_{i}\cdot{\bf F}_{1}^{\star}({\bf U}_{i},{\bf U}_{i+r})+\frac{\phi_{i}}{2}\big{(}B_{1,i}+B_{1,i+r}\big{)} and \Pi_{i,i-r}={\bf W}_{i}\cdot{\bf F}_{1}^{\star}({\bf U}_{i-r},{\bf U}_{i})+\frac{\phi_{i}}{2}\big{(}B_{1,i}+B_{1,i-r}\big{)}. Note that

[TABLE]

where the property (23) of ${\bf F}_{1}^{\star}$ has been used in the penultimate equality sign. Similarly,

[TABLE]

Plugging (50) and (51) into (49) gives

[TABLE]

*which implies the discrete entropy equality (37) for the numerical entropy flux (47). The proof is complete. *

4.2 Entropy Stable Schemes

Entropy is conserved only if the solutions of the RMHD equations (15) are smooth. Entropy conservative schemes preserve the entropy and can work well in smooth regions. However, for solutions containing discontinuity where entropy is dissipated, entropy conservative schemes may produce oscillations; see, e.g., the numerical examples in [55, 22, 62] and Example 1 of the present paper. Consequently, some numerical dissipative mechanism should be added to ensure entropy stability.

The 1D semi-discrete scheme (36) is said to be entropy stable if its computed solutions satisfy a discrete entropy inequality

[TABLE]

for some numerical entropy flux $\widehat{\mathcal{Q}}_{i+\frac{1}{2}}$ consistent with the entropy flux ${\mathcal{Q}}$.

4.2.1 First-order entropy stable scheme

Let us add a numerical dissipation term to the entropy conservative flux ${\bf F}_{1}^{\star}$ and define

[TABLE]

where ${\bf F}_{1}^{\star}({\bf U}_{i},{\bf U}_{i+1})$ is an entropy conservative numerical flux, for example, the one given in (28), and ${\bf D}_{i+\frac{1}{2}}$ is any symmetric positive definite matrix.

Theorem 4.5**.**

The scheme (36) with $\widehat{B}_{1,i+\frac{1}{2}}=\{\kern-3.125pt\{B_{1}\}\kern-3.125pt\}_{i+\frac{1}{2}}$ and numerical flux (53) is entropy stable, and the corresponding numerical entropy flux is given by

[TABLE]

*where $\widetilde{\mathcal{Q}}_{i+\frac{1}{2}}^{\star}$ is defined by (39). *

Proof 4.6**.**

First, one can easily verify that the numerical entropy flux (54) is consistent with the entropy flux $\mathcal{Q}$. Substituting $\widehat{B}_{1,i+\frac{1}{2}}=\{\kern-3.125pt\{B_{1}\}\kern-3.125pt\}_{i+\frac{1}{2}}$ and numerical flux (53) into the scheme (36) and then following the proof of Theorem 4.1, we obtain

[TABLE]

Therefore,

[TABLE]

*which implies the discrete entropy inequality (52) for the numerical entropy flux (54). The proof is complete. *

In the computations, we evaluate the dissipation matrix ${\bf D}_{i+\frac{1}{2}}$ as follows. Let ${\bf U}_{i+\frac{1}{2}}^{\tt ave}$ be an average state at $x=x_{i+\frac{1}{2}}$ with the corresponding primitive variables

[TABLE]

where $\beta=\rho/p$, and $\{\kern-3.125pt\{\cdot\}\kern-3.125pt\}_{i+\frac{1}{2}}^{\rm ln}$ denotes the logarithmic mean defined in (24). Then the dissipation matrix ${\bf D}_{i+\frac{1}{2}}$ is evaluated at the above average state by ${\bf D}_{i+\frac{1}{2}}={\mathbf{R}}_{i+\frac{1}{2}}|{\bf\Lambda}_{i+\frac{1}{2}}|{\mathbf{R}}_{i+\frac{1}{2}}^{\top},$ where $|{\bf\Lambda}|$ is a diagonal matrix to be specified later, ${\mathbf{R}}$ is the matrix formed by the scaled (right) eigenvectors of the Jacobian matrix ${\bf A}_{1}({\bf U}):={\bf F}^{\prime}_{1}({\bf U})+{\bf S}({\bf U})B_{1}^{\prime}({\bf U})$, and it satisfies

[TABLE]

Let $\{\lambda_{\ell}\}_{1\leq\ell\leq 8}$ be the eight eigenvalues of the Jacobian matrix ${\bf A}_{1}$. Then, the diagonal matrix $|{\bf\Lambda}|$ can be chosen as

[TABLE]

which gives the Roe-type dissipation term in (53), or taken as

[TABLE]

which gives the Rusanov (also called generalized Lax-Friedrichs) type dissipation term in (53). We remark that the eigenvector scaling theorem [5] ensures that there exist scaled eigenvalues of ${\bf A}_{1}$ satisfying (56). For the computations of the eigenvalues and eigenvectors of the Jacobian matrix of the RMHD equations, see [2].

4.2.2 High-order entropy stable schemes

The entropy stable scheme (36) with $\widehat{B}_{1,i+\frac{1}{2}}=\{\kern-3.125pt\{B_{1}\}\kern-3.125pt\}_{i+\frac{1}{2}}$ and numerical flux (53) is only first order accurate in space, due to the presence of ${\mathcal{O}}(\Delta x)$ jump $\llbracket{\bf W}\rrbracket_{i+\frac{1}{2}}$ in the dissipation term. Towards achieving higher-order entropy stable schemes, we should use high-order dissipation operators with more accurate estimate of jump at cell interface. In this paper, we consider two approaches to construct high-order dissipation operators: the ENO based approach [23] and WENO based approach [7].

In the ENO based approach, we define high-order entropy stable fluxes as

[TABLE]

where $\widetilde{\bf F}_{i+1/2}^{2k,\star}$ is the $2k$th-order entropy conservative flux defined in (41), $\llangle{\bm{\omega}}\rrangle_{i+1/2}^{{\tt ENO}}:={\bm{\omega}}_{i+1/2}^{+}-{\bm{\omega}}_{i+1/2}^{-}$ with ${\bm{\omega}}_{i+1/2}^{-}$ and ${\bm{\omega}}_{i+1/2}^{+}$ denoting, respectively, the left and right limiting values of the scaled entropy variables $\bm{\omega}:={\bf R}^{\top}_{i+1/2}{\bf W}$ at interface $x_{i+1/2}$, obtained by $2k$th-order ENO reconstruction. The sign preserving property [24] of ENO reconstruction implies that

[TABLE]

We refer the readers to [23, Eq. (3.12)] for a more precise interpretation of the equality (60). In the WENO based approach [7], high-order accurate entropy stable fluxes can be defined as

[TABLE]

where the $\ell$th component of the vector $\llangle{\bm{\omega}}\rrangle_{i+1/2}^{{\tt WENO}}$ is computed by

[TABLE]

with ${\omega}_{\ell,i+1/2}^{-}$ and ${\omega}_{\ell,i+1/2}^{+}$ denoting, respectively, the left and right limiting values of $\omega_{\ell}$ at interface $x_{i+1/2}$ by using $(2k-1)$th-order WENO reconstruction. Although the standard WENO reconstruction may not satisfy the sign stability, the use of switch operator $\theta_{\ell,i+\frac{1}{2}}$ in (62), proposed in [7], ensures that

[TABLE]

Theorem 4.7**.**

The scheme (36), with $\widehat{B}_{1,i+\frac{1}{2}}=\widetilde{B}_{1,i+\frac{1}{2}}^{2k,\star}$ and the ENO-based numerical flux (59) or the WENO-based numerical flux (61), is entropy stable, and the corresponding numerical entropy flux is given by

[TABLE]

*where $\widetilde{\mathcal{Q}}_{i+1/2}^{2k,\star}$ is defined in (47), and $\llangle{\bm{\omega}}\rrangle_{i+1/2}$ is taken as $\llangle{\bm{\omega}}\rrangle_{i+1/2}^{\tt ENO}$ or $\llangle{\bm{\omega}}\rrangle_{i+1/2}^{\tt WENO}$ accordingly. *

Proof 4.8**.**

First, it is evident that the numerical entropy flux (64) is consistent with the entropy flux $\mathcal{Q}$. Substituting $\widehat{B}_{1,i+1/2}=\widetilde{B}_{1,i+1/2}^{2k,\star}$ and numerical flux (59) or (61) into the scheme (36), and then following the proof of Theorem 4.3, we obtain

[TABLE]

It follows that

[TABLE]

*where the last inequality is obtained by using (60) or (63) accordingly. Therefore, the computed solutions of the scheme satisfy a discrete entropy inequality (52) for the numerical entropy flux (64). *

Remark 4.9**.**

*The entropy stability of the scheme in Theorem 4.7 is established only at the semi-discrete level. With explicit time discretization by, for example, a Runge-Kutta method, we cannot prove the entropy stability of the resulting fully discrete schemes. The entropy stability of fully discrete schemes will be only demonstrated by numerical experiments in Sect. 6. *

5 Entropy Conservative Schemes and Entropy Stable Schemes in Two Dimensions

The 1D entropy conservative schemes and entropy stable schemes developed in Sect. 4 can be easily extended to the multidimensional cases on rectangular meshes. This section presents the extension for 2D RMHD equations (15) with $d=2$. To avoid confusing subscripts, we will use $(x,y)$ to denote the 2D spatial coordinates.

Let us consider a uniform 2D Cartesian mesh consisting of grid points $(x_{i},y_{j})=(i\Delta x,j\Delta y)$ for $i,j\in\mathbb{Z}$, where both spatial step-sizes $\Delta x$ and $\Delta y$ are given positive constants. A semi-discrete finite difference scheme for 2D modified RMHD equations (15) can be written as

[TABLE]

where ${\bf U}_{ij}(t)\approx{\bf U}(x_{ij},t)$, and $\widehat{\bf F}_{1,i+1/2,j}$ (resp. $\widehat{\bf F}_{2,i,j+1/2}$) is numerical flux consistent with ${\bf F}_{1}$ (resp. ${\bf F}_{2}$). For convenience, the $t$ dependence of all quantities is suppressed below.

5.1 Entropy Conservative Schemes

The semi-discrete scheme (65) is said to be entropy conservative if its computed solutions satisfy a discrete entropy equality

[TABLE]

for some numerical entropy fluxes $\widetilde{\mathcal{Q}}_{1,i+\frac{1}{2},j}$ and $\widetilde{\mathcal{Q}}_{2,i,j+\frac{1}{2}}$ consistent with the entropy flux ${\mathcal{Q}}_{1}$ and ${\mathcal{Q}}_{2}$, respectively.

Analogously to the one-dimensional case, for any $k\in\mathbb{N}_{+}$ we define

[TABLE]

where the constants $\alpha_{k,r}$ is given by (42).

Theorem 5.1**.**

For any $k\in\mathbb{N}_{+}$, the scheme (65) with

[TABLE]

is a $2k$th-order accurate entropy conservative scheme with the numerical entropy fluxes

[TABLE]

where the function $\widetilde{\mathcal{Q}}_{\ell}$ is defined by

[TABLE]

The proof is similar to those of Theorems 4.1–4.3 and is omitted here.

5.2 Entropy Stable Schemes

The 2D semi-discrete scheme (65) is said to be entropy stable if its computed solutions satisfy a discrete entropy inequality

[TABLE]

for some numerical entropy fluxes $\widehat{\mathcal{Q}}_{1,i+\frac{1}{2},j}$ and $\widehat{\mathcal{Q}}_{2,i,j+\frac{1}{2}}$ consistent with the entropy flux ${\mathcal{Q}}_{1}$ and ${\mathcal{Q}}_{2}$, respectively.

Analogously to the one-dimensional case, we define

[TABLE]

where ${\mathbf{R}}_{1,i+\frac{1}{2},j}$ (resp. ${\mathbf{R}}_{2,i,j+\frac{1}{2}}$) is the matrix formed by the scaled right eigenvectors of the Jacobian matrix ${\bf A}_{1}({\bf U}_{i+\frac{1}{2},j}^{\tt ave}):={\bf F}^{\prime}_{1}({\bf U}_{i+\frac{1}{2},j}^{\tt ave})+{\bf S}({\bf U}_{i+\frac{1}{2},j}^{\tt ave})B_{1}^{\prime}({\bf U}_{i+\frac{1}{2},j}^{\tt ave})$ (resp. ${\bf A}_{2}({\bf U}_{i,j+\frac{1}{2}}^{\tt ave}):={\bf F}^{\prime}_{2}({\bf U}_{i,j+\frac{1}{2}}^{\tt ave})+{\bf S}({\bf U}_{i,j+\frac{1}{2}}^{\tt ave})B_{2}^{\prime}({\bf U}_{i,j+\frac{1}{2}}^{\tt ave})$); the diagonal matrix $|{\bf\Lambda}_{1,i+\frac{1}{2},j}|$ (resp. $|{\bf\Lambda}_{2,i,j+\frac{1}{2}}|$) is defined as (57) or (58) with the eigenvalues of ${\bf A}_{1}({\bf U}_{i+\frac{1}{2},j}^{\tt ave})$ (resp. ${\bf A}_{2}({\bf U}_{i,j+\frac{1}{2}}^{\tt ave})$). Here ${\bf U}_{i+\frac{1}{2},j}^{\tt ave}$ and ${\bf U}_{i,j+\frac{1}{2}}^{\tt ave}$ denote the average states at the corresponding interfaces, analogously to the 1D case defined in (55). The “high-order accurate” jumps $\llangle{\bm{\omega}}\rrangle_{i+1/2,j}$ and $\llangle{\bm{\omega}}\rrangle_{i,j+1/2}$ in (70) are computed by ENO reconstruction or WENO reconstruction using switch operator, which can be performed precisely as in the 1D case, dimension by dimension.

Theorem 5.2**.**

The scheme (65), with $\widehat{B}_{1,i+1/2,j}=\widetilde{B}_{1,i+1/2,j}^{2k,\star},$ $\widehat{B}_{2,i,j+1/2}=\widetilde{B}_{2,i,j+1/2}^{2k,\star}$ and numerical fluxes (70), is entropy stable, and the numerical entropy flux is

[TABLE]

*where $\widetilde{\mathcal{Q}}_{1,i+1/2,j}^{2k,\star}$ and $\widetilde{\mathcal{Q}}_{2,i,j+1/2}^{2k,\star}$ are defined in (68). *

The proof is similar to that of Theorem 4.7 and thus is omitted here.

6 Numerical Experiments

In this section, we conduct numerical experiments on several 1D and 2D benchmark RMHD problems, to demonstrate the performance of the proposed high-order accurate entropy stable schemes and entropy conservative schemes. For convenience, we abbreviate the forth-order and sixth-order accurate ($k=2,3$ respectively) entropy conservative schemes as EC4 and EC6, respectively. The 1D and 2D entropy stable schemes with 1D numerical flux (59) or 2D numerical flux (70), using $k=2$ and fourth-order accurate ENO reconstruction, are abbreviated as ES4. The 1D and 2D entropy stable schemes with 1D numerical flux (61) or 2D numerical flux (70), using $k=3$ and fifth-order accurate WENO reconstruction, are abbreviated as ES5. Unless otherwise stated, all these semi-discrete schemes are equipped with a fourth-order accurate explicit total-variation-diminishing (TVD) Runge-Kutta time discretization to obtain fully discrete schemes, and we use a CFL number of $0.4$ and the ideal equation of state $p=(\Gamma-1)\rho e$ with $\Gamma=5/3$. In all the tests, the Rusanov-type dissipation operator is employed in the entropy stable schemes.

{expl}

[Smooth problem]This test is used to check the accuracy of our schemes. Consider a 1D smooth problem which describes Alfvén waves propagating periodically within the domain $[0,1]$ and has the exact solution

[TABLE]

where the vector $\bf V$ denotes the primitive variables as defined in (9), $\sigma=\sqrt{1+\rho h\gamma^{2}}$, $v_{2}(x,t)=0.2\sin(2\pi(x+t/\sigma))$, and $v_{3}(x,t)=0.2\cos(2\pi(x+t/\sigma))$.

In our computations, the domain $[0,1]$ is divided into $N$ uniform cells, and periodic boundary conditions are specified. The time step-size is taken as $\Delta t=0.4\Delta x^{\frac{6}{4}}$ and $\Delta t=0.4\Delta x^{\frac{5}{4}}$ for $\tt EC6$ and $\tt ES5$, respectively, in order to make the error in spatial discretization dominant in the present case. Tables 1 and 2 list the numerical errors at $t=0.5$ in the numerical velocity component $v_{2}$ and the corresponding convergence rates for the $\tt EC4$, $\tt ES4$, $\tt EC6$ and $\tt ES5$ schemes at different grid resolutions. The convergence behaviors for $v_{3}$, $B_{2}$ and $B_{3}$ are similar and omitted. We clearly observe that the expected convergence orders of the schemes are achieved accordingly. Figure 1 displays the evolution of discrete total entropy $\sum_{i}{\mathcal{E}}({\bf U}_{i}(t))\Delta x$, which approximates $\int_{0}^{1}{\mathcal{E}}({\bf U}(x,t))dx$ that remains conservative for this smooth problem. For the entropy stable schemes, the discrete total entropy is not constant and slowly decays due to the numerical dissipation, but we observe convergence with grid refinement, while for the entropy conservative schemes, it is nearly constant with time as expected.

{expl}

[1D Riemann problems]This example verifies the performance of the proposed entropy conservative and entropy stable schemes in resolving 1D RMHD wave configurations, by simulating four 1D Riemann problems (RPs) [3, 39]. The initial data of each RP comprise two different constant states separated by the initial discontinuity at $x=0$, see Table 3. Their exact solutions were provided in [33].

Figures 2, 3 and 4 give the numerical results obtained by using ES5 with $1000$ uniform cells for RPs I, II and III, respectively. We see that the wave structures including discontinuities are well resolved by ES5 and that the numerical solutions are in good agreement with the exact ones.

The RP IV is a variant of the RMHD shock-tube 2 proposed by Komissarov [39] with $\Gamma=\frac{4}{3}$. The numerical solutions at $t=1$, computed respectively by EC6 and ES5 on $1000$ uniform girds, are displayed in Figure 5. We see that EC6 produces high-frequency oscillations in its numerical solution, because it enforces entropy conservation so that the entropy dissipation at the shock does not take place. Oscillations are not produced in the numerical solution of ES5, owing to its WENO dissipative mechanism. In order to carefully verify the entropy conservative/stable property of the proposed schemes, we investigate the behavior of the discrete total entropy ${\mathcal{E}}_{\rm tot}(t):=\sum_{i}{\mathcal{E}}({\bf U}_{i}(t))\Delta x$. Figure 6 shows the evolution of ${\mathcal{E}}_{\rm tot}(t)$ computed by EC6 and ES5, respectively, with fourth-order TVD explicit Runge-Kutta time discretization and different CFL numbers (correspond to different time step-sizes). For ES5, we see that the entropy dissipates at almost the same rates for different time step-sizes, as expected. Although the semi-discrete scheme EC6 conserves the entropy in theory, the results of the fully discrete scheme show small-magnitude entropy dissipation, which is solely introduced by time discretization. This is evidenced by the observation that the dissipation magnitude for EC6 becomes smaller when smaller CFL numbers are used. To further verify that the entropy conservative property converges in time with the temporal order, we show the logarithmic plots of the decrement of the total entropy in Figure 7. Again, Figure 7(a) clearly shows that the decrement magnitude becomes smaller for smaller time step-sizes. Moreover, Figure 7(b) indicates that the entropy dissipates at a rate of order ${\mathcal{O}}(\Delta t^{5})$, ranging from $\sim 10^{-3}$ with ${\rm CFL}=0.4$ to $\sim 10^{-8}$ with ${\rm CFL}=0.05$. This demonstrates the convergence of the entropy conservative property and confirms that the semi-discrete scheme EC6 does conserve entropy.

{expl}

[Blast problem]Blast problem is a benchmark test for RMHD numerical schemes. Our setup is the same as in [45, 67]. Initially, the computational domain $[-6,6]^{2}$ is filled with a homogeneous gas at rest with adiabatic index $\Gamma=4/3$. The explosion zone ($r<0.8$) has a density of $10^{-2}$ and a pressure of $1$, while the ambient medium ($r>1$) has a low density of $10^{-4}$ and a low pressure of $5\times 10^{-4}$, where $r=\sqrt{x^{2}+y^{2}}$. A linear taper is applied to the density and pressure for $r\in[0.8,1]$. The magnetic field is initialized in the $x$-direction as $0.1$.

Our numerical results at $t=4$, obtained by using ES5 on the mesh of $400\times 400$ uniform cells, are shown in Figure 8. We observe that the wave pattern of the configuration is composed by two main waves, an external fast and a reverse shock waves. The former is almost circular, while the latter is somewhat elliptic. The magnetic field is essentially confined between them, while the inner region is almost devoid of magnetization. Our numerical results agree quite well with those reported in [67, 66]. To validate the entropy stability of ES5, we compute the total entropy $\sum_{i,j}{\mathcal{E}}({\bf U}_{ij}(t))\Delta x\Delta y$, which should decrease with time if the scheme is entropy stable. The evolution of total entropy is displayed in the left of Figure 10 obtained by using ES5 at different gird resolutions. We clearly see a monotonic decay which indicates that the fully discrete scheme ES5 is entropy stable.

{expl}

[Orszag-Tang problem]This test simulates the relativistic version [58] of the classical Orszag-Tang problem [46]. Our setup is the same as in [58]. Initially, the computational domain $[0,2\pi]^{2}$ is filled with hot gas. We set the adiabatic index $\Gamma=4/3$, the initial pressure $p=10$, and the rest-mass density $\rho=1$. The initial velocity field of the fluid is given by ${\bf v}(x,y,0)=(-A\sin(y),A\sin(x),0)^{\top},$ where the parameter $A=0.99/\sqrt{2}$ so that the maximum velocity is $0.99$ (the corresponding Lorentz factor is about $7$). The magnetic field is initialized at ${\bf B}(x,y,0)=(-\sin y,\sin(2x),0)^{\top}$. Periodic conditions are specified at all the boundaries. Although the solution of this problem is smooth initially, complicated wave structures are formed as the time increases, and turbulence behavior will be produced eventually.

Figure 9 gives the numerical results obtained by using ES5 on $600\times 600$ uniform grids. In comparison with the results in [58], the complicated flow structures are well captured by ES5 with high resolution. To validate the entropy stability, the evolution of total entropy is shown in the right of Figure 10 obtained by using ES5 at different gird resolutions. We observe that the total entropy remains constant at initial times because initially the solution is smooth, and it starts to decrease at $t\approx 2$ when discontinuities start to form, as expected.

{expl}

[Shock cloud interaction problem]This problem describes the disruption of a high density cloud by a strong shock wave. Our setup is the same as in [35]. The computational domain is $[-0.2,1.2]\times[0,1]$, with the left boundary specified as inflow condition and the others as outflow conditions. Initially, a shock wave moves to the right from $x=0.05$, with the left and right states

[TABLE]

respectively. There exists a rest circular cloud centred at the point $(0.25,0.5)$ with radius 0.15. The cloud has the same states to the surrounding fluid except for a higher density of 30. Figure 11 displays the schlieren images of rest-mass density logarithm and magnetic pressure logarithm at $t=1.2$ obtained by using ES5 with $560\times 400$ uniform cells. One can see that the discontinuities are captured with high resolution, and the results agree well with those in [35, 68, 66].

7 Conclusions

In this paper, we presented rigorous entropy analysis and developed high-order accurate entropy stable schemes for the RMHD equations. We proved that the conservative RMHD equations (1) are not symmetrizable and thus do not admit an entropy pair. To address this issue, we first proposed a symmetrizable RMHD system (15) by building the divergence-free condition (2) into the RMHD equations through an additional source term. Based on the symmetrizable RMHD system, high-order accurate entropy stable finite difference schemes are developed on Cartesian meshes. These schemes are built on affordable explicit entropy conservative fluxes that are technically derived through carefully selected parameter variables, a special high-order discretization of the source term in the symmetrizable RMHD system, and suitable high-order dissipative operators based on (weighted) essentially non-oscillatory reconstruction to ensure the entropy stability. The accuracy and robustness of the proposed entropy stable schemes were demonstrated by benchmark numerical RMHD examples. The symmetrizable RMHD system is a relativistic analogue to the non-relativistic Godunov–Powell system [31, 49]. Its discovery would also enable one to generalize a large class of related numerical methods, including some other entropy stable schemes (e.g. [43]) and Powell’s eight-wave methods, to RMHDs.