∎

11institutetext: Chetan Singh 22institutetext: Department of Mathematics,
Indian Institute of Technology Delhi, India
22email<EMAIL_ADDRESS>33institutetext: Anshu Yadav
44institutetext: Department of Mathematics, Indian Institute of Technology
Delhi, India
44email<EMAIL_ADDRESS>55institutetext: Deepak Bhoriya
66institutetext: Physics Department, University of Notre Dame, USA
66email<EMAIL_ADDRESS>77institutetext: Harish Kumar 88institutetext:
Department of Mathematics, Indian Institute of Technology Delhi, India
88email<EMAIL_ADDRESS>99institutetext: Dinshaw S. Balsara
1010institutetext: Physics Department, University of Notre Dame, USA
ACMS, University of Notre Dame, USA
1010email<EMAIL_ADDRESS>
# Entropy stable finite difference schemes for Chew, Goldberger & Low
anisotropic plasma flow equations

Chetan Singh  Anshu Yadav111Corresponding Author Deepak Bhoriya  Harish Kumar
Dinshaw S. Balsara

(Received: date / Accepted: date)

###### Abstract

In this article, we consider the Chew, Goldberger & Low (CGL) plasma flow
equations, which is a set of nonlinear, non-conservative hyperbolic PDEs
modelling anisotropic plasma flows. These equations incorporate the double
adiabatic approximation for the evolution of the pressure, making them very
valuable for plasma physics, space physics and astrophysical applications. We
first present the entropy analysis for the weak solutions. We then propose
entropy-stable finite-difference schemes for the CGL equations. The key idea
is to rewrite the CGL equations such that the non-conservative terms do not
contribute to the entropy equations. The conservative part of the rewritten
equations is very similar to the magnetohydrodynamics (MHD) equations. We then
symmetrize the conservative part by following Godunov’s symmetrization process
for MHD. The resulting equations are then discretized by designing entropy
conservative numerical flux and entropy diffusion operator based on the
entropy scaled eigenvectors of the conservative part. We then prove the semi-
discrete entropy stability of the schemes for CGL equations. The schemes are
then tested using several test problems derived from the corresponding MHD
test cases.

###### Keywords:

CGL Equations Symmetrization Entropy stability Higher-order finite difference
schemes

###### MSC:

65M06 65M12 65M20 76P99 76W05

##  1 Introduction

Plasma flows play an important role in numerous applications in astrophysics,
e.g., relativistic jets from active galactic nuclei landau1987relativistic ,
earth’s magnetosphere heinemann1999role ; dumin2002corotation ;
karimabadi2014link ; burch2016magnetospheric , the solar wind
hollweg1976collisionless ; zouganelis2004transonic ; chandran2011incorporating
. Often, to simulate plasma flows, equations of Magnetohydrodynamics (MHD) are
used. It is one of the simplest descriptions of the plasma flows as it treats
plasma as a single fluid having a single density, velocity and pressure.
However, the derivation of MHD equations involves several assumptions, which
may not be true in many cases. Hence, several efforts have been made to
generalize the MHD equations (otto19903d ; zhao_positivity-preserving_2014 ;
hakim_high_2006 ; kumar_entropy_2012 ; abgrall_robust_2014 ).

One of the key assumptions in the derivation of MHD equations is the
assumption of fluid being local thermodynamical equilibrium, which allows the
scalar description of the pressure (or temperature) (see hunana2019brief ;
goedbloed2004principles ). However, for space plasmas, this assumption is
often not true meena2019robust ; brown1995numerical ; dubroca2004magnetic ;
hakim2008extended ; johnson2014gaussian ; sangam2008hllc ; meng2012classical ;
Hirabayashi2016new ; wang_exact_2020  and a tensorial description of the
pressure is needed adequately describe non-local thermodynamical equilibrium
effects. In berthon2006numerical ; meena_robust_2017 ; meena_positivity-
preserving_2017 ; meena_positivity-preserving_2017 ; meena_well-balanced_2017
; biswas2021entropy ; sen_entropy_2018 ; Hirabayashi2016new , a three-
dimensional tensor is used to describe the pressure. However, pressure tensor
in plasmas is often rotated by magnetic field (see meena2019robust ;
hakim2008extended ; sangam2008hllc ) and in this case, pressure tensor can be
simplified by using tensorial components parallel to magnetic field and
perpendicular to the magnetic field (huang2019six ; hunana2019brief ;
meng2012classical ). Furthermore, if two-fluid effects are not important, then
plasma flow can be described using Chew, Goldberger & Low (CGL) equations
(chew1956boltzmann ; gombosi1991transport ; hunana2019brief ;
meng2012classical ).

The CGL equations were first proposed in chew1956boltzmann  and are the
simplest model to describe anisotropic plasma flows. These equations apply to
dilute plasmas where the plasma may be approximated as collisionless. In the
collisionless limit, charged particles gyrating around a magnetic field have
to satisfy certain invariants. The CGL equations are built on the idea that
there are two adiabatic invariants, one for the pressure parallel to the
magnetic field and another for the pressure perpendicular to the magnetic
field. In addition to the mass and momentum conservation equations, the fluid
is described using two scalar pressure components (one parallel to the
magnetic field and another perpendicular to the magnetic field). Furthermore,
the evolution of the magnetic field is modelled via an induction equation. The
CGL equations are a set of non-conservative hyperbolic systems of PDEs
(meng2012classical ; huang2019six ). Due to this, the design of the numerical
schemes for CGL equations is challenging. The weak solution formulations for
these equations require choosing an appropriate path which is often not known
dal1995definition . Furthermore, the numerical solutions are sensitive to the
numerical viscosity abgrall2010comment . In practice, a linear path is often
considered to design approximated Riemann solvers (huang2019six ;
meng2012classical ).

In this article, we consider entropy stable discretization of the equations
following the ideas in chandrashekar2013kinetic ; fjordholm2012arbitrarily ;
ismail2009affordable ; tadmor2003entropy ; yadav2023entropy ;
chandrashekar2016entropy . We proceed as follows:

  1. 1.

We first present the entropy analysis of the CGL system at the continuous
level.

  2. 2.

Following the idea of yadav2023entropy , we then reformulate the CGL equations
such the non-conservative terms do not contribute to the entropy evolution.

  3. 3.

The conservative flux of the reformulated equations is very similar to the MHD
flux. Furthermore, the entropy equation for the conservative part depends on
the divergence of the magnetic field. Hence, we use Godunov’s symmetrization
for the conservative equations following chandrashekar2016entropy .

  4. 4.

We then design entropy conservative numerical flux for the conservative part
of equations. Following fjordholm2012arbitrarily ; we also design a higher-
order numerical diffusion operator based on the entropy-scaled eigenvectors
and sign-preserving reconstruction process.

  5. 5.

Finally, combining the entropy-stable discretization of the conservative part
with a non-conservative part, we prove the semi-discrete entropy stability of
the proposed scheme.

The rest of the article is organized as follows: In the next Section, we
describe the CGL equations and present an entropy analysis of the equations.
In Section 3, we present the reformulated CGL equations and symmetrize the
conservative part. In Section 4, we propose entropy-stable schemes and prove
semi-discrete entropy stability for the CGL equations. Time stepping is
discussed in Section 5, followed by numerical results in Section 6. Final
conclusions are presented in Section 7.

##  2 CGL equations for anisotropic plasma flows

Following huang2019six ; meng2012classical , the CGL equations can be written
as follows:

$\displaystyle\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\textbf{u})=0,$
(1a) $\displaystyle\frac{\partial\rho\textbf{u}}{\partial
t}+\nabla\cdot\left[\rho\textbf{u}\textbf{u}+p_{\perp}\textbf{I}+(p_{\parallel}-p_{\perp})\boldsymbol{b}\boldsymbol{b}-\left(\boldsymbol{B}\boldsymbol{B}-\frac{|\boldsymbol{B}|^{2}}{2}\textbf{I}\right)\right]=0,$
(1b) $\displaystyle\frac{\partial p_{\parallel}}{\partial
t}+\nabla\cdot(p_{\parallel}\textbf{u})=-2p_{\parallel}\boldsymbol{b}\cdot\nabla\textbf{u}\cdot\boldsymbol{b},$
(1c) $\displaystyle\frac{\partial e}{\partial
t}+\nabla\cdot\left[\textbf{u}\left(e+p_{\perp}+\frac{|\boldsymbol{B}|^{2}}{2}\right)+\textbf{u}\cdot\left((p_{\parallel}-p_{\perp})\boldsymbol{b}\boldsymbol{b}-\boldsymbol{B}\boldsymbol{B}\right)\right]=0,$
(1d) $\displaystyle\frac{\partial\boldsymbol{B}}{\partial
t}+\nabla\times[-(\textbf{u}\times\boldsymbol{B})]=0.$ (1e)

Here, $\rho,\textbf{u}=(u_{x},u_{y},u_{z})^{\top}$, and
$\boldsymbol{B}=(B_{x},B_{y},B_{z})^{\top}$ are the density, velocity, and
magnetic field, respectively. The vector $\boldsymbol{b}$ is the unit vector
in the direction of the magnetic field, defined as
$\boldsymbol{b}=\boldsymbol{B}/|\boldsymbol{B}|$. One of the key shortcomings
of the CGL system is that a physics-free study of the CGL equations seems to
require magnetic field $\boldsymbol{B}$ to be non-zero. Otherwise, the
conventional wisdom has been that the magnetic field direction is not defined,
and hence the equations are not valid. However, it is useful to mention that a
recent paper by bhoriya2024  shows that a physics-based approach to the CGL
equations resolves this and many other issues associated with simulating the
CGL equations. However, we will assume that the magnetic field is non-zero.
Integration of the mathematics-based approach here with the physics-based
approach in bhoriya2024  will be the topic of a future paper. The pressure
tensor $\boldsymbol{P}$, rotated by the magnetic field is described using
scalars $p_{\parallel}$ and $p_{\perp}$ via relation,

$\displaystyle\boldsymbol{P}=p_{\parallel}\boldsymbol{b}\boldsymbol{b}+p_{\perp}(\textbf{I}-\boldsymbol{b}\boldsymbol{b})+\boldsymbol{\Pi}$
$\displaystyle=p_{\parallel}\begin{pmatrix}b_{x}b_{x}&b_{x}b_{y}&b_{x}b_{z}\\\
b_{y}b_{x}&b_{y}b_{y}&b_{y}b_{z}\\\
b_{z}b_{x}&b_{z}b_{y}&b_{z}b_{z}\end{pmatrix}+p_{\perp}\begin{pmatrix}1-b_{x}b_{x}&-b_{x}b_{y}&-b_{x}b_{z}\\\
-b_{y}b_{x}&1-b_{y}b_{y}&-b_{y}b_{z}\\\
-b_{z}b_{x}&-b_{z}b_{y}&1-b_{z}b_{z}\end{pmatrix}+\begin{pmatrix}\Pi_{xx}&\Pi_{xy}&\Pi_{xz}\\\
\Pi_{yx}&\Pi_{yy}&\Pi_{yz}\\\ \Pi_{zx}&\Pi_{zy}&\Pi_{zz}\end{pmatrix}$

where the scalars are given by
$p_{\parallel}=\boldsymbol{P}:\boldsymbol{b}\boldsymbol{b}$ and
$p_{\perp}=\boldsymbol{P}:\textbf{I}-\boldsymbol{b}\boldsymbol{b}$, denoting
parallel and perpendicular components of the pressure tensor, respectively.
The CGL model is obtained by neglecting the Finite Larmor Radius pressure
tensor, i.e., we assume $\boldsymbol{\Pi}=0$. The equation of state is given
by,

$\displaystyle
e=\frac{\rho|\textbf{u}|^{2}}{2}+\frac{|\boldsymbol{B}|^{2}}{2}+\frac{3p}{2}$
(2)

where, $p=\frac{2p_{\perp}+p_{\parallel}}{3}$ is the average scalar pressure.
Given the above notations, we define the vector of conservative variables as
$\mathbf{U}=(\rho,\rho\textbf{u},p_{\parallel},e,\boldsymbol{B})^{\top}$.
Equation (1a) represents the conservation of mass, followed by Equation (1b),
representing momentum conservation. The parallel component of the pressure
tensor is governed by (1c). This is the only equation that contains non-
conservative terms. The energy conservation is given by equation (1d),
followed by the induction equation (1e) for the magnetic field. In addition,
we also have divergence-free conditions on the magnetic field given by,

$\nabla\cdot\boldsymbol{B}=0$ (3)

However, this condition is an outcome of the induction equations (1e) and
divergence free initial condition.

The eigenvalues of the CGL system kato1966propagation ; meng2012classical  in
the $x$-direction are given by,

$\boldsymbol{\Lambda}=(u_{x},~{}u_{x},~{}0,~{}u_{x}\pm c_{a},~{}u_{x}\pm
c_{f},~{}u_{x}\pm c_{s})$ (4)

where,

$\displaystyle
c_{a}=\sqrt{{\frac{B_{x}^{2}}{\rho}-\frac{(p_{\parallel}-p_{\perp}){b}_{x}^{2}}{\rho}}},$
$\displaystyle
c_{f}=\frac{1}{\sqrt{2\rho}}\big{[}|\boldsymbol{B}|^{2}+2p_{\bot}+b_{x}^{2}(2p_{\parallel}-p_{\perp})+\\{(|\boldsymbol{B}|^{2}+2p_{\perp}+b_{x}^{2}(2p_{\parallel}-p_{\perp}))^{2}$
$\displaystyle~{}~{}~{}~{}~{}~{}+4(p_{\perp}^{2}b_{x}^{2}(1-b_{x}^{2})-3p_{\parallel}p_{\perp}b_{x}^{2}(2-b_{x}^{2})+3p_{\parallel}^{2}b_{x}^{4}-3B_{x}^{2}p_{\parallel})\\}^{\frac{1}{2}}\big{]}^{\frac{1}{2}},$
$\displaystyle
c_{s}=\frac{1}{\sqrt{2\rho}}\big{[}|\boldsymbol{B}|^{2}+2p_{\perp}+b_{x}^{2}(2p_{\parallel}-p_{\perp})-\\{(|\boldsymbol{B}|^{2}+2p_{\perp}+b_{x}^{2}(2p_{\parallel}-p_{\perp}))^{2}$
$\displaystyle~{}~{}~{}~{}~{}~{}+4(p_{\perp}^{2}b_{x}^{2}(1-b_{x}^{2})-3p_{\parallel}p_{\perp}b_{x}^{2}(2-b_{x}^{2})+3p_{\parallel}^{2}b_{x}^{4}-3B_{x}^{2}p_{\parallel})\\}^{\frac{1}{2}}\big{]}^{\frac{1}{2}}.$

It was shown in kato1966propagation  that eigenvalues of the system can become
imaginary. Hence, some conditions on the admissible domain need to be imposed.
Following kato1966propagation , we define,

$p_{m}=\frac{p_{\bot}^{2}}{6p_{\bot}+3|\boldsymbol{B}|^{2}},\text{ and
}p_{M}={|\boldsymbol{B}|^{2}}+p_{\bot}.$

Then, the admissible domain for the CGL system is given by,

$\Omega=\\{\textbf{U}\in\mathbb{R}^{9}|~{}\rho>0,~{}p_{\parallel}>0,p_{\perp}>0,~{}p_{m}\leq
p_{\parallel}\leq p_{M}\\}.$ (5)

In addition, it was observed in kato1966propagation , that the slow wave is
not necessarily slower than the Alfvén wave, hence the domain $\Omega$ has to
be divided as follows:

  1. 1.

$p_{m}\leq
p_{\parallel}\leq\frac{p_{M}}{4},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\text{if}~{}~{}~{}~{}~{}~{}~{}c_{s}\leq
c_{a}\leq c_{f},$

  2. 2.

$\frac{p_{M}}{4}\leq
p_{\parallel}\leq\frac{p_{M}}{4}+\frac{3p_{m}}{4},~{}~{}\text{if}~{}~{}~{}~{}~{}~{}c_{s}\leq
c_{a}<c_{f},$

  3. 3.

$\frac{p_{M}}{4}+\frac{3p_{m}}{4}\leq p_{\parallel}\leq
p_{M},~{}~{}\text{if}~{}~{}~{}~{}~{}~{}c_{a}\leq c_{s}<c_{f}.$

###  2.1 Entropy analysis

Several hyperbolic systems consider pressure tensor as one of the variables.
Following berthon2006numerical ; biswas2021entropy ; sen_entropy_2018 ;
berthon2015entropy ; yadav2023entropy , we define the entropy $\mathcal{E}$
and the entropy flux $\mathcal{Q}_{x}$ for the CGL system (1), as

$\displaystyle\mathcal{E}=-\rho s,\qquad\mathcal{Q}_{x}=-\rho u_{x}s,\qquad
s=\ln\left(\dfrac{\det\boldsymbol{P}}{\rho^{5}}\right)=\ln\left(\dfrac{p_{\parallel}p_{\perp}^{2}}{\rho^{5}}\right).$
(6)

###### Lemma 1

If U is a smooth solution of the CGL system (1a)-(1e), then it satisfies,

$\partial_{t}\mathcal{E}+\partial_{x}\mathcal{Q}_{x}=0.$ (7)

###### Proof

By the chain rule,

$\displaystyle\partial_{t}(p_{\parallel}p_{\perp}^{2})=p_{\perp}^{2}\partial_{t}p_{\parallel}+2p_{\parallel}p_{\perp}\partial_{t}p_{\perp}.$

Using (1a)-(1e) to derive the evolution of $p_{\perp}$, to get,

$\displaystyle\partial_{t}p_{\perp}=$
$\displaystyle-\partial_{x}(p_{\perp}u_{x})-p_{\perp}\partial_{x}u_{x}+p_{\perp}b_{x}(b_{x}{\partial_{x}u_{x}}+b_{y}{\partial_{x}u_{y}}+b_{z}{\partial_{x}u_{z}})-{(\textbf{u}\cdot\boldsymbol{B})}{\partial_{x}B_{x}}.$
(8)

Combing with (1c), we get,

$\displaystyle\partial_{t}(p_{\parallel}p_{\perp}^{2})$ $\displaystyle=$
$\displaystyle
p_{\perp}^{2}\partial_{t}p_{\parallel}+2p_{\parallel}p_{\perp}\partial_{t}p_{\perp}$
$\displaystyle=$ $\displaystyle-
u_{x}\partial_{x}(det\boldsymbol{P})-5(det(\boldsymbol{P}))\partial_{x}u_{x}-2p_{\parallel}p_{\perp}(\textbf{u}\cdot\boldsymbol{B})\partial_{x}B_{x},$

Using (1a), we have,

$\displaystyle\partial_{t}s+u_{x}\partial_{x}s+2\frac{(\textbf{u}\cdot\boldsymbol{B})}{p_{\perp}}\partial_{x}B_{x}=0.$

which immediately results in

$\partial_{t}\mathcal{E}+\partial_{x}\mathcal{Q}_{x}+\frac{2\rho(\textbf{u}\cdot\boldsymbol{B})}{p_{\perp}}\partial_{x}B_{x}=0.$
(9)

This combining with $\nabla\cdot\boldsymbol{B}=0$, gives (7).

###### Remark 1

In two dimensions, (9) is replaced with,

$\partial_{t}\mathcal{E}+\partial_{x}\mathcal{Q}_{x}+\partial_{y}\mathcal{Q}_{y}+\frac{2\rho(\textbf{u}\cdot\boldsymbol{B})}{p_{\perp}}\nabla\cdot\boldsymbol{B}=0,$
(10)

with $\mathcal{Q}_{y}=-\rho u_{y}s$. Furthermore, the entropy equality (7),
for the weak solutions results in,

$\partial_{t}\mathcal{E}+\partial_{x}\mathcal{Q}_{x}+\partial_{y}\mathcal{Q}_{y}\leq
0.$ (11)

We aim to design the entropy stable schemes so that (11) holds at the semi-
discrete level.

##  3 Reformulation of the CGL system

The set of equations (1) are not suitable for designing the entropy stable
schemes. To do so, motivated by yadav2023entropy , we will rewrite the
equations. We will consider the system (1) in two dimensions. We define
entropy variable $\textbf{V}=\frac{\partial\mathcal{E}}{\partial\textbf{U}}$,
which is given by,

V
$\displaystyle=\left(V_{1},~{}V_{2},~{}V_{3},~{}V_{4},~{}V_{5},~{}V_{6},~{}V_{7},~{}V_{8},~{}V_{9}\right)^{\top}$
$\displaystyle=\Bigg{(}5-s-\beta_{\perp}|\textbf{u}|^{2},~{}2\beta_{\perp}\textbf{u},~{}-\beta_{\parallel}+\beta_{\perp},~{}-2\beta_{\perp},~{}2\beta_{\perp}\boldsymbol{B}\Bigg{)}^{\top},$
(12)

where, $\beta_{\perp}=\frac{\rho}{p_{\perp}},$ and
$\beta_{\parallel}=\frac{\rho}{p_{\parallel}}$. Let us now recall the
following definition from yadav2023entropy :

###### Definition 1

A convex function $\mathcal{E}(\textbf{U})$ is said to be an entropy function
for the system

$\frac{\partial\textbf{U}}{\partial t}+\frac{\partial\textbf{f}_{x}}{\partial
x}+\frac{\partial\textbf{f}_{y}}{\partial
y}+\textbf{C}_{x}(\textbf{U})\frac{\partial\textbf{U}}{\partial
x}+\textbf{C}_{y}(\textbf{U})\frac{\partial\textbf{U}}{\partial y}=0.$

if there exist smooth functions $\mathcal{Q}_{x}(\textbf{U})$ and
$\mathcal{Q}_{y}(\textbf{U})$ such that

${\mathcal{Q}_{x}}^{\prime}(\textbf{U})=\textbf{V}{\textbf{f}_{x}}^{\prime}(\textbf{U}),\qquad{\mathcal{Q}_{y}}^{\prime}(\textbf{U})=\textbf{V}{\textbf{f}_{y}}^{\prime}(\textbf{U})$
(13)

and

$\textbf{V}^{\top}\textbf{C}_{x}(\textbf{U})=\textbf{V}^{\top}\textbf{C}_{y}(\textbf{U})=0.$
(14)

The functions ($\mathcal{E},\mathcal{Q}_{x},\mathcal{Q}_{y}$) form an entropy-
entropy flux pair.

The non-conservative terms in (1) do not satisfy (14). Hence, we rewrite the
equations (1) in the form,

$\displaystyle\frac{\partial\textbf{U}}{\partial
t}+\frac{\partial\textbf{f}_{x}}{\partial
x}+\frac{\partial\textbf{f}_{y}}{\partial
y}+\textbf{C}_{x}(\textbf{U})\frac{\partial\textbf{U}}{\partial
x}+\textbf{C}_{y}(\textbf{U})\frac{\partial\textbf{U}}{\partial y}=0.$ (15)

with flux functions given by,

$\displaystyle\textbf{f}_{x}=\begin{pmatrix}\rho u_{x}\\\ \rho
u_{x}^{2}+p_{\perp}-\left(B_{x}^{2}-\frac{|\boldsymbol{B}|^{2}}{2}\right)\\\
\rho u_{x}u_{y}-B_{x}B_{y}\\\ \rho u_{x}u_{z}-B_{x}B_{z}\\\
p_{\parallel}u_{x}\\\
u_{x}\left(e+p_{\perp}+\frac{|\boldsymbol{B}|^{2}}{2}\right)-{B_{x}}(\boldsymbol{B}\cdot\textbf{u})\\\
0\\\ u_{x}B_{y}-u_{y}B_{x}\\\
u_{x}B_{z}-u_{z}B_{x}\end{pmatrix},~{}\textbf{f}_{y}$
$\displaystyle=\begin{pmatrix}\rho u_{y}\\\ \rho u_{x}u_{y}-B_{x}B_{y}\\\ \rho
u_{y}^{2}+p_{\perp}-\left(B_{y}^{2}-\frac{|\boldsymbol{B}|^{2}}{2}\right)\\\
\rho u_{y}u_{z}-B_{y}B_{z}\\\ p_{\parallel}u_{y}\\\
u_{y}\left(e+p_{\perp}+\frac{|\boldsymbol{B}|^{2}}{2}\right)-B_{y}(\boldsymbol{B}\cdot\textbf{u})\\\
u_{y}B_{x}-u_{x}B_{y}\\\ 0\\\ u_{y}B_{z}-u_{z}B_{y}\end{pmatrix}.$

We note that the above fluxes are very similar to the MHD fluxes. The rest of
the terms are considered in the non-conservative form with matrices
$\textbf{C}_{x}$ and $\textbf{C}_{y}$. The non-conservative form with matrices
$\textbf{C}_{x}$ and $\textbf{C}_{y}$ given by, Matrix
$\textbf{C}_{x}(\textbf{U})$ is given below,

$\displaystyle\begin{pmatrix}0&0&0&0&0&0&0&0&0\\\
-\frac{b_{x}^{2}|\textbf{u}|^{2}}{2}&b_{x}^{2}u_{x}&b_{x}^{2}u_{y}&b_{x}^{2}u_{z}&\frac{3}{2}b_{x}^{2}&-b_{x}^{2}&b_{x}^{2}B_{x}+\frac{2\Gamma_{x}b_{x}}{|\boldsymbol{B}|}&b_{x}^{2}B_{y}-\frac{2\Delta
Pb_{yxx}}{|\boldsymbol{B}|}&b_{x}^{2}B_{z}-\frac{2\Delta
Pb_{zxx}}{|\boldsymbol{B}|}\\\
-\frac{b_{x}b_{y}|\textbf{u}|^{2}}{2}&b_{x}b_{y}u_{x}&b_{x}b_{y}u_{y}&b_{x}b_{y}u_{z}&\frac{3}{2}b_{x}b_{y}&-b_{x}b_{y}&b_{x}b_{y}B_{x}+\frac{\Gamma_{x}b_{y}-\Delta
Pb_{yxx}}{|\boldsymbol{B}|}&b_{x}b_{y}B_{y}+\frac{\Gamma_{y}b_{x}-\Delta
Pb_{xyy}}{|\boldsymbol{B}|}&b_{x}b_{y}B_{z}-\frac{2\Delta
Pb_{xyz}}{|\boldsymbol{B}|}\\\
-\frac{b_{x}b_{z}|\textbf{u}|^{2}}{2}&b_{x}b_{z}u_{x}&b_{x}b_{z}u_{y}&b_{x}b_{z}u_{z}&\frac{3}{2}b_{x}b_{z}&-b_{x}b_{z}&b_{x}b_{z}B_{x}+\frac{\Gamma_{x}b_{z}-\Delta
Pb_{zxx}}{|\boldsymbol{B}|}&b_{x}b_{z}B_{y}-\frac{2\Delta
Pb_{xyz}}{|\boldsymbol{B}|}&b_{x}b_{z}B_{z}+\frac{\Gamma_{z}b_{x}-\Delta
Pb_{xzz}}{|\boldsymbol{B}|}\\\
-\frac{2p_{\parallel}b_{x}}{\rho}(\boldsymbol{b}\cdot\textbf{u})&\frac{2p_{\parallel}b_{x}}{\rho}b_{x}&\frac{2p_{\parallel}b_{x}}{\rho}b_{y}&\frac{2p_{\parallel}b_{x}}{\rho}b_{z}&0&0&0&0&0\\\
\Upsilon_{1}^{x}&\Upsilon_{2}^{x}&\Upsilon_{3}^{x}&\Upsilon_{4}^{x}&\frac{3}{2}b_{x}(\boldsymbol{b}\cdot\textbf{u})&-b_{x}(\boldsymbol{b}\cdot\textbf{u})&b_{x}(\boldsymbol{b}\cdot\textbf{u})B_{x}+\Theta^{x}_{1}&b_{x}(\boldsymbol{b}\cdot\textbf{u})B_{y}+\Theta^{x}_{2}&b_{x}(\boldsymbol{b}\cdot\textbf{u})B_{z}+\Theta^{x}_{3}\\\
0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0\end{pmatrix}$

and matrix $\textbf{C}_{y}(\textbf{U})$ is given below,

$\displaystyle\begin{pmatrix}0&0&0&0&0&0&0&0&0\\\
-\frac{b_{x}b_{y}|\textbf{u}|^{2}}{2}&b_{x}b_{y}u_{x}&b_{x}b_{y}u_{y}&b_{x}b_{y}u_{z}&\frac{3}{2}b_{x}b_{y}&-b_{x}b_{y}&b_{x}b_{y}B_{x}+\frac{\Gamma_{x}b_{y}-\Delta
Pb_{yxx}}{|\boldsymbol{B}|}&b_{x}b_{y}B_{y}+\frac{\Gamma_{y}b_{x}-\Delta
Pb_{xyy}}{|\boldsymbol{B}|}&b_{x}b_{y}B_{z}-\frac{2\Delta
Pb_{xyz}}{|\boldsymbol{B}|}\\\
-\frac{b_{y}^{2}|\textbf{u}|^{2}}{2}&b_{y}^{2}u_{x}&b_{y}^{2}u_{y}&b_{y}^{2}u_{z}&\frac{3}{2}b_{y}^{2}&-b_{y}^{2}&b_{y}^{2}B_{x}-\frac{2\Delta
Pb_{xyy}}{|\boldsymbol{B}|}&b_{y}^{2}B_{y}+\frac{2\Gamma_{y}b_{y}}{|\boldsymbol{B}|}&b_{y}^{2}B_{z}-\frac{2\Delta
Pb_{zyy}}{|\boldsymbol{B}|}\\\
-\frac{b_{y}b_{z}|\textbf{u}|^{2}}{2}&b_{y}b_{z}u_{x}&b_{y}b_{z}u_{y}&b_{y}b_{z}u_{z}&\frac{3}{2}b_{y}b_{z}&-b_{y}b_{z}&b_{y}b_{z}B_{x}-\frac{2\Delta
Pb_{xyz}}{|\boldsymbol{B}|}&b_{y}b_{z}B_{y}+\frac{\Gamma_{y}b_{z}-\Delta
Pb_{zyy}}{|\boldsymbol{B}|}&b_{y}b_{z}B_{z}+\frac{\Gamma_{z}b_{y}-\Delta
Pb_{yzz}}{|\boldsymbol{B}|}\\\
-\frac{2p_{\parallel}b_{y}}{\rho}(\boldsymbol{b}\cdot\textbf{u})&\frac{2p_{\parallel}b_{y}}{\rho}b_{x}&\frac{2p_{\parallel}b_{y}}{\rho}b_{y}&\frac{2p_{\parallel}b_{y}}{\rho}b_{z}&0&0&0&0&0\\\
\Upsilon_{1}^{y}&\Upsilon_{2}^{y}&\Upsilon^{y}_{3}&\Upsilon_{4}^{y}&\frac{3}{2}b_{y}(\boldsymbol{b}\cdot\textbf{u})&-b_{y}(\boldsymbol{b}\cdot\textbf{u})&b_{y}(\boldsymbol{b}\cdot\textbf{u})B_{x}+\Theta_{1}^{y}&b_{y}(\boldsymbol{b}\cdot\textbf{u})B_{y}+\Theta_{2}^{y}&b_{y}(\boldsymbol{b}\cdot\textbf{u})B_{z}+\Theta_{3}^{y}\\\
0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0\end{pmatrix}$

where,
$\Delta P=(p_{\parallel}-p_{\perp}),~{}\Gamma_{i}=\Delta
P(1-b_{i}^{2}),~{}\forall i\in\left\\{x,y,z\right\\}$
$b_{lmn}=b_{l}b_{m}b_{n},~{}\forall~{}l,m,n\in\left\\{x,y,z\right\\}$
$\Theta^{x}_{1}=\\{(\boldsymbol{b}\cdot\textbf{u})+b_{x}u_{x}\\}\left(\frac{\Gamma_{x}}{|\boldsymbol{B}|}\right)-\Delta
P\frac{b^{2}_{x}b_{y}u_{y}}{|\boldsymbol{B}|}-\Delta
P\frac{b^{2}_{x}b_{z}u_{z}}{|\boldsymbol{B}|},$
$\Theta^{x}_{2}=\Gamma_{y}\left(\frac{b_{x}u_{y}}{|\boldsymbol{B}|}\right)-\Delta
P\frac{b_{x}b_{y}b_{z}u_{z}}{|\boldsymbol{B}|}-\\{(\boldsymbol{b}\cdot\textbf{u})+b_{x}u_{x}\\}\frac{\Delta
Pb_{x}b_{y}}{|\boldsymbol{B}|},$
$\Theta^{x}_{3}=\Gamma_{z}\left(\frac{b_{x}u_{z}}{|\boldsymbol{B}|}\right)-\Delta
P\frac{b_{x}b_{y}b_{z}u_{y}}{|\boldsymbol{B}|}-\\{(\boldsymbol{b}\cdot\textbf{u})+b_{x}u_{x}\\}\frac{\Delta
Pb_{x}b_{z}}{|\boldsymbol{B}|}$,
$\Theta_{1}^{y}=\Gamma_{x}\left(\frac{b_{y}u_{x}}{|\boldsymbol{B}|}\right)-\Delta
P\frac{b_{x}b_{y}b_{z}u_{z}}{|\boldsymbol{B}|}-\\{(\boldsymbol{b}\cdot\textbf{u})+b_{y}u_{y}\\}\frac{\Delta
Pb_{x}b_{y}}{|\boldsymbol{B}|},$
$\Theta_{2}^{y}=\\{(\boldsymbol{b}\cdot\textbf{u})+b_{y}u_{y}\\}\left(\frac{\Gamma_{y}}{|\boldsymbol{B}|}\right)-\Delta
P\frac{b^{2}_{y}b_{x}u_{x}}{|\boldsymbol{B}|}-\Delta
P\frac{b^{2}_{y}b_{z}u_{z}}{|\boldsymbol{B}|},$
$\Theta_{3}^{y}=\Gamma_{z}\left(\frac{b_{y}u_{z}}{|\boldsymbol{B}|}\right)-\Delta
P\frac{b_{x}b_{y}b_{z}u_{x}}{|\boldsymbol{B}|}-\\{(\boldsymbol{b}\cdot\textbf{u})+b_{y}u_{y}\\}\frac{\Delta
Pb_{y}b_{z}}{|\boldsymbol{B}|},$
$\Upsilon_{1}^{x}=-b_{x}(\boldsymbol{b}\cdot\textbf{u})\frac{|\textbf{u}|^{2}}{2}-\frac{\Delta
Pb_{x}(\boldsymbol{b}\cdot\textbf{u})}{\rho},~{}\Upsilon_{2}^{x}=b_{x}(\boldsymbol{b}\cdot\textbf{u})u_{x}+\frac{\Delta
Pb^{2}_{x}}{\rho},$
$\Upsilon_{3}^{x}=b_{x}(\boldsymbol{b}\cdot\textbf{u})u_{y}+\frac{\Delta
Pb_{x}b_{y}}{\rho},~{}\Upsilon_{4}^{x}=b_{x}(\boldsymbol{b}\cdot\textbf{u})u_{z}+\frac{\Delta
Pb_{x}b_{z}}{\rho}$
$\Upsilon_{1}^{y}=-b_{y}(\boldsymbol{b}\cdot\textbf{u})\frac{|\textbf{u}|^{2}}{2}-\frac{\Delta
Pb_{y}(\boldsymbol{b}\cdot\textbf{u})}{\rho},~{}\Upsilon_{2}^{y}=b_{y}(\boldsymbol{b}\cdot\textbf{u})u_{x}+\frac{\Delta
Pb_{x}b_{y}}{\rho},$
$\Upsilon_{3}^{y}=b_{y}(\boldsymbol{b}\cdot\textbf{u})u_{y}+\frac{\Delta
Pb_{y}^{2}}{\rho}$ and
$\Upsilon_{4}^{y}=b_{y}(\boldsymbol{b}\cdot\textbf{u})u_{z}+\frac{\Delta
Pb_{y}b_{z}}{\rho}$.

With this choice of the non-conservative matrices $\textbf{C}_{x}$ and
$\textbf{C}_{y}$, we satisfy the equation (14). The complete proof is given in
the Appendix (B).

###### Remark 2

The existence of an entropy pair does not guarantee the symmetrizability of
the system in the case of non-conservative hyperbolic systems. In particular,
the CGL system is not symmetrizable. The proof is given in Appendix (C)

###  3.1 Godunov’s Symmetrization of the conservative part

We will now analyze the conservative part of the equation (15) for entropy and
symmetrization, i.e. we consider,

$\displaystyle\frac{\partial\textbf{U}}{\partial
t}+\frac{\partial\textbf{f}_{x}}{\partial
x}+\frac{\partial\textbf{f}_{y}}{\partial y}=0,$ (16)

###### Definition 2

The conservation law (16) is said to be symmetrizable if there exists a change
of variable $\textbf{U}\rightarrow\textbf{V}$ which symmetrizes it, i.e., (16)
with a change of variable can be written as,

$\frac{\partial\textbf{U}}{\partial\textbf{V}}\frac{\partial\textbf{V}}{\partial{t}}+\frac{\partial\textbf{f}_{x}}{\partial\textbf{U}}\frac{\partial\textbf{U}}{\partial\textbf{V}}\frac{\partial\textbf{V}}{\partial{x}}+\frac{\partial\textbf{f}_{y}}{\partial\textbf{U}}\frac{\partial\textbf{U}}{\partial\textbf{V}}\frac{\partial\textbf{V}}{\partial{y}}=0.$

where $\frac{\partial\textbf{U}}{\partial\textbf{V}}$ is a symmetric positive
definite matrix. Furthermore,
$\frac{\partial\textbf{f}_{x}}{\partial\textbf{U}}\frac{\partial\textbf{U}}{\partial\textbf{V}}$
and
$\frac{\partial\textbf{f}_{y}}{\partial\textbf{U}}\frac{\partial\textbf{U}}{\partial\textbf{V}}$
are symmetric matrices.

The symmetrization of the hyperbolic conservation law and the existence of an
entropy pair are closely related. Let us recall the following Theorem from
godlewski2013numerical .

###### Theorem 3.1

A necessary and sufficient condition for the system to have a strictly convex
entropy $\mathcal{E}(\textbf{U})$ is that there exists a change of dependent
variables $\textbf{U}=\textbf{U}(\textbf{V})$ that symmetrizes system (16).

We note that for the system (16),

$\displaystyle\frac{\partial\textbf{f}_{\alpha}}{\partial\textbf{V}}\neq\bigg{(}\frac{\partial\textbf{f}_{\alpha}}{\partial\textbf{V}}\bigg{)}^{\top},~{}~{}\alpha=x,y.$

Furthermore, the entropy-pair $(\mathcal{E},\mathcal{Q}_{x},\mathcal{Q}_{y})$
does not satisfy (13). In fact, we have,

$\displaystyle\frac{\partial\mathcal{Q}_{\alpha}}{\partial\textbf{U}}=\textbf{V}^{\top}\cdot\frac{\partial\textbf{f}_{\alpha}}{\partial\textbf{U}}+\frac{2\rho}{p_{\perp}}(\textbf{u}\cdot\boldsymbol{B})B_{\alpha}^{\prime}(\textbf{U}).$

This is similar to the case of MHD equations chandrashekar2016entropy . Hence,
the conservative part of the CGL system is not symmetrizable. We follow
Godunov’s approach to attain the symmetrizability of the conservative part of
the CGL system with divergence constraints. We now follow
chandrashekar2016entropy , to symmetrize (16), which based on the earlier
works in godunov ; barth1999numerical ; barth2006role . We will consider the
equations in $x$-direction only. The discussion can be extended to include
$y$-directional flux, similarly.

Consider the modified conservative system (16) in $x$-direction,

$\frac{\partial\textbf{U}}{\partial t}+\frac{\partial\textbf{f}_{x}}{\partial
x}+\phi^{\prime}(\textbf{V})^{\top}\partial_{x}B_{x}=0,$ (17)

where $\phi(\textbf{V})$ is a homogeneous function of degree one, i.e.,

$\textbf{V}\phi^{\prime}(\textbf{V})=\phi(\textbf{V}).$ (18)

The above modification is consistent because of the divergence-free
constraint. If the system (17) is symmetrized by the transformation
$\textbf{U}\rightarrow\textbf{V}$, then it can be inferred that there are
twice differentiable functions $\mathcal{U}(\textbf{V})$ and
$\mathcal{F}_{x}(\textbf{V})$, where $\mathcal{U}(\textbf{V})$ is a strictly
convex function such that

$\textbf{U}=\mathcal{U}^{\prime}(\textbf{V})^{\top},\quad\textbf{f}_{x}=\mathcal{F}_{x}^{\prime}(\textbf{V})^{\top}-\phi^{\prime}(\textbf{V})^{\top}B_{x}.$
(19)

Combining (17) and (19), we get the symmetric form,

$\mathcal{U}^{\prime\prime}(\textbf{V})\frac{\partial\textbf{V}}{\partial
t}+(\mathcal{F}_{x}^{\prime\prime}(\textbf{V})-\phi^{\prime\prime}(\textbf{V})B_{x})\frac{\partial\textbf{V}}{\partial
x}=0.$

Here, $\mathcal{U}^{\prime\prime}(\textbf{V})$ is a symmetric positive
definite matrix and
$(\mathcal{F}_{x}^{\prime\prime}(\textbf{V})-\phi^{\prime\prime}(\textbf{V})B_{x})$
is a symmetric matrix. Let us also define,

$\mathcal{E}(\textbf{U})=\textbf{V}(\textbf{U})\cdot\textbf{U}-\mathcal{U}(\textbf{V}(\textbf{U})),\quad\mathcal{Q}_{x}(\textbf{U})=\textbf{V}(\textbf{U})\cdot\mathcal{F}_{x}^{\prime}(\textbf{V}(\textbf{U}))-\mathcal{F}_{x}(\textbf{V}(\textbf{U}))$
(20)

On differentiating the above equation and using (19) we get

$\mathcal{E}^{\prime}(\textbf{U})=\textbf{V}^{\prime}(\textbf{U})\cdot\textbf{U}+\textbf{V}(\textbf{U})-\mathcal{U}^{\prime}(\textbf{V})\textbf{V}^{\prime}(\textbf{U})=\textbf{V}^{\prime}(\textbf{U})\cdot\textbf{U}+\textbf{V}(\textbf{U})-\textbf{U}\cdot\textbf{V}^{\prime}(\textbf{U})=\textbf{V}$
(21)

i.e.

$\textbf{V}=\mathcal{E}^{\prime}(\textbf{U})^{\top}.$

Hence, the change of variable that symmetrizes the equations is entropy
variable $\textbf{V}(\textbf{U})$. Now differentiating $\mathcal{Q}_{x}$ from
(21),

$\displaystyle\mathcal{Q}_{x}^{\prime}(\textbf{U})$
$\displaystyle=\bigg{(}\textbf{V}(\textbf{U})\cdot\mathcal{F}_{x}^{\prime}(\textbf{V}(\textbf{U}))-\mathcal{F}_{x}(\textbf{V}(\textbf{U}))\bigg{)}^{\prime}$
$\displaystyle=\bigg{(}\textbf{V}(\textbf{U})\cdot(\textbf{f}_{x}(\textbf{V})^{\top}+\phi^{\prime}(\textbf{V})^{\top}B_{x})-\mathcal{F}_{x}(\textbf{V}(\textbf{U}))\bigg{)}^{\prime}$
$\displaystyle=\bigg{(}\textbf{V}(\textbf{U})\cdot\textbf{f}_{x}+\phi(\textbf{V})B_{x}-\mathcal{F}_{x}(\textbf{V}(\textbf{U}))\bigg{)}^{\prime}\quad\text{(from
eq. \eqref{eq:phi})}$
$\displaystyle=\bigg{(}\textbf{V}(\textbf{U})\cdot\textbf{f}_{x}^{\prime}(\textbf{U})+\phi(\textbf{V})B_{x}^{\prime}(\textbf{U})\bigg{)}=\bigg{(}\mathcal{E}^{\prime}(\textbf{U})\cdot\textbf{f}_{x}^{\prime}(\textbf{U})+\phi(\textbf{V})B_{x}^{\prime}(\textbf{U})\bigg{)}.$
(22)

Now, taking the dot product of V with (17), we get,

$\displaystyle\textbf{V}\cdot\bigg{(}\frac{\partial\textbf{U}}{\partial
t}+\frac{\partial\textbf{f}_{x}}{\partial
x}+\phi^{\prime}(\textbf{V})\nabla\cdot\boldsymbol{B}\bigg{)}$
$\displaystyle=\frac{\partial\mathcal{E}}{\partial
t}+\bigg{(}{\mathcal{E}}^{\prime}(\textbf{U})\cdot\frac{\partial\textbf{f}_{x}}{\partial\textbf{U}}+\phi(\textbf{V}){B^{\prime}_{x}(\textbf{U})}\bigg{)}\frac{\partial\textbf{U}}{\partial
x}$ $\displaystyle=\frac{\partial\mathcal{E}}{\partial
t}+\frac{\partial\mathcal{Q}_{x}}{\partial x}=0.\quad\text{ (using
\eqref{eq:cgl_entf'})}$

Now using,(22) with (19), we get,

$\phi(\textbf{V})=\frac{2\rho}{p_{\perp}}(\textbf{u}\cdot\boldsymbol{B})=2\beta_{\perp}(\textbf{u}\cdot\boldsymbol{B}).$
(23)

We also note that the function $\phi(\textbf{V})$ can be written using the
components of V, which results in,

$\phi(\textbf{V})=-\frac{V_{2}V_{7}+V_{3}V_{8}+V_{4}V_{9}}{V_{6}}.$

We note that this is a homogeneous function of degree one. Furthermore,

$\phi^{\prime}(\textbf{V})=[0,\boldsymbol{B},0,\textbf{u}\cdot\boldsymbol{B},\textbf{u}].$

Now, the entropy flux is given by,

$\mathcal{Q}_{x}=\textbf{V}\cdot\textbf{f}_{x}+\phi B_{x}-\mathcal{F}_{x}$
(24)

with

$\mathcal{U}=\textbf{V}\cdot\textbf{U}-\mathcal{E}=2\rho+\beta_{\perp}|\boldsymbol{B}|^{2},$

and

$\mathcal{F}_{x}=\textbf{V}\cdot\textbf{f}_{x}+\phi
B_{x}-\mathcal{Q}_{x}=2\rho
u_{x}+\beta_{\perp}u_{x}|\boldsymbol{B}|^{2}=\mathcal{U}u_{x}.$

Similarly, for the $y$-direction, we get,

$\mathcal{Q}_{y}=\textbf{V}\cdot\textbf{f}_{y}+\phi B_{y}-\mathcal{F}_{y},$
(25)

with

$\mathcal{F}_{y}=\textbf{V}\cdot\textbf{f}_{y}+\phi
B_{y}-\mathcal{Q}_{y}=2\rho
u_{y}+\beta_{\perp}u_{y}|\boldsymbol{B}|^{2}=\mathcal{U}u_{y}.$

We will now construct entropy stable numerical schemes for CGL equations by
adding $\phi^{\prime}(\textbf{V})^{\top}\nabla\cdot\boldsymbol{B}$, i.e. we
will consider,

$\displaystyle\frac{\partial\textbf{U}}{\partial
t}+\frac{\partial\textbf{f}_{x}}{\partial
x}+\frac{\partial\textbf{f}_{y}}{\partial
y}+\textbf{C}_{x}(\textbf{U})\frac{\partial\textbf{U}}{\partial
x}+\textbf{C}_{y}(\textbf{U})\frac{\partial\textbf{U}}{\partial
y}+\phi^{\prime}(\textbf{V})^{\top}(\nabla\cdot\boldsymbol{B})=0.$ (26)

##  4 Semi-discrete entropy stable numerical schemes

In this Section, we will develop semi-discrete numerical schemes for the
system (26), which satisfy a discrete version of the entropy inequality (11),
with modified entropy fluxes (24) and (25).

We consider a uniform mesh of size $\Delta x\times\Delta y$. The cell centers
are represented as $(i,j)$ with cell interfaces for the cell $(i,j)$ is given
by $(i-\frac{1}{2},j)$ and $(i+\frac{1}{2},j)$ in $x$-direction and
$(i,j-\frac{1}{2})$ and $(i,j+\frac{1}{2})$ in $y$-direction. Let us define
the arithmetic averages across the interfaces for the quantity $a$ as

$\bar{a}_{i+\frac{1}{2},j}=\frac{a_{i+1,j}+a_{i,j}}{2}\qquad\text{ and
}\qquad\bar{a}_{i,j+\frac{1}{2}}=\frac{a_{i,j+1}+a_{i,j}}{2}.$

Similarly, let us denote the jumps across the interfaces in $x$ and
$y$-directions as,

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

The semi-discrete finite difference schemes for (26) is given by,

$\displaystyle\frac{d}{dt}\textbf{U}_{i,j}$
$\displaystyle+\frac{\textbf{F}_{x,i+\frac{1}{2},j}-\textbf{F}_{x,i-\frac{1}{2},j}}{\Delta
x}+\frac{\textbf{F}_{y,i,j+\frac{1}{2}}-\textbf{F}_{y,i,j-\frac{1}{2}}}{\Delta
y}$
$\displaystyle+\phi^{\prime}(\textbf{V}_{i,j})^{\top}\bigg{(}\left(\frac{\partial
B_{x}}{\partial x}\right)_{i,j}+\left(\frac{\partial B_{y}}{\partial
y}\right)_{i,j}\bigg{)}$
$\displaystyle+\textbf{C}_{x}(\textbf{U}_{i,j})\left(\frac{\partial\textbf{U}}{\partial
x}\right)_{i,j}+\textbf{C}_{y}(\textbf{U}_{i,j})\left(\frac{\partial\textbf{U}}{\partial
y}\right)_{i,j}=0.$ (27)

Here $\textbf{F}_{x,i+\frac{1}{2},j}$ and $\textbf{F}_{y,i,j+\frac{1}{2}}$ are
the numerical fluxes in $x$ and $y$ directions, respectively. Furthermore,
$\left(\frac{\partial B_{x}}{\partial x}\right)_{i,j},\left(\frac{\partial
B_{y}}{\partial y}\right)_{i,j},\left(\frac{\partial\textbf{U}}{\partial
x}\right)_{i,j}$ and $\left(\frac{\partial\textbf{U}}{\partial
y}\right)_{i,j}$ are discretize using central difference schemes of suitable
orders. We will first construct an entropy conservative scheme for (26) by
designing entropy conservative fluxes.

###  4.1 Entropy conservative schemes

To design entropy conservative schemes of the form (4), we will combine two
key ideas. The first is to design the entropy conservative scheme for the
symmetrized conservative part. We achieve this by following
chandrashekar2016entropy , where a suitable entropy conservative flux is used
to design second-order numerical schemes. Here, we will also extend this
scheme to higher-order schemes as well. The second idea is using the fact that
the nonconservative terms in reformulated equations do not contribute to the
entropy evolution.

Following tadmor2003entropy ; chandrashekar2016entropy , we have the following
result:

###### Theorem 4.1

Assume that the numerical fluxes
$\tilde{\textbf{F}}_{x,i+\frac{1}{2},j},~{}\tilde{\textbf{F}}_{y,i,j+\frac{1}{2}}$
consistent with continuous fluxes in $x$ and $y$ direction satisfy,

$\displaystyle\quad[\\![\textbf{V}]\\!]_{i+\frac{1}{2},j}\cdot\tilde{\textbf{F}}_{x,i+\frac{1}{2},j}$
$\displaystyle=[\\![\mathcal{F}_{x}]\\!]_{i+\frac{1}{2},j}-[\\![\phi]\\!]_{i+\frac{1}{2},j}\bar{B}_{x,i+\frac{1}{2},j},$
(28)
$\displaystyle\quad[\\![\textbf{V}]\\!]_{i,j+\frac{1}{2}}\cdot\tilde{\textbf{F}}_{y,i,j+\frac{1}{2}}$
$\displaystyle=[\\![\mathcal{F}_{y}]\\!]_{i,j+\frac{1}{2}}-[\\![\phi]\\!]_{i,j+\frac{1}{2}}\bar{B}_{y,i,j+\frac{1}{2}}.$
(29)

Furthermore, the non-conservative terms $\left(\frac{\partial B_{x}}{\partial
x}\right)_{i,j},\left(\frac{\partial B_{y}}{\partial
y}\right)_{i,j},\left(\frac{\partial\textbf{U}}{\partial x}\right)_{i,j}$ and
$\left(\frac{\partial\textbf{U}}{\partial y}\right)_{i,j}$ are discretized
with second order central difference formula. Then the semi-discrete finite
difference scheme (4) is second-order accurate and entropy conservative, i.e.,
it satisfies

$\displaystyle\frac{d}{dt}\mathcal{E}\left(\textbf{U}_{i,j}\right)+\frac{\tilde{\mathcal{Q}}_{x,i+\frac{1}{2},j}-\tilde{\mathcal{Q}}_{x,i-\frac{1}{2},j}}{\Delta
x}+\frac{\tilde{\mathcal{Q}}_{y,i,j+\frac{1}{2}}-\tilde{\mathcal{Q}}_{y,i,j-\frac{1}{2}}}{\Delta
y}=0$ (30)

where the numerical entropy fluxes consistent with (24) and (25) are given by

$\displaystyle\tilde{\mathcal{Q}}_{x,i+\frac{1}{2},j}=\tilde{\mathcal{Q}}_{x}(\textbf{U}_{i,j},\textbf{U}_{i+1,j})=\bar{\textbf{V}}_{i+\frac{1}{2},j}\cdot\tilde{\textbf{F}}_{x,i+\frac{1}{2},j}+\bar{\phi}_{i+\frac{1}{2},j}\bar{B}_{x,i+\frac{1}{2},j}-\bar{\mathcal{F}}_{x,i+\frac{1}{2},j},$
(31)
$\displaystyle\tilde{\mathcal{Q}}_{y,i,j+\frac{1}{2}}=\tilde{\mathcal{Q}}_{y}(\textbf{U}_{i,j},\textbf{U}_{i,j+1})=\bar{\textbf{V}}_{i,j+\frac{1}{2}}\cdot\tilde{\textbf{F}}_{y,i,j+\frac{1}{2}}+\bar{\phi}_{i,j+\frac{1}{2}}\bar{B}_{y,i,j+\frac{1}{2}}-\bar{\mathcal{F}}_{y,i,j+\frac{1}{2}}.$
(32)

###### Proof

Following chandrashekar2016entropy , we note that,

$\displaystyle\tilde{\mathcal{Q}}_{x,i+\frac{1}{2},j}-\tilde{\mathcal{Q}}_{x,i-\frac{1}{2},j}$
$\displaystyle=\bar{\textbf{V}}_{i+\frac{1}{2},j}\cdot\tilde{\textbf{F}}_{x,i+\frac{1}{2},j}-\bar{\textbf{V}}_{i-\frac{1}{2},j}\cdot\tilde{\textbf{F}}_{x,i-\frac{1}{2},j}$
$\displaystyle+\bar{\phi}_{i+\frac{1}{2},j}\bar{B}_{x,i+\frac{1}{2},j}-\bar{\phi}_{i-\frac{1}{2},j}\bar{B}_{x,i-\frac{1}{2},j}-\left(\bar{\mathcal{F}}_{x,i+\frac{1}{2},j}-\bar{\mathcal{F}}_{x,i-\frac{1}{2},j}\right)$
(33)

Now,

$\displaystyle\bar{\mathcal{F}}_{x,i+\frac{1}{2},j}-\bar{\mathcal{F}}_{x,i-\frac{1}{2},j}$
$\displaystyle=\frac{1}{2}\left(\mathcal{F}_{x,i+1,j}+\mathcal{F}_{x,i,j}\right)+\frac{1}{2}\left(\mathcal{F}_{x,i,j}+\mathcal{F}_{x,i-1,j}\right).$

Using (28), we get,

$\displaystyle\bar{\mathcal{F}}_{x,i+\frac{1}{2},j}-\bar{\mathcal{F}}_{x,i-\frac{1}{2},j}$
$\displaystyle=\frac{1}{2}\left(\left(\textbf{V}_{i+1,j}-\textbf{V}_{i,j}\right)\cdot\tilde{\textbf{F}}_{x,i+\frac{1}{2},j}+\left(\phi_{i+1,j}-\phi_{i,j}\right)\bar{B}_{x,i+\frac{1}{2},j}\right)$
$\displaystyle+\frac{1}{2}\left(\left(\textbf{V}_{i,j}-\textbf{V}_{i-1,j}\right)\cdot\tilde{\textbf{F}}_{x,i-\frac{1}{2},j}+\left(\phi_{i,j}-\phi_{i-1,j}\right)\bar{B}_{x,1-\frac{1}{2},j}\right).$
(34)

Substituting (4.1) in (4.1) results in,

$\displaystyle\tilde{\mathcal{Q}}_{x,i+\frac{1}{2},j}-\tilde{\mathcal{Q}}_{x,i-\frac{1}{2},j}$
$\displaystyle=\textbf{V}_{i,j}\cdot\left(\tilde{\textbf{F}}_{x,i+\frac{1}{2},j}-\tilde{\textbf{F}}_{x,i-\frac{1}{2},j}\right)+\phi_{i,j}\left(\bar{B}_{x,i+\frac{1}{2},j}-\bar{B}_{x,i-\frac{1}{2},j}\right)$
(35)
$\displaystyle=\textbf{V}_{i,j}\cdot\left(\tilde{\textbf{F}}_{x,i+\frac{1}{2},j}-\tilde{\textbf{F}}_{x,i-\frac{1}{2},j}\right)+\Delta
x\textbf{V}_{i,j}\cdot\phi^{\prime}(\textbf{U}_{i,j})\left(\left(\frac{\partial
B_{x}}{\partial x}\right)_{i,j}\right)$ $\displaystyle+\Delta
x\textbf{V}_{i,j}^{\top}\textbf{C}_{x}(\textbf{U}_{i,j})\left(\frac{\partial\textbf{U}}{\partial
x}\right)_{i,j}.$

Similarly, using (29), we can get,

$\displaystyle\tilde{\mathcal{Q}}_{y,i,j+\frac{1}{2}}-\tilde{\mathcal{Q}}_{y,i,j-\frac{1}{2}}$
$\displaystyle=\textbf{V}_{i,j}\cdot\left(\tilde{\textbf{F}}_{y,i,j+\frac{1}{2}}-\tilde{\textbf{F}}_{y,i,j-\frac{1}{2}}\right)+\Delta
y\textbf{V}_{i,j}\cdot\phi^{\prime}(\textbf{U}_{i,j})\left(\left(\frac{\partial
B_{y}}{\partial y}\right)_{i,j}\right)$ $\displaystyle+\Delta
y\textbf{V}_{i,j}^{\top}\textbf{C}_{y}(\textbf{U}_{i,j})\left(\frac{\partial\textbf{U}}{\partial
y}\right)_{i,j}.$

Combining, we get the discrete entropy equality (30).

####  4.1.1 Entropy conservative fluxes

We now aim to design second order entropy conservative fluxes
$\tilde{\textbf{F}}_{x,i+\frac{1}{2},j}$ and
$\tilde{\textbf{F}}_{y,i,j+\frac{1}{2},}$ which satisfy (28) and (29),
respectively. To simplify the presentation, we will drop the mesh index and
consider left and right states denoted by $\textbf{U}_{l}$ and
$\textbf{U}_{r}$. Using this, let us denote,

$[\\![a]\\!]=a_{r}-a_{l}\qquad\text{ and
}\qquad\bar{a}=\frac{a_{r}+a_{l}}{2}.$

Then we aim to construct
$\tilde{\textbf{F}}_{x}(\textbf{U}_{l},\textbf{U}_{r})$, such that,

$[\\![\textbf{V}]\\!]\cdot\tilde{\textbf{F}}_{x}=[\\![\mathcal{F}_{x}]\\!]-\bar{B}_{x}[\\![\phi]\\!].$

We will first proceed in $x$-direction. The case for $y$-direction is similar.
We also note that the above equation is a scalar. So, we need a vector
$\tilde{\textbf{F}}_{x}$, satisfying the above scalar equations. Hence, in
general, we can have many entropy-conservative numerical fluxes. For Euler
equations of compressible equations, several such fluxes are presented in
tadmor2003entropy ; chandrashekar2013kinetic ; ismail2009affordable .

$[\\![ab]\\!]=\bar{a}[\\![b]\\!]+\bar{b}[\\![a]\\!].$

Following roe2006affordable , let us introduce the logarithmic average.

$a^{\ln}=\frac{[\\![a]\\!]}{[\\![\ln{a}]\\!]}.$

When $a_{l}\approx a_{r},$ we use the Taylor expansion-based formula from
ismail2009affordable . Let us first simplify the $[\\![s]\\!]$, to get,

$\displaystyle[\\![s]\\!]$
$\displaystyle=[\\![\left(\ln\left(p_{\parallel}p_{\perp}^{2}/\rho^{3}\right)-2\ln(\rho)\right)]\\!]$
$\displaystyle=[\\![\left(\ln\left(p_{\parallel}/\rho\right)+\ln\left(p_{\perp}^{2}/\rho^{2}\right)-2\ln(\rho)\right)]\\!]$
$\displaystyle=-[\\![\ln(\beta_{\parallel})]\\!]-2[\\![\ln(\beta_{\perp})]\\!]-2[\\![\ln(\rho)]\\!]$
$\displaystyle=-\dfrac{[\\![\beta_{\parallel}]\\!]}{{\beta}^{\ln}_{\parallel}}-2\dfrac{[\\![\beta_{\perp}]\\!]}{{\beta}^{\ln}_{\perp}}-2\dfrac{[\\![\rho]\\!]}{{\rho}^{\ln}}$

Following the approach in chandrashekar2013kinetic , we expand the jumps of
entropy variables V in terms of the jumps in
$\rho,u_{x},u_{y},u_{z},\beta_{\parallel},\beta_{\perp},B_{x},B_{y},B_{z}$ as

$\displaystyle[\\![\textbf{V}]\\!]=\left[\begin{array}[]{c}\frac{2}{\rho^{\ln}}[\\![\rho]\\!]+\frac{1}{{\beta}^{\ln}_{\parallel}}[\\![\beta_{\parallel}]\\!]-2\bar{\beta}_{\perp}\left(\bar{u}_{x}[\\![u_{x}]\\!]+\bar{u}_{y}[\\![u_{y}]\\!]+\bar{u}_{z}[\\![u_{z}]\\!]\right)+\left[\frac{2}{{\beta}^{\ln}_{\perp}}-\overline{|u|^{2}}\right][\\![\beta_{\perp}]\\!]\\\
2\bar{\beta}_{\perp}[\\![u_{x}]\\!]+2\bar{u}_{x}[\\![\beta_{\perp}]\\!]\\\
2\bar{\beta}_{\perp}[\\![u_{y}]\\!]+2\bar{u}_{y}[\\![\beta_{\perp}]\\!]\\\
2\bar{\beta}_{\perp}[\\![u_{z}]\\!]+2\bar{u}_{z}[\\![\beta_{\perp}]\\!]\\\
-[\\![\beta_{\parallel}]\\!]+[\\![\beta_{\perp}]\\!]\\\
-2[\\![\beta_{\perp}]\\!]\\\
2\bar{\beta}_{\perp}[\\![B_{x}]\\!]+2\bar{B}_{x}[\\![\beta_{\perp}]\\!]\\\
2\bar{\beta}_{\perp}[\\![B_{y}]\\!]+2\bar{B}_{y}[\\![\beta_{\perp}]\\!]\\\
2\bar{\beta}_{\perp}[\\![B_{z}]\\!]+2\bar{B}_{z}[\\![\beta_{\perp}]\\!]\end{array}\right].$

Let us consider the numerical flux,

$\tilde{\textbf{F}}_{x}=\left[\tilde{F}_{x}^{(1)},\tilde{F}_{x}^{(2)},\tilde{F}_{x}^{(3)},\tilde{F}_{x}^{(4)},\tilde{F}_{x}^{(5)},\tilde{F}_{x}^{(6)},\tilde{F}_{x}^{(7)},\tilde{F}_{x}^{(8)},\tilde{F}_{x}^{(9)}\right]^{\top}.$

Then $[\\![\textbf{V}]\\!]\cdot\tilde{\textbf{F}}_{x}$ can be written as,

$\displaystyle[\\![\textbf{V}]\\!]\cdot\tilde{\textbf{F}}_{x}=$
$\displaystyle\frac{2\tilde{F}_{x}^{(1)}}{{\rho}^{\ln}}[\\![\rho]\\!]+\left(-2\bar{u}_{x}\bar{\beta}_{\perp}\tilde{F}_{x}^{(1)}+2\bar{\beta}_{\perp}\tilde{F}_{x}^{(2)}\right)[\\![u_{x}]\\!]$
$\displaystyle+\left(-2\bar{u}_{y}\bar{\beta}_{\perp}\tilde{F}_{x}^{(1)}+2\bar{\beta}_{\perp}\tilde{F}_{x}^{(3)}\right)[\\![u_{y}]\\!]$
$\displaystyle+\left(-2\bar{u}_{z}\bar{\beta}_{\perp}\tilde{F}_{x}^{(1)}+2\bar{\beta}_{\perp}\tilde{F}_{x}^{(4)}\right)[\\![u_{z}]\\!]+\left(\frac{\tilde{F}_{x}^{(1)}}{{\beta}^{\ln}_{\parallel}}-\tilde{F}_{x}^{(5)}\right)[\\![\beta_{\parallel}]\\!]$
$\displaystyle+\Bigg{[}\left(\frac{2}{{\beta}^{\ln}_{\perp}}-\overline{\left|u\right|^{2}}\right)\tilde{F}_{x}^{(1)}+2\bar{u}_{x}\tilde{F}_{x}^{(2)}+2\bar{u}_{y}\tilde{F}_{x}^{(3)}+2\bar{u}_{z}\tilde{F}_{x}^{(4)}+\tilde{F}_{x}^{(5)}-2\tilde{F}_{x}^{(6)}$
$\displaystyle+2\bar{B}_{x}\tilde{F}_{x}^{(7)}+2\bar{B}_{y}\tilde{F}_{x}^{(8)}+2\bar{B}_{z}\tilde{F}_{x}^{(9)}\Bigg{]}[\\![\beta_{\perp}]\\!]$
$\displaystyle+2\bar{\beta}_{\perp}\tilde{F}_{x}^{(7)}[\\![B_{x}]\\!]+2\bar{\beta}_{\perp}\tilde{F}_{x}^{(8)}[\\![B_{y}]\\!]+2\bar{\beta}_{\perp}\tilde{F}_{x}^{(9)}[\\![B_{z}]\\!]$
(36)

Similarly, we simplify the $[\\![\mathcal{F}_{x}]\\!]$ and $[\\![\phi]\\!]$,
to get,

$\displaystyle\\!]=$ $\displaystyle[\\![\left(2\rho
u_{x}\right)]\\!]+[\\![\left(\beta_{\perp}u_{x}|\boldsymbol{B}|^{2}\right)]\\!]$
$\displaystyle=$ $\displaystyle
2\bar{u}_{x}[\\![\rho]\\!]+\left(2\bar{\rho}+\bar{\beta}_{\perp}\overline{|\boldsymbol{B}|^{2}}\right)[\\![u_{x}]\\!]+\bar{u}_{x}\overline{\left|\boldsymbol{B}\right|^{2}}[\\![\beta_{\perp}]\\!]+2\overline{\beta_{\perp}u_{x}}\bar{B}_{x}[\\![B_{x}]\\!]$
$\displaystyle+2\overline{\beta_{\perp}u_{x}}\bar{B}_{y}[\\![B_{y}]\\!]+2\overline{\beta_{\perp}u_{x}}\bar{B}_{z}[\\![B_{z}]\\!],$

and

$\displaystyle\\!]=$ $\displaystyle
2[\\![(\beta_{\perp}(\boldsymbol{u}\cdot\boldsymbol{B}))]\\!]$
$\displaystyle=$ $\displaystyle
2\bar{\beta}_{\perp}\bar{B}_{x}[\\![u_{x}]\\!]+2\bar{\beta}_{\perp}\bar{B}_{y}[\\![u_{y}]\\!]+2\bar{\beta}_{\perp}\bar{B}_{z}[\\![u_{z}]\\!]+2\left(\bar{u}_{x}\bar{B}_{x}+\bar{u}_{y}\bar{B}_{y}+\bar{u}_{z}\bar{B}_{z}\right)[\\![\beta_{\perp}]\\!]$
$\displaystyle+2\overline{\beta u_{x}}[\\![B_{x}]\\!]+2\overline{\beta
u_{y}}[\\![B_{y}]\\!]+2\overline{\beta u_{z}}[\\![B_{z}]\\!].$

Combining, we get,

$\displaystyle[\\![\mathcal{F}_{x}]\\!]-\bar{B}_{x}[\\![\phi]\\!]$
$\displaystyle=2\bar{u}_{x}[\\![\rho]\\!]+\left(2\bar{\rho}+\bar{\beta}_{\perp}\overline{|\boldsymbol{B}|^{2}}-2\bar{\beta}_{\perp}\bar{B}_{x}\bar{B}_{x}\right)[\\![u_{x}]\\!]$
$\displaystyle-2\bar{\beta}_{\perp}\bar{B}_{x}\bar{B}_{y}[\\![u_{y}]\\!]-2\bar{\beta}_{\perp}\bar{B}_{x}\bar{B}_{z}[\\![u_{z}]\\!]$
$\displaystyle+\left[\bar{u}_{x}\overline{|\boldsymbol{B}|^{2}}-2\left(\bar{u}_{x}\bar{B}_{x}+\bar{u}_{y}\bar{B}_{y}+\bar{u}_{z}\bar{B}_{z}\right)\bar{B}_{x}\right][\\![\beta_{\perp}]\\!]$
$\displaystyle+2\left(\overline{\beta_{\perp}u_{x}}\bar{B}_{y}-\overline{\beta_{\perp}u_{y}}\bar{B}_{x}\right)[\\![B_{y}]\\!]$
$\displaystyle+2\left(\overline{\beta_{\perp}u_{x}}\bar{B}_{z}-\overline{\beta_{\perp}u_{z}}\bar{B}_{z}\right)[\\![B_{z}]\\!]$
(37)

Equating the coefficients of
$[\\![\rho]\\!],[\\![u_{x}]\\!],[\\![u_{y}]\\!],[\\![u_{z}]\\!],[\\![\beta_{\parallel}]\\!],[\\![\beta_{\perp}]\\!],[\\![B_{x}]\\!],[\\![B_{y}]\\!]$
and $[\\![B_{z}]\\!]$ in (4.1.1) and (4.1.1), we get,

$\displaystyle\tilde{F}_{x}^{(1)}=$
$\displaystyle{\rho}^{\ln}\bar{u}_{x},~{}\tilde{F}_{x}^{(2)}=\frac{\bar{\rho}}{\bar{\beta}_{\perp}}+\bar{u}_{x}\tilde{F}_{x}^{(1)}+\frac{1}{2}\overline{\left|\boldsymbol{B}\right|^{2}}-\bar{B}_{x}\bar{B}_{x}$
$\displaystyle\tilde{F}_{x}^{(3)}=$
$\displaystyle\bar{u}_{y}\tilde{F}_{x}^{(1)}-\bar{B}_{x}\bar{B}_{y},~{}\tilde{F}_{x}^{(4)}=\bar{u}_{z}\tilde{F}_{x}^{(1)}-\bar{B}_{x}\bar{B}_{z},~{}\tilde{F}_{x}^{(5)}=\frac{\tilde{F}_{x}^{(1)}}{{\beta}^{\ln}_{\parallel}},$
$\displaystyle\tilde{F}_{x}^{(7)}=$ $\displaystyle
0,~{}\tilde{F}_{x}^{(8)}=\frac{1}{\bar{\beta}_{\perp}}\left(\overline{\beta_{\perp}u_{x}}\bar{B}_{y}-\overline{\beta_{\perp}u_{y}}\bar{B}_{x}\right),~{}\tilde{F}_{x}^{(9)}=\frac{1}{\bar{\beta}_{\perp}}\left(\overline{\beta_{\perp}u_{x}}\bar{B}_{z}-\overline{\beta
u_{z}}\bar{B_{x}}\right),$ $\displaystyle\tilde{F}_{x}^{(6)}=$
$\displaystyle\frac{1}{2}\left[\frac{2}{{\beta}^{\ln}_{\perp}}-\overline{|u|^{2}}\right]\tilde{F}_{x}^{(1)}+\bar{u}_{x}\tilde{F}_{x}^{(2)}+\bar{u}_{y}\tilde{F}_{x}^{(3)}+\bar{u}_{z}\tilde{F}_{x}^{(4)}+\frac{\tilde{F}_{x}^{(5)}}{2}$
$\displaystyle+\bar{B}_{x}\tilde{F}_{x}^{(7)}+\bar{B}_{y}\tilde{F}_{x}^{(8)}+\bar{B}_{z}\tilde{F}_{x}^{(9)}-\frac{1}{2}\bar{u}_{x}\overline{\left|\boldsymbol{B}\right|^{2}}+\left(\bar{u}_{x}\bar{B}_{x}+\bar{u}_{y}\bar{B}_{y}+\bar{u}_{z}\bar{B}_{z}\right)\bar{B}_{x},$

Similarly, we obtain the entropy conservative numerical flux in $y$-direction,
$\tilde{\textbf{F}}_{y}=\left[\tilde{F}_{y}^{(1)},\tilde{F}_{y}^{(2)},\tilde{F}_{y}^{(3)},\tilde{F}_{y}^{(4)},\right.$
$\left.\tilde{F}_{y}^{(5)},\tilde{F}_{y}^{(6)},\tilde{F}_{y}^{(7)},\tilde{F}_{y}^{(8)},\tilde{F}_{y}^{(9)}\right]^{\top}$,
where

$\displaystyle\tilde{F}_{y}^{(1)}=$
$\displaystyle{\rho}^{\ln}\bar{u}_{y},~{}\tilde{F}_{y}^{(2)}=\bar{u}_{x}\tilde{F}_{y}^{(1)}-\bar{B}_{x}\bar{B}_{y},$
$\displaystyle\tilde{F}_{y}^{(3)}=$
$\displaystyle\frac{\bar{\rho}}{\bar{\beta}_{\perp}}+\bar{u}_{y}\tilde{F}_{y}^{(1)}+\frac{1}{2}\overline{\left|\boldsymbol{B}\right|^{2}}-\bar{B}_{y}\bar{B}_{y},~{}\tilde{F}_{y}^{(4)}=\bar{u}_{z}\tilde{F}_{y}^{(1)}-\bar{B}_{y}\bar{B}_{z},~{}\tilde{F}_{y}^{(5)}=\frac{\tilde{F}_{y}^{(1)}}{{\beta}^{\ln}_{\parallel}},$
$\displaystyle\tilde{F}_{y}^{(7)}=$
$\displaystyle\frac{1}{\bar{\beta}_{\perp}}\left(\overline{\beta_{\perp}u_{y}}\bar{B}_{x}-\overline{\beta_{\perp}u_{x}}\bar{B}_{y}\right),~{}\tilde{F}_{y}^{(8)}=0,~{}\tilde{F}_{y}^{(9)}=\frac{1}{\bar{\beta}_{\perp}}\left(\overline{\beta_{\perp}u_{y}}\bar{B}_{z}-\overline{\beta_{\perp}u_{z}}\bar{B_{y}}\right),$
$\displaystyle\tilde{F}_{y}^{(6)}=$
$\displaystyle\frac{1}{2}\left[\frac{2}{{\beta}^{\ln}_{\perp}}-\overline{|u|^{2}}\right]\tilde{F}_{y}^{(1)}+\bar{u}_{x}\tilde{F}_{y}^{(2)}+\bar{u}_{y}\tilde{F}_{y}^{(3)}+\bar{u}_{z}\tilde{F}_{y}^{(4)}+\frac{\tilde{F}_{y}^{(5)}}{2}$
$\displaystyle+\bar{B}_{x}\tilde{F}_{y}^{(7)}+\bar{B}_{y}\tilde{F}_{y}^{(8)}+\bar{B}_{z}\tilde{F}_{y}^{(9)}-\frac{1}{2}\bar{u}_{y}\overline{\left|\boldsymbol{B}\right|^{2}}+\left(\bar{u}_{x}\bar{B}_{x}+\bar{u}_{y}\bar{B}_{y}+\bar{u}_{z}\bar{B}_{z}\right)\bar{B}_{y},$

It is easy to check that the above numerical fluxes are consistent. So, the
numerical scheme (4) is now well defined with
$\tilde{\textbf{F}}_{x,i+\frac{1}{2},j}=\tilde{\textbf{F}}_{x}(\textbf{U}_{i,j},\textbf{U}_{i+1,j})$
and
$\tilde{\textbf{F}}_{y,i,j+\frac{1}{2}}=\tilde{\textbf{F}}_{y}(\textbf{U}_{i,j},\textbf{U}_{i,j+1}).$

####  4.1.2 Higher order entropy conservative schemes

In leFloch2002 , authors have constructed $2p$-th ($p$ is any positive
integer) order entropy conservative numerical fluxes by combining second-order
numerical flux. We will use the same formulations and construct a fourth-order
numerical flux. Using the approach we present here, we can further extend to
design $2p$-th order entropy conservative schemes.

Following leFloch2002 , let us define the $4$-th order numerical fluxes,

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

and

$\displaystyle\tilde{\textbf{F}}_{y,i,j+\frac{1}{2}}^{4}$
$\displaystyle=\frac{4}{3}\tilde{\textbf{F}}_{y}(\textbf{U}_{i,j},\textbf{U}_{i,j+1})-\frac{1}{6}\bigg{(}\tilde{\textbf{F}}_{y}(\textbf{U}_{i,j-1},\textbf{U}_{i,j+1})+\tilde{\textbf{F}}_{y}(\textbf{U}_{i,j},\textbf{U}_{i,j+2})\bigg{)}.$
(39)

Using these numerical fluxes, let us consider the following fourth-order
scheme:

$\displaystyle\frac{d}{dt}\textbf{U}_{i,j}$
$\displaystyle+\frac{\tilde{\textbf{F}}^{4}_{x,i+\frac{1}{2},j}-\tilde{\textbf{F}}^{4}_{x,i-\frac{1}{2},j}}{\Delta
x}+\frac{\tilde{\textbf{F}}^{4}_{y,i,j+\frac{1}{2}}-\tilde{\textbf{F}}^{4}_{y,i,j-\frac{1}{2}}}{\Delta
y}$
$\displaystyle+\phi^{\prime}(\textbf{V}_{i,j})^{\top}\bigg{(}\left(\frac{\partial
B_{x}}{\partial x}\right)_{i,j}+\left(\frac{\partial B_{y}}{\partial
y}\right)_{i,j}\bigg{)}$
$\displaystyle+\textbf{C}_{x}(\textbf{U}_{i,j})\left(\frac{\partial\textbf{U}}{\partial
x}\right)_{i,j}+\textbf{C}_{y}(\textbf{U}_{i,j})\left(\frac{\partial\textbf{U}}{\partial
y}\right)_{i,j}=0.$ (40)

where, the nonconservative terms $\left(\frac{\partial B_{x}}{\partial
x}\right)_{i,j},\left(\frac{\partial B_{y}}{\partial
y}\right)_{i,j},\left(\frac{\partial\textbf{U}}{\partial x}\right)_{i,j}$ and
$\left(\frac{\partial\textbf{U}}{\partial y}\right)_{i,j}$ are discretized
using fourth order central difference formula. Then we have the following
results:

###### Theorem 4.2

The numerical scheme (4.1.2) is entropy conservative i.e. it satisfies,

$\displaystyle\frac{d}{dt}\mathcal{E}\left(\textbf{U}_{i,j}\right)+\frac{\tilde{\mathcal{Q}}^{4}_{x,i+\frac{1}{2},j}-\tilde{\mathcal{Q}}^{4}_{x,i-\frac{1}{2},j}}{\Delta
x}+\frac{\tilde{\mathcal{Q}}^{4}_{y,i,j+\frac{1}{2},}-\tilde{\mathcal{Q}}^{4}_{y,i,j-\frac{1}{2}}}{\Delta
y}=0$ (41)

where the consistent numerical entropy fluxes are given by

$\displaystyle\tilde{\mathcal{Q}}^{4}_{x,i+\frac{1}{2},j}=\frac{4}{3}\tilde{\mathcal{Q}}_{x}(\textbf{U}_{i,j},\textbf{U}_{i+1,j})-\frac{1}{6}\bigg{(}\tilde{\mathcal{Q}}_{x}(\textbf{U}_{i-1,j},\textbf{U}_{i+1,j})+\tilde{\mathcal{Q}}_{x}(\textbf{U}_{i,j},\textbf{U}_{i+2,j})\bigg{)},$
$\displaystyle\tilde{\mathcal{Q}}^{4}_{y,i,j+\frac{1}{2}}=\frac{4}{3}\tilde{\mathcal{Q}}_{y}(\textbf{U}_{i,j},\textbf{U}_{i,j+1})-\frac{1}{6}\bigg{(}\tilde{\mathcal{Q}}_{y}(\textbf{U}_{i,j-1},\textbf{U}_{i,j+1})+\tilde{\mathcal{Q}}_{y}(\textbf{U}_{i,j},\textbf{U}_{i,j+2})\bigg{)}.$

###### Proof

We will first consider the $x$-directional case. Note that fourth-order
central discretization of $\frac{\partial B_{x}}{\partial x}$ can be written
as,

$\left(\frac{\partial B_{x}}{\partial
x}\right)_{i,j}=\frac{\bar{B}^{4}_{x,i+\frac{1}{2},j}-\bar{B}^{4}_{x,i,j+\frac{1}{2}}}{\Delta
x}$

where

$\displaystyle\bar{B}_{x,i+\frac{1}{2},j}^{4}$
$\displaystyle=\frac{4}{3}\bar{B}_{x,i+\frac{1}{2},j}-\frac{1}{6}\left(\frac{B_{x,i-1,j}+B_{x,i+1,j}}{2}+\frac{B_{x,i,j}+B_{x,i+2,j}}{2}\right).$

Now multiply the right hand side of (4.1.2) with $\textbf{V}_{i,j}$, we get,

$\displaystyle\textbf{V}_{i,j}\left(\tilde{\textbf{F}}_{x,i+\frac{1}{2},j}^{4}-\tilde{\textbf{F}}_{x,i-\frac{1}{2},j}^{4}\right)$
$\displaystyle+\textbf{V}_{i,j}\phi^{\prime}(\textbf{V}_{i,j})\left(\bar{B}_{x,i+\frac{1}{2},j}^{4}-\bar{B}_{x,i-\frac{1}{2},j}^{4}\right)+$
$\displaystyle\textbf{V}_{i,j}\textbf{C}_{x}(\textbf{U}_{i,j})\left(\frac{\partial\textbf{U}}{\partial
x}\right)_{i,j}$
$\displaystyle=\textbf{V}_{i,j}\left(\tilde{\textbf{F}}_{x,i+\frac{1}{2},j}^{4}-\tilde{\textbf{F}}_{x,i-\frac{1}{2},j}^{4}\right)$
$\displaystyle+\phi(\textbf{V}_{i,j})\left(\bar{B}_{x,i+\frac{1}{2},j}^{4}-\bar{B}_{x,i-\frac{1}{2},j}^{4}\right)$
$\displaystyle=\textbf{V}_{i,j}\frac{4}{3}\left(\tilde{\textbf{F}}_{x}(\textbf{U}_{i,j},\textbf{U}_{i+1,j})-\tilde{\textbf{F}}_{x}(\textbf{U}_{i-1,j},\textbf{U}_{i,j})\right)$
$\displaystyle+\phi(\textbf{V}_{i,j})\frac{4}{3}\left(\bar{B}_{x,i+\frac{1}{2},j}-\bar{B}_{x,i-\frac{1}{2},j}\right)$
$\displaystyle-\textbf{V}_{i,j}\frac{1}{6}\left(\tilde{\textbf{F}}_{x}(\textbf{U}_{i,j},\textbf{U}_{i+2,j})-\tilde{\textbf{F}}_{x}(\textbf{U}_{i-2,j},\textbf{U}_{i,j})\right)$
$\displaystyle-\phi(\textbf{V}_{i,j})\frac{1}{6}\left(\frac{B_{x,i,j}+B_{x,i+2,j}}{2}-\frac{B_{x,i-2,j}+B_{x,i,j}}{2}\right)$
$\displaystyle=\frac{4}{3}\left(\tilde{\mathcal{Q}}_{x}(\textbf{U}_{i,j},\textbf{U}_{i+1,j})-\tilde{\mathcal{Q}}_{x}(\textbf{U}_{i-1,j},\textbf{U}_{i,j})\right)~{}\text{
using \eqref{eq:entf_diff} }$
$\displaystyle-\frac{1}{6}\left(\tilde{\mathcal{Q}}_{x}(\textbf{U}_{i,j},\textbf{U}_{i+2,j})-\tilde{\mathcal{Q}}_{x}(\textbf{U}_{i-2,j},\textbf{U}_{i,j})\right)$
$\displaystyle=\left(\tilde{\mathcal{Q}}^{4}_{x,i+\frac{1}{2},j}-\tilde{\mathcal{Q}}^{4}_{x,i-\frac{1}{2},j}\right).$

Similarly, for $y$-direction,

$\displaystyle\left(\tilde{\mathcal{Q}}^{4}_{y,i,j+\frac{1}{2}}-\tilde{\mathcal{Q}}^{4}_{y,i,j-\frac{1}{2}}\right)$
$\displaystyle=\textbf{V}_{i,j}\left(\tilde{\textbf{F}}_{y,i,j+\frac{1}{2}}^{4}-\tilde{\textbf{F}}_{y,i,j-\frac{1}{2}}^{4}\right)$
$\displaystyle+\textbf{V}_{i,j}\phi^{\prime}(\textbf{V}_{i,j})\left(\frac{\partial
B_{y}}{\partial
y}\right)_{i,j}+\textbf{V}_{i,j}\textbf{C}_{y}(\textbf{U}_{i,j})\left(\frac{\partial\textbf{U}}{\partial
y}\right)_{i,j}.$

Combining $x$ and $y$-direction entropy flux difference, we now get (41).

###  4.2 Entropy stable schemes

Entropy conservative schemes produce oscillations when the solutions contain
shock. To ensure a stable scheme, we need to have entropy stability of the
numerical scheme. We will now modify the numerical flux to ensure the entropy
stability of the scheme. Following tadmor2003entropy , we introduce the
modified numerical flux as follows:

$\displaystyle{\hat{\textbf{F}}}_{x,i+\frac{1}{2},j}=\tilde{\textbf{F}}_{x,i+\frac{1}{2},j}-\frac{1}{2}\textbf{D}_{x,i+\frac{1}{2},j}[\\![\textbf{V}]\\!]_{i+\frac{1}{2},j},~{}{\hat{\textbf{F}}}_{y,i,j+\frac{1}{2}}=\tilde{\textbf{F}}_{y,i,j+\frac{1}{2}}-\frac{1}{2}\textbf{D}_{y,i,j+\frac{1}{2}}[\\![\textbf{V}]\\!]_{i,j+\frac{1}{2}},$
(42)

where $\textbf{D}_{x,i+\frac{1}{2},j}$ and $\textbf{D}_{y,i,j+\frac{1}{2}}$
are symmetric positive definite matrices. Combining Theorem 4.1, with Theorem
3.3 in chandrashekar2016entropy , we have the following result:

###### Theorem 4.3

The numerical scheme (4) with the modified numerical fluxes (42) is entropy
stable, i.e., the solutions satisfies,

$\frac{d}{dt}\mathcal{E}(\textbf{U}_{i,j})+\frac{1}{\Delta
x}\left(\hat{\mathcal{Q}}_{x,i+\frac{1}{2},j}-\hat{\mathcal{Q}}_{x,i-\frac{1}{2},j}\right)+\frac{1}{\Delta
y}\left(\hat{\mathcal{Q}}_{y,i,j+\frac{1}{2}}-\hat{\mathcal{Q}}_{y,i,j-\frac{1}{2}}\right)\leq
0,$

with consistent numerical entropy flux functions,

$\displaystyle\hat{\mathcal{Q}}_{x,i+\frac{1}{2},j}=\tilde{\mathcal{Q}}_{x,i+\frac{1}{2},j}-\frac{1}{2}\bar{\textbf{V}}_{i+\frac{1}{2},j}^{\top}\textbf{D}_{x,i+\frac{1}{2},j}[\\![\textbf{V}]\\!]_{i+\frac{1}{2},j}$

and

$\displaystyle\qquad\hat{\mathcal{Q}}_{y,i,j+\frac{1}{2}}=\tilde{\mathcal{Q}}_{y,i,j+\frac{1}{2}}-\frac{1}{2}\bar{\textbf{V}}_{i,j+\frac{1}{2}}^{\top}\textbf{D}_{y,i,j+\frac{1}{2}}[\\![\textbf{V}]\\!]_{i,j+\frac{1}{2}}.$

To define the diffusion part, we consider the symmetrized conservative part
(17), in quasi-linear form, i.e.

$\frac{\partial\textbf{U}}{\partial
t}+\mathcal{A}_{d}\frac{\partial\textbf{U}}{\partial
x}=0,\qquad\mathcal{A}_{d}=\frac{\partial\textbf{f}_{d}}{\partial\textbf{U}}+\phi^{\prime}(\textbf{V})B^{\prime}_{d}(\textbf{U}).$
(43)

with $d\in\\{x,y\\}$. As the above system is symmetrizable, we follow
barth1999numerical ; chandrashekar2016entropy  to derive entropy scaled
eigenvectors in Appendix A. We denote the matrix of entropy scaled right
eigenvectors as $\tilde{R}_{d}$. Furthermore, let us take
${\Lambda_{d}},\,d\in\\{x,y\\},$ are $9\times 9$ diagonal matrices of the form

${\Lambda_{d}}=\text{diag}\left\\{\left(\max_{1\leq k\leq
9}|\lambda^{k}_{d}|\right)\mathbf{I}_{9\times 9}\right\\},\ \
\,d\in\\{x,y\\}.$

where $\\{\lambda^{k}_{d}:1\leq k\leq 9\\}$ is the set of eigenvalues of the
jacobian $\mathcal{A}_{d}$. Then we define diffusion matrices,

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

With this choice, the numerical scheme with modified fluxes is stable in
entropy. However, the scheme contains the jump term $[\\![\textbf{V}]\\!]$ in
both directions, which results in only first-order accuracy, even with second
or higher-order entropy conservative numerical flux. One simple way to
overcome this is to reconstruct the jumps
$[\\![\textbf{V}]\\!]_{i+\frac{1}{2},j}$ and
$[\\![\textbf{V}]\\!]_{i,j+\frac{1}{2}}$ using higher order reconstruction
process. This indeed leads to higher-order schemes, but it is not possible to
prove the entropy stability of the scheme. Instead, we follow
fjordholm2012arbitrarily  to design higher-order entropy stable schemes.

####  4.2.1 Higher order entropy stable schemes

Following fjordholm2012arbitrarily , we define scaled entropy variables

$\textbf{W}^{\pm}_{x,i,j}={R}^{{\top}}_{x,i\pm\frac{1}{2},j}\textbf{V}_{i,j}.$

The key idea is to reconstruct $\textbf{W}^{\pm}_{x,i,j}$ instead on
$\textbf{V}^{\pm}_{x,i,j}$. Let $\hat{\textbf{W}}^{\pm}_{x,i,j}$, denotes the
sign preserving higher order reconstruction of $\textbf{W}^{\pm}_{x,i,j}$ in
$x$-direction. Then, we define reconstructed $\textbf{V}_{x,i,j}^{\pm}$ as,

$\hat{\textbf{V}}^{\pm}_{x,i,j}=\left\\{\tilde{R}^{{\top}}_{x,i\pm\frac{1}{2},j}\right\\}^{(-1)}\hat{\textbf{W}}_{x,i,j}^{\pm}.$

Hence, now the reconstructed $[\\![\textbf{V}]\\!]_{i+\frac{1}{2},j}$ of
$k$-th order is defined as,

$[\\![\hat{\textbf{V}}]\\!]^{k}_{x,i,j}=\hat{\textbf{V}}^{-}_{x,i+1,j}-\hat{\textbf{V}}^{+}_{x,i,j}.$

Using this, we now define higher-order entropy stable numerical flux in
$x$-direction as,

${\hat{\textbf{F}}}^{k}_{x,i+\frac{1}{2},j}\,=\,\tilde{\textbf{F}}^{2p}_{x,i+\frac{1}{2},j}\,-\,\frac{1}{2}\,\textbf{D}_{x,i+\frac{1}{2},j}[\\![\boldsymbol{\hat{V}}]\\!]_{x,i+\frac{1}{2},j}^{k}.$
(45)

We choose $p=2$ for first and second-order schemes and $p=4$ for third and
fourth-order schemes. We can follow a similar process in $y$-direction. To
ensure the entropy stability condition of the above scheme, the reconstruction
of the scaled entropy variable needs to satisfy sign preserving property for
each component, i.e., we need

$\textrm{sign}\left({\hat{\textbf{W}}}^{-}_{x,i+1,j,l}-{\hat{\textbf{W}}}^{+}_{x,i,j,l}\right)=\textrm{sign}\left(\textbf{W}_{i+1,j,l}-\textbf{W}_{i,j,l}\right)\quad\text{for}\quad
1\leq l\leq 9$ (46)

where $l$ denote the $l$-th component of vector W. This can be ensured by
using MinMod reconstruction for the second order scheme ($k=2$). For third
($k=3$) and fourth ($k=4$) order schemes, we need to use ENO-based
reconstruction fjordholm2013eno  of the third and fourth order. Similarly, in
$y$-direction, we define higher numerical flux as:

${\hat{\textbf{F}}}^{k}_{y,i+\frac{1}{2},j}\,=\,\tilde{\textbf{F}}^{2p}_{y,i+\frac{1}{2},j}\,-\,\frac{1}{2}\,\textbf{D}_{y,i+\frac{1}{2},j}[\\![\boldsymbol{\hat{V}}]\\!]_{y,i+\frac{1}{2},j}^{k}.$
(47)

where reconstruction is performed in $y$-direction.

###  4.3 Semi-discrete entropy stability

We will now prove the entropy stability of the proposed numerical schemes.
Consider the following numerical discretization:

$\displaystyle\frac{d}{dt}\textbf{U}_{i,j}$
$\displaystyle+\frac{\hat{\textbf{F}}^{k}_{x,i+\frac{1}{2},j}-\hat{\textbf{F}}^{k}_{x,i-\frac{1}{2},j}}{\Delta
x}+\frac{\hat{\textbf{F}}^{k}_{y,i,j+\frac{1}{2}}-\hat{\textbf{F}}^{k}_{y,i,j-\frac{1}{2}}}{\Delta
y}$
$\displaystyle+\phi^{\prime}(\textbf{V}_{i,j})^{\top}\bigg{(}\left(\frac{\partial
B_{x}}{\partial x}\right)_{i,j}+\left(\frac{\partial B_{y}}{\partial
y}\right)_{i,j}\bigg{)}$
$\displaystyle+\textbf{C}_{x}(\textbf{U}_{i,j})\left(\frac{\partial\textbf{U}}{\partial
x}\right)_{i,j}+\textbf{C}_{y}(\textbf{U}_{i,j})\left(\frac{\partial\textbf{U}}{\partial
y}\right)_{i,j}=0,$ (48)

where we use the higher order entropy stable numerical flux (45) and (47).
Furthermore, the derivatives $\left(\frac{\partial B_{x}}{\partial
x}\right)_{i,j},~{}\left(\frac{\partial B_{y}}{\partial
y}\right)_{i,j},~{}\left(\frac{\partial\textbf{U}}{\partial x}\right)_{i,j}$
and $\left(\frac{\partial\textbf{U}}{\partial y}\right)_{i,j}$ in the non-
conservative term are discretized using second-order central difference for
$k=1$ and $k=2$ and fourth order central difference for $k=3$ and $k=4$. Then,
we have the following result:

###### Theorem 4.4

The semi-discrete scheme (4.3) is $k$ order accurate and entropy stable, i.e.
it satisfies,

$\frac{d}{dt}\mathcal{E}(\textbf{U}_{i,j})+\frac{1}{\Delta
x}\left(\hat{\mathcal{Q}}^{k}_{x,i+\frac{1}{2},j}-\hat{\mathcal{Q}}^{k}_{x,i-\frac{1}{2},j}\right)+\frac{1}{\Delta
y}\left(\hat{\mathcal{Q}}^{k}_{y,i,j+\frac{1}{2}}-\hat{\mathcal{Q}}^{k}_{y,i,j-\frac{1}{2}}\right)\leq
0,$ (49)

where $\hat{\mathcal{Q}}^{k}_{x}$ and $\hat{\mathcal{Q}}^{k}_{y}$ are given
by,

$\displaystyle\hat{\mathcal{Q}}^{k}_{x,i+\frac{1}{2},j}=\tilde{\mathcal{Q}}^{2p}_{x,i+\frac{1}{2},j}-\frac{1}{2}\bar{\textbf{V}}_{i+\frac{1}{2},j}^{\top}\textbf{D}_{x,i+\frac{1}{2},j}[\\![\hat{\textbf{V}}]\\!]^{k}_{x,i+\frac{1}{2},j}$

and

$\displaystyle\qquad\hat{\mathcal{Q}}^{k}_{y,i,j+\frac{1}{2}}=\tilde{\mathcal{Q}}^{2p}_{y,i,j+\frac{1}{2}}-\frac{1}{2}\bar{\textbf{V}}_{i,j+\frac{1}{2}}^{\top}\textbf{D}_{y,i,j+\frac{1}{2}}[\\![\hat{\textbf{V}}]\\!]^{k}_{y,i,j+\frac{1}{2}}.$

###### Proof

Accuracy follows from fjordholm2012arbitrarily . We also note that from
Theorems (4.1) and (4.2), we get,

$\displaystyle\frac{d}{dt}\mathcal{E}(\textbf{U}_{i,j})$
$\displaystyle+\frac{1}{\Delta
x}\left(\hat{\mathcal{Q}}^{k}_{x,i+\frac{1}{2},j}-\hat{\mathcal{Q}}^{k}_{x,i-\frac{1}{2},j}\right)+\frac{1}{\Delta
y}\left(\hat{\mathcal{Q}}^{k}_{y,i,j+\frac{1}{2}}-\hat{\mathcal{Q}}^{k}_{y,i,j-\frac{1}{2}}\right)$
$\displaystyle=-\frac{1}{2\Delta
x}\left[[\\![\textbf{V}]\\!]_{i+\frac{1}{2},j}\textbf{D}_{x,i+\frac{1}{2},j}[\\![\hat{\textbf{V}}]\\!]^{k}_{x,i+\frac{1}{2},j}+[\\![\textbf{V}]\\!]_{i-\frac{1}{2},j}\textbf{D}_{x,i-\frac{1}{2},j}[\\![\hat{\textbf{V}}]\\!]^{k}_{x,i-\frac{1}{2},j}\right]$
$\displaystyle\quad\;-\frac{1}{2\Delta
y}\left[[\\![\textbf{V}]\\!]_{i,j+\frac{1}{2}}\textbf{D}_{y,i,j+\frac{1}{2}}[\\![\hat{\textbf{V}}]\\!]^{k}_{y,i,j+\frac{1}{2}}+[\\![\textbf{V}]\\!]_{i,j-\frac{1}{2}}\textbf{D}_{y,i,j-\frac{1}{2}}[\\![\hat{\textbf{V}}]\\!]^{k}_{y,i,j-\frac{1}{2}}\right]$
$\displaystyle=-\frac{1}{2\Delta
x}\left[[\\![\textbf{W}]\\!]_{i+\frac{1}{2},j}{\Lambda}_{x,i+\frac{1}{2},j}[\\![\hat{\textbf{W}}]\\!]^{k}_{x,i+\frac{1}{2},j}+[\\![\textbf{W}]\\!]_{i-\frac{1}{2},j}{\Lambda}_{x,i-\frac{1}{2},j}[\\![\hat{\textbf{W}}]\\!]^{k}_{x,i-\frac{1}{2},j}\right]$
$\displaystyle\quad\;-\frac{1}{2\Delta
y}\left[[\\![\textbf{W}]\\!]_{i,j+\frac{1}{2}}\Lambda_{y,i,j+\frac{1}{2}}[\\![\hat{\textbf{W}}]\\!]^{k}_{y,i,j+\frac{1}{2}}+[\\![\textbf{W}]\\!]_{i,j-\frac{1}{2}}\Lambda_{y,i,j-\frac{1}{2}}[\\![\hat{\textbf{W}}]\\!]^{k}_{y,i,j-\frac{1}{2}}\right]$
$\displaystyle\leq 0\qquad\qquad\text{ Using sign preserving property
\eqref{eq:sign_preser}}.$

To compute the solutions close to isotropic (MHD) limits, in some test cases,
we consider an additional source term in the equations (26) by adding,

$\textbf{S}=\Biggl{\\{}0,~{}0,~{}0,~{}0,~{}\frac{p_{\perp}-p_{\parallel}}{\tau},~{}0,~{}0,~{}0,~{}0\Biggr{\\}}^{\top},$
(50)

with $\tau=10^{-5}$. This effects only $p_{\parallel}$ equation (1c). In that
case, (4.3) results in,

$\displaystyle\frac{d}{dt}\textbf{U}_{i,j}$
$\displaystyle+\frac{\hat{\textbf{F}}^{k}_{x,i+\frac{1}{2},j}-\hat{\textbf{F}}^{k}_{x,i-\frac{1}{2},j}}{\Delta
x}+\frac{\hat{\textbf{F}}^{k}_{y,i,j+\frac{1}{2}}-\hat{\textbf{F}}^{k}_{y,i,j-\frac{1}{2}}}{\Delta
y}$
$\displaystyle+\phi^{\prime}(\textbf{V}_{i,j})^{\top}\bigg{(}\left(\frac{\partial
B_{x}}{\partial x}\right)_{i,j}+\left(\frac{\partial B_{y}}{\partial
y}\right)_{i,j}\bigg{)}$
$\displaystyle+\textbf{C}_{x}(\textbf{U}_{i,j})\left(\frac{\partial\textbf{U}}{\partial
x}\right)_{i,j}+\textbf{C}_{y}(\textbf{U}_{i,j})\left(\frac{\partial\textbf{U}}{\partial
y}\right)_{i,j}=\textbf{S}(\textbf{U}_{i,j}).$ (51)

In that case (49), results in,

$\displaystyle\frac{d}{dt}\mathcal{E}(\textbf{U}_{i,j})$
$\displaystyle+\frac{1}{\Delta
x}\left(\hat{\mathcal{Q}}^{k}_{x,i+\frac{1}{2},j}-\hat{\mathcal{Q}}^{k}_{x,i-\frac{1}{2},j}\right)+\frac{1}{\Delta
y}\left(\hat{\mathcal{Q}}^{k}_{y,i,j+\frac{1}{2}}-\hat{\mathcal{Q}}^{k}_{y,i,j-\frac{1}{2}}\right)$
$\displaystyle\leq\textbf{V}_{i,j}^{\top}\cdot\textbf{S}(\textbf{U}_{i,j})=-\frac{\rho_{i,j}}{p_{\parallel,i,j}p_{\perp,i,j}\tau}(p_{\parallel,i,j}-p_{\perp,i,j})^{2}\leq
0.$ (52)

##  5 Fully discrete schemes

The semi-discrete scheme (4.3), can be written as

$\frac{d}{dt}\textbf{U}_{i,j}(t)=\mathcal{L}_{i,j}(\textbf{U}(t))+\textbf{S}(\textbf{U}_{i,j}(t))$
(53)

where,

$\displaystyle\mathcal{L}_{i,j}(\textbf{U}(t))$
$\displaystyle=-\frac{\hat{\textbf{F}}^{k}_{x,i+\frac{1}{2},j}-\hat{\textbf{F}}^{k}_{x,i-\frac{1}{2},j}}{\Delta
x}-\frac{\hat{\textbf{F}}^{k}_{y,i,j+\frac{1}{2}}-\hat{\textbf{F}}^{k}_{y,i,j-\frac{1}{2}}}{\Delta
y}$
$\displaystyle-\phi^{\prime}(\textbf{V}_{i,j})^{\top}\bigg{(}\left(\frac{\partial
B_{x}}{\partial x}\right)_{i,j}+\left(\frac{\partial B_{y}}{\partial
y}\right)_{i,j}\bigg{)}$
$\displaystyle-\textbf{C}_{x}(\textbf{U}_{i,j})\left(\frac{\partial\textbf{U}}{\partial
x}\right)_{i,j}-\textbf{C}_{y}(\textbf{U}_{i,j})\left(\frac{\partial\textbf{U}}{\partial
y}\right)_{i,j}.$

To integrate the ODE (53), we use explicit and Implicit-Explicit (IMEX)
schemes (see gottlieb2001strong ; pareschi2005implicit ; bhoriya2023high ).

###  5.1 Explicit schemes

Following gottlieb2001strong , we use explicit strong stability preserving
Runge Kutta (SSP-RK) methods for the time discretization. The second and
third-order accurate SSP-RK schemes are as follows:

  1. 1.

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

  2. 2.

For $k$ in $\\{1,\dots,m+1\\}$, compute

$\displaystyle\textbf{U}_{i,j}^{(k)}\ =\
\sum_{l=0}^{k-1}\gamma_{kl}\textbf{U}_{i,j}^{(l)}+\delta_{kl}\Delta
t\Big{(}\mathcal{L}_{i,j}(\textbf{U}^{(l)})+\textbf{S}(\textbf{U}_{i,j}^{(l)})\Big{)},$

where $\gamma_{kl}$ and $\delta_{kl}$ are given in the Table (1).

  3. 3.

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

Table 1: Coefficients for Explicit SSP-Runge-Kutta time stepping Order |  $\gamma_{il}$ |  |  |  $\delta_{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

The fourth order SSP-RK scheme gottlieb2001strong  is given by,

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

where
$\mathcal{M}_{i,j}(\textbf{U}^{n})=\mathcal{L}_{i,j}(\textbf{U}^{n})+\textbf{S}(\textbf{U}^{n}_{i,j})$.

We will use explicit schemes to compute CGL solutions only, i.e. S in (50) is
ignored. We are now denoting explicit schemes without source term S, as
follows:

  1. 1.

O2ES-EXP: We consider second-order time stepping with MinMod reconstruction
for scaled entropy variables. We use the second-order central difference for
non-conservative terms. The source term is ignored.

  2. 2.

O3ES-EXP: We consider third-order time stepping with third-order ENO
reconstruction for scaled entropy variables. We use the fourth-order central
difference for non-conservative terms. The source term is ignored.

  3. 3.

O4ES-EXP: We consider fourth-order time stepping with fourth-order ENO
reconstruction for scaled entropy variables. We use the fourth-order central
difference for non-conservative terms. The source term is ignored.

###  5.2 ARK IMEX schemes

The source term S contains the parameter $\tau$, which makes the source term
stiff. To avoid the small time step restriction required for stability of the
explicit scheme, treat the source term implicitly. In particular, we use the
Additive Runge-Kutta (ARK) IMEX scheme. For the second-order discretization,
we use the L-Stable second-order accurate ARK-IMEX scheme from
pareschi2005implicit ; bhoriya2023high . The method has two implicit stages
and is given by

$\displaystyle\mathbf{U}_{i,j}^{(1)}$
$\displaystyle=\mathbf{U}^{n}_{i,j}+\Delta t\
\big{(}\beta\textbf{S}(\mathbf{U}^{(1)}_{i,j}\big{)},$
$\displaystyle\mathbf{U}^{(2)}_{i,j}$
$\displaystyle=\mathbf{U}^{n}_{i,j}+\Delta t\
\big{(}\mathcal{L}_{i,j}(\mathbf{U}^{(1)})+(1-2\beta)\textbf{S}(\mathbf{U}^{(1)}_{i,j})+\beta\textbf{S}(\mathbf{U}^{(2)}_{i,j})\big{)},$
$\displaystyle\mathbf{U}_{i,j}^{(n+1)}$
$\displaystyle=\mathbf{U}^{n}_{i,j}+\Delta t\
\bigg{(}\dfrac{1}{2}\mathcal{L}_{i,j}(\mathbf{U}^{(1)})+\dfrac{1}{2}\mathcal{L}_{i,j}(\mathbf{U}^{(2)})+\dfrac{1}{2}\textbf{S}(\mathbf{U}^{(1)}_{i,j})+\dfrac{1}{2}\textbf{S}(\mathbf{U}^{(2)}_{i,j})\bigg{)}.$

where $\beta=1-\dfrac{1}{\sqrt{2}}$.
For the third and fourth-order ARK scheme, we follow kennedy2003additive . In
the ARK-IMEX schemes, we need to solve a system of nonlinear equations to
solve the Implicit source part. We use the Newton solver, where the
convergence and the next direction search are based on the back-tracing line
search framework. A detailed procedure is given in dennis1996numerical .

We will use IMEX schemes to compute CGL equations with S in (50). We are
denoting these IMEX schemes, with source terms S, as follows:

  1. 1.

O2ES-IMEX: We consider second-order time stepping with MinMod reconstruction
for scaled entropy variables. We use the second-order central difference for
non-conservative terms. Source term S is treated implicitly.

  2. 2.

O3ES-IMEX: We consider third-order time stepping with third-order ENO
reconstruction for scaled entropy variables. We use the fourth-order central
difference for non-conservative terms. Source term S is treated implicitly.

  3. 3.

O4ES-IMEX: We consider fourth-order time stepping with fourth-order ENO
reconstruction for scaled entropy variables. We use the fourth-order central
difference for non-conservative terms. Source term S is treated implicitly.

##  6 Numerical results

In this Section, we present several test cases in one dimension. We also
present a two-dimensional Orzag-Tang vertex test case, which is generalized
from the MHD test case. The time step is computed using

$\Delta
t=\text{CFL}\frac{1}{\max_{i,j}\left(\frac{|\lambda_{x}(\textbf{U}_{i,j})|}{\Delta
x}+\frac{|\lambda_{x}(\textbf{U}_{i,j})|}{\Delta y}\right)},$

where, $\lambda_{x}$ is the maximum eigenvalue in $x$-direction and
$\lambda_{y}$ is the maximum eigenvalue in $y$-direction. To measure the
entropy decay for each test case, we also measure

$\displaystyle\sum_{i,j}\Big{[}\mathcal{E}^{n+1}(\textbf{U}_{i,j})-\mathcal{E}^{n}(\textbf{U}_{i,j})-\frac{\Delta
t}{\Delta
x}\left(\hat{\mathcal{Q}}^{k}_{x,i+\frac{1}{2},j}-\hat{\mathcal{Q}}^{k}_{x,i-\frac{1}{2},j}\right)$
$\displaystyle-\frac{\Delta t}{\Delta
y}\left(\hat{\mathcal{Q}}^{k}_{y,i,j+\frac{1}{2}}-\hat{\mathcal{Q}}^{k}_{y,i,j-\frac{1}{2}}\right)\Big{]},$
(56)

for explicit schemes (see Equation (49)). This measures the total entropy
change after every time step. Similarly, for the implicit scheme, we measure

$\displaystyle\sum_{i,j}\Big{[}\mathcal{E}^{n+1}(\textbf{U}_{i,j})-\mathcal{E}^{n}(\textbf{U}_{i,j})$
$\displaystyle-\frac{\Delta t}{\Delta
x}\left(\hat{\mathcal{Q}}^{k}_{x,i+\frac{1}{2},j}-\hat{\mathcal{Q}}^{k}_{x,i-\frac{1}{2},j}\right)$
$\displaystyle-\frac{\Delta t}{\Delta
y}\left(\hat{\mathcal{Q}}^{k}_{y,i,j+\frac{1}{2}}-\hat{\mathcal{Q}}^{k}_{y,i,j-\frac{1}{2}}\right)-\Delta
t\textbf{V}_{i,j}\textbf{S}(\textbf{U}_{i,j})\Big{]}.$ (57)

which is consistent with Eqn. (52), which also takes source term into
consideration. When computing with the source term, we expect pressure to be
isotropic (i.e. $p_{\parallel}$ and $p_{\perp}$ should be equal); hence, the
solution should be closer to the MHD limit. To compare with the MHD solution,
we have plotted the corresponding MHD solution, which is computed using
Rusanov solver, with MinMod reconstruction. In one-dimensional test cases, we
have used $10000$ cells. For the Orszag-Tang vertex two-dimensional test case,
we have used $1000\times 1000$ cells to compute the reference MHD solution.

###  6.1 One-dimensional test problems

####  6.1.1 Accuracy test

To assess the accuracy of the proposed schemes, we consider a problem with a
smooth solution. The computational domain is taken to be $[0,1]$ with periodic
boundary conditions. Initially, the density profile is assumed to be
$\rho(x,0)=2+\sin{(2\pi x)}$ which is then advected with velocity
$\textbf{u}=(1,0,0)$. The pressures are kept constant at
$p_{\parallel}=p_{\perp}=1.0$. We also take $\boldsymbol{B}=(1,1,0)$. The
exact solution is an advection of the density profile given by
$\rho(x,t)=2+\sin{(2\pi(x-t))}$.

Table 2: Accuracy test: $L^{1}$ errors and order of accuracy for $\rho$. Number of cells | O2ES-EXP | O3ES-EXP | O4ES-EXP
---|---|---|---
|  $L^{1}$ error | Order |  $L^{1}$ error | Order |  $L^{1}$ error | Order
40 | 1.30241e-01 | – | 5.53282e-03 | – | 8.98990e-04 | –
80 | 5.17593e-02 | 1.33 | 7.07952e-04 | 2.96629 | 8.19842e-05 | 3.45489
160 | 1.65341e-02 | 1.65 | 8.84319e-05 | 3.00101 | 5.86817e-06 | 3.80436
320 | 4.57448e-03 | 1.86 | 1.10535e-05 | 3.00007 | 4.11943e-07 | 3.83239
640 | 1.24571e-03 | 1.88 | 1.38170e-06 | 2.99998 | 2.84900e-08 | 3.85391
1280 | 3.35275e-04 | 1.89 | 1.72720e-07 | 2.99993 | 1.92548e-09 | 3.88717
Table 3: Accuracy test: $L^{1}$ errors and order of accuracy for $\rho$. Number of cells | O2ES-IMEX | O3ES-IMEX | O4ES-IMEX
---|---|---|---
|  $L^{1}$ error | Order |  $L^{1}$ error | Order |  $L^{1}$ error | Order
40 | 1.30187e-01 | – | 5.57809e-03 | – | 8.90421e-04 | –
80 | 5.16249e-02 | 1.33 | 7.19948e-04 | 2.95381 | 8.24102e-05 | 3.43359
160 | 1.65536e-02 | 1.64 | 9.01302e-05 | 2.99781 | 5.96265e-06 | 3.7888
320 | 4.58388e-03 | 1.85 | 1.12664e-05 | 2.99999 | 4.19101e-07 | 3.83059
640 | 1.24898e-03 | 1.88 | 1.40835e-06 | 2.99994 | 2.90423e-08 | 3.85107
1280 | 3.36345e-04 | 1.89 | 1.76049e-07 | 2.99996 | 1.95243e-09 | 3.89482

The simulations are conducted till the final time $t=2.0$. Table (2) presents
the $L^{1}$-errors of the density for the O2ES-EXP, O3ES-EXP and O4ES-EXP
schemes. We observe that all the schemes have the expected order of accuracy.
Similarly, in Table (3), we present the $L^{1}$-errors of the density for the
O2ES-IMEX, O3ES-IMEX and O4ES-IMEX schemes. Again, we see that all the schemes
have a theoretical order of accuracy.

####  6.1.2 Riemann problem 1: Brio-Wu shock tube problem

This test case is a generalized form of the MHD test case proposed by Brio and
Wu Brio1988 . We consider the computational domain of $[-1,1]$ with outflow
boundary conditions. The initial discontinuity is centred at $x=0$, with
initial states given by,

$(\rho,\textbf{u},p_{\parallel},p_{\perp},B_{y},B_{z})=\begin{cases}(1,0,0,0,1,1,1,0),&\textrm{if
}x\leq 0\\\ (0.125,0,0,0,0.1,0.1,-1,0),&\textrm{otherwise}\end{cases}$

We consider $B_{x}=0.75$ and perform the computations till final time $t=0.2$
on $2000$ cells.

(a) Anisotropic case: Density for explicit schemes without source term

(b) Anisotropic case: Pressure component $p_{\parallel}$ for explicit schemes
without source term

(c) Anisotropic case: Pressure component $p_{\perp}$ for explicit schemes
without source term

(d) Isotropic case: Density for IMEX schemes with source term

(e) Isotropic case: Pressure component $p_{\parallel}$ for IMEX schemes with
source term

(f) Isotropic case: Pressure component $p_{\perp}$ for IMEX schemes with
source term

(g) Anisotropic case: Total entropy decay at each time step for explicit
schemes

(h) Isotropic case: Total entropy decay at each time step for IMEX schemes
with source term

Figure 1: Riemann problem 1: Brio-Wu shock tube problem: Plots of density,
parallel and perpendicular pressure components, and total entropy decay at
each time step for explicit schemes without source term and IMEX scheme with
source term using $2000$ cells at final time $t=0.2$.

Numerical results are presented in Figure (1). In Figures (1a), (1b) and (1c),
we have presented the density, parallel component of pressure, and
perpendicular component of the pressure using explicit schemes. We observe all
the schemes have resolved all the waves. We also note that both components of
the pressures are different, as expected, because we are computing anisotropic
solutions. The results are similar to the anisotropic solution given in
Hirabayashi2016new . We also note that O2ES-EXP is most diffusive compared to
O3ES-EXP and O4ES-EXP.

The isotropic solutions (with source terms) with IMEX schemes are presented in
Figures (1d), (1e) and (1f). We observe that the profile of $p_{\parallel}$
and $p_{\perp}$ are the same due to the effect of the source term.
Furthermore, the solution for all the components matches the reference MHD
solution. This is similar to the isotropic solution presented in
Hirabayashi2016new . For the IMEX schemes, we again observe that the O2ES-IMEX
is the most diffusive, with O3ES-IMEX and O4ES-IMEX being comparable.

In Figures (1g) and (1h), we have plotted the total entropy change at each
time step for explicit and IMEX schemes, respectively. For schemes to be
entropy stable, we need these values to be less than zero. At the same time,
we desire that the schemes do not decay too much entropy if no new
discontinuity is generated in the solution. We see that all schemes decay lots
of entropy initially, as we have several discontinuities being generated.
After that, as no new discontinuity is generated, all the schemes decay much
less entropy. Again, we see that both second-order schemes, O2ES-EXP and O2ES-
IMEX are most diffusive as they decay most entropy, followed by third-order
and fourth-order schemes. We also see that there is not much difference in the
entropy decay of explicit and IMEX schemes of the same order.

####  6.1.3 Riemann problem 2: Ryu-Jones problem

In this test case, we consider a generalization of the Ryu–Jones test problem
Ryu1995Numerical . The computational domain is taken to be $[-0.5,0.5]$, with
outflow boundary conditions. The initial conditions are given by,

$(\rho,\textbf{u},p_{\parallel},p_{\perp},B_{y},B_{z})=\begin{cases}(1.08,1.2,0,0,0.95,0.95,\frac{3.6}{\sqrt{4\pi}},\frac{2}{\sqrt{4\pi}}),&\textrm{if
}x\leq 0\\\
(1,0,0,0,1,1,\frac{4}{\sqrt{4\pi}},\frac{2}{\sqrt{4\pi}}),&\textrm{otherwise}\end{cases}$

with $B_{x}=\frac{2}{\sqrt{4\pi}}$. We compute the solutions using $2000$
cells till the final time $t=0.2$.

(a) Anisotropic case: Density for explicit schemes without source term

(b) Anisotropic case: Pressure component $p_{\parallel}$ for explicit schemes
without source term

(c) Anisotropic case: Pressure component $p_{\perp}$ for explicit schemes
without source term

(d) Isotropic case: Density for IMEX schemes with source term

(e) Isotropic case: Pressure component $p_{\parallel}$ for IMEX schemes with
source term

(f) Isotropic case: Pressure component $p_{\perp}$ for IMEX schemes with
source term

(g) Anisotropic case: Total entropy decay at each time step for explicit
schemes

(h) Isotropic case: Total entropy decay at each time step for IMEX schemes
with source term

Figure 2: Riemann problem 2: Ryu-Jones problem: Plots of density, parallel and
perpendicular pressure components, and total entropy decay at each time step
for explicit schemes without source term and IMEX scheme with source term
using $2000$ cells at final time $t=0.2$.

Numerical results are presented in Figure (2). Anisotropic solutions (without
source term) using explicit schemes are presented in Figures (2a), (2b) and
(2c). We observe that the profiles of $p_{\parallel}$ and $p_{\perp}$ are
different. We see that all the schemes are able to resolve different waves. We
do see some small oscillations near some discontinuities. However, these
oscillations are stable and do not increase when we refine the mesh. In Figure
(2g), we have plotted the total entropy decay for each scheme. We again see
that O2ES-EXP is most diffusive when compared with O3ES-EXP and O4ES-EXP
schemes.

Isotropic solutions (with source term) using IMEX schemes are presented in
(2d), (2e) and (2f). We see that all the variables have converged to the MHD
solution. Furthermore, both components of pressure have the same profile, and
the solution is isotropic. We again observe that O2ES-IMEX is the most
diffusive, followed by O3ES-IMEX and O4ES-IMEX schemes. In Figure (2h), we
have plotted the total entropy decay at each time step. We see that the
entropy decays are comparable to the explicit schemes.

####  6.1.4 Riemann problem 3: Super-fast expansion

In this test, we consider the generalization of the MHD super-fast expansion
test case, considered in bouchut2007multiwave ; fuchs2009splitting ;
fuchs2011approximate . The solution contains a low-density, low-pressure area
in the middle as the fluid is pushed outward. We consider the computational
domain of $[0,1]$ with outflow boundary conditions. The initial conditions are
given by,

$(\rho,\textbf{u},p_{\parallel},p_{\perp},B_{x},B_{y},B_{z})=\begin{cases}(1,-3.1,0,0,1,1,0,0.5,0),&\textrm{if
}x\leq 0.5\\\ (1,3.1,0,0,1,1,0,0.5,0),&\textrm{otherwise.}\end{cases}$

(a) Anisotropic case: Density for explicit schemes without source term

(b) Anisotropic case: Pressure component $p_{\parallel}$ for explicit schemes
without source term

(c) Anisotropic case: Pressure component $p_{\perp}$ for explicit schemes
without source term

(d) Isotropic case: Density for IMEX schemes with source term

(e) Isotropic case: Pressure component $p_{\parallel}$ for IMEX schemes with
source term

(f) Isotropic case: Pressure component $p_{\perp}$ for IMEX schemes with
source term

(g) Anisotropic case: Total entropy decay at each time step for explicit
schemes

(h) Isotropic case: Total entropy decay at each time step for IMEX schemes
with source term

Figure 3: Riemann problem 3: Super-fast expansion: Plots of density, parallel
and perpendicular pressure components, and total entropy decay at each time
step for explicit schemes without source term and IMEX scheme with source term
using $2000$ cells at final time $t=0.05$.

The numerical results are plotted in Figure (3) at the final time of $t=0.05$
using $2000$ cells. In Figures (3a), (3b) and (3c), we have plotted density,
$p_{\parallel}$ and $p_{\perp}$ for anisotropic case using explicit schemes.
We again see the low-density, low-pressure area is resolved and both
components of the pressure are not identical. Similarly, in Figures (3d), (3e)
and (3f), we have plotted isotropic solutions using IMEX schemes. We observe
that the profiles of pressure components are the same and match the MHD
solution. Entropy decays for explicit and IMEX schemes are plotted in Figures
(3g) and (3h), respectively. We note that the second order schemes O2ES-EXP
and O2ES-IMEX are most diffusive compared to the higher order schemes O3ES-
EXP, O3ES-IMEX, O4ES-EXP and O4ES-IMEX. Furthermore, both explicit and IMEX
schemes of corresponding order have similar entropy decay performance.

####  6.1.5 Riemann problem 4

In this test, we generalize a Riemann problem described in dumbser2016new .
The computational domain is $[-0.5,0.5]$, with outflow boundary conditions.
The initial profile is given by,

$(\rho,\textbf{u},p_{\parallel},p_{\perp},B_{y},B_{z})=\begin{cases}(1,0,0,0,1,1,1,0),&\textrm{if
}x\leq 0\\\ (0.4,0,0,0,0.4,0.4,-1,0),&\textrm{otherwise}\end{cases}$

with $B_{x}=1.3$. We compute the solution using $2000$ cells till the final
time $t=0.15$. The numerical results are presented in Figure (4). Anisotropic
density, $p_{\parallel}$, $p_{\perp}$ and $B_{y}$ using explicit schemes are
plotted in Figure (4a), (4b), (4c) and (4d), respectively. We observe that
both pressure components have significantly different profiles. Isotropic
density, $p_{\parallel}$, $\perp$ and $B_{y}$ using explicit schemes are
plotted in Figure (4e), (4f), (4g) and (4h), respectively. We see that all the
variables are matching the MHD solution. We observe that the second-order
schemes are less accurate when compared with higher-order schemes. We also
observe some small-scale oscillations near some of the discontinuities.
However, they are stable and do not increase with further mesh refinement.

Total entropy decays at each time step are plotted in Figures (4i) and (4j)
for explicit and IMEX schemes, respectively. We again see that both explicit
and IMEX schemes have similar entropy decay performance.

(a) Anisotropic case: Density for explicit schemes without source term

(b) Anisotropic case: Pressure component $p_{\parallel}$ for explicit schemes
without source term

(c) Anisotropic case: Pressure component $p_{\perp}$ for explicit schemes
without source term

(d) Anisotropic case: Magnetic field component $B_{y}$ for explicit schemes
without source term

(e) Isotropic case: Density for IMEX schemes with source term

(f) Isotropic case: Pressure component $p_{\parallel}$ for IMEX schemes with
source term

(g) Isotropic case: Pressure component $p_{\perp}$ for IMEX schemes with
source term

(h) Isotropic case: Magnetic field component $B_{y}$ for explicit schemes
without source term

(i) Anisotropic case: Total entropy decay at each time step for explicit
schemes

(j) Isotropic case: Total entropy decay at each time step for IMEX schemes
with source term

Figure 4: Riemann problem 4: Plots of density, parallel and perpendicular
pressure components, and total entropy decay at each time step for explicit
schemes without source term and IMEX scheme with source term using $2000$
cells at final time $t=0.15$.

####  6.1.6 Riemann problem 5

In this test, we consider another Riemann problem from dumbser2016new . The
initial conditions are given by,

$(\rho,\textbf{u},p_{\parallel},p_{\perp},B_{y},B_{z})=\begin{cases}(1.7,0,0,0,1.7,1.7,\frac{3.544908}{\sqrt{4\pi}},0),&\textrm{if
}x\leq 0\\\
(0.2,0,0,-1.496891,0.2,0.2,\frac{2.785898}{\sqrt{4\pi}},\frac{2.192064}{\sqrt{4\pi}}),&\textrm{otherwise}\end{cases}$

with $B_{x}=\frac{3.899398}{\sqrt{4\pi}}$. We consider the computational
domain $[-0.5,0.5]$ with outflow boundary conditions. We use $2000$ cells and
compute the solutions till the final time $t=0.15$.

(a) Anisotropic case: Density for explicit schemes without source term

(b) Anisotropic case: Pressure component $p_{\parallel}$ for explicit schemes
without source term

(c) Anisotropic case: Pressure component $p_{\perp}$ for explicit schemes
without source term

(d) Anisotropic case: Magnetic field component $B_{y}$ for explicit schemes
without source term

(e) Isotropic case: Density for IMEX schemes with source term

(f) Isotropic case: Pressure component $p_{\parallel}$ for IMEX schemes with
source term

(g) Isotropic case: Pressure component $p_{\perp}$ for IMEX schemes with
source term

(h) Isotropic case: Magnetic field component $B_{y}$ for explicit schemes
without source term

(i) Anisotropic case: Total entropy decay at each time step for explicit
schemes

(j) Isotropic case: Total entropy decay at each time step for IMEX schemes
with source term

Figure 5: Riemann problem 5: Plots of density, parallel and perpendicular
pressure components, and total entropy decay at each time step for explicit
schemes without source term and IMEX scheme with source term using $2000$
cells at final time $t=0.15$.

Numerical results are presented in Figure (5) for anisotrpic solutions using
IMEX schemes O2ES-EXP, O3ES-EXP and O4ES-EXPschemes (see Figures (5a), (5b),
(5c), and (5d) )and for isotropic solutions using O2ES-IMEX, O3ES-IMEX and
O4ES-IMEX scheme (see Figures (5e), (5f), (5g), and (5h)). We have plotted
density, pressure components $p_{\parallel}$ and $p_{\perp}$ and magnetic
filed component $B_{y}$. We observe that all the schemes are able to resolve
different waves. Furthermore, isotropic solutions match MHD solutions for all
the variables. From the entropy decay in Figures (5i) and (5j), we observe
that the second-order schemes are more diffusive than the third and fourth-
order schemes.

####  6.1.7 Riemann problem 6

In this test, we consider the Riemann problem presented in
Balsara2018efficient  for MHD equations. The computational domain is
considered to be $[-0.5,0.5]$ with outflow boundary conditions. The initial
conditions are,

$(\rho,\textbf{u},p_{\parallel},p_{\perp},B_{y},B_{z})=\begin{cases}(\frac{1}{4\pi},-1,1,-1,1,1,-\frac{1}{\sqrt{4\pi}},\frac{1}{\sqrt{4\pi}}),&\textrm{if
}x\leq 0\\\
(\frac{1}{4\pi},-1,-1,-1,1,1,\frac{1}{\sqrt{4\pi}},\frac{1}{\sqrt{4\pi}}),&\textrm{otherwise}\end{cases}$

with magnetic field component $B_{x}=\frac{1}{\sqrt{4\pi}}$. We compute
anisotropic and isotropic solutions using explicit (O2ES-EXP, O3ES-EXP and
O4ES-EXP) and IMEX (O2ES-IMEX, O3ES-IMEX and O4ES-IMEX) schemes, respectively
with $2000$ cells. The numerical results for variable $p_{\parallel}$,
$p_{\perp}$ and $B_{y}$ are presented in Figure (6). We observe that the wave
structures are much more dominant in the anisotropic case than in the
isotropic case. Furthermore, in the case of anisotropic solutions, we clearly
observe that O4ES-EXP is most accurate, followed by O3ES-EXP and O2ES-EXP
schemes (see Figures (6a), (6b) and (6c) ). Furthermore, the isotropic
solutions matched the MHD solution for IMEX schemes see Figures (6d), (6e) and
(6f). We also note that entropy decays are more oscillatory (see Figures (6g)
and (6h)).

(a) Anisotropic case: Pressure component $p_{\parallel}$ for explicit schemes
without source term

(b) Anisotropic case: Pressure component $p_{\perp}$ for explicit schemes
without source term

(c) Anisotropic case: Magnetic field component $B_{y}$ for explicit schemes
without source term

(d) Isotropic case: Pressure component $p_{\parallel}$ for IMEX schemes with
source term

(e) Isotropic case: Pressure component $p_{\perp}$ for IMEX schemes with
source term

(f) Isotropic case: Magnetic field component $B_{y}$ for explicit schemes
without source term

(g) Anisotropic case: Total entropy decay at each time step for explicit
schemes

(h) Isotropic case: Total entropy decay at each time step for IMEX schemes
with source term

Figure 6: Riemann problem 6: Plots of density, parallel and perpendicular
pressure components, and total entropy decay at each time step for explicit
schemes without source term and IMEX scheme with source term using $2000$
cells at final time $t=0.15$.

####  6.1.8 Riemann problem 7

The last one-dimensional test case for CGL system is designed using MHD
Riemann problem from Balsara2018efficient . The initial conditions are given
by,

$(\rho,\textbf{u},p_{\parallel},p_{\perp},B_{y},B_{z})=\begin{cases}(1,0,0,0,1,1,1,0),&\textrm{if
}x\leq 0\\\ (0.2,0,0,0,0.1,0.1,0,0),&\textrm{otherwise},\end{cases}$

with $B_{x}=1$. We take the computational domain of $[-0.5,0.5]$ with outflow
boundary conditions. We compute the solution till the final time $t=0.15$. The
numerical solutions using explicit and IMEX schemes for anisotropic and
isotropic cases are presented in Figure (7) using $2000$ cells. We again
observe that all the schemes are able to resolve wave structures for both
anisotropic and isotropic cases. We also observed a significant change in the
wave structure when comparing anisotropic and isotropic cases. For the
isotropic case, the solutions match with the MHD solution. For the entropy
decay, we observe that explicit and IMEX schemes have similar decays.

(a) Anisotropic case: Density for explicit schemes without source term

(b) Anisotropic case: Pressure component $p_{\parallel}$ for explicit schemes
without source term

(c) Anisotropic case: Pressure component $p_{\perp}$ for explicit schemes
without source term

(d) Anisotropic case: Magnetic field component $B_{y}$ for explicit schemes
without source term

(e) Isotropic case: Density for IMEX schemes with source term

(f) Isotropic case: Pressure component $p_{\parallel}$ for IMEX schemes with
source term

(g) Isotropic case: Pressure component $p_{\perp}$ for IMEX schemes with
source term

(h) Isotropic case: Magnetic field component $B_{y}$ for explicit schemes
without source term

(i) Anisotropic case: Total entropy decay at each time step for explicit
schemes

(j) Isotropic case: Total entropy decay at each time step for IMEX schemes
with source term

Figure 7: Riemann problem 7: Plots of density, parallel and perpendicular
pressure components, and total entropy decay at each time step for explicit
schemes without source term and IMEX scheme with source term using $2000$
cells at final time $t=0.15$.

###  6.2 Two-dimensional test problems

####  6.2.1 Orszag–Tang vortex

To test the stability of the proposed schemes in two dimensions, we generalize
the Orzag-Tang vortex test case of MHD equations, which was first considered
in orszag1979small . Following toth2000b ; chandrashekar2016entropy , we
consider the computational domain $[0,1]\times[0,1]$ with periodic boundary
conditions. The initial conditions are specified as follows:

$\displaystyle\rho=\frac{25}{36\pi},~{}~{}~{}~{}\textbf{u}$
$\displaystyle=(-\sin{(2\pi y)},\sin{(2\pi
x)},0),~{}~{}~{}~{}p_{\perp}=p_{\parallel}=\frac{5}{12\pi},$
$\displaystyle\boldsymbol{B}$ $\displaystyle=\frac{1}{\sqrt{4\pi}}(-\sin{(2\pi
y)},\sin{(4\pi x)},0).$

We present the solutions for the isotropic case using IMEX schemes only, as we
do not have any comparable results for the anisotropic case. The numerical
results for O2ES-IMEX, O3ES-IMEX and O4ES-IMEX schemes at $200\times 200$,
$400\times 400$ and $800\times 800$ cells are presented in Figure (8). We
observe that all the schemes are able to resolve the complicated flow
features. To compare the solution with MHD solution, we have plotted the cut
of variables $p_{\parallel}$ and $p_{\perp}$ for $200\times 200$ cells and
compared the solutions with MHD pressure (see Figures (9a) and (9b) ). We
observe that all the schemes produce solutions similar to the MHD solution,
with O4ES-IMEX being the most accurate and O2ES-IMEX being the most diffusive.

The solution is initially smooth and, over time, develops discontinuities.
Hence, low entropy decay should be expected initially, and as the
discontinuities form, entropy decay should increase. We observe the same
behaviour of the schemes in Figure (9c), (9d) and (9e). We also note that the
second-order scheme has the most entropy decay.

(a) Density using $200\times 200$ cells for O2ES-IMEXscheme

(b) Density using $400\times 400$ cells for O2ES-IMEXscheme

(c) Density using $800\times 800$ cells for O2ES-IMEXscheme

(d) Density using $200\times 200$ cells for O3ES-IMEXscheme

(e) Density using $400\times 400$ cells for O3ES-IMEXscheme

(f) Density using $800\times 800$ cells for O3ES-IMEXscheme

(g) Density using $200\times 200$ cells for O4ES-IMEXscheme

(h) Density using $400\times 400$ cells for O4ES-IMEX4cheme

(i) Density using $800\times 800$ cells for O4ES-IMEXscheme

Figure 8: Orszag–Tang vortex: Isotropic case: Plots of density at different
meshes for IMEX schemes O2ES-IMEX, O3ES-IMEX and O4ES-IMEX at time $t=0.5$.

(a) Isotropic case: Cut of pressure component $p_{\parallel}$ along $y=0.3125$
with $200\times 200$ cells using O2ES-IMEX, O3ES-IMEXand O4ES-IMEXschemes

(b) Isotropic case: Cut of pressure component $p_{\perp}$ along $y=0.3125$
with $200\times 200$ cells using O2ES-IMEX, O3ES-IMEXand O4ES-IMEXschemes

(c) Isotropic case: Total entropy decay at each time step with $200\times 200$
cells using O2ES-IMEX, O3ES-IMEXand O4ES-IMEXschemes

(d) Isotropic case: Total entropy decay at each time step with $400\times 400$
cells using O2ES-IMEX, O3ES-IMEXand O4ES-IMEXschemes

(e) Isotropic case: Total entropy decay at each time step with $800\times 800$
cells using O2ES-IMEX, O3ES-IMEXand O4ES-IMEXschemes

Figure 9: Orszag–Tang vortex: Plots of pressure components and total entropy
decay at each time step for IMEX schemes at time $t=0.5$.

##  7 Conclusion

In this article, we have proposed higher-order finite difference entropy
stable numerical schemes for the CGL equations, which form a system of non-
conservative hyperbolic equations. The proposed schemes are based on the novel
reformulation of the equations, such that the non-conservative terms do not
contribute to the entropy evolution. We then present the Godunov
symmetrization of the conservative part. To design entropy stable schemes for
the complete system, we first consider the entropy conservative numerical
schemes for the conservative part and then extend the same to include non-
conservative terms. To ensure entropy stability, we have considered the
diffusion operator based on the entropy scaled right eigenvector of the
symmetrized conservative part. The time discretization is performed using SSP-
RK schemes for explicit schemes, which are used to calculate anisotropic
solutions. For the isotropic case, we have considered a still source term.
Hence, we have used ARK-IMEX schemes for time discretization. To test the
entropy stability of the proposed numerical schemes, we have designed several
test cases in one dimension by generalizing the one-dimensional MHD Riemann
problems. We demonstrate that for a wide variety of test cases, the numerical
entropy decay holds for each time step for both anisotropic and isotropic
limits. For the isotropic limit, we show that solutions matched the MHD
solution, as expected. For the two-dimensional test case, we have considered
the Orzag-Tang vortex test case and present the isotropic case. We demonstrate
that the schemes decay entropy only when discontinuities start to form in the
solution. Furthermore, the solutions match the MHD solutions.

## Acknowledgments

The work of Harish Kumar is supported in parts by VAJRA grant No.
VJR/2018/000129 by the Dept. of Science and Technology, Govt. of India

## Conflict of interest

The authors declare that they have no conflict of interest.

## Data Availability Declaration

Data will be made available on reasonable request.

## References

  * (1) Abgrall, R., Karni, S.: A comment on the computation of non-conservative products.  Journal of Computational Physics 229(8), 2759–2763 (2010)
  * (2) Abgrall, R., Kumar, H.: Robust Finite Volume Schemes for Two-Fluid Plasma Equations.  Journal of Scientific Computing 60(3), 584–611 (2014)
  * (3) Balsara, D.S.: The eigenstructure of the equations of radiation magnetohydrodynamics.  Journal of Quantitative Spectroscopy and Radiative Transfer 61(5), 637–646 (1999)
  * (4) Balsara, D.S., Li, J., Montecinos, G.I.: An efficient, second-order accurate, universal generalized Riemann problem solver based on the HLLI Riemann solver.  Journal of Computational Physics 375, 1238–1269 (2018)
  * (5) Barth, T.: On the role of involutions in the discontinuous galerkin discretization of maxwell and magnetohydrodynamic systems.  IMA VOLUMES IN MATHEMATICS AND ITS APPLICATIONS 142, 69 (2006)
  * (6) Barth, T.J.: Numerical methods for gasdynamic systems on unstructured meshes.  An Introduction to Recent Developments in Theory and Numerics for Conservation Laws: Proceedings of the International School on Theory and Numerics for Conservation Laws, Freiburg/Littenweiler, October 20–24, 1997 pp. 195–285 (1999)
  * (7) Berthon, C.: Numerical approximations of the 10-moment gaussian closure.  Mathematics of computation 75(256), 1809–1831 (2006)
  * (8) Berthon, C., Dubroca, B., Sangam, A.: An entropy preserving relaxation scheme for ten-moments equations with source terms.  Communications in Mathematical sciences 13(8), 2119–2154 (2015)
  * (9) Bhoriya, D., Balsara, D., Florinski, V., Kumar, H.: Going beyond the mhd approximation: Physics-based numerical solution of the cgl equations.  Submitted (2024)
  * (10) Bhoriya, D., Kumar, H., Chandrashekar, P.: High-order finite-difference entropy stable schemes for two-fluid relativistic plasma flow equations.  Journal of Computational Physics p. 112207 (2023)
  * (11) Biswas, B., Kumar, H., Yadav, A.: Entropy stable discontinuous galerkin methods for ten-moment gaussian closure equations.  Journal of Computational Physics 431, 110148 (2021)
  * (12) Bouchut, F., Klingenberg, C., Waagan, K.: A multi-wave approximate riemann solver for ideal mhd based on relaxation. i: theoretical framework.  Numerische Mathematik 108, 7–42 (2007)
  * (13) Brio, M., Wu, C.: An upwind differencing scheme for the equations of ideal magnetohydrodynamics.  Journal of Computational Physics 75(2), 400–422 (1988)
  * (14) Brown, S., Roe, P., Groth, C.: Numerical solution of a 10-moment model for nonequilibrium gasdynamics.  In: 12th Computational Fluid Dynamics Conference, p. 1677 (1995)
  * (15) Burch, J., Moore, T., Torbert, R., Giles, B.h.: Magnetospheric multiscale overview and science objectives.  Space Science Reviews 199, 5–21 (2016)
  * (16) Chandran, B.D., Dennis, T.J., Quataert, E., Bale, S.D.: Incorporating kinetic physics into a two-fluid solar-wind model with temperature anisotropy and low-frequency alfvén-wave turbulence.  The Astrophysical Journal 743(2), 197 (2011)
  * (17) Chandrashekar, P.: Kinetic energy preserving and entropy stable finite volume schemes for compressible euler and navier-stokes equations.  Communications in Computational Physics 14(5), 1252–1286 (2013)
  * (18) Chandrashekar, P., Klingenberg, C.: Entropy stable finite volume scheme for ideal compressible mhd on 2-d cartesian meshes.  SIAM Journal on Numerical Analysis 54(2), 1313–1340 (2016)
  * (19) Chew, G., Goldberger, M., Low, F.: The boltzmann equation and the one-fluid hydromagnetic equations in the absence of particle collisions.  Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 236(1204), 112–118 (1956)
  * (20) Dal Maso, G., Lefloch, P.G., Murat, F.: Definition and weak stability of nonconservative products.  Journal de mathématiques pures et appliquées 74(6), 483–548 (1995)
  * (21) Dennis Jr, J.E., Schnabel, R.B.: Numerical methods for unconstrained optimization and nonlinear equations.  SIAM (1996)
  * (22) Dubroca, B., Tchong, M., Charrier, P., Tikhonchuk, V., Morreeuw, J.P.: Magnetic field generation in plasmas due to anisotropic laser heating.  Physics of Plasmas 11(8), 3830–3839 (2004)
  * (23) Dumbser, M., Balsara, D.S.: A new efficient formulation of the hllem riemann solver for general conservative and non-conservative hyperbolic systems.  Journal of Computational Physics 304, 275–319 (2016)
  * (24) Dumin, Y.V.: The corotation field in collisionless magnetospheric plasma and its influence on average electric field in the lower atmosphere.  Advances in Space Research 30(10), 2209–2214 (2002)
  * (25) Fjordholm, U.S., Mishra, S., Tadmor, E.: Arbitrarily high-order accurate entropy stable essentially nonoscillatory schemes for systems of conservation laws.  SIAM Journal on Numerical Analysis 50(2), 544–573 (2012)
  * (26) Fjordholm, U.S., Mishra, S., Tadmor, E.: Eno reconstruction and eno interpolation are stable.  Foundations of Computational Mathematics 13(2), 139–159 (2013)
  * (27) Fuchs, F.G., McMurry, A.D., Mishra, S., Risebro, N.H., Waagan, K.: Approximate riemann solvers and robust high-order finite volume schemes for multi-dimensional ideal mhd equations.  Communications in Computational Physics 9(2), 324–362 (2011)
  * (28) Fuchs, F.G., Mishra, S., Risebro, N.H.: Splitting-based finite volume schemes for ideal mhd equations.  Journal of Computational Physics 228(3), 641–660 (2009)
  * (29) Godlewski, E., Raviart, P.A.: Numerical approximation of hyperbolic systems of conservation laws, vol. 118.  Springer Science & Business Media (2013)
  * (30) Godunov, S.K.: Symmetric form of the equations of magnetohydrodynamics.  Numerical Methods for Mechanics of Continuum Medium 1, 26–34 (1972)
  * (31) Goedbloed, J.H., Poedts, S.: Principles of magnetohydrodynamics: with applications to laboratory and astrophysical plasmas.  Cambridge university press (2004)
  * (32) Gombosi, T.I., Rasmussen, C.E.: Transport of gyration-dominated space plasmas of thermal origin: 1. generalized transport equations.  Journal of Geophysical Research: Space Physics 96(A5), 7759–7778 (1991)
  * (33) Gottlieb, S., Shu, C.W., Tadmor, E.: Strong stability-preserving high-order time discretization methods.  SIAM review 43(1), 89–112 (2001)
  * (34) Hakim, A., Loverich, J., Shumlak, U.: A high resolution wave propagation scheme for ideal Two-Fluid plasma equations.  Journal of Computational Physics 219(1), 418–442 (2006)
  * (35) Hakim, A.H.: Extended mhd modelling with the ten-moment equations.  Journal of Fusion Energy 27, 36–43 (2008)
  * (36) Heinemann, M.: Role of collisionless heat flux in magnetospheric convection.  Journal of Geophysical Research: Space Physics 104(A12), 28397–28410 (1999)
  * (37) Hirabayashi, K., Hoshino, M., Amano, T.: A new framework for magnetohydrodynamic simulations with anisotropic pressure.  Journal of Computational Physics 327, 851–872 (2016)
  * (38) Hollweg, J.V.: Collisionless electron heat conduction in the solar wind.  Journal of Geophysical Research 81(10), 1649–1658 (1976)
  * (39) Huang, Z., Tóth, G., van der Holst, B., Chen, Y., Gombosi, T.: A six-moment multi-fluid plasma model.  Journal of Computational Physics 387, 134–153 (2019)
  * (40) Hunana, P., Tenerani, A., Zank, G., Khomenko, E., Goldstein, M., Webb, G., Velli, M., Adhikari, L.: A brief guide to fluid models with anisotropic temperatures part 1-cgl description and collisionless fluid hierarchy.  arXiv preprint arXiv:1901.09354 (2019)
  * (41) Ismail, F., Roe, P.L.: Affordable, entropy-consistent euler flux functions ii: Entropy production at shocks.  Journal of Computational Physics 228(15), 5410–5436 (2009)
  * (42) Johnson, E.A.: Gaussian-moment relaxation closures for verifiable numerical simulation of fast magnetic reconnection in plasma.  arXiv preprint arXiv:1409.6985 (2014)
  * (43) Karimabadi, H., Roytershteyn, V., Vu, H., Omelchenko, Y., Scudder, J., Daughton, W., Dimmock, A., Nykyri, K., Wan, M., Sibeck, D., et al.: The link between shocks, turbulence, and magnetic reconnection in collisionless plasmas.  Physics of Plasmas 21(6) (2014)
  * (44) Kato, Y., Tajiri, M., Taniuti, T.: Propagation of hydromagnetic waves in collisionless plasma. i.  Journal of the Physical Society of Japan 21(4), 765–777 (1966)
  * (45) Kennedy, C.A., Carpenter, M.H.: Additive runge–kutta schemes for convection–diffusion–reaction equations.  Applied numerical mathematics 44(1-2), 139–181 (2003)
  * (46) Kumar, H., Mishra, S.: Entropy Stable Numerical Schemes for Two-Fluid Plasma Equations.  Journal of Scientific Computing 52(2), 401–425 (2012)
  * (47) Landau, L., Lifshitz, E.: Relativistic fluid dynamics (1987)
  * (48) Lefloch, P.G., Mercier, J.M., Rohde, C.: Fully discrete, entropy conservative schemes of arbitraryorder.  SIAM Journal on Numerical Analysis 40(5), 1968–1992 (2002)
  * (49) Meena, A.K., Kumar, H.: Robust MUSCL schemes for ten-moment Gaussian closure equations with source terms.  International Journal on Finite Volumes 14, 34 (2017)
  * (50) Meena, A.K., Kumar, H.: A well-balanced scheme for Ten-Moment Gaussian closure equations with source term.  Zeitschrift für angewandte Mathematik und Physik 69(1), 8 (2017)
  * (51) Meena, A.K., Kumar, H.: Robust numerical schemes for two-fluid ten-moment plasma flow equations.  Zeitschrift für angewandte Mathematik und Physik 70, 1–30 (2019)
  * (52) Meena, A.K., Kumar, H., Chandrashekar, P.: Positivity-preserving high-order discontinuous Galerkin schemes for Ten-Moment Gaussian closure equations.  Journal of Computational Physics 339, 370–395 (2017)
  * (53) Meng, X., Tóth, G., Sokolov, I.V., Gombosi, T.I.: Classical and semirelativistic magnetohydrodynamics with anisotropic ion pressure.  Journal of Computational Physics 231(9), 3610–3622 (2012)
  * (54) Orszag, S.A., Tang, C.M.: Small-scale structure of two-dimensional magnetohydrodynamic turbulence.  Journal of Fluid Mechanics 90(1), 129–143 (1979)
  * (55) Otto, A.: 3d resistive mhd computations of magnetospheric physics.  Computer Physics Communications 59(1), 185–195 (1990)
  * (56) Pareschi, L., Russo, G.: Implicit–explicit runge–kutta schemes and applications to hyperbolic systems with relaxation.  Journal of Scientific computing 25, 129–155 (2005)
  * (57) Roe, P.L.: Affordable, entropy consistent flux functions.  In: Eleventh International Conference on Hyperbolic Problems: Theory, Numerics and Applications, Lyon (2006)
  * (58) Ryu, D., Jones, T.W.: Numerical magetohydrodynamics in astrophysics: Algorithm and tests for one-dimensional flow`.  The Astrophysical Journal 442, 228 (1995)
  * (59) Sangam, A.: An hllc scheme for ten-moments approximation coupled with magnetic field.  International Journal of Computing Science and Mathematics 2(1-2), 73–109 (2008)
  * (60) Sen, C., Kumar, H.: Entropy Stable Schemes For Ten-Moment Gaussian Closure Equations.  Journal of Scientific Computing 75(2), 1128–1155 (2018)
  * (61) Tadmor, E.: Entropy stability theory for difference approximations of nonlinear conservation laws and related time-dependent problems.  Acta Numerica 12, 451–512 (2003)
  * (62) Tóth, G.: The div b= 0 constraint in shock-capturing magnetohydrodynamics codes.  Journal of Computational Physics 161(2), 605–652 (2000)
  * (63) Wang, L., Hakim, A.H., Ng, J., Dong, C., Germaschewski, K.: Exact and locally implicit source term solvers for multifluid-Maxwell systems.  Journal of Computational Physics 415, 109510 (2020)
  * (64) Yadav, A., Bhoriya, D., Kumar, H., Chandrashekar, P.: Entropy stable schemes for the shear shallow water model equations.  Journal of Scientific Computing 97(3), 77 (2023)
  * (65) Zhao, X., Yang, Y., Seyler, C.E.: A positivity-preserving semi-implicit discontinuous Galerkin scheme for solving extended magnetohydrodynamics equations.  Journal of Computational Physics 278, 400–415 (2014)
  * (66) Zouganelis, I., Maksimovic, M., Meyer-Vernet, N., Lamy, H., Issautier, K.: A transonic collisionless model of the solar wind.  The Astrophysical Journal 606(1), 542 (2004)

##  Appendix A Entropy scaled right eigenvectors for conservative part of the
reformulated CGL system

In this section, we present expressions for the right eigenvectors and scaled
right eigenvectors for the conservative system

$\frac{\partial\textbf{U}}{\partial
t}+\mathcal{A}_{x}(\textbf{U})\frac{\partial\textbf{U}}{\partial
x}=0,~{}~{}~{}~{}~{}~{}~{}~{}\text{where,}~{}\mathcal{A}_{x}(\textbf{U})=\frac{\partial\textbf{f}_{x}}{\partial\textbf{U}}+\phi^{\prime}(\textbf{V})^{\top}B_{x}^{\prime}(\textbf{U})$
(58)

The right eigenvectors are described in Section A.1, while the expressions of
scaled eigenvectors have been demonstrated in Section A.2. We will do all the
analysis in the term of primitive variable
$\textbf{W}=\\{\rho,u_{x},u_{y},u_{z},p_{\parallel},p_{\perp},B_{x},B_{y},B_{z}\\}$.

###  A.1 Eigenvalues and Right eigenvectors

To derive the eigenvalues and right eigenvectors, it is useful to transform
the system (58) in terms of the primitive variables W. The set of eigenvalues
$\Lambda_{1}$ of the jacobian matrix $\mathcal{A}_{x}$ are given as

$\displaystyle\lambda_{1,2,3,4,5,6,7,8}=u_{x},~{}u_{x},~{}u_{x}\pm
v_{ax},~{}u_{x}\pm c_{f},~{}u_{x}\pm c_{s},$

where,

$\displaystyle
v^{2}_{ax}=\frac{B_{x}^{2}}{\rho},~{}v^{2}_{a}=\frac{B^{2}}{\rho},~{}a^{2}=\frac{2P_{\perp}}{\rho},$
$\displaystyle
c^{2}_{f,s}=\frac{1}{2}\bigg{[}(v_{a}^{2}+a^{2})\pm\sqrt{(v_{a}^{2}+a^{2})^{2}-4v^{2}_{ax}a^{2}}\bigg{]}.$

The corresponding right eigenvectors are

$\displaystyle R_{\lambda_{1}}=\begin{pmatrix}1\\\ 0\\\ 0\\\ 0\\\ 0\\\ 0\\\
0\\\ 0\\\ \end{pmatrix},R_{\lambda_{2}}=\begin{pmatrix}0\\\ 0\\\ 0\\\ 0\\\
1\\\ 0\\\ 0\\\ 0\end{pmatrix},R_{\lambda_{3,4}}=\begin{pmatrix}0\\\ 0\\\
\pm{\beta_{z}}\\\ \mp{\beta_{y}}\\\ 0\\\ 0\\\
-{\beta_{z}}sgn(B_{x})\sqrt{\rho}\\\ {\beta_{y}}sgn(B_{x})\sqrt{\rho}\\\
\end{pmatrix},$ $\displaystyle
R_{\lambda_{5,6}}=\begin{pmatrix}\alpha_{f}\rho\\\ \pm\alpha_{f}c_{f}\\\
\mp\alpha_{s}c_{s}\beta_{y}sgn(B_{x})\\\
\mp\alpha_{s}c_{s}\beta_{z}sgn(B_{x})\\\ \alpha_{f}p_{\parallel}\\\
\alpha_{f}\rho a^{2}\\\ \alpha_{s}a\beta_{y}\sqrt{\rho}\\\
\alpha_{s}a\beta_{z}\sqrt{\rho}\end{pmatrix},R_{\lambda_{7,8}}=\begin{pmatrix}\alpha_{s}\rho\\\
\pm\alpha_{s}c_{s}\\\ \pm\alpha_{f}c_{f}\beta_{y}sgn(B_{x})\\\
\pm\alpha_{f}c_{f}\beta_{z}sgn(B_{x})\\\ \alpha_{s}p_{\parallel}\\\
\alpha_{s}\rho a^{2}\\\ -\alpha_{f}a\beta_{y}\sqrt{\rho}\\\
-\alpha_{f}a\beta_{z}\sqrt{\rho}\end{pmatrix},$

where,

$\displaystyle\alpha_{f}^{2}=\frac{a^{2}-c_{s}^{2}}{c_{f}^{2}-c_{s}^{2}},~{}\alpha_{s}^{2}=\frac{c_{f}^{2}-a^{2}}{c_{f}^{2}-c_{s}^{2}},~{}\beta_{y}=\frac{B_{y}}{\sqrt{B_{y}^{2}+B_{z}^{2}}},~{}\beta_{z}=\frac{B_{z}}{\sqrt{B_{y}^{2}+B_{z}^{2}}}.$

Above all, eigenvectors are linearly independent. The above-defined right
eigenvectors are singular in a variety of cases as listed by Balsara and Roe
in balsara1999eigenstructure . Similarly, in $y$-direction, the set of
eigenvalues $\Lambda_{2}$ of the jacobian matrix
$\frac{\partial\textbf{f}_{y}}{\partial\textbf{U}}$ are given as

$\displaystyle\lambda_{1,2,3,4,5,6,7,8}=u_{y},~{}u_{y},~{}u_{y}\pm
v_{ay},~{}u_{y}\pm c_{f},~{}u_{x}\pm c_{s},$

where,

$\displaystyle
v^{2}_{ay}=\frac{B_{y}^{2}}{\rho},~{}v^{2}_{a}=\frac{B^{2}}{\rho},~{}a^{2}=\frac{2P_{\perp}}{\rho},$
$\displaystyle
c^{2}_{f,s}=\frac{1}{2}\bigg{[}(v_{a}^{2}+a^{2})\pm\sqrt{(v_{a}^{2}+a^{2})^{2}-4v^{2}_{ay}a^{2}}\bigg{]}.$

the corresponding right eigenvectors are

$\displaystyle R_{\lambda_{1}}=\begin{pmatrix}1\\\ 0\\\ 0\\\ 0\\\ 0\\\ 0\\\
0\\\ 0\\\ \end{pmatrix},R_{\lambda_{2}}=\begin{pmatrix}0\\\ 0\\\ 0\\\ 0\\\
1\\\ 0\\\ 0\\\ 0\end{pmatrix},R_{\lambda_{3,4}}=\begin{pmatrix}0\\\
\pm{\beta_{z}}\\\ 0\\\ \mp{\beta_{x}}\\\ 0\\\ 0\\\
-{\beta_{z}}sgn(B_{y})\sqrt{\rho}\\\ {\beta_{x}}sgn(B_{y})\sqrt{\rho}\\\
\end{pmatrix},$ $\displaystyle
R_{\lambda_{5,6}}=\begin{pmatrix}\alpha_{f}\rho\\\
\mp\alpha_{s}c_{s}\beta_{x}sgn(B_{y})\\\ \pm\alpha_{f}c_{f}\\\
\mp\alpha_{s}c_{s}\beta_{z}sgn(B_{y})\\\ \alpha_{f}p_{\parallel}\\\
\alpha_{f}\rho a^{2}\\\ \alpha_{s}a\beta_{x}\sqrt{\rho}\\\
\alpha_{s}a\beta_{z}\sqrt{\rho}\end{pmatrix},R_{\lambda_{7,8}}=\begin{pmatrix}\alpha_{s}\rho\\\
\pm\alpha_{f}c_{f}\beta_{x}sgn(B_{y})\\\ \pm\alpha_{s}c_{s}\\\
\pm\alpha_{f}c_{f}\beta_{z}sgn(B_{y})\\\ \alpha_{s}p_{\parallel}\\\
\alpha_{s}\rho a^{2}\\\ -\alpha_{f}a\beta_{x}\sqrt{\rho}\\\
-\alpha_{f}a\beta_{z}\sqrt{\rho}\end{pmatrix},$

where,

$\displaystyle\alpha_{f}^{2}=\frac{a^{2}-c_{s}^{2}}{c_{f}^{2}-c_{s}^{2}},~{}\alpha_{s}^{2}=\frac{c_{f}^{2}-a^{2}}{c_{f}^{2}-c_{s}^{2}},~{}\beta_{x}=\frac{B_{x}}{\sqrt{B_{x}^{2}+B_{z}^{2}}},~{}\beta_{z}=\frac{B_{z}}{\sqrt{B_{x}^{2}+B_{z}^{2}}}.$

###  A.2 Entropy scaled right eigenvectors

In this section, we have calculated the entropy scaled right eigenvectors for
the case of $x$-direction. Consider the conservative part of the CGL system
(58). The eigenvalues and Right eigenvectors in terms of primitive of the
Jacobian matrix $\mathcal{A}_{1}$ are described in (A.1). The right
eigenvector matrix $R^{x}$ for the matrix $\mathcal{A}_{1}$ in terms of
conservative variable is given by the relation

$R^{x}=\frac{\partial\textbf{U}}{\partial\textbf{W}}R_{\textbf{W}}^{x}$

where, $\frac{\partial\textbf{U}}{\partial\textbf{W}}$ is the Jacobian matrix
for the change of variable, given by

$\displaystyle\frac{\partial\textbf{U}}{\partial\textbf{W}}=\begin{pmatrix}1&0&0&0&0&0&0&0&0\\\
u_{x}&\rho&0&0&0&0&0&0&0\\\ u_{y}&0&\rho&0&0&0&0&0&0\\\
u_{z}&0&0&\rho&0&0&0&0&0\\\ 0&0&0&0&1&0&0&0&0\\\ \frac{\textbf{u}^{2}}{2}&\rho
u_{x}&\rho u_{y}&\rho u_{z}&\frac{1}{2}&1&B_{x}&B_{y}&B_{z}\\\
0&0&0&0&0&0&1&0&0\\\ 0&0&0&0&0&0&0&1&0\\\ 0&0&0&0&0&0&0&0&1\end{pmatrix}.$

and the matrix $R_{\textbf{W}}^{x}$ is given by,

$\displaystyle\begin{pmatrix}\alpha_{f}\rho&\alpha_{s}\rho&0&1&0&0&0&\alpha_{s}\rho&\alpha_{f}\rho\\\
-\alpha_{f}c_{f}&-\alpha_{s}c_{s}&0&0&0&0&0&\alpha_{s}c_{s}&\alpha_{f}c_{f}\\\
\alpha_{s}\beta_{y}c_{s}\mathcal{S}_{x}&-\alpha_{f}\beta_{y}c_{f}\mathcal{S}_{x}&-\beta_{z}&0&0&0&\beta_{z}&\alpha_{f}\beta_{y}c_{f}\mathcal{S}_{x}&-\alpha_{s}\beta_{y}c_{s}\mathcal{S}_{x}\\\
\alpha_{s}\beta_{z}c_{s}\mathcal{S}_{x}&-\alpha_{f}\beta_{z}c_{f}\mathcal{S}_{x}&\beta_{y}&0&0&0&-\beta_{y}&\alpha_{f}\beta_{z}c_{f}\mathcal{S}_{x}&-\alpha_{s}\beta_{z}c_{s}\mathcal{S}_{x}\\\
\alpha_{f}p_{\parallel}&\alpha_{s}p_{\parallel}&0&0&0&1&0&\alpha_{s}p_{\parallel}&\alpha_{f}p_{\parallel}\\\
a^{2}\alpha_{f}\rho&a^{2}\alpha_{s}\rho&0&0&0&0&0&a^{2}\alpha_{s}\rho&a^{2}\alpha_{f}\rho\\\
0&0&0&0&1&0&0&0&0\\\
a\alpha_{s}\beta_{y}\sqrt{\rho}&-a\alpha_{f}\beta_{y}\sqrt{\rho}&-\beta_{z}\sqrt{\rho}\mathcal{S}_{x}&0&0&0&-\beta_{z}\sqrt{\rho}\mathcal{S}_{x}&-a\alpha_{f}\beta_{y}\sqrt{\rho}&a\alpha_{s}\beta_{y}\sqrt{\rho}\\\
a\alpha_{s}\beta_{z}\sqrt{\rho}&-a\alpha_{f}\beta_{z}\sqrt{\rho}&\beta_{y}\sqrt{\rho}\mathcal{S}_{x}&0&0&0&\beta_{y}\sqrt{\rho}\mathcal{S}_{x}&-a\alpha_{f}\beta_{z}\sqrt{\rho}&a\alpha_{s}\beta_{z}\sqrt{\rho}\\\
\end{pmatrix}.$

where, $\mathcal{S}_{x}=sgn(B_{x})$. We need to find a scaling matrix $T^{x}$
such that the scaled right eigenvector matrix $\tilde{R}^{x}=R^{x}T^{x}$
satisfies

$\frac{\partial\textbf{U}}{\partial\textbf{V}}=\tilde{R}^{x}\tilde{R}^{x\top}$
(59)

where V is the entropy variable vector as in Eqn. (12). We follow the Barth
scaling process barth1999numerical  to scale the right eigenvectors. The
scaling matrix $T^{x}$ is the square root of $Y^{x}$, where $Y^{x}$ has the
expression

$Y^{x}=(R_{\textbf{W}}^{x})^{-1}\frac{\partial\textbf{W}}{\partial\textbf{V}}\left(\frac{\partial\textbf{U}}{\partial\textbf{W}}\right)^{-\top}(R_{\textbf{W}}^{x})^{-\top}$

which results in

$\displaystyle Y^{x}=\begin{pmatrix}\frac{1}{8\rho}&0&0&0&0&0&0&0&0\\\
0&\frac{1}{8\rho}&0&0&0&0&0&0&0\\\
0&0&\frac{p_{\perp}}{4\rho^{2}}&0&0&0&0&0&0\\\
0&0&0&\frac{\rho}{4}&0&\frac{p_{\parallel}}{4}&0&0&0\\\
0&0&0&0&\frac{p_{\perp}}{2\rho}&0&0&0&0\\\
0&0&0&\frac{p_{\parallel}}{4}&0&\frac{5p_{\parallel}^{2}}{4\rho}&0&0&0\\\
0&0&0&0&0&0&\frac{p_{\perp}}{4\rho^{2}}&0&0\\\
0&0&0&0&0&0&0&\frac{1}{8\rho}&0\\\
0&0&0&0&0&0&0&0&\frac{1}{8\rho}\end{pmatrix}.$

Then matrix $T^{x}$ is,

$\displaystyle\begin{pmatrix}\frac{1}{2\sqrt{2\rho}}&0&0&0&0&0&0&0&0\\\
0&\frac{1}{2\sqrt{2\rho}}&0&0&0&0&0&0&0\\\
0&0&\frac{\sqrt{p_{\perp}}}{2\rho}&0&0&0&0&0&0\\\
0&0&0&\frac{\sqrt{\rho}(2p_{\parallel}+\rho)}{2\sqrt{5p_{\parallel}^{2}+4p_{\parallel}\rho+\rho^{2}}}&0&\frac{p_{\parallel}\sqrt{\rho}}{2\sqrt{5p_{\parallel}^{2}+4p_{\parallel}\rho+\rho^{2}}}&0&0&0\\\
0&0&0&0&\sqrt{\frac{p_{\perp}}{2\rho}}&0&0&0&0\\\
0&0&0&\frac{p_{\parallel}\sqrt{\rho}}{2\sqrt{5p_{\parallel}^{2}+4p_{\parallel}\rho+\rho^{2}}}&0&\frac{p_{\parallel}(5p_{\parallel}+2\rho)}{2\sqrt{\rho(5p_{\parallel}^{2}+4p_{\parallel}\rho+\rho^{2})}}&0&0&0\\\
0&0&0&0&0&0&\frac{\sqrt{p_{\perp}}}{2\rho}&0&0\\\
0&0&0&0&0&0&0&\frac{1}{2\sqrt{2\rho}}&0\\\
0&0&0&0&0&0&0&0&\frac{1}{2\sqrt{2\rho}}\end{pmatrix}.$

###### Remark 3

For finding the square root of $3\times 3$ matrix

$\displaystyle M=\begin{pmatrix}a&0&b\\\ 0&e&0\\\ c&0&d\end{pmatrix}$

where $a~{},b~{},c$, $d$ and $e$ are real or complex numbers. we apply the
formula,

$\displaystyle R=\begin{pmatrix}\frac{a+s}{t}&0&\frac{b}{t}\\\ 0&\sqrt{e}&0\\\
\frac{c}{t}&0&\frac{d+s}{t}\end{pmatrix}.$

where, $\delta=ad-bc$, $\tau=a+d$ and $s=\sqrt{\delta}$, $t=\sqrt{\tau+2s}$.
We apply this formula if $t\neq 0$.

###### Remark 4

We proceed similarly in the $y$-direction, the eigenvalues and Right
eigenvectors in terms of primitive of the Jacobian matrix
$\frac{\partial\textbf{f}_{2}}{\partial\textbf{U}}$ describe in (A.1). The
right eigenvector matrix is given by the relation

$R^{y}=\frac{\partial\textbf{U}}{\partial\textbf{W}}R_{\textbf{W}}^{y}$

where, $R_{\textbf{W}}^{y}$ is given below

$\displaystyle\begin{pmatrix}\alpha_{f}\rho&\alpha_{s}\rho&0&1&0&0&0&\alpha_{s}\rho&\alpha_{f}\rho\\\
\alpha_{s}\beta_{x}c_{s}\mathcal{S}_{y}&-\alpha_{f}\beta_{x}c_{f}\mathcal{S}_{y}&-\beta_{z}&0&0&0&\beta_{z}&\alpha_{f}\beta_{x}c_{f}\mathcal{S}_{y}&-\alpha_{s}\beta_{x}c_{s}\mathcal{S}_{y}\\\
-\alpha_{f}c_{f}&-\alpha_{s}c_{s}&0&0&0&0&0&\alpha_{s}c_{s}&\alpha_{f}c_{f}\\\
\alpha_{s}\beta_{z}c_{s}\mathcal{S}_{y}&-\alpha_{f}\beta_{z}c_{f}\mathcal{S}_{y}&\beta_{x}&0&0&0&-\beta_{x}&\alpha_{f}\beta_{z}c_{f}\mathcal{S}_{y}&-\alpha_{s}\beta_{z}c_{s}\mathcal{S}_{y}\\\
\alpha_{f}p_{\parallel}&\alpha_{s}p_{\parallel}&0&0&0&1&0&\alpha_{s}p_{\parallel}&\alpha_{f}p_{\parallel}\\\
a^{2}\alpha_{f}\rho&a^{2}\alpha_{s}\rho&0&0&0&0&0&a^{2}\alpha_{s}\rho&a^{2}\alpha_{f}\rho\\\
a\alpha_{s}\beta_{x}\sqrt{\rho}&-a\alpha_{f}\beta_{x}\sqrt{\rho}&-\beta_{z}\sqrt{\rho}\mathcal{S}_{y}&0&0&0&-\beta_{z}\sqrt{\rho}\mathcal{S}_{y}&-a\alpha_{f}\beta_{x}\sqrt{\rho}&a\alpha_{s}\beta_{x}\sqrt{\rho}\\\
0&0&0&0&1&0&0&0&0\\\
a\alpha_{s}\beta_{z}\sqrt{\rho}&-a\alpha_{f}\beta_{z}\sqrt{\rho}&\beta_{x}\sqrt{\rho}\mathcal{S}_{y}&0&0&0&\beta_{x}\sqrt{\rho}\mathcal{S}_{y}&-a\alpha_{f}\beta_{z}\sqrt{\rho}&a\alpha_{s}\beta_{z}\sqrt{\rho}\\\
\end{pmatrix}$

where, $\mathcal{S}_{y}=sgn(B_{y})$. Accordingly, we obtain the scaling matrix
$T^{y}$ as,

$\displaystyle\begin{pmatrix}\frac{1}{2\sqrt{2\rho}}&0&0&0&0&0&0&0&0\\\
0&\frac{1}{2\sqrt{2\rho}}&0&0&0&0&0&0&0\\\
0&0&\frac{\sqrt{p_{\perp}}}{2\rho}&0&0&0&0&0&0\\\
0&0&0&\frac{\sqrt{\rho}(2p_{\parallel}+\rho)}{2\sqrt{5p_{\parallel}^{2}+4p_{\parallel}\rho+\rho^{2}}}&0&\frac{p_{\parallel}\sqrt{\rho}}{2\sqrt{5p_{\parallel}^{2}+4p_{\parallel}\rho+\rho^{2}}}&0&0&0\\\
0&0&0&0&\sqrt{\frac{p_{\perp}}{2\rho}}&0&0&0&0\\\
0&0&0&\frac{p_{\parallel}\sqrt{\rho}}{2\sqrt{5p_{\parallel}^{2}+4p_{\parallel}\rho+\rho^{2}}}&0&\frac{p_{\parallel}(5p_{\parallel}+2\rho)}{2\sqrt{\rho(5p_{\parallel}^{2}+4p_{\parallel}\rho+\rho^{2})}}&0&0&0\\\
0&0&0&0&0&0&\frac{\sqrt{p_{\perp}}}{2\rho}&0&0\\\
0&0&0&0&0&0&0&\frac{1}{2\sqrt{2\rho}}&0\\\
0&0&0&0&0&0&0&0&\frac{1}{2\sqrt{2\rho}}\end{pmatrix}.$

###  A.3 Entropy scaled right eigenvectors in $x$ and $y$ direction