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

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 possess an entropy pair. To address this issue, a symmetrizable RMHD system, which admits a convex entropy pair, is proposed by adding a source term into the equations. Arbitrarily high-order accurate entropy stable finite difference schemes are then 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 benchmark numerical tests demonstrate the accuracy and robustness of the proposed entropy stable schemes of the symmetrizable RMHD equations.



There are no comments yet.


page 28

page 29

page 31


High-order accurate entropy stable finite difference schemes for the shallow water magnetohydrodynamics

This paper develops the high-order accurate entropy stable (ES) finite d...

High-order accurate entropy stable nodal discontinuous Galerkin schemes for the ideal special relativistic magnetohydrodynamics

This paper studies high-order accurate entropy stable nodal discontinuou...

Entropy stable non-oscillatory fluxes: An optimized wedding of entropy conservative flux with non-oscillatory flux

This work settles the problem of constructing entropy stable non-oscilla...

A Kernel-Based Explicit Unconditionally Stable Scheme for Hamilton-Jacobi Equations on Nonuniform Meshes

In <cit.>, the authors developed a class of high-order numerical schemes...

HOLZ: High-Order Entropy Encoding of Lempel-Ziv Factor Distances

We propose a new representation of the offsets of the Lempel-Ziv (LZ) fa...

High-order Positivity-preserving L2-stable Spectral Collocation Schemes for the 3-D compressible Navier-Stokes equations

This paper extends a new class of positivity-preserving, entropy stable ...

Stable High-Order Cubature Formulas for Experimental Data

In many applications, it is impractical—if not even impossible—to obtain...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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 -dimensional space, the governing equations of special RMHDs can be written as a system of hyperbolic conservation laws


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


where or . Here we employ the geometrized unit system so that the speed of light in vacuum is equal to one. In (1.1

), the conservative vector

, and the flux in the -direction is defined by

with the mass density , the momentum density vector , the energy density , and the vector denoting the -th column of the unit matrix of size . Additionally, is the rest-mass density, denotes the fluid velocity vector, is the Lorentz factor, is the total pressure containing the gas pressure and magnetic pressure , represents the specific enthalpy, and is the specific internal energy. Different from the non-relativistic case, the variables and depend on the magnetic field nonlinearly. The system (1.1) is closed with an equation of state. In this paper, we consider the ideal equation of state , where is constant and denotes the adiabatic index.

The system (1.1) involves strong nonlinearity, making its analytic treatment quite difficult. Numerical simulation is a primary approach to explore physical laws in RMHDs. Since nearly 2000s, numerical study of RMHD has attracted considerable attention, and various numerical methods have been developed for the RMHD equations, for example, the total variation diminishing scheme [2], adaptive mesh methods [49, 31], discontinuous Galerkin methods [57, 58], entropy limited appraoch [28] and the high-order schemes based on subluminal reconstruction [4], etc. Recently, physical-constraints-preserving schemes were proposed for relativistic hydrodynamics [55, 51] and RMHDs [56], A systematic review of numerical RMHD schemes can be found in [39]. Besides the standard difficulty in solving the nonlinear hyperbolic systems, an additional numerical challenge for the RMHD system (1.1) comes from the divergence-free condition (1.2), which is also involved in the non-relativistic ideal MHD system. Numerical preservation of (1.2) is highly non-trivial (for ) but crucial for the robustness of numerical computation. Numerical experiments and analysis in the non-relativistic MHD case indicated that violating the divergence-free condition (1.2) may lead to numerical instability and nonphysical solutions [9, 3, 52]. Various numerical techniques were proposed to reduce or control the effect of divergence error; see, e.g., [19, 44, 48, 15, 3, 36, 37, 58, 53, 54].

Due to the nonlinear hyperbolic nature of the RMHD equations (1.1), solutions of (1.1) can be discontinuous with the presence of shocks or contact discontinuities. This leads to consideration of weak solutions. However, the 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.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 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.1) is meaningful and highly desirable.

Entropy stability analysis was well studied for first-order accurate schemes and scalar conservation laws with all entropy functions in early work [14, 30, 42, 43]. In recent years, there has been much interest in exploring high-order accurate entropy stable schemes for systems of hyperbolic conservation laws, with entropy stability focused on single given entropy function. Tadmor [46, 47] established the framework of entropy conservative fluxes, which conserves entropy locally, and entropy stable fluxes for second-order schemes. Lefloch, Mercier and Rohde [35] proposed a produce to construct higher-order accurate entropy conservative fluxes. Fjordholm, Mishra and Tadmor [21] 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 [22]. High-order entropy stable schemes have also been constructed via the summation-by-parts procedure [20, 10, 24]. Entropy stable space–time discontinuous Galerkin (DG) schemes were studied in [5, 6, 33], where the exact integration is required for the proof of entropy stability. More recently, a unified framework was proposed in [13] for designing provably entropy stable high-order DG methods through suitable numerical quadrature. There are other studies that address various aspects of entropy stability, including but not limited to [8, 23, 7, 32]. 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 flux based on path integration [46, 47], the resulting flux may not have an explicit formula and can be computationally expensive. Explicit entropy conservative fluxes were derived for the Euler equations [34, 11, 45], shallow water equations [25], special relativistic hydrodynamics without magnetic field [18], and ideal non-relativistic MHD equations [12, 50], etc. Different from the Euler equations, the conservative form of ideal non-relativistic MHD equations is not symmetrizable and does not admit an entropy [27, 5, 12]. Entropy symmetrization can be achieved by a modified non-relativistic MHD system with an additional source [27, 5]. Based on the modified formulations, several entropy stable schemes were well developed for non-relativistic MHDs; see, e.g., [12, 16, 50, 38, 17].

This paper aims to present entropy analysis and high-order accurate entropy stable finite difference schemes for RMHD equations on Cartesian meshes. In comparison with the non-relativistic MHD case, the difficulties in this work mainly come from (i) the unclear symmetrizable form of RMHD equations that admits an entropy pair, and (ii) the highly nonlinear coupling between RMHD equations (1.1), which leads to no explicit expression of the primitive variables , entropy variables and the flux in terms of . The effort and findings in this work include the following:

  1. We show that the conservative form (1.1) of RMHD equations is not symmetrizable and thus does not admit an entropy pair. A modified RMHD system is proposed by building the divergence-free condition (1.2) into the equations through adding a source term, which is proportional to . The modified RMHD system is symmetrizable, possesses the natural entropy pair, and can be considered as the relativistic extension of the non-relativistic symmetrizable MHD system proposed by Godunov [27]. As a result, it has almost all good features that the non-relativistic symmetrizable MHD system possesses.

  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 [47] of the entropy conservative flux is not applicable here and should be modified to take the effect of the source term into account.

  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 type discretization to the source term. Higher-order entropy conservative schemes are constructed by suitable linear combinations [35] of the proposed two-point entropy conservative flux. We find that, to ensure the entropy conservative property and high-order accuracy, a particular discretization of the symmetrization source term should be employed, which becomes 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.

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 entropy function.

Definition 2.1.

A convex function is called an entropy for the system (1.1) if there exist entropy fluxes such that


where the gradients and are written as row vectors, and is the Jacobian matrix. The functions form an entropy pair.

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

For non-smooth solutions, the above identity does not hold in general and is replaced with an inequality


which is interpreted in the week sense and known as the entropy condition.

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

Definition 2.2.

The system (1.1) is said to be symmetrizable if there exists a change of variables which symmetrizes it, that is, the equations (1.1) become


where the matrix is symmetric positive definite and is symmetric for all .

Lemma 2.1.

A necessary and sufficient condition for the system (1.1) to possess a strictly convex entropy is that there exists a change of dependent variables which symmetrizes (1.1).

The proof of Lemma 2.1 can be found in [26].

2.1 Entropy Function for RMHD Equations

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

Let us consider the thermodynamic entropy For smooth solutions, the RMHD equations (1.1) can be used to derive an equation for :


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


satisfy an additional conservation law for smooth solution, thus may be an entropy function.

Theorem 2.1.

The entropy variables corresponding to the entropy function are given by


Since cannot be explicitly expressed by , direct derivation of can be quite difficult. Here we consider the following primitive variables


As and can be explicitly formulated in terms of , it is easy to derive that



where , and and

denote the zero square matrix and the identity matrix of size

, respectively. One can verify that the vector satisfies which implies The proof is complete.

However, the change of variable fails to symmetrize the RMHD equations (1.1), and the functions defined in (2.5) do not form an entropy pair for the RMHD equations (1.1); see the following theorem.

Theorem 2.2.

For the RMHD equations (1.1) and the entropy variables in (2.6), one has

  1. the change of variable fails to symmetrize the RMHD equations (1.1), i.e., the matrix is not symmetric in general.

  2. the functions defined in (2.5) satisfy


    which implies that do not satisfy the condition (2.1).


We only prove the conclusions for , as the proofs for are similar.

Let us first show that is not symmetric in general. Because cannot be formulated explicitly in terms of , we calculate the Jacobian matrix with the aid of primitive variables . The Jacobian matrix is computed as


with , , and

Since can be explicitly expressed in terms of , one can derive that


with Then, we obtain by the inverse of the matrix , i.e.



By the chain rule

, we get the expression of and find it is not symmetric in general. For example, the element of the Jacobian matrix is

while the element is

Next, let us prove (2.9). Note that

Using (2.6), (2.8) and (2.10), one can derive that

It follows that


The proof for is complete. Similar arguments yield the conclusions for .

2.2 Modified RMHD Equations and Entropy Symmetrization

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




The system (2.13) can be considered as the relativistic extension of the Godunov–Powell system [27, 44] in the ideal non-relativistic MHDs. The right-hand side term of (2.13) is proportional to . This means, at the continuous level, the modified form (2.13) and conservative form (1.1) are equivalent under the condition (1.2). However, the “source term” modifies the character of the RMHD equations, making the system (2.13) symmetrizable and admit the convex entropy pair , as shown below.

Lemma 2.2.

Let In terms of the entropy variables in (2.6), is a homogeneous function of degree one, i.e.,


In addition, the gradient of with respect to equals , i.e.,


Taking the gradient of with respect to the primitive variables gives

which together with (2.12) imply

Thus the identity (2.16) holds. Based on (2.6), (2.14) and (2.16), one has which gives (2.15).

Theorem 2.3.

For the entropy variables in (2.6), the change of variable symmetrizes the modified RMHD equations (2.13), and the functions defined in (2.5) form an entropy pair of the modified RMHD equations (2.13).




which satisfy In terms of the variables , the gradients of and satisfy the following identities


Substituting (2.19) into (2.13), we can rewrite the modified RMHD equations as


where the Hessian matrices , and are all symmetric. Moreover, is a convex function on and the matrix is positive definite. Hence, the change of variables symmetrizes the modified RMHD equations (2.13). According to Lemma 2.1, the system (2.13) possesses a strictly convex entropy .

We now show that the functions defined in (2.5) form an entropy pair of the modified RMHD system (2.13). The Jacobian matrix of the system (2.13) in -direction is given by

Thanks to (2.9), (2.15) and (2.16), we have

Thus, and the functions form an entropy pair of (2.13).

Remark 2.1.

We note that, in the modified RMHD system (2.13), the induction equation is given by Taking the divergence of this equation gives

Combining the continuity equation of (2.13), it yields


Equation (2.21) implies that the quantity is constant along particle paths. Similar property holds for the Godunov–Powell modified system in the non-relativistic ideal MHDs [27, 44]. Following Powell’s perspective for the non-relativistic ideal MHDs [44], it reasonable to conjecture that suitable discretization of the modified RMHD system (2.13) can provide more stable numerical simulation as the numerical divergence error may be advected away by the flow.

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 (2.13). Similar to the non-relativistic MHD case [12], the standard definition [47] of the entropy conservative flux is not applicable here and should be slightly modified, due to the presence of the term in the symmetrizable RMHD equations (2.13). Here we adopt a definition similar to the one proposed in [12] for the non-relativistic MHD equations.

Definition 3.1.

For a consistent two-point numerical flux is entropy conservative if


where , , and are the entropy variables defined in (2.6), the function defined in Lemma 2.2, the potential fluxes defined in (2.18), and the magnetic field component , respectively. The subscripts and indicate that those quantities are corresponding to the “left” state and the “right” state , respectively.

Now, we would like to construct explicit entropy conservative fluxes satisfying the condition (3.1). For notational convenience, we employ

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

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


which will be frequently used in the following derivation.

Let us introduce the following set of variables


with and . Define

An explicit two-point entropy conservative flux for is given below.

Theorem 3.1.

Let , , , , , and

Then the numerical flux given by


satisfies (3.1) for .


Let and . Then the condition (3.1) for can be rewritten as


To determine the unknown components of , we would like to expand each jump term in (3.7) 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 . There are many options to choose different sets of parameter variables, which may result in different fluxes. Here we take as the parameter variables, to make the resulting formulation of simple.

In terms of the parameter variables , the entropy variables in (2.6) can be explicitly expressed as

and and can be expressed as

Then, using the identities (3.2)–(3.4), we rewrite the jump terms involved in (3.7) as



Substituting the above expressions of jumps into (3.7) gives