# Monolithic parabolic regularization of the MHD equations and entropy
principles111This research is funded by Swedish Research Council (VR) under
grant number 2021-04620.

Tuan Anh Dao<EMAIL_ADDRESS>Murtazo Nazarov<EMAIL_ADDRESS>Department of Information Technology, Uppsala University, Sweden

###### Abstract

We show at the PDE level that the monolithic parabolic regularization of the
equations of ideal magnetohydrodynamics (MHD) is compatible with all the
generalized entropies, fulfills the minimum entropy principle, and preserves
the positivity of density and internal energy. We then numerically investigate
this regularization for the MHD equations using continuous finite elements in
space and explicit strong stability preserving Runge-Kuta methods in time. The
artificial viscosity coefficient of the regularization term is constructed to
be proportional to the entropy residual of MHD. It is shown that the method
has a high order of accuracy for smooth problems and captures strong shocks
and discontinuities accurately for non-smooth problems.

###### keywords:

MHD , artificial viscosity , entropy inequalities , viscous regularization ,
entropy viscosity

##  1 Introduction

Designing numerical methods which produce physically relevant solutions has
been an interesting yet challenging task. In the study of non-conducting
fluids, for example such described by the compressible Euler equations,
entropy principles, and positivity preserving properties have been studied for
many schemes. Those include Godunov scheme [9], Lax scheme [22], and its
variants, see e.g., [28, 34, 32]. Positivity preserving properties often refer
to the fact that density, internal energy, or pressure of the solution should
maintain positive as it evolves in time.  These properties are also required
when solving the MHD equations. The ideal MHD equations describe the coupling
of hydrodynamic motion of plasma fluid with very little resistivity and its
electromagnetic field.  Existing works on the compressible Euler equations are
difficult to extend to the MHD system. Several reasons are: ($i$) the
additional consideration of electromagnetism leads to complex hyperbolicity of
the MHD system; ($ii$) the MHD flux is non-convex; ($iii$) the solenoidal
nature of the magnetic field is numerically important to maintain but it is
not explicitly described by the MHD equations. A common approach to achieving
positivity and entropy principles for hyperbolic systems consists of adding a
vanishing viscous regularization, see e.g., [11, 21]. The viscous
regularization can be chosen such that the resulting system matches the
resistive model of the MHD equations, see e.g., [2, 6], or simply chosen as
Laplacian terms under equal weights to all the conserved quantities, here we
refer to as the “monolithic parabolic regularization”. The resistive MHD
viscous flux suffers from several drawbacks for numerical purposes. One
difficulty is that there is no regularization to the mass equation, which
makes it incompatible with most numerical methods due to the Gibbs phenomenon.
Another problem is that unless the thermal diffusivity is zero, the resistive
MHD flux violates the minimum entropy principle, see [11]. On the other hand,
the non-physically-motivated monolithic viscous flux is employed in numerical
schemes for MHD, e.g., [20, 24] as well as being a continuous analog of the
well-known Lax-Friedrichs or upwind schemes. However, to the best of our
knowledge, investigations regarding entropy principles and other positivity-
preserving properties of viscous regularizations to the MHD equations are
still missing in the literature.

The main focus and contribution of this paper is to investigate the entropy
principles of the monolithic parabolic regularization to the ideal MHD
equations at the PDE level. We prove that the conservative form of the MHD
equations with monolithic parabolic regularization satisfies positivity of
density, minimum entropy principle, positivity of internal energy, and is
compatible with all the generalized entropy inequalities in the manner of
[17]. Our analysis is an extension of the works on compressible Euler
equations of [11] and [17]. The theory encourages the use of the monolithic
viscous flux as a stabilization tool for the numerical approximation of the
ideal MHD equations.

State-of-the-art finite element methods for solving convection-dominated
problems, such as the compressible Euler the MHD equations, are mainly based
on the least-squares argument, see e.g., [18] and references therein. A
secondary contribution of this article is the use of the monolithic parabolic
regularization developed in this work as a finite element stabilization term
without invoking any least-squares argument. First, in Section 4.5 we validate
the developed theory using a first order artificial viscosity and compare the
monolithic viscous flux with the traditional Navier-Stokes resistive viscous
flux. The numerical results reveal several benefits of the monolithic flux
over the resistive MHD flux, which fit well with the theoretical findings.

Then, we introduce the entropy viscosity method for MHD in the spirit of [10,
13, 26], where the regularization coefficients are constructed to be
proportional to the entropy residual of the MHD system. Although the magnetic
field was divergenceless in the continuous case, this is not the case in the
discrete approximation. It is necessary to apply some divergence cleaning
algorithms. In this article, we use the projection method [3] to fix the
divergence error. The tests show that the resulting method can maintain high-
order accuracy of the smooth solutions, capture the shocks accurately, and
remain stable even in the turbulence phase of the numerical solution.

The rest of the paper is organized as follows. In Section 2, we describe the
ideal MHD equations. Section 3 contains the analysis of the monolithic
parabolic regularization at the PDE level. In Section 4, we relate the
theoretical results in the continuous level with numerical results using a CG
discretization. Section 5 contains a summary of our contributions and
concluding remarks.

##  2 The ideal MHD equations

Consider a $d$-dimensional spatial domain
${\boldsymbol{x}}\in{\mathbb{R}}^{d}$, $d\in{\mathbb{N}}^{*}$, and a temporal
domain $t\in[0,+\infty)$. We define the space time domain
$D:={\mathbb{R}}^{d}\times[0,+\infty)$. Let us denote
${\boldsymbol{U}}:=(\rho,{\boldsymbol{m}},E,{\boldsymbol{B}})$ a vector
containing several conserved quantities of a conducting fluid: density
$\rho({\boldsymbol{x}},t):D\to{\mathbb{R}}$, momentum
${\boldsymbol{m}}({\boldsymbol{x}},t):D\to{\mathbb{R}}^{d}$, total energy
$E({\boldsymbol{x}},t):D\to{\mathbb{R}}$, and magnetic field
${\boldsymbol{B}}({\boldsymbol{x}},t):D\to{\mathbb{R}}^{d}$. The velocity
field ${\boldsymbol{u}}({\boldsymbol{x}},t)$ is determined by the relation
${\boldsymbol{u}}({\boldsymbol{x}},t):={\boldsymbol{m}}({\boldsymbol{x}},t)/\rho({\boldsymbol{x}},t)$.
The governing electrodynamics and fluid mechanics of such fluid can be
described in the following ideal MHD equations,

${\partial}_{t}\textsf{{U}}+\nabla{\cdot\,}\textsf{{F}}_{{\mathcal{E}}}(\textsf{{U}})+\nabla{\cdot\,}\textsf{{F}}_{{\mathcal{B}}}(\textsf{{U}})=0,$
(2.1)

where the nonlinear tensor fluxes $\textsf{{F}}_{{\mathcal{E}}}(\textsf{{U}})$
and $\textsf{{F}}_{{\mathcal{B}}}(\textsf{{U}})$ are defined as

$\textsf{{F}}_{{\mathcal{E}}}(\textsf{{U}}):=\begin{pmatrix}{\boldsymbol{m}}\\\
{\boldsymbol{m}}\otimes{\boldsymbol{u}}+p{\mathbb{I}}\\\
{\boldsymbol{u}}(E+p)\\\
0\end{pmatrix},\quad\textsf{{F}}_{{\mathcal{B}}}(\textsf{{U}}):=\begin{pmatrix}0\\\
-{\boldsymbol{\beta}}\\\ -{\boldsymbol{\beta}}{\cdot}{\boldsymbol{u}}\\\
{\boldsymbol{u}}\otimes{\boldsymbol{B}}-{\boldsymbol{B}}\otimes{\boldsymbol{u}}\\\
\end{pmatrix},$

${\mathbb{I}}$ denotes the identity matrix of size $d\times d$, and $p$ is the
thermodynamic pressure. The symmetric term ${\boldsymbol{\beta}}$ is called
the Maxwell stress tensor:

${\boldsymbol{\beta}}:=-\frac{1}{2}({\boldsymbol{B}}{\cdot}{\boldsymbol{B}}){\mathbb{I}}+{\boldsymbol{B}}\otimes{\boldsymbol{B}}.$

In the absence of the magnetic field or the conductivity of the fluid
$\nabla{\cdot\,}\textsf{{F}}_{{\mathcal{B}}}(\textsf{{U}})\equiv 0$, the MHD
equations (2.1) reduce to the compressible Euler equations. Thus, in certain
circumstances, it is useful to extend the knowledge of the Euler equations to
study the MHD equations. The solenoidal nature of the magnetic field
${\boldsymbol{B}}$ is implied from Faraday’s law,

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

which is also known as the “divergence free” constraint. In the rest of the
paper, we use the fact that (2.2) holds in the continuous settings. The
specific internal energy $e$ is defined as

$\rho
e:=E-\frac{1}{2}\rho{\boldsymbol{u}}^{2}-\frac{1}{2}{\boldsymbol{B}}^{2},$
(2.3)

where the term $\frac{1}{2}{\boldsymbol{B}}^{2}$ is often referred to as the
magnetic pressure.

The specific entropy $s(\rho,e)$ and the thermodynamic pressure $p$ are
defined through the following thermodynamic identity, see [31, Section
12.2.4],

$\,\mathrm{d}e:=T\,\mathrm{d}s+\frac{p}{\rho^{2}}\,\mathrm{d}\rho.$ (2.4)

Denote $s_{e}:=\frac{{\partial}s}{{\partial}e}$ and
$s_{\rho}:=\frac{{\partial}s}{{\partial}\rho}$. The following simple chain
rule will be frequently used,

${\partial}_{\alpha}s=s_{e}{\partial}_{\alpha}e+s_{\rho}{\partial}_{\alpha}\rho,$
(2.5)

where $\alpha=\\{{\boldsymbol{x}},t\\}$. The equation (2.4) can be written as

$\,\mathrm{d}s=T^{-1}\,\mathrm{d}e-pT^{-1}\rho^{-2}\,\mathrm{d}\rho.$ (2.6)

Combining (2.6) with the following direct consequence of (2.5),

$\,\mathrm{d}s=\frac{{\partial}s}{{\partial}e}\,\mathrm{d}e+\frac{{\partial}s}{{\partial}\rho}\,\mathrm{d}\rho,$

we deduce that

$s_{e}:=T^{-1},\quad s_{\rho}:=-pT^{-1}\rho^{-2},$

with the second equality being called the “equation of state” and being
equivalent to

$ps_{e}+\rho^{2}s_{\rho}=0.$ (2.7)

Throughout this paper, we have the following assumptions: the temperature $T$
is positive, which is equivalent to

$s_{e}>0,$

and the specific entropy $-s$ is strictly convex with respect to $\rho^{-1}$
and $e$. Note that similar to [11], we do not make the assumption on the sign
of the pressure $p$, as it is allowed to be both positive and negative. The
strict convexity of $-s$ implies that, see Appendix A,

${\partial}_{\rho}(\rho^{2}s_{\rho})<0,\quad
s_{ee}<0,\quad{\partial}_{\rho}(\rho^{2}s_{\rho})s_{ee}-\rho^{2}s_{\rho
e}^{2}>0.$ (2.8)

For ideal gases, it is common to use

$p=(\gamma-1)\rho e=\rho T$ (2.9)

to compute pressure and temperature from the conserved variables, where
$\gamma>1$ is the adiabatic gas constant and

$s=c_{v}\ln\frac{p}{\rho^{\gamma}},$ (2.10)

where $c_{v}$ is the specific heat capacity at constant volume. We note that
the main results of this work, which are derived in Section 3, do not rely on
assuming that the gas is ideal.

##  3 Monolithic parabolic regularization

The ideal MHD equations can be regularized with monolithic parabolic terms as
in the conservative form below,

${\partial}_{t}\textsf{{U}}+\nabla{\cdot\,}\textsf{{F}}_{{\mathcal{E}}}(\textsf{{U}})+\nabla{\cdot\,}\textsf{{F}}_{{\mathcal{B}}}(\textsf{{U}})-\nabla{\cdot\,}\textsf{{F}}_{{\mathcal{V}}}^{m}(\textsf{{U}})=0,$
(3.1)

where

$\textsf{{F}}_{{\mathcal{V}}}^{m}(\textsf{{U}}):=\begin{pmatrix}\epsilon\nabla\rho\\\
\epsilon\nabla{\boldsymbol{m}}\\\ \epsilon\nabla E\\\
\epsilon\nabla{\boldsymbol{B}}\end{pmatrix},$ (3.2)

and $\epsilon$ is a small positive constant. The viscous flux (3.2) is called
“monolithic” because it simply adds the Laplacian terms corresponding to the
conserved variables under an equal weight $\epsilon$. For reference, we
separately write the equations in (3.1) as

$\displaystyle{\partial}_{t}\rho$
$\displaystyle+\nabla{\cdot\,}{\boldsymbol{m}}$
$\displaystyle=\epsilon\Delta\rho$ (3.3)
$\displaystyle{\partial}_{t}{\boldsymbol{m}}$
$\displaystyle+\nabla{\cdot\,}({\boldsymbol{m}}\otimes{\boldsymbol{u}}+p{\mathbb{I}})$
$\displaystyle-\nabla{\cdot\,}{\boldsymbol{\beta}}$
$\displaystyle=\epsilon\Delta{\boldsymbol{m}}$ (3.4)
$\displaystyle{\partial}_{t}E$
$\displaystyle+\nabla{\cdot\,}({\boldsymbol{u}}(E+p))$
$\displaystyle-\nabla{\cdot\,}({\boldsymbol{u}}{\cdot}{\boldsymbol{\beta}})$
$\displaystyle=\epsilon\Delta E$ (3.5)
$\displaystyle{\partial}_{t}{\boldsymbol{B}}$
$\displaystyle+\nabla{\cdot\,}({\boldsymbol{u}}\otimes{\boldsymbol{B}}-{\boldsymbol{B}}\otimes{\boldsymbol{u}})$
$\displaystyle=\epsilon\Delta{\boldsymbol{B}}$ (3.6)

For numerical purposes, the constant $\epsilon$ is often dependent on the
resolution scale of the underlying scheme. The monolithic parabolic viscous
terms can be linked to many classical numerical schemes with attractive
properties. For example, on a uniform grid in one spatial dimension
$\\{\dots,x_{i-1},x_{i},x_{i+1},\dots\\}$, the well-known Lax-Friedrichs
method applied on the unregularized mass conservation equation

${\partial}_{t}\rho+\nabla{\cdot\,}{\boldsymbol{m}}=0$ (3.7)

reads

$\rho_{i}^{n+1}=\frac{1}{2}(\rho_{i-1}^{n}+\rho_{i+1}^{n})-\frac{\delta
t}{2\delta x}\left(m_{i+1}^{n}-m_{i-1}^{n}\right),$ (3.8)

where $\rho_{i}^{n}$ approximates $\rho(x_{i},t_{n})$, $m_{i}^{n}$
approximates ${\boldsymbol{m}}(x_{i},t_{n})$, $x_{i}\in{\mathbb{R}}$,
$t_{n+1}>t_{n}$, $\delta x=x_{i+1}-x_{i}>0$, and $\delta t>0$ is the time
step. Assuming sufficient smoothness, it is known that (3.8) is stable upon
choosing $\delta t=C_{LF}\delta x$, where $C_{LF}$ is a sufficiently small
real number. By adding and subtracting $\rho_{i}^{n}$ from the right-hand-side
of (3.8), we can rewrite (3.8) as

$\rho_{i}^{n+1}=\rho_{i}^{n}-\frac{\delta t}{2\delta
x}\left(m_{i+1}^{n}-m_{i-1}^{n}\right)+\epsilon_{LF}\frac{\delta t}{\delta
x^{2}}(\rho_{i+1}^{n}-2\rho_{i}^{n}+\rho_{i-1}^{n}),$ (3.9)

where $\epsilon_{LF}=\frac{\delta x^{2}}{2\delta t}$. One can see that (3.9)
is a consistent discretization of (3.3) when $\epsilon$ is set to
$\epsilon_{LF}$. If $\delta t$ is chosen to be of order ${\mathcal{O}}(\delta
x)$ following the stability condition, the viscosity coefficient
$\epsilon_{LF}$ is of order ${\mathcal{O}}(\delta x)$ and vanishes as $\delta
x\to\infty$.

Upwind schemes can also be viewed as directly related to (3.3). A regular
upwind scheme approximating the nonlinear equation (3.7) is often known as,

$\rho_{i}^{n+1}=\rho_{i}^{n}-\frac{\delta t}{2\delta
x}\left(m_{i+1}^{n}-m_{i-1}^{n}\right)+\epsilon_{Up}\frac{\delta t}{\delta
x^{2}}(\rho_{i+1}^{n}-2\rho_{i}^{n}+\rho_{i-1}^{n}),$ (3.10)

where the vanishing viscosity coefficient reads
$\epsilon_{Up}=\frac{1}{2}\delta
x\|{\boldsymbol{u}}\|_{L^{\infty}({\mathbb{R}}^{d})}$ with
$\|{\boldsymbol{u}}\|_{L^{\infty}({\mathbb{R}}^{d})}$ being the maximum
density wave speed. Again, it can be seen that (3.10) is a consistent
discretization of (3.3) when $\epsilon$ is set to $\epsilon_{Up}$.

###  3.1 Positivity of density

For positivity of density, one can apply the same proof as [11, Section 3.1]
did for the Euler equations since no equations other than the scalar mass
conservation equation (3.3) are invoked. We conclude that the density to
solution to (3.3) is positive given that the assumptions in the following
theorem hold and the viscous parameter $\epsilon$ in (3.3) is positive.

######  Theorem 3.1 (Positivity of density, [11]).

Upon the following justifiable assumptions:

  1. (i)

${\boldsymbol{u}}$ and $\nabla{\cdot\,}{\boldsymbol{u}}\in L^{1}(D)$;

  2. (ii)

$\int_{0}^{+\infty}({\partial}_{t}\rho+\nabla{\cdot\,}{\boldsymbol{m}})\,\mathrm{d}t$
and $\int_{0}^{+\infty}\epsilon\nabla\rho\,\mathrm{d}t\in
L^{1}({\mathbb{R}}^{d})$;

  3. (iii)

there exists $\rho_{\min}>0$ such that outside of a closed ball in
${\mathbb{R}}^{d}$ centered at the origin with radius
$R_{\widehat{t}}<\infty$, $\rho({\boldsymbol{x}},t)\geq\rho_{\min}$ for all
$t\in(0,\widehat{t}\,),\,\widehat{t}>0$, and
$\int_{0}^{\widehat{t}}(\nabla{\cdot\,}{\boldsymbol{u}}({\boldsymbol{x}},\widehat{t}\,))^{+}\,\mathrm{d}t\in
L^{\infty}({\mathbb{R}}^{d})$,

the density solution to (3.3) satisfies the following positivity property

$\operatorname{essinf}_{{\boldsymbol{x}}\in{\mathbb{R}}^{d}}\rho({\boldsymbol{x}},t)>0,\quad\forall
t>0.$

In simple words, Theorem 3.1 implies that under sufficient smoothness of the
density and the velocity, and that outside of a region of interest, then the
parabolic term $\epsilon\Delta\rho$, for an arbitrary positive value of
$\epsilon$, regularizes (3.3) and guarantees positivity of the density
solution.

###  3.2 Minimum entropy principle

In this section, we investigate the continuous minimum entropy principle of
the monolithic flux. We write (3.3)-(3.6) in the following nonconservative
form,

$\displaystyle{\partial}_{t}\rho+{\boldsymbol{u}}{\cdot}\nabla\rho+\rho\nabla{\cdot\,}{\boldsymbol{u}}$
$\displaystyle=\epsilon\Delta\rho,$ (3.11)
$\displaystyle\rho({\partial}_{t}{\boldsymbol{u}}+{\boldsymbol{u}}{\cdot}\nabla{\boldsymbol{u}})+{\boldsymbol{u}}\epsilon\Delta\rho+\nabla
p$ $\displaystyle-\nabla{\cdot\,}{\boldsymbol{\beta}}$
$\displaystyle=\epsilon\Delta{\boldsymbol{m}},$ (3.12)
$\displaystyle\rho({\partial}_{t}{\mathcal{E}}+{\boldsymbol{u}}{\cdot}\nabla{\mathcal{E}})+{\mathcal{E}}\epsilon\Delta\rho+\nabla{\cdot\,}({\boldsymbol{u}}p)$
$\displaystyle-\nabla{\cdot\,}({\boldsymbol{u}}{\cdot}{\boldsymbol{\beta}})$
$\displaystyle=\epsilon\Delta E,$ (3.13)
$\displaystyle{\partial}_{t}{\boldsymbol{B}}$
$\displaystyle+{\boldsymbol{B}}(\nabla{\cdot\,}{\boldsymbol{u}})+({\boldsymbol{u}}{\cdot}\nabla){\boldsymbol{B}}-({\boldsymbol{B}}{\cdot}\nabla){\boldsymbol{u}}$
$\displaystyle=\epsilon\Delta{\boldsymbol{B}},$ (3.14)

where ${\mathcal{E}}=\rho^{-1}E$.

For convenience, we use the notation ${\boldsymbol{l}}:=\epsilon\nabla(\rho
e)$ which will appear frequently. The term ${\boldsymbol{l}}$ can be written
in another form described by the following lemma.

######  Lemma 1.

The following equality holds

${\boldsymbol{l}}=s_{e}^{-1}\left(es_{e}-\rho
s_{\rho}\right)(\epsilon\nabla\rho)+\epsilon\rho s_{e}^{-1}\nabla s.$ (3.15)

###### Proof.

Using the chain rule (2.5), we obtain $\nabla s=s_{\rho}\nabla\rho+s_{e}\nabla
e$. We then have

$\displaystyle s_{e}\nabla(\rho e)$ $\displaystyle=es_{e}\nabla\rho+\rho
s_{e}\nabla e$ $\displaystyle=es_{e}\nabla\rho+(\rho\nabla s-\rho
s_{\rho}\nabla\rho)$ $\displaystyle=(es_{e}-\rho
s_{\rho})\nabla\rho+\rho\nabla s.$

Since $\epsilon s_{e}^{-1}>0$, multiplying both sides of the above equality
with $\epsilon s_{e}^{-1}$ gives the desired conclusion. ∎

We introduce the following matrices,

$\quad{\mathbb{G}}_{{\boldsymbol{u}}}:=\epsilon\rho\nabla{\boldsymbol{u}}=\epsilon\nabla{\boldsymbol{m}}-(\epsilon\nabla\rho)\otimes{\boldsymbol{u}},\quad{\mathbb{G}}_{{\boldsymbol{B}}}:=\epsilon\nabla{\boldsymbol{B}}.$

An equation governing the specific entropy is derived in the following lemma.

######  Lemma 2.

The following identity holds for any specific entropy $s(\rho,e)$,

$\rho({\partial}_{t}s+{\boldsymbol{u}}{\cdot}\nabla
s)-\nabla{\cdot\,}(\rho\epsilon\nabla
s)-(\epsilon\nabla\rho){\cdot}\nabla(es_{e}-\rho
s_{\rho})+{\boldsymbol{l}}{\cdot}\nabla
s_{e}-s_{e}{\mathbb{G}}_{{\boldsymbol{u}}}:\nabla{\boldsymbol{u}}-s_{e}{\mathbb{G}}_{{\boldsymbol{B}}}:\nabla{\boldsymbol{B}}=0.$

###### Proof.

Denote
${\mathcal{E}}_{{\boldsymbol{B}}}=\frac{1}{2}\rho^{-1}{\boldsymbol{B}}^{2}$,
and ${\mathcal{E}}_{{\boldsymbol{u}}}=\frac{1}{2}{\boldsymbol{u}}^{2}$. From
Definition (2.3), the specific internal energy can be expressed as
$e={\mathcal{E}}-{\mathcal{E}}_{{\boldsymbol{u}}}-{\mathcal{E}}_{{\boldsymbol{B}}}$.
First, we use this equality to derive an equation controlling the internal
energy $e$. Multiplying (3.12) with ${\boldsymbol{u}}$ gives

$\rho\left({\partial}_{t}{\mathcal{E}}_{{\boldsymbol{u}}}+{\boldsymbol{u}}{\cdot}\nabla{\mathcal{E}}_{{\boldsymbol{u}}}\right)+{\boldsymbol{u}}^{2}\epsilon\Delta\rho+{\boldsymbol{u}}{\cdot}\nabla
p-{\boldsymbol{u}}{\cdot}(\nabla{\cdot\,}{\boldsymbol{\beta}})-\epsilon{\boldsymbol{u}}{\cdot}\Delta{\boldsymbol{m}}=0.$
(3.16)

Multiplying (3.14) with ${\boldsymbol{B}}$ and writing
${\boldsymbol{B}}^{2}/2$ as $\rho{\mathcal{E}}_{{\boldsymbol{B}}}$, we end up
with the following equation governing the magnetic energy,

$\rho({\partial}_{t}{\mathcal{E}}_{{\boldsymbol{B}}}+{\boldsymbol{u}}{\cdot}\nabla{\mathcal{E}}_{{\boldsymbol{B}}})+{\mathcal{E}}_{{\boldsymbol{B}}}\epsilon\Delta\rho+{\mathcal{E}}_{{\boldsymbol{B}}}\rho\nabla{\cdot\,}{\boldsymbol{u}}-{\boldsymbol{B}}{\cdot}({\boldsymbol{B}}{\cdot}\nabla){\boldsymbol{u}}-\epsilon{\boldsymbol{B}}{\cdot}\Delta{\boldsymbol{B}}=0.$
(3.17)

An equation describing internal energy balance is obtained by subtracting
(3.16) and (3.17) from (3.13),

$\displaystyle\rho({\partial}_{t}e+{\boldsymbol{u}}{\cdot}\nabla
e)+\left(e-\frac{1}{2}{\boldsymbol{u}}^{2}\right)\epsilon\Delta\rho$
$\displaystyle\hskip
62.59596pt+p(\nabla{\cdot\,}{\boldsymbol{u}})-\nabla{\cdot\,}({\boldsymbol{u}}{\cdot}{\boldsymbol{\beta}})+{\boldsymbol{u}}{\cdot}(\nabla{\cdot\,}{\boldsymbol{\beta}})-{\mathcal{E}}_{{\boldsymbol{B}}}\rho\nabla{\cdot\,}{\boldsymbol{u}}+{\boldsymbol{B}}{\cdot}({\boldsymbol{B}}{\cdot}\nabla){\boldsymbol{u}}$
$\displaystyle\hskip 62.59596pt-\epsilon\Delta
E+\epsilon{\boldsymbol{u}}{\cdot}\Delta{\boldsymbol{m}}+\epsilon{\boldsymbol{B}}{\cdot}\Delta{\boldsymbol{B}}=0.$

By simple manipulation, we can show that

$-\nabla{\cdot\,}({\boldsymbol{u}}{\cdot}{\boldsymbol{\beta}})+{\boldsymbol{u}}{\cdot}(\nabla{\cdot\,}{\boldsymbol{\beta}})-{\mathcal{E}}_{{\boldsymbol{B}}}\rho\nabla{\cdot\,}{\boldsymbol{u}}+{\boldsymbol{B}}{\cdot}({\boldsymbol{B}}{\cdot}\nabla){\boldsymbol{u}}=0.$

Therefore, the equation can be simplified as

$\displaystyle\rho({\partial}_{t}e+{\boldsymbol{u}}{\cdot}\nabla
e)+\left(e-\frac{1}{2}{\boldsymbol{u}}^{2}\right)\epsilon\Delta\rho+p(\nabla{\cdot\,}{\boldsymbol{u}})$
(3.18) $\displaystyle\hskip 62.59596pt-\epsilon\Delta
E+\epsilon{\boldsymbol{u}}{\cdot}\Delta{\boldsymbol{m}}+\epsilon{\boldsymbol{B}}{\cdot}\Delta{\boldsymbol{B}}=0.$

Now, we use (3.18) to derive the desired equality. To utilize the chain rule
(2.5), we multiply (3.18) with $s_{e}$, (3.11) with $\rho s_{\rho}$ and add
them together to obtain

$\displaystyle\rho({\partial}_{t}s+{\boldsymbol{u}}{\cdot}\nabla
s)+(es_{e}-\rho
s_{\rho})\epsilon\Delta\rho+(\rho^{2}s_{\rho}+ps_{e})(\nabla{\cdot\,}{\boldsymbol{u}})$
(3.19)
$\displaystyle+s_{e}\left(-\frac{1}{2}{\boldsymbol{u}}^{2}\epsilon\Delta\rho-\epsilon\Delta
E+\epsilon{\boldsymbol{u}}{\cdot}\Delta{\boldsymbol{m}}+\epsilon{\boldsymbol{B}}{\cdot}\Delta{\boldsymbol{B}}\right)$
$\displaystyle=0,$

where $(\rho^{2}s_{\rho}+ps_{e})(\nabla{\cdot\,}{\boldsymbol{u}})$ is zero due
to the equation of state (2.7). The terms inside the last bracket can be
written as

$-\frac{1}{2}{\boldsymbol{u}}^{2}\Delta\rho-\Delta
E+{\boldsymbol{u}}{\cdot}\Delta{\boldsymbol{m}}+{\boldsymbol{B}}{\cdot}\Delta{\boldsymbol{B}}=-\nabla{\cdot\,}{\boldsymbol{l}}-\nabla{\boldsymbol{m}}:\nabla{\boldsymbol{u}}+\frac{1}{2}\nabla\rho{\cdot}\nabla{\boldsymbol{u}}^{2}-\nabla{\cdot\,}({\boldsymbol{B}}{\cdot}\nabla{\boldsymbol{B}})+{\boldsymbol{B}}{\cdot}\Delta{\boldsymbol{B}}.$

Because
$-\nabla{\cdot\,}({\boldsymbol{B}}{\cdot}\nabla{\boldsymbol{B}})+{\boldsymbol{B}}{\cdot}\Delta{\boldsymbol{B}}=-\nabla{\boldsymbol{B}}:\nabla{\boldsymbol{B}}$,
the equation (3.19) leads to

$\rho({\partial}_{t}s+{\boldsymbol{u}}{\cdot}\nabla s)+(es_{e}-\rho
s_{\rho})\epsilon\Delta\rho-
s_{e}\nabla{\cdot\,}{\boldsymbol{l}}-s_{e}{\mathbb{G}}_{{\boldsymbol{u}}}:\nabla{\boldsymbol{u}}-s_{e}{\mathbb{G}}_{{\boldsymbol{B}}}:\nabla{\boldsymbol{B}}=0.$

Applying the product rules $(es_{e}-\rho
s_{\rho})\epsilon\Delta\rho=-(\epsilon\nabla\rho){\cdot}\nabla(es_{e}-\rho
s_{\rho})+\nabla{\cdot\,}\left[(es_{e}-\rho
s_{\rho})\epsilon\nabla\rho\right]$ and
$s_{e}\nabla{\cdot\,}{\boldsymbol{l}}=\nabla{\cdot\,}({\boldsymbol{l}}s_{e})-{\boldsymbol{l}}{\cdot}\nabla
s_{e}$, we have

$\displaystyle\rho({\partial}_{t}s+{\boldsymbol{u}}{\cdot}\nabla
s)-(\epsilon\nabla\rho){\cdot}\nabla(es_{e}-\rho
s_{\rho})+\nabla{\cdot\,}\left[(es_{e}-\rho
s_{\rho})\epsilon\nabla\rho-{\boldsymbol{l}}s_{e}\right]$
$\displaystyle+{\boldsymbol{l}}{\cdot}\nabla
s_{e}-s_{e}{\mathbb{G}}_{{\boldsymbol{u}}}:\nabla{\boldsymbol{u}}-s_{e}{\mathbb{G}}_{{\boldsymbol{B}}}:\nabla{\boldsymbol{B}}=0.$

Combining with (3.15), we have the statement proved. ∎

The following result is useful for the main theorem of the minimum entropy
principle.

######  Lemma 3.

The quadratic form

$J_{1}(\nabla\rho,\nabla e):=-(\epsilon\nabla\rho)\cdot\nabla\left(es_{e}-\rho
s_{\rho}\right)+\boldsymbol{l}\cdot\nabla s_{e}+\epsilon\nabla\rho\cdot\nabla
s$

is negative definite.

###### Proof.

Using (3.15) and $\nabla s=s_{e}\nabla e+s_{\rho}\nabla\rho$, we can express
$J_{1}$ as

$\displaystyle J_{1}=$ $\displaystyle-\epsilon s_{e}\nabla\rho{\cdot}\nabla
e-\epsilon e\nabla\rho{\cdot}\nabla s_{e}+\epsilon
s_{\rho}|\nabla\rho|^{2}+\epsilon\rho\nabla\rho{\cdot}\nabla s_{\rho}$
$\displaystyle+\epsilon e\nabla\rho{\cdot}\nabla s_{e}-\epsilon s_{e}^{-1}\rho
s_{\rho}\nabla\rho{\cdot}\nabla s_{e}$ $\displaystyle+\epsilon\rho
s_{e}^{-1}\nabla s_{e}{\cdot}\left(s_{\rho}\nabla\rho+s_{e}\nabla
e\right)+\epsilon\nabla\rho{\cdot}\left(s_{\rho}\nabla\rho+s_{e}\nabla
e\right).$

Grouping terms with regards to quadratic terms of $(\nabla\rho,\nabla e)$
gives

$\displaystyle J_{1}=(\epsilon\rho s_{\rho\rho}+2\epsilon
s_{\rho})|\nabla\rho|^{2}+2\epsilon\rho s_{\rho e}\nabla\rho{\cdot}\nabla
e+\epsilon\rho s_{ee}|\nabla e|^{2}.$

That allows us to rewrite $J_{1}$ in matrix form as

$J_{1}=\begin{pmatrix}\nabla\rho\\\ \nabla
e\end{pmatrix}\left(\begin{pmatrix}\epsilon\rho^{-1}{\partial}_{\rho}(\rho^{2}s_{\rho})&\epsilon\rho
s_{\rho e}\\\ \epsilon\rho s_{\rho e}&\epsilon\rho
s_{ee}\end{pmatrix}\otimes{\mathbb{I}}_{d}\right)\begin{pmatrix}\nabla\rho&\nabla
e\end{pmatrix}.$ (3.20)

From the convex entropy inequalities (2.8), we can see that the $2\times 2$
matrix

$\begin{pmatrix}\epsilon\rho^{-1}{\partial}_{\rho}(\rho^{2}s_{\rho})&\epsilon\rho
s_{\rho e}\\\ \epsilon\rho s_{\rho e}&\epsilon\rho s_{ee}\end{pmatrix}$ (3.21)

has positive determinant and a negative trace. They are sufficient conditions
for negative definiteness of (3.21). The lemma is proved. Therefore, we always
have $J_{1}\leq 0$ for all $(\nabla\rho,\nabla e)$. ∎

We proceed as in [11]. The following theorem completes the purpose of this
section.

######  Theorem 3.2 (Minimum entropy principle).

Upon smoothness of the solution to (3.3)-(3.6), and the assumptions of Theorem
3.1, and ${\boldsymbol{B}}\in L^{1}(D)$. We further assume that the density
and the internal energy uniformly converge to constant states $\rho^{*},e^{*}$
and remain constant outside of a compact set of interest. For all $t>0$, the
minimum entropy principle holds:

$\inf_{{\boldsymbol{x}}\in{\mathbb{R}}^{d}}s({\boldsymbol{x}},t)\geq\inf_{\boldsymbol{x}\in\mathbb{R}^{d}}s_{0}(\boldsymbol{x}).$

###### Proof.

From Lemmas 2 and 3, we have

$\rho\left(\partial_{t}s+\boldsymbol{u}\cdot\nabla
s\right)-\nabla\cdot(\epsilon\rho\nabla s)-a\nabla\rho\cdot\nabla
s=-J_{1}+s_{e}{\mathbb{G}}_{{\boldsymbol{u}}}:\nabla\boldsymbol{u}+s_{e}{\mathbb{G}}_{{\boldsymbol{B}}}:\nabla{\boldsymbol{B}}\geq
0.$ (3.22)

If the infimum of $s({\boldsymbol{x}},t)$ is reached outside of the compact
set in the assumption, then the result follows readily since
$\inf_{{\boldsymbol{x}}\in{\mathbb{R}}^{d}}s({\boldsymbol{x}},t)=s^{*}\geq\inf_{{\boldsymbol{x}}\in{\mathbb{R}}^{d}}s_{0}(\boldsymbol{x})$,
with $s({\boldsymbol{x}},t)\to s^{*}$ as ${\boldsymbol{x}}\to\infty$ due to
the uniform convergence and smoothness assumptions on density and internal
energy. Otherwise, at the point $\bar{{\boldsymbol{x}}}(t)$ where
$s({\boldsymbol{x}},t)$ reaches its infimum, by smoothness assumptions, we
have $\nabla s(\bar{{\boldsymbol{x}}}(t),t)=0$ and $\Delta
s(\bar{{\boldsymbol{x}}}(t),t)>0$. From (3.22), it is clear that

$\rho\partial_{t}s\left(\bar{{\boldsymbol{x}}}(t),t\right)-\epsilon\rho\Delta
s\left(\bar{{\boldsymbol{x}}}(t),t\right)\geq 0.$

This says that $\rho\partial_{t}s\left(\bar{{\boldsymbol{x}}}(t),t\right)\geq
0$, which leads to the conclusion since $\rho>0$ by Theorem 3.1. ∎

######  Remark 3.3 (Positivity of internal energy).

For completeness, we repeat this argument by [11]. The minimum entropy
principle and the positivity of density lead to the positivity of specific
internal energy $e>0$ and therefore the positivity of internal energy $\rho
e>0$. From the equation of state (2.7), we can write $s=s(e,\rho^{-1})$ and
$e=e(s,\rho^{-1})$. By [11, Appendix A.1], we have that the assumption of
positive temperature $s_{e}>0$ leads to $e_{s}>0$ if $\rho>0$. Combining this
with the minimum entropy principle, at time $t>0$, we have $e(s,\rho^{-1})\geq
e(\min_{\Omega}s,\rho^{-1})\geq
e(\min_{\Omega}s_{0},\rho^{-1}),\;\forall\rho>0$. This means that if the
specific internal energy is positive at $t=0$, then it remains positive at any
$t>0$. Only in the case of ideal gases, due to (2.9), the positivity of
internal energy implies the positivity of pressure.

###  3.3 Generalized entropy inequalities

Let $f(s)$ be a twice differentiable function. We consider a class of strictly
convex generalized entropies in the form $\rho f(s)$, as derived for the Euler
equations by [17]. The choice of the form $\rho f(s)$ can be motivated by the
unregularized equations of (3.3) and (3.19), which respectively give
${\partial}_{t}\rho+\nabla{\cdot\,}(\rho{\boldsymbol{u}})=0$ and
${\partial}_{t}s+{\boldsymbol{u}}{\cdot}\nabla s=0$. The latter equation leads
to ${\partial}_{t}f(s)+{\boldsymbol{u}}{\cdot}\nabla f(s)=0$ for any
differentiable function $f$. Multiplying
${\partial}_{t}\rho+\nabla{\cdot\,}(\rho{\boldsymbol{u}})=0$ with $f$ and
${\partial}_{t}f(s)+{\boldsymbol{u}}{\cdot}\nabla f(s)=0$ with $\rho$, and
adding them together leads to ${\partial}_{t}(\rho
f(s))+\nabla{\cdot\,}({\boldsymbol{u}}\rho f(s))$ = 0. Therefore, $\rho f(s)$
represents a large class of entropies for the MHD equations. We investigate
some important properties of this class in the following theorem. We note that
the results in Theorem 3.4 are well established for the Euler equations, see
[16, 17].

######  Theorem 3.4 (Inequalities of strictly convex generalized entropies).

The specific heat capacity at constant pressure $c_{p}$ is defined as
$c_{p}=T\frac{\partial s(p,T)}{\partial T}$. The generalized entropy $-\rho
f(s)$ is strictly convex if and only if

$f^{\prime}(s)>0,\quad\frac{f^{\prime}(s)}{c_{p}}-f^{\prime\prime}(s)>0.$
(3.23)

###### Proof.

See Appendix B. ∎

######  Lemma 4.

The following $2\times 2$ matrix is negative definite

$\begin{pmatrix}c_{p}^{-1}\rho
s_{\rho}^{2}+\rho^{-1}{\partial}_{\rho}(\rho^{2}s_{\rho})&c_{p}^{-1}s_{\rho}\rho
s_{e}+\rho s_{\rho e}\\\ c_{p}^{-1}\rho s_{\rho}s_{e}+\rho s_{\rho
e}&c_{p}^{-1}\rho s_{e}^{2}+\rho s_{ee}\end{pmatrix}.$

###### Proof.

The negative definiteness of the above $2\times 2$ matrix is shown in [11,
Lemma A.3] for the Euler equations. The proof remains valid for the MHD
equations because the only invoked necessities are the convexity assumption on
$-s$ and the thermodynamic identity (2.4). ∎

######  Theorem 3.5 (entropy inequalities).

Any smooth solution to the regularized system (3.3)-(3.6) satisfies the
entropy inequality

$\partial_{t}(\rho f(s))+\nabla\cdot(\boldsymbol{u}\rho
f(s)-\epsilon\rho\nabla f(s)-\epsilon f(s)\nabla\rho)\geq 0$

###### Proof.

Multiplying both sides of (3.22) with $f^{\prime}(s)$, we have

$\begin{array}[]{c}\rho\left(\partial_{t}f(s)+\boldsymbol{u}\cdot\nabla
f(s)\right)-\nabla\cdot(\epsilon\rho\nabla f(s))+\epsilon\rho
f^{\prime\prime}(s)|\nabla s|^{2}-\epsilon f^{\prime}(s)\nabla\rho\cdot\nabla
s\\\
+J_{1}f^{\prime}(s)=f^{\prime}(s)s_{e}{\mathbb{G}}_{{\boldsymbol{u}}}:\nabla\boldsymbol{u}+f^{\prime}(s)s_{e}{\mathbb{G}}_{{\boldsymbol{B}}}:\nabla{\boldsymbol{B}}.\end{array}$

Multiplying the density equation (3.3) with $\rho$ and adding it to the above
equation, the product rule for temporal and spatial derivatives gives

$\displaystyle{\partial}_{t}(\rho f(s))$
$\displaystyle+\nabla\cdot(\boldsymbol{u}\rho
f(s))-\nabla\cdot(\epsilon\rho\nabla f(s)+\epsilon f(s)\nabla\rho)$
$\displaystyle+\epsilon\rho f^{\prime\prime}(s)|\nabla
s|^{2}+J_{1}f^{\prime}(s)=f^{\prime}(s)s_{e}{\mathbb{G}}_{{\boldsymbol{u}}}:\nabla{\boldsymbol{u}}+f^{\prime}(s)s_{e}{\mathbb{G}}_{{\boldsymbol{B}}}:\nabla{\boldsymbol{B}},$

which is close to the inequality that we want to prove. Since we know that the
right hand side of the above equation is nonnegative, if $\epsilon\rho
f^{\prime\prime}(s)|\nabla s|^{2}+J_{1}f^{\prime}(s)\leq 0$ holds, then the
proof is complete. Denote $J_{2}:=\epsilon\rho c_{p}^{-1}|\nabla
s|^{2}+J_{1}$. Using the second inequality of (3.23), we can derive an upper
bound of the quantity of interest $\epsilon\rho f^{\prime\prime}(s)|\nabla
s|^{2}+J_{1}f^{\prime}(s)<f^{\prime}(s)\left(\epsilon\rho c_{p}^{-1}|\nabla
s|^{2}+J_{1}\right)=f^{\prime}(s)J_{2}$. Using the chain rule (2.5) on $\nabla
s$ and the matrix form of $J_{1}$ in (3.20), we have

$J_{2}=\epsilon\begin{pmatrix}\nabla\rho\\\ \nabla
e\end{pmatrix}\left(\begin{pmatrix}c_{p}^{-1}\rho
s_{\rho}^{2}+\rho^{-1}{\partial}_{\rho}(\rho^{2}s_{\rho})&c_{p}^{-1}s_{\rho}\rho
s_{e}+\rho s_{\rho e}\\\ c_{p}^{-1}\rho s_{\rho}s_{e}+\rho s_{\rho
e}&c_{p}^{-1}\rho s_{e}^{2}+\rho
s_{ee}\end{pmatrix}\otimes{\mathbb{I}}_{d}\right)\begin{pmatrix}\nabla\rho&\nabla
e\end{pmatrix}.$

Due to Lemma 4, $J_{2}$ is always nonpositive. Therefore, $\epsilon\rho
f^{\prime\prime}(s)|\nabla s|^{2}+J_{1}f^{\prime}(s)\leq 0$, which proves the
theorem. ∎

######  Remark 3.6 (Entropy principles in the fully discrete settings).

In the above analysis, entropy principles such as density positivity, internal
energy, and entropy minimum principles are proved at the PDE level when
$\nabla{\cdot\,}{\boldsymbol{B}}\equiv 0$. However, this condition is usually
violated in fully discrete approximations. The authors of [19, 5] showed that
a slight violation of the divergence-free condition leads to a violation of
the positivity property. Recently [39, 40] proposed to use the Godunov form of
ideal MHD instead, where the so-called Powell terms [29] are added to the
system. A good feature of the Godunov form is that a slight violation of the
divergence-free condition does not prevent the positivity property of the MHD.
Then the authors [39, 40] proved the entropy principles both at the PDE level
and at the fully discrete level for some DG schemes. Extending these works
within finite elements is our ongoing work and will be reported in a separate
article.

##  4 Numerical investigation

In this section, we demonstrate some numerical properties of the monolithic
parabolic flux (3.3)-(3.6) using a robust shock-capturing CG method.

###  4.1 Finite element (CG) discretization

For computation, instead of considering the unbounded spatial domain
${\mathbb{R}}^{d}$ in the main analysis, we consider an open bounded subset
$\Omega\in{\mathbb{R}}^{d}$. The domain $\Omega$ is discretized into $N_{e}$
disjoint triangle elements $K_{i},i=1,\dots,N_{e}$ being open sets in
${\mathbb{R}}^{d}$ such that
$\cup\\{\overline{K_{i}}\\}_{i=1}^{N_{e}}=\overline{\Omega}$, where
$\overline{K_{i}}$ is the closure of $K_{i}$, and all vertices of the polytope
$\overline{\Omega}\approx\Omega$ are contained by the boundary
${\partial}\Omega$ of $\Omega$. For a valid CG discretization, we require that
no hanging nodes are present, i.e., no vertices of any element lie on an edge
of any other element. The set of all vertices and all elements of this
partition constitutes the computational mesh ${\mathcal{T}}_{h}$.

We define a continuous Lagrange finite element function space
${\mathcal{Q}}_{h}$ as

${\mathcal{Q}}_{h}:=\\{v\in
C^{0}(\overline{\Omega})\;\big{|}\;v_{|_{K_{i}}}\in{\mathbb{P}}_{k}(K_{i}),\,\forall
i=1,\dots,N_{e}\\},$

where $C^{0}(\overline{\Omega})$ is the space of continuous functions on
$\overline{\Omega}$, and ${\mathbb{P}}_{k}(K_{i})$ is the space of polynomials
of at most $k$-th degree on $K_{i}$. The corresponding vector function space
${\mathbfcal V}_{h}$ is defined as ${\mathbfcal
V}_{h}:=[{\mathcal{Q}}_{h}]^{d}$. A mesh-size function
$h(x)\in{\mathcal{Q}}_{h}$ is defined through the following projection,

$(h,v_{h})=\left(\frac{h_{K}}{k},v_{h}\right),\,\forall
v_{h}\in{\mathcal{Q}}_{h},$ (4.1)

where the $L^{2}$–inner product
$(u,v)\equiv(u,v)_{\Omega}:=\int_{\Omega}uv\,\mathrm{d}x$ is defined for real-
valued functions $u,v$ in $L^{2}(\Omega)$, and $h_{K}$ is the circumradius of
element $K$.

A weak formulation of the regularized ideal MHD equations (3.1) reads: find
$\textsf{{U}}_{h}(t):=(\rho_{h}(t),{\boldsymbol{m}}_{h}(t),E_{h}(t),{\boldsymbol{B}}_{h}(t))^{\top}\in{\mathcal{C}}^{1}({\mathbb{R}}^{+},{\mathcal{Q}}_{h}\times{\mathbfcal
V}_{h}\times{\mathcal{Q}}_{h}\times{\mathbfcal V}_{h})$ such that

$\displaystyle({\partial}_{t}\textsf{{U}}_{h},\textsf{{V}}_{h})+(\nabla{\cdot\,}\textsf{{F}}_{{\mathcal{E}}}(\textsf{{U}}_{h}),\textsf{{V}}_{h})+(\nabla{\cdot\,}\textsf{{F}}_{{\mathcal{B}}}(\textsf{{U}}_{h}),\textsf{{V}}_{h})$
(4.2)
$\displaystyle+(\textsf{{F}}_{{\mathcal{V}}}^{m}(\textsf{{U}}),\nabla\textsf{{V}}_{h})-({\boldsymbol{n}}{\cdot}\textsf{{F}}_{{\mathcal{V}}}^{m}(\textsf{{U}}),\textsf{{V}}_{h})_{{\partial}\Omega}$
$\displaystyle=0,$

for all test functions $\textsf{{V}}_{h}\in{\mathcal{Q}}_{h}\times{\mathbfcal
V}_{h}\times{\mathcal{Q}}_{h}\times{\mathbfcal V}_{h}$, where the boundary
inner product
$(u,v)_{{\partial}\Omega}=\int_{{\partial}\Omega}uv\,\mathrm{d}s$ is a surface
integral, and the vector ${\boldsymbol{n}}$ is the pointing-outward normal
vector defined at every nodal point on the boundary ${\partial}\Omega$.
Solution of the system (4.2) is often said to be a viscous solution of the
ideal MHD system (3.1). The viscosity coefficient $\epsilon$ is constructed
such that it vanishes with mesh refinement. Therefore, as $h\to 0$, the
viscous solution of (4.2) converges to the weak solution of (3.1).

Our code is implemented in FEniCS, an open source finite element library, see
[23].

###  4.2 Time stepping

To proceed in time, we solve the ODE (4.2) using the strong stability
preserving Runge-Kutta schemes of order 3 when ${\mathbb{P}}_{1}$ elements are
used in space and order 4 when ${\mathbb{P}}_{3}$ elements are used in space,
see [30].

The time step size is adaptively chosen following a CFL condition,

$\Delta
t=\text{CFL}\left(\min_{{\boldsymbol{x}}\in\Omega}h({{\boldsymbol{x}}})\right)\left(\max_{i=1,8}|\Lambda_{i}|\right)^{-1},$

where $\Lambda_{i}$ corresponds to the $i$-th eigenvalue of the MHD system,
see [6, Equation (3.5)]. In the following tests, the CFL number is chosen to
be $0.3$.

###  4.3 Divergence cleaning

Satisfying the divergence-free condition (2.2) has been known as a challenging
task in numerically solving MHD. However, due to being not of the main focus,
we use the simplest divergence cleaning method for the numerical demonstration
in this paper: the projection method [3]. In each Runge-Kutta stage, the
following cleaning procedure is applied:

  1. 1.

Solve the Poisson equation
$\Delta\Psi_{h}-\nabla{\cdot\,}{\boldsymbol{B}}_{h}=0$ for
$\Psi_{h}\in{\mathcal{Q}}_{h}$.

  2. 2.

Calculate the projection ${\boldsymbol{B}}_{h}^{\prime}$ of the magnetic field
${\boldsymbol{B}}_{h}$ onto the divergence-free space:
${\boldsymbol{B}}_{h}^{\prime}={\boldsymbol{B}}_{h}-\nabla\Psi_{h}$.

  3. 3.

Use ${\boldsymbol{B}}_{h}^{\prime}$ as the magnetic field solution to proceed
in the next time step, and update the dependent numerical variables
accordingly to ensure consistency: pressure, temperature, energy, and other
entropy-related variables.

Despite the seemingly ad-hoc nature, in [36], conservation and accuracy
preserving properties of the projection method is proved.

###  4.4 Boundary conditions

In the benchmark tests in this paper, we use two basic types of boundary
conditions: one is Dirichlet, and the other is the periodic boundary. The
Dirichlet boundary conditions are injected into the solution vector in each
time step. The periodic boundary mapping is done via built-in functions in
FEniCS, see [23].

###  4.5 Comparison with the resistive MHD flux

It is natural and reasonable to use the resistive MHD flux, similar to using
the Navier-Stokes flux for the compressible Euler equations. We demonstrate
why using the resistive MHD flux is not suitable for artificial viscosity
methods.

The resistive model of the ideal MHD equations are obtained by replacing the
monolithic flux $\textsf{{F}}_{{\mathcal{V}}}^{m}(\textsf{{U}})$ in (3.1) with
the viscous flux $\textsf{{F}}_{{\mathcal{V}}}^{r}(\textsf{{U}})$ which we
call “resistive MHD flux”, see e.g., [2],

$\textsf{{F}}_{{\mathcal{V}}}^{r}(\textsf{{U}}):=\begin{pmatrix}0\\\
{\boldsymbol{\tau}}\\\ {\boldsymbol{u}}{\cdot}{\boldsymbol{\tau}}+\kappa\nabla
T+\eta{\boldsymbol{B}}{\cdot}\big{(}\nabla{\boldsymbol{B}}-\nabla{\boldsymbol{B}}^{\top}\big{)}\\\
\eta\big{(}\nabla{\boldsymbol{B}}-\nabla{\boldsymbol{B}}^{\top}\big{)}\end{pmatrix},$
(4.3)

where the viscous shear stress tensor ${\boldsymbol{\tau}}$ is

${\boldsymbol{\tau}}=\mu\Big{(}\nabla{\boldsymbol{u}}+\nabla{\boldsymbol{u}}^{\top}\Big{)}-\lambda\nabla{\cdot\,}{\boldsymbol{u}}{\mathbb{I}},$

and $\kappa,\eta,\mu,\lambda$ are different viscosity coefficients. The sign
of $\lambda$ is not determined. In applications, $\lambda$ is often neglected,
or is set as $\lambda=-\frac{2}{3}\mu$.

####  4.5.1 An example: contact waves

We consider a contact line problem as an example: there is a discontinuity in
the density, but the velocity, pressure, and magnetic field are constant
functions. Let $\rho({\boldsymbol{x}},0)=\rho_{0}({\boldsymbol{x}})$ be an
initial density field containing the contact discontinuity. We can verify that
given some uniform fields
${\boldsymbol{u}}({\boldsymbol{x}},t)={\boldsymbol{u}}_{0},{\boldsymbol{B}}({\boldsymbol{x}},t)={\boldsymbol{B}}_{0},p({\boldsymbol{x}},t)=p_{0}$,
where the functions ${\boldsymbol{u}}_{0},{\boldsymbol{B}}_{0},p_{0}$ are
constant, if a density solution $\rho({\boldsymbol{x}},t)$ solve the mass
equation

${\partial}_{t}\rho+\nabla{\cdot\,}(\rho{\boldsymbol{u}}_{0})=\epsilon\Delta\rho,$
(4.4)

then ${\boldsymbol{u}}_{0},{\boldsymbol{B}}_{0},p_{0},\rho$ also solve the
momentum, energy, and magnetic equations. Indeed, when the velocity, pressure,
and magnetic solutions are constant, the momentum equation (3.4), and the
magnetic equation (3.6) follows trivially. Inserting
${\boldsymbol{u}}_{0},{\boldsymbol{B}}_{0},p_{0},\rho$ into the energy
equation (3.5), we have

${\partial}_{t}(\rho e)+\nabla{\cdot\,}({\boldsymbol{u}}_{0}\rho
e)-\epsilon\Delta(\rho
e)+\frac{1}{2}{\boldsymbol{u}}_{0}^{2}\left({\partial}_{t}\rho+\nabla{\cdot\,}(\rho{\boldsymbol{u}}_{0})-\epsilon\Delta\rho\right)=0.$
(4.5)

For ideal gas, the equation of state (2.9) implies that $\rho
e=p_{0}/(\gamma-1)$ is a constant function. Combining this fact with (4.4), we
can see that (4.5) is fulfilled. Therefore, the monolithic parabolic flux is
compatible with contact lines.

The same conclusion cannot be drawn about the resistive MHD flux. We make the
same assumptions as above on velocity, pressure, magnetic field, and initial
density solution. A density solution $\rho({\boldsymbol{x}},t)$ to the MHD
equations regularized by the resistive MHD flux satisfies

${\partial}_{t}\rho+\nabla{\cdot\,}(\rho{\boldsymbol{u}}_{0})=0.$ (4.6)

In a similar manner, the momentum equation (3.4), and the magnetic equation
(3.6) follows trivially. However, inserting
${\boldsymbol{u}}_{0},{\boldsymbol{B}}_{0},p_{0},\rho$ into the energy
equation (3.5) gives

${\partial}_{t}(\rho e)+\nabla{\cdot\,}({\boldsymbol{u}}_{0}\rho
e)+\frac{1}{2}{\boldsymbol{u}}_{0}^{2}\left({\partial}_{t}\rho+\nabla{\cdot\,}(\rho{\boldsymbol{u}}_{0})\right)+\nabla{\cdot\,}(\kappa
p_{0}\rho^{-2}\nabla\rho)=0.$ (4.7)

The equations (2.9) and (4.6) imply that (4.7) only holds if the thermal
diffusivity is zero, i.e., $\kappa=0$. For the compressible Euler equations,
letting $\kappa=0$ in (4.3) is known to lead to Gibbs phenomenon to the
numerical solution, see [26]. This argument suggests that the resistive MHD
flux is not compatible with contact lines.

We numerically demonstrate this argument on a contact solution extracted from
a Riemann solver [35] of the Brio-Wu problem [4]. The spatial domain is one
dimensional $\Omega=[0,1]$. The gas constant is $\gamma=2$. The initial
solution contains a constant velocity field ${\boldsymbol{u}}=(u_{x},u_{y})$,
$u_{x}=0.5915470932$, $u_{y}=-1.5792628803$, constant pressure
$p=0.5122334291$, and a constant magnetic field
${\boldsymbol{B}}=(B_{x},B_{y})$, $B_{x}=0.75$, $B_{y}=-0.5349102426$. There
is a discontinuity in the density

$\rho=\begin{cases}0.7156521382,&{\boldsymbol{x}}\in[0,0.5],\\\
0.2348529760,&{\boldsymbol{x}}\in(0.5,1].\end{cases}$

In one dimension, the divergence-free condition (2.2) reduces to
${\partial}_{x}B_{x}=0$, which means that $B_{x}$ needs to remain constant
across $\Omega$ at all future time. Since the violation of this condition in
the numerical approximations of the Brio-Wu solution is typically negligible,
the divergence of $B_{x}$ is left untreated. At every nodal point, the
viscosity coefficients $\epsilon,\mu,\nu$ are chosen to be
$\frac{1}{2}h\max_{i=1,8}|\Lambda_{i}|$, $h$ is calculated by (4.1), and
$\lambda$ is set to $0$. The monolithic flux with the mentioned choice of
$\epsilon$ resembles the Lax-Friedrichs scheme or the upwind scheme. The
contact line at time $\widehat{t}=0.1$ is captured by the monolithic flux and
the resistive MHD flux in Figure 1 under multiple resolution levels. It can be
seen that the monolithic flux can capture the contact line without undershoots
and overshoots. This is not the case for the resistive MHD flux. Choosing
either $\kappa=0$ or $\kappa=1$ leads to overshoots and undershoots in the
numerical solutions. Increasing or decreasing $\kappa$ from the standard value
$\kappa=1$ does not help with the situation.

(a) Monolithic flux

(b) Resistive MHD flux, $\kappa=0$

(c) Resistive MHD flux, $\kappa=1$

(d) Resistive MHD flux, $\kappa=\frac{1}{2}$

(e) Resistive MHD flux, $\kappa=5$

Figure 1: Single contact solution captured by monolithic flux and resistive
MHD flux with different choices of the thermal diffusivity coefficient
$\kappa$. The solutions are shown under multiple resolution levels: 161, 321,
641 DOFs. First order viscosity. Final time $\widehat{t}=0.1$.

####  4.5.2 Single waves from Brio-Wu problem [4]

The Brio-Wu problem is a one dimensional Riemann problem, $\Omega=[0,1]$. The
initial profile is given by

$(\rho,u,p,B_{x},B_{y})=\begin{cases}(1,0,1,0.75,1),&{\boldsymbol{x}}\in[0,0.5],\\\
(0.125,0,0.1,0.75,-1),&{\boldsymbol{x}}\in(0.5,1].\end{cases}$

The adiabatic constant is $\gamma=2$. The well-known Brio-Wu problem is an
essential yet challenging test to examine if a numerical method can capture
different MHD wave structures accurately: the shocks, the rarefactions, and
the contact lines. The density solutions at the final time $\hat{t}=0.1$
comparing monolithic flux and resistive MHD flux are shown in Figure 2 with
different mesh resolutions. The viscosity coefficients
$\epsilon,\mu,\nu,\lambda$ are calculated at nodal points same to Section
4.5.1. Similar to the previous example, no divergence cleaning procedure is
used for the current one dimensional tests. Visualization of density solution
in Figure 2 shows that using the resistive MHD flux when $\kappa=1$ can be
numerically sufficient to capture the whole compound structure of the Brio-Wu
solution. We show the result by the resistive MHD flux with $\kappa=0$ in
Figure 2(b), which clearly shows that the solution is polluted by spurious
oscillations when the thermal diffusivity is zero. We also test an interesting
setting of the monolithic flux when we drop the regularization to the mass
equation, i.e., set $\epsilon\Delta\rho$ to zero in
$\textsf{{F}}_{{\mathcal{V}}}^{m}(\textsf{{U}})$, and show the result in
Figure 2(c). Under this setting, the compound structure of the Brio-Wu
solution is still captured without the unphysical oscillations, as opposed to
the intuition that mass regularization is the key to eliminating this
phenomenon. In contrast to the distinguishable behaviors resulting from the
two fluxes in Figure 1, due to the combination of different waves, one can
hardly decide whether the monolithic flux in Figure 2(a) or the resistive MHD
flux with $\kappa=1$ in Figure 2(d) is better than the other.

(a) Monolithic flux

(b) Resistive MHD flux, $\kappa=0$

(c) Monolithic flux, no mass regularization

(d) Resistive MHD flux, $\kappa=1$

Figure 2: The Brio-Wu solutions with 161, 321, 641 DOFs comparing the
monolithic flux and the resistive MHD flux with first order viscosity. Final
time $\widehat{t}=0.1$.

Due to the above reason, we separate the waves of the Brio-Wu solution and
look into the behavior of each of the single waves. For this purpose, we use
an exact Riemann solver by [35] to extract the single wave solutions. The
initial profiles to generate them: the contact, fast rarefaction, intermediate
shock-slow rarefaction, and slow shock are given in Table 1. The results are
reported in Figure 3.

Table 1: Initial solutions for the single waves of the Brio-Wu problem. The
solutions are extracted from an exact Riemann solver [35]. Entries with “–”
indicates that the values on the left and the right states are the same.

Left state Right state Contact $\rho$ 0.7156521382 0.2348529760 $u_{x}$
0.5915470932 – $u_{y}$ -1.5792628803 – $p$ 0.5122334291 – $B_{x}$ 0.7500000000
– $B_{y}$ -0.5349102426 – Intermediate shock $\rho$ 0.6799272943 0.2348529760
$u_{x}$ 0.6288155014 0.5915470935 $u_{y}$ -0.2295748706 -1.5792628801 $p$
0.4623011255 0.5122334291 $B_{x}$ 0.7500000000 – $B_{y}$ 0.5900487481
-0.5349102425 Left state Right state Fast rarefaction $\rho$ 1.0000000000
0.6799272943 $u_{x}$ 0.0000000000 0.6288155014 $u_{y}$ 0.0000000000
-0.2295748706 $p$ 1.0000000000 0.4623011255 $B_{x}$ 0.7500000000 – $B_{y}$
1.0000000000 0.5900487481 Slow shock $\rho$ 0.2348529760 0.1168051318 $u_{x}$
0.5915470930 -0.2455906431 $u_{y}$ -1.5792628803 -0.1711653489 $p$
0.5122334291 0.0873180084 $B_{x}$ 0.7500000000 – $B_{y}$ -0.5349102426
-0.9001418247

(a) Monolithic

(b) Resistive, $\kappa=0$

(c) Resistive, $\kappa=1$

(d) Monolithic

(e) Resistive, $\kappa=0$

(f) Resistive, $\kappa=1$

(g) Monolithic

(h) Resistive, $\kappa=0$

(i) Resistive, $\kappa=1$

Figure 3: Simulation of single waves from Brio-Wu problem by the monolithic
flux and the resistive MHD flux. The solutions are shown under multiple
resolution levels: 161, 321, 641 DOFs. Final time $\widehat{t}=0.1$.

In Figure 3(a)-(c), all the viscous fluxes produce undershoot around
$x\in(0.5,0.6)$ for the fast rarefaction solutions. However, the undershoot
produced by the monolithic flux is smaller in both $L^{\infty}$ and $L^{2}$
senses compared to the other two settings by the resistive MHD flux. The same
can be said for the slow shock solutions in Figure 3(g)-(i). The relevance of
the undershoots and overshoots in the invariant domain set is however a
difficult topic on its own, see [12] for this discussion on the compressible
Euler solutions. Figure 3(d)-(f) show an intermediate shock followed by a slow
rarefaction wave by the given initial solution. In this case, however, there
is no clear benefits of using the monolithic flux over the resistive MHD flux
with $\kappa=1$. Notice that the intermediate shock is not present in the
reference solution in Figure 3(d)-(f). The reason is that by default the exact
Riemann solver [35] does not capture this phenomenon. To this day, the
existence of the intermediate shock is still a debatable topic, see e.g., [8,
25, 33, 35].

####  4.5.3 Discrete minimum principle

We investigate the minimum principle of the experimental CG solutions. For
this study, we employ a commonly used entropy function for ideal gas,

$S=\frac{\rho}{\gamma-1}s,$ (4.8)

where the thermodynamic entropy $s=\ln\frac{p}{\rho^{\gamma}}$ is obtained
upon assuming $c_{v}=1$ in (2.10). Because $s$ is linearly scaled by $c_{v}$,
we can investigate the discrete minimum entropy principle using any choice of
$c_{v}$.

######  Remark 4.1.

It has been shown by the authors of [14] that positivity-preserving properties
are impossible to achieve when the consistent mass matrix is used by CG
methods. Positivity-preserving properties are important and are known to be
connected to other invariant-domain preserving properties such as the minimum
entropy principles. Therefore, in this Section 4.5.3, to avoid any possible
effects of the consistent mass matrix, we lump the mass matrix in the linear
system yielded by (4.2).

A history plot of the entropy $s_{h}$ is presented in Figure 4. We note that
$s_{h}\in{\mathcal{Q}}_{h}$ is calculated pointwise from pressure and density,
$s_{h,i}=\ln\frac{p_{h,i}}{\rho_{h,i}^{\gamma}}$ at every nodal point $i$. The
viscosity coefficients $\epsilon,\mu,\nu$ are chosen to be first order, which
is same to the previous test cases. We demonstrate further the dilemma of
using whether $\kappa=0$ or $\kappa>0$ for the resistive MHD flux in Section
4.5.1. The dilemma is that setting $\kappa>0$ makes the flux inconsistent with
contact wave as pointed out in Section 4.5.1, and violation of minimum entropy
principle and entropy inequalities [11], while letting $\kappa=0$ would lead
to Gibbs phenomenon in numerical approximations [11, 26]. It is worth
mentioning that in [2], the authors show that by adding the Powell terms [29]
together with incorporating the generalized Lagrange multiplier (GLM) [7] to
the MHD system, the entropy inequality for the resulting resistive GLM-MHD
system is recovered. Figure 4(a) compares the behavior of $s_{h}$ between
using the monolithic flux $\textsf{{F}}_{{\mathcal{V}}}^{m}(\textsf{{U}})$ and
the resistive MHD flux $\textsf{{F}}_{{\mathcal{V}}}^{r}(\textsf{{U}})$ when
setting $\kappa=0$. Figure 4(b) shows a similar comparison but the parameter
$\kappa$ is set to $1$. In Figure 4(a) and 4(b), one can see the violation of
the minimum principle since $\min_{\Omega}(s_{h})$ is not monotonically
increasing in time. In both cases of $\kappa=0$ and $\kappa=1$, there is a
drastic fall of $\min_{\Omega}(s_{h})$ in the start up phase. It gradually
recovers after that. The portion of time when the minimum principle is
violated seems to be shortened by refining the mesh. Only in the case of the
monolithic flux, the discrete minimum entropy principle is fulfilled. We note
that the minimum principle would not be satisfied if the consistent mass
matrix was used. In the history plots of $\min_{\Omega}(s_{h})$ by the
monolithic flux in Figure 4(a) and 4(b), the minimum principle violation is in
the order of machine epsilon. Even though in this paper we have proved the
minimum principle in the continuum case along with the shown numerical
evidence, much work is needed to design a discretization that provably ensures
this property at the fully discrete level. Constructing such a CG method is
within our ongoing works.

(a) Monolithic vs resistive MHD flux, $\kappa=0$

(b) Monolithic vs resistive MHD flux, $\kappa=1$

Figure 4: Investigation of the discrete minimum entropy principle on the Brio-
Wu problem comparing the monolithic and the resistive MHD flux. The value
$\min_{\Omega}(s_{h})$ are captured at 1000 snapshots in time under multiple
resolution levels: 161, 321, 641 DOFs. The range of $\min_{\Omega}(s_{h})$ for
the monolithic flux is in the order of machine epsilon:
$\min_{\Omega}(s_{h})\in$[-4.6E-15, 0] using 161 DOFs,
$\min_{\Omega}(s_{h})\in$[-5.2E-15,0] using 321 DOFs, and
$\min_{\Omega}(s_{h})\in$[-6.3E-15,0] using 641 DOFs.

###  4.6 Entropy viscosity method

For high-order stabilization, the entropy viscosity method was proposed by
[10, 13]. For systems of conservation laws, entropy is a conserved quantity in
smooth solution regions, which is known as the “entropy equation”. In presence
of shocks and discontinuities, this equation becomes an inequality. The idea
of the entropy viscosity method is to use entropy residual – the violation of
the entropy equation as an indicator of the shock locations. Effectively, it
adds enough viscosity to stabilize the discontinuities, while being negligible
in the smooth regions to preserve high-order accuracy away from the
discontinuities. We employ the entropy function (4.8) to calculate the entropy
residual. In each time step, we seek the entropy residual $R_{h}$ in the
finite element space ${\mathcal{C}}^{1}({\mathbb{R}}^{+},{\mathcal{Q}}_{h})$.
A robust way to construct $R_{h}$ has shown to be by a nodal-based approach
[26],

$R_{h,i}(t):=\sum_{K\in{\mathcal{T}}_{h}}\frac{1}{|K|}\int_{K}|{\partial}_{t}S_{h}+\nabla{\cdot\,}({\boldsymbol{u}}_{h}S_{h})|\varphi_{i}\,\mathrm{d}{\boldsymbol{x}},$

where $|K|$ is the volume/area of the element $K$, and $S_{h}$ is computed as
in Section 4.5.3. In order to construct the necessary amount of artificial
viscosity to stabilize the solution, we need to compute a low-order viscosity
coefficient function $\epsilon_{h}^{L}\in{\mathcal{Q}}_{h}$ proportional to an
approximate maximum wave speed, and a high-order viscosity function
$\epsilon_{h}^{H}\in{\mathcal{Q}}_{h}$ based on the entropy residual $R_{h}$.

The low-order viscosity is calculated as described in Section 4.5.1, for every
node $i$, we compute

$\epsilon_{h,i}^{L}=c_{\max}h_{i}\max_{i=1,8}|\Lambda_{i}|,$

where $\max_{i=1,8}|\Lambda_{i}|$ approximates the maximum wave speed, $h_{i}$
is the corresponding nodal value of $h(x)$ by (4.1), and the parameter
$c_{\max}$ is set to be $0.5$.

The high-order viscosity is set to be proportional to the normalized entropy
residual,

$\epsilon_{h,i}^{H}=c_{E}h_{i}^{2}\frac{|R_{h}|}{\|\overline{R_{h}}-R_{h}\|_{\infty,\Omega}},$

where $\|\overline{R_{h}}-R_{h}\|_{\infty,\Omega}$ normalizes the unit of
$|R_{h}|$, the maximum variance
$\overline{R_{h}}:=\|R_{h}-\text{mean}(R_{h})\|_{\infty,\Omega}$, and $c_{E}$
is set to be 1.

The final artificial viscosity coefficient at each nodal $i$ is assembled as

$\epsilon_{h,i}=\min(\epsilon_{h,i}^{L},\epsilon_{h,i}^{H}).$

In the next sections, we investigate the accuracy and the shock-capturing
capability of the entropy viscosity method when it is incorporated in the
$\epsilon$ coefficient in the monolithic viscous flux (3.2).

####  4.6.1 Accuracy test

Consider a periodic smooth vortex problem on a rectangle domain
$\Omega=[-10,10]\times[-10,10]$. The reference solution is a stationary flow
with a vortex perturbation

$(\rho(t),{\boldsymbol{u}}(t),p(t),{\boldsymbol{B}}(t))=(\rho_{0},{\boldsymbol{u}}_{0}+\delta{\boldsymbol{u}},p_{0}+\delta
p,{\boldsymbol{B}}_{0}+\delta{\boldsymbol{B}}),$

where

$\displaystyle\rho_{0}$ $\displaystyle=1,$ $\displaystyle{\boldsymbol{u}}_{0}$
$\displaystyle=(1,1),$ $\displaystyle\delta{\boldsymbol{u}}$
$\displaystyle=\frac{\mu}{\pi\sqrt{2}}e^{(1-r^{2})/2}(-r_{2},r_{1}),$
$\displaystyle p_{0}$ $\displaystyle=0,$ $\displaystyle\delta p$
$\displaystyle=-\frac{\mu^{2}(1+r^{2})e^{1-r^{2}}}{8\pi^{2}},$
$\displaystyle{\boldsymbol{B}}_{0}$ $\displaystyle=(0.1,0.1),$
$\displaystyle\delta{\boldsymbol{B}}$ $\displaystyle=\frac{\mu
e^{(1-r^{2})/2}}{2\pi}(-r_{2},r_{1}),$

the vortex radius is $r=\sqrt{r_{1}^{2}+r_{2}^{2}}$,
$(r_{1},r_{2})=(x,y)-{\boldsymbol{u}}_{0}t$, and the vortex strength is
$\mu=5.389489439$. The adiabatic constant is $\gamma=\frac{5}{3}$. The errors
measured at final time $\widehat{t}=0.05$.

Table 2: Smooth vortex problem. L1, L2 error of the finite element solution
and convergence rates of velocity ${\boldsymbol{u}}$ and magnetic field
${\boldsymbol{B}}$ at final time $\widehat{t}=0.05$ using ${\mathbb{P}}_{1}$
elements.

#DOFs Entropy viscosity solution ${\boldsymbol{u}}$ Unregularized Galerkin
solution ${\boldsymbol{u}}$ L1 Rate L2 Rate L1 Rate L2 Rate 7442 6.11E-04 –
3.37E-03 – 5.89E-04 – 3.24E-03 – 29282 1.51E-04 2.01 8.36E-04 2.01 1.48E-04
1.99 8.18E-04 1.99 116162 3.75E-05 2.01 2.07E-04 2.01 3.71E-05 2.00 2.05E-04
2.00 462722 9.33E-06 2.01 5.16E-05 2.01 9.28E-06 2.00 5.13E-05 2.00 #DOFs
Entropy viscosity solution ${\boldsymbol{B}}$ Unregularized Galerkin solution
${\boldsymbol{B}}$ L1 Rate L2 Rate L1 Rate L2 Rate 7442 2.47E-02 – 2.77E-02 –
2.35E-02 – 2.60E-02 – 29282 6.10E-03 2.04 6.80E-03 2.05 5.90E-03 2.01 6.56E-03
2.01 116162 1.50E-03 2.03 1.68E-03 2.03 1.48E-03 2.01 1.64E-03 2.01 462722
3.73E-04 2.02 4.16E-04 2.02 3.70E-04 2.01 4.11E-04 2.01

Table 3: Smooth vortex problem. L1, L2 error of the finite element solution
and convergence rates of velocity ${\boldsymbol{u}}$ and magnetic field
${\boldsymbol{B}}$ at final time $\widehat{t}=0.05$ using ${\mathbb{P}}_{2}$
elements.

#DOFs Entropy viscosity solution ${\boldsymbol{u}}$ Unregularized Galerkin
solution ${\boldsymbol{u}}$ L1 Rate L2 Rate L1 Rate L2 Rate 7442 1.92E-04 –
1.14E-03 – 1.91E-04 – 1.13E-03 – 29282 3.41E-05 2.50 2.01E-04 2.50 3.41E-05
2.48 2.02E-04 2.49 116162 7.96E-06 2.10 4.69E-05 2.10 7.98E-06 2.10 4.70E-05
2.10 462722 1.99E-06 2.00 1.17E-05 2.00 2.00E-06 2.00 1.17E-05 2.00 #DOFs
Entropy viscosity solution ${\boldsymbol{B}}$ Unregularized Galerkin solution
${\boldsymbol{B}}$ L1 Rate L2 Rate L1 Rate L2 Rate 7442 7.95E-03 – 9.23E-03 –
7.83E-03 – 9.08E-03 – 29282 1.40E-03 2.53 1.54E-03 2.61 1.40E-03 2.51 1.54E-03
2.59 116162 3.18E-04 2.16 3.32E-04 2.23 3.18E-04 2.15 3.32E-04 2.22 462722
7.75E-05 2.04 8.01E-05 2.06 7.76E-05 2.04 8.03E-05 2.06

Table 4: Smooth vortex problem. L1, L2 error of the finite element solution
and convergence rates of velocity ${\boldsymbol{u}}$ and magnetic field
${\boldsymbol{B}}$ at final time $\widehat{t}=0.05$ using ${\mathbb{P}}_{3}$
elements.

#DOFs Entropy viscosity solution ${\boldsymbol{u}}$ Unregularized Galerkin
solution ${\boldsymbol{u}}$ L1 Rate L2 Rate L1 Rate L2 Rate 7442 1.13E-04 –
5.70E-04 – 1.14E-04 – 5.87E-04 – 29282 8.09E-06 3.80 4.89E-05 3.54 8.11E-06
3.82 4.95E-05 3.57 116162 5.29E-07 3.93 3.95E-06 3.63 5.18E-07 3.97 3.89E-06
3.67 462722 3.93E-08 3.75 3.96E-07 3.32 3.42E-08 3.92 3.67E-07 3.41 #DOFs
Entropy viscosity solution ${\boldsymbol{B}}$ Unregularized Galerkin solution
${\boldsymbol{B}}$ L1 Rate L2 Rate L1 Rate L2 Rate 7442 4.43E-03 – 4.43E-03 –
4.44E-03 – 4.51E-03 – 29282 2.82E-04 4.02 2.94E-04 3.96 2.81E-04 4.03 2.96E-04
3.98 116162 1.84E-05 3.96 2.35E-05 3.67 1.79E-05 3.99 2.26E-05 3.73 462722
1.60E-06 3.54 2.64E-06 3.16 1.36E-06 3.73 2.15E-06 3.41

Figure 5: Accuracy test on the smooth vortex problem: convergence history for
velocity ${\boldsymbol{u}}_{h}$ against the exact solution. The entropy
viscosity method using the monolithic flux does not destroy the accuracy of
high-order solutions.

####  4.6.2 Brio-Wu problem with the EV method

We again consider the Brio-Wu problem described in Section 4.5.2. It can be
seen that due to being a first order method, the solutions presented in
Section 4.5.2 are over dissipative for the given number of degrees of freedom
(DOFs). In this section, we want to demonstrate the convergence of the EV
method and its efficiency as a shock capturing technique. The problem settings
are kept same but the viscosity coefficients are now calculated by the EV
method. The result is presented in Figure 6. For comparison, we show a
solution under 641 nodes as it has been shown in Section 4.5.2. One can see
that under the same modest number of nodes, the EV solution is a much better
approximation of the expected solution compared to the first order solution.

Figure 6: EV solutions to the Brio-Wu problem with 641, 1281 DOFs. Final time
$\widehat{t}=0.1$. For comparison, the first order solution with 641 DOFs in
Section 4.5.2 is presented by the blue dashed line. The red line presents the
EV solution obtained under the same number of 641 DOFs. The green line
presents the EV solution on a twice finer mesh.

####  4.6.3 Orszag-Tang problem, [27]

The popular 2D Orszag-Tang benchmark is investigated in this section. We
consider the unit square domain $\Omega=[0,1]\times[0,1]$ with periodic
boundaries in both $x$\- and $y$\- directions. The initial solution is given
as follows,

$(\rho_{0},{\boldsymbol{u}}_{0},p_{0},{\boldsymbol{B}}_{0})=\left(\frac{25}{36\pi},(-\sin(2\pi
y),\sin(2\pi x)),\frac{5}{12\pi},\left(-\frac{\sin(2\pi
y)}{\sqrt{4\pi}},\frac{\sin(4\pi x)}{\sqrt{4\pi}}\right)\right).$

The adiabatic constant is $\gamma=\frac{5}{3}$. We compute the solution with
the entropy viscosity method using the monolithic parabolic flux. The density
solution $\rho_{h}$ and the artificial viscosity $\mu_{h}$ using
${\mathbb{P}}_{3}$ elements are respectively shown at time $t=0.5$ and $t=1.0$
in Figure 7 and 8. It can be seen that the artificial viscosity locally tracks
the shocks and effectively adds enough viscosity to stabilize the solution. At
$t=1.0$, the behavior of the Orszag-Tang solution is considered turbulence
[37]. It is worth noting that simulation of the Orszag-Tang problem after
$t=0.5$ is a challenging task that is not straightforwardly achieved by many
numerical methods due to the solution being prone to divergence blowups, see
[6, 15]. Nonetheless, the solution at $t=1.0$ is still well captured by the
described viscosity method.

(a) Density $\rho_{h}$

(b) Entropy viscosity $\mu_{h}$

(c) Magnetic pressure $\frac{1}{2}{\boldsymbol{B}}_{h}^{2}$

(d) Thermodynamic pressure $p_{h}$

Figure 7: Entropy viscosity solution of the Orszag-Tang problem at time
$\widehat{t}=0.5$, 90000 ${\mathbb{P}}_{3}$ nodes

(a) Density $\rho_{h}$

(b) Entropy viscosity $\mu_{h}$

(c) Magnetic pressure $\frac{1}{2}{\boldsymbol{B}}_{h}^{2}$

(d) Thermodynamic pressure $p_{h}$

Figure 8: Entropy viscosity solution of the Orszag-Tang problem at time
$\widehat{t}=1.0$, 90000 ${\mathbb{P}}_{3}$ nodes

The next two benchmarks in Section 4.6.4 and 4.6.5 are more challenging for
numerical MHD solvers. They contain strong discontinuities while having a low
plasma-beta number $\beta:=\frac{2p}{{\boldsymbol{B}}^{2}}$, small density,
and small gas pressure. This condition may easily lead to breakdowns of the
solvers due to the negativity of the pressure and density of the numerical
solutions.

####  4.6.4 MHD Rotor problem, [1]

The computational domain is the unit square $\Omega=[0,1]\times[0,1]$. The
initial pressure and magnetic fields are uniform,
$p_{0}=1,{\boldsymbol{B}}_{0}=\left(\frac{5}{\sqrt{4\pi}},0\right)$. For
${\boldsymbol{x}}=(x,y)\in\Omega$, we define a radius $r$ as
$r:=\sqrt{(x-0.5)^{2}+(y-0.5)^{2}}$. The density and velocity are defined as

$(\rho_{0},{\boldsymbol{u}}_{0})=\begin{cases}\left(10,\left(\frac{2}{r_{0}}\left(\frac{1}{2}-y\right),\frac{2}{r_{0}}\left(x-\frac{1}{2}\right)\right)^{\top}\right)&\text{
if }r<r_{0},\\\
\left(1+9f,\left(f\frac{2}{r}\left(\frac{1}{2}-y\right),f\frac{2}{r}\left(x-\frac{1}{2}\right)\right)^{\top}\right)&\text{
if }r_{0}\leq r<r_{1},\\\ (1,(0,0)^{\top})&\text{ otherwise},\end{cases}$

where $r_{0}=0.1,r_{1}=0.115$, and $f=\frac{r_{1}-r}{r_{1}-r_{0}}$. We use
periodic boundaries in both $x-$ and $y-$ directions. A ${\mathbb{P}}_{3}$
solution by the EV method at the final time $\widehat{t}=0.15$ with 90000
${\mathbb{P}}_{3}$ nodes is shown in Figure 9. The computed solution agrees
with the existing results, e.g., [1, 37].

(a) Density $\rho_{h}$

(b) Thermodynamic pressure $p_{h}$

(c) Contour plot of Mach number

(d) Magnetic pressure $\frac{1}{2}{\boldsymbol{B}}_{h}^{2}$

Figure 9: Entropy viscosity solution to the Rotor problem at time
$\widehat{t}=0.15$, 90000 ${\mathbb{P}}_{3}$ nodes

####  4.6.5 MHD Blast problem, [1]

We consider the same setup as [39]. See [39] for a more detailed description
of this problem.

The computational domain is $\Omega=[-0.5,0.5]$. Initially, the density is
uniform $\rho_{0}=1$; the magnetic field is uniform
${\boldsymbol{B}}_{0}=\left(\frac{100}{\sqrt{4\pi}},0\right)$; and the
velocity ${\boldsymbol{u}}_{0}$ is uniformly zero. For
${\boldsymbol{x}}=(x,y)\in\Omega$, we define a radius $r:=\sqrt{x^{2}+y^{2}}$.
There is a sharp jump in the initial pressure,

$p_{0}=\begin{cases}1000&\text{ if }r<0.1,\\\ 0.1&\text{
otherwise}.\end{cases}$

We use periodic boundaries in both $x-$ and $y-$ directions. We show a
numerical solution to this problem in Figure 10. The computed solution agrees
with the existing results, e.g., [1, 39].

(a) Density $\rho_{h}$

(b) Thermodynamic pressure $p_{h}$

(c) Velocity $|{\boldsymbol{u}}_{h}|$

(d) Magnetic pressure $\frac{1}{2}{\boldsymbol{B}}_{h}^{2}$

Figure 10: Contour plots of the entropy viscosity solution to the MHD Blast
wave problem at time $\widehat{t}=0.01$, 90000 ${\mathbb{P}}_{3}$ nodes

##  5 Conclusion

The purpose of this paper has been to investigate the continuous and numerical
properties of the monolithic parabolic regularization to the ideal MHD
equations. Despite having no physical motivation, this regularization has been
shown to be compatible with all the generalized entropies, fulfill minimum
entropy principles, and positivity of pressure and density. As demonstrated in
the CG context, the monolithic parabolic regularization also holds attractive
numerical properties which can be related to the continuous analysis.

A known shortcoming with this regularization is that it is not Galilean and
rotational invariant. This means that the regularization changes with a shift
or rotation of the reference frame. However, since the viscosity coefficient
to scale the regularization is artificial and vanishing, it is unclear how
this violation would affect the numerical solutions. During the investigation
in this paper, we have also made attempts to demonstrate the downsides of the
monolithic flux but so far we have not witnessed any significant pitfalls.

In future works, we may investigate Galilean invariant viscous regularizations
to the ideal MHD equations, or add control of nonzero divergence in the
regularization. Allowing small violations of the divergence is numerically
important, and has been done by several other methods such as [38, 39].

##  Appendix A Strict convexity of $-s$

We prove that the strict convexity of $-s$ with respect to $\rho^{-1}$ and $e$
implies (2.8). The properties (2.8) are borrowed from the Euler equations, see
[11] but for completeness, we want to give a short proof to demonstrate that
many entropy properties still hold in the presence of electromagnetism.

The second inequality $s_{ee}<0$ is rather obvious. We have
${\partial}_{\rho^{-1}}s=-\rho^{2}s_{\rho}$. Since we know that
${\partial}_{\rho^{-1}}({\partial}_{\rho^{-1}}s)<0$, this leads to
${\partial}_{\rho}(\rho^{2}s_{\rho})<0$. The Hessian matrix of $-s$ with
respect to $\rho^{-1}$ and $e$ is also positive definite. Hence, the
determinant of this matrix is positive, which says
$\rho^{2}{\partial}_{\rho}(\rho^{2}s_{\rho})s_{ee}-\rho^{4}s_{\rho e}^{2}>0$.
This inequality is equivalent to the third property in (2.8).

##  Appendix B Proof of Theorem 3.4

We include a proof of Theorem 3.4 since this result, to the best of our
knowledge, has not existed in the literature. Although the two inequalities in
(3.23) have been proved for the compressible Euler equation [17], it is
certainly nontrivial to show that they also hold for the MHD equations. In
this appendix, we follow the same procedure of [17] to show (3.23) for ideal
MHD.

###### Proof.

Recall that our solution vector is
${\boldsymbol{U}}:=(\rho,{\boldsymbol{m}},E,{\boldsymbol{B}})^{\top}$. An
entropy $S=-\rho f(s)$ is said to be strictly convex if the Hessian matrix
$S_{{\boldsymbol{U}}{\boldsymbol{U}}}$ is positive definite. The minus sign in
$S$ differs from [17] because the entropy $s$ in the thermodynamic identity
(2.4) is defined with opposite sign to $s$ in [17].

We express ${\boldsymbol{u}}$ and $e$ under the components of
${\boldsymbol{U}}$

${\boldsymbol{u}}=\frac{{\boldsymbol{m}}}{\rho},\quad
e=\frac{E}{\rho}-\frac{{\boldsymbol{m}}^{2}}{2\rho^{2}}-\frac{{\boldsymbol{B}}^{2}}{2\rho}.$

We have

$\rho_{{\boldsymbol{U}}}=(1,{\boldsymbol{0}}^{\top},0,{\boldsymbol{0}}^{\top})^{\top},\quad\rho_{{\boldsymbol{U}}{\boldsymbol{U}}}={\mathbb{O}},$

where ${\boldsymbol{0}}$ is the zero column vector and ${\mathbb{O}}$ is the
zero square matrix of appropriate sizes. Because
$(\rho_{{\boldsymbol{U}}{\boldsymbol{U}}})_{ij}=0,\;\forall i,j$, by
differentiation rules,

$\displaystyle\left(S_{{\boldsymbol{U}}{\boldsymbol{U}}}\right)_{ij}$
$\displaystyle=\frac{{\partial}^{2}(-\rho
f(s))}{{\partial}{\boldsymbol{U}}_{i}{\partial}{\boldsymbol{U}}_{j}}$
$\displaystyle=\frac{{\partial}}{{\partial}{\boldsymbol{U}}_{i}}\left(-f(s){\partial}_{{\boldsymbol{U}}_{j}}\rho-\rho
f^{\prime}(s){\partial}_{{\boldsymbol{U}}_{j}}s\right)$
$\displaystyle=-f^{\prime}(s)\left({\partial}_{{\boldsymbol{U}}_{i}}\rho{\partial}_{{\boldsymbol{U}}_{j}}s+{\partial}_{{\boldsymbol{U}}_{j}}\rho{\partial}_{{\boldsymbol{U}}_{i}}s+\rho\frac{{\partial}^{2}s}{{\partial}_{{\boldsymbol{U}}_{i}}{\partial}_{{\boldsymbol{U}}_{j}}}\right)-\rho
f^{\prime\prime}(s){\partial}_{{\boldsymbol{U}}_{i}}s{\partial}_{{\boldsymbol{U}}_{j}}s.$

Therefore, in matrix form, $S_{{\boldsymbol{U}}{\boldsymbol{U}}}$ can be
written as

$S_{{\boldsymbol{U}}{\boldsymbol{U}}}=-f^{\prime}(s){\boldsymbol{H}}-f^{\prime\prime}(s)\rho
s_{\boldsymbol{U}}s_{\boldsymbol{U}}^{\top},$

where
${\boldsymbol{H}}=\rho_{{\boldsymbol{U}}}s_{\boldsymbol{U}}^{\top}+s_{{\boldsymbol{U}}}\rho_{\boldsymbol{U}}^{\top}+\rho
s_{{\boldsymbol{U}}{\boldsymbol{U}}}$. Using the chain rule (2.5), we can
write $s_{{\boldsymbol{U}}}$ as

$s_{{\boldsymbol{U}}}=s_{\rho}\rho_{\boldsymbol{U}}+s_{e}e_{\boldsymbol{U}}.$

Consequently, we have

${\boldsymbol{H}}=(\rho
s_{\rho\rho}+2s_{\rho})\rho_{\boldsymbol{U}}\rho_{\boldsymbol{U}}^{\top}+\rho
s_{e\rho}(\rho_{\boldsymbol{U}}e_{\boldsymbol{U}}^{\top}+e_{\boldsymbol{U}}\rho_{\boldsymbol{U}}^{\top})+\rho
s_{ee}e_{\boldsymbol{U}}e_{\boldsymbol{U}}^{\top}-s_{e}{\boldsymbol{C}},$

where ${\boldsymbol{C}}=-(\rho
e_{{\boldsymbol{U}}{\boldsymbol{U}}}+\rho_{{\boldsymbol{U}}}e_{{\boldsymbol{U}}}^{\top}+e_{{\boldsymbol{U}}}\rho_{{\boldsymbol{U}}}^{\top})$.
We now compute $e_{\boldsymbol{U}}$, $e_{{\boldsymbol{U}}{\boldsymbol{U}}}$,
and ${\boldsymbol{C}}$ to determine ${\boldsymbol{H}}$ and
$s_{\boldsymbol{U}}$. Straightforward calculation gives

$e_{{\boldsymbol{U}}}=\begin{pmatrix}-\frac{E}{\rho^{2}}+\frac{{\boldsymbol{m}}^{2}}{\rho^{3}}+\frac{{\boldsymbol{B}}^{2}}{2\rho^{2}}\\\
-\frac{{\boldsymbol{m}}}{\rho^{2}}\\\ \frac{1}{\rho}\\\
-\frac{{\boldsymbol{B}}}{\rho}\end{pmatrix}=\frac{1}{\rho}\begin{pmatrix}-\frac{E}{\rho}+\frac{{\boldsymbol{m}}^{2}}{\rho^{2}}+\frac{{\boldsymbol{B}}^{2}}{2\rho}\\\
-\frac{{\boldsymbol{m}}}{\rho}\\\ 1\\\
-{\boldsymbol{B}}\end{pmatrix}=\frac{1}{\rho}\begin{pmatrix}\frac{1}{2}{\boldsymbol{u}}^{2}-e\\\
-{\boldsymbol{u}}\\\ 1\\\ -{\boldsymbol{B}}\end{pmatrix},$
$e_{{\boldsymbol{U}}{\boldsymbol{U}}}=\begin{pmatrix}\frac{2e}{\rho^{3}}-\frac{3{\boldsymbol{m}}^{2}}{\rho^{4}}-\frac{{\boldsymbol{B}}^{2}}{\rho^{3}}&\frac{2{\boldsymbol{m}}^{\top}}{\rho^{3}}&-\frac{1}{\rho^{2}}&\frac{{\boldsymbol{B}}^{\top}}{\rho^{2}}\\\
\frac{2{\boldsymbol{m}}}{\rho^{3}}&-\frac{1}{\rho^{2}}{\mathbb{I}}&{\boldsymbol{0}}&{\mathbb{O}}\\\
-\frac{1}{\rho^{2}}&{\boldsymbol{0}}^{\top}&0&{\boldsymbol{0}}^{\top}\\\
\frac{{\boldsymbol{B}}}{\rho^{2}}&{\mathbb{O}}&{\boldsymbol{0}}&-\frac{1}{\rho}{\mathbb{I}}\end{pmatrix},$

and

${\boldsymbol{C}}=\frac{1}{\rho}\begin{pmatrix}{\boldsymbol{u}}^{2}&-{\boldsymbol{u}}^{\top}&0&{\boldsymbol{0}}^{\top}\\\
-{\boldsymbol{u}}&{\mathbb{I}}&{\boldsymbol{0}}&{\mathbb{O}}\\\
0&{\boldsymbol{0}}^{\top}&0&{\boldsymbol{0}}^{\top}\\\
{\boldsymbol{0}}&{\mathbb{O}}&{\boldsymbol{0}}&\rho{\mathbb{I}}\end{pmatrix}.$

The key to this proof is to introduce the following invertible matrix

${\boldsymbol{P}}:=\begin{pmatrix}1&{\boldsymbol{u}}^{\top}&\frac{1}{2}{\boldsymbol{u}}^{2}+e&{\boldsymbol{0}}^{\top}\\\
{\boldsymbol{0}}&\rho{\mathbb{I}}&\rho{\boldsymbol{u}}&{\mathbb{O}}\\\
0&{\boldsymbol{0}}^{\top}&\rho&{\boldsymbol{0}}^{\top}\\\
{\boldsymbol{0}}&{\mathbb{O}}&{\boldsymbol{B}}&{\mathbb{I}}\end{pmatrix}.$

Applying a left multiplication with ${\boldsymbol{P}}$ and a right
multiplication with ${\boldsymbol{P}}^{\top}$ simplifies the structure of
$S_{{\boldsymbol{U}}{\boldsymbol{U}}}$. We have

${\boldsymbol{P}}\rho_{{\boldsymbol{U}}}=\begin{pmatrix}1\\\
{\boldsymbol{0}}\\\ 0\\\
{\boldsymbol{0}}\end{pmatrix},\quad{\boldsymbol{P}}e_{{\boldsymbol{U}}}=\begin{pmatrix}0\\\
{\boldsymbol{0}}\\\ 1\\\
{\boldsymbol{0}}\end{pmatrix},\quad{\boldsymbol{P}}{\boldsymbol{C}}{\boldsymbol{P}}^{\top}=\begin{pmatrix}0&{\boldsymbol{0}}^{\top}&0&{\boldsymbol{0}}^{\top}\\\
{\boldsymbol{0}}&\rho{\mathbb{I}}&{\boldsymbol{0}}&{\mathbb{O}}\\\
0&{\boldsymbol{0}}^{\top}&0&{\boldsymbol{0}}^{\top}\\\
{\boldsymbol{0}}&{\mathbb{O}}&{\boldsymbol{0}}&{\mathbb{I}}\end{pmatrix},$
${\boldsymbol{P}}{\boldsymbol{H}}{\boldsymbol{P}}^{\top}=\begin{pmatrix}2s_{\rho}+\rho
s_{\rho\rho}&{\boldsymbol{0}}^{\top}&\rho s_{e\rho}&{\boldsymbol{0}}^{\top}\\\
{\boldsymbol{0}}&-\rho s_{e}{\mathbb{I}}&{\boldsymbol{0}}&{\mathbb{O}}\\\ \rho
s_{e\rho}&{\boldsymbol{0}}^{T}&\rho s_{ee}&{\boldsymbol{0}}^{\top}\\\
{\boldsymbol{0}}&{\mathbb{O}}&{\boldsymbol{0}}&-s_{e}{\mathbb{I}}\end{pmatrix}.$

Combining all the above derivations gives

${\boldsymbol{P}}S_{{\boldsymbol{U}}{\boldsymbol{U}}}{\boldsymbol{P}}^{\top}=-f^{\prime}(s)\begin{pmatrix}2s_{\rho}+\rho
s_{\rho\rho}&{\boldsymbol{0}}^{\top}&\rho s_{e\rho}&{\boldsymbol{0}}^{\top}\\\
{\boldsymbol{0}}&-\rho s_{e}{\mathbb{I}}&{\boldsymbol{0}}&{\mathbb{O}}\\\ \rho
s_{e\rho}&{\boldsymbol{0}}^{\top}&\rho s_{ee}&{\boldsymbol{0}}^{\top}\\\
{\boldsymbol{0}}&{\mathbb{O}}&{\boldsymbol{0}}&-s_{e}{\mathbb{I}}\end{pmatrix}-f^{\prime\prime}(s)\rho\begin{pmatrix}s_{\rho}\\\
{\boldsymbol{0}}\\\ s_{e}\\\
{\boldsymbol{0}}\end{pmatrix}\begin{pmatrix}s_{\rho}&{\boldsymbol{0}}^{\top}&s_{e}&{\boldsymbol{0}}^{\top}\end{pmatrix}.$

This reveals that two eigenvalues of
${\boldsymbol{P}}S_{{\boldsymbol{U}}{\boldsymbol{U}}}{\boldsymbol{P}}^{\top}$
are $f^{\prime}(s)\rho s_{e}$ and $f^{\prime}(s)s_{e}$. Since
$\rho>0,s_{e}>0$, and the invertability of ${\boldsymbol{P}}$, we conclude
that $f^{\prime}(s)>0$ is a necessary condition for strict positivity of
$S_{{\boldsymbol{U}}{\boldsymbol{U}}}$. Disregarding the two detached
dimensions corresponding to the two known eigenvalues, the remaining
dimensions require that the following matrix ${\boldsymbol{M}}$ of size
$2\times 2$ must be positive definite

${\boldsymbol{M}}:=-f^{\prime}(s)\begin{pmatrix}2s_{\rho}+\rho
s_{\rho\rho}&\rho s_{e\rho}\\\ \rho s_{e\rho}&\rho s_{ee}\\\
\end{pmatrix}-f^{\prime\prime}(s)\rho\begin{pmatrix}s_{\rho}\\\
s_{e}\end{pmatrix}\begin{pmatrix}s_{\rho}&s_{e}\end{pmatrix}.$

We then can proceed as in [17, from (3.13)] to see that the positive
definiteness of ${\boldsymbol{M}}$ is equivalent to
$\frac{f^{\prime}(s)}{c_{p}}-f^{\prime\prime}(s)>0$. This is because the same
thermodynamic identity (2.4) holds for the MHD equations as it holds for the
Euler equations. ∎

## Acknowledgement

Some computations were performed on UPPMAX provided by the Swedish National
Infrastructure for Computing (SNIC) under project number SNIC 2021/22-233.

## References

  * Balsara and Spicer [1999] D. S. Balsara and D. S. Spicer.  A staggered mesh algorithm using high order Godunov fluxes to ensure solenoidal magnetic fields in magnetohydrodynamic simulations.  _J. Comput. Phys._ , 149(2):270–292, 1999.  ISSN 0021-9991.  doi: 10.1006/jcph.1998.6153.  URL https://doi.org/10.1006/jcph.1998.6153.
  * Bohm et al. [2018] M. Bohm, A. R. Winters, G. J. Gassner, D. Derigs, F. Hindenlang, and J. Saur.  An entropy stable nodal discontinuous Galerkin method for the resistive MHD equations. Part I: Theory and numerical verification.  _J. Comput. Phys._ , 422:108076, 35, 2018.  ISSN 0021-9991.  doi: 10.1016/j.jcp.2018.06.027.  URL https://doi.org/10.1016/j.jcp.2018.06.027.
  * Brackbill and Barnes [1980] J. U. Brackbill and D. C. Barnes.  The effect of nonzero $\nabla\cdot{\bf B}$ on the numerical solution of the magnetohydrodynamic equations.  _J. Comput. Phys._ , 35(3):426–430, 1980.  ISSN 0021-9991.  doi: 10.1016/0021-9991(80)90079-0.  URL https://doi.org/10.1016/0021-9991(80)90079-0.
  * Brio and Wu [1988] M. Brio and C. C. Wu.  An upwind differencing scheme for the equations of ideal magnetohydrodynamics.  _J. Comput. Phys._ , 75(2):400–422, 1988.  ISSN 0021-9991.  doi: 10.1016/0021-9991(88)90120-9.  URL https://doi.org/10.1016/0021-9991(88)90120-9.
  * Christlieb et al. [2015] A. J. Christlieb, Y. Liu, Q. Tang, and Z. Xu.  Positivity-preserving finite difference weighted ENO schemes with constrained transport for ideal magnetohydrodynamic equations.  _SIAM J. Sci. Comput._ , 37(4):A1825–A1845, 2015\.  ISSN 1064-8275.  doi: 10.1137/140971208.  URL https://doi.org/10.1137/140971208.
  * Dao and Nazarov [2021] T. A. Dao and M. Nazarov.  A high-order residual-based viscosity finite element method for the ideal mhd equations, 2021.  URL https://arxiv.org/abs/2112.08885.
  * Dedner et al. [2002] A. Dedner, F. Kemm, D. Kröner, C.-D. Munz, T. Schnitzer, and M. Wesenberg.  Hyperbolic divergence cleaning for the MHD equations.  _J. Comput. Phys._ , 175(2):645–673, 2002.  ISSN 0021-9991.  doi: 10.1006/jcph.2001.6961.  URL https://doi.org/10.1006/jcph.2001.6961.
  * Freistühler and Szmolyan [1995] H. Freistühler and P. Szmolyan.  Existence and bifurcation of viscous profiles for all intermediate magnetohydrodynamic shock waves.  _SIAM J. Math. Anal._ , 26(1):112–128, 1995.  ISSN 0036-1410.  doi: 10.1137/S0036141093247366.  URL https://doi.org/10.1137/S0036141093247366.
  * Godunov [1959] S. K. Godunov.  A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics.  _Mat. Sb. (N.S.)_ , 47 (89):271–306, 1959.
  * Guermond and Pasquetti [2008] J.-L. Guermond and R. Pasquetti.  Entropy-based nonlinear viscosity for Fourier approximations of conservation laws.  _C. R. Math. Acad. Sci. Paris_ , 346(13-14):801–806, 2008.  ISSN 1631-073X.  doi: 10.1016/j.crma.2008.05.013.  URL https://doi.org/10.1016/j.crma.2008.05.013.
  * Guermond and Popov [2014] J.-L. Guermond and B. Popov.  Viscous regularization of the Euler equations and entropy principles.  _SIAM J. Appl. Math._ , 74(2):284–305, 2014.  ISSN 0036-1399.  doi: 10.1137/120903312.  URL https://doi.org/10.1137/120903312.
  * Guermond and Popov [2016] J.-L. Guermond and B. Popov.  Invariant domains and first-order continuous finite element approximation for hyperbolic systems.  _SIAM J. Numer. Anal._ , 54(4):2466–2489, 2016\.  ISSN 0036-1429.  doi: 10.1137/16M1074291.  URL https://doi.org/10.1137/16M1074291.
  * Guermond et al. [2011] J.-L. Guermond, R. Pasquetti, and B. Popov.  Entropy viscosity method for nonlinear conservation laws.  _J. Comput. Phys._ , 230(11):4248–4267, 2011\.  ISSN 0021-9991.  doi: 10.1016/j.jcp.2010.11.043.  URL https://doi.org/10.1016/j.jcp.2010.11.043.
  * Guermond et al. [2017] J.-L. Guermond, B. Popov, and Y. Yang.  The effect of the consistent mass matrix on the maximum-principle for scalar conservation equations.  _J. Sci. Comput._ , 70(3):1358–1366, 2017.  ISSN 0885-7474.  doi: 10.1007/s10915-016-0285-7.  URL https://doi.org/10.1007/s10915-016-0285-7.
  * Guillet et al. [2019] T. Guillet, R. Pakmor, V. Springel, P. Chandrashekar, and C. Klingenberg.  High-order magnetohydrodynamics for astrophysics with an adaptive mesh refinement discontinuous galerkin scheme.  _Mon. Notices Royal Astron. Soc._ , 485(3):4209–4246, 2019.  ISSN 1365-2966.
  * Harten [1983] A. Harten.  On the symmetric form of systems of conservation laws with entropy.  _J. Comput. Phys._ , 49(1):151–164, 1983.  ISSN 0021-9991.  doi: 10.1016/0021-9991(83)90118-3.  URL https://doi.org/10.1016/0021-9991(83)90118-3.
  * Harten et al. [1998] A. Harten, P. D. Lax, C. D. Levermore, and W. J. Morokoff.  Convex entropies and hyperbolicity for general Euler equations.  _SIAM J. Numer. Anal._ , 35(6):2117–2127, 1998\.  ISSN 0036-1429.  doi: 10.1137/S0036142997316700.  URL https://doi.org/10.1137/S0036142997316700.
  * Hughes et al. [2010] T. J. R. Hughes, G. Scovazzi, and T. E. Tezduyar.  Stabilized methods for compressible flows.  _J. Sci. Comput._ , 43(3):343–368, 2010.  ISSN 0885-7474.  doi: 10.1007/s10915-008-9233-5.  URL https://doi.org/10.1007/s10915-008-9233-5.
  * Janhunen [2000] P. Janhunen.  A positive conservative method for magnetohydrodynamics based on HLL and Roe methods.  _J. Comput. Phys._ , 160(2):649–661, 2000.  ISSN 0021-9991.  doi: 10.1006/jcph.2000.6479.  URL https://doi.org/10.1006/jcph.2000.6479.
  * Kuzmin and Klyushnev [2020] D. Kuzmin and N. Klyushnev.  Limiting and divergence cleaning for continuous finite element discretizations of the MHD equations.  _J. Comput. Phys._ , 407:109230, 18, 2020.  ISSN 0021-9991.  doi: 10.1016/j.jcp.2020.109230.  URL https://doi.org/10.1016/j.jcp.2020.109230.
  * Lax [1971] P. Lax.  Shock waves and entropy.  In _Contributions to nonlinear functional analysis (Proc. Sympos., Math. Res. Center, Univ. Wisconsin, Madison, Wis., 1971)_ , pages 603–634, 1971.
  * Lax [1954] P. D. Lax.  Weak solutions of nonlinear hyperbolic equations and their numerical computation.  _Comm. Pure Appl. Math._ , 7:159–193, 1954.  ISSN 0010-3640.  doi: 10.1002/cpa.3160070112.  URL https://doi.org/10.1002/cpa.3160070112.
  * Logg et al. [2012] A. Logg, K.-A. Mardal, and G. N. Wells, editors.  _Automated solution of differential equations by the finite element method_ , volume 84 of _Lecture Notes in Computational Science and Engineering_.  Springer, Heidelberg, 2012.  ISBN 978-3-642-23098-1; 978-3-642-23099-8.  doi: 10.1007/978-3-642-23099-8.  URL https://doi.org/10.1007/978-3-642-23099-8.  The FEniCS book.
  * Mabuza et al. [2020] S. Mabuza, J. N. Shadid, E. C. Cyr, R. P. Pawlowski, and D. Kuzmin.  A linearity preserving nodal variation limiting algorithm for continuous Galerkin discretization of ideal MHD equations.  _J. Comput. Phys._ , 410:109390, 28, 2020.  ISSN 0021-9991.  doi: 10.1016/j.jcp.2020.109390.  URL https://doi.org/10.1016/j.jcp.2020.109390.
  * Myong and Roe [1997] R. S. Myong and P. L. Roe.  Shock waves and rarefaction waves in magnetohydrodynamics. part 2. the mhd system.  _J. Plasma Phys._ , 58(3):485–519, 1997.  ISSN 0022-3778.  doi: 10.1006/jcph.1998.6101.  URL https://doi.org/10.1006/jcph.1998.6101.
  * Nazarov and Larcher [2017] M. Nazarov and A. Larcher.  Numerical investigation of a viscous regularization of the Euler equations by entropy viscosity.  _Comput. Methods Appl. Mech. Engrg._ , 317:128–152, 2017\.  ISSN 0045-7825.
  * Orszag and Tang [1979] S. A. Orszag and C.-M. Tang.  Small-scale structure of two-dimensional magnetohydrodynamic turbulence.  _J. Fluid Mech._ , 90(1):129–143, 1979.  ISSN 0022-1120.  doi: 10.1017/S002211207900210X.  URL https://doi.org/10.1006/10.1017/S002211207900210X.
  * Perthame and Shu [1996] B. Perthame and C.-W. Shu.  On positivity preserving finite volume schemes for Euler equations.  _Numer. Math._ , 73(1):119–130, 1996.  ISSN 0029-599X.  doi: 10.1007/s002110050187.  URL https://doi.org/10.1007/s002110050187.
  * Powell et al. [1999] K. G. Powell, P. L. Roe, T. J. Linde, T. I. Gombosi, and D. L. De Zeeuw.  A solution-adaptive upwind scheme for ideal magnetohydrodynamics.  _J. Comput. Phys._ , 154(2):284–309, 1999.  ISSN 0021-9991.  doi: 10.1006/jcph.1999.6299.  URL https://doi.org/10.1006/jcph.1999.6299.
  * Ruuth [2006] S. J. Ruuth.  Global optimization of explicit strong-stability-preserving Runge-Kutta methods.  _Math. Comp._ , 75(253):183–207, 2006.  ISSN 0025-5718.  doi: 10.1090/S0025-5718-05-01772-2.  URL https://doi.org/10.1090/S0025-5718-05-01772-2.
  * Somov [2012] B. V. Somov.  _Plasma Astrophysics, Part I: Fundamentals and Practice_ , volume 391\.  Springer Science & Business Media, 2012.
  * Tadmor [1986] E. Tadmor.  A minimum entropy principle in the gas dynamics equations.  _Appl. Numer. Math._ , 2(3-5):211–219, 1986.  ISSN 0168-9274.  doi: 10.1016/0168-9274(86)90029-2.  URL https://doi.org/10.1016/0168-9274(86)90029-2.
  * Takahashi and Yamada [2014] K. Takahashi and S. Yamada.  Exact riemann solver for ideal magnetohydrodynamics that can handle all types of intermediate shocks and switch-on/off waves.  _Journal of Plasma Physics_ , 80(2):255–287, 2014\.
  * Tang and Xu [2000] H.-Z. Tang and K. Xu.  Positivity-preserving analysis of explicit and implicit Lax-Friedrichs schemes for compressible Euler equations.  _J. Sci. Comput._ , 15(1):19–28, 2000.  ISSN 0885-7474.  doi: 10.1023/A:1007593601466.  URL https://doi.org/10.1023/A:1007593601466.
  * Torrilhon [2002] M. Torrilhon.  Exact solver and uniqueness conditions for riemann problems of ideal magnetohydrodynamics.  Technical report, Zurich: Seminar for Applied Mathematics, ETH, 2002.
  * Tóth [2000] G. Tóth.  The $\nabla\cdot B=0$ constraint in shock-capturing magnetohydrodynamics codes.  _J. Comput. Phys._ , 161(2):605–652, 2000.  ISSN 0021-9991.  doi: 10.1006/jcph.2000.6519.  URL https://doi.org/10.1006/jcph.2000.6519.
  * Tricco et al. [2016] T. S. Tricco, D. J. Price, and M. R. Bate.  Constrained hyperbolic divergence cleaning in smoothed particle magnetohydrodynamics with variable cleaning speeds.  _J. Comput. Phys._ , 322:326–344, 2016.  ISSN 0021-9991.  doi: 10.1016/j.jcp.2016.06.053.  URL https://doi.org/10.1016/j.jcp.2016.06.053.
  * Winters et al. [2017] A. R. Winters, D. Derigs, G. J. Gassner, and S. Walch.  A uniquely defined entropy stable matrix dissipation operator for high Mach number ideal MHD and compressible Euler simulations.  _J. Comput. Phys._ , 332:274–289, 2017.  ISSN 0021-9991.  doi: 10.1016/j.jcp.2016.12.006.  URL https://doi.org/10.1016/j.jcp.2016.12.006.
  * Wu and Shu [2018] K. Wu and C.-W. Shu.  A provably positive discontinuous Galerkin method for multidimensional ideal magnetohydrodynamics.  _SIAM J. Sci. Comput._ , 40(5):B1302–B1329, 2018\.  ISSN 1064-8275.  doi: 10.1137/18M1168042.  URL https://doi.org/10.1137/18M1168042.
  * Wu and Shu [2021] K. Wu and C.-W. Shu.  Provably physical-constraint-preserving discontinuous Galerkin methods for multidimensional relativistic MHD equations.  _Numer. Math._ , 148(3):699–741, 2021.  ISSN 0029-599X.  doi: 10.1007/s00211-021-01209-4.  URL https://doi.org/10.1007/s00211-021-01209-4.