# High-order finite-difference entropy stable schemes for two-fluid
relativistic plasma flow equations

Deepak Bhoriya  Harish Kumar  Praveen Chandrashekar

###### Abstract

In this article, we propose high-order finite-difference entropy stable
schemes for the two-fluid relativistic plasma flow equations. This is achieved
by exploiting the structure of the equations, which consists of three
independent flux components. The first two components describe the ion and
electron flows, which are modeled using the relativistic hydrodynamics
equation. The third component is Maxwell’s equations, which are linear
systems. The coupling of the ion and electron flows, and electromagnetic
fields is via source terms only. Furthermore, we also show that the source
terms do not affect the entropy evolution. To design semi-discrete entropy
stable schemes, we extend the RHD entropy stable schemes in [1] to three
dimensions. This is then coupled with entropy stable discretization of the
Maxwell’s equations. Finally, we use SSP-RK schemes to discretize in time. We
also propose ARK-IMEX schemes to treat the stiff source terms; the resulting
nonlinear set of algebraic equations is local (at each discretization point).
These equations are solved using the Newton’s Method, which results in an
efficient method. The proposed schemes are then tested using various test
problems to demonstrate their stability, accuracy and efficiency.

###### Contents

  1. 1 Introduction
  2. 2 Two-fluid relativistic plasma flow equations
  3. 3 Analysis of continuous problem
  4. 4 Semi-discrete entropy stable schemes
    1. 4.1 Entropy stable numerical schemes for the fluid equations
      1. 4.1.1 High-order entropy conservative schemes for fluid equations
      2. 4.1.2 High-order entropy stable schemes for fluid equations
    2. 4.2 Discretization of Maxwell’s equations
  5. 5 Fully discrete scheme
    1. 5.1 Explicit schemes
    2. 5.2 ARK-IMEX schemes
      1. 5.2.1 Second order time discretization.
      2. 5.2.2 Third and fourth-order time discretization.
  6. 6 Numerical results
    1. 6.1 One-dimensional test cases
      1. 6.1.1 Accuracy test
      2. 6.1.2 Relativistic Brio–Wu test problem with finite plasma skin depth
      3. 6.1.3 Self-similar current sheet with finite resistivity
    2. 6.2 Two-dimensional test cases
      1. 6.2.1 Relativistic Orzag-Tang test case
      2. 6.2.2 Relativistic two-fluid blast problem
      3. 6.2.3 Relativistic two-fluid GEM challenge problem
  7. 7 Conclusion
  8. A Barth scaling of right eigenvectors
    1. A.1 Right eigenvectors
    2. A.2 Barth scaling and entropy scaled right eigenvectors
  9. B ARK3-IMEX coefficients table

Keywords:  Finite-difference entropy stable schemes, two-fluid relativistic
plasma flows, balance laws, IMEX-schemes

##  1 Introduction

Relativistic plasma flows play a central role in astrophysics. In these flows,
fluids moving with the speed comparable to the speed of light interact with
electric and magnetic fields. Some examples are pulsar winds, gamma-ray
bursts, relativistic jets from active galactic nuclei, and quasars
[2],[3],[4],[5]. To model relativistic plasma flows, often, equations of
relativistic Magneto-hydrodynamics (RMHD) are used [6],[7],[8],[9],[10],[11].
However, RMHD has limitations when modeling several astrophysical systems,
like pulsar, gamma-ray burst, etc., see [12] for a detailed discussion.
Following the development of non-relativistic two-fluid plasma flow equations
[13],[14],[15],[16],[17],[18],[19], several authors have considered an
analogous relativistic two-fluid plasma flow equations
[20],[21],[22],[23],[24].

Two-fluid relativistic plasma flow equations describe each fluid (ion and
electron) in plasma using the equations of special relativistic hydrodynamics
(RHD). The fluid components are coupled via electromagnetic quantities using
Lorentz force terms. Finally, the electric and magnetic fields are evolved
using Maxwell’s equations, where current and charges are described using fluid
variables. The resulting set of PDEs is a Hyperbolic system of Balance laws
with nonlinear flux and stiff source. Due to nonlinearity in the flux, the
solutions will exhibit discontinuities [25]. Hence, we need to consider the
weak solutions, which can be characterized using the Rankine-Hugoniot
condition across the discontinuities. As the weak solutions are non-unique,
entropy inequality is imposed to avoid nonphysical solutions. More recently,
even the entropy solutions for the systems have been shown to be non-unique
[26], still, entropy stability is one of the few nonlinear stability estimates
for the solutions, and hence it is desirable to have a numerical scheme which
replicates this stability at the discrete level.

As the flux and source terms are nonlinear, in general, it is not possible to
find an analytical solution for the two-fluid relativistic plasma flow
equations. Hence, the development of stable, efficient and accurate numerical
methods is highly desirable. However, in addition to the well-known
difficulties in designing stable numerical schemes for hyperbolic balance laws
[27],[25], there are several additional challenges in designing such a scheme
for two-fluid relativistic plasma flow equations. First, we need to compute
the primitive variables from the conservative variables as the analytical
expressions are not available. Second, we must satisfy additional constraints
for the magnetic and electric fields. Third, we have nonlinear stiff source
terms, and a suitable time implicit treatment of such source is highly
desirable. Furthermore, as discussed above, it is also desirable to have
entropy stability of the scheme.

Given these difficulties, there are few numerical methods for the two-fluid
relativistic plasma flow model. In [28], authors have simulated the
relativistic magnetic reconnection problem by assuming the symmetric motions
of the fluid species, which was further improved in [20] to include the
independent motion of the two-fluid species under a guide field. In [22],
authors used one-dimensional two-fluid relativistic plasma flow system to
investigate the importance of the electric-field-dominated regime in the
dynamics of a strongly magnetized plasma by simulating a standing shock front.
In [23], authors presented a third-order accurate explicit multidimensional
numerical scheme for special relativistic two-fluid plasma flow equations.
Here, the authors used the method of generalized Lagrange multiplier to keep
the divergence error small. Divergence-free and constraint preserving
multidimensional scheme has been designed in [12, 29].

In this article, we propose arbitrarily high-order finite-difference entropy
stable schemes for the relativistic two-fluid plasma flow model. High order
finite-difference entropy stable schemes for general hyperbolic systems were
first proposed in [30, 31]. Recently, several authors have developed entropy
stable finite-difference scheme for non-relativistic fluids and plasma flows
[15],[32],[33],[34]. For the relativistic fluids, Duan et.al [35], Biswas
[36], and Bhoriya et.al [1] have proposed entropy stable schemes.

We proceed as follows.

  * •

We first present the entropy framework for the two-fluid relativistic plasma
flow equations. We also prove that the source terms do not affect entropy
evolution.

  * •

To design entropy stable scheme, we exploit the structure of the flux. Similar
to the case of [15], the flux consists of three independent blocks. For the
fluid parts, we extend the entropy conservative schemes of [1] to include a
three-dimensional velocity field for each fluid component. This is then
coupled with a Rusanov’s solver based on higher-order Maxwell’s
discretization. The resultant scheme is then shown to be entropy conservative.

  * •

To make the scheme entropy stable, following [30], [15], [1], we design high-
order diffusion operator, for each component. As the source does not affect
entropy evolution, the semi-discrete scheme is then shown to be entropy
stable.

  * •

As the source terms are stiff, in addition to the explicit Runge-Kutta method,
we also propose IMEX schemes for the model, where flux contributions are
treated explicitly, and the source is treated implicitly. The resulting set of
nonlinear algebraic equations is local in each cell and are solved using
Newton’s method for nonlinear systems.

The rest of the article is organized as follows. In Section (2), we describe
the two-fluid relativistic plasma flow equations. In Section (3), we discuss
hyperbolicity and the entropy framework for the model. In Section (4), we
propose semi-discrete entropy stable numerical schemes for the model. This
includes the derivation of the fluid entropy conservative numerical flux and
the construction of the high-order entropy diffusion operators. We also
describe the discretization of the Maxwell’s equations. In Section (5), we
present the fully discrete scheme using time explicit and IMEX schemes. In
Section (6), we provide detailed numerical tests in one and two dimensions.

##  2 Two-fluid relativistic plasma flow equations

The two-fluid relativistic plasma flows consists of electrons and ions, whose
flows are governed by the equations of special relativistic flows. The
macroscopic variables for the fluids parts are density, fluid velocity and
pressure, denoted by $\rho_{\alpha},\ \mathbf{u}_{\alpha}$, and $\
p_{\alpha}$, respectively with $\alpha\in\\{i,\ e\\}$. Here, the subscripts
$\\{i,\ e\\}$ denote the variables corresponding to the ion and electron
species. These are coupled to Maxwell’s equations via source terms. The
complete system is given as follows:

$\displaystyle\frac{\partial(\rho_{i}\Gamma_{i})}{\partial
t}+\nabla\cdot(\rho_{i}\Gamma_{i}\mathbf{u}_{i})$ $\displaystyle=0,$ (1a)
$\displaystyle\frac{\partial(\rho_{i}h_{i}\Gamma_{i}^{2}\mathbf{u}_{i})}{\partial
t}+\nabla\cdot(\rho_{i}h_{i}\Gamma_{i}^{2}\textbf{u}_{i}\mathbf{u}_{i}^{\top}+p_{i}\mathbf{I})$
$\displaystyle=r_{i}\Gamma_{i}\rho_{i}(\mathbf{E}+\mathbf{u}_{i}\times\mathbf{B}),$
(1b) $\displaystyle\frac{\partial(\mathcal{E}_{i})}{\partial
t}+\nabla\cdot((\mathcal{E}_{i}+p_{i})\mathbf{u}_{i})$
$\displaystyle=r_{i}\Gamma_{i}\rho_{i}(\mathbf{u}_{i}\cdot\mathbf{E}),$ (1c)
$\displaystyle\frac{\partial(\rho_{e}\Gamma_{e})}{\partial
t}+\nabla\cdot(\rho_{e}\Gamma_{e}\mathbf{u}_{e})$ $\displaystyle=0,$ (1d)
$\displaystyle\
\frac{\partial(\rho_{e}h_{e}\Gamma_{e}^{2}\mathbf{u}_{e})}{\partial
t}+\nabla\cdot(\rho_{e}h_{e}\Gamma_{e}^{2}\textbf{u}_{e}\mathbf{u}_{e}^{\top}+p_{e}\mathbf{I})$
$\displaystyle=r_{e}\Gamma_{e}\rho_{e}(\mathbf{E}+\mathbf{u}_{e}\times\mathbf{B}),$
(1e) $\displaystyle\frac{\partial(\mathcal{E}_{e})}{\partial
t}+\nabla\cdot((\mathcal{E}_{e}+p_{e})\mathbf{u}_{e})$
$\displaystyle=r_{e}\Gamma_{e}\rho_{e}(\mathbf{u}_{e}\cdot\mathbf{E}),$ (1f)
$\displaystyle\dfrac{\partial\mathbf{B}}{\partial t}+\nabla\times\mathbf{E}$
$\displaystyle=0,$ (1g) $\displaystyle\dfrac{\partial\mathbf{E}}{\partial
t}-\nabla\times\mathbf{B}$ $\displaystyle=-\mathbf{j},$ (1h)
$\displaystyle\nabla\cdot\mathbf{E}$ $\displaystyle=\rho_{c},$ (1i)
$\displaystyle\nabla\cdot\mathbf{B}$ $\displaystyle=0,$ (1j)

where,
$\mathcal{E}_{\alpha}=\rho_{\alpha}h_{\alpha}\Gamma_{\alpha}^{2}-p_{\alpha}$
are the energy densities, $h_{\alpha}$ are the specific enthalpy,
$\Gamma_{\alpha}$ are the Lorentz factors, given by
$\Gamma_{\alpha}=\dfrac{1}{\sqrt{1-\mathbf{u}_{\alpha}^{2}}}$,
$\mathbf{B}=(B_{x},B_{y},B_{z})$ is the magnetic field vector, and
$\mathbf{E}=(E_{x},E_{y},E_{z})$ is the electric field vector. We have assumed
the speed of light to be unity. We use notations, $q_{\alpha}$ for the
particle charge and $m_{\alpha}$ for the particle mass, while the charge to
mass ratios is denoted by $r_{\alpha}=\dfrac{q_{\alpha}}{m_{\alpha}}$,
$\alpha\in\\{i,\ e\\}.$ The total charge density $\rho_{c}$ and the current
density vector $\mathbf{j}=(j_{x},j_{y},j_{z})$ are given by

$\rho_{c}=r_{i}\rho_{i}\Gamma_{i}+r_{e}\rho_{e}\Gamma_{e},\qquad\mathbf{j}=r_{i}\rho_{i}\Gamma_{i}\mathbf{u}_{i}+r_{e}\rho_{e}\Gamma_{e}\mathbf{u}_{e}.$

The above system of equations is closed using the equation of state,
$h_{\alpha}=h_{\alpha}(p_{\alpha},\rho_{\alpha}),$ where we use the ideal
equation of state for $h_{\alpha}$, which is given by

$h_{\alpha}=1+\frac{\gamma_{\alpha}}{\gamma_{\alpha}-1}\frac{p_{\alpha}}{\rho_{\alpha}},$

where, $\gamma_{\alpha}=c_{p_{\alpha}}/c_{V_{\alpha}}$ is the ratio of
specific heats. Accordingly, the polytropic index $n_{\alpha}$, and the sound
speed $c_{\alpha}=c_{s_{\alpha}}$ can be written as

$n_{\alpha}=k_{\alpha}-1\qquad\text{ and }\qquad
c_{\alpha}^{2}=\frac{k_{\alpha}p_{\alpha}}{n_{\alpha}\rho_{\alpha}h_{\alpha}},$

where $k_{\alpha}=\dfrac{\gamma_{\alpha}}{\gamma_{\alpha}-1}$ are constants.

Equations (1a) and (1d) are mass conservation laws for ion and electrons.
Similarly, Equations (1b) and (1e) represent momentum balance, where source
terms are due to Lorentz force acting on the fluid due to the electromagnetic
variables. The energy conservation equations for ion and electron fluids are
(1c) and (1f), respectively. Here, source terms are kinetic energy
contributions due to electric and magnetic fields. Equations (1g)-(1j) are
Maxwell’s equations, describing the constraint evolution of the electric and
magnetic fields.

Usually, to deal with the constraints (1i) and (1j) in plasma flow models
[15], a hyperbolic form of Maxwell’s equations [37] is considered. However, we
did not observe any difference in the numerical results if we consider the
Hyperbolic formulation of the equations. Furthermore, more suitable methods
are presented in [12, 29]. Hence, we have ignored these constraint equations.

##  3 Analysis of continuous problem

Let us introduce the notations, $D_{\alpha}=\Gamma_{\alpha}\rho_{\alpha}$ and
$\mathbf{M}_{\alpha}=(M_{x_{\alpha}},M_{y_{\alpha}},M_{z_{\alpha}})=\rho_{\alpha}h_{\alpha}\Gamma^{2}_{\alpha}\mathbf{u}_{\alpha}$
for $\alpha\in\\{i,e\\}$ . The vector of conservative variables $\mathbf{U}$
can be written as,

$\mathbf{U}=(\mathbf{U}_{i}^{\top},\mathbf{U}_{e}^{\top},\mathbf{U}_{m}^{\top})^{\top},$

where $\mathbf{U}_{i}=(D_{i},\mathbf{M}_{i},\mathcal{E}_{i})^{\top}$ are ion
fluid variables,
$\mathbf{U}_{e}=(D_{e},\mathbf{M}_{e},\mathcal{E}_{e})^{\top}$ are electron
fluid variables, and $\mathbf{U}_{m}=(\mathbf{B},\ \mathbf{E})^{\top}$ are
Maxwell’s variables. In two dimensions, let us also denote $x$\- and
$y$-directional fluxes as $\mathbf{f}^{x}$ and $\mathbf{f}^{y}$, respectively.
From (1), expressions for the fluxes can be written as,

$\mathbf{f}^{x}=\begin{pmatrix}\mathbf{f}^{x}_{i}(\mathbf{U}_{i})\\\
\mathbf{f}^{x}_{e}(\mathbf{U}_{e})\\\
\mathbf{f}^{x}_{m}(\mathbf{U}_{m})\end{pmatrix},\qquad\text{ and
}\qquad\mathbf{f}^{y}=\begin{pmatrix}\mathbf{f}^{y}_{i}(\mathbf{U}_{i})\\\
\mathbf{f}^{y}_{e}(\mathbf{U}_{e})\\\
\mathbf{f}^{y}_{m}(\mathbf{U}_{m})\end{pmatrix},$

where

$\displaystyle\mathbf{f}_{\alpha}^{x}=\begin{pmatrix}D_{\alpha}u_{x_{\alpha}}\\\
M_{x_{\alpha}}u_{x_{\alpha}}+p_{\alpha}\\\ M_{y_{\alpha}}u_{x_{\alpha}}\\\
M_{z_{\alpha}}u_{x_{\alpha}}\\\ M_{x_{\alpha}}\end{pmatrix},\qquad\text{ and
}\qquad\mathbf{f}_{\alpha}^{y}=\begin{pmatrix}D_{\alpha}u_{y_{\alpha}}\\\
M_{x_{\alpha}}u_{y_{\alpha}}\\\ M_{y_{\alpha}}u_{y_{\alpha}}+p_{\alpha}\\\
M_{z_{\alpha}}u_{y_{\alpha}}\\\ M_{y_{\alpha}}\end{pmatrix},$ (2) for the flow
parts and $\displaystyle\mathbf{f}_{m}^{x}=\begin{pmatrix}0\\\ -E_{z}\\\
E_{y}\\\ 0\\\ B_{z}\\\ -B_{y}\end{pmatrix},\qquad\text{ and
}\qquad\mathbf{f}_{m}^{y}=\begin{pmatrix}E_{z}\\\ 0\\\ -E_{x}\\\ -B_{z}\\\
0\\\ B_{x}\end{pmatrix},$ (3)

for the Maxwell’s part. We note that the flux for the whole system contains
three independent parts, two fluid parts modeled using the equations of
special relativistic flows, and the third part is linear Maxwell’s flux. These
three parts are coupled via source terms, which are given as,

$\mathbf{s}=\begin{pmatrix}\mathbf{s}_{i}(\mathbf{U}_{i},\mathbf{U}_{m})\\\
\mathbf{s}_{e}(\mathbf{U}_{e},\mathbf{U}_{m})\\\
\mathbf{s}_{m}(\mathbf{U}_{i},\mathbf{U}_{e})\end{pmatrix},$

where,

$\displaystyle\mathbf{s}_{\alpha}=\begin{pmatrix}0\\\
r_{\alpha}D_{\alpha}(E_{x}+u_{y_{\alpha}}B_{z}-u_{z_{\alpha}}B_{y})\\\
r_{\alpha}D_{\alpha}(E_{y}+u_{z_{\alpha}}B_{x}-u_{x_{\alpha}}B_{z})\\\
r_{\alpha}D_{\alpha}(E_{z}+u_{x_{\alpha}}B_{y}-u_{y_{\alpha}}B_{x})\\\ \
r_{\alpha}D_{\alpha}(u_{x_{\alpha}}E_{x}+u_{y_{\alpha}}E_{y}+u_{z_{\alpha}}E_{z})\end{pmatrix},\qquad\text{
and }\qquad\mathbf{s}_{m}=\begin{pmatrix}0\\\ 0\\\ 0\\\ -j_{x}\\\ -j_{y}\\\
-j_{z}\end{pmatrix},$

for $\alpha\in\\{i,e\\}$. Using the above notations, the system (1), for the
two-dimensional case can be written in the conservation form as follows:

$\frac{\partial\mathbf{U}}{\partial t}+\frac{\partial\mathbf{f}^{x}}{\partial
x}+\frac{\partial\mathbf{f}^{y}}{\partial y}=\mathbf{s}.$ (4)

Let us also introduce the vector of primitive variables
$\mathbf{W}=(\mathbf{W}_{i}^{\top},\mathbf{W}_{e}^{\top},\mathbf{W}_{m}^{\top})^{\top}$,
where
$\mathbf{W}_{\alpha}=(\rho_{\alpha},\mathbf{u}_{\alpha},{p}_{\alpha})^{\top}$
and $\mathbf{W}_{m}=(\mathbf{B,E})^{\top}$. We follow the procedure given in
[38, 1] to extract the primitive variables from the conservative variables. We
define the set of admissible solution space

$\Omega=\\{\mathbf{U}\in\mathbb{R}^{16}:\ \rho_{i}>0,\rho_{e}>0,\
p_{i}>0,p_{e}>0,\ |\mathbf{u}_{i}|<1,|\mathbf{u}_{e}|<1\\}.$

The eigenvalues of the system in $x$-direction are,

$\displaystyle\Lambda^{x}=\biggl{\\{}$
$\displaystyle\frac{(1-c_{i}^{2})u_{x_{i}}-(c_{i}/\Gamma_{i})\sqrt{Q_{i}^{x}}}{1-c_{i}^{2}|\mathbf{u}_{i}|^{2}},\
u_{x_{i}},\ u_{x_{i}},\
u_{x_{i}},\frac{(1-c_{i}^{2})u_{x_{i}}+(c_{i}/\Gamma_{i})\sqrt{Q_{i}^{x}}}{1-c_{i}^{2}|\mathbf{u}_{i}|^{2}},$
$\displaystyle\frac{(1-c_{e}^{2})u_{x_{e}}-(c_{e}/\Gamma_{e})\sqrt{Q_{e}^{x}}}{1-c_{e}^{2}|\mathbf{u}_{e}|^{2}},\
u_{x_{e}},\ u_{x_{e}},\
u_{x_{e}},\frac{(1-c_{e}^{2})u_{x_{e}}+(c_{e}/\Gamma_{e})\sqrt{Q_{e}^{x}}}{1-c_{e}^{2}|\mathbf{u}_{e}|^{2}},$
$\displaystyle-1,\ -1,\ 0,\ 0,\ 1,\ 1\ \biggr{\\}},$

where,
$Q_{\alpha}^{x}=1-u_{x_{\alpha}}^{2}-c_{\alpha}^{2}(u_{y_{\alpha}}^{2}+u_{z_{\alpha}}^{2}),\
\alpha\in\\{i,\ e\\}$. For $\mathbf{U}\in\Omega$ and $\gamma_{\alpha}\in(1,2]$
we have $c_{\alpha}<1$ which implies that $Q_{\alpha}^{x}>0$, hence all the
eigenvalues are real. The complete set of right eigenvectors for the Jacobian
matrix $\dfrac{\partial\mathbf{f}^{x}}{\partial\mathbf{U}}$ is given in
Appendix (A.1). The Appendix also contains the eigenvalues and right
eigenvectors for the jacobian matrix
$\dfrac{\partial\mathbf{f}^{y}}{\partial\mathbf{U}}$. Consequently, we can
state the following Lemma.

######  Lemma 3.1.

The system (4) is hyperbolic for the states $\mathbf{U}\in\Omega$, with real
eigenvalues and a complete set of eigenvectors.

We now introduce entropy functions for the two-fluid relativistic plasma flow
equations. The entropy functions $\mathcal{U}_{\alpha}$ and associated entropy
fluxes $\mathcal{F}^{d}_{\alpha}$ for the fluid part of the system (4) are
given by,

$\mathcal{U}_{\alpha}=-\frac{\rho_{\alpha}\Gamma_{\alpha}s_{\alpha}}{\gamma_{\alpha}-1}\qquad\text{
and
}\qquad\mathcal{F}^{d}_{\alpha}=-\frac{\rho_{\alpha}\Gamma_{\alpha}s_{\alpha}u_{d_{\alpha}}}{\gamma_{\alpha}-1},\qquad\alpha\in\\{i,e\\},\quad
d=x,y.$

Here, $s_{\alpha}=\ln(p_{\alpha}\rho_{\alpha}^{-\gamma_{\alpha}})$. The pair
$(\mathcal{U}_{\alpha},\mathcal{F}^{d}_{\alpha})$ is called an entropy-entropy
flux pair. For simplicity, we only consider the one-dimensional case, i.e.,

$\frac{\partial\mathbf{U}}{\partial t}+\frac{\partial\mathbf{f}^{x}}{\partial
x}=\mathbf{s}.$ (5)

The extension of entropy framework to higher dimensions is straightforward. We
will now prove the following result:

######  Proposition 3.1.

The smooth solutions of (5) satisfy the entropy equality,

$\partial_{t}s_{\alpha}+u_{x_{\alpha}}\partial_{x}s_{\alpha}=0,\ \
\alpha\in\\{i,e\\}.$

As a consequence, every smooth function $H(s_{\alpha})$ of $s_{\alpha}$
satisfies

$\partial_{t}(\rho_{\alpha}\Gamma_{\alpha}H(s_{\alpha}))+\partial_{x}(\rho_{\alpha}\Gamma_{\alpha}u_{x_{\alpha}}H(s_{\alpha}))=0.$
(6)

In particular, we have the following entropy equality for the smooth
solutions,

$\partial_{t}\mathcal{U}_{\alpha}+\partial_{x}\mathcal{F}^{x}_{\alpha}=0.$ (7)

We will give the proof of Proposition (3.1) in two steps. Firstly, we will
prove two lemmas, and then we give the proof of proposition using these two
lemmas. The proof is similar to those in [1] for relativistic hydrodynamics
except for the additional source terms and three-dimensional velocity field.

######  Lemma 3.2.

The smooth solutions of the system (5) satisfies the following identity

$\displaystyle\frac{1}{\Gamma_{\alpha}^{2}}\partial_{t}p_{\alpha}=\partial_{t}p_{\alpha}+u_{x\alpha}\partial_{x}p_{\alpha}+\left(\frac{\rho_{\alpha}h_{\alpha}}{\Gamma_{\alpha}}\right)(\partial_{t}\Gamma_{\alpha}+\partial_{x}{(\Gamma_{\alpha}u_{x\alpha})})-\rho_{\alpha}h_{\alpha}\partial_{x}u_{x,\alpha}-\mathbf{u}_{\alpha}\cdot(\mathbf{s}_{\mathbf{M}_{\alpha}}-\mathbf{u}_{\alpha}\mathbf{s}_{\mathcal{E}_{\alpha}}),$
(8)

where

$\mathbf{s}_{\mathbf{M}_{\alpha}}=r_{\alpha}\Gamma_{\alpha}\rho_{\alpha}(\mathbf{E}+\mathbf{u}_{\alpha}\times\mathbf{B}),\qquad\mathbf{s}_{\mathcal{E}_{\alpha}}=r_{\alpha}\Gamma_{\alpha}\rho_{\alpha}(\mathbf{u}_{\alpha}\cdot\mathbf{E}).$

are the source terms in the momentum and energy equations.

###### Proof.

Let $\mathbf{U}$ be the smooth solution of the system (5). Then we can write,

$\displaystyle\partial_{t}(\rho_{\alpha}\Gamma_{\alpha})+\partial_{x}(\rho_{\alpha}\Gamma_{\alpha}u_{x_{\alpha}})=0,$
(9)
$\displaystyle\partial_{t}(\rho_{\alpha}h_{\alpha}\Gamma_{\alpha}^{2}u_{x_{\alpha}})+\partial_{x}(\rho_{\alpha}h_{\alpha}\Gamma_{\alpha}^{2}u_{x_{\alpha}}^{2})+\partial_{x}p_{\alpha}=\mathbf{s}_{M_{x_{\alpha}}},$
(10)
$\displaystyle\partial_{t}(\rho_{\alpha}h_{\alpha}\Gamma_{\alpha}^{2}u_{y_{\alpha}})+\partial_{x}(\rho_{\alpha}h_{\alpha}\Gamma_{\alpha}^{2}u_{x_{\alpha}}u_{y_{\alpha}})=\mathbf{s}_{M_{y_{\alpha}}},$
(11)
$\displaystyle\partial_{t}(\rho_{\alpha}h_{\alpha}\Gamma_{\alpha}^{2}u_{z_{\alpha}})+\partial_{x}(\rho_{\alpha}h_{\alpha}\Gamma_{\alpha}^{2}u_{x_{\alpha}}u_{z_{\alpha}})=\mathbf{s}_{M_{z_{\alpha}}},$
(12)
$\displaystyle\partial_{t}(\rho_{\alpha}h_{\alpha}\Gamma_{\alpha}^{2})-\partial_{t}(p_{\alpha})+\partial_{x}(\rho_{\alpha}h_{\alpha}\Gamma_{\alpha}^{2}u_{x_{\alpha}})=\mathbf{s}_{\mathcal{E}_{\alpha}},$
(13)

where, we have used
$\mathbf{s}_{\mathbf{M}_{\alpha}}=(\mathbf{s}_{M_{x_{\alpha}}},\mathbf{s}_{M_{y_{\alpha}}},\mathbf{s}_{M_{z_{\alpha}}})$,
as the momentum source vector, i.e., source terms of the Eqns. (1b) and (1e).
Also, $\mathbf{s}_{\mathcal{E}_{\alpha}}$ is the source term of the energy
equations (1c), and (1f). We apply the product rule on Equations (10)-(12) and
substitute the value of
$\partial_{x}(\rho_{\alpha}h_{\alpha}\Gamma^{2}_{\alpha}u_{x_{\alpha}})$ from
Eqn. (13) in each of the three equations to obtain the following set of three
identities:

$\displaystyle(\rho_{\alpha}h_{\alpha}\Gamma^{2}_{\alpha}u_{x_{\alpha}})(\partial_{t}u_{x_{\alpha}}+u_{x_{\alpha}}\partial_{x}u_{x_{\alpha}})+u_{x_{\alpha}}^{2}\partial_{t}p_{\alpha}+u_{x_{\alpha}}\partial_{x}p_{\alpha}=u_{x_{\alpha}}(\mathbf{s}_{M_{x_{\alpha}}}-u_{x_{\alpha}}\mathbf{s}_{\mathcal{E}_{\alpha}})$
(14a)
$\displaystyle(\rho_{\alpha}h_{\alpha}\Gamma^{2}_{\alpha}u_{y_{\alpha}})(\partial_{t}u_{y_{\alpha}}+u_{y_{\alpha}}\partial_{x}u_{y_{\alpha}})+u_{y_{\alpha}}^{2}\partial_{t}p_{\alpha}=u_{y_{\alpha}}(\mathbf{s}_{M_{y_{\alpha}}}-u_{y_{\alpha}}\mathbf{s}_{\mathcal{E}_{\alpha}})$
(14b)
$\displaystyle(\rho_{\alpha}h_{\alpha}\Gamma^{2}_{\alpha}u_{z_{\alpha}})(\partial_{t}u_{z_{\alpha}}+u_{z_{\alpha}}\partial_{x}u_{z_{\alpha}})+u_{z_{\alpha}}^{2}\partial_{t}p_{\alpha}=u_{z_{\alpha}}(\mathbf{s}_{M_{z_{\alpha}}}-u_{z_{\alpha}}\mathbf{s}_{\mathcal{E}_{\alpha}})$
(14c)

We observe that
$\partial_{t}\Gamma_{\alpha}=\Gamma^{3}_{\alpha}\boldsymbol{u}_{\alpha}\cdot\partial_{t}\boldsymbol{u}_{\alpha}$,
and
$\partial_{x}\Gamma_{\alpha}=\Gamma^{3}_{\alpha}\boldsymbol{u}_{\alpha}\cdot\partial_{x}\boldsymbol{u}_{\alpha}.$
Consequently, Eqn. (14a) simplifies to

$\displaystyle\left(\frac{\rho_{\alpha}h_{\alpha}}{\Gamma_{\alpha}}\right)(\partial_{t}\Gamma_{\alpha}+u_{x_{\alpha}}\partial_{x}\Gamma_{\alpha})-\rho_{\alpha}h_{\alpha}\Gamma_{\alpha}^{2}u_{y_{\alpha}}(\partial_{t}u_{y_{\alpha}}+u_{x_{\alpha}}\partial_{x}u_{y_{\alpha}})-\rho_{\alpha}h_{\alpha}\Gamma_{\alpha}^{2}u_{z_{\alpha}}(\partial_{t}u_{z_{\alpha}}+u_{x_{\alpha}}\partial_{x}u_{z_{\alpha}})$
$\displaystyle+u_{x_{\alpha}}^{2}\partial_{t}p_{\alpha}+u_{x_{\alpha}}\partial_{x}p_{\alpha}=u_{x_{\alpha}}(\mathbf{s}_{M_{x_{\alpha}}}-u_{x_{\alpha}}\mathbf{s}_{\mathcal{E}_{\alpha}}),$
(15)

From the definition of Lorentz factor $\Gamma_{\alpha}$, we have
$u_{x_{\alpha}}^{2}=1-\dfrac{1}{\Gamma_{\alpha}^{2}}-u_{y_{\alpha}}^{2}-u_{z_{\alpha}}^{2}$.
We substitute this value of $u_{x_{\alpha}}^{2}$ in Eqn. (15), followed by the
substitution of value of $u_{y_{\alpha}}^{2}\partial_{t}p_{\alpha}$ and
$u_{z_{\alpha}}^{2}\partial_{t}p_{\alpha}$ from the Eqn. (14b)-(14c), to
obtain the required identity (8). ∎

######  Lemma 3.3.

For smooth solutions of the system (5), we have the following identity

$(\partial_{t}p_{\alpha}+u_{x_{\alpha}}\partial_{x}p_{\alpha})+\frac{p_{\alpha}\gamma_{\alpha}}{\Gamma_{\alpha}}(\partial_{t}\Gamma_{\alpha}+\partial_{x}{(\Gamma_{\alpha}u_{x_{\alpha}})})=0.$
(16)

###### Proof.

We expand Eqn. (13) using the product rule and simplify further using the
substitution of variables from Eqn. (9) to obtain,

$\displaystyle\partial_{t}p_{\alpha}+u_{x_{\alpha}}\partial_{x}p_{\alpha}=\left(\dfrac{\gamma_{\alpha}-1}{\Gamma_{\alpha}^{2}\gamma_{\alpha}}\right)\mathbf{s}_{\mathcal{E}_{\alpha}}+\left(\frac{\gamma_{\alpha}-1}{\Gamma_{\alpha}^{2}\gamma_{\alpha}}\right)\partial_{t}p_{\alpha}-\left(\frac{\rho_{\alpha}}{\Gamma_{\alpha}}\frac{\gamma_{\alpha}-1}{\gamma_{\alpha}}\right)(\partial_{t}\Gamma_{\alpha}+\partial_{x}{(\Gamma_{\alpha}u_{x_{\alpha}})})$
$\displaystyle-2(\partial_{t}\Gamma_{\alpha}+\partial_{x}{(\Gamma_{\alpha}u_{x_{\alpha}})})\frac{p_{\alpha}}{\Gamma_{\alpha}}+\rho_{\alpha}\left(\frac{\gamma_{\alpha}-1}{\gamma_{\alpha}}\right)\partial_{x}u_{x_{\alpha}}+p_{\alpha}\partial_{x}u_{x_{\alpha}}.$
(17)

Substituting the value of
$\dfrac{1}{\Gamma_{\alpha}^{2}}\partial_{t}p_{\alpha}$ from the identity of
Lemma (3.2), and simplifying further we get

$\displaystyle\partial_{t}p_{\alpha}+u_{x_{\alpha}}\partial_{x}p_{\alpha}$
$\displaystyle=\left(\frac{\gamma_{\alpha}-1}{\gamma_{\alpha}}\right)\left(-\mathbf{u}_{\alpha}\cdot\mathbf{s}_{\mathbf{M}_{\alpha}}+\left(\frac{1+\mathbf{u}_{\alpha}^{2}}{\Gamma_{\alpha}^{2}}\right)\mathbf{s}_{\mathcal{E}_{\alpha}}\right)$
$\displaystyle+\left(\frac{\gamma_{\alpha}-1}{\gamma_{\alpha}}\right)(\partial_{t}p_{\alpha}+u_{x_{\alpha}}\partial_{x}p_{\alpha})$
$\displaystyle+\left(\frac{\gamma_{\alpha}-1}{\gamma_{\alpha}}\right)\left[\left(\frac{\rho_{\alpha}h_{\alpha}}{\Gamma_{\alpha}}\right)(\partial_{t}\Gamma_{\alpha}+\partial_{x}(\Gamma_{\alpha}u_{x_{\alpha}}))-\rho_{\alpha}h_{\alpha}\partial_{x}u_{x_{\alpha}}\right]$
$\displaystyle-\left(\frac{\rho_{\alpha}}{\Gamma_{\alpha}}\frac{\gamma_{\alpha}-1}{\gamma_{\alpha}}\right)(\partial_{t}\Gamma_{\alpha}+\partial_{x}{(\Gamma_{\alpha}u_{x_{\alpha}})})-\frac{2p_{\alpha}}{\Gamma_{\alpha}}(\partial_{t}\Gamma_{\alpha}+\partial_{x}{(\Gamma_{\alpha}u_{x_{\alpha}})})$
$\displaystyle+\left(\frac{\gamma_{\alpha}-1}{\gamma_{\alpha}}\right)\rho_{\alpha}\partial_{x}u_{x_{\alpha}}+p_{\alpha}\partial_{x}u_{x_{\alpha}},$

We observe that
$\mathbf{u}_{\alpha}\cdot\mathbf{s_{M}}_{\alpha}-\mathbf{s}_{\mathcal{E}_{\alpha}}=0$,
thus, after further simplifications and cancellations [1], we obtain (16). ∎

###### Proof of Proposition (3.1).

We find the partial derivatives of
$s_{\alpha}=\ln(p_{\alpha}\rho_{\alpha}^{-\gamma_{\alpha}})$ with respect to
$t$ and $x$, to obtain

$\partial_{t}s_{\alpha}=\frac{1}{p_{\alpha}}\partial_{t}p_{\alpha}-\frac{\gamma_{\alpha}}{\rho_{\alpha}}\partial_{t}\rho_{\alpha}\quad\text{
and
}\quad\partial_{x}s_{\alpha}=\frac{1}{p_{\alpha}}\partial_{x}p_{\alpha}-\frac{\gamma_{\alpha}}{\rho_{\alpha}}\partial_{x}\rho_{\alpha}.$
(18)

Once we have the partial derivatives, we apply the product rule on the mass
density Eqn. (9) to obtain an identity in terms of the partial derivatives of
$\rho_{\alpha}$

$\frac{1}{\rho_{\alpha}}(\partial_{t}\rho_{\alpha}+u_{x_{\alpha}}\partial_{x}\rho_{\alpha})=-\frac{1}{\Gamma_{\alpha}}(\partial_{t}\Gamma_{\alpha}+\partial_{x}{(\Gamma_{\alpha}u_{x_{\alpha}})}).$
(19)

Combining Eqn. (16) and Eqn. (19) we obtain the expression

$\frac{1}{p_{\alpha}}(\partial_{t}p_{\alpha}+u_{x_{\alpha}}\partial_{x}p_{\alpha})-\frac{\gamma_{\alpha}}{\rho_{\alpha}}(\partial_{t}\rho_{\alpha}+u_{x_{\alpha}}\partial_{x}\rho_{\alpha})=0,$

which, on using Eqn. (18), is equivalent to
$\partial_{t}s_{\alpha}+u_{x_{\alpha}}\partial_{x}s_{\alpha}=0$. The identity
(6) is a direct consequence of the product rule. To obtain the final entropy
equality (7), we simplify the identity (6) by choosing particular value of
$H(s_{\alpha})$ as $H(s_{\alpha})=\dfrac{-s_{\alpha}}{\gamma_{\alpha}-1}$. ∎

######  Remark 3.1.

The entropy equality (7) in Proposition (3.1) gives rise to the entropy
inequality

$\partial_{t}\mathcal{U}_{\alpha}+\partial_{x}\mathcal{F}^{x}_{\alpha}\leq 0,$
(20)

in the sense of distributions for non-smooth solutions.

In the above analysis, we also note that the Maxwell’s flux does not affect
the fluid entropies. This fact will be exploited to design semi-discrete high-
order numerical schemes which will satisfy the entropy inequality (20).

##  4 Semi-discrete entropy stable schemes

In this section, we will design semi-discrete entropy stable numerical schemes
for the two-dimensional two-fluid relativistic plasma flow Eqns. (4). The
extension to three dimensions is fairly easy. Let us consider the domain
$D=I_{x}\times I_{y}$, where $I_{x}=(x_{min},x_{max})$ and
$I_{y}=(y_{min},y_{max})$. The domain $D$ is discretized using the uniform
mesh with cells of size $\Delta x\times\Delta y$. We define
$x_{i}=x_{min}+i\Delta x$ for $0\leq i\leq N_{x}$ and $y_{j}=y_{min}+j\Delta
y$ for $0\leq j\leq N_{y}$. The cell $(i,j)$ is given by
$I_{i,j}=[x_{i-1/2},x_{i+1/2}]\times[y_{j-1/2},y_{j+1/2}]$, where
$x_{i+1/2}=\frac{x_{i}+x_{i+1}}{2}$ and $y_{j+1/2}=\frac{y_{j}+y_{j+1}}{2}$.

The semi-discrete finite difference scheme for the evolution of solution
$\mathbf{U}_{i,j}$ at $(x_{i},y_{j})$ of system (4) is given by,

$\frac{d}{dt}\mathbf{U}_{i,j}(t)+\frac{1}{\Delta
x}\left(\mathbf{F}_{i+\frac{1}{2},j}^{x}(t)-\mathbf{F}_{i-\frac{1}{2},j}^{x}(t)\right)+\frac{1}{\Delta
y}\left(\mathbf{F}_{i,j+\frac{1}{2}}^{y}(t)-\mathbf{F}_{i,j-\frac{1}{2}}^{y}(t)\right)=\mathbf{s}(\mathbf{U}_{i,j}(t)),$
(21)

where $\mathbf{F}_{i+\frac{1}{2},j}^{x}$, $\mathbf{F}_{i,j+\frac{1}{2}}^{y}$
are the numerical fluxes consistent with the continuous fluxes
$\mathbf{f}^{x}$, $\mathbf{f}^{y}$, respectively.

###  4.1 Entropy stable numerical schemes for the fluid equations

To simplify the discussion, we will first present the discretization of the
fluid part, i.e., we consider the following differential form:

$\frac{\partial\mathbf{U}_{\alpha}}{\partial
t}+\frac{\partial\mathbf{f}^{x}_{\alpha}}{\partial
x}+\frac{\partial\mathbf{f}^{y}_{\alpha}}{\partial
y}=\mathbf{s}_{\alpha}(\mathbf{U}_{\alpha},\mathbf{U}_{m}),$ (22)

for $\alpha\in\\{i,e\\}$. Let us now define:

######  Definition 4.1.

The numerical scheme

$\displaystyle\frac{d}{dt}\mathbf{U}_{\alpha,i,j}(t)+\frac{1}{\Delta
x}\left(\mathbf{F}_{\alpha,i+\frac{1}{2},j}^{x}(t)-\mathbf{F}_{\alpha,i-\frac{1}{2},j}^{x}(t)\right)+\frac{1}{\Delta
y}\left(\mathbf{F}_{\alpha,i,j+\frac{1}{2}}^{y}(t)-\mathbf{F}_{\alpha,i,j-\frac{1}{2}}^{y}(t)\right)$
$\displaystyle=\mathbf{s}_{\alpha}(\mathbf{U}_{\alpha,i,j}(t),\mathbf{U}_{m,i,j}(t)),$
(23)

is said to be entropy conservative if the following entropy equality is
satisfied

$\frac{d}{dt}\mathcal{U}_{\alpha}(\mathbf{U}_{ij})+\frac{1}{\Delta
x}\left(\hat{\mathcal{F}}_{\alpha,i+\frac{1}{2},j}^{x}-\hat{\mathcal{F}}_{\alpha,i-\frac{1}{2},j}^{x}\right)+\frac{1}{\Delta
y}\left(\hat{\mathcal{F}}_{\alpha,i,j+\frac{1}{2}}^{y}-\hat{\mathcal{F}}_{\alpha,i,j-\frac{1}{2}}^{y}\right)=0,\
\ \ \alpha\in\\{i,e\\},$ (24)

where $\hat{\mathcal{F}}_{\alpha,i+\frac{1}{2},j}^{x}$,
$\hat{\mathcal{F}}_{\alpha,i,j+\frac{1}{2}}^{y}$ are some numerical entropy
flux functions consistent with the entropy fluxes
${\mathcal{F}}_{\alpha}^{x}$, ${\mathcal{F}}_{\alpha}^{y}$, respectively. Here
$\mathbf{F}_{\alpha,i+\frac{1}{2},j}^{x}$,
$\mathbf{F}_{\alpha,i,j+\frac{1}{2}}^{y}$ are the numerical fluxes consistent
with $\mathbf{f}^{x}$, $\mathbf{f}^{y}$, respectively.

The scheme is said to be entropy stable if the following entropy inequality is
satisfied:

$\frac{d}{dt}\mathcal{U}_{\alpha}(\mathbf{U}_{\alpha,i,j})+\frac{1}{\Delta
x}\left(\hat{\mathcal{F}}_{\alpha,i+\frac{1}{2},j}^{x}-\hat{\mathcal{F}}_{\alpha,i-\frac{1}{2},j}^{x}\right)+\frac{1}{\Delta
y}\left(\hat{\mathcal{F}}_{\alpha,i,j+\frac{1}{2}}^{y}-\hat{\mathcal{F}}_{\alpha,i,j-\frac{1}{2}}^{y}\right)\leq
0\ \ \ \alpha\in\\{i,e\\},$ (25)

To obtain the entropy stable schemes, we will first construct the entropy
conservative schemes and then add dissipative fluxes to obtain entropy
inequality.

####  4.1.1 High-order entropy conservative schemes for fluid equations

Let us introduce entropy variables,

$\mathbf{V}_{\alpha}(\mathbf{U}_{\alpha})=\frac{\partial\mathcal{U}_{\alpha}}{\partial\mathbf{U}_{\alpha}}$

and entropy potentials,

$\psi_{\alpha}^{x}(\mathbf{U}_{\alpha})=\mathbf{V}_{\alpha}^{\top}(\mathbf{U}_{\alpha})\cdot\mathbf{f}_{\alpha}^{x}(\mathbf{U}_{\alpha})-\mathcal{F}_{\alpha}^{x}(\mathbf{U}_{\alpha}),\qquad\psi_{\alpha}^{y}(\mathbf{U}_{\alpha})=\mathbf{V}_{\alpha}^{\top}(\mathbf{U}_{\alpha})\cdot\mathbf{f}_{\alpha}^{y}(\mathbf{U}_{\alpha})-\mathcal{F}_{\alpha}^{y}(\mathbf{U}_{\alpha})$

for $\alpha\in\\{i,e\\}$. A simple calculation results in the following
expressions

$\mathbf{V}_{\alpha}=\begin{pmatrix}\dfrac{\gamma_{\alpha}-s_{\alpha}}{\gamma_{\alpha}-1}+{\beta_{\alpha}}\\\
{u_{x_{\alpha}}\Gamma_{\alpha}\beta_{\alpha}}\\\
{u_{y_{\alpha}}\Gamma_{\alpha}\beta_{\alpha}}\\\
{u_{z_{\alpha}}\Gamma_{\alpha}\beta_{\alpha}}\\\
-{\Gamma_{\alpha}\beta_{\alpha}}\end{pmatrix},\qquad\text{with
}\qquad\beta_{\alpha}=\dfrac{\rho_{\alpha}}{p_{\alpha}},$ (26)

and

${\psi}_{\alpha}^{x}=\rho_{\alpha}\Gamma_{\alpha}u_{x_{\alpha}},\qquad{\psi}_{\alpha}^{y}=\rho_{\alpha}\Gamma_{\alpha}u_{y_{\alpha}}.$
(27)

We also introduce the following notations for jumps and averaging operations
over the cell interfaces

$[\\![a]\\!]_{i+\frac{1}{2},j}=a_{i+1,j}-a_{i,j},\qquad\bar{a}_{i+\frac{1}{2},j}=\frac{1}{2}(a_{i+1,j}+a_{i,j}),$
$[\\![a]\\!]_{i,j+\frac{1}{2}}=a_{i,j+1}-a_{i,j},\qquad\bar{a}_{i,j+\frac{1}{2}}=\frac{1}{2}(a_{i,j+1}+a_{i,j}).$

Let us first consider the homogeneous fluid part, i.e., we consider

$\frac{\partial\mathbf{U}_{\alpha}}{\partial
t}+\frac{\partial\mathbf{f}^{x}_{\alpha}}{\partial
x}+\frac{\partial\mathbf{f}^{y}_{\alpha}}{\partial y}=0.$ (28)

and the corresponding semi-discrete scheme

$\displaystyle\frac{d}{dt}\mathbf{U}_{\alpha,i,j}(t)+\frac{1}{\Delta
x}\left(\mathbf{F}_{\alpha,i+\frac{1}{2},j}^{x}(t)-\mathbf{F}_{\alpha,i-\frac{1}{2},j}^{x}(t)\right)+\frac{1}{\Delta
y}\left(\mathbf{F}_{\alpha,i,j+\frac{1}{2}}^{y}(t)-\mathbf{F}_{\alpha,i,j-\frac{1}{2}}^{y}(t)\right)=0$
(29)

Then we have the following result from [39].

######  Theorem 4.1 (Tadmor [39]).

Let $\tilde{\mathbf{F}}^{x}_{\alpha}$ and $\tilde{\mathbf{F}}^{y}_{\alpha}$ be
the consistent numerical fluxes which satisfy

$[\\![\mathbf{V}_{\alpha}]\\!]^{\top}_{i+\frac{1}{2},j}\,\tilde{\mathbf{F}}_{\alpha,i+\frac{1}{2},j}^{x}=[\\![\psi^{x}_{\alpha}]\\!]_{i+\frac{1}{2},j},\
\quad\
[\\![\mathbf{V}_{\alpha}]\\!]^{\top}_{i,j+\frac{1}{2}}\,\tilde{\mathbf{F}}_{\alpha,i,j+\frac{1}{2}}^{y}=[\\![\psi^{y}_{\alpha}]\\!]_{i,j+\frac{1}{2}},$
(30)

then the scheme (29) with the numerical fluxes
$\tilde{\mathbf{F}}^{x}_{\alpha}$ and $\tilde{\mathbf{F}}^{y}_{\alpha}$ is
second order accurate and entropy conservative, i.e., the entropy equality

$\frac{d}{dt}\mathcal{U}_{\alpha}(\mathbf{U}_{ij})+\frac{1}{\Delta
x}\left(\tilde{\mathcal{F}}_{\alpha,i+\frac{1}{2},j}^{x}-\tilde{\mathcal{F}}_{\alpha,i-\frac{1}{2},j}^{x}\right)+\frac{1}{\Delta
y}\left(\tilde{\mathcal{F}}_{\alpha,i,j+\frac{1}{2}}^{y}-\tilde{\mathcal{F}}_{\alpha,i,j-\frac{1}{2}}^{y}\right)=0,$

is satisfied with the consistent entropy numerical fluxes,

$\tilde{\mathcal{F}}_{\alpha,i+\frac{1}{2},j}^{x}=\bar{\mathbf{V}}_{\alpha,i+\frac{1}{2},j}^{\top}\tilde{\mathbf{F}}_{\alpha,i+\frac{1}{2},j}^{x}-\bar{\psi}_{\alpha,i+\frac{1}{2},j}^{x},\qquad\text{
and
}\qquad\tilde{\mathcal{F}}_{\alpha,i,j+\frac{1}{2}}^{y}=\bar{\mathbf{V}}_{\alpha,i,j+\frac{1}{2}}^{\top}\tilde{\mathbf{F}}_{\alpha,i,j+\frac{1}{2}}^{y}-\bar{\psi}_{\alpha,i,j+\frac{1}{2}}^{y}.$

We consider the $x$-directional identity of Eqn. (30) and observe that the
equation (30) contains five unknowns,
$\tilde{\mathbf{F}}^{x}_{\alpha}=(F^{x}_{\alpha,1},\ F^{x}_{\alpha,2},\
F^{x}_{\alpha,3},\ F^{x}_{\alpha,4},\ F^{x}_{\alpha,5})^{\top}$. Hence, in
general, we will have non-uniqueness for the solutions of this algebraic
equation. Several authors [40, 32] have presented different techniques to find
an affordable entropy conservative flux. We will follow the procedure of [32].

We will follow [1] to derive the following expression for the fluxes. Consider
the flux at $(i+\frac{1}{2},j)$. To simplify the notation, we suppress the
cell indices $i,j$, and let

$[\\![\cdot]\\!]=(\cdot)_{i+1,j}-(\cdot)_{i,j},\qquad(\bar{\cdot})=\frac{(\cdot)_{i,j}+(\cdot)_{i+1,j}}{2}$

Define $a^{\ln}=\frac{[\\![a]\\!]}{[\\![\log a]\\!]}$ as the logarithmic
average for a strictly positive scalar $a$.

Applying the jump condition
$[\\![ab]\\!]=\bar{a}[\\![b]\\!]+\bar{b}[\\![a]\\!]$ on the equation
$\Gamma_{\alpha}=\dfrac{1}{\sqrt{1-\mathbf{u}^{2}_{\alpha}}}$ written as
$\Gamma_{\alpha}^{2}=1+m_{x_{\alpha}}^{2}+m_{y_{\alpha}}^{2}+m_{z_{\alpha}}^{2}$,
we obtain

$[\\![\Gamma_{\alpha}]\\!]=\frac{1}{\overline{\Gamma_{\alpha}}}(\overline{m_{x_{\alpha}}}[\\![m_{x_{\alpha}}]\\!]+\overline{m_{y_{\alpha}}}[\\![m_{y_{\alpha}}]\\!]+\overline{m_{z_{\alpha}}}[\\![m_{z_{\alpha}}]\\!]),$

where, $m_{x_{\alpha}}=\Gamma_{\alpha}u_{x_{\alpha}}$,
$m_{y_{\alpha}}=\Gamma_{\alpha}u_{y_{\alpha}}$ and
$m_{z_{\alpha}}=\Gamma_{\alpha}u_{z_{\alpha}}$. Using this, we can write the
jump in $\mathbf{V}_{\alpha}$ in terms of the jump in $\rho_{\alpha}$,
$\beta_{\alpha}$, $m_{x_{\alpha}}$, $m_{y_{\alpha}}$ and $m_{z_{\alpha}}$ as
follows:

${[\\![\mathbf{V}_{\alpha}]\\!]}=\begin{pmatrix}\frac{[\\![\rho_{\alpha}]\\!]}{{\rho_{\alpha}}^{\ln}}+k_{\alpha}[\\![\beta_{\alpha}]\\!]\\\
\overline{m_{x_{\alpha}}}[\\![\beta_{\alpha}]\\!]+\bar{\beta_{\alpha}}[\\![m_{x_{\alpha}}]\\!]\\\
\overline{m_{y_{\alpha}}}[\\![\beta_{\alpha}]\\!]+\bar{\beta_{\alpha}}[\\![m_{y_{\alpha}}]\\!]\\\
\overline{m_{z_{\alpha}}}[\\![\beta_{\alpha}]\\!]+\bar{\beta_{\alpha}}[\\![m_{z_{\alpha}}]\\!]\\\
-\bar{\Gamma_{\alpha}}[\\![\beta_{\alpha}]\\!]-\frac{\bar{\beta_{\alpha}}\overline{m_{x_{\alpha}}}}{\bar{\Gamma_{\alpha}}}[\\![m_{x_{\alpha}}]\\!]-\frac{\bar{\beta_{\alpha}}\overline{m_{y_{\alpha}}}}{\bar{\Gamma_{\alpha}}}[\\![m_{y_{\alpha}}]\\!]-\frac{\bar{\beta_{\alpha}}\overline{m_{z_{\alpha}}}}{\bar{\Gamma_{\alpha}}}[\\![m_{z_{\alpha}}]\\!]\end{pmatrix},$

where $\beta_{\alpha}=\frac{\rho_{\alpha}}{p_{\alpha}}$, and
$k_{\alpha}=\left(\dfrac{1}{\gamma-1}\dfrac{1}{{\beta^{\ln}_{\alpha}}}+1\right)$.
Furthermore, $[\\![\psi^{x}_{\alpha}]\\!]$ can be written as

$[\\![\psi^{x}_{\alpha}]\\!]=\overline{\rho_{\alpha}}[\\![m_{x_{\alpha}}]\\!]+\overline{m_{x_{\alpha}}}[\\![\rho_{\alpha}]\\!].$

and Eqn. (30) becomes

$F_{\alpha,1}^{x}[\\![V_{\alpha,1}]\\!]+F_{\alpha,2}^{x}[\\![V_{\alpha,2}]\\!]+F_{\alpha,3}^{x}[\\![V_{\alpha,3}]\\!]+F_{\alpha,4}^{x}[\\![V_{\alpha,4}]\\!]+F_{\alpha,5}^{x}[\\![V_{\alpha,5}]\\!]=[\\![\psi^{x}_{\alpha}]\\!],$

which simplifies to,

$\displaystyle\left(\frac{F_{\alpha,1}^{x}}{{\rho_{\alpha}}^{\ln}}\right)[\\![\rho_{\alpha}]\\!]+\left(\overline{\beta_{\alpha}}F_{\alpha,2}^{x}-\frac{\bar{\beta_{\alpha}}\overline{m_{x_{\alpha}}}}{\bar{\Gamma_{\alpha}}}F_{\alpha,5}^{x}\right)[\\![m_{x_{\alpha}}]\\!]+\left(\overline{\beta_{\alpha}}F_{\alpha,3}^{x}-\frac{\bar{\beta_{\alpha}}\overline{m_{y_{\alpha}}}}{\bar{\Gamma_{\alpha}}}F_{\alpha,5}^{x}\right)[\\![m_{y_{\alpha}}]\\!]$
$\displaystyle+\left(\overline{\beta_{\alpha}}F_{\alpha,4}^{x}-\frac{\bar{\beta_{\alpha}}\overline{m_{z_{\alpha}}}}{\bar{\Gamma_{\alpha}}}F_{\alpha,5}^{x}\right)[\\![m_{z_{\alpha}}]\\!]+\bigg{(}k_{\alpha}F_{\alpha,1}^{x}+\overline{m_{x_{\alpha}}}F_{\alpha,2}^{x}+\overline{m_{y_{\alpha}}}F_{\alpha,3}^{x}+\overline{m_{z_{\alpha}}}F_{\alpha,4}^{x}-F_{\alpha,5}^{x}\bar{\Gamma_{\alpha}}\bigg{)}[\\![\beta_{\alpha}]\\!]$
$\displaystyle=\bar{\rho_{\alpha}}[\\![m_{x_{\alpha}}]\\!]+\overline{m_{x_{\alpha}}}[\\![\rho_{\alpha}]\\!].$

We want this equation to hold for all possible values of the jumps. This is
possible if the coefficients of $[\\![\rho_{\alpha}]\\!]$,
$[\\![\beta_{\alpha}]\\!]$, $[\\![m_{x_{\alpha}}]\\!]$,
$[\\![m_{y_{\alpha}}]\\!]$, and $[\\![m_{z_{\alpha}}]\\!]$, agree on both
sides, which yields the following fluxes,

$\displaystyle\mathbf{\tilde{F}}^{x}_{\alpha,i+\frac{1}{2},j}=\mathbf{\tilde{F}}^{x}_{\alpha}(\mathbf{U}_{\alpha,i,j},\mathbf{U}_{\alpha,i+1,j})=\begin{pmatrix}{\rho_{\alpha}^{\ln}}\overline{m_{x_{\alpha}}}\vspace{0.2cm}\\\
\frac{1}{\overline{\beta_{\alpha}}}\big{(}\frac{\overline{\beta_{\alpha}}\overline{m_{x_{\alpha}}}}{\bar{\Gamma_{\alpha}}}F_{\alpha,5}^{x}+\overline{\rho_{\alpha}}\big{)}\\\
\frac{\overline{m_{y_{\alpha}}}}{\overline{\Gamma_{\alpha}}}F_{\alpha,5}^{x}\\\
\frac{\overline{m_{z_{\alpha}}}}{\overline{\Gamma_{\alpha}}}F_{\alpha,5}^{x}\vspace{0.2cm}\\\
\frac{-\overline{\Gamma_{\alpha}}{\big{(}k_{\alpha}{\rho_{\alpha}^{\ln}}\overline{m_{x_{\alpha}}}+\frac{\overline{m_{x_{\alpha}}}\overline{\rho_{\alpha}}}{\overline{\beta_{\alpha}}}}\big{)}}{\big{(}\overline{m_{x_{\alpha}}}^{2}+\overline{m_{y_{\alpha}}}^{2}+\overline{m_{z_{\alpha}}}^{2}-{(\overline{\Gamma_{\alpha}})}^{2}\big{)}}\end{pmatrix}_{i+\frac{1}{2},j}$
(31)

It is easy to verify that $\tilde{{\mathbf{F}}}_{\alpha}^{x}$ is consistent
with the flux $\mathbf{f}_{\alpha}^{x}$. Similarly, we can derive the
expression for $\mathbf{\tilde{F}}_{\alpha}^{y}$ to get

$\displaystyle\mathbf{\tilde{F}}^{y}_{\alpha,i,j+\frac{1}{2}}=\mathbf{\tilde{F}}^{y}_{\alpha}(\mathbf{U}_{\alpha,i,j},\mathbf{U}_{\alpha,i,j+1})=\begin{pmatrix}{\rho_{\alpha}^{\ln}}\overline{m_{y_{\alpha}}}\vspace{0.2cm}\\\
\frac{\overline{m_{x_{\alpha}}}}{\overline{\Gamma_{\alpha}}}F_{{\alpha,5}}^{y}\\\
\frac{1}{\overline{\beta_{\alpha}}}\bigg{(}\frac{\bar{\beta_{\alpha}}\overline{m_{y_{\alpha}}}}{\overline{\Gamma_{\alpha}}}F_{{\alpha,5}}^{y}+\overline{\rho_{\alpha}}\bigg{)}\\\
\frac{\overline{m_{z_{\alpha}}}}{\overline{\Gamma_{\alpha}}}F_{{\alpha,5}}^{y}\vspace{0.2cm}\\\
\frac{-\overline{\Gamma_{\alpha}}{\left(k_{\alpha}{\rho_{\alpha}^{\ln}}\overline{m_{y_{\alpha}}}+\frac{\overline{m_{y_{\alpha}}}\overline{\rho_{\alpha}}}{\bar{\beta_{\alpha}}}\right)}}{\left(\overline{m_{x_{\alpha}}}^{2}+\overline{m_{y_{\alpha}}}^{2}+\overline{m_{z_{\alpha}}}^{2}-(\overline{\Gamma_{\alpha}})^{2}\right)}\end{pmatrix}_{i,j+\frac{1}{2}},$
(32)

which is consistent with ${\mathbf{f}}^{y}_{\alpha}$.

######  Theorem 4.2.

The numerical scheme (23), with the numerical fluxes (31) and (32) is second-
order accurate and entropy conservative, i.e., it satisfy (24).

###### Proof.

Following [39] and [30], multiplying (23) with
$\mathbf{V}_{\alpha,i,j}^{\top}$, we get,

$\displaystyle\frac{d}{dt}\mathcal{U}_{\alpha}(\mathbf{U}_{\alpha,i,j})+\frac{1}{\Delta
x}\left(\tilde{\mathcal{F}}_{\alpha,i+\frac{1}{2},j}^{x}-\tilde{\mathcal{F}}_{\alpha,i-\frac{1}{2},j}^{x}\right)+\frac{1}{\Delta
y}\left(\tilde{\mathcal{F}}_{\alpha,i,j+\frac{1}{2}}^{y}-\tilde{\mathcal{F}}_{\alpha,i,j-\frac{1}{2}}^{y}\right)$
$\displaystyle=\mathbf{V}_{\alpha,i,j}^{\top}\cdot\mathbf{s}_{\alpha}(\mathbf{U}_{\alpha,i,j}(t),\mathbf{U}_{m,i,j}(t))$

Now, we note that

$\mathbf{V}_{\alpha,i,j}^{\top}\cdot\mathbf{s}_{\alpha}(\mathbf{U}_{\alpha,i,j}(t),\mathbf{U}_{m,i,j}(t))=0$

to arrive at the equality (24). ∎

The entropy conservative fluxes derived above are only second-order accurate.
Following [41], we can use second-order fluxes to construct $2p^{th}$-order
accurate fluxes for any positive integer $p.$ In particular, the
$4^{th}$-order ($p=2$) $x$-directional entropy conservative flux
$\tilde{\mathbf{{F}}}^{x,4}_{\alpha,i+\frac{1}{2},j}$ is given by,

$\tilde{\mathbf{{F}}}^{x,4}_{\alpha,i+\frac{1}{2},j}=\frac{4}{3}\tilde{\mathbf{{F}}}^{x}_{\alpha}(\mathbf{U}_{\alpha,i,j},\mathbf{U}_{\alpha,i+1,j})-\frac{1}{6}\bigg{(}\tilde{\mathbf{{F}}}^{x}_{\alpha}(\mathbf{U}_{\alpha,i-1,j},\mathbf{U}_{\alpha,i+1,j})+\tilde{\mathbf{{F}}}^{x}_{\alpha}(\mathbf{U}_{\alpha,i,j},\mathbf{U}_{\alpha,i+2,j})\bigg{)}.$

Similarly, we can obtain the fourth-order entropy conservative flux
$\tilde{\mathbf{{F}}}^{y,4}_{\alpha,i+\frac{1}{2},j}$ in $y$-direction.
Combining this we have the following remarks.

######  Remark 4.1.

Replacing second-order fluxes with the higher-order fluxes in (23), we arrive
at high-order, entropy conservative schemes.

######  Remark 4.2.

The result from [41] gives only even order accurate fluxes. So, for
$q^{th}$-order (q is an odd integer) accurate scheme, we use a
$(q+1)^{th}$-order (an even number) accurate flux.

####  4.1.2 High-order entropy stable schemes for fluid equations

The entropy conservative schemes presented above will produce high-frequency
oscillations near the shocks as they are central fluxes which do not dissipate
entropy. Following [39], we introduce modified fluxes which will ensure
entropy dissipation at the shocks. We consider a modified numerical flux of
the form,

$\displaystyle\mathbf{F}_{\alpha,i+\frac{1}{2},j}^{x}=\tilde{\mathbf{{F}}}_{\alpha,i+\frac{1}{2},j}^{x}-\frac{1}{2}\textbf{D}_{\alpha,i+\frac{1}{2},j}^{x}[\\![\mathbf{V}_{\alpha}]\\!]_{i+\frac{1}{2},j},$
(33)
$\displaystyle\mathbf{F}_{\alpha,i,j+\frac{1}{2}}^{y}=\tilde{\mathbf{{F}}}_{\alpha,i,j+\frac{1}{2}}^{y}-\frac{1}{2}\textbf{D}_{\alpha,i,j+\frac{1}{2}}^{y}[\\![\mathbf{V}_{\alpha}]\\!]_{i,j+\frac{1}{2}},$

where $\textbf{D}^{x}_{\alpha,i+\frac{1}{2},j}$ and
$\textbf{D}^{y}_{\alpha,i,j+\frac{1}{2}}$ are symmetric positive definite
matrices. Then we have the following Lemma:

######  Lemma 4.1 (Tadmor[39]).

The numerical scheme (21) with the numerical fluxes (33) is entropy stable,
i.e., the entropy inequality

$\frac{d}{dt}\mathcal{U}_{\alpha}(\mathbf{U}_{ij})+\frac{1}{\Delta
x}\left(\hat{\mathcal{F}}_{\alpha,i+\frac{1}{2},j}^{x}-\hat{\mathcal{F}}_{\alpha,i-\frac{1}{2},j}^{x}\right)+\frac{1}{\Delta
y}\left(\hat{\mathcal{F}}_{\alpha,i,j+\frac{1}{2}}^{y}-\hat{\mathcal{F}}_{\alpha,i,j-\frac{1}{2}}^{y}\right)\leq
0,\ \ \ \alpha\in\\{i,e\\},$

is satisfied with the consistent numerical entropy flux functions,

$\hat{\mathcal{F}}^{x}_{\alpha,i+\frac{1}{2},j}=\tilde{\mathcal{F}}^{x}_{\alpha,i+\frac{1}{2},j}+\frac{1}{2}(\bar{\mathbf{V}}_{\alpha})^{\top}_{i+\frac{1}{2},j}\mathbf{D}_{\alpha,i+\frac{1}{2},j}^{x}[\\![\mathbf{V}_{\alpha}]\\!]_{i+\frac{1}{2},j},$
$\hat{\mathcal{F}}^{y}_{\alpha,i,j+\frac{1}{2}}=\tilde{\mathcal{F}}^{y}_{\alpha,i,j+\frac{1}{2}}+\frac{1}{2}(\bar{\mathbf{V}}_{\alpha})^{\top}_{i,j+\frac{1}{2}}\mathbf{D}_{\alpha,i,j+\frac{1}{2}}^{y}[\\![\mathbf{V}_{\alpha}]\\!]_{i,j+\frac{1}{2}}.$

We will use Rusanov type diffusion operators for the matrix
$\textbf{D}_{\alpha}$, which is given by,

$\mathbf{D}_{\alpha,i+\frac{1}{2},j}^{x}={\mathbf{\tilde{R}}}_{\alpha,i+\frac{1}{2},j}^{x}\Lambda_{\alpha,i+\frac{1}{2},j}^{x}{\mathbf{\tilde{R}}}_{\alpha,i+\frac{1}{2},j}^{x\top},\text{
\ \ and \ \
}\textbf{D}_{\alpha,i,j+\frac{1}{2}}^{y}={\mathbf{\tilde{R}}}_{\alpha,i,j+\frac{1}{2}}^{y}\Lambda_{\alpha,i,j+\frac{1}{2}}^{y}{\mathbf{\tilde{R}}}_{\alpha,i,j+\frac{1}{2}}^{y\top}.$
(34)

where ${\mathbf{\tilde{R}}^{d}}_{\alpha},\,d\in\\{x,y\\}$, are matrices of the
entropy scaled right eigenvectors and ${\Lambda^{d}_{\alpha}}$ are $5\times 5$
diagonal matrices of the form

${\Lambda^{d}_{\alpha}}=\left(\max_{1\leq k\leq
5}|\Lambda_{\alpha_{k}}^{d}|\right)\mathbf{I}_{5\times 5},$

where $\\{\Lambda_{\alpha_{k}}^{d}:1\leq k\leq 5\\},\,d\in\\{x,y\\}$, are the
eigenvalues. Following [42], we have
$\partial_{\mathbf{V}_{\alpha}}\mathbf{U}_{\alpha}=\mathbf{\tilde{R}}_{\alpha}^{d}(\mathbf{\tilde{R}}_{\alpha}^{d})^{\top}$
for $d\in\\{x,y\\}$, hence,

$\mathbf{\tilde{R}}^{d}_{\alpha}\Lambda_{\alpha}^{d}(\mathbf{\tilde{R}}^{d}_{\alpha})^{(-1)}[\\![\mathbf{U}_{\alpha}]\\!]\approx\mathbf{\tilde{R}}^{d}_{\alpha}\Lambda^{d}_{\alpha}(\mathbf{\tilde{R}}^{d}_{\alpha})^{(-1)}\partial_{\mathbf{V}_{\alpha}}\mathbf{U}_{\alpha}[\\![\mathbf{V}_{\alpha}]\\!]=\mathbf{\tilde{R}}^{d}_{\alpha}\Lambda^{d}_{\alpha}(\mathbf{\tilde{R}}^{d}_{\alpha})^{\top}[\\![\mathbf{V}_{\alpha}]\\!],$

i.e., the entropy diffusion operator used here is similar to the Roe diffusion
operator when the entropy scaled right eigenvectors are used. The complete
expressions for $\mathbf{\tilde{R}}^{d}_{\alpha}$ are derived in Appendix
(A.2).

With the choice of diffusion operator (34), the numerical scheme (21) with the
numerical flux (33) is entropy stable. However, the scheme is only first-order
accurate due to the presence of the first-order jump terms
$[\\![\mathbf{V}_{\alpha}]\\!]_{i+\frac{1}{2},j}$ and
$[\\![\mathbf{V}_{\alpha}]\\!]_{i,j+\frac{1}{2}}$. A straightforward way to
increase the order of accuracy is to approximate the jumps in the diffusive
term of (33) using a higher-order reconstruction process. However, proving the
entropy stability of the resulting scheme is not possible. Instead, we follow
the process prescribed in [30] and introduce the scaled entropy variables. We
illustrate the procedure for the $x$-direction only, as the $y$-directional
case is similar. Define the change of variables

$\mathcal{V}_{\alpha,{m,j}}^{x,\pm}\,=\,(\mathbf{\tilde{R}}^{x}_{\alpha,{i\pm\frac{1}{2},j}})^{\top}\mathbf{V}_{\alpha,{m,j}},\qquad\textrm{$m$
are neighbours of cell $(i,j)$ along $x$}$

Using ENO procedure we choose a stencil of cells and construct the polynomials
$P^{x,\pm}_{i,j}(x)$ of degree $k$ and evaluate it at the faces of cell
$(i,j)$

$\mathcal{\tilde{V}}_{\alpha,i,j}^{x,\pm}=P_{i,j}^{x,\pm}(x_{i\pm\frac{1}{2}})$

Converting back to the entropy variables,

$\mathbf{\tilde{V}}_{\alpha,{i,j}}^{x,\pm}\,=\,\left\\{(\mathbf{\tilde{R}}^{x}_{\alpha,{i\pm\frac{1}{2},j}})^{\top}\right\\}^{(-1)}\tilde{\mathcal{V}}_{\alpha,{i,j}}^{x,\pm}$

are the corresponding $k$-th order reconstructed values for
$\mathbf{V}_{\alpha}$. Hence the high-order entropy stable numerical flux is
given by,

$\mathbf{F}_{\alpha,i+\frac{1}{2},j}^{x,k}\,=\,\tilde{\mathbf{F}}_{\alpha,i+\frac{1}{2},j}^{x,2p}\,-\,\frac{1}{2}\,\mathbf{D}_{\alpha,i+\frac{1}{2},j}^{x}[\\![\mathbf{\tilde{V}}^{x}_{\alpha}]\\!]_{i+\frac{1}{2},j}$
(35)

where $[\\![\mathbf{\tilde{V}}^{x}_{\alpha}]\\!]_{i+\frac{1}{2},j}$ stands
for,

$[\\![\mathbf{\tilde{V}}^{x}_{\alpha}]\\!]_{i+\frac{1}{2},j}\,=\,\mathbf{\tilde{V}}^{x,-}_{\alpha,{i+1,j}}\,-\,\mathbf{\tilde{V}}^{x,+}_{\alpha,{i,j}}.$

As only even order entropy conservative fluxes
$\tilde{\mathbf{F}}_{\alpha,i+\frac{1}{2},j}^{x,2p}$ are available and we
choose $p\in\mathbb{N}$ as

  * •

$p=k/2$ if $k$ is even,

  * •

$p=(k+1)/2$ if $k$ is odd,

where $k$ is the order of the scheme. Following [30], a sufficient condition
for the numerical flux (35) to be entropy stable is that the reconstruction
process for ${\mathcal{V}}$ must satisfy the sign preserving property. For the
second-order reconstruction, we use the min-mod reconstruction, which
satisfies this property. Following [31], we use ENO reconstruction as it
satisfies the sign preserving property. For third and fourth-order schemes, we
use fourth-order entropy conservative flux along with the third and fourth-
order ENO reconstruction of the scaled entropy variables. We now have the
following result:

######  Theorem 4.3.

The numerical scheme (23) with entropy stable flux (35) satisfies the entropy
inequality (25).

###### Proof.

The proof follows from [30] and the fact that,

$\mathbf{V}_{\alpha,i,j}^{\top}\cdot\mathbf{s}_{\alpha}(\mathbf{U}_{\alpha,i,j}(t),\mathbf{U}_{m,i,j}(t))=0$

∎

######  Remark 4.3.

The above result holds for any discretization of Maxwell’s equations.

###  4.2 Discretization of Maxwell’s equations

To spatially discretize the Maxwell’s equations (1g)-(1j), we use the finite-
difference scheme based on the Rusanov solver. The $x$\- and $y$-directional
numerical fluxes are given by

$\displaystyle\mathbf{F}^{x}_{m,i+\frac{1}{2},j}=\dfrac{1}{2}(\mathbf{f}^{x}_{m}(\mathbf{U}_{m,i,j})+\mathbf{f}^{x}_{m}(\mathbf{U}_{m,i+1,j}))-\dfrac{c}{2}(\mathbf{U}_{m,i+1,j}-\mathbf{U}_{m,i,j}),$
$\displaystyle\mathbf{F}^{y}_{m,i,j+\frac{1}{2}}=\dfrac{1}{2}(\mathbf{f}^{y}_{m}(\mathbf{U}_{m,i,j})+\mathbf{f}^{y}_{m}(\mathbf{U}_{m,i,j+1}))-\dfrac{c}{2}(\mathbf{U}_{m,i,j+1}-\mathbf{U}_{m,i,j}).$

where $c=1$ is the speed of light. To obtain higher accuracy, we use min-mod
limiter [25] to reconstruct the variables $\mathbf{U}_{m}$ at the cell faces
to get the second-order scheme. Furthermore, we use third and fifth-order WENO
reconstructions [43] to obtain the third- and fourth-order Rusanov’s solver
for the Maxwell’s equations.

##  5 Fully discrete scheme

Let $\mathbf{U}^{n}$ be the discrete solution at time $t^{n}$, where the time
step is $\Delta t=t^{n+1}-t^{n}$. The semi-discrete scheme (21) can be
expressed as

$\frac{d}{dt}\mathbf{U}_{i,j}(t)=\mathcal{L}_{i,j}(\mathbf{U}(t))+\mathbf{s}(\mathbf{U}_{i,j}(t)),\qquad\mathbf{U}_{i,j}(t^{n})=\mathbf{U}^{n}$
(37)

where,

$\mathcal{L}_{i,j}(\mathbf{U}(t))=-\frac{1}{\Delta
x}\left(\mathbf{F}_{i+\frac{1}{2},j}^{x}(t)-\mathbf{F}_{i-\frac{1}{2},j}^{x}(t)\right)-\frac{1}{\Delta
y}\left(\mathbf{F}_{i,j+\frac{1}{2}}^{y}(t)-\mathbf{F}_{i,j-\frac{1}{2}}^{y}(t)\right).$

We propose fully explicit and IMEX schemes for the above system.

###  5.1 Explicit schemes

For the explicit time discretizations, we use explicit strong stability
preserving Runge Kutta (SSP-RK) methods [44]. The second and third-order
accurate SSP-RK schemes are given as follows.

  1. 1.

Set $\mathbf{U}^{0}\ =\ \mathbf{U}^{n}$.

  2. 2.

For $m=1,\dots,k+1$, compute,

$\mathbf{U}_{i,j}^{(m)}\ =\
\sum_{l=0}^{m-1}\alpha_{ml}\mathbf{U}_{i,j}^{(l)}+\beta_{ml}\Delta
t\big{(}\mathcal{L}_{i,j}(\mathbf{U}^{(l)})+\mathbf{s}(\mathbf{U}_{i,j}^{(l)})\big{)},$

where, $\alpha_{ml}$ and $\beta_{ml}$ follows from the Table (1).

  3. 3.

Finally, $\mathbf{U}_{i,j}^{n+1}\ =\ \mathbf{U}_{i,j}^{(k+1)}$.

Order | $\alpha_{il}$ | $\beta_{il}$
---|---|---
2 | 1 |  |  | 1 |  |
1/2 | 1/2 |  | 0 | 1/2 |
3 | 1 |  |  | 1 |  |
3/4 | 1/4 |  | 0 | 1/4 |
1/3 | 0 | 2/3 | 0 | 0 | 2/3
Table 1: Coefficients for SSP Runge-Kutta time stepping

The fourth order RK-SSP is given by,

$\displaystyle\mathbf{U}^{(1)}_{i,j}$
$\displaystyle=\textbf{U}^{n}_{i,j}+0.39175222700392\;(\Delta
t)\left(\mathcal{L}_{i,j}(\mathbf{U}^{n})+\mathbf{s}(\mathbf{U}^{n}_{i,j})\right),$
$\displaystyle\mathbf{U}^{(2)}_{i,j}$
$\displaystyle=0.44437049406734\;\textbf{U}^{n}_{i,j}+0.55562950593266\;\mathbf{U}^{(1)}_{i,j}$
$\displaystyle+0.36841059262959\;(\Delta
t)\left(\mathcal{L}_{i,j}(\mathbf{U}^{(1)})+\mathbf{s}(\mathbf{U}^{(1)}_{i,j})\right),$
$\displaystyle\mathbf{U}^{(3)}_{i,j}$
$\displaystyle=0.62010185138540\;\textbf{U}^{n}_{i,j}+0.37989814861460\;\mathbf{U}^{(2)}_{i,j}$
$\displaystyle+0.25189177424738\;(\Delta
t)\left(\mathcal{L}_{i,j}(\mathbf{U}^{(2)})+\mathbf{s}(\mathbf{U}^{(2)}_{i,j})\right),$
$\displaystyle\mathbf{U}^{(4)}_{i,j}$
$\displaystyle=0.17807995410773\;\textbf{U}^{n}_{i,j}+0.82192004589227\;\mathbf{U}^{(3)}_{i,j}$
$\displaystyle+0.54497475021237\;(\Delta
t)\left(\mathcal{L}_{i,j}(\mathbf{U}^{(3)})+\mathbf{s}(\mathbf{U}^{(3)}_{i,j})\right),$
$\displaystyle\textbf{U}^{n+1}_{i,j}$
$\displaystyle=0.00683325884039\;\textbf{U}^{n}+0.51723167208978\;\mathbf{U}^{(2)}_{i,j}+0.12759831133288\;\mathbf{U}^{(3)}_{i,j}$
$\displaystyle+0.34833675773694\;\mathbf{U}^{(4)}_{i,j}+0.08460416338212\;(\Delta
t)\left(\mathcal{L}_{i,j}(\mathbf{U}^{(3)})+\mathbf{s}(\mathbf{U}^{(3)}_{i,j})\right)$
$\displaystyle+0.22600748319395\;(\Delta
t)\left(\mathcal{L}_{i,j}(\mathbf{U}^{(4)})+\mathbf{s}(\mathbf{U}^{(4)}_{i,j})\right).$

We use a second-order spatial accurate scheme with a second-order explicit
time scheme, which is denoted by, O2-ES-Exp. Similarly, the third-order
spatial accurate scheme is time-discretized with the third-order RK method
above and denoted by, O3-ES-Exp. The fourth-order scheme consists of fourth-
order spatial discretization with the fourth-order time update above and
denoted by O4-ES-Exp.

###  5.2 ARK-IMEX schemes

In several cases, the values of charge to mass ratios $r_{\alpha}$ can be very
large. In those cases, the source terms will be stiff. This will result in a
time-step restriction when using the explicit scheme discussed above. To
overcome this difficulty, we propose Additive Runge Kutta Implicit Explicit
(ARK IMEX) time updates from[45]. Here, the flux terms will be treated
explicitly, and the source terms will be treated implicitly. We write the
semi-discrete scheme (21) in a simplified form as follows

$\dfrac{\partial\mathbf{U}}{\partial
t}=F_{NS}(\mathbf{U}(t))+F_{S}(\mathbf{U}(t)),$

where $F_{NS}(\mathbf{U}(t))=\mathcal{L}_{i,j}(\mathbf{U}(t))$ denotes the
non-stiff flux terms, and
$F_{S}(\mathbf{U}(t))=\mathbf{s}(\mathbf{U}_{i,j}(t))$ denotes the stiff
source terms. Note that the source terms depend only on the solution value in
the parent cell and do not depend on the solutions in neighboring cells.

####  5.2.1 Second order time discretization.

For the second-order discretization, we use L-Stable second-order accurate ARK
IMEX scheme from [45]. The method has two implicit stages and is given by

$\displaystyle\mathbf{U}^{(1)}$ $\displaystyle=\mathbf{U}^{n}+\Delta t\
\big{(}\beta F_{S}(\mathbf{U}^{(1)}(t^{n}))\big{)},$
$\displaystyle\mathbf{U}^{(2)}$ $\displaystyle=\mathbf{U}^{n}+\Delta t\
\big{(}F_{NS}(\mathbf{U}^{(1)}(t^{n}))+(1-2\beta)F_{S}(\mathbf{U}^{(1)}(t^{n}))+\beta
F_{S}(\mathbf{U}^{(2)}(t^{n}))\big{)},$ $\displaystyle\mathbf{U}^{(n+1)}$
$\displaystyle=\mathbf{U}^{n}+\Delta t\
\bigg{(}\dfrac{1}{2}F_{NS}(\mathbf{U}^{(1)}(t^{n}))+\dfrac{1}{2}F_{NS}(\mathbf{U}^{(2)}(t^{n}))+\dfrac{1}{2}F_{S}(\mathbf{U}^{(1)}(t^{n}))+\dfrac{1}{2}F_{S}(\mathbf{U}^{(2)}(t^{n}))\bigg{)},$

where $\beta=1-\dfrac{1}{\sqrt{2}}$. Note that, to evaluate $\mathbf{U}^{1}$
and $\mathbf{U}^{2}$, we need to solve a nonlinear system of algebraic
equations. However, the resulting equations are only local, i.e., we get an
independent set of equations for each cell. We solve these equations using the
Newton’s method, where the convergence and the next direction search is based
on the backtracing line search framework. A detailed procedure is given in
[46]. In the higher-order method, we will use the same process to update the
internal IMEX steps.

####  5.2.2 Third and fourth-order time discretization.

For the higher-order discretization, we use L-Stable ARK IMEX methods from
[47]. We ignore the subscripts $\\{i,j\\}$ for simplicity. The ARK scheme from
[47] has the following steps.

  1. 1.

Set $\mathbf{U}^{0}\ =\ \mathbf{U}^{n}$.

  2. 2.

For $m=1,\dots,k$, compute,

$\mathbf{U}^{(m)}=\mathbf{U}^{n}+\Delta
t\sum_{l=0}^{m-1}a_{ml}^{[NS]}F_{NS}(\mathbf{U}^{l})+\Delta
t\sum_{l=0}^{m}a_{ml}^{[S]}F_{S}(\mathbf{U}^{l}),$

where, $\mathbf{U}^{(m)}$ approximates $\mathbf{U}(t^{n}+c_{m}\Delta t).$

  3. 3.

Finally,

$\mathbf{U}^{n+1}=\mathbf{U}^{n}+\Delta
t\sum_{l=0}^{k}b_{l}^{[NS]}F_{NS}(\mathbf{U}^{l})+\Delta
t\sum_{l=0}^{k}b_{l}^{[S]}F_{S}(\mathbf{U}^{l}),$

where $k$ is the order of the desired time discretization ($k$ ranges from 3
to 5), and coefficients $a_{ml}^{[NS]}$, $a_{ml}^{[S]}$, $b_{l}^{[NS]}$,
$b_{l}^{[S]}$, $c_{l}$ are defined accordingly by the order of accuracy and
stability considerations.

We present coefficients for the third-order scheme in Appendix (B). The
coefficients for the fourth-order time update can be found in [47]. We denote
the second-order IMEX scheme with O2-ES-IMEX, third-order IMEX scheme with
O3-ES-IMEX and fourth-order IMEX scheme with O4-ES-IMEX.

##  6 Numerical results

In this section, we present numerical results for various test cases. We set
specific heat constants $\gamma_{i}=\gamma_{e}=4/3$, unless stated otherwise.
In some test cases, we need to consider the resistive effects. To take them
into consideration, following [12], we modify the momentum ((1e),(1b)) and
energy equations ((1c),(1f)) for the fluid parts by adding resistive terms as
follows,

$\displaystyle\frac{\partial(\rho_{\alpha}h_{\alpha}\Gamma_{\alpha}^{2}\mathbf{u}_{\alpha})}{\partial
t}+\nabla\cdot(\rho_{\alpha}h_{\alpha}\Gamma_{\alpha}^{2}\textbf{u}_{\alpha}\mathbf{u}_{\alpha}^{\top}+p_{\alpha}\mathbf{I})$
$\displaystyle=r_{\alpha}\Gamma_{\alpha}\rho_{\alpha}(\mathbf{E}+\mathbf{u}_{\alpha}\times\mathbf{B})+\boldsymbol{R}_{\alpha},$
(40a) $\displaystyle\frac{\partial(\mathcal{E}_{\alpha})}{\partial
t}+\nabla\cdot((\mathcal{E}_{\alpha}+p_{\alpha})\mathbf{u}_{\alpha})$
$\displaystyle=r_{\alpha}\Gamma_{\alpha}\rho_{\alpha}(\mathbf{u}_{\alpha}\cdot\mathbf{E})+R_{\alpha}^{0},$
(40b)

for $\alpha\in\\{i,e\\}$. Following[12], an anti-symmetry relationship,
$(\boldsymbol{R}_{e},{R}_{e}^{0})=(-\boldsymbol{R}_{i},-R_{i}^{0})$, has been
assumed for the conservation of the total momentum and energy density. The
expressions for the terms $\boldsymbol{R}_{i}$ and $R_{i}^{0}$ are given by

$\displaystyle\boldsymbol{R}_{i}=-\eta\dfrac{\omega_{p}^{2}}{r_{i}-r_{e}}(\mathbf{j}-\rho_{0}\boldsymbol{\Phi}),\qquad{R}_{i}^{0}=-\eta\dfrac{\omega_{p}^{2}}{r_{i}-r_{e}}(\rho_{c}-\rho_{0}{\Lambda}),$

where,

$\displaystyle\omega_{p}^{2}$
$\displaystyle=r_{i}^{2}\rho_{i}+r_{e}^{2}\rho_{e},\qquad\boldsymbol{\Phi}=\dfrac{r_{i}\rho_{i}\Gamma_{i}\mathbf{u}_{i}+r_{e}\rho_{e}\Gamma_{e}\mathbf{u}_{e}}{\omega_{p}^{2}},$
(41) $\displaystyle\Lambda$
$\displaystyle=\dfrac{r_{i}^{2}\rho_{i}\Gamma_{i}+r_{e}^{2}\rho_{e}\Gamma_{e}}{\omega_{p}^{2}},\qquad\rho_{0}=\Lambda\rho_{c}-\mathbf{j}\cdot\boldsymbol{\Phi}$

and $\eta$ is the resistivity constant. Here, $\omega_{p}$ is the total plasma
frequency. The total plasma skin depth $d_{p}$ is defined as
$d_{p}=\frac{1}{\omega_{p}}$. For the time update using ARK IMEX scheme, we
couple the resistive terms with the other source terms and treat them
implicitly, i.e., we modify $F_{S}$ to include the resistive terms.

To compare the numerical results with published work consistently, for all the
test cases taken from [29], we multiply the source term for Maxwell’s
equations by 4$\pi$, i.e., we consider,

$\mathbf{s}_{m}=4\pi\begin{pmatrix}0\\\ 0\\\ 0\\\ -j_{x}\\\ -j_{y}\\\
-j_{z}\end{pmatrix}.$ (42)

Consequently, in the test cases from [29], the total plasma skin depth is
given by $d_{p}=\frac{\sqrt{2}}{\omega_{p}}$ with
$\omega_{p}^{2}=4\pi(r_{i}^{2}\rho_{i}+r_{e}^{2}\rho_{e}).$

###  6.1 One-dimensional test cases

For the one-dimensional test cases, the time step is chosen using

${\Delta t}=\text{CFL}\cdot\min\left\\{\dfrac{\Delta
x}{\Lambda_{max}^{x}(\mathbf{U}_{i})}:1\leq i\leq N_{x}\right\\},\textrm{
where
}\Lambda_{max}^{x}(\mathbf{U}_{i})=\max\\{|\Lambda^{x}_{k}(\mathbf{U}_{i})|:1\leq
k\leq 16\\}.$

We take CFL to be $0.8$, for the explicit (O2-ES-Exp, O3-ES-Exp and O4-ES-
Exp), and ARK IMEX schemes (O2-ES-IMEX, O3-ES-IMEX and O4-ES-IMEX), unless
stated otherwise.

####  6.1.1 Accuracy test

Number of cells | O2-ES-Exp | O3-ES-Exp | O4-ES-Exp
---|---|---|---
– |  $L^{1}$ error | Order |  $L^{1}$ error | Order |  $L^{1}$ error | Order
50 | 4.65354e-02 | – | 1.45685e-03 | – | 3.72752e-05 | –
100 | 1.46178e-02 | 1.6706 | 2.10073e-04 | 2.7939 | 2.83586e-06 | 3.7164
200 | 4.06539e-03 | 1.8463 | 2.68670e-05 | 2.9670 | 2.00245e-07 | 3.8239
400 | 1.10156e-03 | 1.8838 | 3.38951e-06 | 2.9867 | 1.38798e-08 | 3.8507
800 | 2.96426e-04 | 1.8938 | 4.12557e-07 | 3.0384 | 9.36376e-10 | 3.8898
1600 | 7.83484e-05 | 1.9197 | 4.74611e-08 | 3.1198 | 6.20992e-11 | 3.9144
3200 | 2.04441e-05 | 1.9382 | 5.12937e-09 | 3.2099 | 4.42086e-12 | 3.8122
Table 2: Accuracy test: $L^{1}$ errors and order of convergence for $\rho_{i}$ using the explicit entropy stable schemes O2-ES-Exp, O3-ES-Exp and O4-ES-Exp. Number of cells | O2-ES-IMEX | O3-ES-IMEX | O4-ES-IMEX
---|---|---|---
– |  $L^{1}$ error | Order |  $L^{1}$ error | Order |  $L^{1}$ error | Order
50 | 4.61786e-02 | – | 1.47485e-03 | – | 1.49887e-04 | –
100 | 1.44393e-02 | 1.6772 | 2.11249e-04 | 2.8036 | 1.13242e-05 | 3.7264
200 | 4.00815e-03 | 1.8490 | 2.70271e-05 | 2.9665 | 7.99400e-07 | 3.8243
400 | 1.09588e-03 | 1.8708 | 3.40274e-06 | 2.9896 | 5.54233e-08 | 3.8504
800 | 2.94580e-04 | 1.8954 | 4.13678e-07 | 3.0401 | 3.73687e-09 | 3.8906
1600 | 7.77739e-05 | 1.9213 | 4.75938e-08 | 3.1197 | 2.47401e-10 | 3.9169
3200 | 2.02876e-05 | 1.9387 | 5.14414e-09 | 3.2098 | 1.63538e-11 | 3.9192
Table 3: Accuracy test: $L^{1}$ errors and order of convergence for $\rho_{i}$
using the ARK IMEX entropy stable schemes O2-ES-IMEX, O3-ES-IMEX and O4-ES-
IMEX.

To test the accuracy and order of convergence of the proposed scheme, we
consider a test case with smooth solution. We follow a forced solution
approach of [15] to modify the smooth test case from [1]. We add forcing term
$\mathcal{S}(x,t)$ to write the one-dimensional modified system (4) as

$\frac{\partial\mathbf{U}}{\partial t}+\frac{\partial\mathbf{f}^{x}}{\partial
x}=\mathbf{s}+\mathcal{S}(x,t)$

with

$\mathcal{S}(x,t)=\left(\mathbf{0}_{13},-\dfrac{1}{\sqrt{3}}(2+\sin(2\pi(x-0.5t))),0,-3\pi\cos(2\pi(x-0.5t))\right)^{\top}.$

We set initial ion and electron densities as $\rho_{i}=\rho_{e}=2+\sin(2\pi
x)$, with initial velocities $u_{x_{i}}=u_{x_{e}}=0.5$ and initial pressures
$p_{i}=p_{e}=1$. The $y-$magnetic component is $B_{y}=2\sin(2\pi(x))$ and the
$z-$electric field component is $E_{z}=-\sin(2\pi(x)$. All other primitive
variables are set to zero. We use computational domain $I=[0,1]$ with the
periodic boundary conditions. We also set charge to mass ratios $r_{i}=1$ and
$r_{e}=-2$, consequently we have non-zero electric field source term
$\mathbf{j}$. For the present test, we use ion-electron adiabatic index
$\gamma=5/3$. Using the above conditions, it is easy to verify that the exact
solution is $\rho_{i}=\rho_{e}=2+\sin(2\pi(x-0.5t))$.

In Table (2) and Table (3), we have presented the results for explicit (O2-ES-
Exp, O3-ES-Exp and O4-ES-Exp) and ARK-IMEX schemes (O2-ES-IMEX, O3-ES-IMEX,
and O4-ES-IMEX) at the final time $t=2.0$. We have presented $L^{1}$-errors
for $\rho_{i}$ at different resolutions. For both explicit and IMEX schemes,
we observe a consistent order of accuracy at various resolutions. Furthermore,
we also note that the errors for both explicit and IMEX schemes of same order
are comparable at a given resolution.

####  6.1.2 Relativistic Brio–Wu test problem with finite plasma skin depth

In this test case, we consider one-dimensional Brio-Wu shock tube problem with
finite plasma skin depth from [29, 12] . This problem is a two-fluid extension
of the relativistic analogue of the Brio-Wu shock tube problem given in [7].
In this case, we consider the Maxwell’s source (42), and no resistive effects
are considered. We take the computational domain to be $[-0.5,0.5]$ with
Neumann boundary conditions at both boundaries. The initial discontinuity is
assumed to be placed at $x=0.0$. Following [29], the initial Riemann data is
given by,

$\displaystyle\mathbf{W_{L}}=\begin{pmatrix}\rho_{i}\\\ p_{i}\\\ \rho_{e}\\\
p_{e}\\\ B_{x}\\\ B_{y}\end{pmatrix}_{L}=\begin{pmatrix}0.5\\\ 0.5\\\ 0.5\\\
0.5\\\ \sqrt{\pi}\\\
\sqrt{4\pi}\end{pmatrix},\qquad\mathbf{W_{R}}=\begin{pmatrix}\rho_{i}\\\
p_{i}\\\ \rho_{e}\\\ p_{e}\\\ B_{x}\\\
B_{y}\end{pmatrix}_{R}=\begin{pmatrix}0.0625\\\ 0.05\\\ 0.0625\\\ 0.05\\\
\sqrt{\pi}\\\ -\sqrt{4\pi}\end{pmatrix}.$

All other variables are set to zero. Following [29], the specific heat ratios
are set to $\gamma_{i}=\gamma_{e}=2.0$, and we compute the solution till final
time $t=0.4$.

First, we consider charge to mass ratios as $r_{i}=-r_{e}=10^{3}/\sqrt{4\pi}$,
thus the plasma skin depth is $10^{-3}/\sqrt{\rho_{i}}$. Consequently, we
expect solutions to be close to the RMHD results whenever the resolution of
the cells is larger than the plasma skin depth.

(a) Plots of total density $\rho_{i}+\rho_{e}$ for second-order schemes O2-ES-
Exp and O2-ES-IMEX using $400$, $1600$ and $6400$ cells.

(b) Plots of total density $\rho_{i}+\rho_{e}$, for second (O2-ES-Exp and
O2-ES-IMEX), third (O3-ES-Exp and O3-ES-IMEX), and fourth (O4-ES-Exp and
O4-ES-IMEX) order schemes using $400$ cells.

(c) Plots of total ion entropy for second (O2-ES-Exp and O2-ES-IMEX), third
(O3-ES-Exp and O3-ES-IMEX), and fourth (O4-ES-Exp and O4-ES-IMEX) order
schemes using $400$ cells.

(d) Plots of total electron entropy for second (O2-ES-Exp and O2-ES-IMEX),
third (O3-ES-Exp and O3-ES-IMEX), and fourth (O4-ES-Exp and O4-ES-IMEX) order
schemes using $400$ cells.

Figure 1: Numerical results for Relativistic Brio–Wu test problem with finite
plasma skin depth using different schemes and grid resolutions.

(a) Plots of total density $\rho_{i}+\rho_{e}$, for second (O2-ES-Exp and
O2-ES-IMEX), third (O3-ES-Exp and O3-ES-IMEX), and fourth (O4-ES-Exp and
O4-ES-IMEX) order schemes with $r_{i}=-r_{e}=10^{4}/\sqrt{4\pi}$ using $400$
cells.

(b) Plots of total density $\rho_{i}+\rho_{e}$, for second (O2-ES-Exp and
O2-ES-IMEX), third (O3-ES-Exp and O3-ES-IMEX), and fourth (O4-ES-Exp and
O4-ES-IMEX) order schemes with $r_{i}=-r_{e}=10^{5}/\sqrt{4\pi}$ using $400$
cells.

(c) Plots of total density $\rho_{i}+\rho_{e}$, for O4-ES-IMEX scheme for
various stiff cases at $400$ cells.

Figure 2: Numerical results for Relativistic Brio–Wu test problem with finite
plasma skin depth using various stiff cases at 400 cells.

To compare with the results of [29], in Figure (1(a)), we plot the total
density $\rho_{i}+\rho_{e}$ for O2-ES-Exp and O2-ES-IMEX schemes using $400$,
$1600$ and $6400$ cells. To compare the solution with RMHD solution, we have
computed “RMHD” solution using $r_{i}=-r_{e}=10^{5}/\sqrt{4\pi}$ at $6400$
cells and O2-ES-IMEX scheme. We observe that on a resolution of $400$ cells,
the cell size $\Delta x=0.0025$ is slightly larger than the plasma skin depth
in the entire domain; thus, the profile roughly matches with the RMHD results
[7], while at the higher resolutions, $1600$ and $6400$, the cell size is much
smaller than the plasma skin depth, consequently, appearance of the two-fluid
effect in the form of dispersive waves is observed. We also note that both
explicit (O2-ES-Exp) and IMEX (O2-ES-IMEX) schemes have similar performance in
capturing all the waves and additional effects at the various resolutions.

To compare second-order schemes, with the third and fourth-order schemes, in
Figure (1(b)), we have plotted total density for explicit (O2-ES-Exp, O3-ES-
Exp and O4-ES-Exp) and ARK-IMEX schemes (O2-ES-IMEX, O3-ES-IMEX, and O4-ES-
IMEX) at $400$ cells. We observe that the third and fourth-order schemes are
able to capture the additional dispersive effects even on the coarser mesh. We
also note that fourth-order schemes are more effective in capturing these
effects compared to the third-order scheme. Furthermore, both third-order
schemes (O3-ES-Exp and O3-ES-IMEX) have similar performance. Similarly, both
fourth-order schemes (O4-ES-Exp and O4-ES-IMEX) have similar performance.

We can also observe this from the time evolution of total entropy plots for
ion flow in Figure (1(c)) and electron flow in Figure (1(d)). We note that the
same order schemes have almost the same entropy decay behaviour. Furthermore,
third and fourth-order schemes decay less entropy compared to the second-order
schemes. Similarly, fourth-order schemes are less entropy decaying than the
third-order schemes.

| $r_{i}=-r_{e}=10^{4}/\sqrt{4\pi}$ | $r_{i}=-r_{e}=10^{5}/\sqrt{4\pi}$
---|---|---
Number of cells | $100$ | $200$ | $400$ | $800$ | $100$ | $200$ | $400$ | $800$
O2-ES-Exp | 3.31s | 6.0s | 11.04s | 20.81s | 76.99s | 157.17s | 288.66s | 506.90s
O2-ES-IMEX | 0.45s | 1.06s | 4.17s | 11.27s | 0.61s | 1.41s | 4.40s | 15.51s
O3-ES-Exp | 2.19s | 4.34s | 7.90s | 13.62s | 29.86s | 54.09s | 98.92s | 190.71s
O3-ES-IMEX | 0.86s | 2.64s | 6.57s | 24.79s | 1.13s | 2.67s | 6.59s | 25.13s
O4-ES-Exp | 5.84s | 10.02s | 18.03s | 34.61s | 56.53s | 102.79s | 198.36s | 359.40s
O4-ES-IMEX | 1.28s | 3.85s | 13.96s | 54.91s | 1.32s | 4.23s | 14.31s | 49.10s
Table 4: Relativistic Brio–Wu test problem with finite plasma skin depth:
Computational times using $r_{i}=-r_{e}=10^{4}/\sqrt{4\pi}$ and
$r_{i}=-r_{e}=10^{5}/\sqrt{4\pi}$ for the second (O2-ES-Exp and O2-ES-IMEX),
third (O3-ES-Exp and O3-ES-IMEX), and fourth (O4-ES-Exp and O4-ES-IMEX) order
schemes.

We also consider the stiff source terms with $r_{i}=-r_{e}=10^{4}/\sqrt{4\pi}$
and $r_{i}=-r_{e}=10^{5}/\sqrt{4\pi}.$ We remark that in this case plasma skin
depth is very small, hence we do not expect any scheme to resolve the
dispersive effects and we will need substantial more number of cells to
resolve them. In Figure (2), we have plotted solutions for different schemes
using 400 cells. We note that for both the value of $r_{i}$, all the schemes
produce solutions close to RMHD solution, with third and fourth-order schemes
performing better than the second-order schemes. Also, there is no visible
difference in Explicit and IMEX scheme of same order. To further show the
effects of higher values of $r_{i}$, in Figure (2(c)), we have plotted
solutions using O4-ES-IMEX scheme on $400$ cells, for different values of
$r_{i}$. For $r_{i}=10^{3}/\sqrt{4\pi}$, we can see the dispersive effects are
resolved. However, for $r_{i}=10^{4}/\sqrt{4\pi}$ and
$r_{i}=10^{5}/\sqrt{4\pi}$, scheme was not able to resolve the small scale
oscillation; nevertheless, the scheme produce a stable solution thanks to the
built in entropy stability property of the scheme.

To demonstrate the effectiveness of ARK-IMEX schemes over the explicit schemes
for the stiff source terms, in Table (4), we present the computational time
for the cases $r_{i}=-r_{e}=10^{4}/\sqrt{4\pi}$ and
$r_{i}=-r_{e}=10^{5}/\sqrt{4\pi}$ on $100$, $200$, $400$, and $800$ cells. We
use the MPI-parallelized code on 10 cores. At the lower resolutions, the
source terms will govern the time step as time step given by the CFL
restriction is large. Hence, we observe consistently that IMEX schemes
outperform the explicit schemes on the coarser mesh. In the case of
$r_{i}=-r_{e}=10^{4}/\sqrt{4\pi}$, we see that IMEX outperform the explicit
schemes up to $400$ cells. On $800$ cells, we observe that the IMEX schemes
are more expensive than the explicit schemes as at this resolution, stable
time step is governed by the CFL condition, not the source terms. In the case
of $r_{i}=-r_{e}=10^{5}/\sqrt{4\pi}$, we observe that IMEX schemes
consistently outperform explicit scheme even on $800$ cells as the stiffness
in the source is very strong and hence it decides the stable time step for
explicit schemes.

Given the similar accuracy of the explicit and IMEX schemes of the same order,
and the ability of IMEX schemes to overcome time-step restriction due to stiff
source terms, especially at the lower resolutions, in the remaining test
cases, we will present numerical results only for the ARK-IMEX schemes. We
have also computed the results using explicit schemes (whenever source terms
are not stiff) also, and the solutions are similar to those of the IMEX
schemes of the same order.

####  6.1.3 Self-similar current sheet with finite resistivity

In this test case, we consider a self-similar current sheet which includes the
finite resistivity effects. The test case was first considered in [10], and
the two-fluid relativistic extension is presented in [12], which is used here.
We consider the modified Equations (40a) and (40b) to capture the resistive
effects.

For resistive RMHD solution [10], only the $B_{y}$-component has non-zero
variation, whose time evolution includes the resistive effects and is governed
by the diffusion equation,

$\dfrac{\partial B_{y}}{\partial t}-D\dfrac{\partial^{2}B_{y}}{\partial
x^{2}}=0.$ (43)

The relation between the diffusion coefficient and resistivity constant is
given by the expression $D=\eta c^{2}$. An exact self similar solution of Eqn.
(43) in the resistive RMHD regime [10] is given by

$B_{y}(x,t)=B_{0}\text{ erf}\left(\dfrac{x}{2\sqrt{Dt}}\right),$

where “erf” is the error function.

For two-fluid relativistic plasma equations, we consider the initial $B_{y}$
at the time $t=1.0$, using $B_{0}=1.0$ and $\eta c^{2}=0.01$. The domain is
$[-1.5,1.5]$ with Neumann boundary conditions. The charge to mass ratios are
$r_{i}=-r_{e}=10^{3}$. The ion and electron density are 0.5, and the ion and
electron pressure are 25.0. The $z$-component of the ion and electron velocity
is given by

$u_{z_{i}}=-u_{z_{e}}=\dfrac{B_{0}}{r_{i}\rho_{i}\sqrt{\pi
D}}\exp\left(-\dfrac{x^{2}}{4D}\right).$

All other variables are set to zero. We evolve the solution from the initial
time $t=1$ to final time $t=9.0$ using ARK-IMEX schemes O2-ES-IMEX, O3-ES-IMEX
and O4-ES-IMEX. The solutions using 400 cells are plotted in Figure (3). We
observe that, in the presence of friction terms as in Eqn. (40a), (40b), the
obtained $B_{y}$ profile for all the schemes matches with the resistive RMHD
analytical expression of the $B_{y}$-component. Furthermore, all the schemes
are highly accurate and there is no visible difference in the solution
profiles .

Figure 3: Self-similar current sheet with finite resistivity: Comparison of
the $B_{y}$ profile for ARK-IMEX schemes (O2-ES-IMEX, O3-ES-IMEX, and O4-ES-
IMEX) using $400$ cells.

###  6.2 Two-dimensional test cases

For two-dimensional test cases, the time step is chosen using,

${\Delta t}=\text{CFL}\cdot\min\left\\{\left(\dfrac{\Delta
x}{\Lambda_{max}^{x}(\mathbf{U}_{i,j})}+\dfrac{\Delta
y}{\Lambda_{max}^{y}(\mathbf{U}_{i,j})}\right):1\leq i\leq N_{x},\ 1\leq y\leq
N_{y}\right\\}$

where
$\Lambda_{max}^{x}(\mathbf{U}_{i,j})=\max\\{|\Lambda^{x}_{k}(\mathbf{U}_{i,j})|:1\leq
k\leq 16\\}$ and
$\Lambda_{max}^{y}(\mathbf{U}_{i,j})=\max\\{|\Lambda^{y}_{k}(\mathbf{U}_{i,j})|:1\leq
k\leq 16\\}$. We take CFL to be $0.45$ and present the numerical results for
the ARK-IMEX schemes (O3-ES-IMEX, and O4-ES-IMEX). We have not presented the
results for second-order scheme O2-ES-IMEX as the scheme is less accurate when
compared to the third and fourth order schemes.

####  6.2.1 Relativistic Orzag-Tang test case

Orzag-Tang test case was first proposed for the non-relativistic MHD systems
in [48]. Here, we consider the two-fluid relativistic extension of the test
case, presented in [29]. To have consistent comparison with the results in
[29], we consider the Maxwell’s equations with scaled source (42). We consider
the computational domain $[0.0,1.0]\times[0.0,1.0]$ with periodic boundary
conditions. The initial conditions are given by,

$\displaystyle\begin{pmatrix}[c]\rho_{i}\\\ u_{x_{i}}\\\ u_{y_{i}}\\\
p_{i}\end{pmatrix}=\begin{pmatrix}[c]\frac{25}{72\pi}\\\ -\frac{\sin(2\pi
y)}{2}\\\ \frac{\sin(2\pi x)}{2}\\\
\frac{5}{24\pi}\end{pmatrix},\qquad\begin{pmatrix}[c]\rho_{e}\\\ u_{x_{e}}\\\
u_{y_{e}}\\\ p_{e}\end{pmatrix}=\begin{pmatrix}[c]\frac{25}{72\pi}\\\
-\frac{\sin(2\pi y)}{2}\\\ \frac{\sin(2\pi x)}{2}\\\
\frac{5}{24\pi}\end{pmatrix},\qquad\begin{pmatrix}[c]B_{x}\\\
B_{y}\end{pmatrix}=\begin{pmatrix}[c]-\sin(2\pi y)\\\ \sin(4\pi
x)\end{pmatrix}$

The initial electric field is given by $-\mathbf{u}_{i}\times\mathbf{B}$ and
all other variables are set to zero. We use the charge to mass ratios of
$r_{i}=-r_{e}=10^{3}/\sqrt{4\pi}$, thus the plasma skin depth approximates to
$3.0\times 10^{-3}$. We choose adiabatic indices $\gamma_{i}=\gamma_{e}=5/3$
and simulate till time $t=1.0$ on a mesh with $200\times 200$ cells.

(a) Plot of total density $\rho_{i}+\rho_{e}$.

(b) Plot of total pressure $p_{i}+p_{e}$.

(c) Plot of ion Lorentz factor $\Gamma_{i}$.

(d) Plot of magnitude of the Magnetic field, $\dfrac{|\mathbf{B}|^{2}}{2}$.

Figure 4: Relativistic Orzag-Tang test case: Plots of total density, total
pressure, Ion Lorentz factor, and magnitude of magnetic field
$\dfrac{|\mathbf{B}|^{2}}{2}$ using O3-ES-IMEX scheme and $200\times 200$
cells at time $t=0.1$. We have plotted 30 contours for each variable.

(a) Plot of total density $\rho_{i}+\rho_{e}$.

(b) Plot of total pressure $p_{i}+p_{e}$.

(c) Plot of ion Lorentz factor $\Gamma_{i}$.

(d) Plot of magnitude of the Magnetic field, $\dfrac{|\mathbf{B}|^{2}}{2}$.

Figure 5: Relativistic Orzag-Tang test case: Plots of total density, total
pressure, Ion Lorentz factor, and magnitude of magnetic field
$\dfrac{|\mathbf{B}|^{2}}{2}$ using O4-ES-IMEX scheme and $200\times 200$
cells at time $t=0.1$. We have plotted 30 contours for each variable. Figure
6: Relativistic Orzag-Tang test case: Total fluid entropy
$\mathcal{U}_{i}+\mathcal{U}_{e}$ evolution for the schemes ES-O3-IMEX and
ES-O4-IMEX.

In Figure (4) and (5), we have plotted the total density, total pressure, ion
Lorentz factor and magnitude of the magnetic field for O3-ES-IMEX and O4-ES-
IMEX schemes, respectively. The results for both the schemes are comparable to
those in [29]. Furthermore, both the schemes have similar performance. As the
cell size is higher than the plasma skin depth, the solution is comparable to
those of the RMHD case.

In Figure (6), we plot the time evolution of the total fluid entropy
$\mathcal{U}_{i}+\mathcal{U}_{e}$ for O3-ES-IMEX and O4-ES-IMEX schemes.
Initially, we do not see significant decay as discontinuities are not present
in the solution. However, at time approximately $t=0.4$, the solutions start
to develop discontinuities and we see sharp decay in the entropies. We also
observe that third-order scheme decays slightly more entropy than the fourth-
order scheme, however, the difference is very small.

####  6.2.2 Relativistic two-fluid blast problem

(a) Plot of $\log_{10}(\rho_{i}+\rho_{e})$.

(b) Plot of $\log_{10}(p_{i}+p_{e})$.

(c) Ion Lorentz factor $\Gamma_{i}$.

(d) Magnitude of the Magnetic field $\dfrac{|\mathbf{B}|^{2}}{2}$.

Figure 7: Relativistic two-fluid blast problem: Plot for the weakly magnetized
medium $B_{0}=0.1$, using O3-ES-IMEX scheme with $200\times 200$ cells. We
have plotted 30 contours for each variable.

(a) Plot of $\log_{10}(\rho_{i}+\rho_{e})$.

(b) Plot of $\log_{10}(p_{i}+p_{e})$.

(c) Plot of ion Lorentz factor $\Gamma_{i}$.

(d) Plot of magnitude of the Magnetic field $\dfrac{|\mathbf{B}|^{2}}{2}$.

Figure 8: Relativistic two-fluid blast problem: Plot for the weakly magnetized
medium $B_{0}=0.1$, using O4-ES-IMEX scheme with $200\times 200$ cells. We
have plotted 30 contours for each variable.

(a) Plot of $\log_{10}(\rho_{i}+\rho_{e})$.

(b) Plot of $\log_{10}(p_{i}+p_{e})$.

(c) Plot of ion Lorentz factor $\Gamma_{i}$.

(d) Plot of magnitude of the Magnetic field $\dfrac{|\mathbf{B}|^{2}}{2}$.

Figure 9: Relativistic two-fluid blast problem: Plot for the strongly
magnetized medium $B_{0}=1.0$, using O3-ES-IMEX scheme with $200\times 200$
cells. We have plotted 30 contours for each variable.

(a) Plot of $\log_{10}(\rho_{i}+\rho_{e})$.

(b) Plot of $\log_{10}(p_{i}+p_{e})$.

(c) Plot of ion Lorentz factor $\Gamma_{i}$.

(d) Plot magnitude of the Magnetic field $\dfrac{|\mathbf{B}|^{2}}{2}$.

Figure 10: Relativistic two-fluid blast problem: Plot for the strongly
magnetized medium $B_{0}=1.0$, using O4-ES-IMEX scheme with $200\times 200$
cells. We have plotted 30 contours for each variable.

(a) Total entropy fluid decay for weakly magnetized medium, $B_{0}=0.1$.

(b) Total fluid entropy decay for strongly magnetized medium, $B_{0}=1.0$.

Figure 11: Relativistic two-fluid blast problem: Total fluid entropy
$\mathcal{U}_{i}+\mathcal{U}_{e}$ evolution of the schemes ES-O3-IMEX and
ES-O4-IMEX for the weakly and strongly magnetized medium.

We consider the two-fluid relativistic extension of the strong cylindrical
RMHD explosion described in [6]. We follow [12, 29] to consider a
computational domain of $[-6,6]\times[-6,6]$ with Neumann boundary conditions.
To describe the initial fluid profile, consider $\rho_{in}=10^{-2}$,
$p_{in}=1.0$, $\rho_{out}=10^{-4}$ and $p_{out}=5\times 10^{-4}$. Also, let us
denote the radial distance $r=\sqrt{x^{2}+y^{2}}$ from the center of the
computational domain. Inside the disc of radius 0.8, i.e., $r<0.8$, we define
the densities $\rho_{i}=\rho_{e}=0.5\times\rho_{in}$ and the pressures
$p_{i}=p_{e}=0.5\times p_{in}$. Also, outside the disc of radius 1.0, i.e.,
$r>1$, we define the densities $\rho_{i}=\rho_{e}=0.5\times\rho_{out}$ and the
pressures $p_{i}=p_{e}=0.5\times p_{out}$. In the range $0.8\leq r\leq 1.0$,
the densities and pressures are defined using a linear profile, such that the
densities and the pressures decreases with increase in the radius $r$ and at
the boundaries $r=0.8$ and $r=1.0$ they matches the prescribed constant
states. The magnetic field has been initialized in the $x$-direction only, as
$B_{x}=B_{0}$. We set all remaining variables to zero. The charge to mass
ratios are taken to be $r_{i}=-r_{e}=10^{3}$ and $\gamma_{i}=\gamma_{e}=4/3$.

We compute solutions on a mesh of $200\times 200$ cells using the O3-ES-IMEX
and O4-ES-IMEX schemes till final time $t=4.0.$ We consider two values of
parameter $B_{0}$, such that $B_{0}=0.1$ denoting the case of the weakly
magnetized medium, and $B_{0}=1.0$ denoting the case of the strongly
magnetized medium. In Figures (7) and (8), we plot the results for weakly
magnetized medium for O3-ES-IMEX and O4-ES-IMEX schemes, respectively. We have
plotted the log to the base 10 of total density and total pressure; we also
plot the ion Lorentz factor and magnitude of the magnetic field. We observe
that both the schemes produce stable results and are able to capture the waves
accurately. Furthermore, the results are comparable to those in [12] and [29].

Similarly, in the case of strongly magnetized medium we plot the results for
O3-ES-IMEX and O4-ES-IMEX schemes in Figures (9) and (10), respectively. We
observe a substantial change in the solution profile. Again both the schemes
are able to capture the blast waves accurately.

In Figure (11), we plot the evolution of the total fluid entropy
$\mathcal{U}_{i}+\mathcal{U}_{e}$ for the O3-ES-IMEX and O4-ES-IMEX schemes in
weakly and strongly magnetized mediums. As the solutions contain many
discontinuities, we see a strong decay in the entropies. Furthermore, O4-ES-
IMEX scheme decays less entropy than the O3-ES-IMEX scheme.

####  6.2.3 Relativistic two-fluid GEM challenge problem

In this test, we consider two-fluid relativistic Geospace Environment Modeling
(GEM) magnetic reconnection problem from [12], which is an extension of the
non-relativistic GEM magnetic reconnection given in [49]. We consider the
modified Equations (40a) and (40b) to capture the resistive effects. The
computational domain is $[-L_{x}/2,L_{x}/2]\times[-L_{y}/2,L_{y}/2]$ where
$L_{x}=8\pi$ and $L_{y}=4\pi$, with periodic boundary conditions at $x=\pm
L_{x}/2$ and conducting wall boundary at $y=\pm L_{y}/2$ boundary. The ion-
electron mass ratio is taken to be $m_{i}/m_{e}=25$ with $m_{i}=1$.
Accordingly, we have $r_{i}=1.0$ and $r_{e}=-25.0.$ The unperturbed
$x$-component of the magnetic field is taken to be $B_{x}(y)=B_{0}\tan(y/d)$
with $B_{0}=1.0$, where $d=1.0$ is the thickness of the current sheet. On the
other hand, unperturbed $y$\- and $z$-components of the magnetic field are
assumed to be zero. After the perturbation initial conditions become

$\displaystyle\begin{pmatrix}\rho_{i}\\\ u_{z_{i}}\\\ p_{i}\\\ \\\ \rho_{e}\\\
u_{z_{e}}\\\ p_{e}\\\ \\\ B_{x}\\\ B_{y}\end{pmatrix}=\begin{pmatrix}n\\\
\frac{c}{2d}\frac{B_{0}\text{sech}^{2}(y/d)}{n}\\\
0.2+\frac{B_{0}^{2}\text{sech}^{2}(y/d)}{4}\frac{5}{24\pi}\\\ \\\
\frac{m_{e}}{m_{p}}n\\\ -(u_{z})_{i}\\\ -p_{i}\\\ \\\
B_{0}\tan(y/d)-B_{0}\psi_{0}\frac{\pi}{L_{y}}\cos(\frac{\pi
x}{L_{x}})\sin(\frac{\pi y}{L_{y}})\\\
B_{0}\psi_{0}\frac{\pi}{L_{x}}\sin(\frac{\pi x}{L_{x}})\cos(\frac{\pi
y}{L_{y}})\end{pmatrix}$

where $n=\mathrm{sech}^{2}(y/d)+0.2$. All other variables are set to zero. We
set resistivity constant as $\eta=0.01$, and adiabatic indices as
$\gamma_{i}=\gamma_{e}=4.0/3.0$.

We plot total density, $z$-component of magnetic field, $x$-component of ion
velocity and $x$-component of electron velocity on a mesh of $512\times 256$
cells. The plots also contain field line for $(B_{x},B_{y})$. We plot the
solution at $t=40$ and $t=80$ for O3-ES-IMEX and O4-ES-IMEX. We observe that
both the solutions for both the schemes are similar and agrees with the
solutions presented in [12].

In Figure (16), we also plot time evolution of the reconnected magnetic flux,

$\psi(t)=\dfrac{1}{2B_{0}}\int_{-L_{x}/2}^{L_{x}/2}|B_{y}(x,y=0,t)|dx.$

at various resolutions. We also compare it with the results presented in [12].
We observe that both schemes are able to capture the reconnected flux at
various resolutions. Only for the time more than $t=80$, we see a difference
in the reconnection flux for the different schemes, but at this time, flow is
already quite turbulent.

(a) Total Density ($\rho_{i}+\rho_{e}$)

(b) $B_{z}$

(c) Ion x-velocity

(d) Electron x-velocity

Figure 12: Relativistic two-fluid GEM challenge problem: Plots for the total
density, $B_{z}$-component, Ion $x$-velocity, and Electron $x$-velocity on the
mess $512\times 256$, using scheme O3-ES-IMEX, at time t=40.0.

(a) Total Density ($\rho_{i}+\rho_{e}$)

(b) $B_{z}$

(c) Ion x-velocity

(d) Electron x-velocity

Figure 13: Relativistic two-fluid GEM challenge problem: Plots for the total
density, $B_{z}$-component, Ion $x$-velocity, and Electron $x$-velocity on the
mess $512\times 256$, using scheme O3-ES-IMEX, at time t=80.0.

(a) Total Density ($\rho_{i}+\rho_{e}$)

(b) $B_{z}$

(c) Ion x-velocity

(d) Electron x-velocity

Figure 14: Relativistic two-fluid GEM challenge problem: Plot for the total
density, $B_{z}$-component, Ion $x$-velocity, and Electron $x$-velocity on the
mess $512\times 256$, using scheme O4-ES-IMEX, at time t=40.0.

(a) Total Density ($\rho_{i}+\rho_{e}$)

(b) $B_{z}$

(c) Ion x-velocity

(d) Electron x-velocity

Figure 15: Relativistic two-fluid GEM challenge problem: Plot for the total
Density, $B_{z}$-component, Ion $x$-velocity, and Electron $x$-velocity on the
mess $512\times 256$, using scheme O4-ES-IMEX, at time t=80.0.  Figure 16:
Relativistic two-fluid GEM challenge problem: Time development of the
reconnected magnetic flux for $N_{x}=128,\ 256,\ 512$, using schemes O3-ES-
IMEX and O4-ES-IMEX. We overlay the plot on Amano’s GEM data [12] (solid
lines).

##  7 Conclusion

The two-fluid relativistic plasma flow equations are a multi-physics model
that couples relativistic hydrodynamic conservation laws with Maxwell’s
equations. This model is non-linear and exhibits non-smooth and multi-scale
phenomena. In this article, we have designed arbitrary high-order finite-
difference entropy stable schemes for the model. This was achieved by
exploiting the structure of the flux and by demonstrating that the source
terms do not contribute to the fluid entropies. Furthermore, we have also
presented ARK-IMEX schemes to overcome the time-step restriction imposed by
the stiff source terms. The implicit step in the ARK-IMEX scheme is efficient
as it involves solving only one set of nonlinear equations in each cell and
there is no need for any global nonlinear solver.

The proposed schemes are then tested on various test cases in one and two
dimensions. First, we demonstrate that the proposed schemes have formal order
of accuracy. Then we test the schemes on the Brio-Wu test problem, where we
have finite plasma skin depth. We note that the proposed schemes can resolve
the dispersive effects, and higher-order schemes can capture the finite skin
depth effects even on coarser meshes. We also demonstrate that on the coarser
meshes with the stiff source terms, IMEX schemes are highly efficient. The
SSP-RK and ARK-IMEX schemes have similar accuracy and entropy decay
performance. We then test the scheme on a current sheet problem where
resistive effects are considered. We observe that the proposed schemes
accurately capture the RMHD solution. We compute the Orzag-Tang vertex, blast
problem, and GEM challenge problem in two-dimensional test cases. We again
show that the schemes are able to capture solution features very well and
ensure entropy stability of the solutions. For the GEM problem, the schemes
are also able to predict magnetic reconnection flux.

## Acknowledgments

The work of Harish Kumar is supported by VAJRA grant No. VJR/2018/000129 by
the Dept. of Science and Technology, Govt. of India. The work of Praveen
Chandrashekar was supported by SERB-DST, Govt. of India, under the MATRICS
grant No. MTR/2018/000006, and by the Department of Atomic Energy, Government
of India, under project no. 12-R&D-TFR-5.01-0520. The authors would also like
to thank Prof. Dinshaw Balsara, Univ. of Notre-Dame, USA, for his several
suggestions.

## References

  * [1] D. Bhoriya, H. Kumar, Entropy-stable schemes for relativistic hydrodynamics equations, Zeitschrift fur Angewandte Mathematik und Physik 71 (1) (2020).  doi:10.1007/s00033-020-1250-8.
  * [2] Y. A. Gallant, J. Arons, Structure of relativistic shocks in pulsar winds: A model of the wisps in the Crab Nebula, The Astrophysical Journal 435 (1994) 230\.  doi:10.1086/174810.
  * [3] R. Mochkovitch, V. Maitia, R. Marques, Internal shocks in a relativistic wind as a source for gamma-ray bursts?, Astrophysics and Space Science 231 (1-2) (1995) 441–444.  doi:10.1007/BF00658666.
  * [4] L. D. Landau, E. M. Lifshitz, Relativistic Fluid Dynamics, in: Fluid Mechanics, Elsevier, 1987, pp. 505–514.  doi:10.1016/b978-0-08-033933-7.50023-4.
  * [5] J. F. Wardle, D. C. Homan, R. Ojha, D. H. Roberts, Electron-positron jets associated with the quasar 3C279, Nature 395 (6701) (1998) 457–461.  doi:10.1038/26675.
  * [6] S. S. Komissarov, A Godunov-type scheme for relativistic magnetohydrodynamics, Monthly Notices of the Royal Astronomical Society 303 (2) (1999) 343–366.  doi:10.1046/j.1365-8711.1999.02244.x.
  * [7] D. Balsara, Total Variation Diminishing Scheme for Relativistic Magnetohydrodynamics, The Astrophysical Journal Supplement Series 132 (1) (2001) 83–101.  doi:10.1086/318941.
  * [8] L. Del Zanna, N. Bucciantini, P. Londrillo, An efficient shock-capturing central-type scheme for multidimensional relativistic flows II. Magnetohydrodynamics, Astronomy and Astrophysics 400 (2) (2003) 397–413.  doi:10.1051/0004-6361:20021641.
  * [9] A. Mignone, G. Bodo, An HLLC Riemann solver for relativistic flows – II. Magnetohydrodynamics, Monthly Notices of the Royal Astronomical Society 368 (3) (2006) 1040–1054.  doi:10.1111/J.1365-2966.2006.10162.X.
  * [10] S. S. Komissarov, Multidimensional numerical scheme for resistive relativistic magnetohydrodynamics, Monthly Notices of the Royal Astronomical Society 382 (3) (2007) 995–1004.  doi:10.1111/j.1365-2966.2007.12448.x.
  * [11] D. S. Balsara, J. Kim, A subluminal relativistic magnetohydrodynamics scheme with ADER-WENO predictor and multidimensional Riemann solver-based corrector, Journal of Computational Physics 312 (2016) 357–384.  doi:10.1016/j.jcp.2016.02.001.
  * [12] T. Amano, A second-order divergence-constrained multidimensional numerical scheme for Relativistic Two-Fluid Electrodynamics, The Astrophysical Journal 831 (1) (2016) 100.  doi:10.3847/0004-637x/831/1/100.
  * [13] U. Shumlak, J. Loverich, Approximate Riemann solver for the two-fluid plasma model, Journal of Computational Physics 187 (2) (2003) 620–638.  doi:10.1016/S0021-9991(03)00151-7.
  * [14] A. Hakim, J. Loverich, U. Shumlak, A high resolution wave propagation scheme for ideal Two-Fluid plasma equations, Journal of Computational Physics 219 (1) (2006) 418–442.  doi:10.1016/j.jcp.2006.03.036.
  * [15] H. Kumar, S. Mishra, Entropy stable numerical schemes for two-fluid plasma equations, Journal of Scientific Computing 52 (2) (2012) 401–425.  doi:10.1007/s10915-011-9554-7.
  * [16] R. Abgrall, H. Kumar, Robust Finite Volume Schemes for Two-Fluid Plasma Equations, Journal of Scientific Computing 60 (3) (2014) 584–611.  doi:10.1007/s10915-013-9809-6.
  * [17] D. M. Bond, V. Wheatley, R. Samtaney, Plasma flow simulation using the two-fluid model, in: Proceedings of the 20th Australasian Fluid Mechanics Conference, 2016.
  * [18] Y. Li, R. Samtaney, D. Bond, V. Wheatley, Richtmyer-Meshkov instability of an imploding flow with a two-fluid plasma model, Physical Review Fluids 5 (11) (2020) 113701.  doi:10.1103/PhysRevFluids.5.113701.
  * [19] A. K. Meena, H. Kumar, Robust numerical schemes for Two-Fluid Ten-Moment plasma flow equations, Zeitschrift fur Angewandte Mathematik und Physik 70 (1) (2019) 1–30.  doi:10.1007/s00033-018-1061-3.
  * [20] S. Zenitani, M. Hesse, A. Klimas, Relativistic two-fluid simulations of guide field magnetic reconnection, Astrophysical Journal 705 (1) (2009) 907–913.  doi:10.1088/0004-637X/705/1/907.
  * [21] S. Zenitani, M. Hesse, A. Klimas, Resistive magnetohydrodynamic simulations of relativistic magnetic reconnection, Astrophysical Journal Letters 716 (2) (2010) L214–L218.  doi:10.1088/2041-8205/716/2/L214.
  * [22] T. Amano, J. G. Kirk, The role of superluminal electromagnetic waves in pulsar wind termination shocks, Astrophysical Journal 770 (1) (2013) 18.  doi:10.1088/0004-637X/770/1/18.
  * [23] M. Barkov, S. S. Komissarov, V. Korolev, A. Zankovich, A multidimensional numerical scheme for two-fluid relativistic magnetohydrodynamics, Monthly Notices of the Royal Astronomical Society 438 (1) (2014) 704–716.  doi:10.1093/mnras/stt2247.
  * [24] M. V. Barkov, S. S. Komissarov, Relativistic tearing and drift-kink instabilities in two-fluid simulations, Monthly Notices of the Royal Astronomical Society 458 (2) (2016) 1939–1947.  doi:10.1093/mnras/stw384.
  * [25] R. J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Cambridge University Press, 2002.  doi:10.1017/cbo9780511791253.
  * [26] E. Chiodaroli, C. De Lellis, O. Kreml, Global Ill-Posedness of the Isentropic System of Gas Dynamics, Communications on Pure and Applied Mathematics 68 (7) (2015) 1157–1190.  doi:10.1002/CPA.21537.
  * [27] E. Godlewski, P.-A. Raviart, Numerical Approximation of Hyperbolic Systems of Conservation Laws, Vol. 118 of Applied Mathematical Sciences, Springer New York, New York, NY, 1996.  doi:10.1007/978-1-4612-0713-9.
  * [28] S. Zenitani, M. Hesse, A. Klimas, Two-fluid magnetohydrodynamic simulations of relativistic magnetic reconnection, Astrophysical Journal 696 (2) (2009) 1385–1401.  doi:10.1088/0004-637X/696/2/1385.
  * [29] D. S. Balsara, T. Amano, S. Garain, J. Kim, A high-order relativistic two-fluid electrodynamic scheme with consistent reconstruction of electromagnetic fields and a multidimensional Riemann solver for electromagnetism, Journal of Computational Physics 318 (2016) 169–200.  doi:10.1016/j.jcp.2016.05.006.
  * [30] U. S. Fjordholm, S. Mishra, E. Tadmor, Arbitrarily high-order accurate entropy stable essentially nonoscillatory schemes for systems of conservation laws, SIAM Journal on Numerical Analysis 50 (2) (2012) 544–573.  doi:10.1137/110836961.
  * [31] U. S. Fjordholm, S. Mishra, E. Tadmor, ENO Reconstruction and ENO Interpolation Are Stable, Foundations of Computational Mathematics 13 (2) (2013) 139–159.  doi:10.1007/s10208-012-9117-9.
  * [32] P. Chandrashekar, Kinetic energy preserving and entropy stable finite volume schemes for compressible euler and Navier-Stokes equations, Communications in Computational Physics 14 (5) (2013) 1252–1286.  doi:10.4208/cicp.170712.010313a.
  * [33] C. Sen, H. Kumar, Entropy Stable Schemes For Ten-Moment Gaussian Closure Equations, Journal of Scientific Computing 75 (2) (2018) 1128–1155.  doi:10.1007/S10915-017-0579-4/FIGURES/8.
  * [34] J. Duan, H. Tang, High-order accurate entropy stable finite difference schemes for the shallow water magnetohydrodynamics, Journal of Computational Physics 431 (2021) 110136.  doi:10.1016/J.JCP.2021.110136.
  * [35] J. Duan, H. Tang, High-Order Accurate Entropy Stable Finite Difference Schemes for One- And Two-Dimensional Special Relativistic Hydrodynamics, Advances in Applied Mathematics and Mechanics 12 (1) (2020) 1–29.  doi:10.4208/AAMM.OA-2019-0124.
  * [36] B. Biswas, H. Kumar, D. Bhoriya, Entropy stable discontinuous Galerkin schemes for the special relativistic hydrodynamics equations, Computers & Mathematics with Applications 112 (September 2021) (2022) 55–75.  doi:10.1016/j.camwa.2022.02.019.
  * [37] C. D. Munz, P. Omnes, R. Schneider, E. Sonnendrücker, U. Voß, Divergence Correction Techniques for Maxwell Solvers Based on a Hyperbolic Model, Journal of Computational Physics 161 (2) (2000) 484–511.  doi:10.1006/jcph.2000.6507.
  * [38] V. Schneider, U. Katscher, D. H. Rischke, B. Waldhauser, J. A. Maruhn, C. D. Munz, New Algorithms for Ultra-relativistic Numerical Hydrodynamics, Journal of Computational Physics 105 (1) (1993) 92–107.  doi:10.1006/jcph.1993.1056.
  * [39] E. Tadmor, The Numerical Viscosity of Entropy Stable Schemes for Systems of Conservation Laws. I, Mathematics of Computation 49 (179) (1987) 91.  doi:10.2307/2008251.
  * [40] F. Ismail, P. L. Roe, Affordable, entropy-consistent Euler flux functions II: Entropy production at shocks, Journal of Computational Physics 228 (15) (2009) 5410–5436.  doi:10.1016/j.jcp.2009.04.021.
  * [41] P. G. Lefloch, J. M. Mercier, C. Rohde, Fully discrete, entropy conservative schemes of arbitrary order, SIAM Journal on Numerical Analysis 40 (5) (2002) 1968–1992.  doi:10.1137/S003614290240069X.
  * [42] T. J. Barth, Numerical methods for gas-dynamics systems on unstructured meshes, An introduction to recent developments in theory and numerics of conservation Laws, Lecture notes in computational science and engineering volume {5}, Springer, Berlin. Eds: D. Kroner, M., Springer, Berlin, Heidelberg, 1999.
  * [43] G. S. Jiang, C. W. Shu, Efficient implementation of weighted ENO schemes, Journal of Computational Physics 126 (1) (1996) 202–228.  doi:10.1006/jcph.1996.0130.
  * [44] S. Gottlieb, C. W. Shu, E. Tadmor, Strong stability-preserving high-order time discretization methods, SIAM Review 43 (1) (2001) 89–112.  doi:10.1137/S003614450036757X.
  * [45] L. Pareschi, G. Russo, Implicit-explicit Runge-Kutta schemes and applications to hyperbolic systems with relaxation, Journal of Scientific Computing 25 (1) (2005) 129–155.  doi:10.1007/s10915-004-4636-4.
  * [46] J. E. Dennis, R. B. Schnabel, Numerical Methods for Unconstrained Optimization and Nonlinear Equations, Society for Industrial and Applied Mathematics, 1996\.  doi:10.1137/1.9781611971200.
  * [47] C. A. Kennedy, M. H. Carpenter, Additive Runge-Kutta schemes for convection-diffusion-reaction equations, Applied Numerical Mathematics 44 (1-2) (2003) 139–181.  doi:10.1016/S0168-9274(02)00138-1.
  * [48] S. A. Orszag, C. M. Tang, Small-scale structure of two-dimensional magnetohydrodynamic turbulence, Journal of Fluid Mechanics 90 (1) (1979) 129–143.  doi:10.1017/S002211207900210X.
  * [49] J. Birn, J. F. Drake, M. A. Shay, B. N. Rogers, R. E. Denton, M. Hesse, M. Kuznetsova, Z. W. Ma, A. Bhattacharjee, A. Otto, P. L. Pritchett, Geospace Environmental Modeling (GEM) Magnetic Reconnection Challenge, Journal of Geophysical Research: Space Physics 106 (A3) (2001) 3715–3719.  doi:10.1029/1999ja900449.

##  Appendix A Barth scaling of right eigenvectors

Here, we present expressions for the right eigenvectors and scaled right
eigenvectors for the conservative system (4). The right eigenvectors are
described in Section A.1, while the expressions and procedure to obtain the
scaled eigenvectors is presented in Section A.2

###  A.1 Right eigenvectors

The set of right eigenvectors of the matrix
$\mathbf{A}^{x}=\frac{\partial\mathbf{f}^{x}}{\partial\mathbf{U}}$
corresponding to the eigenvalues $\Lambda^{x}$ of the system (4) is obtained
by taking $d=x$ in the ordered set
$\mathbf{R}^{d}_{\Lambda^{d}}=\left\\{\left(\mathbf{R}_{\Lambda^{d}}^{d}\right)_{n}:\
n=1,\ 2,\ 3,\dots,\ 16\right\\}$ where the vectors
$\left(\mathbf{R}_{\Lambda^{d}}^{d}\right)_{n}$ are defined as

$\displaystyle\left(\mathbf{R}_{\Lambda^{d}}^{d}\right)_{n}=\begin{cases}\Big{(}(\mathbf{R}_{i,k}^{d})_{1\times
5},\mathbf{0}_{1\times 5},\mathbf{0}_{1\times 6}\Big{)}^{\top},&1\leq n\leq
5,\ k=n\\\ \Big{(}\mathbf{0}_{1\times 5},(\mathbf{R}_{e,k}^{d})_{1\times
5},\mathbf{0}_{1\times 6}\Big{)}^{\top},&6\leq n\leq 10,\ k=n-5\\\
\Big{(}\mathbf{0}_{1\times 5},\mathbf{0}_{1\times
5},(\mathbf{R}_{m,k}^{d})_{1\times 6}\Big{)}^{\top},&11\leq n\leq 16,\
k=n-10,\end{cases}$ (44)

where $\mathbf{R}^{d}_{\alpha,k}$, $\alpha\in\\{i,e\\}$, is the $k^{th}$
column vector of the $5\times 5$ right eigenvector matrices
$\mathbf{R}^{d}_{\alpha}$ of the flux jacobians
$\dfrac{\partial\mathbf{f}^{d}_{\alpha}}{\partial\mathbf{U}_{\alpha}}$, and
$\mathbf{R}^{d}_{m,k}$ is the $k^{th}$ column vector of the right eigenvector
matrix $\mathbf{R}^{d}_{m}$ of the flux jacobian matrix
$\dfrac{\partial\mathbf{f}^{d}_{m}}{\partial\mathbf{U}_{m}}$. The matrices
$\mathbf{R}^{k}_{\alpha}$ and $\mathbf{R}^{k}_{m}$ have the following
expressions.

  * •

For $d=x,y$, and $\alpha\in\\{i,e\\}$, the right eigenvector matrix
$\mathbf{R}^{d}_{\alpha}$ is given by the relation,

$\mathbf{R}_{\alpha}^{d}=\left(\dfrac{\partial\mathbf{U}_{\alpha}}{\partial\mathbf{W}_{\alpha}}\right)\mathbf{R}_{\alpha,\mathbf{W}}^{d},$
(45)

where the matrix $\mathbf{R}_{\alpha,\mathbf{W}}^{d}$ for $d=x$ is given by,

$\displaystyle\mathbf{R}_{\alpha,\mathbf{W}}^{x}=\begin{pmatrix}\frac{1}{c_{\alpha}^{2}h_{\alpha}}&1&0&0&\frac{1}{c_{\alpha}^{2}h_{\alpha}}\\\
\frac{-\sqrt{Q^{x}_{\alpha}}}{c_{\alpha}h_{\alpha}\Gamma_{\alpha}\rho_{\alpha}}&0&0&0&\frac{+\sqrt{Q^{x}_{\alpha}}}{c_{\alpha}h_{\alpha}\Gamma_{\alpha}\rho_{\alpha}}\\\
\frac{\left(c_{\alpha}-\Gamma_{\alpha}\sqrt{Q^{x}_{\alpha}}{u_{x_{\alpha}}}\right){u_{y_{\alpha}}}}{c_{\alpha}h_{\alpha}\Gamma^{2}_{\alpha}\rho_{\alpha}\left(u_{x_{\alpha}}^{2}-1\right)}&0&1&0&\frac{\left(c_{\alpha}+\Gamma_{\alpha}\sqrt{Q^{x}_{\alpha}}{u_{x_{\alpha}}}\right){u_{y_{\alpha}}}}{c_{\alpha}h_{\alpha}\Gamma^{2}_{\alpha}\rho_{\alpha}\left(u_{x_{\alpha}}^{2}-1\right)}\\\
\frac{\left(c_{\alpha}-\Gamma_{\alpha}\sqrt{Q^{x}_{\alpha}}{u_{x_{\alpha}}}\right){u_{z_{\alpha}}}}{c_{\alpha}h_{\alpha}\Gamma^{2}_{\alpha}\rho_{\alpha}\left(u_{x_{\alpha}}^{2}-1\right)}&0&1&0&\frac{\left(c_{\alpha}+\Gamma_{\alpha}\sqrt{Q^{x}_{\alpha}}{u_{x_{\alpha}}}\right){u_{z_{\alpha}}}}{c_{\alpha}h_{\alpha}\Gamma^{2}_{\alpha}\rho_{\alpha}\left(u_{x_{\alpha}}^{2}-1\right)}\\\
1&0&0&0&1\end{pmatrix}.$ (46)

  * •

The eigenvector matrix $\mathbf{R}^{d}_{m}$ for $d=x$ is given by,

$\displaystyle\mathbf{R}_{m}^{x}=\begin{pmatrix}0&0&1&0&0&0\\\ 0&1&0&0&0&-1\\\
-1&0&0&0&1&0\\\ 0&0&0&1&0&0\\\ 1&0&0&0&1&0\\\ 0&1&0&0&0&1\end{pmatrix}.$

######  Remark A.1.

For the y-directional flux, $\mathbf{f}^{y}$, we proceed similarly to get the
set of eigenvalues $\Lambda^{y}$ of the jacobian matrix
$\dfrac{\partial\mathbf{f}^{y}}{\partial\mathbf{U}}$ as,

$\displaystyle\Lambda^{y}=\biggl{\\{}$
$\displaystyle\frac{(1-c_{i}^{2})u_{y_{i}}-(c_{i}/\Gamma_{i})\sqrt{Q_{i}^{y}}}{1-c_{i}^{2}|\mathbf{u}_{i}|^{2}},\
u_{y_{i}},\ u_{y_{i}},\
u_{y_{i}},\frac{(1-c_{i}^{2})u_{y_{i}}+(c_{i}/\Gamma_{i})\sqrt{Q_{i}^{y}}}{1-c_{i}^{2}|\mathbf{u}_{i}|^{2}},$
$\displaystyle\frac{(1-c_{e}^{2})u_{y_{e}}-(c_{e}/\Gamma_{e})\sqrt{Q_{e}^{y}}}{1-c_{e}^{2}|\mathbf{u}_{e}|^{2}},\
u_{y_{e}},\ u_{y_{e}},\
u_{y_{e}},\frac{(1-c_{e}^{2})u_{y_{e}}+(c_{e}/\Gamma_{e})\sqrt{Q_{e}^{y}}}{1-c_{e}^{2}|\mathbf{u}_{e}|^{2}},$
$\displaystyle-1,\ -1,\ 0,\ 0,\ 1,\ 1\ \biggr{\\}},$

where,
$Q_{\alpha}^{y}=1-u_{y_{\alpha}}^{2}-c_{\alpha}^{2}(u_{x_{\alpha}}^{2}+u_{z_{\alpha}}^{2}),\
\alpha\in\\{i,\ e\\}$. The corresponding right eigenvectors are given by
fixing $d=y$ in the ordered set
$\mathbf{R}^{d}_{\Lambda^{d}}=\left\\{\left(\mathbf{R}_{\Lambda^{d}}^{d}\right)_{n}:\
n=1,\ 2,\ 3,\dots,\ 16\right\\}$ where
$\left(\mathbf{R}_{\Lambda^{d}}^{d}\right)_{n}$ is defined by Eqn. (44) and
the updated $y$-directional matrices $\mathbf{R}_{\alpha,\mathbf{W}}^{y}$ and
$\mathbf{R}^{y}_{m}$ are given by

$\displaystyle\mathbf{R}_{\alpha,\mathbf{W}}^{y}=\begin{pmatrix}\frac{1}{c_{\alpha}^{2}h_{\alpha}}&1&0&0&\frac{1}{c_{\alpha}^{2}h_{\alpha}}\\\
\frac{\left(c_{\alpha}-\Gamma_{\alpha}\sqrt{Q^{y}_{\alpha}}{u_{y_{\alpha}}}\right){u_{x_{\alpha}}}}{c_{\alpha}h_{\alpha}\Gamma^{2}_{\alpha}\rho_{\alpha}\left(u_{y_{\alpha}}^{2}-1\right)}&0&1&0&\frac{\left(c_{\alpha}+\Gamma_{\alpha}\sqrt{Q^{y}_{\alpha}}{u_{y_{\alpha}}}\right){u_{x_{\alpha}}}}{c_{\alpha}h_{\alpha}\Gamma^{2}_{\alpha}\rho_{\alpha}\left(u_{y_{\alpha}}^{2}-1\right)}\\\
\frac{-\sqrt{Q^{y}_{\alpha}}}{c_{\alpha}h_{\alpha}\Gamma_{\alpha}\rho_{\alpha}}&0&0&0&\frac{+\sqrt{Q^{y}_{\alpha}}}{c_{\alpha}h_{\alpha}\Gamma_{\alpha}\rho_{\alpha}}\\\
\frac{\left(c_{\alpha}-\Gamma_{\alpha}\sqrt{Q^{y}_{\alpha}}{u_{y_{\alpha}}}\right){u_{z_{\alpha}}}}{c_{\alpha}h_{\alpha}\Gamma^{2}_{\alpha}\rho_{\alpha}\left(u_{y_{\alpha}}^{2}-1\right)}&0&0&1&\frac{\left(c_{\alpha}+\Gamma_{\alpha}\sqrt{Q^{y}_{\alpha}}{u_{y_{\alpha}}}\right){u_{z_{\alpha}}}}{c_{\alpha}h_{\alpha}\Gamma^{2}_{\alpha}\rho_{\alpha}\left(u_{y_{\alpha}}^{2}-1\right)}\\\
1&0&0&0&1\end{pmatrix},$ (47)
$\displaystyle\mathbf{R}_{m}^{y}=\begin{pmatrix}0&-1&0&0&0&1\\\ 0&0&1&0&0&0\\\
1&0&0&0&-1&0\\\ 1&0&0&0&1&0\\\ 0&0&0&1&0&0\\\ 0&1&0&0&0&1\end{pmatrix}.$ (48)

###  A.2 Barth scaling and entropy scaled right eigenvectors

We ignore the superscript $x$ to derive the entropy scaled right eigenvectors
for the $x$-directional part; for the $y$-directional flux, we proceed
similarly. We have the entropy variable vector $\mathbf{V}_{\alpha}$ as

$\mathbf{V}_{\alpha}=\begin{pmatrix}\dfrac{\gamma_{\alpha}-s_{\alpha}}{\gamma_{\alpha}-1}+{\beta_{\alpha}}\\\
{u_{x_{\alpha}}\Gamma_{\alpha}\beta_{\alpha}}\\\
{u_{y_{\alpha}}\Gamma_{\alpha}\beta_{\alpha}}\\\
{u_{z_{\alpha}}\Gamma_{\alpha}\beta_{\alpha}}\\\
-{\Gamma_{\alpha}\beta_{\alpha}}\end{pmatrix},\qquad\text{ where
$\beta_{\alpha}=\frac{\rho_{\alpha}}{p_{\alpha}}$.}$

We want to find scaling matrices $\mathbf{T}^{d}_{\alpha}$
($d\in\\{x,y\\}\text{ and }\alpha\in\\{i,e\\}$) such that the scaled right
eigenvector matrices
${\mathbf{\tilde{R}}}^{d}_{\alpha}=\mathbf{R}^{d}_{\alpha}\mathbf{T}^{d}_{\alpha}$,
with $\mathbf{R}^{d}_{\alpha}$ as in Eqn. (45), satisfy

$\frac{\partial\mathbf{U}_{\alpha}}{\partial\mathbf{V}_{\alpha}}={\mathbf{\tilde{R}}^{d}_{\alpha}}({\mathbf{\tilde{R}}^{d}_{\alpha}})^{\top}.$

The matrices $\mathbf{T}^{d}_{\alpha}$ are known as Barth scaling matrices.
Following the procedure of [42], the scaling matrices
$\mathbf{T}^{d}_{\alpha}$ are given by the square root of the matrices
$\mathbf{Y}^{d}_{\alpha}$ where the expression for matrices
$\mathbf{Y}^{d}_{\alpha}$ is given by the formulae

$\mathbf{Y}^{d}_{\alpha}=({\mathbf{R}^{d}_{\alpha,\mathbf{W}}})^{-1}\frac{\partial\mathbf{W}_{\alpha}}{\partial\mathbf{V}_{\alpha}}\frac{\partial\mathbf{U}_{\alpha}}{\partial\mathbf{W}_{\alpha}}^{-\top}{\mathbf{R}^{d}_{\alpha,\mathbf{W}}}^{-\top},$

where the matrices $\mathbf{R}_{\alpha,\mathbf{W}}^{d}$ are given by Eqn. (46)
for $d=x$, and for the case $d=y$ the expressions for the matrices follows
from Eqn. (47). Assuming $d=x$ for the clarity. A long simplification leads to
the expressions

$\mathbf{Y}^{x}_{\alpha}=\begin{pmatrix}\frac{c_{\alpha}^{2}h_{\alpha}p_{\alpha}\left(1+\frac{c_{\alpha}{u_{x_{\alpha}}}}{\Gamma_{\alpha}\sqrt{Q^{x}_{\alpha}}}\right)}{2\Gamma_{\alpha}\left(1-c_{\alpha}^{2}|\mathbf{u}|_{\alpha}^{2}\right)}&0&0&0&0\\\
0&\frac{\rho_{\alpha}}{\Gamma_{\alpha}}\frac{\gamma_{\alpha}-1}{\gamma_{\alpha}}&0&0&0\\\
0&0&Y^{x}_{\alpha_{33}}&Y^{x}_{\alpha_{34}}&0\\\
0&0&Y^{x}_{\alpha_{43}}&Y^{x}_{\alpha_{44}}&0\\\
0&0&0&0&\frac{c_{\alpha}^{2}h_{\alpha}p_{\alpha}\left(1-\frac{c_{\alpha}{u_{x_{\alpha}}}}{\Gamma_{\alpha}\sqrt{Q^{x}_{\alpha}}}\right)}{2\Gamma_{\alpha}\left(1-c_{\alpha}^{2}|\mathbf{u}|_{\alpha}^{2}\right)}\end{pmatrix}$

where,
$Q_{\alpha}^{x}=1-u_{x_{\alpha}}^{2}-c_{\alpha}^{2}(u_{y_{\alpha}}^{2}+u_{z_{\alpha}}^{2})$,
$Y^{x}_{\alpha_{33}}=C^{x}_{\alpha}{(1-u_{x_{\alpha}}^{2}-u_{y_{\alpha}}^{2})}$,
$Y^{x}_{\alpha_{44}}=C^{x}_{\alpha}{(1-u_{x_{\alpha}}^{2}-u_{z_{\alpha}}^{2})}$,
$Y^{x}_{\alpha_{34}}=Y^{x}_{\alpha_{43}}=-C^{x}_{\alpha}u_{y_{\alpha}}u_{z_{\alpha}}$
and

$C^{x}_{\alpha}=\frac{p_{\alpha}}{h_{\alpha}\Gamma_{\alpha}^{3}\rho_{\alpha}^{2}\left(1-u_{x_{\alpha}}^{2}\right)}.$

Consequently, the matrix
$\mathbf{T}^{x}_{\alpha}=\sqrt{\mathbf{Y}^{x}_{\alpha}}$ is given by

$\mathbf{T}^{x}_{\alpha}=\begin{pmatrix}\sqrt{\frac{c_{\alpha}^{2}h_{\alpha}p_{\alpha}\left(1+\frac{c_{\alpha}{u_{x_{\alpha}}}}{\Gamma_{\alpha}\sqrt{Q^{x}_{\alpha}}}\right)}{2\Gamma_{\alpha}\left(1-c_{\alpha}^{2}|\mathbf{u}|_{\alpha}^{2}\right)}}&0&0&0&0\\\
0&\sqrt{\frac{\rho_{\alpha}}{\Gamma_{\alpha}}\frac{\gamma_{\alpha}-1}{\gamma_{\alpha}}}&0&0&0\\\
0&0&T^{x}_{\alpha_{33}}&T^{x}_{\alpha_{34}}&0\\\
0&0&T^{x}_{\alpha_{43}}&T^{x}_{\alpha_{44}}&0\\\
0&0&0&0&\sqrt{\frac{c_{\alpha}^{2}h_{\alpha}p_{\alpha}\left(1-\frac{c_{\alpha}{u_{x_{\alpha}}}}{\Gamma_{\alpha}\sqrt{Q^{x}_{\alpha}}}\right)}{2\Gamma_{\alpha}\left(1-c_{\alpha}^{2}|\mathbf{u}|_{\alpha}^{2}\right)}}\end{pmatrix}$

where,

$\scriptstyle\begin{pmatrix}T^{x}_{\alpha_{33}}&T^{x}_{\alpha_{34}}\\\
T^{x}_{\alpha_{43}}&T^{x}_{\alpha_{44}}\end{pmatrix}=\scriptstyle\begin{cases}\begin{pmatrix}\sqrt{Y^{x}_{\alpha_{33}}}&0\\\
0&\sqrt{Y^{x}_{\alpha_{44}}}\end{pmatrix}&\text{ if
}u_{y_{\alpha}}=u_{z_{\alpha}}=0,\\\ \\\
\sqrt{C^{x}_{\alpha}}\scriptstyle\begin{pmatrix}\frac{u_{z_{\alpha}}^{2}\sqrt{1-u_{x_{\alpha}}^{2}}+\frac{u_{y_{\alpha}}^{2}}{\Gamma_{\alpha}}}{u_{y_{\alpha}}^{2}+u_{z_{\alpha}}^{2}}&\frac{u_{y_{\alpha}}u_{z_{\alpha}}}{u_{y_{\alpha}}^{2}+u_{z_{\alpha}}^{2}}\left(\frac{1}{\Gamma_{\alpha}}-\sqrt{1-u_{x_{\alpha}}^{2}}\right)\\\
\frac{u_{y_{\alpha}}u_{z_{\alpha}}}{u_{y_{\alpha}}^{2}+u_{z_{\alpha}}^{2}}\left(\frac{1}{\Gamma_{\alpha}}-\sqrt{1-u_{x_{\alpha}}^{2}}\right)&\frac{u_{y_{\alpha}}^{2}\sqrt{1-u_{x_{\alpha}}^{2}}+\frac{u_{z_{\alpha}}^{2}}{\Gamma_{\alpha}}}{u_{y_{\alpha}}^{2}+u_{z_{\alpha}}^{2}}\end{pmatrix}&\text{
otherwise.}\end{cases}$

######  Remark A.2.

Proceeding similarly in the y-direction with the matrix
$\mathbf{R}_{\alpha,\mathbf{W}}^{y}$ as in Eqn. (47), the scaling matrix is
given by

$\mathbf{T}^{y}_{\alpha}=\begin{pmatrix}\sqrt{\frac{c_{\alpha}^{2}h_{\alpha}p_{\alpha}\left(1+\frac{c_{\alpha}{u_{y_{\alpha}}}}{\Gamma_{\alpha}\sqrt{Q^{y}_{\alpha}}}\right)}{2\Gamma_{\alpha}\left(1-c_{\alpha}^{2}|\mathbf{u}|_{\alpha}^{2}\right)}}&0&0&0&0\\\
0&\sqrt{\frac{\rho_{\alpha}}{\Gamma_{\alpha}}\frac{\gamma_{\alpha}-1}{\gamma_{\alpha}}}&0&0&0\\\
0&0&T^{y}_{\alpha_{33}}&T^{y}_{\alpha_{34}}&0\\\
0&0&T^{y}_{\alpha_{43}}&T^{y}_{\alpha_{44}}&0\\\
0&0&0&0&\sqrt{\frac{c_{\alpha}^{2}h_{\alpha}p_{\alpha}\left(1-\frac{c_{\alpha}{u_{y_{\alpha}}}}{\Gamma_{\alpha}\sqrt{Q^{y}_{\alpha}}}\right)}{2\Gamma_{\alpha}\left(1-c_{\alpha}^{2}|\mathbf{u}|_{\alpha}^{2}\right)}}\end{pmatrix}$

where,
$Q_{\alpha}^{y}=1-u_{y_{\alpha}}^{2}-c_{\alpha}^{2}(u_{x_{\alpha}}^{2}+u_{z_{\alpha}}^{2})$,
and the block entries are given by the matrices

$\scriptstyle\begin{pmatrix}T^{y}_{\alpha_{33}}&T^{y}_{\alpha_{34}}\\\
T^{y}_{\alpha_{43}}&T^{y}_{\alpha_{44}}\end{pmatrix}=\scriptstyle\begin{cases}\begin{pmatrix}\sqrt{Y^{y}_{\alpha_{33}}}&0\\\
0&\sqrt{Y^{y}_{\alpha_{44}}}\end{pmatrix}&\text{ if
}u_{x_{\alpha}}=u_{z_{\alpha}}=0,\\\ \\\
\sqrt{C^{y}_{\alpha}}\scriptstyle\begin{pmatrix}\frac{u_{z_{\alpha}}^{2}\sqrt{1-u_{y_{\alpha}}^{2}}+\frac{u_{x_{\alpha}}^{2}}{\Gamma_{\alpha}}}{u_{x_{\alpha}}^{2}+u_{z_{\alpha}}^{2}}&\frac{u_{x_{\alpha}}u_{z_{\alpha}}}{u_{x_{\alpha}}^{2}+u_{z_{\alpha}}^{2}}\left(\frac{1}{\Gamma_{\alpha}}-\sqrt{1-u_{y_{\alpha}}^{2}}\right)\\\
\frac{u_{x_{\alpha}}u_{z_{\alpha}}}{u_{x_{\alpha}}^{2}+u_{z_{\alpha}}^{2}}\left(\frac{1}{\Gamma_{\alpha}}-\sqrt{1-u_{y_{\alpha}}^{2}}\right)&\frac{u_{x_{\alpha}}^{2}\sqrt{1-u_{y_{\alpha}}^{2}}+\frac{u_{z_{\alpha}}^{2}}{\Gamma_{\alpha}}}{u_{x_{\alpha}}^{2}+u_{z_{\alpha}}^{2}}\end{pmatrix}&\text{
otherwise,}\end{cases}$

with
$Y^{y}_{\alpha_{33}}=C^{y}_{\alpha}{(1-u_{y_{\alpha}}^{2}-u_{x_{\alpha}}^{2})}$,
$Y^{y}_{\alpha_{44}}=C^{y}_{\alpha}{(1-u_{y_{\alpha}}^{2}-u_{z_{\alpha}}^{2})}$
and

$C^{y}_{\alpha}=\frac{p_{\alpha}}{h_{\alpha}\Gamma_{\alpha}^{3}\rho_{\alpha}^{2}\left(1-u_{y_{\alpha}}^{2}\right)}.$

##  Appendix B ARK3-IMEX coefficients table

The coefficients for the third order ARK-IMEX time update are given as
follows,

|  | $l=0$ | $l=1$ | $l=2$ | $l=3$
---|---|---|---|---|---
| $m=1$ | $\frac{1767732205903}{2027836641118}$ | 0 | 0 | 0
$a_{ml}^{[NS]}$ | $m=2$ | $\frac{5535828885825}{10492691773637}$ | $\frac{788022342437}{10882634858940}$ | 0 | 0
| $m=3$ | $\frac{6485989280629}{16251701735622}$ | $-\frac{4246266847089}{9704473918619}$ | $\frac{10755448449292}{10357097424841}$ | 0
| $m=1$ | $\frac{1767732205903}{4055673282236}$ | $\frac{1767732205903}{4055673282236}$ | 0 | 0
$a_{ml}^{[S]}$ | $m=2$ | $\frac{2746238789719}{10658868560708}$ | $-\frac{640167445237}{6845629431997}$ | $\frac{1767732205903}{4055673282236}$ | 0
| $m=3$ | $\frac{1471266399579}{7840856788654}$ | $-\frac{4482444167858}{7529755066697}$ | $\frac{11266239266428}{11593286722821}$ | $\frac{1767732205903}{4055673282236}$
| $b_{l}^{[NS]}$ | $\frac{1471266399579}{7840856788654}$ | $-\frac{4482444167858}{7529755066697}$ | $\frac{11266239266428}{11593286722821}$ | $\frac{1767732205903}{4055673282236}$
| $b_{l}^{[S]}$ | $\frac{2756255671327}{12835298489170}$ | $-\frac{10771552573575}{22201958757719}$ | $\frac{9247589265047}{10645013368117}$ | $\frac{2193209047091}{5459859503100}$
| $c_{l}$ | 0 | $\frac{1767732205903}{2027836641118}$ | $\frac{3}{5}$ | 1
Table 5: The coefficients for ARK3-IMEX time update scheme.