Monolithic parabolic regularization of the MHD equations and entropy principles

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.


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.

The ideal MHD equations
Consider a d-dimensional spatial domain x ∈ R d , d ∈ N * , and a temporal domain t ∈ [0, +∞). We define the space time domain D := R d × [0, +∞). Let us denote U := (ρ, m, E, B) a vector containing several conserved quantities of a conducting fluid: density ρ(x, t) : The velocity field u(x, t) is determined by the relation u(x, t) := m(x, t)/ρ(x, t). The governing electrodynamics and fluid mechanics of such fluid can be described in the following ideal MHD equations, where the nonlinear tensor fluxes F E (U) and F B (U) are defined as I denotes the identity matrix of size d × d, and p is the thermodynamic pressure. The symmetric term β is called the Maxwell stress tensor: In the absence of the magnetic field or the conductivity of the fluid ∇· F B (U) ≡ 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 B is implied from Faraday's law, 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 where the term 1 2 B 2 is often referred to as the magnetic pressure. The specific entropy s(ρ, e) and the thermodynamic pressure p are defined through the following thermodynamic identity, see [31,Section 12.2.4], Denote s e := ∂s ∂e and s ρ := ∂s ∂ρ . The following simple chain rule will be frequently used, ∂ α s = s e ∂ α e + s ρ ∂ α ρ, (2.5) where α = {x, t}. The equation (2.4) can be written as Combining (2.6) with the following direct consequence of (2.5), we deduce that s e := T −1 , s ρ := −pT −1 ρ −2 , 3 with the second equality being called the "equation of state" and being equivalent to 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 ρ −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, For ideal gases, it is common to use to compute pressure and temperature from the conserved variables, where γ > 1 is the adiabatic gas constant and 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.

Monolithic parabolic regularization
The ideal MHD equations can be regularized with monolithic parabolic terms as in the conservative form below, and 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 . For reference, we separately write the equations in (3.1) as For numerical purposes, the constant is often dependent on the resolution scale of the underlying scheme. The monolithic parabolic viscous terms can be linked to many classical numerical 4 schemes with attractive properties. For example, on a uniform grid in one spatial dimension {. . . , x i−1 , x i , x i+1 , . . . }, the well-known Lax-Friedrichs method applied on the unregularized mass conservation equation , and δt > 0 is the time step. Assuming sufficient smoothness, it is known that (3.8) is stable upon choosing δt = C LF δx, where C LF is a sufficiently small real number. By adding and subtracting ρ n i from the right-hand-side of (3.8), we can rewrite (3.8) as where LF = δx 2 2δt . One can see that (3.9) is a consistent discretization of (3.3) when is set to LF . If δt is chosen to be of order O(δx) following the stability condition, the viscosity coefficient LF is of order O(δx) and vanishes as δx → ∞.
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, where the vanishing viscosity coefficient reads U p = 1 2 δx u L ∞ (R d ) with u L ∞ (R d ) being the maximum density wave speed. Again, it can be seen that (3.10) is a consistent discretization of (3.3) when is set to U p .

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 in (3.3) is positive. Theorem 3.1 (Positivity of density, [11]). Upon the following justifiable assumptions: (i) u and ∇· u ∈ L 1 (D); (iii) there exists ρ min > 0 such that outside of a closed ball in R d centered at the origin with the density solution to (3.3) satisfies the following positivity property 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 ∆ρ, for an arbitrary positive value of , regularizes (3.3) and guarantees positivity of the density solution.
Since s −1 e > 0, multiplying both sides of the above equality with s −1 e gives the desired conclusion.
We introduce the following matrices, An equation governing the specific entropy is derived in the following lemma. Proof. Denote E B = 1 2 ρ −1 B 2 , and E u = 1 2 u 2 . From Definition (2.3), the specific internal energy can be expressed as e = E − E u − E B . First, we use this equality to derive an equation controlling the internal energy e. Multiplying (3.12) with u gives Multiplying (3.14) with B and writing B 2 /2 as ρE B , we end up with the following equation governing the magnetic energy, An equation describing internal energy balance is obtained by subtracting (3.16) and (3.17) from (3.13), By simple manipulation, we can show that Therefore, the equation can be simplified as 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 ρs ρ and add them together to obtain ρ(∂ t s + u·∇s) + (es e − ρs ρ ) ∆ρ + (ρ 2 s ρ + ps e )(∇· u) where (ρ 2 s ρ + ps e )(∇· u) is zero due to the equation of state (2.7). The terms inside the last bracket can be written as Combining with (3.15), we have the statement proved.
The following result is useful for the main theorem of the minimum entropy principle.
That allows us to rewrite J 1 in matrix form as From the convex entropy inequalities (2.8), we can see that the 2 × 2 matrix is positive definite and has a negative trace. They are sufficient conditions for negative definiteness of (3.21). The lemma is proved. Therefore, we always have J 1 ≤ 0 for all (∇ρ, ∇e).
We proceed as in [11]. The following theorem completes the purpose of this section. . We further assume that the density and the internal energy uniformly converge to constant states ρ * , e * and remain constant outside of a compact set of interest. For all t > 0, the minimum entropy principle holds: Proof. From Lemmas 2 and 3, we have If the infimum of s(x, t) is reached outside of the compact set in the assumption, then the result follows readily since inf x∈R d s(x, t) = s * ≥ inf x∈R d s 0 (x), with s(x, t) → s * as x → ∞ due to the uniform convergence and smoothness assumptions on density and internal energy. Otherwise, at the pointx(t) where s(x, t) reaches its infimum, by smoothness assumptions, we have ∇s(x(t), t) = 0 and ∆s(x(t), t) > 0. From (3.22), it is clear that This says that ρ∂ t s (x(t), t) ≥ 0, which leads to the conclusion since ρ > 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 ρe > 0. From the equation of state (2.7), we can write s = s(e, ρ −1 ) and e = e(s, ρ −1 ). By [11, Appendix A.1], we have that the assumption of positive temperature s e > 0 leads to e s > 0 if ρ > 0. Combining this with the minimum entropy principle, at time t > 0, we have e(s, 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.

Generalized entropy inequalities
Let f (s) be a twice differentiable function. We consider a class of strictly convex generalized entropies in the form ρ f (s), as derived for the Euler equations by [17]. The choice of the form ρ f (s) can be motivated by the unregularized equations of (3.3) and (3.19), which respectively give ∂ t ρ + ∇· (ρu) = 0 and ∂ t s + u·∇s = 0. The latter equation leads to ∂ t f (s) + u·∇ f (s) = 0 for any differentiable function f . Multiplying ∂ t ρ + ∇· (ρu) = 0 with f and ∂ t f (s) + u·∇ f (s) = 0 with ρ, and adding them together leads to ∂ t (ρ f (s)) + ∇· (uρ f (s)) = 0. Therefore, ρ 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].
Lemma 4. The following 2 × 2 matrix is negative definite Proof. The negative definiteness of the above 2 × 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).
Proof. Multiplying both sides of (3.22) with f (s), we have Multiplying the density equation (3.3) with ρ and adding it to the above equation, the product rule for temporal and spatial derivatives gives 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 ρ f (s)|∇s| 2 + J 1 f (s) ≤ 0 holds, then the proof is complete. Denote J 2 := ρc −1 p |∇s| 2 + J 1 . Using the second inequality of (3.23), we can derive an upper bound of the quantity of interest ρ f (s)|∇s| 2 Using the chain rule (2.5) on ∇s and the matrix form of J 1 in (3.20), we have Due to Lemma 4, J 2 is always nonpositive. Therefore, ρ f (s)|∇s| 2 + J 1 f (s) ≤ 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 ∇· B ≡ 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.

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.

Finite element (CG) discretization
For computation, instead of considering the unbounded spatial domain R d in the main analysis, we consider an open bounded subset Ω ∈ R d . The domain Ω is discretized into N e disjoint triangle elements K i , i = 1, . . . , N e being open sets in R d such that ∪{K i } N e i=1 = Ω, where K i is the closure of K i , and all vertices of the polytope Ω ≈ Ω are contained by the boundary ∂Ω of Ω. 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 T h .
We define a continuous Lagrange finite element function space Q h as where C 0 (Ω) is the space of continuous functions on Ω, and P k (K i ) is the space of polynomials of at most k-th degree on K i . The corresponding vector function space V h is defined as where the boundary inner product (u, v) ∂Ω = ∂Ω uv ds is a surface integral, and the vector n is the pointing-outward normal vector defined at every nodal point on the boundary ∂Ω. Solution of the system (4.2) is often said to be a viscous solution of the ideal MHD system (3.1). The viscosity coefficient is constructed such that it vanishes with mesh refinement. Therefore, as h → 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].

Time stepping
To proceed in time, we solve the ODE (4.2) using the strong stability preserving Runge-Kutta schemes of order 3 when P 1 elements are used in space and order 4 when P 3 elements are used in space, see [30].
The time step size is adaptively chosen following a CFL condition, where Λ 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.

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: 2. Calculate the projection B h of the magnetic field B h onto the divergence-free space: 3. Use B h 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.

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].

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 F m V (U) in (3.1) with the viscous flux F r V (U) which we call "resistive MHD flux", see e.g., [2], where the viscous shear stress tensor τ is τ = µ ∇u + ∇u − λ∇· uI, and κ, η, µ, λ are different viscosity coefficients. The sign of λ is not determined. In applications, λ is often neglected, or is set as λ = − 2 3 µ.

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 ρ(x, 0) = ρ 0 (x) be an initial density field containing the contact discontinuity. We can verify that given some uniform fields For ideal gas, the equation of state (2.9) implies that ρe = p 0 /(γ − 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 ρ(x, t) to the MHD equations regularized by the resistive MHD flux satisfies In a similar manner, the momentum equation (3.4), and the magnetic equation (3.6) follows trivially. However, inserting u 0 , B 0 , p 0 , ρ into the energy equation (3.5) gives The equations (2.9) and (4.6) imply that (4.7) only holds if the thermal diffusivity is zero, i.e., κ = 0. For the compressible Euler equations, letting κ = 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 Ω = [0, 1]. The gas constant is γ = 2. The initial solution contains a constant velocity field u = (u x , u y ), u x = 0.5915470932, u y = −1.5792628803, constant pressure p = 0.5122334291, and a constant magnetic field B = (B x , B y ), B x = 0.75, B y = −0.5349102426. There is a discontinuity in the density In one dimension, the divergence-free condition (2.2) reduces to ∂ x B x = 0, which means that B x needs to remain constant across Ω 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 , µ, ν are chosen to be  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 κ = 0 or κ = 1 leads to overshoots and undershoots in the numerical solutions. Increasing or decreasing κ from the standard value κ = 1 does not help with the situation.

Single waves from Brio-Wu problem [4]
The Brio-Wu problem is a one dimensional Riemann problem, Ω = [0, 1]. The initial profile is given by The adiabatic constant is γ = 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 timet = 0.1 comparing monolithic flux and resistive MHD flux are shown in Figure 2 with different mesh resolutions. The viscosity coefficients , µ, ν, λ 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 κ = 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 κ = 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 ∆ρ to zero in F m V (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 13 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 κ = 1 in Figure 2(d) is better than the other. 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.
In Figure 3(a)-(c), all the viscous fluxes produce undershoot around x ∈ (0.5, 0.6) for the fast rarefaction solutions. However, the undershoot produced by the monolithic flux is smaller in both L ∞ 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 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 κ = 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].

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, where the thermodynamic entropy s = ln p ρ γ is obtained upon assuming c v = 1 in (2.10). Because  [14] that positivity-preserving properties are impossible to achieve when the consistent mass matrix is used by CG methods. Positivitypreserving 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 ∈ Q h is calculated pointwise from pressure and density, s h,i = ln p h,i ρ γ h,i at every nodal point i. The viscosity coefficients , µ, ν are chosen to be first order, which is same to the previous test cases. We demonstrate further the dilemma of using whether κ = 0 or κ > 0 for the resistive MHD flux in Section 4.5.1. The dilemma is that setting κ > 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 κ = 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 F m V (U) and the resistive MHD flux F r V (U) when setting κ = 0. Figure 4(b) shows a similar comparison but the parameter κ is set to 1. In Figure 4(a) and 4(b), one can see the violation of the minimum principle since min Ω (s h ) is not monotonically increasing in time. In both cases of κ = 0 and κ = 1, there is a drastic fall of min Ω (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 Ω (s h ) by the monolithic flux in Figure 4 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.

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 residualthe 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 17 (4.8) to calculate the entropy residual. In each time step, we seek the entropy residual R h in the finite element space C 1 (R + , Q h ). A robust way to construct R h has shown to be by a nodal-based approach [26], 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 L h ∈ Q h proportional to an approximate maximum wave speed, and a high-order viscosity function H h ∈ 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 where max i=1,8 |Λ 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, and c E is set to be 1. The final artificial viscosity coefficient at each nodal i is assembled as In the next sections, we investigate the accuracy and the shock-capturing capability of the entropy viscosity method when it is incorporated in the coefficient in the monolithic viscous flux (3.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.   [27] The popular 2D Orszag-Tang benchmark is investigated in this section. We consider the unit square domain Ω = [0, 1] × [0, 1] with periodic boundaries in both x-and y-directions. The initial solution is given as follows, The adiabatic constant is γ = 5 3 . We compute the solution with the entropy viscosity method using the monolithic parabolic flux. The density solution ρ h and the artificial viscosity µ h using 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. 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 β := 2p 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.  shown in Figure 9. The computed solution agrees with the existing results, e.g., [1,37].

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 Ω = [−0.5, 0.5]. Initially, the density is uniform ρ 0 = 1; the magnetic field is uniform B 0 = 100 √ 4π , 0 ; and the velocity u 0 is uniformly zero. For x = (x, y) ∈ Ω, we define a radius r := x 2 + y 2 . There is a sharp jump in the initial pressure, 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].

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. 24 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].
Consequently, we have H = (ρs ρρ + 2s ρ )ρ U ρ U + ρs eρ (ρ U e U + e U ρ U ) + ρs ee e U e U − s e C, where C = −(ρe UU + ρ U e U + e U ρ U ). We now compute e U , e UU , and C to determine H and s U . Straightforward calculation gives The key to this proof is to introduce the following invertible matrix Applying a left multiplication with P and a right multiplication with P simplifies the structure of S UU . We have Combining all the above derivations gives

27
This reveals that two eigenvalues of PS UU P are f (s)ρs e and f (s)s e . Since ρ > 0, s e > 0, and the invertability of P, we conclude that f (s) > 0 is a necessary condition for strict positivity of S UU . Disregarding the two detached dimensions corresponding to the two known eigenvalues, the remaining dimensions require that the following matrix M of size 2 × 2 must be positive definite M := − f (s) 2s ρ + ρs ρρ ρs eρ ρs eρ ρs ee − f (s)ρ s ρ s e s ρ s e .
We then can proceed as in [17, from (3.13)] to see that the positive definiteness of M is equivalent to f (s) c p − f (s) > 0. This is because the same thermodynamic identity (2.4) holds for the MHD equations as it holds for the Euler equations.