# The Virtual Element Method for the 3D Resistive Magnetohydrodynamic Model

Lourenço Beirão da Veiga, Franco Dassi11footnotemark: 1, Gianmarco Manzini,
Lorenzo Mascotto11footnotemark: 1 33footnotemark: 3 Dipartimento di Matematica
e Applicazioni, Università degli Studi di Milano-Bicocca, Italy
<EMAIL_ADDRESS><EMAIL_ADDRESS>lorenzo.mascotto@unimib.it)IMATI-CNR, Pavia, Italy
(marco.manzini@imati.cnr.it)Faculty of Mathematics, University of Vienna, 1090
Vienna, Austria (lorenzo.mascotto@univie.ac.at)Corresponding author

###### Abstract

We present a four-field Virtual Element discretization for the time-dependent
resistive Magnetohydrodynamics equations in three space dimensions, focusing
on the semi-discrete formulation. The proposed method employs general
polyhedral meshes and guarantees velocity and magnetic fields that are
divergence free up to machine precision. We provide a full convergence
analysis under suitable regularity assumptions, which is validated by some
numerical tests.

AMS subject classification: 65N12; 65N15; 76E25; 76W05.

Keywords: Resistive Magnetohydrodynamics; Navier-Stokes equations; Maxwell’s
equations; Virtual Element Method; Polyhedral Meshes.

##  1 Introduction

Magnetohydrodynamics (MHD) is the physical-mathematical framework that
describes the dynamics of magnetic fields in electrically conducting fluids,
e.g., plasmas and ionized gases [16, 17]. Indeed, the electrically charged
particles moving in a plasma generate an electromagnetic field that self-
consistently interacts with the fluid motion through the Lorentz force acting
on such particles. The theoretical framework of MHD has widely been used to
develop predictive mathematical models, for example, in astrophysics (solar
wind and space weather [40, Chapter 1], galactic jets [54]), geophysics
(dynamo theory [55]), and nuclear fusion (design of fusion reactors, ignition
in _inertial confinement fusion_ [32, 53, 42]). Roughly speaking, MHD resorts
on a strongly nonlinear coupling of a fluid flow submodel and an
electromagnetic one. In a general setting, as the one considered in this work,
the fluid flow submodel is based on the incompressible Navier-Stokes equations
and the electromagnetic submodel on the Maxwell equations. A thorough
description of the MHD models, their variants, derivation, and physical and
mathematical properties can be found in many textbooks and review papers,
e.g., [35, 63, 16, 17] just to mention a few.

Despite its ubiquity in scientific applications, which is reflected by the
number of papers and books published in the last few decades, the computer
resolution of the system of the MHD equations is still a formidable and
challenging task and may pose pitfalls to the numerical scientists. Here, we
shortly review some of the more critical points, with particular attention on
the usefulness of polyhedral meshes and related schemes, always referring all
interested readers to the rich technical literature for specific and detailed
expositions.

First, vorticial and/or shear flows may lead to large domain deformations,
i.e., severe mesh deformations, also with the occurence of non-planar mesh
faces, in the Lagrangian and Arbitrary Lagrangian-Eulerian (ALE) numerical
frameworks [48, 67]. Such mesh deformations are treated by remeshing the
computational domain; then, by remapping the unknowns to the new mesh [66]. An
optimal remap algorithm should be accurate, cheap, and feature-preserving,
i.e., it should preserve the divergence-free constraints, the solution
positivity, and conserve mass, momentum, and energy, cf. [21]. However,
current remap algorithms can be too expensive, prone to significant accuracy
losses, and not feature-preserving. Therefore, we infer that a good
approximation method should be accurate and stable even in presence of severe
mesh deformations in order to reduce the number of remeshing steps in a
numerical simulation.

Next, low-order accurate algorithms can be affected by an excessive numerical
diffusion and dispersion. For example, numerical diffusion can smear physical
vorticity and shocks at an unacceptable level and even force the magnetic
reconnection in ideal MHD models, which is clearly an artificial and
physically meaningless effect. Numerical diffusion and dispersion can be
reduced by increasing the order of the approximation and by using carefully
selected polygonal or polyhedral meshes. In [37, 52], the Finite Difference
Time Domain (FDTD) method [74] is applied to a grid of hexagonal prisms and
yields much less numerical dispersion and anisotropy than on regular
hexahedral grids where such method is normally considered.

The final major point that we want to consider is related to the divergence-
free nature of both the fluid velocity and the magnetic flux fields. If such
constraints are not accurately satisfied, possibly at the machine precision
level, unreliable or even unphysical solutions may result from a numerical
simulation. The consequence of the violation of the divergence-free constraint
for the magnetic flux field has widely been investigated in the literature. It
was seen that the numerical simulations are prone to significant errors, see,
e.g., [17, 18, 33, 71], as fictitious forces and an unphysical behavior may
appear, cf. [33]. For this reason, in the last two decades, a great effort has
been devoted to the development of divergence-free numerical approximations.
Just to mention a few possible approaches, the divergence-free constraint of
the magnetic flux field can explicitly be enforced by introducing a Lagrange
multiplier in the set of the unknowns, cf. [36]; by using a special flux
limiter in the formulation of the numerical scheme, cf. [56]; by minimizing a
special energy functional in a least squares finite element formulation, cf.
[51]. However, such fixing strategies can result in costly inefficient
schemes, they can be in conflict with each other or not suitable to polyhedral
meshes. A major breakthrough was provided by realizing that “classical”
numerical discretizations fail to reproduce the divergence-free nature of
these fields if the discrete version of the divergence applied to the discrete
version of the rotational operator does not annihilate at the zero machine
precision. In fact, a small but not zero remainder can significantly
accumulate over the many cycles of a long-time calculation, thus breaking the
constraint. This understanding has led to the design of numerical
approximation techniques in the framework of compatible/mimetic methods, as
for example in [66, 67, 58]. A major drawback of compatible/mimetic
discretizations is that they are only low-order accurate and can be too
dispersive and diffusive. The discontinuous Galerkin (DG) method offers high-
order accurate discretizations suitable to general polyhedral meshes that can
reduce the overall numerical diffusion and dispersion. Nonetheless, DG methods
are not compatible and may require a costly divergence-cleaning procedure
based on solving an additional global equation, thus resulting in very
expensive calculations. Finally, we note that the FEM literature for MHD
models on tetrahedral and hexahedral meshes is very broad. Here, we limit
ourselves to recalling the divergence-free schemes in [49, 44, 47], see the
related work [68] as well, and the partial convergence analysis results from
[45, 64]

In this work, we start exploring and presenting the design of numerical
approximations of the MHD equations based on the Virtual Element Method (VEM).
Our goal is to address at least part of the issues discussed so far, but in a
framework that makes it possible to address all of them in the future
developments of this research project. The VEM, initially proposed in [4] for
the numerical approximation of elliptic problems, is a Galerkin projection
method such as the Finite Element Method (FEM). The major difference between
the VEM and the FEM is that in the formulation and practical implementation of
the VEM we do not need an explicit knowledge of the basis functions that
generate the finite element approximation spaces. In fact, all the bilinear
forms and linear functionals of the discrete variational formulation are built
using suitable polynomial projections that are always computable from a
careful choice of the degrees of freedom. For these reasons, such
approximation spaces and the method itself are dubbed as “virtual”. Local
polynomial consistency and an additional stability term provide the well-
posedness of the final discretization. Since the formulation of the method
does not need a closed form of the basis functions, the resulting
computational framework is extremely powerful and offers important advantages
with respect to the FEM. We can build approximation spaces that are in
principle of any order of accuracy and global regularity, suitable to general
polyhedral meshes, even including nonconvex or nonconforming elements in 2D
and 3D, and, importantly, satifying other additional properties such as being
part of a discrete de Rham chain. This last property has the major consequence
that both the numerical approximation of the fluid velocity and the magnetic
flux field are intrinsically divergence free, a fact that we can verify
numerically up to the machine precision.

The discrete divergence-free property in the VEM has been firstly analyzed in
[6, 7] (see also [3]) in the framework of the Stokes and Navier-Stokes
problems. Later, it has been extended to other fluid flow models, a short list
of representatives being [25, 60, 30, 59, 22, 24, 43, 73, 72, 23, 9, 14, 50,
34]. On the other hand, De-Rahm complexes in the virtual element setting have
been introduced in [12] and then improved (and applied to the analysis of
magnetostatic problems) in [5, 11, 10]. The Maxwell’s equations were the
target of [13], see also [27, 28], and their analysis hinges upon technical
tools from [8]. Other interesting results along these directions are contained
in [29] that paves the way for a VEM discretization of the Bianchi-Einstein
equations.

In the present contribution we present a low-order four-field VEM formulation
of the time-dependent MHD equations in three dimensions, in the spirit of
[49]. In order to discretize the unknown fields (pressure, and velocity,
magnetic, and electric fields), we exploit the compatible VEM discrete spaces
outlined above. Since the novelty of the method is in the space discretization
approach, we focus on the semi-discrete version of the scheme, and leave the
development of different time advancing techniques, together with the
associated linear and nonlinear solvers, to future publications. The ensuing
VEM scheme employs general polyhedral meshes and guarantees that both the
velocity and magnetic fields are divergence free (up to machine precision). An
important achievement of the present contribution is that, under suitable
regularity conditions on the exact solution, we are able to prove the linear
convergence of the scheme in the “natural” norms of the problem. To the best
of our knowledge, there are few convergence analyses of FEM schemes for the
MHD model in the literature, and none for the 4-field formulation here
presented. Thus, a major contribution of this work is that we also provide the
details for the convergence analysis for the fully nonlinear model; such an
analysis can be simplified to the FE setting, thus leading to possibly useful
results also for FEM. Clearly, our analysis has the drawback of requiring a
sufficiently regular exact solution; thus, our theoretical results should be
intended as a guarantee that, at least in good conditions, the method delivers
an accurate solution. In the final part of the article, we present some
numerical tests that show the good performance of the proposed method.

The paper is organized as follows. After discussing some basic notation and
definitions in Section 2, we review the continuous problem in Section 3. In
Section 4, we introduce the velocity/pressure discrete fields and some
theoretical results essential to the approximation of the convection term. In
Section 5, after presenting the discrete spaces adopted for the
electromagnetic fields and some related property, we state the discrete
method. In Section 6, we show and prove the theoretical convergence results,
which we validate in Section 7 with some numerical tests. Finally, in Section
8, we draw some conclusions and outline future developments.

##  2 Notation

###  2.1 Mesh notation and functional spaces

Let $\mathcal{T}=\\{\Omega_{h}\\}_{h}$ be a family of mesh decompositions
$\Omega_{h}$ of the computational domain $\Omega$ uniquely identified by the
value of the mesh size parameter $h$. Every mesh $\Omega_{h}$ is a finite
collection of polytopal elements $K$ forming a finite covering of $\Omega$,
i.e., $\overline{\Omega}=\bigcup_{K\in\Omega_{h}}K$, planar faces
$\mathcal{F}_{h}$, straight edges $\mathcal{E}_{h}$, and vertices
$\mathcal{V}_{h}$. The mesh elements are nonoverlapping in the sense that the
intersection in $\mathbb{R}^{3}$ of any pair of them, e.g., $K$ and
$K^{\prime}$, has Lebesgue measure (volume) equal to zero, i.e., $\left|K\cap
K^{\prime}\right|=0$ if $K\neq K^{\prime}$. Accordingly, the intersection of
the elemental boundaries of $K$ and $K^{\prime}$ is either the empty set, or a
set of common vertices, edges, or faces.

For every element $K$, we denote its volume, center, and diameter by $|K|$,
$\mathbf{b}_{K}=(x_{K},y_{K},z_{K})^{T}$, and
$h_{K}=\max_{\mathbf{x},\mathbf{y}\in K}\left|\mathbf{x}-\mathbf{y}\right|$,
and its set of faces, edges, and vertices by $\mathcal{F}^{K}_{h}$,
$\mathcal{E}^{K}_{h}$, and $\mathcal{V}^{K}_{h}$, respectively. As usual, the
maximum of the diameters $h_{K}$ for $K\in\Omega_{h}$ is the mesh size, i.e.,
$h=\max_{K\in\Omega_{h}}h_{K}$. The boundary of every element $K$ is denoted
by $\partial K$ and formed by a finite set of nonintersecting planar faces
$F\in\mathcal{F}^{K}_{h}$ so that $\partial
K=\bigcup_{F\in\mathcal{F}^{K}_{h}}F$. For every face $F$, we denote its area,
center, and diameter by $|F|$, $\mathbf{b}_{F}=(x_{F},y_{F},z_{F})^{T}$, and
$h_{F}=\max_{\mathbf{x},\mathbf{y}\in F}\left|\mathbf{x}-\mathbf{y}\right|$,
and its set of edges and vertices by $\mathcal{E}^{F}_{h}$ and
$\mathcal{V}^{K}_{h}$, respectively. The boundary of every face $F$ is the
planar polygon $\partial F$, which is formed by a finite number of
nonintersecting straigh segments $e\in\mathcal{E}^{F}_{h}$ so that $\partial
F=\bigcup_{e\in\mathcal{E}^{F}_{h}}e$. We denote the length of edge $e$ by
$h_{e}$ and its center, which we shall also refer to as the edge mid-point, by
$\mathbf{b}_{e}=(x_{e},y_{e},z_{e})$.

We assume that all the meshes $\Omega_{h}$ of a given
sequence$\\{\Omega_{h}\\}_{h}$ satisfy these conditions for $h\to 0$: there
exists a real constant factor $\gamma\in(0,1)$ that is independent of
$\Omega_{h}$ and $h_{K}\in\Omega_{h}$ such that

$\mathbf{({M}_{1}})$

shape-regularity: all the elements $K\in\Omega_{h}$ and faces
$F\in\mathcal{F}_{h}$ are “$\gamma$-shape” regular [31];

$\mathbf{({M}_{2}})$

uniform scaling: $\gamma h_{K}\leq h_{F}$ for every face
$F\in\mathcal{F}^{K}_{h}$ of every element $K\in\Omega_{h}$, and, analogously,
$\gamma h_{F}\leq h_{e}$ for every edge $e\in\mathcal{E}^{F}_{h}$ of every
$F\in\mathcal{F}^{K}_{h}$.

Throughout the paper, we shall refer to
$\mathbf{({M}_{1}})$-$\mathbf{({M}_{2}})$ as the _mesh regularity
assumptions_. Such assumptions could be weakened [19, 20, 15, 26]; however,
for the sake of presentation, we stick to a simpler setting.

On every face $F\in\mathcal{F}_{h}$, we define a local coordinate system
$(\xi_{1},\xi_{2})$ and denote the differentiation along $\xi_{1}$ and
$\xi_{2}$ as $\partial_{\xi_{1}}$ and $\partial_{\xi_{2}}$, respectively. The
corresponding second order derivatives are denoted by
$\partial^{2}_{\xi_{1}\xi_{1}}$ and $\partial^{2}_{\xi_{2}\xi_{2}}$. Then, we
consider the two-dimensional vector-valued field
$\mathbf{v}=(v_{1},v_{2}):F\subseteq\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$
and the scalar field $v:F\subseteq\mathbb{R}^{2}\rightarrow\mathbb{R}$, and
let the face-based differential operators $\operatorname{div_{F}}$,
$\operatorname{rot_{F}}$, $\operatorname{\mathbf{curl}_{F}}$ and
$\operatorname{\Delta_{F}}$ applied to $\mathbf{v}(\xi_{1},\xi_{2})$ and
$v(\xi_{1},\xi_{2})$ be defined as

$\displaystyle\operatorname{div_{F}}\mathbf{v}$
$\displaystyle:=\partial_{\xi_{1}}v_{1}+\partial_{\xi_{2}}v_{2},\qquad\operatorname{rot_{F}}\mathbf{v}:=\partial_{\xi_{2}}v_{1}-\partial_{\xi_{1}}v_{2},$
(1a) $\displaystyle\operatorname{\mathbf{curl}_{F}}v$
$\displaystyle:=\big{(}\partial_{\xi_{2}}v,-\partial_{\xi_{1}}v\big{)}^{T},\qquad\operatorname{\Delta_{F}}v:=\partial^{2}_{\xi_{1}\xi_{1}}v+\partial^{2}_{\xi_{2}\xi_{2}}v.$
(1b)

We denote the partial derivatives along the directions $x$, $y$, and $z$ by
$\partial_{x}$, $\partial_{y}$, and $\partial_{z}$, respectively, and the
corresponding second derivatives by $\partial^{2}_{xx}$, $\partial^{2}_{yy}$,
and $\partial^{2}_{zz}$. We define the Laplace, divergence, and curl operators
of the three-dimensional field $\mathbf{v}=(v_{1},v_{2},v_{3})^{T}$ as
follows:

$\displaystyle\Delta\mathbf{v}:=\partial^{2}_{xx}v_{1}+\partial^{2}_{yy}v_{2}+\partial^{2}_{zz}v_{3},\qquad\operatorname{div}\mathbf{v}:=\partial_{x}v_{1}+\partial_{y}v_{2}+\partial_{z}v_{3},$
$\displaystyle\operatorname{\mathbf{curl}}\mathbf{v}:=\big{(}\partial_{y}v_{3}-\partial_{z}v_{2},\partial_{z}v_{1}-\partial_{x}v_{3},\partial_{x}v_{2}-\partial_{y}v_{1}\big{)}^{T}.$

###  2.2 Functional spaces

Consider the polygonal/polyhedral domain $D\subset\mathbb{R}^{d}$, $d=2,3$,
with (Lipschitz) boundary $\partial D$. Throughout the paper, $D$ can be
either a mesh face $F$, a mesh element $K$, or the whole domain $\Omega$.

##### Sobolev spaces

According to [1], $L^{2}(D)$ denotes the Lebesgue space of real-valued square
integrable functions defined on $D$; $L^{2}_{0}(D)$ is the subspace of the
functions in $L^{2}(D)$ having zero average on $D$; for any integer number
$s>0$, $H^{s}(D)$ is the Sobolev space $s$ of the real-valued functions in
$L^{2}(D)$ with all weak partial derivatives of order up to $s$ in $L^{2}(D)$;
$\big{[}L^{2}(D)\big{]}^{d}$, $\big{[}L^{2}_{0}(D)\big{]}^{d}$, and
$\big{[}H^{s}(D)\big{]}^{d}$ are the vector version of these spaces. The
Sobolev spaces of noninteger orders are constructed by using the interpolation
theory and the Sobolev spaces of negative orders by duality [1].

We denote the Sobolev space on $\partial D$ by $H^{s}(\partial D)$ and recall
that the functions in $H^{s}(D)$ admit a trace on $\partial D$ when $s>1/2$.
Since $D$ can be either a polygonal or a polyhedral domain, the bound $s<3/2$
must also be valid. We denote the inner product in $L^{2}(D)$ and $H^{s}(D)$
by $(v,w)_{D}$ and $(v,w)_{s,D}$, respectively; we also denote the
corresponding induced norms by $\|v\|_{D}$ and $\|v\|_{s,D}$, and the seminorm
in $H^{s}(D)$ by $|v|_{s,D}$. When $D$ is the whole computational domain, we
prefer to omit the subindex $\Omega$ and rather use the notation $(v,w)$,
$\|v\|$, $\|v\|_{s}$, etc. instead of $(v,w)_{\Omega}$, $\|v\|_{\Omega}$,
$\|v\|_{s,\Omega}$, etc.

In the light of these definitions, for a given $s>0$, we introduce the
functional spaces

$\displaystyle H^{s}(\operatorname{div},D)$
$\displaystyle:=\Big{\\{}\mathbf{v}\in\big{[}H^{s}(D)\big{]}^{3}\mid\operatorname{div}\mathbf{v}\in
H^{s}(D)\Big{\\}},$ (2) $\displaystyle H^{s}(\operatorname{\mathbf{curl}},D)$
$\displaystyle:=\Big{\\{}\mathbf{v}\in\big{[}H^{s}(D)\big{]}^{3}\mid\operatorname{\mathbf{curl}}\mathbf{v}\in\big{[}H^{s}(D)\big{]}^{3}\big{\\}}.$
(3)

If $s=0$, we write $H(\operatorname{div},D)$ and
$H(\operatorname{\mathbf{curl}},D)$ instead of $H^{0}(\operatorname{div},D)$
and $H^{0}(\operatorname{\mathbf{curl}},D)$. Let $\mathbf{n}_{D}$ be the unit
vector orthogonal to the boundary $\partial D$ and pointing out of $D$.
According, e.g., to [62, Section $3.5$], for the spaces (2) and (3), we can
define the trace operators

$\displaystyle\operatorname{\operatorname{tr}_{\operatorname{div}}}:H(\operatorname{div},D)\to
H^{-\frac{1}{2}}(\partial
D),\qquad\operatorname{\operatorname{tr}_{\operatorname{\mathbf{curl}}}}:H(\operatorname{\mathbf{curl}},D)\to\big{[}H^{-\frac{1}{2}}(\partial
D)\big{]}^{d},\qquad d=2,3,$

which are such that

$\displaystyle\operatorname{\operatorname{tr}_{\operatorname{div}}}(\mathbf{v}):=\mathbf{n}_{D}\cdot\mathbf{v},\qquad\operatorname{\operatorname{tr}_{\operatorname{\mathbf{curl}}}}(\mathbf{v}):=\mathbf{n}_{D}\times\mathbf{v},$

for all sufficiently smooth vector-valued field $\mathbf{v}$. These trace
operators allow us to define the subspaces of $H^{s}(\operatorname{div},D)$
and $H^{s}(\operatorname{\mathbf{curl}},D)$

$\displaystyle H^{s}_{0}(\operatorname{div},D)$
$\displaystyle:=\big{\\{}\mathbf{v}\in
H^{s}(\operatorname{div},D)\mid\operatorname{\operatorname{tr}_{\operatorname{div}}}\mathbf{v}=0\big{\\}},$
$\displaystyle H^{s}_{0}(\operatorname{\mathbf{curl}},D)$
$\displaystyle:=\big{\\{}\mathbf{v}\in
H^{s}(\operatorname{\mathbf{curl}},D)\mid\operatorname{\operatorname{tr}_{\operatorname{\mathbf{curl}}}}\mathbf{v}=\mathbf{0}\big{\\}}.$

These subspaces incorporate the homogeneous boundary conditions in their
definition. Finally, we define the norms for div and curl spaces: for every
sufficiently smooth field $\mathbf{v}$,

$\|\mathbf{v}\|_{\operatorname{\mathbf{curl}}}^{2}:=\|\mathbf{v}\|_{0}^{2}+\|\operatorname{\mathbf{curl}}\mathbf{v}\|_{0}^{2},\qquad\qquad\|\mathbf{v}\|_{\operatorname{div}}^{2}:=\|\mathbf{v}\|_{0}^{2}+\|\operatorname{div}\mathbf{v}\|_{0}^{2}.$

##### Bochner spaces

Let $T>0$ be a real number and $(X,\|\cdot\|_{X})$ a normed space, where $X$
can be either $L^{2}(\Omega)$ or $H^{s}(\Omega)$, $s\geq 0$. According to
[39], the Bochner space $L^{p}(0,T;X)$ is the space of functions $v$ such that
the sublinear functional

$\displaystyle\|v\|_{L^{p}(0,T;X)}=\begin{cases}\displaystyle\left(\int_{0}^{T}\|v(t)\|_{X}^{p}\,dt\right)^{1/{p}}&1\leq
p<\infty,\\\ \textrm{ess sup}_{t\in[0,T]}\|v(t)\|_{X}&p=\infty,\end{cases}$

is a _finite_ norm for almost every $t\in[0,T]$. According to this notation,
$C(0,T;X)$ is the space of the continuous functions from $[0,T]$ to $X$.

##### Polynomial spaces

We denote the space of polynomials of degree $\ell=0,1$ defined on the element
$K$, the face $F$ and the edge $e$ by $\mathbbm{P}_{\ell}(K)$,
$\mathbbm{P}_{\ell}(F)$, and $\mathbbm{P}_{\ell}(e)$, respectively. We set
$\mathbbm{P}_{-1}(K)=\mathbbm{P}_{-1}(F)=\mathbbm{P}_{-1}(e)=\\{0\\}$. The
space $\mathbbm{P}_{1}(K)$ is the span of the _scaled monomials_ defined as:

$\displaystyle
m_{0}(\mathbf{x})=1,\;m_{1}(\mathbf{x})=\frac{x-x_{K}}{h_{K}},\;m_{2}(\mathbf{x})=\frac{y-y_{K}}{h_{K}},\;m_{3}(\mathbf{x})=\frac{z-z_{K}}{h_{K}}\quad\forall\mathbf{x}=(x,y,z)^{T}\in
K.$

The bases of $\mathbbm{P}_{1}(F)$ and $\mathbbm{P}_{1}(e)$ are defined in a
similar way. We let $\mathbbm{P}_{\ell}(\Omega_{h})$ denote the space of the
piecewise discontinuous polynomials of degree $\ell=0,1$ that are globally
defined on $\Omega$ and such that ${q}_{|K}\in\mathbbm{P}_{1}(K)$ for all
elements $K\in\Omega_{h}$.

##### Orthogonal projections onto polynomial spaces

In the forthcoming discrete formulation, given any $\omega$ either in
$\Omega_{h}$, $\mathcal{F}_{h}$, or $\mathcal{E}_{h}$, we shall use the
polynomial projectors
$\Pi^{0,\omega}_{\ell}:L^{2}(\omega)\rightarrow\mathbbm{P}_{\ell}(\omega)$,
$\ell=0,1$, and
$\Pi^{\nabla,\omega}_{1}:H^{1}(\omega)\rightarrow\mathbbm{P}_{1}(\omega)$
defined on all the mesh objects $\omega$. The operator $\Pi^{0,\omega}_{\ell}$
is the orthogonal projection onto constant ($\ell=0$) and linear $(\ell=1)$
polynomials with respect to the inner product in $L^{2}(\omega)$. The operator
$\Pi^{\nabla,\omega}_{1}$ is the orthogonal projection onto linear polynomials
with respect to the semi-inner product in $H^{1}(\omega)$; we call it the
_elliptic projection_. The elliptic projection $\Pi^{\nabla,\omega}_{\ell}v$
of a function $v\in H^{1}(\omega)$ is the linear polynomial solving the
variational problem

$\displaystyle\big{(}\nabla(\Pi^{\nabla,\omega}_{1}v-v),\nabla
q\big{)}_{\omega}=0\quad\forall
q\in\mathbbm{P}_{1}(\omega)\quad\textrm{and}\quad\int_{\partial\omega}\big{(}\Pi^{\nabla,\omega}_{1}v-v\big{)}=0.$
(4)

The second condition in (4) is considered to ensure the uniqueness of the
projection $\Pi^{\nabla,K}_{1}v$.

With an abuse of notation, we extend these definitions in a component-wise
manner to the multidimensional projection operators
$\Pi^{0,\omega}_{\ell}:\big{[}L^{2}(\omega)\big{]}^{d}\to\big{[}\mathbbm{P}_{\ell}(\omega)\big{]}^{d}$
and
$\Pi^{\nabla,\omega}_{1}:\big{[}H^{1}(\omega)\big{]}^{d}\to\big{[}\mathbbm{P}_{1}(\omega)\big{]}^{d}$,
$d=2,3$. In particular, we shall use the orthogonal projection
$\Pi^{0,K}_{0}\nabla v$ of the gradient of a function $v\in H^{1}(K)$, which
is the constant vector polynomial solving the variational problem:

$\displaystyle\big{(}\Pi^{0,K}_{0}\nabla v-\nabla
v,\mathbf{q}\big{)}_{K}=0\quad\forall\mathbf{q}\in\big{[}\mathbbm{P}_{0}(K)\big{]}^{3}.$

For $\ell=0,1$, we also define the global projection operators
$\Pi^{0}_{\ell}:L^{2}(\Omega)\to\mathbbm{P}_{\ell}(\Omega_{h})$ and
$\Pi^{0}_{\ell}:\big{[}L^{2}(\Omega)\big{]}^{d}\to\big{[}\mathbbm{P}_{\ell}(\Omega_{h})\big{]}^{d}$
as the operators respectively satisfying
${\big{(}\Pi^{0}_{\ell}v\big{)}}_{|K}=\Pi^{0,K}_{\ell}\big{(}{v}_{|K}\big{)}$
and
${\big{(}\Pi^{0}_{\ell}\mathbf{v}\big{)}}_{|K}=\Pi^{0,K}_{\ell}\big{(}{\mathbf{v}}_{|K}\big{)}$
for all mesh elements $K\in\Omega_{h}$.

Similarly, $\Pi^{\nabla}_{1}:H^{1}(\Omega)\to\mathbbm{P}_{1}(\Omega_{h})$ and
$\Pi^{\nabla}_{1}:\big{[}H^{1}(\Omega)\big{]}^{d}\to\big{[}\mathbbm{P}_{1}(\Omega_{h})\big{]}^{d}$
are the global projection operators respectively satisfying
${\big{(}\Pi^{\nabla}_{1}v\big{)}}_{|K}=\Pi^{\nabla,K}_{1}\big{(}{v}_{|K}\big{)}$
and
${\big{(}\Pi^{\nabla}_{1}\mathbf{v}\big{)}}_{|K}=\Pi^{\nabla,K}_{1}\big{(}{\mathbf{v}}_{|K}\big{)}$
for all mesh elements $K\in\Omega_{h}$.

##### Some names for constants

In some of the forthcoming estimates, we shall occasionally write explicit
constants and denote them with different symbols, depending on their meaning.
Notably, we shall use the following notation:

  * •

$C_{S}$ denotes a constant depending on a Sobolev embedding;

  * •

$C_{D}$ denotes a constant depending on the shape of a domain;

  * •

$C_{P}$ denotes a constant appearing in a Poincaré inequality;

  * •

$C_{I}$ denotes a constant depending on an interpolation estimate with respect
to functions in virtual element spaces;

  * •

$C_{\operatorname{inv}}$ denotes a constant depending on an inverse estimate;

  * •

$C_{\operatorname{appr}}$ denotes a constant depending on a polynomial
approximation estimate.

We deem that this notation might help the reader to better follow some steps
in the forthcoming proofs.

Henceforth, we use the letter “$C$” to denote a strictly positive constant
that can take a different value at any occurrence. The constant $C$ is
independent of the mesh size parameter $h$ but may depend on the other
parameters of the differential problem and virtual element discretization such
as the domain shape, the mesh regularity constant $\gamma$, and the coercivity
and continuity constants of the bilinear forms used in the variational
formulations that will be introduced in the next sections.

##  3 The continuous problem

Let $\Omega\subset\mathbb{R}^{3}$ be a polyhedral domain with Lipschitz
continuous boundary $\partial\Omega$. We model the interaction of an
electrically charged, incompressible fluid having velocity $\mathbf{u}$ and
pressure $p$ with the self-consistently generated electric field $\mathbf{E}$
and magnetic flux field $\mathbf{B}$. We denote the viscous Reynolds number by
$Re$, the magnetic Reynolds number by $Re_{m}$, and the Hartman number by $s$.
Furthermore, we define the electric current density

$\mathbf{j}:=\mathbf{E}+\mathbf{u}\times\mathbf{B}.$ (5)

Henceforth, the subscript $t$ as in $\mathbf{u}_{t}$ and $\mathbf{B}_{t}$
denotes the time derivative and $[0,T]$ is the time integration interval for a
given final time $T>0$.

The MHD problem reads as follows: _For all times $t\in(0,T]$, find
$\mathbf{u}$, $p$, $\mathbf{E}$, and $\mathbf{B}$, such that_

$\displaystyle\mathbf{u}_{t}+(\nabla\mathbf{u})\mathbf{u}-Re^{-1}\Delta\mathbf{u}-s\mathbf{j}\times\mathbf{B}+\nabla
p$ $\displaystyle=\mathbf{f}\phantom{0\mathbf{0}}\qquad\text{in }\Omega,$ (6a)
$\displaystyle\mathbf{j}-Re_{m}^{-1}\operatorname{\mathbf{curl}}\mathbf{B}$
$\displaystyle=\mathbf{0}\phantom{0\mathbf{f}}\qquad\text{in }\Omega,$ (6b)
$\displaystyle\mathbf{B}_{t}+\operatorname{\mathbf{curl}}\mathbf{E}$
$\displaystyle=\mathbf{0}\phantom{0\mathbf{f}}\qquad\text{in }\Omega,$ (6c)
$\displaystyle\operatorname{div}\mathbf{B}$
$\displaystyle=0\phantom{\mathbf{0}\mathbf{f}}\qquad\text{in }\Omega,$ (6d)
$\displaystyle\operatorname{div}\mathbf{u}$
$\displaystyle=0\phantom{\mathbf{0}\mathbf{f}}\qquad\text{in }\Omega,$ (6e)

where $(\nabla\mathbf{u})\mathbf{u}=\sum_{i,j}(\partial_{j}u_{i})u_{j}$. The
MHD equations (6) are completed by the set of initial conditions for the
velocity and the magnetic flux fields

$\mathbf{u}(\mathbf{x},0)=\mathbf{u}_{0}(\mathbf{x}),\qquad\mathbf{B}(\mathbf{x},0)=\mathbf{B}_{0}(\mathbf{x}),\qquad\forall\mathbf{x}\in\Omega,$
(7)

and the homogeneous boundary conditions

$\mathbf{u}(\mathbf{x},t)=\mathbf{0},\qquad\mathbf{B}(\mathbf{x},t)\cdot\mathbf{n}_{\Omega}=0,\qquad\mathbf{E}(\mathbf{x},t)\times\mathbf{n}_{\Omega}=\mathbf{0}\qquad\forall\mathbf{x}\in\partial\Omega,\quad\forall
t\in[0,T].$ (8)

Equations (6a) and (6e) describe the hydrodynamic behavior of an electrically
charged, incompressible fluid under the action of an external force
$\mathbf{f}$ and the electromagnetic force $\mathbf{j}\times\mathbf{B}$
multiplied by the Hartmann number, which acts as a coupling coefficient. The
electromagnetic force is self-consistently generated by the electromagnetic
fields $\mathbf{E}$ and $\mathbf{B}$ satisfying equations (6b), (6c), and
(6d). These last three equations describe the electromagnetic submodel in the
magneto-hydrodynamics approximation [16]. The incompressibility of the fluid
velocity and the solenoidal nature of the magnetic field require that both
$\mathbf{u}_{0}$ and $\mathbf{B}_{0}$ in (7) are divergence free.

The weak formulation of problem (6)-(8) reads as follows: _For almost every
$t\in[0,T]$, find $(\mathbf{u}(t),p(t)$,
$\mathbf{E}(t),\mathbf{B}(t))\in\big{[}H^{1}_{0}(\Omega)\big{]}^{3}\times
L^{2}_{0}(\Omega)\times H_{0}(\operatorname{\mathbf{curl}},\Omega)\times
H_{0}(\operatorname{div},\Omega)$ such that_

$\displaystyle(\mathbf{u}_{t},\mathbf{v})+Re^{-1}a(\mathbf{u},\mathbf{v})+b(\mathbf{v},p)+c(\mathbf{u};\mathbf{u},\mathbf{v})-s(\mathbf{j}\times\mathbf{B},\mathbf{v})$
$\displaystyle=(\mathbf{f},\mathbf{v})\phantom{0}\qquad\forall\mathbf{v}\in\big{[}H^{1}_{0}(\Omega)\big{]}^{3},$
(9a)
$\displaystyle(\mathbf{j},\mathbf{F})-Re_{m}^{-1}(\mathbf{B},\operatorname{\mathbf{curl}}\mathbf{F})$
$\displaystyle=0\phantom{(\mathbf{f},\mathbf{v})}\qquad\forall\mathbf{F}\in
H_{0}(\operatorname{\mathbf{curl}},\Omega),$ (9b)
$\displaystyle(\mathbf{B}_{t},\mathbf{C})+(\operatorname{\mathbf{curl}}\mathbf{E},\mathbf{C})$
$\displaystyle=0\phantom{(\mathbf{f},\mathbf{v})}\qquad\forall\mathbf{C}\in
H_{0}(\operatorname{div},\Omega),$ (9c) $\displaystyle b(\mathbf{u},q)$
$\displaystyle=0\phantom{(\mathbf{f},\mathbf{v})}\qquad\forall q\in
L_{0}(\Omega),$ (9d)

where $\mathbf{j}$ is defined in (5). The bilinear forms $a(\cdot,\cdot)$ and
$b(\cdot,\cdot)$, and the trilinear form $c(\cdot;\cdot,\cdot)$ are defined as

$\displaystyle
a:\big{[}H^{1}(\Omega)\big{]}^{3}\times\big{[}H^{1}(\Omega)\big{]}^{3}\to\mathbb{R}:$
$\displaystyle\quad
a(\mathbf{u},\mathbf{v}):=(\nabla\mathbf{u},\nabla\mathbf{v})=\displaystyle\int_{\Omega}\nabla\mathbf{u}:\nabla\mathbf{v},$
$\displaystyle b:\big{[}H^{1}(\Omega)\big{]}^{3}\times
L^{2}_{0}(\Omega)\to\mathbb{R}:$ $\displaystyle\quad
b(\mathbf{v},q):=-(\operatorname{div}\mathbf{v},q)=\displaystyle-\int_{\Omega}q\operatorname{div}\mathbf{v},$
$\displaystyle
c:\big{[}H^{1}(\Omega)\big{]}^{3}\times\big{[}H^{1}(\Omega)\big{]}^{3}\times\big{[}H^{1}(\Omega)\big{]}^{3}\to\mathbb{R}:$
$\displaystyle\quad
c(\mathbf{w};\mathbf{u},\mathbf{v}):=((\nabla\mathbf{u})\mathbf{w},\mathbf{v})=\displaystyle\int_{\Omega}(\nabla\mathbf{u})\mathbf{w}\cdot\mathbf{v},$

where $\mathbf{A}:\mathbf{B}=\sum_{ij}A_{ij}B_{ij}$ and
$(\mathbf{A})\mathbf{b}\cdot\mathbf{c}=\sum_{ij}A_{ij}b_{j}c_{i}$ for the
matrices $\mathbf{A}=(A_{ij})$ and $\mathbf{B}=(B_{ij})$, and the vectors
$\mathbf{b}=(b_{i})$ and $\mathbf{c}=(c_{i})$. We also consider the local
forms that are defined by splitting the above forms on the mesh elements
$K\in\Omega_{h}$:

$\displaystyle a(\mathbf{u},\mathbf{v})$
$\displaystyle=\sum_{K\in\Omega_{h}}a^{K}(\mathbf{u},\mathbf{v})\quad\textrm{and}\quad
a^{K}(\mathbf{u},\mathbf{v})=(\nabla\mathbf{u},\nabla\mathbf{v})_{K},$
$\displaystyle b(\mathbf{v},q)$
$\displaystyle=\sum_{K\in\Omega_{h}}b^{K}(\mathbf{u},q)\quad\textrm{and}\quad
b^{K}(\mathbf{u},q)=-(\operatorname{div}\mathbf{v},q)_{K},$ $\displaystyle
c(\mathbf{v};\mathbf{u},\mathbf{w})$
$\displaystyle=\sum_{K\in\Omega_{h}}c^{K}(\mathbf{v};\mathbf{u},\mathbf{w})\quad\textrm{and}\quad
c^{K}(\mathbf{v};\mathbf{u},\mathbf{w})=((\nabla\mathbf{u})\mathbf{v},\mathbf{w})_{K}.$

Several (partial) results are available from the technical literature about
the well-posedness of the MHD model in the strong and weak formulation, e.g.,
problems (6) and (9) and their variants, cf. [38, 70, 61, 65, 46, 41] and the
citations therein. Some of these results concerning existence, uniqueness and
stability of the solution fields $(\mathbf{u},p,\mathbf{E},\mathbf{B})$ have
been derived under specific assumptions that may reduce the generality of the
model. To the best of our knowledge, the full mathematical understanding of
the MHD model is still an open issue and an active research area. This topic
is beyond the scope of our work; thus, we shall simply assume that the MHD
model is well-posed, at least in the setting that we are using in this paper.

##  4 Virtual element spaces for the incompressible flow equations

###  4.1 The discrete pressure space

We approximate the pressure unknown in the space of piecewise discontinuous,
constant functions with zero average on $\Omega$:

$\displaystyle Q_{h}:=\Big{\\{}q_{h}\in
L^{2}_{0}(\Omega)\,\mid\,{q_{h}}_{|K}\in\mathbbm{P}_{0}(K)\;\forall
K\in\Omega_{h}\Big{\\}}=\mathbbm{P}_{0}(\Omega_{h})\cap L^{2}_{0}(\Omega).$

Every $q_{h}\in Q_{h}$ is uniquely determined by the set of constant values
$\big{(}{q_{h}}_{|K}\big{)}$. Accordingly, we can approximate any scalar
function $q\in L^{2}_{0}(\Omega)$ by its averages over the mesh elements.

###  4.2 The discrete velocity space

Here, we consider the “lowest order” version of the spaces introduced in [7,
14]. This includes an “enhancement” procedure that allows for the computation
of an $L^{2}$ linear polynomial projector. We refer to such papers for a
better understanding of the motivations behind this construction.

We define the nodal velocity space on a face $F\in\mathcal{F}_{h}$ as

$\displaystyle\mathbf{W}_{h}(F):=\Big{\\{}\mathbf{w}_{h}\in\big{[}H^{1}(F)\big{]}^{3}\,\mid\,$
$\displaystyle{\mathbf{w}_{h}}_{|\partial F}\in\big{[}C^{0}(\partial
F)\big{]}^{3},\;{\mathbf{w}_{h}}_{|e}\in\big{[}\mathbbm{P}_{1}(e)\big{]}^{3}\,\forall
e\in\mathcal{E}^{F}_{h},\quad$
$\displaystyle\operatorname{\Delta_{F}}\mathbf{w}_{h}\in\big{[}\mathbbm{P}_{2}(F)\big{]}^{3},$
$\displaystyle\big{(}\mathbf{w}_{h,\bm{\tau}}-\bm{\Pi}^{\nabla,F}_{1}\mathbf{w}_{h,\bm{\tau}},\mathbf{q}\big{)}_{0,F}=0\quad\forall\mathbf{q}\in\big{[}\mathbbm{P}_{2}(F)\big{]}^{2},$
$\displaystyle\big{(}w_{h,\mathbf{n}}-\Pi^{\nabla,F}_{1}w_{h,\mathbf{n}},q\big{)}_{0,F}=0\quad\forall
q\in\mathbbm{P}_{2}(F)\setminus\mathbb{R}\Big{\\}}.$

Above, we denoted the tangential and normal components of a given field
$\mathbf{w}$ by
$\mathbf{w}_{h,\bm{\tau}}=(\mathbf{n}_{F}\times\mathbf{w})\times\mathbf{n}_{F}$
and $w_{h,\mathbf{n}}=\mathbf{n}_{F}\cdot\mathbf{w}$, respectively. We define
the nodal velocity space on the whole boundary $\partial K$ of a mesh element
$K$ as

$\mathbf{W}_{h}(\partial K):=\Big{\\{}\mathbf{w}_{h}\in\big{[}C^{0}(\partial
K)\big{]}^{3}\,\mid\,{\mathbf{w}_{h}}_{|F}\in\mathbf{W}_{h}(F)\quad\forall
F\in\mathcal{F}^{K}_{h}\Big{\\}}.$

Finally, we introduce the nodal velocity space on the element $K\in\Omega_{h}$
as follows:

$\displaystyle\mathbf{W}_{h}(K):=\Big{\\{}\,\mathbf{w}_{h}\in\big{[}H^{1}(K)\big{]}^{3}\,\mid\,$
$\displaystyle{\mathbf{w}_{h}}_{|\partial K}\in\mathbf{W}_{h}(\partial
K),\quad$ $\displaystyle\begin{cases}\Delta\mathbf{w}_{h}+\nabla
s\in\mathbf{x}\times\big{[}\mathbbm{P}_{0}(K)\big{]}^{3}\text{\leavevmode\nobreak\
for\leavevmode\nobreak\ some\leavevmode\nobreak\ }s\in L^{2}_{0}(K),\\\
\operatorname{div}\mathbf{w}_{h}\in\mathbbm{P}_{0}(K),\end{cases}$
$\displaystyle\big{(}\mathbf{w}_{h}-\bm{\Pi}^{\nabla,K}_{1}\mathbf{w}_{h},\mathbf{x}\times\mathbf{q}\big{)}_{0,K}=0\quad\forall\mathbf{q}\in\big{[}\mathbbm{P}_{0}(K)\big{]}^{3}\,\Big{\\}}.$

Every virtual element vector field $\mathbf{w}_{h}\in\mathbf{W}_{h}(K)$ is
uniquely determined by the set of values
$\big{(}\,(\mathbf{w}_{\textsf{v}})_{\textsf{v}\in\mathcal{V}^{K}_{h}},(w_{F})_{F\in\mathcal{F}^{K}_{h}}\,\big{)}$,
where

  * •

$\mathbf{w}_{\textsf{v}}=\mathbf{w}_{h}(\mathbf{x}_{\textsf{v}})$ is the value
of $\mathbf{w}_{h}$ at the vertex $\textsf{v}\in\mathcal{V}^{K}_{h}$ (with
position vector $\mathbf{x}_{\textsf{v}}$);

  * •

$w_{F}$ is the scaled zero-th order moment of the normal component of
$\mathbf{w}_{h}$ on the face $F\in\mathcal{F}^{K}_{h}$, given by

$\displaystyle w_{F}=\frac{1}{|F|}\int_{F}w_{h,\mathbf{n}}.$

The unisolvence of these degrees of freedom can be proved as in [3, 14]. For
any $\mathbf{w}_{h}\in\mathbf{W}_{h}(K)$, the elliptic projection
$\Pi^{\nabla,K}_{1}\mathbf{w}_{h}\in\big{[}\mathbbm{P}_{1}(K)\big{]}^{3}$, the
$L^{2}$-orthogonal projection
$\bm{\Pi}^{0,K}_{1}\mathbf{w}_{h}\in\big{[}\mathbbm{P}_{1}(K)\big{]}^{3}$, and
the divergence of $\mathbf{w}_{h}$ are directly computable from such degrees
of freedom. Finally, we define the global nodal velocity space by an
$H^{1}$-conforming coupling of the local degrees of freedom:

$\displaystyle\mathbf{W}_{h}:=\Big{\\{}\,\mathbf{w}_{h}\in\big{[}H^{1}_{0}(\Omega)\big{]}^{3}\,\mid\,{\mathbf{w}_{h}}_{|K}\in\mathbf{W}_{h}(K)\quad\forall
K\in\Omega_{h}\,\Big{\\}}.$

######  Remark 4.1.

As in [7, 14], the couple of spaces $\mathbf{W}_{h}\times Q_{h}$ satisfies a
discrete inf-sup condition, i.e., there exists $\beta>0$ independent of $h$
such that

$\beta\leq\inf_{q_{h}\in
Q_{h}}\sup_{\mathbf{v}_{h}\in\mathbf{W}_{h}}\frac{b(\mathbf{v}_{h},q_{h})}{\|\mathbf{v}_{h}\|_{1}\|q_{h}\|}.$

Furthermore, we also have the important property

$\operatorname{div}\mathbf{W}_{h}\subseteq Q_{h}.$

###  4.3 The virtual element bilinear forms $a_{h}(\cdot,\cdot)$ and
$m_{h}(\cdot,\cdot)$

We use the elliptic projection operator $\Pi^{\nabla,K}_{1}$ to define the
elemental bilinear form
$a_{h}^{K}:\mathbf{W}_{h}(K)\times\mathbf{W}_{h}(K)\to\mathbb{R}$, which
mimics the $H^{1}$-inner product on the element $K$:

$a_{h}^{K}(\mathbf{u}_{h},\mathbf{v}_{h}):=\big{(}\nabla\Pi^{\nabla,K}_{1}\mathbf{u}_{h},\nabla\Pi^{\nabla,K}_{1}\mathbf{v}_{h}\big{)}_{K}+S^{K}\big{(}(I-\Pi^{\nabla,K}_{1})\mathbf{u}_{h},(I-\Pi^{\nabla,K}_{1})\mathbf{v}_{h}\big{)}\\\
\quad\forall\mathbf{u}_{h},\,\mathbf{v}_{h}\in\mathbf{W}_{h}(K).$ (10)

Here, $S^{K}:\mathbf{W}_{h}(K)\times\mathbf{W}_{h}(K)\to\mathbb{R}$ can be any
computable symmetric bilinear form such that there exist two positive
constants $\sigma_{*}$ and $\sigma^{*}$ independent of $h_{K}$ satisfying

$\displaystyle\sigma_{*}|\mathbf{v}_{h}|^{2}_{1,K}\leq
S^{K}(\mathbf{v}_{h},\mathbf{v}_{h})\leq\sigma^{*}|\mathbf{v}_{h}|^{2}_{1,K}\quad\quad\forall\mathbf{v}_{h}\in\ker\big{(}\Pi^{\nabla,K}_{1}\big{)}\cap\mathbf{W}_{h}(K).$

A “standard” choice for the stabilization is given by [14]

$\displaystyle
S^{K}(\mathbf{u}_{h},\mathbf{v}_{h}):=h_{K}\sum_{i}\operatorname{\operatorname{dof}^{K}_{i}}(\mathbf{u}_{h})\,\operatorname{\operatorname{dof}^{K}_{i}}(\mathbf{v}_{h}),\qquad\forall\mathbf{u}_{h},\,\mathbf{v}_{h}\in\mathbf{W}_{h}(K),$

where the summation on the index $i$ is carried over all the elemental degrees
of freedom and $\operatorname{\operatorname{dof}^{K}_{i}}(\mathbf{w}_{h})$ is
the bounded, linear functional defined on $\mathbf{W}_{h}(K)$ providing the
$i$-th degree of freedom of $\mathbf{w}_{h}\in\mathbf{W}_{h}(K)$. The
summation term in $S^{K}(\cdot,\cdot)$ is multiplied by $h_{K}$ to have a
consistent scaling for both terms of $a_{h}^{K}(\cdot,\cdot)$ on the right-
hand side of (10) with respect to the element size.

The local bilinear form $a_{h}^{K}(\cdot,\cdot)$ satisfies the two fundamental
properties of consistency and stability:

  * •

consistency: for all $\mathbf{q}\in\big{[}\mathbbm{P}_{1}(K)\big{]}^{3}$ and
$\mathbf{w}_{h}\in\mathbf{W}_{h}(K)$,

$\displaystyle
a_{h}^{K}(\mathbf{q},\mathbf{w}_{h})=a^{K}(\mathbf{q},\mathbf{w}_{h});$ (11)

  * •

stability: for all $\mathbf{w}_{h}\in\mathbf{W}_{h}(K)$,

$\displaystyle\alpha_{*}\|\mathbf{w}_{h}\|_{1,K}^{2}\leq
a_{h}^{K}(\mathbf{w}_{h},\mathbf{w}_{h})\leq\alpha^{*}\|\mathbf{w}_{h}\|_{1,K}^{2},$
(12)

with $\alpha_{*}=\min(1,\sigma_{*})$ and $\alpha^{*}=\max(1,\sigma^{*})$.

Finally, we define the global bilinear form
$a_{h}:\mathbf{W}_{h}\times\mathbf{W}_{h}\to\mathbb{R}$ that we use in the
virtual element formulation of the MHD model:

$\displaystyle
a_{h}(\mathbf{u}_{h},\mathbf{v}_{h}):=\sum_{K\in\Omega_{h}}a_{h}^{K}(\mathbf{u}_{h},\mathbf{v}_{h})\quad\quad\forall\mathbf{u}_{h},\mathbf{v}_{h}\in\mathbf{W}_{h}.$

In what follows, we use the discrete norm

$\|\mathbf{v}_{h}\|_{\mathbf{W}_{h}}^{2}:=a_{h}(\mathbf{v}_{h},\mathbf{v}_{h})\qquad\forall\mathbf{v}_{h}\in\mathbf{W}_{h}.$

We use the $L^{2}$-projection operator $\bm{\Pi}^{0,K}_{1}$ to define the
bilinear form
$m_{h}^{K}:\mathbf{W}_{h}(K)\times\mathbf{W}_{h}(K)\to\mathbb{R}$, which
mimics the local $L^{2}$-inner product on the element $K$:

$\displaystyle
m_{h}^{K}(\mathbf{u}_{h},\mathbf{v}_{h}):=\big{(}\bm{\Pi}^{0,K}_{1}\mathbf{u}_{h},\bm{\Pi}^{0,K}_{1}\mathbf{v}_{h}\big{)}_{K}+\tilde{S}^{K}\Big{(}(I-\bm{\Pi}^{0,K}_{1})\mathbf{u}_{h},(I-\bm{\Pi}^{0,K}_{1})\mathbf{v}_{h}\Big{)}\quad\forall\mathbf{u}_{h},\,\mathbf{v}_{h}\in\mathbf{W}_{h}(K).$
(13)

Here, $\tilde{S}^{K}:\mathbf{W}_{h}(K)\times\mathbf{W}_{h}(K)\to\mathbb{R}$
can be any computable bilinear form such that there exist two positive
constants $\widetilde{\sigma}_{*}$ and $\widetilde{\sigma}^{*}$ independent of
$h_{K}$ satisfying

$\widetilde{\sigma}_{*}\|\mathbf{v}_{h}\|^{2}_{0,K}\leq\tilde{S}^{K}(\mathbf{v}_{h},\mathbf{v}_{h})\leq\widetilde{\sigma}^{*}\|\mathbf{v}_{h}\|^{2}_{0,K}\quad\quad\forall\mathbf{v}_{h}\in\ker\big{(}\bm{\Pi}^{0,K}_{1}\big{)}\cap\mathbf{W}_{h}(K).$

A “standard” choice for the stabilization term $\tilde{S}^{K}(\cdot,\cdot)$ is
given by:

$\displaystyle\tilde{S}^{K}(\mathbf{u}_{h},\mathbf{v}_{h}):=h_{K}^{3}\sum_{i}\operatorname{\operatorname{dof}^{K}_{i}}(\mathbf{u}_{h})\,\operatorname{\operatorname{dof}^{K}_{i}}(\mathbf{v}_{h}),\qquad\forall\mathbf{u}_{h},\,\mathbf{v}_{h}\in\mathbf{W}_{h}(K),$

where the summation on the index $i$ is again on all the elemental degrees of
freedom and $\operatorname{\operatorname{dof}^{K}_{i}}(\cdot)$ is the same
functional used in the definition of the stabilization term
$S^{K}(\cdot,\cdot)$. The summation term in $\tilde{S}^{K}(\cdot,\cdot)$ is
multiplied by $h_{K}^{3}$ to have a consistent scaling for both terms of
$m_{h}^{K}(\cdot,\cdot)$ on the right-hand side of (13) with respect to the
element size.

The local bilinear form $m_{h}^{K}(\cdot,\cdot)$ satisfies the two fundamental
properties of consistency and stability:

  * •

consistency: for all $\mathbf{q}\in\big{[}\mathbbm{P}_{1}(K)\big{]}^{3}$ and
$\mathbf{w}_{h}\in\mathbf{W}_{h}(K)$,

$\displaystyle
m_{h}^{K}(\mathbf{q},\mathbf{w}_{h})=(\mathbf{q},\mathbf{w}_{h})_{0,K};$ (14)

  * •

stability: for all $\mathbf{w}_{h}\in\mathbf{W}_{h}(K)$,

$\displaystyle\mu_{*}\|\mathbf{w}_{h}\|^{2}_{0,K}\leq
m_{h}^{K}(\mathbf{w}_{h},\mathbf{w}_{h})\leq\mu^{*}\|\mathbf{w}_{h}\|^{2}_{0,K},$
(15)

with $\mu_{*}=\min(1,\widetilde{\sigma}_{*})$ and
$\mu^{*}=\max(1,\widetilde{\sigma}^{*})$.

Finally, we define the global inner product
$m_{h}:\mathbf{W}_{h}\times\mathbf{W}_{h}\to\mathbb{R}$ that we use in the
virtual element formulation of the MHD model:

$\displaystyle
m_{h}(\mathbf{u}_{h},\mathbf{v}_{h}):=\sum_{K\in\Omega_{h}}m_{h}^{K}(\mathbf{u}_{h},\mathbf{v}_{h}).$

Associated with $m_{h}(\cdot,\cdot)$, we define the corresponding discrete
norm

$\|\mathbf{v}_{h}\|_{m_{h}}^{2}:=m_{h}(\mathbf{v}_{h},\mathbf{v}_{h}).$

####  4.3.1 Interpolation in $\mathbf{W}_{h}$ and a stability result

We define the (energy) interpolant $\mathbf{u}_{I}\in\mathbf{W}_{h}$ of a
sufficiently smooth vector-valued field $\mathbf{u}$ as the unique function in
$\mathbf{W}_{h}$ sharing the same degrees of freedom of $\mathbf{u}$. We have
the following local and global interpolation property. We omit the proof as it
would be a fairly standard modification of the analogous results in [7].

######  Lemma 4.2.

Let $\mathbf{u}\in\big{[}H^{2}(\Omega)\big{]}^{3}$ and
$\mathbf{u}_{I}\in\mathbf{W}_{h}$ be its degrees of freedom interpolant. Then,
a real, positive constant $C_{I}$ independent of $h$ exists such that

$\|\mathbf{u}-\mathbf{u}_{I}\|_{0,K}+h_{K}|\mathbf{u}-\mathbf{u}_{I}|_{1,K}\leq
C_{I}h_{K}^{2}|\mathbf{u}|_{2,K}\quad\forall
K\in\Omega_{h};\quad\|\mathbf{u}-\mathbf{u}_{I}\|+h\|\mathbf{u}-\mathbf{u}_{I}\|_{1}\leq
C_{I}h^{2}|\mathbf{u}|_{2}.$

The constant $C_{I}$ depends on the mesh parameter $\gamma$.

An important consequence of the definition of the interpolant $\mathbf{u}_{I}$
is that, if $\mathbf{u}$ is divergence free, then (see also Remark 4.1)

$b(\mathbf{u}_{I},q_{h})=b(\mathbf{u},q_{h})=0\qquad\forall q_{h}\in
Q_{h}\quad\Longrightarrow\quad\operatorname{div}\mathbf{u}_{I}=0.$ (16)

We conclude this section with a technical lemma that will be useful in Section
6.

######  Lemma 4.3.

We have a stability bound for the $L^{\infty}$-norm of the elemental $L^{2}$
polynomial projection of $\mathbf{u}_{I}$: for all $\ell\in\mathbb{N}$,

$\quad\|\bm{\Pi}^{0}_{\ell}\mathbf{u}_{I}\|_{L^{\infty}}\leq
C_{\operatorname{inv}}C_{D}\big{(}C_{I}|\mathbf{u}|_{W^{1,3}}+\|\mathbf{u}\|_{L^{\infty}}).$

###### Proof.

On each element $K\in\Omega_{h}$, standard manipulations imply

$\displaystyle\|\bm{\Pi}^{0}_{\ell}\mathbf{u}_{I}\|_{L^{\infty}(K)}$
$\displaystyle\leq\|\bm{\Pi}^{0}_{\ell}(\mathbf{u}-\mathbf{u}_{I})\|_{L^{\infty}(K)}+\|\bm{\Pi}^{0}_{\ell}\mathbf{u}\|_{L^{\infty}(K)}\leq
C_{\operatorname{inv}}h_{K}^{-\frac{3}{2}}\big{(}\|\mathbf{u}-\mathbf{u}_{I}\|_{0,K}+\|\mathbf{u}\|_{0,K}\big{)}$
$\displaystyle\leq
C_{\operatorname{inv}}\big{(}C_{I}h_{K}^{-\frac{1}{2}}|\mathbf{u}|_{1,K}+h_{K}^{-\frac{3}{2}}\|\mathbf{u}\|_{0,K}\Big{)}\leq
C_{\operatorname{inv}}\big{(}C_{I}C_{D}|\mathbf{u}|_{W^{1,3}(K)}+C_{D}\|\mathbf{u}\|_{L^{\infty}(K)}\big{)}.$

Taking the maximum over all elements gives the assertion. ∎

###  4.4 The discrete trilinear forms

We introduce the local continuous and discrete trilinear forms on the element
$K$:

$\displaystyle c^{K}(\mathbf{v};\mathbf{u},\mathbf{w})$
$\displaystyle:=\int_{K}(\nabla\mathbf{u})\mathbf{v}\cdot\mathbf{w}\qquad\forall\mathbf{u},\mathbf{v},\mathbf{w}\in\big{[}H^{1}(K)\big{]}^{3},$
$\displaystyle c_{h}^{K}(\mathbf{v}_{h};\mathbf{u}_{h},\mathbf{w}_{h})$
$\displaystyle:=\int_{K}(\bm{\Pi}^{0,K}_{0}\nabla\mathbf{u}_{h})\bm{\Pi}^{0,K}_{1}\mathbf{v}_{h}\cdot\bm{\Pi}^{0,K}_{1}\mathbf{w}_{h}\qquad\forall\mathbf{u}_{h},\mathbf{v}_{h},\mathbf{w}_{h}\in\mathbf{W}_{h}(K).$

We also define their skew-symmetric counterparts

$\displaystyle\widetilde{c}^{K}(\mathbf{v};\mathbf{u},\mathbf{w})$
$\displaystyle:=\frac{1}{2}\left(c^{K}(\mathbf{v};\mathbf{u},\mathbf{w})-c^{K}(\mathbf{v};\mathbf{w},\mathbf{u})\right)\qquad\forall\mathbf{u},\mathbf{v},\mathbf{w}\in\big{[}H^{1}(K)\big{]}^{3},$
(17)
$\displaystyle\widetilde{c}_{h}^{K}(\mathbf{v}_{h};\mathbf{u}_{h},\mathbf{w}_{h})$
$\displaystyle:=\frac{1}{2}\left(c_{h}^{K}(\mathbf{v}_{h};\mathbf{u}_{h},\mathbf{w}_{h})-c_{h}^{K}(\mathbf{v}_{h};\mathbf{w}_{h},\mathbf{u}_{h})\right)\qquad\forall\mathbf{u}_{h},\mathbf{v}_{h},\mathbf{w}_{h}\in\mathbf{W}_{h}(K),$
(18)

and the associated global skew-symmetric trilinear forms

$\displaystyle\widetilde{c}(\mathbf{v};\mathbf{u},\mathbf{w})$
$\displaystyle=\sum_{K\in\Omega_{h}}\widetilde{c}^{K}(\mathbf{v};\mathbf{u},\mathbf{w})\qquad\forall\mathbf{u},\mathbf{v},\mathbf{w}\in\big{[}H^{1}(\Omega)\big{]}^{3},$
$\displaystyle\widetilde{c}_{h}(\mathbf{v}_{h};\mathbf{u}_{h},\mathbf{w}_{h})$
$\displaystyle=\sum_{K\in\Omega_{h}}\widetilde{c}_{h}^{K}(\mathbf{v}_{h};\mathbf{u}_{h},\mathbf{w}_{h})\qquad\forall\mathbf{u}_{h},\mathbf{v}_{h},\mathbf{w}_{h}\in\mathbf{W}_{h}.$

If $\mathbf{v}\in[H^{1}_{0}(\Omega)]^{3}$ is divergence free, then

$\widetilde{c}(\mathbf{v};\mathbf{u},\mathbf{w})=c(\mathbf{v};\mathbf{u},\mathbf{w}).$

In the remainder of the section, we discuss three properties of these
trilinear forms. For the sake of exposition, we consider the extension of the
virtual element bilinear and trilinear forms to $\big{[}H^{1}(K)\big{]}^{3}$
and $\big{[}H^{1}(\Omega)\big{]}^{3}$.

The first property is the continuity of the trilinear form, which can be
proved as in [7, Proposition $3.3$].

######  Lemma 4.4.

The following continuity bound is valid:

$\displaystyle\left|\widetilde{c}_{h}(\mathbf{v}_{h};\mathbf{u}_{h},\mathbf{w}_{h})\right|\lesssim|\mathbf{u}_{h}|_{1}\,|\mathbf{v}_{h}|_{1}\,|\mathbf{w}_{h}|_{1},$

where the hidden constant is independent of $h$.

The second property is the measure of the variational crime perpetrated in the
discretization of the trilinear form. The proof is a modification of the proof
of [7, Lemma $4.3$].

######  Lemma 4.5.

Given $\mathbf{w}\in\big{[}H^{1}_{0}(\Omega)\big{]}^{3}$ and
$\mathbf{v}\in\big{[}H^{1}_{0}(\Omega)\big{]}^{3}\cap\big{[}H^{2}(\Omega)\big{]}^{3}$,
we have the following bound

$\displaystyle\left|\widetilde{c}(\mathbf{v};\mathbf{v},\mathbf{w})-\widetilde{c}_{h}(\mathbf{v};\mathbf{v},\mathbf{w})\right|\leq
C_{C}h\|\mathbf{v}\|_{2}^{2}\,|\mathbf{w}|_{1},$

where the constant $C_{C}$ depends only on the constant of inverse estimates,
approximation bounds, and Sobolev embedding inequalities, but is independent
of $h_{K}$.

###### Proof.

A straightforward manipulation of relations (17)-(18) defining the skew-
symmetric trilinear forms yields

$\begin{split}\left|\widetilde{c}(\mathbf{v};\mathbf{v},\mathbf{w})-\widetilde{c}_{h}(\mathbf{v};\mathbf{v},\mathbf{w})\right|&\leq\frac{1}{2}\left(\left|c(\mathbf{v};\mathbf{v},\mathbf{w})-c_{h}(\mathbf{v};\mathbf{v},\mathbf{w})\right|+\left|c(\mathbf{v};\mathbf{w},\mathbf{v})-c_{h}(\mathbf{v};\mathbf{w},\mathbf{v})\right|\right)\\\
&=\frac{1}{2}\left(\textrm{T}_{1}+\textrm{T}_{2}\right).\end{split}$ (19)

We estimate the two terms on the right-hand side separately. We begin with the
first one. Using the definition of the orthogonal projection
$\bm{\Pi}^{0}_{1}\mathbf{v}$, and adding and subtracting
$(\nabla\mathbf{v})\mathbf{v}\cdot\bm{\Pi}^{0}_{1}\mathbf{w}$ to the integral
argument, we write

$\displaystyle\begin{split}\textrm{T}_{1}&=\left|\int_{\Omega}(\nabla\mathbf{v})\mathbf{v}\cdot\mathbf{w}-\int_{\Omega}(\bm{\Pi}^{0}_{0}\nabla\mathbf{v})\bm{\Pi}^{0}_{1}\mathbf{v}\cdot\bm{\Pi}^{0}_{1}\mathbf{w}\right|\\\\[5.0pt]
&=\left|\int_{\Omega}(\nabla\mathbf{v})\mathbf{v}\cdot\mathbf{w}-\int_{\Omega}(\bm{\Pi}^{0}_{0}\nabla\mathbf{v})\mathbf{v}\cdot\bm{\Pi}^{0}_{1}\mathbf{w}\right|\\\\[5.0pt]
&\leq\left|\int_{\Omega}(\nabla\mathbf{v})\mathbf{v}\cdot(\mathbf{w}-\bm{\Pi}^{0}_{1}\mathbf{w})\right|+\left|\int_{\Omega}(\nabla\mathbf{v}-\bm{\Pi}^{0}_{0}\nabla\mathbf{v})\mathbf{v}\cdot\bm{\Pi}^{0}_{1}\mathbf{w}\right|=:\textrm{T}_{1,1}+\textrm{T}_{1,2}.\end{split}$

Again, we focus on the two terms on the right-hand side separately. Using the
Hölder inequality, polynomial approximation properties, and the Sobolev
embedding theorem gives

$\displaystyle\begin{split}\textrm{T}_{1,1}&\leq\|\nabla\mathbf{v}\|_{L^{4}}\,\|\mathbf{v}\|_{L^{4}}\,\|\mathbf{w}-\bm{\Pi}^{0}_{1}\mathbf{w}\|\leq
C_{\operatorname{appr}}h\|\nabla\mathbf{v}\|_{L^{4}}\,\|\mathbf{v}\|_{L^{4}}\|\mathbf{w}\|_{1}\leq
C_{S}^{2}C_{\operatorname{appr}}h\|\mathbf{v}\|_{1}\,\|\mathbf{v}\|_{2}\|\mathbf{w}\|_{1}.\end{split}$

Next, for each element $K\in\Omega_{h}$, we apply a polynomial inverse
inequality, the Hölder inequality, and standard manipulations:

$\|\bm{\Pi}^{0}_{1}\mathbf{w}\|_{L^{4}(K)}\leq
C_{\operatorname{inv}}h_{K}^{-\frac{3}{4}}\|\mathbf{w}\|_{0,K}\leq
C_{\operatorname{inv}}h_{K}^{-\frac{3}{4}}\|1\|_{L^{4}(K)}\,\|\mathbf{w}\|_{L^{4}(K)}\leq
C_{\operatorname{inv}}C_{D}\|\mathbf{w}\|_{L^{4}(K)}.$

This entails

$\displaystyle\textrm{T}_{1,2}$
$\displaystyle\leq\sum_{K\in\Omega_{h}}\|\nabla\mathbf{v}-\bm{\Pi}^{0}_{0}\nabla\mathbf{v}\|_{0,K}\,\|\mathbf{v}\|_{L^{4}(K)}\,\|\bm{\Pi}^{0}_{1}\mathbf{w}\|_{L^{4}(K)}$
$\displaystyle\leq
C_{\operatorname{inv}}C_{D}C_{\operatorname{appr}}h\sum_{K\in\Omega_{h}}|\mathbf{v}|_{2,K}\,\|\mathbf{v}\|_{L^{4}(K)}\,\|\mathbf{w}\|_{L^{4}(K)}.$

We apply a sequential $\ell^{2}-\ell^{4}-\ell^{4}$ Hölder inequality and the
Sobolev embedding theorem twice, and deduce

$\textrm{T}_{1,2}\leq
C_{\operatorname{inv}}C_{D}C_{\operatorname{appr}}h|\mathbf{v}|_{2}\,\|\mathbf{v}\|_{L^{4}}\,\|\mathbf{w}\|_{L^{4}}\leq
C_{\operatorname{appr}}C_{\operatorname{inv}}C_{D}C_{S}^{2}h\|\mathbf{v}\|_{2}\,\|\mathbf{v}\|_{1}\,\|\mathbf{w}\|_{1}.$

This step concludes the derivation of a bound on the term $\textrm{T}_{1}$ in
(19).

Next, we focus on the term $\textrm{T}_{2}$. Adding and subtracting
$\big{(}\bm{\Pi}^{0}_{0}\nabla\mathbf{w}\big{)}\mathbf{v}\cdot\mathbf{v}$
yields

$\displaystyle\textrm{T}_{2}$
$\displaystyle=\left|\int_{\Omega}(\nabla\mathbf{w})\mathbf{v}\cdot\mathbf{v}-\int_{\Omega}(\bm{\Pi}^{0}_{0}\nabla\mathbf{w})\mathbf{v}\cdot\bm{\Pi}^{0}_{1}\mathbf{v}\right|$
$\displaystyle\leq\left|\int_{\Omega}(\nabla\mathbf{w}-\bm{\Pi}^{0}_{0}\nabla\mathbf{w})\mathbf{v}\cdot\mathbf{v}\right|+\left|\int_{\Omega}(\bm{\Pi}^{0}_{0}\nabla\mathbf{w})\mathbf{v}\cdot(\mathbf{v}-\bm{\Pi}^{0}_{1}\mathbf{v})\right|=:\textrm{T}_{2,1}+\textrm{T}_{2,2}.$

We derive the bounds on the terms $\textrm{T}_{2,1}$ and $\textrm{T}_{2,2}$ by
using the approximation properties of the orthogonal projections, and
theoretical tools similar to those used for the bound on $\textrm{T}_{1}$:

$\displaystyle\textrm{T}_{2,1}$
$\displaystyle=\left|\int_{\Omega}(\nabla\mathbf{w})(\mathbf{v}\cdot\mathbf{v}-\bm{\Pi}^{0}_{0}\mathbf{v}\cdot\mathbf{v})\right|\leq|\mathbf{w}|_{1}\
\|\mathbf{v}\cdot\mathbf{v}-\bm{\Pi}^{0}_{0}(\mathbf{v}\cdot\mathbf{v})\|\leq
C_{\operatorname{appr}}h|\mathbf{w}|_{1}\ |\mathbf{v}\cdot\mathbf{v}|_{1}$
$\displaystyle\leq C_{\operatorname{appr}}h|\mathbf{w}|_{1}\
\|\mathbf{v}\|_{W^{1,4}}^{2}\leq
C_{\operatorname{appr}}C_{S}^{2}h\|\mathbf{v}\|_{2}^{2}|\mathbf{w}|_{1}.$

On the other hand, given $\mathbf{v}_{\pi}$ the best piecewise linear
approximant in $L^{4}$ of $\mathbf{v}$ over $\Omega_{h}$, we also have

$\displaystyle\textrm{T}_{2,2}$
$\displaystyle\leq\sum_{K\in\Omega_{h}}|\mathbf{w}|_{1,K}\,\|\mathbf{v}\|_{L^{4}(K)}\,\|\mathbf{v}-\bm{\Pi}^{0}_{1}\mathbf{v}\|_{L^{4}(K)}$
$\displaystyle\leq\sum_{K\in\Omega_{h}}|\mathbf{w}|_{1,K}\,\|\mathbf{v}\|_{L^{4}(K)}\left(\|\mathbf{v}-\mathbf{v}_{\pi}\|_{L^{4}(K)}+\|\mathbf{v}_{\pi}-\bm{\Pi}^{0}_{1}\mathbf{v}\|_{L^{4}(K)}\right).$

Observe that

$\displaystyle\|\mathbf{v}_{\pi}-\bm{\Pi}^{0,K}_{1}\mathbf{v}\|_{L^{4}(K)}$
$\displaystyle\leq
C_{\operatorname{inv}}h_{K}^{-\frac{3}{4}}(\|\mathbf{v}-\mathbf{v}_{\pi}\|_{0,K}+\|\mathbf{v}-\bm{\Pi}^{0,K}_{1}\mathbf{v}\|_{0,K})$
$\displaystyle\leq
2C_{\operatorname{inv}}h_{K}^{-\frac{3}{4}}\|\mathbf{v}-\mathbf{v}_{\pi}\|_{0,K}\leq
2C_{\operatorname{inv}}C_{D}\|\mathbf{v}-\mathbf{v}_{\pi}\|_{L^{4}(K)}$
$\displaystyle\leq
2C_{\operatorname{inv}}C_{D}C_{\operatorname{appr}}h_{K}|\mathbf{v}|_{W^{1,4}(K)}.$

Inserting this bound above and using polynomial approximation estimates yield

$\textrm{T}_{2,2}\leq(1+2C_{\operatorname{inv}}C_{D})C_{\operatorname{appr}}h\sum_{K\in\Omega_{h}}|\mathbf{w}|_{1,K}\|\mathbf{v}\|_{L^{4}(K)}\|\mathbf{v}\|_{W^{1,4}(K)}.$

We apply a sequential $\ell^{2}-\ell^{4}-\ell^{4}$ Hölder inequality and the
Sobolev embedding theorem, and deduce

$\textrm{T}_{2,2}\leq(1+2C_{\operatorname{inv}}C_{D})C_{\operatorname{appr}}h|\mathbf{w}|_{1}\|\mathbf{v}\|_{L^{4}}\|\mathbf{v}\|_{W^{1,4}}\leq(1+2C_{\operatorname{inv}}C_{D})C_{\operatorname{appr}}C_{S}^{2}h\|\mathbf{v}\|_{1}\|\mathbf{v}\|_{2}\|\mathbf{w}\|_{1}.$

Collecting all the bounds together yields the assertion. ∎

The last property measures the distance between the continuous and the
discrete solutions through the discrete trilinear form. The proof differs from
that of [7, Lemma $4.4$], since here we are interested in the analysis of the
time-dependent case. For this reason, we also use some techniques from [69,
Lemma $4.1$]. Further, we require some extra regularity on the exact velocity
solution $\mathbf{u}$ to (6).

######  Lemma 4.6.

Let $\mathbf{u}\in W^{2,3}(\Omega)$ and
$\mathbf{e}_{h}^{\mathbf{u}}:=\mathbf{u}_{h}-\mathbf{u}_{I}$, where
$\mathbf{u}_{I}$ is the degrees of freedom interpolant of $\mathbf{u}$; see
Lemma 4.2. Then, for every positive $\varepsilon$, there exist two positive
constants $C_{1}$ and $C_{2}$ independent of $h$ such that

$\displaystyle\left|\widetilde{c}_{h}(\mathbf{u};\mathbf{u},\mathbf{e}_{h}^{\mathbf{u}})-\widetilde{c}_{h}(\mathbf{u}_{h};\mathbf{u}_{h},\mathbf{e}_{h}^{\mathbf{u}})\right|\leq\varepsilon|\mathbf{e}_{h}^{\mathbf{u}}|_{1}^{2}+\left(1+\frac{1}{\varepsilon}\right)C_{1}R_{1}(\mathbf{u})\|\mathbf{e}_{h}^{\mathbf{u}}\|^{2}+\left(1+\frac{1}{\varepsilon}\right)C_{2}R_{2}(\mathbf{u})h^{2},$

where

$R_{1}(\mathbf{u}):=\|\mathbf{u}\|_{W^{2,3}}^{2}+\|\mathbf{u}\|_{W^{2,3}}+\|\mathbf{u}\|_{W^{1,\infty}}+1,\qquad
R_{2}(\mathbf{u}):=\|\mathbf{u}\|_{W^{2,3}}^{4}+\|\mathbf{u}\|_{W^{2,3}}^{2}|\mathbf{u}|_{W^{1,\infty}}^{2}.$

###### Proof.

The triangle inequality implies

$\displaystyle\left|\widetilde{c}_{h}(\mathbf{u};\mathbf{u},\mathbf{e}_{h}^{\mathbf{u}})-\widetilde{c}_{h}(\mathbf{u}_{h};\mathbf{u}_{h},\mathbf{e}_{h}^{\mathbf{u}})\right|$
$\displaystyle\leq\left|\widetilde{c}_{h}(\mathbf{u};\mathbf{u},\mathbf{e}_{h}^{\mathbf{u}})-\widetilde{c}_{h}(\mathbf{u}_{h};\mathbf{u},\mathbf{e}_{h}^{\mathbf{u}})\right|+\left|\widetilde{c}_{h}(\mathbf{u}_{h};\mathbf{u},\mathbf{e}_{h}^{\mathbf{u}})-\widetilde{c}_{h}(\mathbf{u}_{h};\mathbf{u}_{h},\mathbf{e}_{h}^{\mathbf{u}})\right|=:\textrm{T}_{1}+\textrm{T}_{2}.$

We estimate the two terms on the right-hand side separately and begin by
splitting $\textrm{T}_{1}$ into the sum of two further contributions, due to
the definition of the skew symmetric form
$\widetilde{c}_{h}(\cdot;\cdot,\cdot)$:

$\textrm{T}_{1}\leq\frac{1}{2}(\left|c_{h}(\mathbf{u};\mathbf{u},\mathbf{e}_{h}^{\mathbf{u}})-c_{h}(\mathbf{u}_{h};\mathbf{u},\mathbf{e}_{h}^{\mathbf{u}})\right|+\left|c_{h}(\mathbf{u};\mathbf{e}_{h}^{\mathbf{u}},\mathbf{u})-c_{h}(\mathbf{u}_{h};\mathbf{e}_{h}^{\mathbf{u}},\mathbf{u})\right|)=:\frac{1}{2}(\textrm{T}_{1,1}+\textrm{T}_{1,2}).$

Using standard manipulations, we get

$\displaystyle\textrm{T}_{1,1}$
$\displaystyle=\left|\int_{\Omega}(\bm{\Pi}^{0}_{0}\nabla\mathbf{u})(\mathbf{u}-\mathbf{u}_{h})\cdot\bm{\Pi}^{0}_{1}\mathbf{e}_{h}^{\mathbf{u}}\right|\leq\|\bm{\Pi}^{0}_{0}\nabla\mathbf{u}\|_{L^{\infty}}\,\|\mathbf{u}-\mathbf{u}_{h}\|\,\|\bm{\Pi}^{0}_{1}\mathbf{e}_{h}^{\mathbf{u}}\|$
$\displaystyle\leq|\mathbf{u}|_{W^{1,\infty}}\left(\|\mathbf{u}-\mathbf{u}_{I}\|+\|\mathbf{e}_{h}^{\mathbf{u}}\|\right)\|\mathbf{e}_{h}^{\mathbf{u}}\|\leq
C_{I}h|\mathbf{u}|_{W^{1,\infty}}\,|\mathbf{u}|_{1}\,\|\mathbf{e}_{h}^{\mathbf{u}}\|+|\mathbf{u}|_{W^{1,\infty}}\,\|\mathbf{e}_{h}^{\mathbf{u}}\|^{2}$
$\displaystyle\leq
C_{I}h^{2}|\mathbf{u}|_{W^{1,\infty}}^{2}\,\|\mathbf{u}\|_{1}^{2}+(|\mathbf{u}|_{W^{1,\infty}}+C_{I})\|\mathbf{e}_{h}^{\mathbf{u}}\|^{2}$
$\displaystyle\leq
C_{I}C_{S}^{2}h^{2}|\mathbf{u}|_{W^{1,\infty}}^{2}\|\mathbf{u}\|_{W^{2,3}}^{2}+(|\mathbf{u}|_{W^{1,\infty}}+C_{I})\|\mathbf{e}_{h}^{\mathbf{u}}\|^{2}.$

As for the term $\textrm{T}_{1,2}$, we write

$\displaystyle\textrm{T}_{1,2}$
$\displaystyle=\left|\int_{\Omega}(\bm{\Pi}^{0}_{0}\nabla\mathbf{e}_{h}^{\mathbf{u}})(\bm{\Pi}^{0}_{1}(\mathbf{u}-\mathbf{u}_{h}))\cdot\mathbf{u}\right|\leq|\mathbf{e}_{h}^{\mathbf{u}}|_{1}\,\|\mathbf{u}-\mathbf{u}_{h}\|\,\|\mathbf{u}\|_{L^{\infty}}$
$\displaystyle\leq|\mathbf{e}_{h}^{\mathbf{u}}|_{1}\,\left(\|\mathbf{u}-\mathbf{u}_{I}\|+\|\mathbf{e}_{h}^{\mathbf{u}}\|\right)\|\mathbf{u}\|_{L^{\infty}}\leq
C_{I}h|\mathbf{e}_{h}^{\mathbf{u}}|_{1}\,|\mathbf{u}|_{1}\,\|\mathbf{u}\|_{L^{\infty}}+|\mathbf{e}_{h}^{\mathbf{u}}|_{1}\,\|\mathbf{e}_{h}^{\mathbf{u}}\|\,\|\mathbf{u}\|_{L^{\infty}}$
$\displaystyle\leq\frac{\varepsilon}{2}|\mathbf{e}_{h}^{\mathbf{u}}|_{1}^{2}+\frac{1}{\varepsilon}C_{I}^{2}h^{2}|\mathbf{u}|_{1}^{2}\,\|\mathbf{u}\|_{L^{\infty}}^{2}+\frac{1}{\varepsilon}\|\mathbf{u}\|_{L^{\infty}}^{2}\,\|\mathbf{e}_{h}^{\mathbf{u}}\|^{2}$
$\displaystyle\leq\frac{\varepsilon}{2}|\mathbf{e}_{h}^{\mathbf{u}}|_{1}^{2}+\frac{1}{\varepsilon}C_{I}^{2}C_{S}^{2}h^{2}\|\mathbf{u}\|_{W^{2,3}}^{4}+\frac{1}{\varepsilon}\|\mathbf{u}\|_{L^{\infty}}^{2}\,\|\mathbf{e}_{h}^{\mathbf{u}}\|^{2}.$

Next, we focus on the term $\textrm{T}_{2}$. Note that

$\widetilde{c}_{h}^{K}(\mathbf{u}_{h};\mathbf{u}_{h},\mathbf{e}_{h}^{\mathbf{u}})=\widetilde{c}_{h}^{K}(\mathbf{u}_{h};\mathbf{u}_{I},\mathbf{e}_{h}^{\mathbf{u}}),$
(20)

which is a consequence of the skew symmetry of
$\widetilde{c}_{h}^{K}(\cdot;\cdot,\cdot)$ and (16).

We split the term $\textrm{T}_{2}$ into two contributions, due to the
definition of the skew symmetric form $\widetilde{c}_{h}(\cdot;\cdot,\cdot)$.
Recalling (20), we write

$\textrm{T}_{2}\leq\frac{1}{2}\left(\left|c_{h}(\mathbf{u}_{h};\mathbf{u},\mathbf{e}_{h}^{\mathbf{u}})-c_{h}(\mathbf{u}_{h};\mathbf{u}_{I},\mathbf{e}_{h}^{\mathbf{u}})\right|+\left|c_{h}(\mathbf{u}_{h};\mathbf{e}_{h}^{\mathbf{u}},\mathbf{u})-c_{h}(\mathbf{u}_{h};\mathbf{e}_{h}^{\mathbf{u}},\mathbf{u}_{I})\right|\right)=:\frac{1}{2}(\textrm{T}_{2,1}+\textrm{T}_{2,2}).$
(21)

Using standard manipulations as above, we get

$\displaystyle\textrm{T}_{2,1}\leq\left|c_{h}(\mathbf{e}_{h}^{\mathbf{u}};\mathbf{u}-\mathbf{u}_{I},\mathbf{e}_{h}^{\mathbf{u}})\right|+\left|c_{h}(\mathbf{u}_{I};\mathbf{u}-\mathbf{u}_{I},\mathbf{e}_{h}^{\mathbf{u}})\right|=:\textrm{T}_{2,1,1}+\textrm{T}_{2,1,2}.$

Preliminary, for all $K\in\Omega_{h}$, we observe that

$\begin{split}\|\bm{\Pi}^{0}_{0}\nabla(\mathbf{u}-\mathbf{u}_{I})\|_{L^{\infty}(K)}&\leq
C_{D}h_{K}^{-\frac{3}{2}}|\mathbf{u}-\mathbf{u}_{I}|_{1,K}\leq
C_{D}C_{I}h_{K}^{-\frac{1}{2}}|\mathbf{u}|_{2,K}\leq
C_{D}^{2}C_{I}|\mathbf{u}|_{W^{2,3}(K)}.\end{split}$

By taking the maximum over all the elements on both sides entails

$\|\bm{\Pi}^{0}_{0}\nabla(\mathbf{u}-\mathbf{u}_{I})\|_{L^{\infty}}\leq
C_{D}^{2}C_{I}|\mathbf{u}|_{W^{2,3}}.$

Then, we have

$\displaystyle\textrm{T}_{2,1,1}$
$\displaystyle\leq\int_{\Omega}\left|\bm{\Pi}^{0}_{0}\nabla(\mathbf{u}-\mathbf{u}_{I})\right|\,\left|\bm{\Pi}^{0}_{1}\mathbf{e}_{h}^{\mathbf{u}}\right|\,\left|\bm{\Pi}^{0}_{1}\mathbf{e}_{h}^{\mathbf{u}}\right|\leq\|\bm{\Pi}^{0}_{0}\nabla(\mathbf{u}-\mathbf{u}_{I})\|_{L^{\infty}}\,\|\mathbf{e}_{h}^{\mathbf{u}}\|^{2}$
$\displaystyle\leq
C_{D}^{2}C_{I}|\mathbf{u}|_{W^{2,3}}\,\|\mathbf{e}_{h}^{\mathbf{u}}\|^{2}.$

Next, we apply Lemma 4.3 and write

$\|\bm{\Pi}^{0,K}_{1}\mathbf{u}_{I}\|_{L^{\infty}(K)}\leq
C_{\operatorname{inv}}C_{D}\Big{(}C_{I}|\mathbf{u}|_{W^{1,3}}+\|\mathbf{u}\|_{L^{\infty}}\Big{)}.$

By means of the above inequality, we deduce the bound

$\displaystyle\textrm{T}_{2,1,2}$
$\displaystyle\leq\int_{\Omega}\left|\bm{\Pi}^{0}_{0}\nabla(\mathbf{u}-\mathbf{u}_{I})\right|\,\left|\bm{\Pi}^{0}_{1}\mathbf{u}_{I}\right|\,\left|\mathbf{e}_{h}^{\mathbf{u}}\right|\leq\|\bm{\Pi}^{0}_{0}\nabla(\mathbf{u}-\mathbf{u}_{I})\|\,\|\bm{\Pi}^{0}_{1}\mathbf{u}_{I}\|_{L^{\infty}}\,\|\mathbf{e}_{h}^{\mathbf{u}}\|$
$\displaystyle\leq
C_{I}h|\mathbf{u}|_{2}C_{\operatorname{inv}}C_{D}(C_{I}+1)(\|\mathbf{u}\|_{W^{1,3}}+\|\mathbf{u}\|_{L^{\infty}})\|\mathbf{e}_{h}^{\mathbf{u}}\|$
$\displaystyle\leq
C_{I}^{2}C_{\operatorname{inv}}^{2}C_{D}^{2}(C_{I}+1)^{2}h^{2}|\mathbf{u}|_{2}^{2}(\|\mathbf{u}\|_{W^{1,3}}^{2}+\|\mathbf{u}\|_{L^{\infty}}^{2})+\frac{1}{4}\|\mathbf{e}_{h}^{\mathbf{u}}\|^{2}$
$\displaystyle\leq
C_{I}^{2}C_{\operatorname{inv}}^{2}C_{D}^{2}(C_{I}+1)^{2}(C_{S}^{2}+1)h^{2}\|\mathbf{u}\|_{W^{2,3}}^{4}+\frac{1}{4}\|\mathbf{e}_{h}^{\mathbf{u}}\|^{2}.$

Next, we focus on the term $\textrm{T}_{2,2}$ appearing on the right-hand side
of (21):

$\displaystyle\textrm{T}_{2,2}$
$\displaystyle\leq\left|\int_{\Omega}(\bm{\Pi}^{0}_{0}\nabla\mathbf{e}_{h}^{\mathbf{u}})\mathbf{u}_{h}\cdot\bm{\Pi}^{0}_{1}(\mathbf{u}-\mathbf{u}_{I})\right|$
$\displaystyle\leq\left|\int_{\Omega}(\bm{\Pi}^{0}_{0}\nabla\mathbf{e}_{h}^{\mathbf{u}})\mathbf{e}_{h}^{\mathbf{u}}\cdot\bm{\Pi}^{0}_{1}(\mathbf{u}-\mathbf{u}_{I})\right|+\left|\int_{\Omega}(\bm{\Pi}^{0}_{0}\nabla\mathbf{e}_{h}^{\mathbf{u}})\mathbf{u}_{I}\cdot\bm{\Pi}^{0}_{1}(\mathbf{u}-\mathbf{u}_{I})\right|=:\textrm{T}_{2,2,1}+\textrm{T}_{2,2,2}.$

We have the bounds

$\displaystyle\textrm{T}_{2,2,1}$
$\displaystyle\leq\sum_{K\in\Omega_{h}}\|\bm{\Pi}^{0,K}_{0}\nabla\mathbf{e}_{h}^{\mathbf{u}}\|_{0,K}\,\|\mathbf{e}_{h}^{\mathbf{u}}\|_{0,K}\,\|\bm{\Pi}^{0,K}_{1}(\mathbf{u}-\mathbf{u}_{I})\|_{L^{\infty}(K)}$
$\displaystyle\leq\sum_{K\in\Omega_{h}}|\mathbf{e}_{h}^{\mathbf{u}}|_{1,K}\,\|\mathbf{e}_{h}^{\mathbf{u}}\|_{0,K}\,C_{\operatorname{inv}}h_{K}^{-\frac{3}{2}}\,\|\mathbf{u}-\mathbf{u}_{I}\|_{0,K}$
$\displaystyle\leq\sum_{K\in\Omega_{h}}|\mathbf{e}_{h}^{\mathbf{u}}|_{1,K}\|\mathbf{e}_{h}^{\mathbf{u}}\|_{0,K}C_{\operatorname{inv}}C_{I}h_{K}^{-\frac{1}{2}}|\mathbf{u}|_{1,K}\leq
C_{\operatorname{inv}}C_{I}C_{D}|\mathbf{u}|_{W^{1,3}}\sum_{K\in\Omega_{h}}|\mathbf{e}_{h}^{\mathbf{u}}|_{1,K}\|\mathbf{e}_{h}^{\mathbf{u}}\|_{0,K}$
$\displaystyle\leq
C_{\operatorname{inv}}C_{I}C_{D}|\mathbf{u}|_{W^{1,3}}|\mathbf{e}_{h}^{\mathbf{u}}|_{1}\|\mathbf{e}_{h}^{\mathbf{u}}\|\leq\frac{\varepsilon}{4}|\mathbf{e}_{h}^{\mathbf{u}}|_{1}^{2}+\frac{1}{\varepsilon}C_{\operatorname{inv}}^{2}C_{I}^{2}C_{D}^{2}\,|\mathbf{u}|_{W^{2,3}}^{2}\,\|\mathbf{e}_{h}^{\mathbf{u}}\|^{2}$

and, with similar computations,

$\displaystyle\begin{split}\textrm{T}_{2,2,2}&\leq\sum_{K\in\Omega_{h}}|\mathbf{e}_{h}^{\mathbf{u}}|_{1,K}\,\|\mathbf{u}_{I}\|_{0,K}\,\|\bm{\Pi}^{0,K}_{1}(\mathbf{u}-\mathbf{u}_{I})\|_{L^{\infty}(K)}\\\
&\leq\sum_{K\in\Omega_{h}}|\mathbf{e}_{h}^{\mathbf{u}}|_{1,K}(1+C_{I}h_{K})\|\mathbf{u}\|_{1,K}\,C_{\operatorname{inv}}C_{I}h_{K}^{\frac{1}{2}}|\mathbf{u}|_{2,K}\\\
&\leq\sum_{K\in\Omega_{h}}|\mathbf{e}_{h}^{\mathbf{u}}|_{1,K}(1+C_{I}h_{K})\|\mathbf{u}\|_{1,K}\,C_{\operatorname{inv}}C_{I}C_{D}h_{K}|\mathbf{u}|_{W^{2,3}(K)}\\\
&\leq\frac{\varepsilon}{4}|\mathbf{e}_{h}^{\mathbf{u}}|_{1}^{2}+\frac{1}{\varepsilon}(1+C_{I}h_{K})^{2}C_{\operatorname{inv}}^{2}C_{I}^{2}C_{D}^{2}h^{2}\|\mathbf{u}\|_{1}^{2}\,|\mathbf{u}|_{W^{2,3}}^{2}\\\
&\leq\frac{\varepsilon}{4}|\mathbf{e}_{h}^{\mathbf{u}}|_{1}^{2}+\frac{1}{\varepsilon}(1+C_{I}h_{K})^{2}C_{\operatorname{inv}}^{2}C_{I}^{2}C_{D}^{2}h^{2}\|\mathbf{u}\|_{W^{2,3}}^{4}.\end{split}$

We collect the bounds on all the terms above and the assertion follows. ∎

##  5 Virtual element spaces for the electromagnetic equations

###  5.1 A pivot nodal space

Here, we introduce a nodal virtual element space that will be instrumental in
defining polynomial projectors for edge spaces. Consider the mesh face
$F\in\mathcal{F}_{h}$ and the translated position vector
$\mathbf{x}_{F}=\mathbf{x}-\mathbf{b}_{F}$ for all $\mathbf{x}\in F$ (we
recall that $\mathbf{b}_{F}$ is the center of $F$). We define the nodal
virtual element space on $F$ as

$\displaystyle V_{h}^{\operatorname{node}}(F):=\bigg{\\{}v_{h}\in
C^{0}(\overline{F})\,\mid\,\Delta
v_{h}\in\mathbbm{P}_{0}(F),\;{v_{h}}_{|e}\in\mathbbm{P}_{1}(e)\,\,\forall
e\in\mathcal{E}^{F}_{h},\;\int_{F}\nabla
v_{h}\cdot\mathbf{x}_{F}=0\bigg{\\}}.$

Then, for all elements $K\in\Omega_{h}$, we define the local nodal space as:

$\displaystyle V_{h}^{\operatorname{node}}(K):=\Big{\\{}v_{h}\in
C^{0}(\overline{K})\,\mid\,\Delta v_{h}=0,\;{v_{h}}_{|F}\in
V_{h}^{\operatorname{node}}(F)\,\,\forall F\in\mathcal{F}^{K}_{h}\Big{\\}}.$

The virtual element functions in the local space
$V_{h}^{\operatorname{node}}(K)$ for all $K\in\Omega_{h}$ are uniquely
determined by the set of their vertex values over $\partial K$; see, e.g.,
[11] for a proof of their unisolvence. We introduce the global nodal space

$V_{h}^{\operatorname{node}}:=\\{v_{h}\in\mathcal{C}^{0}(\overline{\Omega})\mid
v_{h}{}_{|K}\in V_{h}^{\operatorname{node}}(K)\;\forall K\in\Omega_{h}\\},$
(22)

which will be used in the definitions of the edge interpolant
$\mathbf{E}_{I\\!I}$ in the forthcoming Section 5.4.

###  5.2 The edge space

We define the virtual element edge space on the mesh face
$F\in\mathcal{F}_{h}$ as

$\displaystyle\mathbf{V}_{h}^{\operatorname{edge}}(F):=\bigg{\\{}\mathbf{E}_{h}\in\big{[}L^{2}(F)\big{]}^{2}\,\mid\,$
$\displaystyle\operatorname{div_{F}}\mathbf{E}_{h}\in\mathbbm{P}_{0}(F),\;\operatorname{rot_{F}}\mathbf{E}_{h}\in\mathbbm{P}_{0}(F),$
$\displaystyle\hskip
2.84526pt\mathbf{E}_{h}\cdot\mathbf{t}_{e}\in\mathbbm{P}_{0}(e)\,\,\forall
e\in\mathcal{E}^{F}_{h},\;\int_{F}\mathbf{E}_{h}\cdot\mathbf{x}_{F}=0\bigg{\\}},$

where $\operatorname{div_{F}}$ and $\operatorname{rot_{F}}$ are defined in
(1), and again $\mathbf{x}_{F}=\mathbf{x}-\mathbf{b}_{F}$.

The local edge space on an element $K\in\Omega_{h}$ is defined as

$\displaystyle\mathbf{V}_{h}^{\operatorname{edge}}(K):=\bigg{\\{}\mathbf{E}_{h}\in\big{[}L^{2}(K)\big{]}^{3}\,\mid\,\,$
$\displaystyle\mathbf{E}_{h,{\bm{\tau}}_{F}}\in\mathbf{V}_{h}^{\operatorname{edge}}(F)\,\,\forall
F\in\mathcal{F}^{K}_{h},\,\,\mathbf{E}_{h}\cdot\mathbf{t}_{e}\text{\leavevmode\nobreak\
continuous at each edge\leavevmode\nobreak\ $e$,\leavevmode\nobreak\ }$
$\displaystyle\begin{cases}\operatorname{\mathbf{curl}}\operatorname{\mathbf{curl}}\mathbf{E}_{h}\in\big{[}\mathbbm{P}_{0}(K)\big{]}^{3}\\\
\operatorname{div}\mathbf{E}_{h}=0\end{cases},$ $\displaystyle\hskip
2.84526pt\int_{K}\operatorname{\mathbf{curl}}\mathbf{E}_{h}\cdot(\mathbf{x}_{K}\times\mathbf{q})=0\,\,\forall\mathbf{q}\in\big{[}\mathbbm{P}_{0}(K)\big{]}^{3}\bigg{\\}},$

where $\mathbf{x}_{K}=\mathbf{x}-\mathbf{b}_{K}$ for all $\mathbf{x}\in K$,
$\mathbf{b}_{K}$ is the center of $K$, and
$\mathbf{E}_{h,{\bm{\tau}}_{F}}=(\mathbf{n}_{F}\times\mathbf{v}_{h})\times\mathbf{n}_{F}$.

Each element $\mathbf{E}_{h}\in\mathbf{V}_{h}^{\operatorname{edge}}(K)$ is
uniquely determined by the set of constant values
$(E_{e})_{e\in\mathcal{E}^{K}_{h}}$, where
$E_{e}=\mathbf{E}_{h}\cdot\mathbf{t}_{e}$ is the tangential component of
$\mathbf{E}_{h}$ along the elemental edge $e$. We take these values as the
degrees of freedom of $\mathbf{E}_{h}$ and refer the reader to [11] for the
proof of their unisolvence in $\mathbf{V}_{h}^{\operatorname{edge}}(K)$. A
major consequence of such a space definition is that
$\bm{\Pi}^{0,K}_{0}\mathbf{E}_{h}$ is computable from the degrees of freedom
of $\mathbf{E}_{h}$ for all
$\mathbf{E}_{h}\in\mathbf{V}_{h}^{\operatorname{edge}}(K)$.

We define the global edge space $\mathbf{V}_{h}^{\operatorname{edge}}$ by an
$H(\operatorname{\mathbf{curl}},\Omega)$-conforming coupling of the degrees of
freedom, so that

$\mathbf{V}_{h}^{\operatorname{edge}}:=\Big{\\{}\mathbf{E}_{h}\in
H_{0}(\operatorname{\mathbf{curl}},\Omega)\,\mid\,{\mathbf{E}_{h}}_{|K}\in\mathbf{V}_{h}^{\operatorname{edge}}(K)\,\,\forall
K\in\Omega_{h}\Big{\\}}.$ (23)

This definition includes the homogeneous tangential boundary conditions on
$\partial\Omega$.

####  5.2.1 Virtual element inner product

According to [11, 13], we equip the virtual element edge space
$\mathbf{V}_{h}^{\operatorname{edge}}$ with an inner product mimicking the
$L^{2}$ inner product. Notably, we first introduce the local bilinear form

$\big{[}\mathbf{E}_{h},\mathbf{F}_{h}\big{]}_{\operatorname{edge},K}:=\Big{(}\bm{\Pi}^{0,K}_{0}\mathbf{E}_{h},\bm{\Pi}^{0,K}_{0}\mathbf{F}_{h}\Big{)}_{K}+S^{K}_{\operatorname{edge}}\Big{(}(I-\bm{\Pi}^{0,K}_{0})\mathbf{E}_{h},(I-\bm{\Pi}^{0,K}_{0})\mathbf{F}_{h}\Big{)}\\\
\forall\mathbf{E}_{h},\,\mathbf{F}_{h}\in\mathbf{V}_{h}^{\operatorname{edge}}(K),$
(24)

where $S^{K}_{\operatorname{edge}}(\cdot,\cdot)$ can be any computable,
symmetric bilinear form such that there exist two positive constants
$\gamma_{*}$ and $\gamma^{*}$ independent of $h_{K}$ satisfying

$\displaystyle\gamma_{*}|\mathbf{E}_{h}|_{0,K}^{2}\leq
S^{K}_{\operatorname{edge}}(\mathbf{E}_{h},\mathbf{E}_{h})\leq\gamma^{*}|\mathbf{E}_{h}|_{0,K}^{2}\quad\quad\forall\mathbf{E}_{h}\in\ker(\bm{\Pi}^{0,K}_{0})\cap\mathbf{V}_{h}^{\operatorname{edge}}(K).$
(25)

An explicit stabilization satisfying (25) was introduced in [11, formula
(4.8)] and analyzed in [8, Proposition 5.5], and reads

$S^{K}_{\operatorname{edge}}(\mathbf{E}_{h},\mathbf{F}_{h}):=h_{K}^{2}\sum_{F\in\mathcal{F}^{K}_{h}}\sum_{e\in\mathcal{E}^{F}_{h}}(\mathbf{E}_{h}\cdot\mathbf{t}_{e},\mathbf{F}_{h}\cdot\mathbf{t}_{e})_{e}\quad\quad\forall\mathbf{E}_{h},\,\mathbf{F}_{h}\in\mathbf{V}_{h}^{\operatorname{edge}}(K).$

The summation term in $S^{K}_{\operatorname{edge}}(\cdot,\cdot)$ is multiplied
by $h_{K}^{2}$ to have a consistent scaling for both terms of the edge inner
product (24) with respect to the element size. The bilinear form (24) is
computable on all elements $K\in\Omega_{h}$ using the degrees of freedom of
$\mathbf{E}_{h},\mathbf{F}_{h}\in\mathbf{V}_{h}^{\operatorname{edge}}(K)$, and
has the two crucial properties of consistency and stability:

  * •

consistency: for all $\mathbf{q}\in\big{[}\mathbbm{P}_{0}(K)\big{]}^{3}$ and
all $\mathbf{E}_{h}\in\mathbf{V}_{h}^{\operatorname{edge}}(K)$,

$\displaystyle\big{[}\mathbf{q},\mathbf{E}_{h}\big{]}_{\operatorname{edge},K}=(\mathbf{q},\mathbf{E}_{h})_{K};$
(26)

  * •

stability: for all $\mathbf{E}_{h}\in\mathbf{V}_{h}^{\operatorname{edge}}(K)$,

$\displaystyle\eta_{*}\|\mathbf{E}_{h}\|^{2}_{0,K}\leq\big{[}\mathbf{E}_{h},\mathbf{E}_{h}\big{]}_{\operatorname{edge},K}\leq\eta^{*}\|\mathbf{E}_{h}\|^{2}_{0,K},$
(27)

with $\eta_{*}=\min(1,\gamma_{*})$ and $\eta^{*}=\max(1,\gamma^{*})$.

Finally, the global inner product over the virtual element edge space
$\mathbf{V}_{h}^{\operatorname{edge}}(K)$ is given by adding all the elemental
contributions:

$\displaystyle\big{[}\mathbf{E}_{h},\mathbf{F}_{h}\big{]}_{\operatorname{edge}}:=\sum_{K\in\Omega_{h}}\big{[}\mathbf{E}_{h},\mathbf{F}_{h}\big{]}_{\operatorname{edge},K}.$

In what follows, we shall use the following discrete norm:

$\|\cdot\|_{\mathbf{V}_{h}^{\operatorname{edge}}}^{2}:=\big{[}\cdot,\cdot\big{]}_{\operatorname{edge}}.$

###  5.3 The face space

We define the virtual element face space on an element $K\in\Omega_{h}$ as

$\displaystyle\mathbf{V}_{h}^{\operatorname{face}}(K):=\bigg{\\{}\mathbf{B}_{h}$
$\displaystyle\in[L^{2}(K)]^{3}\,\mid\,\mathbf{B}_{h}\cdot\mathbf{n}_{F}\in\mathbbm{P}_{0}(F)\,\forall
F\in\mathcal{F}^{K}_{h},$
$\displaystyle\begin{cases}\operatorname{div}\mathbf{B}_{h}\in\mathbbm{P}_{0}(K),\\\
\operatorname{\mathbf{curl}}\mathbf{B}_{h}\in\big{[}\mathbbm{P}_{0}(K)\big{]}^{3}\end{cases},\qquad\int_{K}\mathbf{B}_{h}\cdot(\mathbf{x}_{K}\times\mathbf{q})=0\
\forall\mathbf{q}\in\big{[}\mathbbm{P}_{0}(K)\big{]}^{3}\bigg{\\}},$

where $\mathbf{x}_{K}=\mathbf{x}-\mathbf{b}_{K}$.

Each element $\mathbf{B}_{h}\in\mathbf{V}_{h}^{\operatorname{face}}(K)$ is
uniquely determined by its normal components on the elemental faces, i.e., the
set of constant values of
$B_{F}=\big{(}\mathbf{B}_{h}\cdot\mathbf{n}_{F}\big{)}_{F\in\mathcal{F}^{K}_{h}}$.
We take these values as the degrees of freedom of $\mathbf{B}_{h}$ and refer
the reader to [11] for the proof of their unisolvence. A major consequence of
such a space definition is that $\bm{\Pi}^{0,K}_{0}\mathbf{B}_{h}$ and the
divergence of $\mathbf{B}_{h}$ are computable from the degrees of freedom of
$\mathbf{B}_{h}$ for all
$\mathbf{B}_{h}\in\mathbf{V}_{h}^{\operatorname{face}}(K)$.

We define the global face space $\mathbf{V}_{h}^{\operatorname{face}}$ by an
$H{}(\operatorname{div})$-coupling of the degrees of freedom:

$\mathbf{V}_{h}^{\operatorname{face}}:=\Big{\\{}\mathbf{B}_{h}\in
H_{0}(\operatorname{div},\Omega)\,\mid\,{\mathbf{B}_{h}}_{|K}\in\mathbf{V}_{h}^{\operatorname{face}}(K)\,\,\forall
K\in\Omega_{h}\Big{\\}}.$

This definition includes the homogeneous normal boundary conditions on
$\partial\Omega$.

####  5.3.1 Virtual element inner product

From [12, 13], we recall the discretization of the $L^{2}$-inner product for
the local face virtual element space, which reads as

$\big{[}\mathbf{B}_{h},\mathbf{C}_{h}\big{]}_{\operatorname{face},K}:=\Big{(}\bm{\Pi}^{0,K}_{0}\mathbf{B}_{h},\bm{\Pi}^{0,K}_{0}\mathbf{C}_{h}\Big{)}_{0,K}+S^{K}_{\operatorname{face}}((I-\bm{\Pi}^{0,K}_{0})\mathbf{B}_{h},(I-\bm{\Pi}^{0,K}_{0})\mathbf{C}_{h})\\\
\quad\forall\mathbf{B}_{h},\,\mathbf{C}_{h}\in\mathbf{V}_{h}^{\operatorname{face}}(K),$
(28)

where $S^{K}_{\operatorname{face}}(\cdot,\cdot)$ can be any computable,
symmetric bilinear form such that there exist two positive constants
$\widetilde{\gamma}_{*}$ and $\widetilde{\gamma}^{*}$ independent of $h_{K}$
satisfying

$\widetilde{\gamma}_{*}\|\mathbf{B}_{h}\|^{2}_{0,K}\leq
S^{K}_{\operatorname{face}}(\mathbf{B}_{h},\mathbf{B}_{h})\leq\widetilde{\gamma}^{*}\|\mathbf{B}_{h}\|^{2}_{0,K}\quad\quad\forall\mathbf{B}_{h}\in\ker(\bm{\Pi}^{0,K}_{0})\cap\mathbf{V}_{h}^{\operatorname{face}}(K).$
(29)

An explicit stabilization satisfying (29) was introduced in [11, formula
(4.8)] analyzed in [8, Proposition 5.2], and reads

$S^{K}_{\operatorname{face}}(\mathbf{B}_{h},\mathbf{C}_{h}):=h_{K}\sum_{F\in\mathcal{F}^{K}_{h}}\big{(}\mathbf{n}_{F}\cdot\mathbf{B}_{h},\mathbf{n}_{F}\cdot\mathbf{C}_{h}\big{)}_{0,F}\quad\quad\forall\mathbf{B}_{h},\,\mathbf{C}_{h}\in\mathbf{V}_{h}^{\operatorname{face}}(K).$

The bilinear form (28) is computable on all elements $K\in\Omega_{h}$ using
the degrees of freedom of
$\mathbf{B}_{h},\mathbf{C}_{h}\in\mathbf{V}_{h}^{\operatorname{face}}(K)$, and
has the two crucial properties of consistency and stability:

  * •

consistency: for all $\mathbf{q}\in\big{[}\mathbbm{P}_{0}(K)\big{]}^{3}$ and
all $\mathbf{B}_{h}\in\mathbf{V}_{h}^{\operatorname{face}}(K)$,

$\big{[}\mathbf{q},\mathbf{B}_{h}\big{]}_{\operatorname{face},K}=(\mathbf{q},\mathbf{B}_{h})_{K};$

  * •

stability: for all $\mathbf{B}_{h}\in\mathbf{V}_{h}^{\operatorname{face}}(K)$:

$\displaystyle\chi_{*}\|\mathbf{B}_{h}\|^{2}_{0,K}\leq\big{[}\mathbf{B}_{h},\mathbf{B}_{h}\big{]}_{\operatorname{face},K}\leq\chi^{*}\|\mathbf{B}_{h}\|^{2}_{0,K},$
(30)

with $\chi_{*}=\min(1,\widetilde{\gamma}_{*})$ and
$\chi^{*}=\max(1,\widetilde{\gamma}^{*})$.

Finally, the inner product over the global virtual element space
$\mathbf{V}_{h}^{\operatorname{face}}$ is given by summing all the elemental
contributions:

$\displaystyle\big{[}\mathbf{B}_{h},\mathbf{C}_{h}\big{]}_{\operatorname{face}}:=\sum_{K\in\Omega_{h}}\big{[}\mathbf{B}_{h},\mathbf{C}_{h}\big{]}_{\operatorname{face},K}.$

In what follows, we shall use the following discrete norm:

$\|\cdot\|_{\mathbf{V}_{h}^{\operatorname{face}}}^{2}:=\big{[}\cdot,\cdot\big{]}_{\operatorname{face}}.$

###  5.4 Interpolations and exact sequence properties

$H(\operatorname{\mathbf{curl}},\Omega)$$H(\operatorname{div},\Omega)$$\operatorname{\mathbf{curl}}$$\mathbf{V}_{h}^{\operatorname{edge}}$$\widetilde{\mathbf{V}}_{h}^{\operatorname{face}}$$\operatorname{\mathbf{curl}}$$\operatorname{\mathcal{I}^{\mathcal{E}_{h}}}$$\operatorname{\mathcal{I}^{\mathcal{F}_{h}}}$
Figure 1: Commuting diagram for the edge and face virtual element functions.

As thoroughly discussed, e.g., in [11], the above spaces satisfy fundamental
exact sequence properties. Introduce an edge approximation operator
$\operatorname{\mathcal{I}^{\mathcal{E}_{h}}}:H(\operatorname{\mathbf{curl}},\Omega)\to\mathbf{V}_{h}^{\operatorname{edge}}$
as follows. The edge approximant
$\operatorname{\mathcal{I}^{\mathcal{E}_{h}}}(\mathbf{E}):=\mathbf{E}_{I\\!I}\in\mathbf{V}_{h}^{\operatorname{edge}}$
of a given vector field $\mathbf{E}\in H(\operatorname{\mathbf{curl}},\Omega)$
is defined as the solution to the variational problem

$\displaystyle\big{[}\operatorname{\mathbf{curl}}\mathbf{E}_{I\\!I},\operatorname{\mathbf{curl}}\mathbf{F}_{h}\big{]}_{\operatorname{face}}$
$\displaystyle=\big{(}\operatorname{\mathbf{curl}}\mathbf{E},\operatorname{\mathbf{curl}}\mathbf{F}_{h}\big{)}$
$\displaystyle\forall\mathbf{F}_{h}\in\mathbf{V}_{h}^{\operatorname{edge}},$
(31a) $\displaystyle\big{[}\mathbf{E}_{I\\!I},\nabla
r_{h}\big{]}_{\operatorname{edge}}$ $\displaystyle=\big{(}\mathbf{E},\nabla
r_{h}\big{)}$ $\displaystyle\forall r_{h}\in V_{h}^{\operatorname{node}},$
(31b)

where the spaces $\mathbf{V}_{h}^{\operatorname{edge}}$ and
$V_{h}^{\operatorname{node}}$ are defined in (23) and (22), respectively. The
following approximation bound is valid [13, Proposition 5].

######  Lemma 5.1.

Let $\mathbf{E}\in H(\operatorname{\mathbf{curl}},\Omega)$ and consider the
approximant $\mathbf{E}_{I\\!I}\in\mathbf{V}_{h}^{\operatorname{edge}}$
defined as in (31). Then, there exists a real, positive constant $C_{I}$
independent of $h$ such that

$\|\mathbf{E}-\mathbf{E}_{I\\!I}\|_{\operatorname{\mathbf{curl}}}\leq
C_{I}h\big{(}|\mathbf{E}|_{1,\Omega}+\|\operatorname{\mathbf{curl}}\mathbf{E}\|_{0,\Omega}+h|\operatorname{\mathbf{curl}}\mathbf{E}|_{1,\Omega}\big{)}.$

Introduce a face approximation operator
$\operatorname{\mathcal{I}^{\mathcal{F}_{h}}}:H(\operatorname{div},\Omega)\to\mathbf{V}_{h}^{\operatorname{face}}$
as follows. The face approximant

$\operatorname{\mathcal{I}^{\mathcal{F}_{h}}}(\mathbf{B}):=\mathbf{B}_{I\\!I}\in\widetilde{\mathbf{V}}_{h}^{\operatorname{face}}:=\left\\{\mathbf{C}_{h}\in\mathbf{V}_{h}^{\operatorname{face}}\mid\operatorname{div}\mathbf{C}_{h}=0\right\\}$

of a given vector field $\mathbf{B}\in H(\operatorname{div},\Omega)$ is
defined as the solution to the variational problem

$\displaystyle\big{[}\mathbf{B}_{I\\!I},\mathbf{C}_{h}\big{]}_{\operatorname{face}}=\big{(}\mathbf{B},\mathbf{C}_{h}\big{)}\qquad\forall\mathbf{C}_{h}\in\widetilde{\mathbf{V}}_{h}^{\operatorname{face}}.$
(32)

We have that
$\widetilde{\mathbf{V}}_{h}^{\operatorname{face}}=\operatorname{\mathbf{curl}}(\mathbf{V}_{h}^{\operatorname{edge}})$;
see [11, equation (4.35)].

We recall the following approximation result [13, Lemma 7].

######  Lemma 5.2.

Let $\mathbf{B}\in[H^{1}(\Omega)]^{3}$ be a divergence-free magnetic field and
$\mathbf{B}_{I\\!I}\in\mathbf{V}_{h}^{\operatorname{face}}$ be its approximant
defined as in (32). Then, there exists a real, positive constant $C_{I}$
independent of $h$ such that

$\|\mathbf{B}-\mathbf{B}_{I\\!I}\|\leq C_{I}h\|\mathbf{B}\|_{1}.$

Additionally, we have a commuting property involving the operators defined in
(31) and (32), and the curl operator, which we report in the next lemma, whose
proof can be found in [13, Proposition 6]. The interpolation diagram is shown
in Figure 1.

######  Lemma 5.3.

Let $\mathbf{E}_{I\\!I}\in\mathbf{V}_{h}^{\operatorname{edge}}$ and
$(\operatorname{\mathbf{curl}}\mathbf{E})_{I\\!I}\in\mathbf{V}_{h}^{\operatorname{face}}$
be the approximant of $\mathbf{E}$ and
$\operatorname{\mathbf{curl}}\mathbf{E}\in H(\operatorname{div},\Omega)$
defined in (31) and (32), respectively. Then, the following commuting property
is valid:

$\operatorname{\mathbf{curl}}(\mathbf{E}_{I\\!I})=(\operatorname{\mathbf{curl}}\mathbf{E})_{I\\!I}\quad\quad\forall\mathbf{E}\in
H(\operatorname{\mathbf{curl}},\Omega).$

###  5.5 The semi-discrete MHD model

First, note that
$-(\mathbf{j}\times\mathbf{B},\mathbf{v})=(\mathbf{j},\mathbf{v}\times\mathbf{B})$
in equation (9a). To discretize this term, we consider the elemental bilinear
operator
${\bm{\chi}}_{h,K}:\mathbf{W}_{h}(K)\times\mathbf{V}_{h}^{\operatorname{face}}(K)\to\mathbf{V}_{h}^{\operatorname{edge}}(K)$.
This operator denotes the discretization of the cross product
“$\mathbf{v}_{h}\times\mathbf{B}_{h}$” for
$\mathbf{v}_{h}\in\mathbf{W}_{h}(K)$ and
$\mathbf{B}_{h}\in\mathbf{V}_{h}^{\operatorname{face}}(K)$ through the
elemental edge function
${\bm{\chi}}_{h,K}(\mathbf{v}_{h},\mathbf{B}_{h})\in\mathbf{V}_{h}^{\operatorname{edge}}(K)$.
For the sake of notation, we also introduce the global operator

$\displaystyle{\bm{\chi}}_{h}:\mathbf{W}_{h}\times\mathbf{V}_{h}^{\operatorname{face}}\to\Pi_{K\in\Omega_{h}}\mathbf{V}_{h}^{\operatorname{edge}}(K)\quad\textrm{such
that}\quad{{\bm{\chi}}_{h}(\mathbf{v}_{h},\mathbf{B}_{h})}_{|K}={\bm{\chi}}_{h,K}({\mathbf{v}_{h}}_{|K},{\mathbf{B}_{h}}_{|K})$

for the global fields $\mathbf{v}_{h}\in\mathbf{W}_{h}$ and
$\mathbf{B}_{h}\in\mathbf{V}_{h}^{\operatorname{face}}$. We emphasize that
${\bm{\chi}}_{h}(\mathbf{v}_{h},\mathbf{B}_{h})$ does not belong to
$\mathbf{V}_{h}^{\operatorname{edge}}$ since it is not
$H(\operatorname{\mathbf{curl}},\Omega)$-conforming. However, to simplify
several equations below, with an abuse of notation, we let

$\displaystyle\big{[}\mathbf{F}_{h},{\bm{\chi}}_{h}(\mathbf{v}_{h},\mathbf{B}_{h})\big{]}_{\operatorname{edge}}:=\sum_{K\in\Omega_{h}}\big{[}\mathbf{F}_{h},{\bm{\chi}}_{h}(\mathbf{v}_{h},\mathbf{B}_{h})\big{]}_{\operatorname{edge},K}.$

Using this bilinear operator, we define the “broken” virtual element edge
function

$\displaystyle\mathbf{j}_{h}=\mathbf{E}_{h}+{\bm{\chi}}_{h}(\mathbf{u}_{h},\mathbf{B}_{h}),$

which can be considered as the electric current density of the discrete MHD
model. For the sake of simplicity, we take $s=1$ in (6a) and (9a).

Several choices for ${\bm{\chi}}_{h}$ are possible. Among them, we pick

${\bm{\chi}}_{h}(\mathbf{v}_{h},\mathbf{B}_{h})_{|K}:=\Pi^{0,K}_{0}\mathbf{v}_{h}\times\Pi^{0,K}_{0}\mathbf{B}_{h}.$
(33)

The semi-discrete virtual element method for the MHD equations reads as: _For
every $t\in(0,T)$, find
$(\mathbf{u}_{h},p_{h},\mathbf{E}_{h},\mathbf{B}_{h})\in\mathbf{W}_{h}\times
Q_{h}\times\mathbf{V}_{h}^{\operatorname{edge}}\times\mathbf{V}_{h}^{\operatorname{face}}$
such that_

$\displaystyle
m_{h}(\mathbf{u}_{h,t},\mathbf{v}_{h})+Re^{-1}a_{h}(\mathbf{u}_{h},\mathbf{v}_{h})+b(\mathbf{v}_{h},p_{h})$
$\displaystyle+\widetilde{c}_{h}(\mathbf{u}_{h};\mathbf{u}_{h},\mathbf{v}_{h})+\big{[}\mathbf{j}_{h},{\bm{\chi}}_{h}(\mathbf{v}_{h},\mathbf{B}_{h})\big{]}_{\operatorname{edge}}$
$\displaystyle\qquad\qquad\qquad=(\mathbf{f},\Pi^{0}_{1}\mathbf{v}_{h})$
$\displaystyle\phantom{=\leavevmode\nobreak\
0}\qquad\forall\mathbf{v}_{h}\in\mathbf{W}_{h}\subset\big{[}H^{1}_{0}(\Omega)\big{]}^{3},$
(34a)
$\displaystyle\big{[}\mathbf{j}_{h},\mathbf{F}_{h}\big{]}_{\operatorname{edge}}-Re_{m}^{-1}\big{[}\mathbf{B}_{h},\operatorname{\mathbf{curl}}\mathbf{F}_{h}\big{]}_{\operatorname{face}}$
$\displaystyle=0\qquad\forall\mathbf{F}_{h}\in\mathbf{V}_{h}^{\operatorname{edge}}\subset
H_{0}(\operatorname{\mathbf{curl}},\Omega),$ (34b)
$\displaystyle\big{[}\mathbf{B}_{h,t},\mathbf{C}_{h}\big{]}_{\operatorname{face}}+\big{[}\operatorname{\mathbf{curl}}\mathbf{E}_{h},\mathbf{C}_{h}\big{]}_{\operatorname{face}}$
$\displaystyle=0\qquad\forall\mathbf{C}_{h}\in\mathbf{V}_{h}^{\operatorname{face}}\subset
H_{0}(\operatorname{div},\Omega),$ (34c) $\displaystyle
b(\mathbf{u}_{h},q_{h})$ $\displaystyle=0\qquad\forall q_{h}\in Q_{h}\subset
L_{0}(\Omega).$ (34d)

We provide the semi-discrete scheme (34) with the discrete initial conditions
$\mathbf{u}_{h}(0)=\mathbf{u}_{h,0}$ and $\mathbf{B}_{h}(0)=\mathbf{B}_{h,0}$,
where $\mathbf{u}_{h,0}$ and $\mathbf{B}_{h,0}$ are the degrees of freedom
interpolants in $\mathbf{W}_{h}$ and $\mathbf{V}_{h}^{\operatorname{face}}$ of
the continuous initial conditions $\mathbf{u}_{0}$ and $\mathbf{B}_{0}$ in
(7). Provided that $\mathbf{B}$ is sufficiently regular, from [8, Proposition
$4.5$], the degrees of freedom interpolant of $\mathbf{B}$ in the face space
$\mathbf{V}_{h}^{\operatorname{face}}(K)$ satisfies

$\|\mathbf{B}(0)-\mathbf{B}_{h,0}\|_{0,K}\lesssim h_{K}.$ (35)

######  Remark 5.4.

Recalling that $\operatorname{div}\mathbf{u}_{h}\in Q_{h}$ and
$\operatorname{\mathbf{curl}}\mathbf{E}_{h}\in\mathbf{V}_{h}^{\operatorname{face}}$,
the third and fourth equations provide

$\mathbf{B}_{h,t}+\operatorname{\mathbf{curl}}\mathbf{E}_{h}=\mathbf{0},\qquad\operatorname{div}\mathbf{u}_{h}=0\qquad\text{in
}\Omega.$ (36)

As a consequence, by setting the initial data such that
$\operatorname{div}\mathbf{B}_{h,0}=0$, the first equation in (36) implies
$\operatorname{div}\mathbf{B}_{h}=0$ at all times. Thus, the proposed scheme
satisfies both solenoidal constraints exactly.

We conclude this section by presenting the following stability result on the
solution to the semi-discrete scheme (34) and briefly commenting on the
existence of a (unique) solution to (34).

######  Proposition 5.5.

Let $\mathbf{u}_{h}$, $\mathbf{E}_{h}$, and $\mathbf{B}_{h}$ be the solutions
to (34). Then, the following stability estimates are valid:

$\begin{split}\frac{1}{2}\frac{\text{d}}{\text{dt}}m_{h}(\mathbf{u}_{h},\mathbf{u}_{h})+\frac{1}{2}Re_{m}^{-1}\frac{\text{d}}{\text{dt}}\big{[}\mathbf{B}_{h},\mathbf{B}_{h}\big{]}_{\operatorname{face}}+Re^{-1}a_{h}(\mathbf{u}_{h},\mathbf{u}_{h})+\big{[}\mathbf{j}_{h},\mathbf{j}_{h}\big{]}_{\operatorname{edge}}=(\mathbf{f},\Pi^{0}_{1}\mathbf{u}_{h})\quad\forall
t\in(0,T)\end{split}$

and

$\begin{split}&\max_{0\leq t\leq
T}\left(m_{h}(\mathbf{u}_{h},\mathbf{u}_{h})+Re_{m}^{-1}\big{[}\mathbf{B}_{h},\mathbf{B}_{h}\big{]}_{\operatorname{face}}\right)+\text{Re}^{-1}\int_{0}^{T}a_{h}(\mathbf{u}_{h},\mathbf{u}_{h})+2\int_{0}^{T}\big{[}\mathbf{j}_{h},\mathbf{j}_{h}\big{]}_{\operatorname{edge}}\\\
&\quad\leq
m_{h}(\mathbf{u}_{h,0},\mathbf{u}_{h,0})+Re_{m}^{-1}\big{[}\mathbf{B}_{h,0},\mathbf{B}_{h,0}\big{]}_{\operatorname{face}}+\alpha_{*}^{-1}C_{P}Re\int_{0}^{T}\|\mathbf{f}\|_{0}.\end{split}$

###### Proof.

The first estimate follows picking $\mathbf{v}_{h}=\mathbf{u}_{h}$ in (34a),
$\mathbf{F}_{h}=\mathbf{E}_{h}$ in (34b), and $\mathbf{C}_{h}=\mathbf{B}_{h}$
in (34c), multiplying (34b) by $Re_{m}$, and summing up the three resulting
equations.

As for the second one, we observe that

$\begin{split}(\mathbf{f},\Pi^{0}_{1}\mathbf{u}_{h})&\leq\frac{1}{2}\alpha_{*}^{-1}C_{P}Re\|\mathbf{f}\|_{0}^{2}+\frac{1}{2}\alpha_{*}C_{P}^{-1}Re^{-1}\|\Pi^{0}_{1}\mathbf{u}_{h}\|^{2}\leq\frac{1}{2}\alpha_{*}^{-1}C_{P}Re\|\mathbf{f}\|_{0}^{2}+\frac{1}{2}\alpha_{*}Re^{-1}|\mathbf{u}_{h}|_{1}^{2}\\\
&\leq\frac{1}{2}\alpha_{*}^{-1}C_{P}Re\|\mathbf{f}\|_{0}^{2}+\frac{1}{2}Re^{-1}a_{h}(\mathbf{u}_{h},\mathbf{u}_{h}).\end{split}$

Next, we absorb the second term on the right-hand side in the left-hand side
of the first bound. The assertion follows by integrating in time the first
bound between $0$ and any time $0<t\leq T$, multiplying by two both sides, and
taking the maximum over the time $t$. ∎

The existence of a unique solution to the semi-discrete problem (34) follows
from the standard ODE theory in finite dimensions, in the spirit of Section
9.2 of [57]. Here, we only sketch the main steps of the proof. First, we
eliminate the pressure variable and equation $\eqref{eq:VEM:weak:D}$ by
restricting the (test and trial) velocity space to divergence-free functions.
Next, (34b) at the initial time $t=0$ together with the definition of
$\mathbf{j}_{h}$ allows us to define an initial condition $\mathbf{E}_{h,0}$
for the electric field, depending on the initial conditions $\mathbf{u}_{h,0}$
and $\mathbf{B}_{h,0}$ of the velocity and magnetic fields. Then, we can
substitute $\eqref{eq:VEM:weak:B}$ with the same equation derived in time;
upon further applying (34c), we obtain

$\big{[}{\mathbf{E}}_{h,t},\mathbf{F}_{h}\big{]}_{\operatorname{edge}}+Re_{m}^{-1}\big{[}\operatorname{\mathbf{curl}}\mathbf{E}_{h},\operatorname{\mathbf{curl}}\mathbf{F}_{h}\big{]}_{\operatorname{face}}+\big{[}\partial_{t}{\bm{\chi}}_{h}(\mathbf{v}_{h},\mathbf{B}_{h}),\mathbf{F}_{h}\big{]}_{\operatorname{edge}}=0$
(37)

for all $\mathbf{F}_{h}$ in $\mathbf{V}_{h}^{\operatorname{edge}}\subset
H_{0}(\operatorname{\mathbf{curl}},\Omega)$. Note that the term

$\displaystyle\big{[}\partial_{t}{\bm{\chi}}_{h}(\mathbf{v}_{h},\mathbf{B}_{h}),\mathbf{F}_{h}\big{]}_{\operatorname{edge}}=\big{[}\Pi^{0,K}_{0}\mathbf{u}_{h,t}\times\Pi^{0,K}_{0}\mathbf{B}_{h},\mathbf{F}_{h}\big{]}_{\operatorname{edge}}+\big{[}\Pi^{0,K}_{0}\mathbf{u}_{h}\times\Pi^{0,K}_{0}\mathbf{B}_{h,t},\mathbf{F}_{h}\big{]}_{\operatorname{edge}}$

can be interpreted as a (quartic) function of
$(\mathbf{u}_{h},\mathbf{B}_{h},\mathbf{E}_{h})$ using the algebraic form of
equations (34a), (34c) and inverting the corresponding “mass” matrices.
Therefore, substituting (34b) with equation (37), and recalling that ${\bf
E}_{h,0}$ is known, yields an equivalent system. The final system is a first
order Cauchy problem on the finite-dimensional variable
$(\mathbf{u}_{h},\mathbf{B}_{h},\mathbf{E}_{h})$, the nonlinearity being of
quartic type. Thanks to Proposition 5.5, asserting that the potential discrete
solutions to the system cannot blow up in finite time, the nonlinearity is
uniformly Lipschitz. The existence and uniqueness result (for finite times)
follows from the standard ODE theory in finite dimension.

##  6 Convergence analysis of the virtual element approximation

Introduce the “discrete errors”

$\displaystyle\mathbf{e}_{h}^{\mathbf{u}}:=\mathbf{u}_{h}-\mathbf{u}_{I},\qquad\mathbf{e}_{h}^{\mathbf{E}}:=\mathbf{E}_{h}-\mathbf{E}_{I\\!I},\qquad\mathbf{e}_{h}^{\mathbf{B}}:=\mathbf{B}_{h}-\mathbf{B}_{I\\!I},$
(38)

and the virtual element edge fields

$\displaystyle\mathbf{j}_{I}$
$\displaystyle=\mathbf{E}_{I\\!I}+{\bm{\chi}}_{h}(\mathbf{u}_{I},\mathbf{B}_{h})$
(39)

and

$\mathbf{e}_{h}^{\mathbf{j}}=\mathbf{j}_{h}-\mathbf{j}_{I}=\big{(}\mathbf{E}_{h}+{\bm{\chi}}_{h}(\mathbf{u}_{h},\mathbf{B}_{h})\big{)}-\big{(}\mathbf{E}_{I\\!I}+{\bm{\chi}}_{h}(\mathbf{u}_{I},\mathbf{B}_{h})\big{)}=\mathbf{e}_{h}^{\mathbf{E}}+{\bm{\chi}}_{h}(\mathbf{e}_{h}^{\mathbf{u}},\mathbf{B}_{h}).$
(40)

These quantities are bounded as in the following theorem, whose proof can be
found in Section 6.1.

######  Theorem 6.1.

Let $(\mathbf{u}(t),p(t)$, $\mathbf{E}(t),\mathbf{B}(t))$ in
$\big{[}H^{1}_{0}(\Omega)\big{]}^{3}\times L^{2}_{0}(\Omega)\times
H_{0}(\operatorname{\mathbf{curl}},\Omega)\times
H_{0}(\operatorname{div},\Omega)$ for almost every $t\in(0,T)$ be the solution
to the variational formulation (9) under the assumptions in Section 3.
Furthermore, assume that the following regularity properties

$\displaystyle\mathbf{u}\in[W^{2,\infty}(\Omega)]^{3},\qquad\mathbf{u}_{t}\in[H^{1}(\Omega)]^{3},\qquad\mathbf{B}\in[H^{1}(\Omega)\cap
L^{\infty}(\Omega)]^{3}\quad\textrm{and}\quad\mathbf{j}\in[L^{\infty}(\Omega)]^{3}$
(41)

hold uniformly in time.

Let $(\mathbf{u}_{h}(t),p_{h}(t),\mathbf{E}_{h}(t),\mathbf{B}_{h}(t))$ in
$\mathbf{W}_{h}\times
Q_{h}\times\mathbf{V}_{h}^{\operatorname{edge}}\times\mathbf{V}_{h}^{\operatorname{face}}$
for every $t\in(0,T)$ be the solution to the virtual element method (34).
Then, a positive constant $C$ exists that is independent of $h$ such that

$\|\mathbf{e}_{h}^{\mathbf{u}}(t)\|+\|\mathbf{e}_{h}^{\mathbf{B}}(t)\|+\left(\int_{0}^{t}\|\mathbf{e}_{h}^{\mathbf{j}}(s)\|^{2}\right)^{\frac{1}{2}}\leq
C\big{(}\|\mathbf{e}_{h}^{\mathbf{u}}(0)\|+\|\mathbf{e}_{h}^{\mathbf{B}}(0)\|+h\big{)}\qquad\text{for
every }t\in(0,T],$

where $\mathbf{e}_{h}^{\mathbf{u}}$, $\mathbf{e}_{h}^{\mathbf{B}}$, and
$\mathbf{e}_{h}^{\mathbf{j}}$ are defined in (38)-(40). The constant $C$
depends on the parameters of the discretization, the final time $T$, and the
regularity of the MHD solution fields
$\mathbf{u}(t),\mathbf{E}(t),\mathbf{B}(t)$, and $\mathbf{B}_{h}(t)$.

In view of Theorem 6.1, we have the following convergence result.

######  Theorem 6.2.

Let $(\mathbf{u}(t),p(t)$, $\mathbf{E}(t),\mathbf{B}(t))$ in
$\big{[}H^{1}_{0}(\Omega)\big{]}^{3}\times L^{2}_{0}(\Omega)\times
H_{0}(\operatorname{\mathbf{curl}},\Omega)\times
H_{0}(\operatorname{div},\Omega)$ for almost every $t\in(0,T)$ be the solution
to the variational formulation (9) under the assumptions in Section 3. Let the
regularity assumptions of Theorem 6.1 be valid and assume that the initial
conditions $\mathbf{u}_{0}$, $\mathbf{B}_{0}$, and $\mathbf{E}_{0}$ of the
fields $\mathbf{u}$, $\mathbf{B}$, and $\mathbf{E}$ belong to
$[H^{1}(\Omega)]^{3}$. Let
$(\mathbf{u}_{h}(t),p_{h}(t),\mathbf{E}_{h}(t),\mathbf{B}_{h}(t))$ in
$\mathbf{W}_{h}\times
Q_{h}\times\mathbf{V}_{h}^{\operatorname{edge}}\times\mathbf{V}_{h}^{\operatorname{face}}$
for every $t\in(0,T)$ be the solution to the virtual element method (34).
Then, a positive constant $C$ exists, independent of $h$, such that

$\|\mathbf{u}(t)-\mathbf{u}_{h}(t)\|+\|\mathbf{B}(t)-\mathbf{B}_{h}(t)\|+\left(\int_{0}^{t}\|\mathbf{E}(s)-\mathbf{E}_{h}(s)\|^{2}\right)^{\frac{1}{2}}\leq
Ch$ (42)

for almost every $t\in(0,T]$. The constant $C$ depends on the parameters of
the discretization, the final time $T$, and the regularity of the MHD solution
fields $\mathbf{u}(t),\mathbf{E}(t),\mathbf{B}(t)$, and $\mathbf{B}_{h}(t)$.

###### Proof.

We add and subtract $\mathbf{u}_{I}$ (see Lemma 4.2), $\mathbf{E}_{I\\!I}$
(see (31)), and $\mathbf{B}_{I\\!I}$ (see (32)) to the left-hand side of (42)
and get

$\begin{split}&\|\mathbf{u}(t)-\mathbf{u}_{h}(t)\|+\|\mathbf{B}(t)-\mathbf{B}_{h}(t)\|+\left(\int_{0}^{t}\|\mathbf{E}(s)-\mathbf{E}_{h}(s)\|^{2}\right)^{\frac{1}{2}}\\\
&\leq\left[\|\mathbf{u}(t)-\mathbf{u}_{I}(t)\|+\|\mathbf{B}(t)-\mathbf{B}_{I\\!I}(t)\|+\left(\int_{0}^{t}\|\mathbf{E}(s)-\mathbf{E}_{I\\!I}(s)\|^{2}\right)^{\frac{1}{2}}\right]\\\
&\quad+\left[\|\mathbf{e}_{h}^{\mathbf{u}}(t)\|+\|\mathbf{e}_{h}^{\mathbf{B}}(t)\|+\left(\int_{0}^{t}\|\mathbf{e}_{h}^{\mathbf{E}}(s)\|^{2}\right)^{\frac{1}{2}}\right]:=S_{1}+S_{2}.\end{split}$

An upper bound on the term $S_{1}$ can be shown using the interpolation
properties of Lemmas 4.2, 5.2, and 5.1.

Observe that $\mathbf{e}_{h}^{\mathbf{u}}(0)=0$. Furthermore, the triangle
inequality, (35), and Lemma 5.2 give

$\|\mathbf{e}_{h}^{\mathbf{B}}(0)\|\leq\|\mathbf{B}(0)-\mathbf{B}_{h,0}\|+\|\mathbf{B}(0)-\mathbf{B}_{I\\!I}(0)\|\lesssim
h.$

With this at hand, Theorem 6.1 allows us to show a direct bound on the terms
involving the velocity and magnetic fields appearing in $S_{2}$. For the term
involving the electric field we proceed as follows: using the triangle
inequality, the definition of ${\bm{\chi}}_{h}$ in (33), the bound on
$\|\mathbf{e}_{h}^{\mathbf{u}}\|$ from (41), the stability estimates on the
discrete face element bilinear forms (30), and the bound on
$\big{[}\mathbf{B}_{h},\mathbf{B}_{h}\big{]}_{\operatorname{face}}$ in
Proposition 5.5, we arrive at

$\begin{split}&\left(\int_{0}^{t}\|\mathbf{e}_{h}^{\mathbf{E}}(s)\|^{2}\right)^{\frac{1}{2}}\lesssim\left(\int_{0}^{t}\left(\|\mathbf{e}_{h}^{\mathbf{j}}(s)\|^{2}+\|{\bm{\chi}}_{h}(\mathbf{e}_{h}^{\mathbf{u}}(s),\mathbf{B}_{h}(s))\|^{2}\right)\right)^{\frac{1}{2}}\\\
&\quad\leq
h+\left(\int_{0}^{t}\|\mathbf{e}_{h}^{\mathbf{u}}(s)\|^{2}\|\mathbf{B}_{h}(s)\|^{2}\right)^{\frac{1}{2}}\lesssim
h+\|\mathbf{e}_{h}^{\mathbf{u}}\|_{L^{\infty}(0,t;L^{2}(\Omega))}\|\mathbf{B}_{h}\|_{L^{\infty}(0,t;\mathbf{V}_{h}^{\operatorname{face}})}\lesssim
h.\end{split}$

The final hidden constant depends on the square root of the instant time $t$.
∎

###  6.1 Proof of Theorem 6.1

###### Proof.

Throughout, for presentation’s sake, we assume that $Re_{m}=1$ and $Re=1$. The
general assertion can be proved by means of minor modifications but note that
the final constant in Theorem 6.1 will depend on $Re_{m}$ and $Re$; treating
convection-dominated cases is beyond the scope of this work.

First, we add and subtract $\mathbf{E}_{I\\!I}$ and
${\bm{\chi}}_{h}(\mathbf{u}_{I},\mathbf{B}_{h})$ to
$\mathbf{j}_{h}=\mathbf{E}_{h}+{\bm{\chi}}_{h}(\mathbf{u}_{h},\mathbf{B}_{h})$
so that

$\displaystyle\mathbf{j}_{h}=\mathbf{E}_{h}+{\bm{\chi}}_{h}(\mathbf{u}_{h},\mathbf{B}_{h})=(\mathbf{E}_{h}-\mathbf{E}_{I\\!I})+{\bm{\chi}}_{h}(\mathbf{u}_{h}-\mathbf{u}_{I},\mathbf{B}_{h})+\mathbf{E}_{I\\!I}+{\bm{\chi}}_{h}(\mathbf{u}_{I},\mathbf{B}_{h}),$
(43)

where $\mathbf{E}_{I\\!I}$ and $\mathbf{u}_{I}$ are defined in (31) and
Section 4.3.1 (see also Lemma 4.2), respectively. We substitute (43) in (34a),
add and subtract $\mathbf{u}_{I}$, and find that

$\displaystyle
m_{h}\big{(}(\mathbf{u}_{h}-\mathbf{u}_{I})_{t},\mathbf{v}_{h}\big{)}+a_{h}(\mathbf{u}_{h}-\mathbf{u}_{I},\mathbf{v}_{h})+b(\mathbf{v}_{h},p_{h})$
$\displaystyle\hskip
170.71652pt+\big{[}\mathbf{E}_{h}-\mathbf{E}_{I\\!I}+{\bm{\chi}}_{h}(\mathbf{u}_{h}-\mathbf{u}_{I}),{\bm{\chi}}_{h}(\mathbf{v}_{h},\mathbf{B}_{h})\big{]}_{\operatorname{edge}}$
$\displaystyle\qquad=(\mathbf{f},\Pi^{0}_{1}\mathbf{v}_{h})-m_{h}(\mathbf{u}_{I,t},\mathbf{v}_{h}\big{)}-a_{h}(\mathbf{u}_{I},\mathbf{v}_{h})-\widetilde{c}_{h}(\mathbf{u}_{h};\mathbf{u}_{h},\mathbf{v}_{h})$
$\displaystyle\qquad\phantom{=}-\big{[}\mathbf{E}_{I\\!I}+{\bm{\chi}}_{h}(\mathbf{u}_{I},\mathbf{B}_{h}),{\bm{\chi}}_{h}(\mathbf{v}_{h},\mathbf{B}_{h})\big{]}_{\operatorname{edge}}.$
(44)

Similarly, we use (43) in (34b), add and subtract $\mathbf{B}_{I\\!I}$ defined
in (32), and find that

$\big{[}\mathbf{E}_{h}-\mathbf{E}_{I\\!I}+{\bm{\chi}}_{h}(\mathbf{u}_{h}-\mathbf{u}_{I},\mathbf{B}_{h}),\mathbf{F}_{h}\big{]}_{\operatorname{edge}}-\big{[}\mathbf{B}_{h}-\mathbf{B}_{I\\!I},\operatorname{\mathbf{curl}}\mathbf{F}_{h}\big{]}_{\operatorname{face}}\\\\[5.0pt]
=-\big{[}\mathbf{E}_{I\\!I}+{\bm{\chi}}_{h}(\mathbf{u}_{I},\mathbf{B}_{h}),\mathbf{F}_{h}\big{]}_{\operatorname{edge}}+\big{[}\mathbf{B}_{I\\!I},\operatorname{\mathbf{curl}}\mathbf{F}_{h}\big{]}_{\operatorname{face}}.$
(45)

Then, we use Lemma 5.3, (32) twice and (9c), and find that, for all
$\mathbf{C}_{h}\in\widetilde{\mathbf{V}}_{h}^{\operatorname{face}}$,

$\displaystyle\big{[}\operatorname{\mathbf{curl}}\mathbf{E}_{I\\!I},\mathbf{C}_{h}\big{]}_{\operatorname{face}}=\big{[}(\operatorname{\mathbf{curl}}\mathbf{E})_{I\\!I},\mathbf{C}_{h}\big{]}_{\operatorname{face}}=(\operatorname{\mathbf{curl}}\mathbf{E},\mathbf{C}_{h})=-(\mathbf{B}_{t},\mathbf{C}_{h})=-\big{[}(\mathbf{B}_{t})_{I\\!I},\mathbf{C}_{h}\big{]}_{\operatorname{face}}.$

Since the derivation in time and the interpolation commute, e.g.,
$(\mathbf{B}_{t})_{I\\!I}=\mathbf{B}_{I\\!I,t}$, we find

$\displaystyle\big{[}\mathbf{B}_{I\\!I,t},\mathbf{C}_{h}\big{]}_{\operatorname{face}}+\big{[}\operatorname{\mathbf{curl}}\mathbf{E}_{I\\!I},\mathbf{C}_{h}\big{]}_{\operatorname{face}}=0.$

Finally, we subtract this equation to (34c) and obtain, for all
$\mathbf{C}_{h}\in\widetilde{\mathbf{V}}_{h}^{\operatorname{face}}$,

$\displaystyle\big{[}(\mathbf{B}_{h}-\mathbf{B}_{I\\!I})_{t},\mathbf{C}_{h}\big{]}_{\operatorname{face}}+\big{[}\operatorname{\mathbf{curl}}(\mathbf{E}_{h}-\mathbf{E}_{I\\!I}),\mathbf{C}_{h}\big{]}_{\operatorname{face}}=0.$
(46)

Next, we take $\mathbf{v}_{h}=\mathbf{e}_{h}^{\mathbf{u}}$ in (44),
$\mathbf{F}_{h}=\mathbf{e}_{h}^{\mathbf{E}}$ in (45), and
$\mathbf{C}_{h}=\mathbf{e}_{h}^{\mathbf{B}}$ in (46). This choice of
$\mathbf{C}_{h}$ is admissible since
$\operatorname{div}\mathbf{e}_{h}^{\mathbf{B}}=0$ for every $t\geq 0$ thanks
to Remark 5.4 and the fact that $\mathbf{B}_{I\\!I}$ has zero divergence by
definition. We further observe $b(\mathbf{e}_{h}^{\mathbf{u}},p_{h})=0$ since

$\displaystyle\mathbf{e}_{h}^{\mathbf{u}}\in\mathbf{Z}_{h}:=\big{\\{}\mathbf{v}_{h}\in\mathbf{W}_{h}\,\mid\,(\operatorname{div}\mathbf{v}_{h},q_{h})=0\quad\forall
q_{h}\in Q_{h}\big{\\}};$

see also (16). Equations (44), (45), and (46) thus become

$\displaystyle
m_{h}\big{(}\mathbf{e}^{\mathbf{u}}_{h,t},\mathbf{e}_{h}^{\mathbf{u}}\big{)}+a_{h}(\mathbf{e}_{h}^{\mathbf{u}},\mathbf{e}_{h}^{\mathbf{u}})+\big{[}\mathbf{e}_{h}^{\mathbf{E}}+{\bm{\chi}}_{h}(\mathbf{e}_{h}^{\mathbf{u}},\mathbf{B}_{h}),{\bm{\chi}}_{h}(\mathbf{e}_{h}^{\mathbf{u}},\mathbf{B}_{h})\big{]}_{\operatorname{edge}}$
$\displaystyle\qquad=(\mathbf{f},\Pi^{0}_{1}\mathbf{e}_{h}^{\mathbf{u}})-m_{h}(\mathbf{u}_{I,t},\mathbf{e}_{h}^{\mathbf{u}}\big{)}-a_{h}(\mathbf{u}_{I},\mathbf{e}_{h}^{\mathbf{u}})-\widetilde{c}_{h}(\mathbf{u}_{h};\mathbf{u}_{h},\mathbf{e}_{h}^{\mathbf{u}})$
$\displaystyle\qquad\phantom{=}-\big{[}\mathbf{E}_{I\\!I}+{\bm{\chi}}_{h}(\mathbf{u}_{I},\mathbf{B}_{h}),{\bm{\chi}}_{h}(\mathbf{e}_{h}^{\mathbf{u}},\mathbf{B}_{h})\big{]}_{\operatorname{edge}};$
(47)
$\displaystyle\big{[}\mathbf{e}_{h}^{\mathbf{E}}+{\bm{\chi}}_{h}(\mathbf{e}_{h}^{\mathbf{u}},\mathbf{B}_{h}),\mathbf{e}_{h}^{\mathbf{E}}\big{]}_{\operatorname{edge}}-\big{[}\mathbf{e}_{h}^{\mathbf{B}},\operatorname{\mathbf{curl}}\mathbf{e}_{h}^{\mathbf{E}}\big{]}_{\operatorname{face}}=-\big{[}\big{(}\mathbf{E}_{I\\!I}+{\bm{\chi}}_{h}(\mathbf{u}_{I},\mathbf{B}_{h})\big{)},\mathbf{e}_{h}^{\mathbf{E}}\big{]}_{\operatorname{edge}}$
$\displaystyle\phantom{\big{[}\mathbf{e}_{h}^{\mathbf{E}}+{\bm{\chi}}_{h}(\mathbf{e}_{h}^{\mathbf{u}},\mathbf{B}_{h}),\mathbf{e}_{h}^{\mathbf{E}}\big{]}_{\operatorname{edge}}-\big{[}\mathbf{e}_{h}^{\mathbf{B}},\operatorname{\mathbf{curl}}\mathbf{e}_{h}^{\mathbf{E}}\big{]}_{\operatorname{face}}=}+\big{[}\mathbf{B}_{I\\!I},\operatorname{\mathbf{curl}}\mathbf{e}_{h}^{\mathbf{E}}\big{]}_{\operatorname{face}};$
(48)
$\displaystyle\big{[}\mathbf{e}^{\mathbf{B}}_{h,t},\mathbf{e}_{h}^{\mathbf{B}}\big{]}_{\operatorname{face}}+\big{[}\operatorname{\mathbf{curl}}\mathbf{e}_{h}^{\mathbf{E}},\mathbf{e}_{h}^{\mathbf{B}}\big{]}_{\operatorname{face}}=0.$
(49)

We add (49) to (48), so that

$\displaystyle\big{[}\mathbf{e}_{h}^{\mathbf{E}}+{\bm{\chi}}_{h}(\mathbf{e}_{h}^{\mathbf{u}},\mathbf{B}_{h}),\mathbf{e}_{h}^{\mathbf{E}}\big{]}_{\operatorname{edge}}+\big{[}\mathbf{e}^{\mathbf{B}}_{h,t},\mathbf{e}_{h}^{\mathbf{B}}\big{]}_{\operatorname{face}}$
$\displaystyle=-\big{[}\big{(}\mathbf{E}_{I\\!I}+{\bm{\chi}}_{h}(\mathbf{u}_{I},\mathbf{B}_{h})\big{)},\mathbf{e}_{h}^{\mathbf{E}}\big{]}_{\operatorname{edge}}$
$\displaystyle\phantom{=\,}+\big{[}\mathbf{B}_{I\\!I},\operatorname{\mathbf{curl}}\mathbf{e}_{h}^{\mathbf{E}}\big{]}_{\operatorname{face}}.$

Then, we add this equation to (47) and write the resulting equation as

LHS
$\displaystyle:=m_{h}(\mathbf{e}^{\mathbf{u}}_{h,t},\mathbf{e}_{h}^{\mathbf{u}}\big{)}+a_{h}(\mathbf{e}_{h}^{\mathbf{u}},\mathbf{e}_{h}^{\mathbf{u}})+\big{[}\mathbf{e}^{\mathbf{B}}_{h,t},\mathbf{e}_{h}^{\mathbf{B}}\big{]}_{\operatorname{face}}+\big{[}\mathbf{e}_{h}^{\mathbf{E}}+{\bm{\chi}}_{h}(\mathbf{e}_{h}^{\mathbf{u}},\mathbf{B}_{h}),\mathbf{e}_{h}^{\mathbf{E}}+{\bm{\chi}}_{h}(\mathbf{e}_{h}^{\mathbf{u}},\mathbf{B}_{h})\big{]}_{\operatorname{edge}}$
$\displaystyle=(\mathbf{f},\Pi^{0}_{1}\mathbf{e}_{h}^{\mathbf{u}})-m_{h}(\mathbf{u}_{I,t},\mathbf{e}_{h}^{\mathbf{u}}\big{)}-a_{h}(\mathbf{u}_{I},\mathbf{e}_{h}^{\mathbf{u}})-\widetilde{c}_{h}(\mathbf{u}_{h};\mathbf{u}_{h},\mathbf{e}_{h}^{\mathbf{u}})$
$\displaystyle\phantom{:=}-\big{[}\mathbf{E}_{I\\!I}+{\bm{\chi}}_{h}(\mathbf{u}_{I},\mathbf{B}_{h}),{\bm{\chi}}_{h}(\mathbf{e}_{h}^{\mathbf{u}},\mathbf{B}_{h})\big{]}_{\operatorname{edge}}$
$\displaystyle\phantom{:=}-\big{[}\big{(}\mathbf{E}_{I\\!I}+{\bm{\chi}}_{h}(\mathbf{u}_{I},\mathbf{B}_{h}),\mathbf{e}_{h}^{\mathbf{E}}\big{]}_{\operatorname{edge}}+\big{[}\mathbf{B}_{I\\!I},\operatorname{\mathbf{curl}}\mathbf{e}_{h}^{\mathbf{E}}\big{]}_{\operatorname{face}}=:\textrm{RHS}.$
(50)

We reformulate LHS and RHS in (50) as

$\begin{split}\textrm{LHS}&=\frac{1}{2}\frac{d}{dt}m_{h}(\mathbf{e}_{h}^{\mathbf{u}},\mathbf{e}_{h}^{\mathbf{u}}\big{)}+\frac{1}{2}\frac{d}{dt}\big{[}\mathbf{e}_{h}^{\mathbf{B}},\mathbf{e}_{h}^{\mathbf{B}}\big{]}_{\operatorname{face}}+a_{h}(\mathbf{e}_{h}^{\mathbf{u}},\mathbf{e}_{h}^{\mathbf{u}})+\big{[}\mathbf{e}_{h}^{\mathbf{j}},\mathbf{e}_{h}^{\mathbf{j}}\big{]}_{\operatorname{edge}},\\\\[5.0pt]
\textrm{RHS}&=(\mathbf{f},\Pi^{0}_{1}\mathbf{e}_{h}^{\mathbf{u}})-m_{h}(\mathbf{u}_{I,t},\mathbf{e}_{h}^{\mathbf{u}}\big{)}-a_{h}(\mathbf{u}_{I},\mathbf{e}_{h}^{\mathbf{u}})-\widetilde{c}_{h}(\mathbf{u}_{h};\mathbf{u}_{h},\mathbf{e}_{h}^{\mathbf{u}})-\big{[}\mathbf{j}_{I},\mathbf{e}_{h}^{\mathbf{j}}\big{]}_{\operatorname{edge}}\\\\[5.0pt]
&\phantom{=}+\big{[}\mathbf{B}_{I\\!I},\operatorname{\mathbf{curl}}\mathbf{e}_{h}^{\mathbf{E}}\big{]}_{\operatorname{face}}.\end{split}$
(51)

Next, we observe that

  * •

equation (9a) with
$\mathbf{v}=\mathbf{e}_{h}^{\mathbf{u}}\in\big{[}H^{1}_{0}(\Omega)\big{]}^{3}$
can be rewritten as

$\displaystyle(\mathbf{f},\mathbf{e}_{h}^{\mathbf{u}})=(\mathbf{u}_{t},\mathbf{e}_{h}^{\mathbf{u}})+a(\mathbf{u},\mathbf{e}_{h}^{\mathbf{u}})+\widetilde{c}(\mathbf{u};\mathbf{u},\mathbf{e}_{h}^{\mathbf{u}})+(\mathbf{j},\mathbf{e}_{h}^{\mathbf{u}}\times\mathbf{B}),$
(52)

since $b(\mathbf{e}_{h}^{\mathbf{u}},p)=0$ as
$\operatorname{div}\mathbf{e}_{h}^{\mathbf{u}}=0$ and
$\widetilde{c}(\mathbf{u};\mathbf{u},\mathbf{e}_{h}^{\mathbf{u}})=c(\mathbf{u};\mathbf{u},\mathbf{e}_{h}^{\mathbf{u}})$
since $\operatorname{div}\mathbf{u}=0$ and
$\mathbf{u}\in\big{[}H^{1}_{0}(\Omega)\big{]}$;

  * •

using the definition of $\mathbf{B}_{I\\!I}$ in (32), and (9b) with
$\mathbf{F}=\mathbf{e}_{h}^{\mathbf{E}}\in
H_{0}(\operatorname{\mathbf{curl}},\Omega)$ yields

$\displaystyle\big{[}\mathbf{B}_{I\\!I},\operatorname{\mathbf{curl}}\mathbf{e}_{h}^{\mathbf{E}}\big{]}_{\operatorname{face}}=(\mathbf{B},\operatorname{\mathbf{curl}}\mathbf{e}_{h}^{\mathbf{E}})=(\mathbf{j},\mathbf{e}_{h}^{\mathbf{E}}).$
(53)

We substitute (52) and (53) in RHS defined in (51), add and subtract
$(\mathbf{j},\mathbf{e}_{h}^{\mathbf{u}}\times\mathbf{B}_{h})$, use (38), and
finally obtain

RHS
$\displaystyle:=(\mathbf{f},\Pi^{0}_{1}\mathbf{e}_{h}^{\mathbf{u}})-m_{h}(\mathbf{u}_{I,t},\mathbf{e}_{h}^{\mathbf{u}}\big{)}-a_{h}(\mathbf{u}_{I},\mathbf{e}_{h}^{\mathbf{u}})-\widetilde{c}_{h}(\mathbf{u}_{h};\mathbf{u}_{h},\mathbf{e}_{h}^{\mathbf{u}})-\big{[}\mathbf{j}_{I},\mathbf{e}_{h}^{\mathbf{j}}\big{]}_{\operatorname{edge}}+(\mathbf{j},\mathbf{e}_{h}^{\mathbf{E}})$
$\displaystyle=\Big{[}(\mathbf{u}_{t},\mathbf{e}_{h}^{\mathbf{u}})-m_{h}(\mathbf{u}_{I,t},\mathbf{e}_{h}^{\mathbf{u}}\big{)}\Big{]}+\Big{[}a(\mathbf{u},\mathbf{e}_{h}^{\mathbf{u}})-a_{h}(\mathbf{u}_{I},\mathbf{e}_{h}^{\mathbf{u}})\Big{]}+\Big{[}\widetilde{c}(\mathbf{u};\mathbf{u},\mathbf{e}_{h}^{\mathbf{u}})-\widetilde{c}_{h}(\mathbf{u}_{h};\mathbf{u}_{h},\mathbf{e}_{h}^{\mathbf{u}})\Big{]}$
$\displaystyle\phantom{:=}+(\mathbf{j},\mathbf{e}_{h}^{\mathbf{u}}\times\mathbf{B})-\big{[}\mathbf{j}_{I},\mathbf{e}_{h}^{\mathbf{j}}\big{]}_{\operatorname{edge}}+(\mathbf{j},\mathbf{e}_{h}^{\mathbf{E}})+(\mathbf{f},\Pi^{0}_{1}\mathbf{e}_{h}^{\mathbf{u}})-(\mathbf{f},\mathbf{e}_{h}^{\mathbf{u}})$
$\displaystyle=\Big{[}(\mathbf{u}_{t},\mathbf{e}_{h}^{\mathbf{u}})-m_{h}(\mathbf{u}_{I,t},\mathbf{e}_{h}^{\mathbf{u}}\big{)}\Big{]}+\Big{[}\widetilde{c}(\mathbf{u};\mathbf{u},\mathbf{e}_{h}^{\mathbf{u}})-\widetilde{c}_{h}(\mathbf{u}_{h};\mathbf{u}_{h},\mathbf{e}_{h}^{\mathbf{u}})\Big{]}+(\mathbf{j},\mathbf{e}_{h}^{\mathbf{E}}+\mathbf{e}_{h}^{\mathbf{u}}\times\mathbf{B})$
$\displaystyle\phantom{:=}-\big{[}\mathbf{j}_{I},\mathbf{e}_{h}^{\mathbf{j}}\big{]}_{\operatorname{edge}}+\Big{[}a(\mathbf{u},\mathbf{e}_{h}^{\mathbf{u}})-a_{h}(\mathbf{u}_{I},\mathbf{e}_{h}^{\mathbf{u}})\Big{]}+\Big{[}(\mathbf{f},\Pi^{0}_{1}\mathbf{e}_{h}^{\mathbf{u}})-(\mathbf{f},\mathbf{e}_{h}^{\mathbf{u}})\Big{]}$
$\displaystyle=\Big{[}(\mathbf{u}_{t},\mathbf{e}_{h}^{\mathbf{u}})-m_{h}(\mathbf{u}_{I,t},\mathbf{e}_{h}^{\mathbf{u}}\big{)}\Big{]}+\Big{[}\widetilde{c}(\mathbf{u};\mathbf{u},\mathbf{e}_{h}^{\mathbf{u}})-\widetilde{c}_{h}(\mathbf{u}_{h};\mathbf{u}_{h},\mathbf{e}_{h}^{\mathbf{u}})\Big{]}+\Big{[}\big{(}\mathbf{j},\mathbf{e}_{h}^{\mathbf{u}}\times(\mathbf{B}-\mathbf{B}_{h})\big{)}\Big{]}$
$\displaystyle\phantom{=}+\Big{[}\big{(}\mathbf{j},\mathbf{e}_{h}^{\mathbf{E}}+\mathbf{e}_{h}^{\mathbf{u}}\times\mathbf{B}_{h}\big{)}-\big{[}\mathbf{j}_{I},\mathbf{e}_{h}^{\mathbf{j}}\big{]}_{\operatorname{edge}}\Big{]}+\Big{[}a(\mathbf{u},\mathbf{e}_{h}^{\mathbf{u}})-a_{h}(\mathbf{u}_{I},\mathbf{e}_{h}^{\mathbf{u}})\Big{]}+\Big{[}(\mathbf{f},\Pi^{0}_{1}\mathbf{e}_{h}^{\mathbf{u}})-(\mathbf{f},\mathbf{e}_{h}^{\mathbf{u}})\Big{]}$
$\displaystyle=:\textrm{T}_{1}+\textrm{T}_{2}+\textrm{T}_{3}+\textrm{T}_{4}+\textrm{T}_{5}+\textrm{T}_{6}.$
(54)

We estimate the six terms $\textrm{T}_{j}$, $j=1,\dots,6$, separately, and we
eventually apply Grönwall’s inequality to derive the assertion of the theorem.

##### Estimate of $\textrm{T}_{1}$

We split $\textrm{T}_{1}$ into a summation on the elemental contributions;
note that $\mathbf{u}_{I,t}=\mathbf{u}_{t,I}$; use the consistency property
(14) with the piecewise linear discontinuous polynomial approximation
$\mathbf{u}_{t,\pi}$ of $\mathbf{u}_{t}$; use the stability property (15), the
Cauchy-Schwarz inequality, and the Young inequality, and obtain:

$\displaystyle\textrm{T}_{1}$
$\displaystyle=\sum_{K\in\Omega_{h}}\left|\big{(}\mathbf{u}_{t},\mathbf{e}_{h}^{\mathbf{u}}\big{)}_{K}-m_{h}^{K}(\mathbf{u}_{I,t},\mathbf{e}_{h}^{\mathbf{u}})\right|=\sum_{K\in\Omega_{h}}\left|\big{(}\mathbf{u}_{t}-\mathbf{u}_{t,\pi},\mathbf{e}_{h}^{\mathbf{u}}\big{)}_{K}-m_{h}^{K}(\mathbf{u}_{t,I}-\mathbf{u}_{t,\pi},\mathbf{e}_{h}^{\mathbf{u}})\right|$
$\displaystyle\leq\sum_{K\in\Omega_{h}}\big{(}\|\mathbf{u}_{t}-\mathbf{u}_{t,\pi}\|_{K}+\mu^{*}\|\mathbf{u}_{t,I}-\mathbf{u}_{t,\pi}\|_{K}\big{)}\,\|\mathbf{e}_{h}^{\mathbf{u}}\|_{K}\leq\big{(}\|\mathbf{u}_{t}-\mathbf{u}_{t,\pi}\|+\mu^{*}\|\mathbf{u}_{t,I}-\mathbf{u}_{t,\pi}\|\big{)}\,\|\mathbf{e}_{h}^{\mathbf{u}}\|$
$\displaystyle\leq\big{(}(1+\mu^{*})\|\mathbf{u}_{t}-\mathbf{u}_{t,\pi}\|+\mu^{*}\|\mathbf{u}_{t}-\mathbf{u}_{t,I}\|\big{)}\,\|\mathbf{e}_{h}^{\mathbf{u}}\|$
$\displaystyle\leq\frac{1}{2\varepsilon}\Big{(}\big{(}1+\mu^{*}\big{)}^{2}\|\mathbf{u}_{t}-\mathbf{u}_{t,\pi}\|^{2}+(\mu^{*})^{2}\|\mathbf{u}_{t}-\mathbf{u}_{t,I}\|^{2}\Big{)}+\varepsilon\|\mathbf{e}_{h}^{\mathbf{u}}\|^{2}$
$\displaystyle\leq
h^{2}\,C_{\textrm{T}_{1}}(\varepsilon)+\varepsilon\|\mathbf{e}_{h}^{\mathbf{u}}\|^{2},$

with
$C_{\textrm{T}_{1}}(\varepsilon):=\big{(}C_{\pi}^{2}(1+\mu^{*})^{2}+C_{I}^{2}(\mu^{*})^{2}\big{)}|\mathbf{u}_{t}|_{1}^{2}/{2\varepsilon}$.

##### Estimate of $\textrm{T}_{2}$

We add and subtract
$\widetilde{c}_{h}(\mathbf{u};\mathbf{u},\mathbf{e}_{h}^{\mathbf{u}})$, use
the triangle inequality, and write

$\displaystyle\textrm{T}_{2}$
$\displaystyle=\left|\widetilde{c}(\mathbf{u};\mathbf{u},\mathbf{e}_{h}^{\mathbf{u}})-\widetilde{c}_{h}(\mathbf{u}_{h};\mathbf{u}_{h},\mathbf{e}_{h}^{\mathbf{u}})\right|\leq\left|\widetilde{c}(\mathbf{u};\mathbf{u},\mathbf{e}_{h}^{\mathbf{u}})-\widetilde{c}_{h}(\mathbf{u};\mathbf{u},\mathbf{e}_{h}^{\mathbf{u}})\right|+\left|\widetilde{c}_{h}(\mathbf{u};\mathbf{u},\mathbf{e}_{h}^{\mathbf{u}})-\widetilde{c}_{h}(\mathbf{u}_{h};\mathbf{u}_{h},\mathbf{e}_{h}^{\mathbf{u}})\right|$
$\displaystyle=:\textrm{T}_{2,1}+\textrm{T}_{2,2}.$

We control $\textrm{T}_{2,1}$ by using Lemma 4.5 and the Young inequality:

$\displaystyle\textrm{T}_{2,1}\leq
C_{C}h\|\mathbf{u}\|_{2}^{2}\,|\mathbf{e}_{h}^{\mathbf{u}}|_{1}\leq\frac{C_{C}^{2}}{4\varepsilon}h^{2}\|\mathbf{u}\|_{2}^{4}+\varepsilon|\mathbf{e}_{h}^{\mathbf{u}}|_{1}^{2},$

where the constant $C_{C}$ is independent of $h$. Similarly, we control
$\textrm{T}_{2,2}$ using Lemma 4.6:

$\displaystyle\textrm{T}_{2,2}$
$\displaystyle\leq\varepsilon|\mathbf{e}_{h}^{\mathbf{u}}|_{1}^{2}+\Big{(}1+\frac{1}{\varepsilon}\Big{)}C_{1}R_{1}(\mathbf{u})\|\mathbf{e}_{h}^{\mathbf{u}}\|^{2}+\Big{(}1+\frac{1}{\varepsilon}\Big{)}C_{2}R_{2}(\mathbf{u})h^{2},$

where the constants $C_{1}$ and $C_{2}$ are independent of $h$. We end up with

$\textrm{T}_{2}\leq
h^{2}\,C_{\textrm{T}_{2}}(\varepsilon)+\Big{(}1+\frac{1}{\varepsilon}\Big{)}C_{1}R_{1}(\mathbf{u})\|\mathbf{e}_{h}^{\mathbf{u}}\|^{2}+2\varepsilon|\mathbf{e}_{h}^{\mathbf{u}}|_{1}^{2},$

where

$C_{\textrm{T}_{2}}(\varepsilon):=\frac{C_{C}^{2}}{4\varepsilon}\|\mathbf{u}\|_{2}^{4}+\Big{(}1+\frac{1}{\varepsilon}\Big{)}C_{2}R_{2}(\mathbf{u}).$

##### Estimate of $\textrm{T}_{3}$

By adding and subtracting $\mathbf{B}_{I\\!I}$ defined in (32), we obtain

$\displaystyle\textrm{T}_{3}$
$\displaystyle=\big{(}\mathbf{j},\mathbf{e}_{h}^{\mathbf{u}}\times(\mathbf{B}-\mathbf{B}_{h})\big{)}=\big{(}\mathbf{j},\mathbf{e}_{h}^{\mathbf{u}}\times(\mathbf{B}-\mathbf{B}_{I\\!I})\big{)}+\big{(}\mathbf{j},\mathbf{e}_{h}^{\mathbf{u}}\times(\mathbf{B}_{I\\!I}-\mathbf{B}_{h})\big{)}=:\textrm{T}_{3,1}+\textrm{T}_{3,2}.$

We estimate the two terms $\textrm{T}_{3,1}$ and $\textrm{T}_{3,2}$
separately. Notably, we write

$\displaystyle\left|\textrm{T}_{3,1}\right|$
$\displaystyle\leq\|\mathbf{j}\|_{L^{\infty}}\,\|\mathbf{e}_{h}^{\mathbf{u}}\|\,\|\mathbf{B}-\mathbf{B}_{I\\!I}\|\leq\frac{\|\mathbf{j}\|_{L^{\infty}}}{2}\,\bigg{(}\|\mathbf{e}_{h}^{\mathbf{u}}\|^{2}+\|\mathbf{B}-\mathbf{B}_{I\\!I}\|^{2}\bigg{)}$
$\displaystyle\leq\frac{\|\mathbf{j}\|_{L^{\infty}}}{2}\,\bigg{(}\|\mathbf{e}_{h}^{\mathbf{u}}\|^{2}+C_{I\\!I}^{2}h^{2}|\mathbf{B}|_{1}^{2}\bigg{)}$

and

$\displaystyle\left|\textrm{T}_{3,2}\right|\leq\|\mathbf{j}\|_{L^{\infty}}\,\|\mathbf{e}_{h}^{\mathbf{u}}\|\,\|\mathbf{B}_{I\\!I}-\mathbf{B}_{h}\|\leq\|\mathbf{j}\|_{L^{\infty}}\,\|\mathbf{e}_{h}^{\mathbf{u}}\|\,\|\mathbf{e}_{h}^{\mathbf{B}}\|\leq\frac{1}{2}\|\mathbf{j}\|_{L^{\infty}}\,\big{(}\|\mathbf{e}_{h}^{\mathbf{u}}\|^{2}+\|\mathbf{e}_{h}^{\mathbf{B}}\|^{2}\big{)}.$

Collecting the two above estimates, we find that

$\displaystyle\textrm{T}_{3}\leq
h^{2}\,C_{\textrm{T}_{3}}+\|\mathbf{j}\|_{L^{\infty}}\,\|\mathbf{e}_{h}^{\mathbf{u}}\|^{2}+\frac{1}{2}\|\mathbf{j}\|_{L^{\infty}}\|\mathbf{e}_{h}^{\mathbf{B}}\|^{2},$

where

$C_{\textrm{T}_{3}}:=\frac{C_{I\\!I}^{2}}{2}\,\|\mathbf{j}\|_{L^{\infty}}|\mathbf{B}|_{1}^{2}.$

##### Estimate of $\textrm{T}_{4}$

Recalling that
$\mathbf{e}_{h}^{\mathbf{j}}=\mathbf{e}_{h}^{\mathbf{E}}+\mathbf{e}_{h}^{\mathbf{u}}\times\mathbf{B}_{h}$,
we use (by now) standard manipulations based on the consistency and stability
properties (26) and (27), and get

$\displaystyle\textrm{T}_{4}$
$\displaystyle=\big{(}\mathbf{j},\mathbf{e}_{h}^{\mathbf{j}}\big{)}-\big{[}\mathbf{j}_{I},\mathbf{e}_{h}^{\mathbf{j}}\big{]}_{\operatorname{edge}}=\sum_{K\in\Omega_{h}}\Big{(}\big{(}\mathbf{j},\mathbf{e}_{h}^{\mathbf{j}}\big{)}_{K}-\big{[}\mathbf{j}_{I},\mathbf{e}_{h}^{\mathbf{j}}\big{]}_{\operatorname{edge},K}\Big{)}$
$\displaystyle=\sum_{K\in\Omega_{h}}\Big{(}\big{(}\mathbf{j}-\mathbf{j}_{\pi},\mathbf{e}_{h}^{\mathbf{j}}\big{)}_{K}-\big{[}\mathbf{j}_{I}-\mathbf{j}_{\pi},\mathbf{e}_{h}^{\mathbf{j}}\big{]}_{\operatorname{edge},K}\Big{)}\leq\sum_{K\in\Omega_{h}}\Big{(}\|\mathbf{j}-\mathbf{j}_{\pi}\|_{0,K}+\eta^{*}\|\mathbf{j}_{I}-\mathbf{j}_{\pi}\|_{0,K}\Big{)}\|\mathbf{e}_{h}^{\mathbf{j}}\|_{0,K}$
$\displaystyle\leq\sum_{K\in\Omega_{h}}\Big{(}\big{(}1+\eta^{*}\big{)}\|\mathbf{j}-\mathbf{j}_{\pi}\|_{0,K}+\eta^{*}\|\mathbf{j}-\mathbf{j}_{I}\|_{0,K}\Big{)}\|\mathbf{e}_{h}^{\mathbf{j}}\|_{0,K}\leq\Big{(}\big{(}1+\eta^{*}\big{)}\|\mathbf{j}-\mathbf{j}_{\pi}\|+\eta^{*}\|\mathbf{j}-\mathbf{j}_{I}\|\Big{)}\|\mathbf{e}_{h}^{\mathbf{j}}\|$
$\displaystyle\leq\frac{1}{2\varepsilon}\Big{(}\big{(}1+\eta^{*}\big{)}^{2}\|\mathbf{j}-\mathbf{j}_{\pi}\|^{2}+(\eta^{*})^{2}\|\mathbf{j}-\mathbf{j}_{I}\|^{2}\Big{)}+\varepsilon\|\mathbf{e}_{h}^{\mathbf{j}}\|^{2}$
$\displaystyle\leq\frac{\big{(}1+\eta^{*}\big{)}^{2}C_{\operatorname{appr}}^{2}}{2\varepsilon}h^{2}|\mathbf{j}|_{1}^{2}+\frac{(\eta^{*})^{2}}{2\varepsilon}\|\mathbf{j}-\mathbf{j}_{I}\|^{2}+\varepsilon\|\mathbf{e}_{h}^{\mathbf{j}}\|^{2}.$
(55)

Above, $\mathbf{j}_{\pi}$ stands for the piecewise constant best approximant
of $\mathbf{j}$.

Then, we transform $\|\mathbf{j}-\mathbf{j}_{I}\|$ by using the triangle
inequality and adding and subtracting $\bm{\Pi}^{0,K}_{h}\mathbf{u}_{I}$ and
$\bm{\Pi}^{0,K}_{h}\mathbf{B}$:

$\displaystyle\|\mathbf{j}-\mathbf{j}_{I}\|$
$\displaystyle=\|\mathbf{E}+\mathbf{u}\times\mathbf{B}-\big{(}\mathbf{E}_{I\\!I}+\bm{\Pi}^{0,K}_{h}\mathbf{u}_{I}\times\bm{\Pi}^{0,K}_{h}\mathbf{B}_{h}\big{)}\|$
$\displaystyle\leq\|\mathbf{E}-\mathbf{E}_{I\\!I}\|+\|\mathbf{u}\times\mathbf{B}-\bm{\Pi}^{0,K}_{h}\mathbf{u}_{I}\times\bm{\Pi}^{0,K}_{h}\mathbf{B}_{h}\|$
$\displaystyle\leq
C_{I}h|\mathbf{E}|+\|\big{(}\mathbf{u}-\bm{\Pi}^{0,K}_{h}\mathbf{u}_{I}\big{)}\times\mathbf{B}\|+\|\bm{\Pi}^{0,K}_{h}\mathbf{u}_{I}\times\big{(}\mathbf{B}-\bm{\Pi}^{0,K}_{h}\mathbf{B}\big{)}\|$
$\displaystyle\hskip
35.85031pt+\|\bm{\Pi}^{0,K}_{h}\mathbf{u}_{I}\times\big{(}\bm{\Pi}^{0,K}_{h}\mathbf{B}-\bm{\Pi}^{0,K}_{h}\mathbf{B}_{h}\big{)}\|$
$\displaystyle=:C_{I}h|\mathbf{E}|+\textrm{T}_{4,1}+\textrm{T}_{4,2}+\textrm{T}_{4,3}.$
(56)

Next, we estimate separately the three terms $\textrm{T}_{4,1}$,
$\textrm{T}_{4,2}$, and $\textrm{T}_{4,3}$. We control $\textrm{T}_{4,1}$ by
adding and substracting $\bm{\Pi}^{0,K}_{h}\mathbf{u}$, and using the triangle
inequality:

$\displaystyle\textrm{T}_{4,1}$
$\displaystyle=\|\big{(}\mathbf{u}-\bm{\Pi}^{0,K}_{h}\mathbf{u}_{I}\big{)}\times\mathbf{B}\|\leq\|\mathbf{B}\|_{L^{\infty}}\|\mathbf{u}-\bm{\Pi}^{0,K}_{h}\mathbf{u}_{I}\|$
$\displaystyle\leq\|\mathbf{B}\|_{L^{\infty}}\Big{(}\|\mathbf{u}-\bm{\Pi}^{0,K}_{h}\mathbf{u}\|+\|\bm{\Pi}^{0,K}_{h}\big{(}\mathbf{u}-\mathbf{u}_{I}\big{)}\|\Big{)}\leq\|\mathbf{B}\|_{L^{\infty}}\,h(C_{\operatorname{appr}}+C_{I})\|\mathbf{u}\|_{1}.$

We control $\textrm{T}_{4,2}$ by using Lemma 4.3:

$\displaystyle\textrm{T}_{4,2}$
$\displaystyle=\|\bm{\Pi}^{0,K}_{h}\mathbf{u}_{I}\times\big{(}\mathbf{B}-\bm{\Pi}^{0,K}_{h}\mathbf{B}\big{)}\|\leq\|\bm{\Pi}^{0,K}_{h}\mathbf{u}_{I}\|_{L^{\infty}}\,\|\mathbf{B}-\bm{\Pi}^{0,K}_{h}\mathbf{B}\|$
$\displaystyle\leq(C_{\operatorname{inv}}C_{I}h_{K}^{\frac{1}{2}}+C_{\operatorname{inv}}C_{D}C_{S})\|\mathbf{u}\|_{2}\,hC_{\operatorname{appr}}|\mathbf{B}|_{1}.$

We control $\textrm{T}_{4,3}$ by adding and subtracting $\mathbf{B}_{I\\!I}$
defined in (32) and using again Lemma 4.3:

$\displaystyle\textrm{T}_{4,3}$
$\displaystyle=\|\bm{\Pi}^{0,K}_{h}\mathbf{u}_{I}\times\big{(}\bm{\Pi}^{0,K}_{h}(\mathbf{B}-\mathbf{B}_{h})\|\leq\|\bm{\Pi}^{0,K}_{h}\mathbf{u}_{I}\|_{L^{\infty}}\,\|\mathbf{B}-\mathbf{B}_{h}\|$
$\displaystyle\leq\|\bm{\Pi}^{0,K}_{h}\mathbf{u}_{I}\|_{L^{\infty}}\,\Big{(}\|\mathbf{B}-\mathbf{B}_{I\\!I}\|+\|\mathbf{B}_{I}-\mathbf{B}_{h}\|\Big{)}$
$\displaystyle\leq(C_{\operatorname{inv}}C_{I}h_{K}^{\frac{1}{2}}+C_{\operatorname{inv}}C_{D}C_{S})\|\mathbf{u}\|_{2}(hC_{I}|\mathbf{B}|_{1}+\|\mathbf{e}_{h}^{\mathbf{B}}\|).$

We collect the three above estimates in (6.1):

$\displaystyle\|\mathbf{j}-\mathbf{j}_{I}\|$ $\displaystyle\leq
h\Big{(}C_{I}|\mathbf{E}|_{1}+(C_{\operatorname{appr}}+C_{I})\|\mathbf{B}\|_{L^{\infty}}\|\mathbf{u}\|_{1}+(C_{\operatorname{inv}}C_{I}h_{K}^{\frac{1}{2}}+C_{\operatorname{inv}}C_{D}C_{S})(C_{\operatorname{appr}}+C_{I})\|\mathbf{u}\|_{2}|\mathbf{B}|_{1}\Big{)}$
$\displaystyle\qquad+(C_{\operatorname{inv}}C_{I}h_{K}^{\frac{1}{2}}+C_{\operatorname{inv}}C_{D}C_{S})\|\mathbf{u}\|_{2}\|\mathbf{e}_{h}^{\mathbf{B}}\|,$

and this inequality in (55), and obtain the following upper bound for
$\textrm{T}_{4}$:

$\textrm{T}_{4}\leq
h^{2}C_{\textrm{T}_{4}}(\varepsilon)+\frac{(\eta^{*})^{2}}{\varepsilon}(C_{\operatorname{inv}}C_{I}h_{K}^{\frac{1}{2}}+C_{\operatorname{inv}}C_{D}C_{S})^{2}\|\mathbf{u}\|_{2}^{2}\|\mathbf{e}_{h}^{\mathbf{B}}\|^{2}+\varepsilon\|\mathbf{e}_{h}^{\mathbf{j}}\|^{2},$

where $C_{\textrm{T}_{4}}(\varepsilon)\approx\varepsilon^{-1}$ collects all
the constants above and depends also on $|\mathbf{B}|$, $|\mathbf{j}|_{1}$,
$|\mathbf{E}|_{1}$, $\|\mathbf{B}\|_{L^{\infty}}$, and $\|\mathbf{u}\|_{2}$.

##### Estimate of $\textrm{T}_{5}$

We proceed with (by now) standard techniques based on the consistency and
stability properties (11) and (12): given $\mathbf{u}_{\pi}$ the piecewise
linear best approximant of $\mathbf{u}$,

$\displaystyle\textrm{T}_{5}$
$\displaystyle=\sum_{K\in\Omega_{h}}\left[a^{K}(\mathbf{u},\mathbf{e}_{h}^{\mathbf{u}})-a_{h}^{K}(\mathbf{u}_{I},\mathbf{e}_{h}^{\mathbf{u}})\right]=\sum_{K\in\Omega_{h}}\left[a^{K}(\mathbf{u}-\mathbf{u}_{\pi},\mathbf{e}_{h}^{\mathbf{u}})-a_{h}^{K}(\mathbf{u}_{I}-\mathbf{u}_{\pi},\mathbf{e}_{h}^{\mathbf{u}})\right]$
$\displaystyle\leq\sum_{K\in\Omega_{h}}\big{[}|\mathbf{u}-\mathbf{u}_{\pi}|_{1,K}+\sigma^{*}\left(|\mathbf{u}-\mathbf{u}_{I}|_{1,K}+|\mathbf{u}-\mathbf{u}_{\pi}|_{1,K}\right)\big{]}|\mathbf{e}_{h}^{\mathbf{u}}|_{1,K}$
$\displaystyle\leq
h\left((1+\sigma^{*})C_{\operatorname{appr}}+\sigma^{*}C_{I})\right)|\mathbf{u}|_{2}|\mathbf{e}_{h}^{\mathbf{u}}|_{1}\leq
h^{2}\
C_{\textrm{T}_{5}}(\varepsilon)+\varepsilon|\mathbf{e}_{h}^{\mathbf{u}}|_{1}^{2},$

where

$C_{\textrm{T}_{5}}(\varepsilon):=\frac{1}{4\varepsilon}\left((1+\sigma^{*})C_{\operatorname{appr}}+\sigma^{*}C_{I})\right)^{2}|\mathbf{u}|_{1}^{2}.$

##### Estimate of $\textrm{T}_{6}$

Standard manipulations based on the Poincaré and the Young inequalities yield

$\begin{split}\textrm{T}_{6}&=(\mathbf{f},\Pi^{0}_{1}\mathbf{e}_{h}^{\mathbf{u}})-(\mathbf{f},\mathbf{e}_{h}^{\mathbf{u}})=(\mathbf{f}-\Pi^{0}_{1}\mathbf{f},\Pi^{0}_{1}\mathbf{e}_{h}^{\mathbf{u}}-\mathbf{e}_{h}^{\mathbf{u}})\\\
&\leq\|\mathbf{f}-\Pi^{0}_{1}\mathbf{f}\|\|\mathbf{e}_{h}^{\mathbf{u}}-\Pi^{0}_{1}\mathbf{e}_{h}^{\mathbf{u}}\|\leq
C_{\operatorname{appr}}^{2}h^{2}\|\mathbf{f}\|_{1}|\mathbf{e}_{h}^{\mathbf{u}}|_{1}\leq
h^{2}C_{\textrm{T}_{6}}+\varepsilon|\mathbf{e}_{h}^{\mathbf{u}}|_{1}^{2},\end{split}$

where

$C_{\textrm{T}_{6}}:=\frac{1}{4}h^{2}C_{\operatorname{appr}}^{2}\|\mathbf{f}\|_{1}^{2}.$

##### Use of the Grönwall’s inequality and final step

Combining the above bounds on the term $\textrm{T}_{i}$, $i=1,\dots,6$, we get

$\displaystyle\sum_{i=1}^{6}\textrm{T}_{i}$ $\displaystyle\leq
S_{1}(\varepsilon)h^{2}+S_{2}(\varepsilon)\|\mathbf{e}_{h}^{\mathbf{u}}\|^{2}+5\varepsilon|\mathbf{e}_{h}^{\mathbf{u}}|_{1}^{2}+S_{3}(\varepsilon)\|\mathbf{e}_{h}^{\mathbf{B}}\|^{2}+\varepsilon\|\mathbf{e}_{h}^{\mathbf{j}}\|^{2},$

where

$\displaystyle S_{1}(\varepsilon)$
$\displaystyle=\sum_{j=1}^{6}C_{\textrm{T}_{j}}(\varepsilon),\qquad
S_{3}(\varepsilon)=\frac{1}{2}\|\mathbf{j}\|_{L^{\infty}}+\frac{(\eta^{*})^{2}}{\varepsilon}(C_{\operatorname{inv}}C_{I}h_{K}^{\frac{1}{2}}+C_{\operatorname{inv}}C_{D}C_{S})^{2}\|\mathbf{u}\|_{2}^{2},$
$\displaystyle S_{2}(\varepsilon)$
$\displaystyle=\varepsilon+(1+1/\varepsilon)C_{1}R_{1}(\mathbf{u})+\|\mathbf{j}\|_{L^{\infty}}.$

The terms $S_{j}(\varepsilon)$, $j=1,2,3$, depend on the time instant $t$,
since they depend on $\mathbf{u}$, $\mathbf{B}$, and $\mathbf{j}$. For the
presentation’s sake, we do not express this dependence explicitly in our
notation.

We substitute the above bound in (54) and the resulting inequality in (50).
Next, we pick

$\varepsilon=\widetilde{\varepsilon}=\min(\alpha_{*}^{-1},\eta_{*}^{-1})/10$

and, at every time instant $t\in(0,T]$, we obtain

$\begin{split}\textrm{LHS}{}&\leq
2\left[S_{1}(\widetilde{\varepsilon})h^{2}+S_{2}(\widetilde{\varepsilon})\|\mathbf{e}_{h}^{\mathbf{u}}\|^{2}+S_{3}(\widetilde{\varepsilon})\,\|\mathbf{e}_{h}^{\mathbf{B}}\|^{2}\right]\\\\[5.0pt]
&\leq