{'fraction_non_alphanumeric': 0.058459081396629586, 'fraction_numerical': 0.021478136358629805, 'mean_word_length': 3.6945449844881075, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 31, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Mixed DG method and an HDG method for incompressible magnetohydrodynamics

In this paper we propose and analyze a mixed DG method and an HDG method for the stationary Magnetohydrodynamics (MHD) equations with two types of boundary (or constraint) conditions. The mixed DG method is based a recent work proposed by Houston et. al. for the linearized MHD. With two novel discrete Sobolev embedding type estimates for the discontinuous polynomials, we provide a priori error estimates for the method on the nonlinear MHD equations. In the smooth case, we have optimal convergence rate for the velocity, magnetic field and pressure in the energy norm, the Lagrange multiplier only has suboptimal convergence order. With the minimal regularity assumption on the exact solution, the approximation is optimal for all unknowns. To the best of our knowledge, this is the first a priori error estimates of DG methods for nonlinear MHD equations. In addition, we also propose and analyze the first divergence-free HDG method for the problem with several unique features comparing with the mixed DG method.


Introduction
MHD describes the interaction of electrically conducting fluids and electromagnetic fields [15,20,34]. Examples of such magneto-fluids include plasmas, liquid metals, and salt water or electrolytes. We refer to [20,33,27] for a more comprehensive discussion on the applications of the MHD system. The physical model is based on two principles: first, the motion of a conducting fluid in the presence of a magnetic field induces an electric current which also interacts with the existing electromagnetic field. Second, the Lorentz force generated by the current and the magnetic field also affects the motion of the fluid. The governing equations of the stationary incompressible MHD system can be written as: The domain Ω is a simply connected, bounded Lipschitz polyhedron in R 3 . We denote by n the unit outward normal vector on ∂Ω. The unknowns are the velocity u, the pressure p, magnetic field b and the Lagrange multiplier r associated with the divergence constraint on the magnetic field b. The functions f , g are external force terms. These equations are characterized by three dimensionless parameters: the hydrodynamic Reynolds number Re = ν −1 , the magnetic Reynolds number Rm = ν −1 m and the coupling number κ. We refer to [1,15,20] for further discussion of these parameters and their typical values. We consider two types of boundary (or constraint) conditions for the magnetic field b and the Lagrange multiplier r. The first type is Due to its significant role in applications, there have been many studies on the MHD equations see [5,7,8,2,3,18,37,20,21,22,35,32,38,39] and the references therein. Designing and analyzing numerical methods to solve this system is in general a challenging task due to multiple vector and scalar unknowns, to the various differential operators involved, and to the nonlinearities of the PDEs. To the best of the our knowledge, all existing finite element methods mentioned above for (1.1) use conforming elements to approximate the magnetic field b. As a consequence, local conservation does not hold for (1.1d). For the detailed explanation of importance of local conservation, we refer to [16].
In this paper, we propose and analyze a mixed DG method for the stationary incompressible MHD with two types of boundary (or constraint) conditions (1.2) and (1.3), which provides optimal convergent approximation to the velocity, pressure and magnetic field. Due to the nature of DG methods, the local conservation for both velocity and magnetic field is preserved. Our method is based on the IP-DG scheme proposed in [25] for the linearized MHD equation. The bottle-neck to extend the analysis in [25] to nonlinear MHD system comes from the nonlinear coupling term between the fluid and magnetic field. Namely, a key ingredient in the analysis is a Sobolev embedding like estimate for the L 3 -norm of the magnetic field. In Section 5 & 6, we develop the estimates for the discrete nonconforming magnetic field with the first type of boundary (or constraint) conditions (1.2) and the second type of boundary (or constraint) conditions (1.3), respectively. In fact, these L 3 -norm estimates of the discrete magnetic field help to show that our DG method has optimal approximation to all unknowns with the minimal regularity assumption on the exact solution. [38] is the first paper which tried to obtain optimal approximation under the minimal regularity assumption on the exact solution. However, according to [39,Section 1], the author of [38] failed to prove [38,Proposition 3.2] and [38,Corollary 3.1], which are indispensable to give an error estimate. In [39], a correct error analysis is given to show that with the first type of boundary (or constraint) conditions (1.2), optimal convergence is achieved by the conforming method in [38] under the minimal regularity assumption. Up to our knowledge, with the second type of boundary (or constraint) conditions (1.3), our DG method is the first numerical method shown to achieve optimal convergence of all unknowns under the minimal regularity assumption.
Comparing with conforming mixed methods, DG approach has several attractive features such as local conservation, high order accuracy, hp-adaptivity, easy implementation. Nevertheless DG methods are also critized with much more degrees of freedoms comparing with conforming methods. This disadvantage becomes more severe in problems involving multiple vector and scalar unkowns such as MHD. In order to make the DG approach more competitive, in Section 9 we propose a new hybridizable DG method (HDG) for the MHD problem. Thanks to the nature of the HDG framework, the scheme can be hybridized so that the global degree of freedoms can be reduced significantly and is more efficient than existing mixed methods when high order polynomial spaces are employed. In addition, the proposed HDG method provides exactly divergence-free velocity field while maintaining all existing features of HDG framework. As a consequence, the errors of the velocity and magnetic fields are independent of the pressure. Violation of divergence-free constraint can cause large errors in practice even for stable elements, we refer a review paper [28] for more discussions on this issue.
The rest of the paper is organized as follows: in Section 2 we present the mixedDG scheme for the MHD system and introduce the notations and definitions; in Section 3 we present our main results; Section 4 provides several auxiliary results which needed for the proofs; Section 5-8 are the detail proofs for the main results; In Section, we propose and analyze a new HDG method for the MHD system; concluding remarks are in Section 10.

Mixed Discontinuous Galerkin method
To define the DG method for the problem, we adopt notations and norms in [25]. We consider a family of conforming triangulations T h made of shape-regular tetrahedra. We denote by F I h the set of all interior faces of T h , and by F B h the set of all boundary faces. We define F h := F I h ∪ F B h . h K denotes the diameter of the element K, and h F is the diameter of the face F . The mesh size of T h is defined as h := max K∈T h h K . We denote by n K the unit outward normal vector on ∂K. We also introduce the average and jump operators. Let F = ∂K ∩ ∂K ′ be an interior face shared by K and K ′ . Let φ be a generic piecewise smooth function (scalar-, vector-or tensor-valued). We define the average of φ on F as { {φ} } := 1 2 (φ + φ ′ ) where φ and φ ′ denote the trace of φ from the interior of K and K ′ . Furthermore, let u be a piecewise smooth function and u a piecewise smooth vector-valued field. Analogously, we define the following jumps on F : On a boundary face F = ∂K ∩ ∂Ω, we set accordingly { {φ} } := φ, u := un, u := u ⊗ n, u T := n × u and u N := u · n.
Throughout this paper, we assume the integer k ≥ 1. Here P k (T h ; R 3 ) denotes the space contains vector-valued piecewise polynomials of degree no more than k on T h . Similarly, P k (T h ) denotes the space contains piecewise polynomials of degree no more than k on T h . In addition, standard inner product notations are used throughout the paper. Namely, (f, g) D = D f gdx for D ∈ R 3 and f, g D := D f gds for D ∈ R 2 . We use the standard notations for Sobolev norms. In addition, we use the following notations and spaces: In order to define mixed discontinuous Galerkin method, we introduce the finite dimensional spaces: . V h is for the approximation to the velocity field u. Q h is for the approximation to the pressure p. C h is for the approximation to the magnetic field b. S h is for the approximation to the Lagrange multiplier r with the first type of boundary (or constraint) conditions (1.2). S h ∩ L 2 0 (Ω) is for the approximation to the Lagrange multiplier r with the second type of boundary (constraint) conditions (1.3).
2.1. Mixed Discontinuous Galerkin method for the first type of boundary (or constraint) conditions. The mixed DG method for the first type of boundary (or constraint) conditions (1.2) for the magnetic field and the Lagrange multiplier seeks an approximation the exact solution (u, b, p, r) of (1.1) with (1.2). The method determines the approximate solution by requiring that it solves the following weak formulation: The postprocessing operator P from is defined on the element K by the following equations, see [6]: The parameter a 0 > 0 is a sufficiently large stabilization parameter; see [25,Proposition 2.4].
For the convective form, we take the usual upwind form defined by Here, u e is the value of the trace of u taken from the exterior of K, Next, we define the forms for the discretization of the Maxwell operator. The form M h for the curl-curl operator is given by As for the diffusion form, the stabilization parameter m 0 > 0 must be chosen large enough; see [25,Proposition 2.4]. The form D h for the divergence-free constraint on b is given by The form J h is the stabilization term that ensures the H 1 -conformity of the multiplier r h in a weak sense. It is given by with s 0 > 0 denoting a third stabilization parameter. Finally, the coupling form C h is defined by Mixed Discontinuous Galerkin method for the second type of boundary (or constraint) conditions. We can obtain the mixed DG formulation for the second type of boundary conditions (1.3) with slight modification of the above method. Namely, we use the same spaces for all unknowns except the Lagrange multiplier for which we use S h ∩ L 2 0 (Ω) instead. Now the formulation becomes: where the forms M I h , D I h and J I h are defined as:

Main results
We first present two novel discrete Sobolev embedding results which is the key ingredient for the analysis. For any b h ∈ C h , we define the discrete divergence of b h , denoted by ∇ h · b h to be the unique function in H 1 0 (Ω) ∩ S h satisfying: Similarly, for the second type of boundary condition, we need to modify the above discrete divergence operator as: In addition, we provide an analogue of Theorem 3.1 for the second type of boundary conditions as follows, Theorem 3.2. There is a positive constant C such that for any b h ∈ C h , we have . The norms that we are going to use in the error analysis are defined as follows: . Finally, we use the standard L 2 -norm for the pressure p.
then we have ∇ h · c h = 0 and therefore by Theorem 3.1, we have Now we are ready to present the well-posedness result: The proofs of the above results are in Section 5 and Section 6.
The norms that we are going to use in the error analysis are defined as follows: . Finally, we use the standard L 2 -norm for the pressure p.
Now we are ready to present the well-posedness results for the method: If the following quantities are all small enough, the system (1.1) with the first type of boundary (or constraint) conditions (1.2) has a unique weak solution (u, b, p, r) ∈ H 1 0 (Ω; Theorem 3.4. If the following quantities are all small enough, the system (1.1) with the second type of boundary (or constraint) conditions (1.3) has a unique weak solution (u, b, p, r) ∈ H 1 0 (Ω; R 3 ) × H(curl; Ω) × L 2 0 (Ω) × H 1 (Ω) ∩ L 0 2 (Ω) and the DG method (2.4) has a unique solution The last two results are for the convergence of the numerical solutions. We make minimal regularity assumptions on the exact solutions, see Remark 4.1 in [21]. Namely, we assume the exact solution (u, b, p, r) of (1.1) possesses the smoothness for σ, σ m > 1 2 . Now we are ready to state our main convergence results: Theorem 3.5. Let (u, b, p, r) be the exact solution of the system (1.1) with the first type of boundary (or constraint) conditions (1.2), and (u h , b h , p h , r h ) be the solution of the DG method (2.1). With the same assumption as in Theorem 3.3, in addition with the regularity assumption (3.7) and that 1 min(ν, ν m ) u H 1 (Ω) and 1 √ νκν m ∇ × b L 2 (Ω) are small enough, then we have here C depends on the physical parameters κ, ν, ν m and the external forces f , g but is independent of mesh size h.
here C depends on the physical parameters κ, ν, ν m and the external forces f , g but is independent of mesh size h.
Remark 3.2. The above results indicate that our method obtains same convergence rate as existing conforming methods [39,21]. Namely, if the exact solution is sufficient smooth, we have optimal convergence for u, b, p in the energy norms. The convergence rate for the Lagrange multiplier r is suboptimal in the discrete H 1 -norm. With minimal regularity assumption (3.7) the method is optimal for all unknowns.

Auxiliary results
In this section we gather some auxiliary results needed to carry out the error estimates in the next section.
Here, Π RT is the k-th order Raviart-Thomas (RT) projection. In addition, for For any We refer [10] for the proof of above result. The next result is from [25,Proposition 2.4].
Here C is independent of the mesh size, ν, ν m and κ.
Proof. We recall that Due to discrete trace inequalities, we get that

Proof of L 3 -norm control of discrete magnetic field: Homogeneous tangential components boundary condition
In this section, we prove Theorem 3.1 which provides L 3 -norm control of discrete magnetic field b h . We begin by the following result: Proof. We first recall the inverse inequality for discrete functions [?, Lemma 4.5.3 ]. For any

By (5.2), we have that
, the last step in the above estimate is due to the Sobolev inequality. Consequently, we have Thus we have that Next we present an intermediate result: Proof. Due to [24,Proposition 4.5 According to (5.2) and (5.5a), we have that Next, by the definition of b h , σ h we notice that The last inequality above is due to (5.5b). Finally the proof is complete by combining (5.6), Now we are ready to prove Theorem 3.1.
. We consider the auxiliary Poisson equation: Find φ satisfying: On a polygonal domain Ω, we have the regularity result [26]: there exists δ 0 > 0 such that Let φ h be the numerical solution of (5.9) in the Lagrange space H 1 0 (Ω) ∩ S h . i.e. It solves the system: The last step is due the fact that ∇φ h ∈ C h ∩ H 0 (curl; Ω). Next we present a bound for ∇φ h L 3 (Ω) . To this end, by triangle inequality, we have

by (5.2) and (5.4), we further have
Here we used the approximation property of φ h and the regularity property (5.10). Finally, the proof is complete by combining above estimate with (5.12).
6. Proof of L 3 -norm control of discrete magnetic field: homogeneous normal component boundary condition In this section, we will present detailed proof for Theorem 3.2 which plays a crucial role for the analysis for the second type of boundary condition. Analogous to the proof of Theorem 3.1. We first derive the estimate for functions in C h which are discretely divergence free: . Similar as for the Lemma 5.2, we begin by the following Lemma 6.2, which is similar to [31,Lemma 3.6]. We provide the proof of Lemma 6.2 in Appendix A.
We also need Lemma 6.3, which is similar to [24,Proposition 4.5]. The proof is presented in Appendix B.
. With above two lemmas, the proof of Lemma 6.1 is almost the same as that of Lemma 5.2. We only need to use Lemma 6.2 and Lemma 6.3 to replace Lemma 5.1 and [24, Proposition 4.5] in the proof of Theorem 3.1. Finally, the proof of Theorem 3.2 can be obtained by replacing the auxillary Poisson problem (5.9) with homogeneous Neumann boundary condition in the proofs in the Theorem 3.1. According to [14,Corollary 3.9], there exists δ 0 > 0 such that So, the proof of Theorem 3.2 can be carried out in the same way as that of Theorem 3.1.

Proof of the existence, uniqueness and boundedness of the approximate solution
In this section, we prove Theorem 3.3 on the existence, uniqueness and boundedness of the approximate solution of the DG method. We skip the proof of Theorem 3.4 since it is almost the same as that of Theorem 3.3 (we only need to use Theorem 3.2 to replace Theorem 3.1). The counterpart for the exact solution was provided in [21]. We first define a mapping F on 1) We will show that the mapping is a contraction on a subset of Z h and apply the Brower fixed point theorem for the existence of the solution. Finally the uniqueness follows easily.
Step 1: Definition of the operator F. We start by defining F.
The above system is the original mixed DG scheme for the linearized MHD equations in [25], we refer [25] for the existence and uniqueness of the solutions.
Step 2: Proof of the upper bound of the approximate solution. Next, we establish the boundedness result of the mapping F.

4)
On the other hand, by Lemma 4.2 we get that . In addition, since δ b satisfies (3.3), by Theorem 3.1 and a simple scaling argument, we have Combining the above inequality with (7.3)-(7.5), we get that C . By virtue of (3.6a), F is a contraction on K h if the following quantities are all small enough such that for some constant ρ ∈ [0, 1). As a consequence, by the Brower's fixed point theorem, F has a unique fixed point in K h which is a solution of (2.1). The uniqueness of the solution is trivial since if u h , b h , p h , r h is a solution of (2.1), by (3.6a) u h , b h must be a fixed point of F in K h which is unique.

Proof of the error estimates
In this section, we prove the error estimates of Theorem 3.5. We skip the proof of Theorem 3.6 since it is almost the same as that of Theorem 3.5 (we only need to use Theorem 3.2 and the projection Π C I defined in (4.2) to replace Theorem 3.1 and Π C ). To do that, we proceed in the following steps to give estimates of the projection of the approximation errors, where Π RT is the k-th order Raviart-Thomas (RT) projection, Π C is defined in (4.2), Π Q is the L 2 -orthogonal projection onto Q h and Π S is the Lagrange interpolation onto H 1 0 (Ω) ∩ S h . Since ∇ · u = 0, then ∇ · Π RT u = 0 due to the communication between Π RT and the divergence operator. Thus Π RT u ∈ V h such that e u ∈ V h . Furthermore, due to (4.3a) and (2.1d), e b satisfies (3.3). So we have: Step 1: The error equations. We start our error analysis by obtaining the equations satisfied by the projections of the errors.
Lemma 8.1. The projection of the error (e u , e b , e p , e r ) satisfies Proof. Notice that the exact solution (u, b, p, r) of the equations (1.1) satisfies By the definition of Π RT , Π Q and Π S , we have that for any (v, c, q, r) Subtracting the above equations by (2.1) gives the result.
Step 2: The energy identity derived from the error equations. Now we derive the energy identity which states as follows: We have the energy identity: :=T 1 + · · · + T 9 .

2) and adding all equations, we get that
we denote by Then we get that By (4.1e), we have that Thus we obtain (8.3).
Step 3: Lower bound of the left hand side of the energy identity. Now we provide the lower bound of the left hand side of the energy identity (8.3).

Lemma 8.3.
There is a positive constant γ 0 independent of mesh size such that Proof. By Lemma 4.1, we get that O h (P(u h , { {u h } }); e u , e u ) ≥ 0. By (8.1) and Lemma 4.2, we get that Thus we obtain (8.5).
Step 4: Upper bound of the right hand side of the energy identity. Now we provide upper bound of the right hand side of the energy identity (8.3). We denote by Π V the L 2orthogonal projection onto V h . Since V h = C h , Π V is also the L 2 -orthogonal projection onto C h . We bound T 1 − T 9 as follows: For T 1 , T 2 , we apply Cauchy-Schwarz inequality with the approximation properties of the projections to have: Similarly, for T 4 we have: The last inequality above is due to (4.3c).
With respect to the term T 5 , we chooseẽ r ∈ H 1 0 (Ω) ∩ S h (see [ By the definition of D h , the fact that ∇ · b = 0 and (4.3a), we have that On the other hand, The last inequality above is due to (4.3b). Thus we get that For the convection term T 6 , according to [10, Proposition 4.2], we get that The last inequality is due to (4.1c) and the approximation property of Π RT . The estimates for the last three terms are more involved which we state as following lemmas.
Proof. We recall that Obviously, T ∂K∩F I h . By scaling argument, it is easy to see that Π V u L 6 (Ω) ≤ C u L 6 (Ω) ≤ C u H 1 (Ω) . So by discrete trace inequality and (8.1), On the other hand, by scaling argument and discrete trace inequality to e b , we get that We utilized (8.1) in the last inequality above. Thus we get that T 81 ≤ Cκ u H 1 (Ω) e b 2 C . By generalized Hölder's inequality, Sobolev imbedding, [17,Theorem 5.3] and (8.1), we get that Thus we obtain (8.8).
Proof. We recall that By discrete trace inequality and (4.3b), By (4.3c), it is east to see that
Finally, the estimate of u − u h , b − b h in Theorem 3.5 can be obtained by combining the estimates for T 1 −T 9 together with Theorem 3.3 and the assumption that 1 min(ν, ν m ) u H 1 (Ω) and 1 √ νκν m ∇ × b L 2 (Ω) are small enough.
Step 5: Estimates for r − r h . By a triangle inequality we have with the approximation property of the interpolant we have Therefore it suffice to bound e r S . To this end, by (8.6), consideringẽ r ∈ H 1 0 (Ω) ∩ S h we have The last term is in the lower bound in Lemma 8.5 which has the same estimate as u − u h in Theorem 3.5. For ∇e r L 2 (T h ) , taking c = ∇ẽ r in the error equation (8.2b) we have: By the definition of D h (·, ·) and the fact thatẽ r ∈ H 1 0 (Ω) ∩ S h we have: On the other hand, due to the fact thatẽ r ∈ H 1 0 (Ω) ∩ S h , with the definition of C h (·; ·, ·) we have Inserting above two identities into (8.11), rearranging terms we have For the rest terms we apply Cauchy Schwarz inequality to obtain: the last step is due to a discrete trace inequality and approximation property of Π S . Finally the estimate is complete by (8.10), Lemma 8.5 and estimates for u − u h , b − b h .
Step 6: Estimates for p − p h . For the pressure we apply a standard inf − sup argument, see [4,25,12]. Since e p ∈ L 2 (Ω), there exists a w ∈ H 1 0 (Ω : R 3 ) satisfies It suffices to estimate the term on the right hand side. To this end, let Π BDM w denote the BDM projection of w into v h ∩ H(div; Ω). Due to the orthogonal property of the BDM projection, we have Taking v = Π BDM w in the error equation (8.2a), after rearranging terms, we arrive at: For Y 1 we simply apply Cauchy-Schwarz inequality and discrete trace inequality to have For Y 2 , with a similar estimate as T 6 we can have For Y 3 , we can apply similar estimates as for T 7 − T 9 in Lemma 8.7 -Lemma 8.9 to have Finally, it is well-known that we have If we combine the estimates (8.12) -(8.16), we have ). Finally, the estimate for p − p h is complete by the above estimate together with estimates for u − u h , b − b h and Theorem 3.3.

A divergence-free HDG method for MHD
In this section, we present the hybridizable discontinuous Galerkin (HDG) method for the MHD problem (1.1). There are two main advantages comparing with the mixed DG method proposed in the previous sections: it provides exactly divergence-free velocity field. This feature makes HDG methods more robust in the sense that the convergence of the velocity and magnetic fields is independent of the pressure. Secondly, like all existing HDG methods for different problems, the only global unknowns for the system are the traces of some of the interior unknowns. The hybridization can significantly reduce the size of the global system.
To present the HDG scheme, we first write the original problem (1.1) into a system of first order equations by introducing two more auxiliary unknowns L, w: can be treated and analyzed in a similar manner.
in the interior of each element, we also need to introduce new unknowns (L h , w h , λ h , u h , b Here we define the space P k (∂T h ) := {w ∈ L 2 (∂T h )|w| F ∈ P k (F ), for all F ∈ ∂K, K ∈ T h }. This space, together with the Lagrange multiplier λ h was introduced in [11,19] for HDG methods with exact divergence-free velocity fields. It is worth to mention that on each interior face F shared by two elements, λ h is "double-valued" and eventually it approximates 0. Finally, we adopt the standard integral notation for HDG methods [12], for scalar-valued functions φ and ψ, we write The above notation also applies for vector-valued functions.
Following the standard procedure of devising HDG methods [12], we first first multiplying the test functions (G, v, z, c, q, s) 1a) -(9.1f) and applying integrating by parts to have: Next we replace the exact solutions with the numerical approximations. The HDG method then seeks the approximation magnetic fields, we adopt the choices used in [9]: With the above fluxes, the HDG scheme is complete. Motivated by the work in [6,11,19], we define the convection flux to be We recall that λ h ∈ Λ h is a Lagrange multiplier which approximates 0 on ∂T h . Finally, the constraint equation due to the Lagrange multiplier is given by: This completes the HDG scheme for the problem. Hence, we can formulate our HDG method as to seek together with the transmission conditions: On each face F ∈ ∂K for each K ∈ T h , the stablization parameter S u is defined as: (see [6,36]) S u = max{β · n, 0}.
Finally, we set the vector fields β, d as: Remark 9.1. From (9.2e) we can see that u h ∈ H(div; Ω). Further, since ∇ · u h ∈ Q h we can take q = ∇ · u h in (9.2g) to conclude that u h is exactly divergence-free.
Remark 9.2. In practice, to solve the numerical solution of the above nonlinear system we need to apply Picard iteration. Namely, each iteration we need to solve a linearized system by setting (β, d) is the corresponding solution from the previous iteration. Thanks to the exactly divergence-free feature, we don't need to construct a divergence-free convection field β as for DG scheme in (2.3).
For the purpose of presenting the error estimates, we first rewrite the scheme into a more compact form. To this end, we group (9.2a), (9.2b), (9.2g) and (9.2h) together; group (9.2c), (9.2d) and (9.2i) together; group (9.2f) and (9.2j) together, we can rewrite the system as: Here the operators are defined as: 9.1. Well-posedness and hybridization. The well-posedness of the scheme (9.2) can be validated by two steps: first we show that the linearized scheme (β, d given data) is wellposed. Then we can apply a similar Brower fixed point argument as in Section 7 for the DG scheme. Here we only present the first step in detail. The second step is very similar as in Section 7 and we leave it for readers. We define the norms used in the analysis: Theorem 9.1. If (β, d) in the system (9.2) are given data and they satisfy β ∈ H(div 0 ; Ω), d ∈ L 2 (Ω; R 3 ) and d| ∂T h ∈ L 2 (∂T h ; R 3 ), then the linear system (9.2) has a unique solution In addition, we have u h ∈ H(div 0 ; Ω).
Proof. With (β, d) given, the system (9.2) is a square linear system. Therefore, it is suffice to show that we only have zero solution if the right hand side (f , g) = (0, 0). Let t h , r h ) be a solution, next we show that all components vanish under this assumption.
By Remark 9.1 we know that u h ∈ H(div 0 , Ω) and (u h − u h ) · n = 0 on F h . Therefore, by Theorem 2.1 in [19] we have 2) and adding all the equations, after some algebraic manipulation of the terms, we obtain the energy identity as: This implies that Together with (9.6) we have Next, by (9.7) and (9.2c) we have b h ∈ H 0 (curl 0 ; Ω). The (9.2f), (9.2j) reduces to for all (s, s) ∈ S h × N h . This implies that b h ∈ H 0 (div 0 , Ω). Therefore, by Lemma 5.1 we have b h = 0, and then b t h = 0 by (9.7). Now the equation (9.2d) reduces to: (∇r h , c) T h = 0, for all c ∈ C h , Simply taking c = ∇r h implies that ∇r h = 0. This shows that r h is piecewise constant over T h . By the fact that r h = r h on F h and r h vanishes on ∂Ω we conclude that Finally, we need to show that p h , λ h also vanish. To this end, we use the equations (9.2b), (9.2h) which reduce to: Since p h ∈ L 2 0 (Ω), there exists a γ ∈ H(div; Ω) such that ∇ · γ = p h in Ω, γ · n = 0 on ∂Ω.
Taking v = Π BDM γ in (9.8a) we have by the property of the BDM-projection and (9.8b) we have p h = 0. Now (9.8a) becomes: On each element K we can always find a v ∈ P k (K; R 3 ) such that v · n| ∂K = λ h since this is part of the degree of freedoms of BDM element of degree k. With this v in the above equation we can conclude that λ h = 0. This completes the proof.
Like all HDG methods, the above scheme can be hybridized so that the only globally coupled unknowns are ( u h , b t h , r h ,p h ) wherep h approximates the average pressure within each element. To be more specific, let us decompose the space In practice, we use (9.2a) -(9.2g) as local solvers so that L h , u h , p ⊥ h , b h , w h , λ h are elinimated from the glocal system. The global unknowns are ( u h , b t h , r h ,p h ) coupled by the transmission conditions (9.2h) -(9.2j). In other word, the size of the global system is the dimension of the space M h × M T h × N h ×Q h . The hybridized system reduces the global degrees of freedom significantly, which in turn to make the method more efficient and compatetive with mixed DG and conforming mixed methods. For more details on hybridization of the method we refer to [13,19] and references therein. 9.2. Error equations. In this section, we will derive the error equations based on (9.4). We begin by some notations that will be used in the analysis. For a generic unknown U with its numerical counterpart U h we can split the error as: Here ΠU is some projection/interpolant of U into the discrete polynomial space which we specify next. For (u, p, c, r) we use the same projections (Π RT , Π Q , Π C , Π S ) as for the mixed DG method in Section 8. For the unknowns (L, w) we use the standard L 2 −projections (Π G L, Π W w). For the boundary unknowns (u| F h , b t | F h , r| F h ) we simply use the traces of the projections so that on F h we define: Now we are ready to derive the error equations for the method. Notice that the exact solution (u, L, p, b, w, r, 0, u| F h , b t | F h , r| F h ) also satisfies the discrete system (9.4), if we subtract these two systems and use the splitting discussed above, we can obtain the error equations as: (c, c t )) − N h (c; (e r , e r )) (9.9b) 9.2.1. Energy identity. Similar as the analysis for the mixed DG method in Section 8, we have the following energy identity and several auxiliary estimates for the final error estimates: Lemma 9.2. The projection of the errors satisfy: (1) e u ∈ H(div 0 ; Ω) ∩ V h and (e u − e u ) · n = 0 on ∂T h .
(2) we have (e u , e u ) 1,h ≤ C e L L 2 (Ω) . (3) (e b , ∇s) T h = 0 for all s ∈ H 1 0 (Ω) ∩ S h , and it holds: In addition, we have the energy identity: Proof. (1) is direct consequence of the fact that u h ∈ H(div 0 ; Ω) and (u h − u h )·n = 0 on ∂T h and the conforming property of the Raviart-Thomas projection. (2) is a direct consequence of Theorem 2.1 in [19] and (1) in this lemma.
For (3), in (9.4d) if we restrict s ∈ H 1 0 (Ω) ∩ S h and s = s on F h , we have (b h , ∇s) T h = 0, ∀s ∈ H 1 0 (Ω) ∩ S h , then we also have (e b , ∇s) T h = 0, ∀s ∈ H 1 0 (Ω) ∩ S h , by the property of projection Π C in Lemma 4.3. Further, this means that e b satisfies the condition in Theorem 3.1, therefore, we have The last inequality is due to the fact that e b,t is single valued on ∂T h and a triangle inequality.
We only need to show that the first two terms vanish in the last equality. By (1), we have D h (δ p ; (e u , e u )) = 0 and D h (e p ; (δ u , δ u )) = (∇ · δ u , e p ) T h = 0, by the commuting property: ∇ · Π RT u = Π Q (∇ · u) and e p ∈ Q h . This completes the proof.

Error estimates.
Comparing energy identity for the HDG method in 9.10 and the one for the mixed DG method in (8.3), we can see the the right hand side of 9.10 is independent of pressure p. This is due to the fact that the HDG method provides exactly divergence-free velocity. Thanks to this feature, the error estimate for the energy norm is independent of the regularity of the pressure. Consequently, for the error estimates, we can further relax the regularity assumption (3.7) to be: for σ, σ m > 1 2 , σ p > 0. Now we are ready to state our main convergence results for the HDG method: Theorem 9.3. Let (u, L, p, b, w, r) be the exact solution of the system (9.1), and (L h , u h , p h , b h , w h , r h , λ h , u h , b t h , r h ) be the solution of the HDG method (9.2). With the same assumption as in Theorem 3.3, in addition with the regularity assumption (9.11) and that 1 min(ν, ν m ) u H 1 (Ω) and 1 √ νκν m ∇ × b L 2 (Ω) are small enough, then we have here C depends on the physical parameters κ, ν, ν m and the external forces f , g but is independent of mesh size h.
The proof of the result is based on Lemma 9.2 and mimicking the proofs for the mixed DG method in Section 8. We leave the details for readers who are interested.

Concluding Remarks
In this paper we rigorously analyzed a DG scheme for MHD problem. With standard regularity assumption on the exact solution, we proved that the numerical solution converges to the exact solution optimally for all unknowns in the energy norm. To the best of our knowledge, it is the first analysis dedicated to DG methods for nonlinear MHD problems. In order to make the method more attractive and competitive, we also derive and analyze the first HDG scheme for the problem with several unique features in addition to those for the mixed DG method, including but not limited to (1) It reduces the size of the global system significantly by the hybridization technique, (2) It provides exactly divergence-free velocity fields, (3) The errors for the velocity and magnetic fields are independent of the regularity of the pressure. The issues related with implementation of the mixed DG and HDG method are subjected to ongoing work.
By scaling argument and (A.2), we have that Again by scaling argument, It is easy to see there is a constant C > 0 such that ∀v ∈ L 3 (Ω).
Appendix B. Proof of Lemma 6.3 In this section, we give the proof for Lemma 6.3, which highly mimics that of [24,Proposition 4.5].
For any element K ∈ T h , we use E(K) and F(K) to denote the sets of edges and faces of K. For any face f ∈ F h , we denote by E(f ) the set of the edges of f . The functions {ϕ i K,e } Ne i=1 , are Lagrange basis functions on P k (K; R 3 ), which are introduced in step 2 of the proof of [24,Proposition 4.5]. In fact, {ϕ i K,e } Ne i=1 are basis functions associated with edge degrees of freedom for any edge edge e ∈ E(K); {ϕ i K,f } N f i=1 are basis functions associated with face degrees of freedom for any face f ∈ F(K); and {ϕ i K,b } N b i=1 are basis functions associated with volume degrees of freedom of K.
Like (A.1) in step 2 of the proof of [24,Proposition 4.5], for any K ∈ T h , we have that For any edge e ∈ E(K), we denote by F(e) the set of the faces sharing the edge e. For f ∈ F(e), we denote by K f and K ′ f the elements sharing the face f . It is easy to see that K f = K ′ f if f is a face on ∂Ω. By the construction ofb h , Cauchy-Schwarz inequality and shape-regularity assumption, we have that Then, by (A.3) in step 4 of the proof of [24,Proposition 4.5], we have that It is easy to see that for any face f ∈ F(K), we have that Combing the above estimates yields Summing over all elements, taking into account the shape-regularity of the mesh, we finish the proof.