# Virtual elements for Maxwell’s equations

L. Beirão da Veiga  , F. Dassi  , G. Manzini  , and L. Mascotto  Dip. di
Matematica e Applicazioni, Università degli Studi di Milano-Bicocca, Italy
<EMAIL_ADDRESS> Dip. di Matematica e Applicazioni, Università
degli Studi di Milano-Bicocca, Italy<EMAIL_ADDRESS> Istituto di
Matematica Applicata e Tecnologie Informatiche, CNR, Pavia, Italy
<EMAIL_ADDRESS> Fakultät für Mathematik, Universität Wien,
Austria<EMAIL_ADDRESS>
###### Abstract

We present a low order virtual element discretization for time dependent
Maxwell’s equations, which allow for the use of general polyhedral meshes.
Both the semi- and fully-discrete schemes are considered. We derive optimal a
priori estimates and validate them on a set of numerical experiments. As pivot
results, we discuss some novel inequalities for de Rahm sequences of nodal,
edge, and face virtual element spaces.
AMS subject classification: 65N12; 65N15

###### keywords:

polyhedral meshes; virtual element method; Maxwell’s equations

##  1 Introduction

The Virtual Element Method (VEM) was introduced in [8] as a generalization of
the Finite Element Method (FEM) that allows for the use of general polygonal
and polyhedral meshes. Since its introduction, the VEM has shared a wide
success in the numerical analysis and engineering communities. After the
introduction of $H^{1}$ conforming spaces in [8, 2, 14], also $H({\rm div})$
and $H({\rm curl})$ conforming spaces in both two and three space dimensions
were proposed. Mixed finite elements for the diffusion problem in mixed form
in 2D were introduced in [19, 13], while in [12, 9, 11, 10] various families
of discrete exact VEM complexes of $H^{1}-H({\rm div})-H({\rm curl})-L^{2}$
type were introduced in 2D and 3D. In the above contributions, all such
families of spaces are applied to the Kikuchi formulation of the magnetostatic
equations, used as a simple model problem to showcase the proposed discrete
construction. A recent application for permanent magnet simulations can be
found in [23].

On the other hand, finite elements have been widely used for numerical
modelling of Maxwell’s equations, a very short representative list being [29,
16, 36, 37, 39, 34, 45, 21, 6]. Important applications involve, for instance,
the analysis and design of microwave devices [22], cavity resonators [32, 43,
40], coaxial cables and waveguides [44], antennas and high-power amplifiers
[42, 41, 28, 30], electromagnetic scattering [31, 35].

Due to the complex geometries that are often faced in many applicative areas
of electromagnetism, the additional flexibility of general polytopal grids is
an important asset, not only in generating an efficient mesh to partition the
domain of interest, but also in handling/gluing/adapting existing meshes.
Among the other polytopal technologies, in the realm of electromagnetism it is
possible to find (in a nonexhaustive list) polygonal finite elements [26],
mimetic finite differences [33], hybrid high-order methods [20], and discrete
exact sequences [25].

The aim of the present paper is to use the discrete spaces introduced in [11]
to develop a virtual element discretization of the full time-dependent
Maxwell’s equations. In order to ease the reader’s understanding, we restrict
the presentation and analysis to the lowest order case; the generalization of
the scheme and the analysis to the general order case, see, e.g., [10], would
follow the same steps.

#### Structure of the paper.

After introducing several Sobolev spaces at the end of this introduction, we
present the model problem in Section 2. The virtual element schemes for the
semi- and fully discrete Maxwell’s equations are detailed in Section 3; here,
we also address the approximation properties in virtual element and polynomial
spaces, as well as the design of suitable stabilization terms. We develop
convergence estimates for the semi-discrete and the fully discrete cases in
the spirit of [45], the latter restricted to the backward Euler case, in
Sections 4 and 5. The error estimates show the optimal behaviour of the
proposed method. In order to investigate the practical performance of the
scheme, we develop a set of academic numerical tests in Section 6. Eventually,
we state some conclusions in Section 7.

#### Notation and functional spaces.

We employ the standard definitions and notation for Hilbert and Sobolev spaces
[1]. Given $s\in\mathbb{R}$ and a Lipschitz domain $D$, we denote the Hilbert
space of order $s$ by $H^{s}(D)$. We endow $H^{s}(D)$ with the standard inner
product, norm, and seminorms, which we indicate as $(\cdot,\cdot)_{s,D}$,
$|\cdot|_{s,D}$, and $\|\cdot\|_{s,D}$. The special case $s=0$ consists of the
Lebesgue space $L^{2}(D)$ of real-valued, square integrable functions defined
on $D$. We define Sobolev spaces of noninteger order by interpolation and
Sobolev spaces of negative order by duality. We analogously consider Sobolev
spaces $H^{s}(\partial D)$ on the boundary $\partial D$ of $D$.

We recall the definition of some differential operators that we shall use in
the paper. Let $\partial_{x}$, $\partial_{y}$, and $\partial_{z}$ denote the
partial derivative along $x$, $y$, and $z$. Given a two-dimensional vector-
valued field
$\mathbf{v}=(v_{1},v_{2}):F\subseteq\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$
and a scalar field $v:F\subseteq\mathbb{R}^{2}\rightarrow\mathbb{R}$, we
consider

$\displaystyle\operatorname{div_{F}}\mathbf{v}:=\partial_{x}v_{1}+\partial_{y}v_{2},\quad\quad\operatorname{rot_{F}}\mathbf{v}:=\partial_{y}v_{1}-\partial_{x}v_{2},\quad\quad\operatorname{\mathbf{curl}_{F}}v:=\big{(}\partial_{y}v,-\partial_{x}v\big{)}^{T}.$

In turn, given a three-dimensional vector-valued field
$\mathbf{v}=(v_{1},v_{2},v_{3}):K\subseteq\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$,
we consider

$\displaystyle\operatorname{div}\mathbf{v}:=\partial_{x}v_{1}+\partial_{y}v_{2}+\partial_{z}v_{3},\quad\quad\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}.$

For Lipschitz domains $D\subset\mathbb{R}^{3}$, we introduce the Sobolev
$\operatorname{\mathbf{curl}}$ and $\operatorname{div}$ spaces of order $s>0$

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

If $s=0$, we write $H(\operatorname{\mathbf{curl}},D)$ and
$H(\operatorname{div},D)$. We denote the unit vector that is orthogonal to the
boundary $\partial D$ and pointing out of $D$ by $\mathbf{n}_{D}$.
Furthermore, we recall the existence of the two trace operators
$\text{tr}_{\textbf{curl}}:H(\operatorname{\mathbf{curl}},D)\rightarrow[H^{-\frac{1}{2}}(\partial
D)]^{2}$ and $\text{tr}_{\text{div}}:H(\operatorname{div},D)\rightarrow
H^{-\frac{1}{2}}(\partial D)$ such that
$\text{tr}_{\textbf{curl}}(\mathbf{v})=\mathbf{v}\times\mathbf{n}_{D}$ and
$\text{tr}_{\text{div}}(\mathbf{v})={\bm{\psi}}\cdot\mathbf{n}_{D}$ for all
$\mathbf{v}\in H(\operatorname{\mathbf{curl}},D)$ and ${\bm{\psi}}\in
H(\operatorname{div},D)$, respectively; see, e.g., [36, Section 3.5].
According to the standard notation, $L^{\infty}(D)$ is the Sobolev space of
functions that are bounded almost everywhere and $W^{1,\infty}(D)$ the Sobolev
space of functions in $L^{\infty}(D)$ whose first weak derivatives are also in
$L^{\infty}(D)$. We shall also consider $\operatorname{\mathbf{curl}}$\- and
$\operatorname{div}$-spaces with zero boundary conditions such as

$\displaystyle H_{0}(\operatorname{\mathbf{curl}},D)$
$\displaystyle:=\left\\{\mathbf{v}\in
H(\operatorname{\mathbf{curl}},D)\mid\text{tr}_{\textbf{curl}}(\mathbf{v})=0\text{~{}on~{}}\partial
D\right\\},$ $\displaystyle H_{0}(\operatorname{div},D)$
$\displaystyle:=\left\\{\mathbf{v}\in
H(\operatorname{div},D)\mid\text{tr}_{\text{div}}(\mathbf{v})=0\text{~{}on~{}}\partial
D\right\\}.$

Let $X$ denote a scalar or vector Sobolev space of any order over the domain
$D\subset\mathbb{R}^{3}$; $(a,b)$ an open, connected subset of $\mathbb{R}$,
and $p$ a real number in the interval $[1,\infty]$. The Bochner space [27]
$L^{p}((a,b),X)$ is the vector space of functions $\mathbf{v}$ with finite
norm

$\displaystyle\|\mathbf{v}\|_{L^{p}((a,b),X)}:=\begin{cases}\left(\int_{a}^{b}\|\mathbf{v}(t)\|_{X}^{p}\text{dt}\right)^{\frac{1}{p}}&\text{if
}p\in[1,\infty)\\\
\text{essSup}_{(a,b)}\|\mathbf{v}(\cdot)\|_{X}&\text{otherwise}.\end{cases}$

Finally, for any two positive quantities $a$ and $b$, we write $a\lesssim b$
and $a\gtrsim b$ if there exists a positive constant $c$ such that $a\leq
c\,b$ and $a\geq c\,b$, respectively. We also write $a\approx b$ if $a\lesssim
b$ and $b\lesssim a$. We require the constant $c$ to be independent of the
discretization parameters. In the following proofs, the explanation of the
identities and upper and lower bounds will appear either in the preceeding
text or as an equation reference above the equality symbol “$=$” or the
inequality symbols “$\leq$”, “$\geq$” etc, whichever we believe it is easier
for the reader.

##  2 The continuous problem

We consider the strong form of Maxwell’s equations on a polyhedral domain
$\Omega\subset\mathbb{R}^{3}$ with Lipschitz boundary $\partial\Omega$: _Given
the initial data $\mathbf{E}^{0}$ and $\mathbf{B}^{0}$, find the electric
field $\mathbf{E}$ and the magnetic induction field $\mathbf{B}$ such that_

$\begin{cases}\varepsilon\mathbf{E}_{t}+\sigma\mathbf{E}-\operatorname{\mathbf{curl}}(\mu^{-1}\mathbf{B})=\mathbf{J}&\quad\text{in}~{}\Omega,\,\forall
t\in[0,T],\\\\[5.0pt]
\mathbf{B}_{t}+\operatorname{\mathbf{curl}}\mathbf{E}=\mathbf{0}&\quad\text{in}~{}\Omega,\,\forall
t\in[0,T],\\\\[5.0pt]
\mathbf{E}(0)=\mathbf{E}^{0},\quad\mathbf{B}(0)=\mathbf{B}^{0}\quad&\quad\text{in}~{}\Omega,\\\\[5.0pt]
\mathbf{E}\times\mathbf{n}_{\Omega}=0,\quad\mathbf{B}\cdot\mathbf{n}_{\Omega}=0&\quad\text{on}~{}\partial\Omega,\end{cases}$
(1)

where the subscript $t$ denotes the first derivative in time (so throughout
the paper we use $q_{t}$ instead of $\partial q/{\partial t}$ for a given
time-dependent quantity $q(t)$). Above, $\mathbf{J}$, $\varepsilon$, $\sigma$,
and $\mu$ denote the electric current density that is externally applied to
the system, the electric permittivity, the electric conductivity, and the
magnetic permeability. We consider homogeneous boundary conditions to ease the
exposition, since the nonhomogeneous boundary case presents further
complications. We assume that the initial magnetic induction is a solenoidal
field, i.e.,

$\operatorname{div}\mathbf{B}^{0}=0.$ (2)

Taking the divergence of the second equation in (1), we readily deduce that

$\displaystyle\operatorname{div}\mathbf{B}(t)=0\quad\quad\forall t\in[0,T].$
(3)

The weak formulation of Maxwell’s equations reads as follows:

$\begin{cases}\text{Find}~{}(\mathbf{E},\mathbf{B})\in
H_{0}(\operatorname{\mathbf{curl}},\Omega)\times
H_{0}(\operatorname{div},\Omega)\text{~{}such~{}that}\\\\[5.0pt]
(\varepsilon\mathbf{E}_{t}+\sigma\mathbf{E},\mathbf{w})_{0,\Omega}-(\mu^{-1}\mathbf{B},\operatorname{\mathbf{curl}}\mathbf{w})_{0,\Omega}=(\mathbf{J},\mathbf{w})_{0,\Omega}&\quad\forall\mathbf{w}\in
H_{0}(\operatorname{\mathbf{curl}},\Omega),\\\\[5.0pt]
(\mu^{-1}\mathbf{B}_{t},{\bm{\psi}})_{0,\Omega}+(\mu^{-1}{\bm{\psi}},\operatorname{\mathbf{curl}}\mathbf{E})_{0,\Omega}=0&\quad\forall{\bm{\psi}}\in
H_{0}(\operatorname{div},\Omega).\end{cases}$ (4)

In the sequel, we shall assume that there exist strictly positive constants
$\sigma^{\star},\varepsilon_{\star},\varepsilon^{\star},\mu_{\star},\mu^{\star}$
such that, for all $x\in\Omega$, the material parameters satisfy

$0\leq\sigma(x)\leq\sigma^{\star}\
,\quad\varepsilon_{\star}\leq\varepsilon\leq\varepsilon^{\star}\
,\quad\mu_{\star}\leq\mu(x)\leq\mu^{\star}\ .$ (5)

#### On the regularity of the solutions to (4).

Under suitable assumptions on the regularity of the data, problem (4) admits a
unique solution; see, e.g., [45, Theorem 2.1] and the references therein. We
here recall sufficient conditions from [45] leading to extra smoothness in
space for the solutions to Maxwell’s equations that will be needed in the
following derivations. To the aim, given
$(\mathbf{u},\mathbf{v})\in\mathcal{H}:=[L^{2}(\Omega)]^{3}\times[L^{2}(\Omega)]^{3}$,
we first introduce an operator $\mathcal{A}$ with domain

$\mathcal{D}(\mathcal{A})=\left\\{\begin{pmatrix}\mathbf{u}\\\
\mathbf{v}\end{pmatrix}\in\mathcal{H}\text{ such that }\mathbf{u}\in
H_{0}(\operatorname{\mathbf{curl}},\Omega)\text{ and
}\operatorname{\mathbf{curl}}(\mu^{-1}\mathbf{v})\in[L^{2}(\Omega)]^{3}\right\\},$

where the operator $\mathcal{A}$ is given by

$\mathcal{A}\begin{pmatrix}\mathbf{u}\\\
\mathbf{v}\end{pmatrix}=\begin{pmatrix}-\varepsilon^{-1}\operatorname{\mathbf{curl}}(\mu^{-1}\mathbf{v})+\varepsilon^{-1}\sigma\mathbf{u}\\\
\operatorname{\mathbf{curl}}\mathbf{u}\end{pmatrix}\qquad\forall\begin{pmatrix}\mathbf{u}\\\
\mathbf{v}\end{pmatrix}\in\mathcal{D}(\mathcal{A}).$

Introduce

$\mathbf{X_{N}}(\Omega,\varepsilon):=\\{\mathbf{u}\in
H_{0}(\operatorname{\mathbf{curl}},\Omega)\mid\varepsilon\mathbf{u}\in
H(\operatorname{div},\Omega)\\},\qquad\mathbf{X_{T}}(\Omega,\mu):=\\{\mathbf{u}\in
H(\operatorname{\mathbf{curl}},\Omega)\mid\mu\mathbf{u}\in
H_{0}(\operatorname{div},\Omega)\\}.$

Let $\mathbf{E}$ and $\mathbf{B}$ be the solutions to (4). Assume that $\mu$
or $\varepsilon$ is a constant function, $\sigma$, $\mu$, and $\varepsilon$
are continuous, and $\sigma/\varepsilon\in W^{1,p}(\Omega)$ with $p>3$.
Further assume

$\begin{split}&\mathbf{J}\in\mathcal{C}^{3}([0,T],[L^{2}(\Omega)]^{3}),\qquad\operatorname{div}\mathbf{J}\in
C^{2}([0,T],[L^{2}(\Omega)]^{3}),\\\
&(\mathbf{J}(0),0)^{T},\;\mathcal{A}(\mathbf{J}(0),0)^{T}+(\mathbf{J}_{t}(0),0)^{T}\in\mathbf{X_{N}}(\Omega,\varepsilon)\times\mu\mathbf{X_{T}}(\Omega,\mu).\end{split}$

Then, as in [45, Theorem $2.3$], we have that $\mathbf{E}(t)$,
$\mathbf{B}(t)$, $\mathbf{E}_{t}(t)$,
$\operatorname{\mathbf{curl}}\mathbf{E}(t)$, and
$\operatorname{\mathbf{curl}}\mathbf{E}_{t}$ belong to $H^{s}(\Omega)$, for
some $s>1/2$, for all $t\in[0,T]$.

##  3 The virtual element method

In this section, we construct the virtual element method for the variational
formulation of Maxwell’s equations (4) and discuss its main properties. We
formulate the VEM on sequences of polyhedral meshes, whose properties are
discussed in Section 3.1. In Sections 3.2–3.4, we briefly review the
definitions of the lowest-order nodal, edge and face virtual element spaces
and the design of computable discrete bilinear forms. In Section 3.5, we
recall from [11] that these spaces form an exact de Rham sequence and review
some related property. The design of the virtual element spaces follows the
guidelines of [11]; see also [12, 9, 10]. In Sections 3.7 and 3.8, we present
the semi-discrete and fully-discrete method.

###  3.1 Polyhedral meshes and mesh assumptions

Let $\\{\Omega_{h}\\}_{h}$ be a sequence of mesh partitionings of the
computational domain $\Omega$ labeled by the subscript $h$, which stands for
the _mesh size parameter_. Every mesh $\Omega_{h}$ is a collection of open,
bounded, simply connected polyhedral elements $K$ such that
$\overline{\Omega}=\bigcup_{K\in\Omega_{h}}\overline{K}$. The mesh elements
are nonoverlapping in the sense that the intersection of any possible pair of
them can only be either the empty set, a set of common vertices, or a shared
portion of their boundaries (which is a union of edges). The mesh size
parameter is defined as $h=\max_{K\in\Omega_{h}}h_{K}$, where
$h_{K}=\sup_{\mathbf{x},\mathbf{x}^{\prime}\in
K}|\mathbf{x}-\mathbf{x}^{\prime}|$ is the _diameter_ of $K$. Other
characteristic lengths are the face diameters
$h_{F}=\sup_{\mathbf{x},\mathbf{x}^{\prime}\in
F}|\mathbf{x}-\mathbf{x}^{\prime}|$, which are defined for any mesh face $F$,
and the edge lengths $h_{e}$, which are defined for any mesh edge $e$. For all
$h$, we denote the set of faces and edges by $\mathcal{F}_{h}$ and
$\mathcal{E}_{h}$. Moreover, we denote the set of faces of an element
$K\in\Omega_{h}$ by $\mathcal{F}^{K}$ and the set of edges of a face
$F\in\mathcal{F}_{h}$ by $\mathcal{E}^{F}$. Consistently with our previous
notation, $\mathbf{n}_{K}$ is the unit vector pointing out of element $K$, and
$\mathbf{b}_{K}$ and $\mathbf{b}_{F}$ are the centroids of $K$ and $F$.

Let $\gamma>0$. A face $F$ is said to be _$\gamma$ -shape regular_ if there
exists a a two-dimensional ball $B$ with diameter $h_{B}$ in the interior of
$F$ such that $h_{B}>\gamma h_{F}$. Similarly, an element $K$ is said to be
_$\gamma$ -shape regular_ if there exists a three-dimensional ball $B$ with
diameter $h_{B}$ in the interior of $K$ such that $h_{B}>\gamma h_{K}$.

In the rest of the manuscript we assume that all the meshes $\Omega_{h}$ of a
given sequence $\\{\Omega_{h}\\}$ satisfy these conditions uniformly: there
exists a real constant factor $\gamma\in(0,1)$ independent of $h$ such that

  * •

all the elements $K\in\Omega_{h}$ and faces $F\in\mathcal{F}_{h}$ are
$\gamma$-shape regular;

  * •

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

We assume that the (scalar and real valued) problem coefficients
$\varepsilon$, $\sigma$, and $\mu$ in (1) and (4) are piecewise continuous
over $\Omega_{h}$. As a consequence, we can approximate them by the three
piecewise constant functions $\widehat{\varepsilon}$, $\widehat{\sigma}$, and
$\widehat{\mu}$ given by, in every mesh element $K\in\Omega_{h}$,

$\widehat{\varepsilon}_{K}:=\varepsilon(\mathbf{b}_{K}),\quad\widehat{\sigma}_{K}:=\sigma(\mathbf{b}_{K}),\quad\widehat{\mu}_{K}:=\mu(\mathbf{b}_{K}).$
(6)

To perform the analysis of the method, we also need the additional regularity
condition that, for every element $K\in\Omega_{h}$,

${\mu^{-1}}_{|K},\,\,\,{\sigma}_{|K},\,\,{\varepsilon}_{|K}\in
W^{1,\infty}(K).$ (7)

On every mesh $\Omega_{h}$, we consider the broken Sobolev space of order
$s\geq 0$

$\displaystyle
H^{s}(\Omega_{h}):=\Big{\\{}\mathbf{v}\in[L^{2}(\Omega)]^{3}\mid{\mathbf{v}}_{|K}\in[H^{s}(K)]^{3}\,\,\forall
K\in\Omega_{h}\Big{\\}},$

endowed with the seminorm

$\displaystyle|\mathbf{v}|_{s,\Omega_{h}}^{2}:=\sum_{K\in\Omega_{h}}|\mathbf{v}|^{2}_{s,K}.$

For all elements $K$, we define the local $L^{2}$-orthogonal projector
$\bm{\Pi}^{0}_{h}:[L^{2}(K)]^{3}\rightarrow[\mathbb{P}_{0}(K)]^{3}$ onto
constant vectors as

$\begin{split}(\bm{\Pi}^{0}_{h}\mathbf{E},\mathbf{c})_{0,K}=(\mathbf{E},\mathbf{c})_{0,K}\quad\quad\forall\mathbf{E}\in[L^{2}(K)]^{3},\,\forall\mathbf{c}\in[\mathbb{P}_{0}(K)]^{3}.\end{split}$
(8)

Given a function $\mathbf{E}\in[H^{s}(K)]^{3}$, $0<s\leq 1$, we have the
standard approximation property

$\|\mathbf{E}-\bm{\Pi}^{0}_{h}\mathbf{E}\|_{0,K}\lesssim
h_{K}^{s}|\mathbf{E}|_{s,K}.$

###  3.2 Nodal virtual element spaces

Consider a mesh face $F\in\mathcal{F}_{h}$ and set

$\mathbf{x}_{F}=\mathbf{x}-\mathbf{b}_{F}\quad\quad\forall\mathbf{x}\in F.$
(9)

The nodal virtual element space on face $F$ is

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

We use $V^{\operatorname{node}}_{h}(F)$ in the definition of the nodal virtual
element space on element $K\in\Omega_{h}$, which is given by

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

Every virtual element function $v_{h}\in V^{\operatorname{node}}_{h}(K)$ is
uniquely characterized by the set of its values at the vertices of $K$, which
we take as the degrees of freedom . This unisolvence property is proved in
[11]. Then, we introduce the global virtual element space of the functions
that are globally defined on the computational domain $\Omega$ and have zero
trace on $\partial\Omega$:

$\displaystyle V^{\operatorname{node}}_{h}:=\Big{\\{}v_{h}\in
H^{1}_{0}(\Omega)\mid{v_{h}}_{|K}\in V^{\operatorname{node}}_{h}(K)\,\,\forall
K\in\Omega_{h}\Big{\\}}.$

The degrees of freedom of $V^{\operatorname{node}}_{h}$ are given by an
$H^{1}$-conforming coupling of the local degrees of freedom, i.e., collecting
all the internal vertex values.

For future reference, we also define the interpolant $v_{I}\in
V^{\operatorname{node}}_{h}$ of a function
$v\in\mathcal{C}^{0}(\overline{\Omega})$ as the unique nodal virtual element
function with the same vertex values as $v$. Formally, this definition reads
as

$\displaystyle v_{I}(\nu):=v(\nu)\quad\quad\text{for all vertices~{}$\nu$
in~{}$\Omega_{h}$.}$ (10)

Upper bounds for the approximation error $v-v_{I}$ are available in the
literature; see, e.g., [38, 15, 18].

###  3.3 Edge virtual element spaces

#### Space definitions.

The edge virtual element space on face $F\in\mathcal{F}_{h}$ is

$\begin{split}\mathbf{V}^{\operatorname{edge}}_{h}(F):=\bigg{\\{}\mathbf{v}_{h}\in[L^{2}(F)]^{2}\mid&\operatorname{div_{F}}\mathbf{v}_{h}\in\mathbb{P}_{0}(F),\;\operatorname{rot_{F}}\mathbf{v}_{h}\in\mathbb{P}_{0}(F),\\\
&\,\mathbf{v}_{h}\cdot\mathbf{t}_{e}\in\mathbb{P}_{0}(e)\,\,\forall
e\in\mathcal{E}^{F},\;\int_{F}\mathbf{v}_{h}\cdot\mathbf{x}_{F}=0\bigg{\\}},\end{split}$
(11)

where $\mathbf{x}_{F}$ is defined as in (9). Next, we define the edge virtual
element space on an element $K\in\Omega_{h}$ as

$\begin{split}\mathbf{V}^{\operatorname{edge}}_{h}(K):=\bigg{\\{}\mathbf{v}_{h}&\in[L^{2}(K)]^{3}\mid\operatorname{div}\mathbf{v}_{h}=0,\;\operatorname{\mathbf{curl}}\operatorname{\mathbf{curl}}\mathbf{v}_{h}\in[\mathbb{P}_{0}(K)]^{3},\\\
&(\mathbf{n}_{F}\times\mathbf{v}_{h}{}_{|F})\times\mathbf{n}_{F}\in\mathbf{V}^{\operatorname{edge}}_{h}(F)\,\,\forall
F\in\mathcal{F}^{K},\quad\mathbf{v}_{h}\cdot\mathbf{t}_{e}\text{~{}continuous
at each edge~{}$e$,~{}}\\\\[10.00002pt]
&\int_{K}\operatorname{\mathbf{curl}}\mathbf{v}_{h}\cdot(\mathbf{x}_{K}\times\mathbf{p}_{0})=0\,\,\forall\mathbf{p}_{0}\in[\mathbb{P}_{0}(K)]^{3}\bigg{\\}},\end{split}$
(12)

where $\mathbf{x}_{K}=\mathbf{x}-\mathbf{b}_{K}$ for all $\mathbf{x}\in K$. We
note that $(\mathbf{n}_{F}\times\mathbf{v}_{h}{}_{|F})\times\mathbf{n}_{F}$
corresponds to the projection of $\mathbf{v}_{h}{}_{|F}$ onto the tangent
plane to $F$. The last integral constraints in (11) and (12) are required to
allow for the computation of the $L^{2}$ projector onto vector constant
functions defined in (8); see [11, Section 4.1.2] for more details.

Every virtual element function
$\mathbf{v}_{h}\in\mathbf{V}^{\operatorname{edge}}_{h}(K)$ is uniquely
characterized by the constant values of $\mathbf{v}_{h}\cdot\mathbf{t}_{e}$ on
the edges $e$ of $K$, which we take as the degrees of freedom. The unisolvence
of this set of degrees of freedom for the space
$\mathbf{V}^{\operatorname{edge}}_{h}(K)$ is proved, e.g., in [12]. A
noteworthy property of the local edge virtual element space is that the
$L^{2}$-orthogonal projector $\bm{\Pi}^{0}_{h}$ defined in (8) is computable
from the degrees of freedom of the edge virtual element functions; see, e.g.,
[12]. We define the global edge virtual element space as

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

This definition includes the homogeneous boundary conditions on
$\partial\Omega$. The global set of degrees of freedom of
$\mathbf{V}^{\operatorname{edge}}_{h}$ is obtained via an
$H(\operatorname{\mathbf{curl}})$-conforming coupling of the local ones.

#### Bilinear forms.

As customary in the virtual element framework, we introduce local computable
bilinear forms

$[\cdot,\cdot]_{\operatorname{edge},K}:\mathbf{V}^{\operatorname{edge}}_{h}(K)\times\mathbf{V}^{\operatorname{edge}}_{h}(K)\rightarrow\mathbb{R},$

mimicking the $L^{2}$ inner product $(\cdot,\cdot)_{0,K}$. In particular, we
first introduce the stabilizing bilinear form
$S^{K}_{\operatorname{edge}}:\mathbf{V}^{\operatorname{edge}}_{h}(K)\times\mathbf{V}^{\operatorname{edge}}_{h}(K)\rightarrow\mathbb{R}$
satisfying

$S^{K}_{\operatorname{edge}}(\mathbf{v}_{h},\mathbf{v}_{h})\approx\|\mathbf{v}_{h}\|^{2}_{0,K}\quad\quad\forall\mathbf{v}_{h}\in\ker(\bm{\Pi}^{0}_{h})\cap\mathbf{V}^{\operatorname{edge}}_{h}(K).$
(13)

Then, we define the local discrete counterpart of the $L^{2}$ inner product as

$\big{[}\mathbf{u}_{h},\mathbf{v}_{h}\big{]}_{\operatorname{edge},K}:=(\bm{\Pi}^{0}_{h}\mathbf{u}_{h},\bm{\Pi}^{0}_{h}\mathbf{v}_{h})_{0,K}+S^{K}_{\operatorname{edge}}((\mathbf{I}-\bm{\Pi}^{0}_{h})\mathbf{u}_{h},(\mathbf{I}-\bm{\Pi}^{0}_{h})\mathbf{v}_{h})\quad\forall\mathbf{u}_{h},\,\mathbf{v}_{h}\in\mathbf{V}^{\operatorname{edge}}_{h}(K).$
(14)

The local discrete bilinear forms $[\cdot,\cdot]_{\operatorname{edge},K}$
satisfies the stability condition

$\big{[}\mathbf{v}_{h},\mathbf{v}_{h}\big{]}_{\operatorname{edge},K}\approx\|\mathbf{v}_{h}\|^{2}_{0,K}\quad\quad\forall\mathbf{v}_{h}\in\mathbf{V}^{\operatorname{edge}}_{h}(K)$
(15)

and the consistency property

$\big{[}\mathbf{p}_{0},\mathbf{v}_{h}\big{]}_{\operatorname{edge},K}=(\mathbf{p}_{0},\mathbf{v}_{h})_{0,K}\quad\quad\forall\mathbf{p}_{0}\in[\mathbb{P}_{0}(K)]^{3},\,\mathbf{v}_{h}\in\mathbf{V}^{\operatorname{edge}}_{h}(K).$
(16)

Whereas property (16) follows from definition (14), property (15) requires the
design of a suitable stabilization satisfying (13). If we consider the
stabilization

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

proposed in [11, formula (4.8)], then the stability bounds (15) are proven in
[7, Proposition 5.5].

We introduce the global discrete bilinear forms

$\displaystyle[\widehat{\varepsilon}\mathbf{u}_{h},\mathbf{v}_{h}]_{\operatorname{edge}}$
$\displaystyle:=\sum_{K\in\Omega_{h}}\widehat{\varepsilon}_{|K}[\mathbf{u}_{h},\mathbf{v}_{h}]_{\operatorname{edge},K}\quad\quad\forall\mathbf{u}_{h},\,\mathbf{v}_{h}\in\mathbf{V}^{\operatorname{edge}}_{h},$
(18) and
$\displaystyle[\widehat{\sigma}\mathbf{u}_{h},\mathbf{v}_{h}]_{\operatorname{edge}}$
$\displaystyle:=\sum_{K\in\Omega_{h}}\widehat{\sigma}_{|K}[\mathbf{u}_{h},\mathbf{v}_{h}]_{\operatorname{edge},K}\quad\quad\forall\mathbf{u}_{h},\,\mathbf{v}_{h}\in\mathbf{V}^{\operatorname{edge}}_{h}.$
(19)

In these definitions, we scale the local bilinear forms in the right-hand side
of (14) by $\widehat{\varepsilon}$ and $\widehat{\sigma}$. Moreover, in the
forthcoming analysis, we shall employ the mesh-dependent norm

$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathbf{v}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}^{2}:=[\mathbf{v}_{h},\mathbf{v}_{h}]_{\operatorname{edge}},$

which is induced by the scalar product defined in (18) by setting
$\widehat{\varepsilon}=1$.

#### Interpolation properties.

We denote the interpolation in $\mathbf{V}^{\operatorname{edge}}_{h}$ of a
given vector-valued field $\mathbf{v}\in
H^{s}(\operatorname{\mathbf{curl}},\Omega)$, $1/2<s\leq 1$ by
$\mathbf{v}_{I}$. By definition, $\mathbf{v}_{I}$ is the only function in
$\mathbf{V}^{\operatorname{edge}}_{h}$ such that

$\int_{e}(\mathbf{v}-\mathbf{v}_{I})\cdot\mathbf{t}_{e}=0\quad\quad\forall
e\in\mathcal{E}_{h}.$ (20)

We recall the following interpolation result; see [7, Proposition 4.5;
Corollary 4.6].

###### Proposition 1

Let $\mathbf{v}\in H^{s}(\operatorname{\mathbf{curl}},\Omega)$, $1/2<s\leq 1$,
and $\mathbf{v}_{I}\in\mathbf{V}^{\operatorname{edge}}_{h}$ be its interpolant
as defined in (20). Then, for all $K\in\Omega_{h}$, we have that that

$\begin{split}\|\mathbf{v}-\mathbf{v}_{I}\|_{0,K}&\lesssim
h_{K}^{s}|\mathbf{v}|_{s,K}+h_{K}\|\operatorname{\mathbf{curl}}\mathbf{v}\|_{0,K}+h_{K}^{s+1}|\operatorname{\mathbf{curl}}\mathbf{v}|_{s,K},\\\\[5.0pt]
\|\operatorname{\mathbf{curl}}(\mathbf{v}-\mathbf{v}_{I})\|_{0,K}&\lesssim
h_{K}^{s}|\operatorname{\mathbf{curl}}\mathbf{v}|_{s,K}.\end{split}$ (21)

We shall use inequalities (21) in the forthcoming analysis of Sections 4 and
5.

###  3.4 Face virtual element spaces

#### Space definitions.

Given an element $K\in\Omega_{h}$, we define the face virtual element space as

$\displaystyle\begin{split}\mathbf{V}^{\operatorname{face}}_{h}(K):=\bigg{\\{}\bm{\psi}_{h}&\in[L^{2}(K)]^{3}\mid\operatorname{div}\bm{\psi}_{h}\in\mathbb{P}_{0}(K),\;\operatorname{\mathbf{curl}}\bm{\psi}_{h}\in[\mathbb{P}_{0}(K)]^{3},\\\
&\bm{\psi}_{h}\cdot\mathbf{n}_{F}\in\mathbb{P}_{0}(F)\,\forall
F\in\mathcal{E}^{K},\\\
&\int_{K}\bm{\psi}_{h}\cdot(\mathbf{x}_{K}\times\mathbf{p}_{0})=0\
\forall\mathbf{p}_{0}\in[\mathbb{P}_{0}(K)]^{3}\bigg{\\}},\end{split}$

where we recall that $\mathbf{x}_{K}$ is defined as
$\mathbf{x}-\mathbf{b}_{K}$ for all $\mathbf{x}\in K$. Every virtual element
function $\bm{\psi}_{h}\in\mathbf{V}^{\operatorname{face}}_{h}(K)$ is uniquely
characterized by the constant values of $\bm{\psi}_{h}\cdot\mathbf{n}_{F}$ on
the faces $F$ of $\mathcal{F}^{K}$, which we take as the degrees of freedom.
The unisolvence of this set of degrees of freedom for the space
$\mathbf{V}^{\operatorname{face}}_{h}(K)$ is proved in [12]. A noteworthy
property of the local face virtual element space is that the
$L^{2}$-orthogonal projector $\bm{\Pi}^{0}_{h}$ defined in (8) is computable
from the degrees of freedom of the face virtual element functions; see, e.g.,
[12]. We define the global face virtual element space as

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

This definition includes homogeneous boundary conditions on $\partial\Omega$.
The set of degrees of freedom of $\mathbf{V}^{\operatorname{face}}_{h}$ is
obtained via an $H(\operatorname{div})$-conforming coupling of the local ones,
i.e., by collecting together the internal degrees of freedom.

#### Bilinear forms.

As in the edge element case, we introduce local computable bilinear forms

$\big{[}\cdot,\cdot\big{]}_{\operatorname{face},K}:\mathbf{V}^{\operatorname{face}}_{h}(K)\times\mathbf{V}^{\operatorname{face}}_{h}(K)\rightarrow\mathbb{R},$

mimicking the $L^{2}$ inner product $(\cdot,\cdot)_{0,K}$ on
$\mathbf{V}^{\operatorname{face}}_{h}(K)$. In particular, we first introduce
the stabilizing bilinear form
$S^{K}_{\operatorname{face}}:\mathbf{V}^{\operatorname{face}}_{h}(K)\times\mathbf{V}^{\operatorname{face}}_{h}(K)\rightarrow\mathbb{R}$
satisfying

$S^{K}_{\operatorname{face}}(\bm{\psi}_{h},\bm{\psi}_{h})\approx\|\bm{\psi}_{h}\|^{2}_{0,K}\quad\quad\forall\bm{\psi}_{h}\in\ker(\bm{\Pi}^{0}_{h})\cap\mathbf{V}^{\operatorname{face}}_{h}(K).$
(22)

Then, we define the local discrete counterpart of the $L^{2}$ inner product as

$\big{[}\bm{\psi}_{h},\bm{\phi}_{h}\big{]}_{\operatorname{face},K}:=(\bm{\Pi}^{0}_{h}\bm{\psi}_{h},\bm{\Pi}^{0}_{h}\bm{\phi}_{h})_{0,K}+S^{K}_{\operatorname{face}}((\mathbf{I}-\bm{\Pi}^{0}_{h})\bm{\psi}_{h},(\mathbf{I}-\bm{\Pi}^{0}_{h})\bm{\phi}_{h})\quad\forall\bm{\psi}_{h},\,\bm{\phi}_{h}\in\mathbf{V}^{\operatorname{face}}_{h}(K).$
(23)

The local discrete bilinear forms $[\cdot,\cdot]_{\operatorname{face},K}$
satisfies the stability condition

$\big{[}\bm{\psi}_{h},\bm{\psi}_{h}\big{]}_{\operatorname{face},K}\approx\|\bm{\psi}_{h}\|^{2}_{0,K}\quad\quad\forall\bm{\psi}_{h}\in\mathbf{V}^{\operatorname{face}}_{h}(K)$
(24)

and the consistency property

$\big{[}\mathbf{p}_{0},\bm{\psi}_{h}\big{]}_{\operatorname{face},K}=(\mathbf{p}_{0},\bm{\psi}_{h})_{0,K}\quad\quad\forall\mathbf{p}_{0}\in[\mathbb{P}_{0}(K)]^{3},\,\bm{\psi}_{h}\in\mathbf{V}^{\operatorname{face}}_{h}(K).$
(25)

Whereas (25) follows from definition (23), property (24) requires the design
of a suitable stabilization satisfying (22). We consider the stabilization,
cf. [11, (4.17)],

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

The stability bounds (24) are proven in [7, Proposition 5.2]. Finally, we
introduce the global discrete bilinear form

$\big{[}\widehat{\mu}^{-1}\bm{\psi}_{h},\bm{\phi}_{h}\big{]}_{\operatorname{face}}:=\sum_{K\in\Omega_{h}}\widehat{\mu}^{-1}_{|K}[\bm{\psi}_{h},\bm{\phi}_{h}]_{\operatorname{face},K}\quad\quad\forall\bm{\psi}_{h},\,\bm{\phi}_{h}\in\mathbf{V}^{\operatorname{face}}_{h}.$
(27)

In the forthcoming analysis, we shall employ the mesh-dependent norm

$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\bm{\psi}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{face}}^{2}:=\big{[}\bm{\psi}_{h},\bm{\psi}_{h}\big{]}_{\operatorname{face}}\quad\forall\bm{\psi}_{h}\in\mathbf{V}^{\operatorname{face}}_{h},$

which is induced by the scalar product defined in (27) by setting
$\widehat{\mu}=1$.

#### Interpolation properties.

We denote the interpolation in $\mathbf{V}^{\operatorname{face}}$ of a given
vector-valued field ${\bm{\psi}}\in[H^{s}(\Omega)]^{3}$, $1/2<s\leq 1$ by
$\bm{\psi}_{I}$. By definition, $\bm{\psi}_{I}$ is the only function in
$\mathbf{V}^{\operatorname{face}}$ such that

$\int_{F}({\bm{\psi}}-\bm{\psi}_{I})\cdot\mathbf{n}_{F}=0\quad\quad\forall
F\in\mathcal{F}_{h}.$ (28)

We recall the following interpolation result; see [7, Proposition 3.2;
Corollary 3.3].

###### Proposition 2

Let ${\bm{\psi}}\in[H^{s}(\operatorname{div},\Omega)]^{3}$, $1/2<s\leq 1$, and
$\bm{\psi}_{I}$ be defined as in (28). Then, for all $K\in\Omega_{h}$, we find
that

$\begin{split}&\|{\bm{\psi}}-\bm{\psi}_{I}\|_{0,K}\lesssim
h_{K}^{s}|{\bm{\psi}}|_{s,K},\quad\quad\|\operatorname{div}({\bm{\psi}}-\bm{\psi}_{I})\|_{0,K}\lesssim
h_{K}^{s}|\operatorname{div}\psi|_{s,K}.\end{split}$ (29)

We shall use inequalities (29) in Sections 4 and 5.

###  3.5 Exact sequence properties

We set

$\mathbf{\widetilde{V}}_{h}^{\operatorname{edge}}:=\left\\{\mathbf{v}_{h}\in\mathbf{V}^{\operatorname{edge}}_{h}\mid\operatorname{\mathbf{curl}}(\mathbf{v}_{h})=0\right\\}.$
(30)

As observed in [11, equation (4.33)], the following identity is valid:

$\mathbf{\widetilde{V}}_{h}^{\operatorname{edge}}=\nabla(V^{\operatorname{node}}_{h}).$
(31)

Analogously, we set

$\mathbf{\widetilde{V}}_{h}^{\operatorname{face}}:=\left\\{\bm{\psi}_{h}\in\mathbf{V}^{\operatorname{face}}_{h}\mid\operatorname{div}\bm{\psi}_{h}=0\right\\}.$
(32)

As observed in [11, equation (4.35)], the following identity is also valid:

$\mathbf{\widetilde{V}}_{h}^{\operatorname{face}}=\operatorname{\mathbf{curl}}(\mathbf{V}^{\operatorname{edge}}_{h}).$
(33)

In particular, the spaces $V^{\operatorname{node}}_{h}$,
$\mathbf{V}^{\operatorname{edge}}_{h}$, and
$\mathbf{V}^{\operatorname{face}}_{h}$ form an exact sequence; see, e.g., [12,
11].

###### Remark 1

As shown in [11], the following commuting diagram properties are valid:

  * •

$\operatorname{div}\mathbf{B}_{\mathbf{I}}=\Pi^{0}_{h}(\operatorname{div}\mathbf{B})$
for a given vector-valued field $\mathbf{B}\in
H^{s}(\operatorname{div},\Omega)$, $1/2<s\leq 1$, where $\Pi^{0}_{h}$ is the
scalar version of the $L^{2}$ projector in (8) and $\mathbf{B}_{\mathbf{I}}$
is the face interpolation of $\mathbf{B}$ defined in (28);

  * •

$\operatorname{\mathbf{curl}}\mathbf{E}_{\mathbf{I}}=(\operatorname{\mathbf{curl}}\mathbf{E})_{I}$
for a given vector-valued field $\mathbf{E}\in
H^{s}(\operatorname{\mathbf{curl}},\Omega)$, $1/2<s\leq 1$, where on the left-
and right-hand sides we consider the edge and face interpolations of
$\mathbf{E}$ and $\operatorname{\mathbf{curl}}(\mathbf{E})$ that are
respectively defined in (20) and (28);

  * •

$\nabla v_{I}=(\nabla v)_{I}$ for a given scalar field $v\in H^{1}(\Omega)$,
where on the left- and right-hand sides we consider the nodal and edge
interpolations of $v$ and $\nabla v$ that are respectively defined in (10) and
(20).

###  3.6 Two novel operators

We introduce two novel operators on the spaces
$\mathbf{V}^{\operatorname{edge}}_{h}$ and
$\mathbf{V}^{\operatorname{face}}_{h}$, which also satisfy a commuting diagram
property; see Proposition 6 below. We begin by defining the weighted, global
projector
$\bm{\Pi}_{h}^{\text{curl}}:H(\operatorname{\mathbf{curl}},\Omega)\rightarrow\mathbf{V}^{\operatorname{edge}}_{h}$
as

$\begin{split}&\begin{cases}[\widehat{\mu}^{-1}\operatorname{\mathbf{curl}}(\bm{\Pi}_{h}^{\text{curl}}\mathbf{E}),\operatorname{\mathbf{curl}}\mathbf{w}_{h}]_{\operatorname{face}}=(\mu^{-1}\operatorname{\mathbf{curl}}\mathbf{E},\operatorname{\mathbf{curl}}\mathbf{w}_{h})_{0,\Omega}\quad\forall\mathbf{w}_{h}\in\mathbf{V}^{\operatorname{edge}}_{h}\\\\[5.0pt]
[\bm{\Pi}_{h}^{\text{curl}}\mathbf{E},\nabla
s_{h}]_{\operatorname{edge}}=(\mathbf{E},\nabla
s_{h})_{0,\Omega}\quad\forall\mathbf{E}\in
H(\operatorname{\mathbf{curl}},\Omega),\,\forall s_{h}\in
V^{\operatorname{node}}_{h}.\end{cases}\end{split}$ (34)

The discrete bilinear form $[\cdot,\cdot]_{\operatorname{face}}$ appearing in
the first equation of (34) is well defined. In fact, thanks to (33),
$\operatorname{\mathbf{curl}}(\bm{\Pi}_{h}^{\text{curl}}(\mathbf{E}))$ belongs
to $\mathbf{V}^{\operatorname{face}}_{h}(K)$. An analogous observation applies
for the form $[\cdot,\cdot]_{\operatorname{edge}}$ appearing in the second
equation.

In order to prove the approximation properties of the projector
$\bm{\Pi}_{h}^{\text{curl}}$, we need two preliminary technical results.

###### Lemma 3

Consider $K\in\Omega_{h}$. Then it exists a real parameter $1/2<s\leq 1$,
depending on the shape regularity constant of $K$, such that the following
inverse inequality is valid:

$|\bm{\psi}_{h}|_{s,K}\lesssim
h_{K}^{-s}\|\bm{\psi}_{h}\|_{0,K}\quad\quad\bm{\psi}_{h}\in\mathbf{V}^{\operatorname{face}}_{h}(K).$
(35)

Proof. Define $\Psi$ as the solution to

$\displaystyle\begin{cases}\Delta\Psi=0&\text{in }K\\\
\mathbf{n}_{K}\cdot\nabla\Psi=\mathbf{n}_{K}\cdot\bm{\psi}_{h}&\text{on
}\partial K\\\ \int_{K}\Psi=0.\end{cases}$

Standard regularity results for elliptic problems, see, e.g., [24, Corollary
$23.5$], entail that there exists $1/2<s\leq 1$ such that

$|\Psi|_{s+1,K}\lesssim\|\mathbf{n}_{K}\cdot\bm{\psi}_{h}\|_{s-\frac{1}{2},\partial
K}.$ (36)

Define $\mathbf{z}:=\bm{\psi}_{h}-\nabla\Psi$ and observe that
$\mathbf{n}_{K}\cdot\mathbf{z}=0$ on $\partial K$. We recall from [5,
Proposition 3.7] that for all $\mathbf{z}\in
H(\operatorname{\mathbf{curl}},K)\cap H(\operatorname{div},K)$ with
$\mathbf{n}_{K}\cdot\mathbf{z}{}_{|\partial K}=0$, there exists $1/2<s\leq 1$
such that

$|\mathbf{z}|_{s,K}\lesssim
h_{K}^{-s}\|\mathbf{z}\|_{0,K}+h_{K}^{1-s}\|\operatorname{\mathbf{curl}}\mathbf{z}\|_{0,K}+h_{K}^{1-s}\|\operatorname{div}\mathbf{z}\|_{0,K}.$
(37)

Recalling that $\mathbf{n}_{K}\cdot\mathbf{z}=0$ on $\partial K$,
$\operatorname{div}\mathbf{z}=\operatorname{div}\bm{\psi}_{h}\in\mathbb{P}_{0}(K)$,
and $\Psi$ has vanishing integral on $K$, an integration by parts yields

$(\mathbf{z},\nabla\Psi)_{0,K}=-(\operatorname{div}\mathbf{z},\Psi)_{0,K}+(\mathbf{n}_{K}\cdot\mathbf{z},\Psi)_{0,\partial
K}=0\ .$

Therefore

$\|\mathbf{z}\|_{0,K}^{2}=(\mathbf{z},\mathbf{z})_{0,K}=(\mathbf{z},\bm{\psi}_{h})_{0,K}\leq\|\mathbf{z}\|_{0,K}\|\bm{\psi}_{h}\|_{0,K}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color<EMAIL_ADDRESS>
Thus, (37) becomes

$|\mathbf{z}|_{s,K}\lesssim
h_{K}^{-s}\|\bm{\psi}_{h}\|_{0,K}+h_{K}^{1-s}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{(\|\operatorname{div}\mathbf{z}\|_{0,K}+\|\operatorname{\mathbf{curl}}\mathbf{z}\|_{0,K})}}.$
(38)

We take the minimum $s$ such that (36) and (38) are valid. Using the triangle
inequality, (36), (38), the fact that
$\operatorname{\mathbf{curl}}(\nabla\Psi)=\mathbf{0}$ and
$\operatorname{div}(\nabla\Psi)=0$, and the fact that
$\mathbf{n}_{K}\cdot\bm{\psi}_{h}{}_{|F}\in\mathbb{P}_{0}(F)$ for all
$F\in\mathcal{E}^{K}$, we easily obtain that

$\displaystyle\begin{split}|\bm{\psi}_{h}|_{s,K}&\lesssim|\mathbf{z}|_{s,K}+|\Psi|_{s+1,K}\\\
&\lesssim
h_{K}^{1-s}(\|\operatorname{div}\bm{\psi}_{h}\|_{0,K}+\|\operatorname{\mathbf{curl}}\bm{\psi}_{h}\|_{0,K})+\|\mathbf{n}_{K}\cdot\bm{\psi}_{h}\|_{s-\frac{1}{2},\partial
K}+h_{K}^{-s}\|\bm{\psi}_{h}\|_{0,K}.\end{split}$

We are left to show a bound for each term on the right-hand side in terms of
$h_{K}^{-s}\|\bm{\psi}_{h}\|_{0,K}$. We can prove such bounds based on
employing polynomial inverse inequalities as, e.g., in the proofs of [7,
Proposition $4.1$, Proposition $4.2$]. The main ingredients are the fact that
$\operatorname{div}\bm{\psi}_{h}\in\mathbb{P}_{0}(K)$,
$\operatorname{\mathbf{curl}}\bm{\psi}_{h}\in[\mathbb{P}_{0}(K)]^{3}$,
$\mathbf{n}_{K}\cdot\bm{\psi}_{h}{}_{|F}\in\mathbb{P}_{0}(F)$ for all
$F\in\mathcal{E}^{K}$, and inverse inequalities involving bubbles. ∎

The second auxiliary result is a coercivity property on the kernel of edge
functions

$\displaystyle\mathbf{Z}_{h}:=\left\\{\mathbf{v}_{h}\in\mathbf{V}^{\operatorname{edge}}_{h}\mid[\mathbf{v}_{h},\nabla
v_{h}]_{\operatorname{edge}}=0\;\forall v_{h}\in
V^{\operatorname{node}}_{h}\right\\}.$

###### Lemma 4

The following coercivity property on the kernel $\mathbf{Z}_{h}$ is valid:

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\operatorname{\mathbf{curl}}\mathbf{v}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{face}}\gtrsim\|\mathbf{v}_{h}\|_{0,\Omega}\quad\quad\forall\mathbf{v}_{h}\in\mathbf{Z}_{h}.$
(39)

Proof. Given $\mathbf{v}_{h}\in\mathbf{Z}_{h}$, let $p$ be the solution to the
following problem:

$\displaystyle\begin{cases}\text{find }p\in H^{1}(\Omega)\text{ such that
}\int_{\Omega}p=0\text{ and}\\\ (\nabla
p,\nabla\xi)_{0,\Omega}=(\mathbf{v}_{h},\nabla\xi)_{0,\Omega}\quad\forall\xi\in
H^{1}(\Omega).\end{cases}$

Set $\mathbf{v}:=\mathbf{v}_{h}-\nabla p$. We clearly have that

$\operatorname{div}\mathbf{v}=0,\quad\quad\operatorname{\mathbf{curl}}\mathbf{v}=\operatorname{\mathbf{curl}}\mathbf{v}_{h},\quad\quad\mathbf{n}_{\Omega}\cdot\mathbf{v}_{|\partial\Omega}=0.$
(40)

For a $1/2<s\leq 1$ depending on the shape of $\Omega$, using [5, Proposition
$3.7$] gives

$\|\mathbf{v}\|_{s,\Omega}\lesssim\|\operatorname{\mathbf{curl}}\mathbf{v}\|_{0,\Omega}\overset{\eqref{properties:function:v}}{=}\|\operatorname{\mathbf{curl}}\mathbf{v}_{h}\|_{0,\Omega}.$
(41)

Denote the nodal interpolant of $p$ by $p_{I}$; see (10). As in Remark 1, the
edge interpolant of $\nabla p$ in the sense of (20) is such that $(\nabla
p)_{I}=\nabla p_{I}.$ Therefore, the edge interpolant of $\mathbf{v}$ in the
sense of (20) satisfies

$\displaystyle\mathbf{v}_{I}=\mathbf{v}_{h}-(\nabla
p)_{I}=\mathbf{v}_{h}-\nabla p_{I}.$

Next, recalling that $\mathbf{v}_{h}\in\mathbf{Z}_{h}$, we observe that

$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathbf{v}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}_{\operatorname{edge}}=[\mathbf{v}_{h},\mathbf{v}_{h}]_{\operatorname{edge}}=[\mathbf{v}_{h},\mathbf{v}_{I}+\nabla
p_{I}]_{\operatorname{edge}}=[\mathbf{v}_{h},\mathbf{v}_{I}]_{\operatorname{edge}}\leq{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathbf{v}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathbf{v}_{I}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}.$

We deduce that

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathbf{v}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}\leq{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathbf{v}_{I}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}\overset{\eqref{stability:edge}}{\lesssim}\|\mathbf{v}_{I}\|_{0,\Omega}.$
(42)

We estimate from above the right-hand side of (42) elementwise. Let
$K\in\Omega_{h}$. Using the triangle inequality and (21), we write:

$\begin{split}\|\mathbf{v}_{I}\|_{0,K}&\lesssim\|\mathbf{v}\|_{0,K}+h_{K}^{s}|\mathbf{v}|_{s,K}+h_{K}\|\operatorname{\mathbf{curl}}\mathbf{v}\|_{0,K}+h_{K}^{s+1}|\operatorname{\mathbf{curl}}\mathbf{v}|_{s,K}\\\
&\overset{\eqref{properties:function:v}}{=}\|\mathbf{v}\|_{0,K}+h_{K}^{s}|\mathbf{v}|_{s,K}+h_{K}\|\operatorname{\mathbf{curl}}\mathbf{v}_{h}\|_{0,K}+h_{K}^{s+1}|\operatorname{\mathbf{curl}}\mathbf{v}_{h}|_{s,K}.\\\
\end{split}$

We know that $\operatorname{\mathbf{curl}}\mathbf{v}_{h}$ belongs to
$\mathbf{V}^{\operatorname{face}}(K)$; see (33). Therefore, we can apply the
inverse estimate (35), possibly taking the minimum among the scalar $s$ in
(41) and the minimum over all elements of the parameter $s$ appearing in Lemma
3, and find that

$\|\mathbf{v}_{I}\|_{0,K}\lesssim\|\mathbf{v}\|_{0,K}+h_{K}^{s}|\mathbf{v}|_{s,K}+h_{K}\|\operatorname{\mathbf{curl}}\mathbf{v}_{h}\|_{0,K}.$
(43)

Inserting (43) in (42) and summing over all mesh elements yield

$\begin{split}\|\mathbf{v}_{h}\|_{0,\Omega}\overset{\eqref{stability:edge}}{\lesssim}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathbf{v}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}\lesssim\|\mathbf{v}\|_{0,\Omega}+|\mathbf{v}|_{s,\Omega}+\|\operatorname{\mathbf{curl}}\mathbf{v}_{h}\|_{0,\Omega}\overset{\eqref{using:Amrouche-
Friedrichs}}{\lesssim}\|\operatorname{\mathbf{curl}}\mathbf{v}_{h}\|_{0,\Omega}\overset{\eqref{stability:face}}{\lesssim}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\operatorname{\mathbf{curl}}\mathbf{v}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{face}},\end{split}$

which is the assertion. ∎

We are in the position of proving the approximation properties of the
projector $\bm{\Pi}_{h}^{\text{curl}}$.

###### Proposition 5

Let $\mathbf{E}\in H^{s}(\operatorname{\mathbf{curl}},\Omega)$, $1/2<s\leq 1$,
and $\bm{\Pi}_{h}^{\text{curl}}$ be the projector defined in (34). Then, the
following inequality is valid:

$\|\mathbf{E}-\bm{\Pi}_{h}^{\text{curl}}\mathbf{E}\|_{\operatorname{\mathbf{curl}},\Omega}\lesssim
h^{s}(|\mathbf{E}|_{s,\Omega}+h^{1-s}\|\operatorname{\mathbf{curl}}\mathbf{E}\|_{0,\Omega}+h|\operatorname{\mathbf{curl}}\mathbf{E}|_{s,\Omega}).$
(44)

Proof. Let $\mathbf{E}_{\mathbf{I}}$ denote the interpolant of $\mathbf{E}$ in
the edge space $\mathbf{V}^{\operatorname{edge}}_{h}$; cf. (20). Then, for a
given $\mathbf{E}\in H(\operatorname{\mathbf{curl}},\Omega)$,
$\bm{\Pi}_{h}^{\text{curl}}\mathbf{E}$ is the solution to the following mixed
variational problem:

$\begin{cases}[\widehat{\mu}^{-1}\operatorname{\mathbf{curl}}(\bm{\Pi}_{h}^{\text{curl}}\mathbf{E}),\operatorname{\mathbf{curl}}\mathbf{w}_{h}]_{\operatorname{face}}+[\mathbf{w}_{h},\nabla
p_{h}]_{\operatorname{edge}}=(\mu^{-1}\operatorname{\mathbf{curl}}\mathbf{E},\operatorname{\mathbf{curl}}\mathbf{w}_{h})_{0,\Omega},\\\\[5.0pt]
[\bm{\Pi}_{h}^{\text{curl}}\mathbf{E},\nabla
s_{h}]_{\operatorname{edge}}=(\mathbf{E},\nabla
s_{h})_{0,\Omega}\qquad\forall\mathbf{w}_{h}\in\mathbf{V}^{\operatorname{edge}}_{h},\;\;\forall
s_{h}\in V^{\operatorname{node}}_{h}.\end{cases}$ (45)

Indeed, it can be easily shown that $p_{h}=0$ and the coercivity of the
bilinear form
$[\operatorname{\mathbf{curl}}\cdot,\operatorname{\mathbf{curl}}\cdot]_{\operatorname{face}}$
on the kernel $\mathbf{Z}_{h}$ is shown in Lemma 4. On the other hand, the
discrete inf-sup condition for the bilinear form
$[\cdot,\nabla\cdot]_{\operatorname{edge}}:\mathbf{V}^{\operatorname{edge}}_{h}\times
V^{\operatorname{node}}_{h}\rightarrow\mathbb{R}$ is a trivial consequence of
the fact that the virtual element spaces under consideration form an exact
sequence. Therefore, we can use the standard analysis for mixed problems; see,
e.g., [17]. Notably, there exist $\mathbf{w}_{h}$ and $s_{h}$ in
$\mathbf{V}^{\operatorname{edge}}_{h}\times V^{\operatorname{node}}_{h}$ with

$\|\mathbf{w}_{h}\|_{\operatorname{\mathbf{curl}},\Omega}+\|s_{h}\|_{1,\Omega}\leq
1,$ (46)

such that

$\begin{split}&\|\bm{\Pi}_{h}^{\text{curl}}\mathbf{E}-\mathbf{E}_{\mathbf{I}}\|_{\operatorname{\mathbf{curl}},\Omega}=\|\bm{\Pi}_{h}^{\text{curl}}\mathbf{E}-\mathbf{E}_{\mathbf{I}}\|_{\operatorname{\mathbf{curl}},\Omega}+\|p_{h}\|_{1,\Omega}\\\\[5.0pt]
&\qquad\lesssim[\widehat{\mu}^{-1}\operatorname{\mathbf{curl}}(\bm{\Pi}_{h}^{\text{curl}}\mathbf{E}-\mathbf{E}_{\mathbf{I}}),\operatorname{\mathbf{curl}}\mathbf{w}_{h}]_{\operatorname{face}}+[\mathbf{w}_{h},\nabla
p_{h}]_{\operatorname{edge}}+[\bm{\Pi}_{h}^{\text{curl}}\mathbf{E}-\mathbf{E}_{\mathbf{I}},\nabla
s_{h}]_{\operatorname{edge}}\\\\[5.0pt]
&\qquad\overset{p_{h}=0}{=}[\widehat{\mu}^{-1}\operatorname{\mathbf{curl}}(\bm{\Pi}_{h}^{\text{curl}}\mathbf{E}-\mathbf{E}_{\mathbf{I}}),\operatorname{\mathbf{curl}}\mathbf{w}_{h}]_{\operatorname{face}}+[\bm{\Pi}_{h}^{\text{curl}}\mathbf{E}-\mathbf{E}_{\mathbf{I}},\nabla
s_{h}]_{\operatorname{edge}}\\\\[5.0pt] &\qquad\overset{\eqref{curl-projector-
edge}}{=}(\mu^{-1}\operatorname{\mathbf{curl}}\mathbf{E},\operatorname{\mathbf{curl}}\mathbf{w}_{h})_{0,\Omega}-[\widehat{\mu}^{-1}\operatorname{\mathbf{curl}}\mathbf{E}_{\mathbf{I}},\operatorname{\mathbf{curl}}\mathbf{w}_{h}]_{\operatorname{face}}+(\mathbf{E},\nabla
s_{h})_{0,\Omega}-[\mathbf{E}_{\mathbf{I}},\nabla
s_{h}]_{\operatorname{edge}}=:\big{[}I\big{]}+\big{[}II\big{]}.\end{split}$

Since the bounds for the two terms on the right-hand side follow using
standard VE calculations, we address them briefly. Recall the definition of
the projector $\bm{\Pi}^{0}_{h}$ in (8). As for the term $\big{[}I\big{]}$, we
get

$\displaystyle\big{[}I\big{]}$ $\displaystyle\overset{\hskip
8.53581pt\eqref{consistency:face}\hskip
8.53581pt}{=}\sum_{K\in\Omega_{h}}\bigg{(}(\mu^{-1}(\operatorname{\mathbf{curl}}\mathbf{E}-\bm{\Pi}^{0}_{h}(\operatorname{\mathbf{curl}}\mathbf{E})),\operatorname{\mathbf{curl}}\mathbf{w}_{h})_{0,K}-[\widehat{\mu}^{-1}(\operatorname{\mathbf{curl}}\mathbf{E}_{\mathbf{I}}-\bm{\Pi}^{0}_{h}(\operatorname{\mathbf{curl}}\mathbf{E})),\operatorname{\mathbf{curl}}\mathbf{w}_{h}]_{\operatorname{face},K}$
$\displaystyle\quad\quad\quad\quad\quad\quad+((\mu^{-1}-\widehat{\mu}^{-1})\bm{\Pi}^{0}_{h}(\operatorname{\mathbf{curl}}\mathbf{E}),\operatorname{\mathbf{curl}}\mathbf{w}_{h})_{0,K}\bigg{)}$
$\displaystyle\overset{\eqref{stability:face},\eqref{living-on-the-
ball}}{\lesssim}\Bigg{(}\sum_{K\in\Omega_{h}}\bigg{(}\|\operatorname{\mathbf{curl}}\mathbf{E}-\bm{\Pi}^{0}_{h}(\operatorname{\mathbf{curl}}\mathbf{E})\|_{0,K}^{2}+\|\operatorname{\mathbf{curl}}(\mathbf{E}-\mathbf{E}_{\mathbf{I}})\|_{0,K}^{2}$
$\displaystyle\quad\quad\quad\quad\quad\quad+\|\mu^{-1}-\widehat{\mu}^{-1}\|_{L^{\infty}(K)}\|\bm{\Pi}^{0}_{h}(\operatorname{\mathbf{curl}}\mathbf{E})\|_{0,K}^{2}\bigg{)}\Bigg{)}^{\frac{1}{2}}$
$\displaystyle\overset{\hskip
5.69054pt\eqref{assumption:coefficients},\eqref{interpolation-
estimates:edge}}{\lesssim}h^{s}|\operatorname{\mathbf{curl}}\mathbf{E}|_{s,\Omega}+h\max_{K\in\Omega_{h}}|\mu^{-1}|_{W^{1,\infty}(K)}^{2}\|\operatorname{\mathbf{curl}}\mathbf{E}\|_{0,\Omega}.$

We proceed similarly for the term $\big{[}II\big{]}$:

$\displaystyle\big{[}II\big{]}\overset{\hskip
8.53581pt\eqref{consistency:edge}\hskip
8.53581pt}{=}\sum_{K\in\Omega_{h}}\left((\mathbf{E}-\bm{\Pi}^{0}_{h}\mathbf{E},\nabla
s_{h})_{0,\Omega}-[\mathbf{E}_{\mathbf{I}}-\bm{\Pi}^{0}_{h}\mathbf{E},\nabla
s_{h}]_{\operatorname{edge}}\right)$
$\displaystyle\overset{\eqref{stability:face},\eqref{living-on-the-
ball}}{\lesssim}\left(\sum_{K\in\Omega_{h}}\left(\|\mathbf{E}-\bm{\Pi}^{0}_{h}\mathbf{E}\|_{0,K}^{2}+\|\mathbf{E}-\mathbf{E}_{\mathbf{I}}\|_{0,K}^{2}\right)\right)^{\frac{1}{2}}\overset{\hskip
8.53581pt\eqref{interpolation-estimates:edge}\hskip
8.53581pt}{\lesssim}h^{s}\left(|\mathbf{E}|_{s,\Omega}+h^{1-s}\|\operatorname{\mathbf{curl}}\mathbf{E}\|_{0,\Omega}+h|\operatorname{\mathbf{curl}}\mathbf{E}|_{s,\Omega}\right).$

The assertion follows easily collecting the bounds on the terms
$\big{[}I\big{]}$, $\big{[}II\big{]}$ and by the triangle inequality. ∎

Next, define the weighted, global projector
$\mathbf{P}_{h}:[L^{2}(\Omega)]^{3}\rightarrow\mathbf{\widetilde{V}}_{h}^{\operatorname{face}}$
as

$[\widehat{\mu}^{-1}\mathbf{P}_{h}\mathbf{B},\bm{\psi}_{h}]_{\operatorname{face}}=(\mu^{-1}\mathbf{B},\bm{\psi}_{h})_{0,\Omega}\quad\quad\forall\mathbf{B}\in[L^{2}(\Omega)]^{3},\,\bm{\psi}_{h}\in\mathbf{\widetilde{V}}_{h}^{\operatorname{face}}.$
(47)

As, e.g., in [45], a crucial point in the analysis of the semi-discrete scheme
in Section 4 below is the following commuting diagram result.

###### Proposition 6

Let $\bm{\Pi}_{h}^{\text{curl}}$ and $\mathbf{P}_{h}$ be the two projectors
introduced in (34) and (47), respectively. Then, the following identity is
valid:

$\operatorname{\mathbf{curl}}(\bm{\Pi}_{h}^{\text{curl}}\mathbf{E})=\mathbf{P}_{h}(\operatorname{\mathbf{curl}}\mathbf{E})\quad\quad\forall\mathbf{E}\in
H(\operatorname{\mathbf{curl}},\Omega).$ (48)

Proof. Recall that
$\operatorname{\mathbf{curl}}(\bm{\Pi}_{h}^{\text{curl}}\mathbf{E})$ belongs
to $\mathbf{\widetilde{V}}_{h}^{\operatorname{face}}$; see (48). Using (33),
(34), and (47), we get

$\displaystyle[\widehat{\mu}^{-1}\operatorname{\mathbf{curl}}(\bm{\Pi}_{h}^{\text{curl}}\mathbf{E}),\bm{\psi}_{h}]_{\operatorname{face}}=(\mu^{-1}\operatorname{\mathbf{curl}}\mathbf{E},\bm{\psi}_{h})_{0,\Omega}=[\widehat{\mu}^{-1}\mathbf{P}_{h}(\operatorname{\mathbf{curl}}\mathbf{E}),\bm{\psi}_{h}]_{\operatorname{face}}\quad\forall\bm{\psi}_{h}\in\mathbf{\widetilde{V}}_{h}^{\operatorname{face}}.$

The assertion follows using the stability property (24). ∎

In the light of commuting diagram (48), the projector $\mathbf{P}_{h}$
satisfies the following property.

###### Lemma 7

Let $\mathbf{P}_{h}$ be defined in (47). Then, for all divergence free
${\bm{\psi}}\in H^{s}(\Omega)$, $1/2<s\leq 1$, the following bound is valid:

$\|{\bm{\psi}}-\mathbf{P}_{h}{\bm{\psi}}\|_{0,\Omega}\lesssim
h\|{\bm{\psi}}\|_{0,\Omega}+h^{s}|{\bm{\psi}}|_{s,\Omega}.$ (49)

Proof. Let $\bm{\Pi}^{0}_{h}$ be the projector defined in (8) and
$\bm{\psi}_{I}$ the interpolant in $\mathbf{V}^{\operatorname{edge}}_{h}$ of
${\bm{\psi}}$; see (28). Then, we can write

$\displaystyle\|{\bm{\psi}}-\mathbf{P}_{h}{\bm{\psi}}\|_{0,\Omega}^{2}\lesssim\|{\bm{\psi}}-\bm{\psi}_{I}\|_{0,\Omega}^{2}+\|\bm{\psi}_{I}-\mathbf{P}_{h}{\bm{\psi}}\|_{0,\Omega}^{2}.$

Since $\operatorname{div}{\bm{\psi}}=0$, we have that
$\bm{\psi}_{I}\in\mathbf{\widetilde{V}}_{h}^{\operatorname{face}}$ and thus
also
$(\bm{\psi}_{I}-\mathbf{P}_{h}{\bm{\psi}})\in\mathbf{\widetilde{V}}_{h}^{\operatorname{face}}$;
see definition (47). We focus on the second term on the right-hand side:

$\begin{split}\|\bm{\psi}_{I}-&\mathbf{P}_{h}{\bm{\psi}}\|_{0,\Omega}^{2}\overset{\eqref{stability:face}}{\lesssim}\sum_{K\in\Omega_{h}}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\mu}^{-\frac{1}{2}}(\bm{\psi}_{I}-\mathbf{P}_{h}{\bm{\psi}})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{face},K}^{2}\\\
&\overset{\eqref{L2-projector-
face}}{=}\sum_{K\in\Omega_{h}}\\{[\widehat{\mu}^{-1}\bm{\psi}_{I},\bm{\psi}_{I}-\mathbf{P}_{h}{\bm{\psi}}]_{\operatorname{face},K}-(\mu^{-1}{\bm{\psi}},\bm{\psi}_{I}-\mathbf{P}_{h}{\bm{\psi}})_{0,K}\\}\\\
&\overset{\eqref{consistency:face}}{=}\sum_{K\in\Omega_{h}}\\{[\widehat{\mu}^{-1}(\bm{\psi}_{I}-\bm{\Pi}^{0}_{h}{\bm{\psi}}),\bm{\psi}_{I}-\mathbf{P}_{h}{\bm{\psi}}]_{\operatorname{face},K}-(\mu^{-1}({\bm{\psi}}-\bm{\Pi}^{0}_{h}{\bm{\psi}}),\bm{\psi}_{I}-\mathbf{P}_{h}{\bm{\psi}})_{0,K}\\\
&\quad\quad\quad\quad\quad+((\widehat{\mu}^{-1}-\mu^{-1})\bm{\Pi}^{0}_{h}{\bm{\psi}},\bm{\psi}_{I}-\mathbf{P}_{h}{\bm{\psi}})_{0,K}\\}\\\
&\overset{\eqref{stability:face}}{\lesssim}\left(\sum_{K\in\Omega_{h}}\\{\|{\bm{\psi}}-\bm{\Pi}^{0}_{h}{\bm{\psi}}\|_{0,K}^{2}+\|{\bm{\psi}}-\bm{\psi}_{I}\|_{0,K}^{2}+\|\widehat{\mu}^{-1}-\mu^{-1}\|_{L^{\infty}(K)}^{2}\|{\bm{\psi}}\|_{0,K}^{2}\\}\right)^{\frac{1}{2}}\|\bm{\psi}_{I}-\mathbf{P}_{h}{\bm{\psi}}\|_{0,\Omega}.\end{split}$

Using (29) and collecting the two above bounds, we get the assertion. ∎

###  3.7 The semi-discrete scheme

We denote the virtual element interpolant of the density current vector
$\mathbf{J}$ in $\mathbf{V}^{\operatorname{edge}}_{h}$ by $\mathbf{J}_{I}$. In
other words, $\mathbf{J}_{I}$ is the unique function in
$\mathbf{V}^{\operatorname{edge}}_{h}$ satisfying (20). Similarly, we define
the interpolants $\mathbf{E}_{0,I}$ and $\mathbf{B}_{0,I}$ of the initial
vector field $\mathbf{E}^{0}$ and $\mathbf{B}^{0}$ in
$\mathbf{V}^{\operatorname{edge}}_{h}$ and
$\mathbf{V}^{\operatorname{face}}_{h}$; cf. (20) and (28). The semi-discrete
scheme reads

$\begin{cases}\text{find
}(\mathbf{E}_{h},\mathbf{B}_{h})\in\mathbf{V}^{\operatorname{edge}}_{h}\times\mathbf{V}^{\operatorname{face}}_{h}\text{
such that }\\\\[5.0pt]
[\widehat{\varepsilon}\mathbf{E}_{h,t}+\widehat{\sigma}\mathbf{E}_{h},\mathbf{w}_{h}]_{\operatorname{edge}}-[\widehat{\mu}^{-1}\mathbf{B}_{h},\operatorname{\mathbf{curl}}\mathbf{w}_{h}]_{\operatorname{face}}=[\mathbf{J}_{I},\mathbf{w}_{h}]_{\operatorname{edge}}\quad\forall\mathbf{w}_{h}\in\mathbf{V}^{\operatorname{edge}}_{h}\\\\[5.0pt]
[\widehat{\mu}^{-1}\mathbf{B}_{h,t},\bm{\psi}_{h}]_{\operatorname{face}}+[\widehat{\mu}^{-1}\bm{\psi}_{h},\operatorname{\mathbf{curl}}\mathbf{E}_{h}]_{\operatorname{face}}=0\quad\forall\bm{\psi}_{h}\in\mathbf{V}^{\operatorname{face}}_{h},\end{cases}$
(50)

where we recall that the subscript $t$ stands for a derivative in time.

###  3.8 The fully-discrete scheme

We consider the fully-discrete approximation of (4) that is obtained by
applying the backward Euler time-stepping scheme to the semi-discrete problem
in (50). Higher order schemes in time can be built analogously. As customary,
we start by splitting the time integration interval $[0,T]$ in $N$ equally
spaced subintervals $[t^{n-1},t^{n}]$ with size $\tau=t^{n}-t^{n-1}$ for all
$n=1,\dots,N$. Moreover, let
$\mathbf{E}_{h}^{0}=\mathbf{E}_{\mathbf{I}}^{0}\in\mathbf{V}^{\operatorname{edge}}_{h}$
and
$\mathbf{B}_{h}^{0}=\mathbf{B}_{\mathbf{I}}^{0}\in\mathbf{V}^{\operatorname{face}}_{h}$
be the virtual element interpolations of $\mathbf{E}^{0}$ and
$\mathbf{B}^{0}$, cf. (20) and (28), respectively satisfying inequalities (15)
and (24).

Let $\mathbf{E}_{h}^{m}\in\mathbf{V}^{\operatorname{edge}}_{h}$ and
$\mathbf{B}_{h}^{m}\in\mathbf{V}^{\operatorname{face}}_{h}$ be the discrete
solutions at steps $m=0,\dots,N-1$. We compute the two discrete vector fields
at the time step $t^{m+1}$ using the implicit Euler scheme: _find
$\mathbf{E}_{h}^{m+1}\in\mathbf{V}^{\operatorname{edge}}_{h}$ and
$\mathbf{B}_{h}^{m+1}\in\mathbf{V}^{\operatorname{face}}_{h}$ such that, for
all $\mathbf{w}_{h}\in\mathbf{V}^{\operatorname{edge}}_{h}$ and
$\bm{\psi}_{h}\in\mathbf{V}^{\operatorname{face}}_{h}$,_

$\begin{cases}\frac{1}{\tau}[\widehat{\varepsilon}(\mathbf{E}_{h}^{m+1}-\mathbf{E}_{h}^{m}),\mathbf{w}_{h}]_{\operatorname{edge}}+[\widehat{\sigma}\mathbf{E}_{h}^{m+1},\mathbf{w}_{h}]_{\operatorname{edge}}-[\widehat{\mu}^{-1}\mathbf{B}_{h}^{m+1},\operatorname{\mathbf{curl}}\mathbf{w}_{h}]_{\operatorname{face}}=[\mathbf{J}_{I}^{m+1},\mathbf{w}_{h}]_{\operatorname{edge}}\\\\[10.00002pt]
\frac{1}{\tau}[\widehat{\mu}^{-1}(\mathbf{B}_{h}^{m+1}-\mathbf{B}_{h}^{m},\bm{\psi}_{h})]_{\operatorname{face}}+[\widehat{\mu}^{-1}\bm{\psi}_{h},\operatorname{\mathbf{curl}}\mathbf{E}_{h}^{m+1}]_{\operatorname{face}}=0.\end{cases}$
(51)

The existence and uniqueness of a solution for problem (51) follows with
standard arguments.

We can simplify (51) by rewriting the second equation as

$\displaystyle[\widehat{\mu}^{-1}(\mathbf{B}_{h}^{m+1}-\mathbf{B}_{h}^{m}+\tau\operatorname{\mathbf{curl}}\mathbf{E}_{h}^{m+1}),\bm{\psi}_{h}]_{\operatorname{face}}=0\quad\quad\forall\bm{\psi}_{h}\in\mathbf{V}^{\operatorname{face}}_{h}.$

Since
$\operatorname{\mathbf{curl}}(\mathbf{V}^{\operatorname{edge}}_{h})=\mathbf{\widetilde{V}}_{h}^{\operatorname{face}}$,
cf. (33), we deduce that

$\mathbf{B}_{h}^{m+1}=\mathbf{B}_{h}^{m}-\tau\operatorname{\mathbf{curl}}\mathbf{E}_{h}^{m+1}.$
(52)

Then, we substitute (52) in the first equation of (51) and find that

$\hskip
56.9055pt[(\widehat{\varepsilon}+\tau\,\widehat{\sigma})\mathbf{E}_{h}^{m+1},\mathbf{w}_{h}]_{\operatorname{edge}}+[\tau^{2}\widehat{\mu}^{-1}\operatorname{\mathbf{curl}}\mathbf{E}_{h}^{m+1},\operatorname{\mathbf{curl}}\mathbf{w}_{h}]_{\operatorname{face}}\\\\[5.0pt]
=[\tau\mathbf{J}_{I}^{m+1}+\widehat{\varepsilon}\mathbf{E}_{h}^{m},\mathbf{w}_{h}]_{\operatorname{edge}}+[\tau\,\widehat{\mu}^{-1}\mathbf{B}_{h}^{m},\operatorname{\mathbf{curl}}\mathbf{w}_{h}]_{\operatorname{face}}.\hskip
56.9055pt$

This reformulation allows us to reduce the computational effort in solving
(51).

In view of Remark 1 and assumption (2), we find that

$\operatorname{div}\mathbf{B}_{h}^{0}=\operatorname{div}\mathbf{B}_{\mathbf{I}}^{0}=(\operatorname{div}\mathbf{B}^{0})_{I}=0.$
(53)

We use (52) and (53), and apply the divergence operator to derive the discrete
counterpart of (3):

$\displaystyle\operatorname{div}\mathbf{B}_{h}^{m}=0\quad\quad\forall
m=0,\dots,M,$

which implies that our scheme provides an approximation of $\mathbf{B}$ that
naturally satisfies the divergence-free constraint.

###### Remark 2

The proposed scheme can be immediately extended to the case of general order
$k>1$ in space by substituting the above low order spaces
$\mathbf{V}^{\operatorname{edge}}_{h}$ and
$\mathbf{V}^{\operatorname{face}}_{h}$ with the corresponding ones from [10].
The theoretical analysis would follow along the same lines as that shown below
for the lowest order case. Yet, interpolation and stability properties in high
order edge and face virtual elements are work in progress.

##  4 Analysis of the semi-discrete scheme

In this section, we prove the convergence of the semi-discrete scheme (50).

###### Theorem 8

Let $(\mathbf{E},\mathbf{B})$ and $(\mathbf{E}_{h},\mathbf{B}_{h})$ be the
solutions to (4) and (50) under the geometric assumptions of Section 3.1 and
assumption (7). For all $t\in[0,T]$, we assume that $\mathbf{E}$,
$\mathbf{E}_{t}$ and $\mathbf{J}$ belong to
$L^{1}((0,T),H^{s}(\operatorname{\mathbf{curl}},\Omega))$, $1/2<s\leq 1$.
Furthermore, we recall that the initial vector-valued fields
$\mathbf{E}_{\mathbf{I}}^{0}\in\mathbf{V}^{\operatorname{edge}}_{h}$ and
$\mathbf{B}_{\mathbf{I}}^{0}\in\mathbf{V}^{\operatorname{face}}_{h}$
interpolate $\mathbf{E}^{0}$ and $\mathbf{B}^{0}$ in the sense of (20) and
(28), respectively. Then, the following error estimate is valid:

$\|\mathbf{E}(t)-\mathbf{E}_{h}(t)\|_{0,\Omega}+\|\mathbf{B}(t)-\mathbf{B}_{h}(t)\|_{0,\Omega}\lesssim
h^{s}\quad\quad\forall t\in[0,T].$ (54)

Proof. For all $t\in[0,T]$, we introduce

$\mathbf{e}_{h}=\mathbf{e}_{h}(t):=\bm{\Pi}_{h}^{\text{curl}}\mathbf{E}(t)-\mathbf{E}_{h}(t),\quad\mathbf{b}_{h}=\mathbf{b}_{h}(t):=\mathbf{P}_{h}\mathbf{B}(t)-\mathbf{B}_{h}(t).$
(55)

We recall that $\operatorname{div}\mathbf{B}_{\mathbf{I}}^{0}=0$. For all
$t\in[0,T]$, the definition of quantities (55), the projector $\mathbf{P}_{h}$
in (47), and the first equation in (50) allow us to prove that

$\begin{split}&[\widehat{\varepsilon}\mathbf{e}_{h,t}+\widehat{\sigma}\mathbf{e}_{h},\mathbf{w}_{h}]_{\operatorname{edge}}-[\widehat{\mu}^{-1}\mathbf{b}_{h},\operatorname{\mathbf{curl}}\mathbf{w}_{h}]_{\operatorname{face}}\\\\[5.0pt]
&=-[\mathbf{J}_{I},\mathbf{w}_{h}]_{\operatorname{edge}}+[\widehat{\varepsilon}\bm{\Pi}_{h}^{\text{curl}}\mathbf{E}_{t}+\widehat{\sigma}\bm{\Pi}_{h}^{\text{curl}}\mathbf{E},\mathbf{w}_{h}]_{\operatorname{edge}}-[\widehat{\mu}^{-1}\mathbf{P}_{h}\mathbf{B},\operatorname{\mathbf{curl}}\mathbf{w}_{h}]_{\operatorname{face}}\\\\[5.0pt]
&=-[\mathbf{J}_{I},\mathbf{w}_{h}]_{\operatorname{edge}}+[\widehat{\varepsilon}\bm{\Pi}_{h}^{\text{curl}}\mathbf{E}_{t}+\widehat{\sigma}\bm{\Pi}_{h}^{\text{curl}}\mathbf{E},\mathbf{w}_{h}]_{\operatorname{edge}}-(\mu^{-1}\mathbf{B},\operatorname{\mathbf{curl}}\mathbf{w}_{h})_{0,\Omega}\quad\mathbf{w}_{h}\in\mathbf{V}^{\operatorname{edge}}_{h}.\end{split}$
(56)

Moreover, for all
$\bm{\psi}_{h}\in\mathbf{\widetilde{V}}_{h}^{\operatorname{face}}$, we apply
(50), (47), the commuting property (48), and (4), and obtain

$\displaystyle\begin{split}&[\widehat{\mu}^{-1}\mathbf{b}_{h,t},\bm{\psi}_{h}]_{\operatorname{face}}+[\widehat{\mu}^{-1}\operatorname{\mathbf{curl}}\mathbf{e}_{h},\bm{\psi}_{h}]_{\operatorname{face}}=[\widehat{\mu}^{-1}\mathbf{P}_{h}\mathbf{B}_{t},\bm{\psi}_{h}]_{\operatorname{face}}+[\widehat{\mu}^{-1}\operatorname{\mathbf{curl}}\bm{\Pi}_{h}^{\text{curl}}(\mathbf{E}),\bm{\psi}_{h}]_{\operatorname{face}}\\\
&\hskip
28.45274pt=(\mu^{-1}\mathbf{B}_{t},\bm{\psi}_{h})_{0,\Omega}+[\widehat{\mu}^{-1}\mathbf{P}_{h}(\operatorname{\mathbf{curl}}\mathbf{E}),\bm{\psi}_{h}]_{\operatorname{face}}=(\mu^{-1}(\mathbf{B}_{t}+\operatorname{\mathbf{curl}}\mathbf{E}),\bm{\psi}_{h})_{0,\Omega}=0.\end{split}$

Since
$\mathbf{b}_{h,t}+\operatorname{\mathbf{curl}}\mathbf{e}_{h}\in\mathbf{\widetilde{V}}_{h}^{\operatorname{face}}$,
the above equation implies that

$\mathbf{b}_{h,t}+\operatorname{\mathbf{curl}}\mathbf{e}_{h}=\mathbf{0}.$ (57)

We set $\mathbf{w}_{h}=\mathbf{e}_{h}$ in (56), use (57), and deduce that

$\displaystyle\begin{split}&[\widehat{\varepsilon}\mathbf{e}_{h,t}+\widehat{\sigma}\mathbf{e}_{h},\mathbf{e}_{h}]_{\operatorname{edge}}+[\widehat{\mu}^{-1}\mathbf{b}_{h},\mathbf{b}_{h,t}]_{\operatorname{face}}\\\\[5.0pt]
&\hskip
56.9055pt=-[\mathbf{J}_{I},\mathbf{e}_{h}]_{\operatorname{edge}}+[\widehat{\varepsilon}\bm{\Pi}_{h}^{\text{curl}}\mathbf{E}_{t}+\widehat{\sigma}\bm{\Pi}_{h}^{\text{curl}}\mathbf{E},\mathbf{e}_{h}]_{\operatorname{edge}}-(\mu^{-1}\mathbf{B},\operatorname{\mathbf{curl}}\mathbf{e}_{h})_{0,\Omega}.\end{split}$

Next, we substitute
$(\mu^{-1}\mathbf{B},\operatorname{\mathbf{curl}}\mathbf{e}_{h})_{0,\Omega}$
with the expression given by the first equation (4):

$\begin{split}&\frac{1}{2}\partial_{t}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\mathbf{e}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}_{\operatorname{edge}}+[\widehat{\sigma}\mathbf{e}_{h},\mathbf{e}_{h}]_{\operatorname{edge}}+\frac{1}{2}\partial_{t}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\mu}^{-1/2}\mathbf{b}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}_{\operatorname{face}}\\\
&\quad=\frac{1}{2}\partial_{t}[\widehat{\varepsilon}\mathbf{e}_{h},\mathbf{e}_{h}]_{\operatorname{edge}}+[\widehat{\sigma}\mathbf{e}_{h},\mathbf{e}_{h}]_{\operatorname{edge}}+\frac{1}{2}\partial_{t}[\widehat{\mu}^{-1}\mathbf{b}_{h},\mathbf{b}_{h}]_{\operatorname{face}}\\\
&\quad=[-\mathbf{J}_{I},\mathbf{e}_{h}]_{\operatorname{edge}}-(\widehat{\mu}^{-1}\mathbf{B},\operatorname{\mathbf{curl}}\mathbf{e}_{h})_{0,\Omega}+[\widehat{\varepsilon}\bm{\Pi}_{h}^{\text{curl}}\mathbf{E}_{t}+\widehat{\sigma}\bm{\Pi}_{h}^{\text{curl}}\mathbf{E},\mathbf{e}_{h}]_{\operatorname{edge}}\\\
&\quad=\underbrace{(\mathbf{J},\mathbf{e}_{h})_{0,\Omega}-[\mathbf{J}_{I},\mathbf{e}_{h}]_{\operatorname{edge}}}_{=:T_{1}}+\underbrace{[\widehat{\varepsilon}\bm{\Pi}_{h}^{\text{curl}}\mathbf{E}_{t},\mathbf{e}_{h}]_{\operatorname{edge}}-(\varepsilon\mathbf{E}_{t},\mathbf{e}_{h})_{0,\Omega}}_{=:T_{2}}+\underbrace{[\widehat{\sigma}\bm{\Pi}_{h}^{\text{curl}}\mathbf{E},\mathbf{e}_{h}]_{\operatorname{edge}}-(\sigma\mathbf{E},\mathbf{e}_{h})_{0,\Omega}}_{=:T_{3}}.\end{split}$
(58)

We derive an upper bound for the three terms $T_{1}$, $T_{2}$ and $T_{3}$ on
the right-hand side of (58) separately. To estimate the term $T_{1}$, we use
the stability properties (13) and (15) of the bilinear form
$S^{K}_{\operatorname{face}}(\cdot,\cdot)$, employ standard polynomial
approximation results, use the interpolation property (21), and obtain

$\begin{split}T_{1}&=(\mathbf{J},\mathbf{e}_{h})_{0,\Omega}-(\bm{\Pi}^{0}_{h}\mathbf{J}_{I},\mathbf{e}_{h})_{0,\Omega}-\sum_{K\in\Omega_{h}}S^{K}_{\operatorname{edge}}((\mathbf{I}-\bm{\Pi}^{0}_{h})\mathbf{J}_{I},(\mathbf{I}-\bm{\Pi}^{0}_{h})\mathbf{e}_{h})\\\\[5.0pt]
&\lesssim(\|\mathbf{J}-\bm{\Pi}^{0}_{h}\mathbf{J}_{I}\|_{0,\Omega}+\|(\mathbf{I}-\bm{\Pi}^{0}_{h})\mathbf{J}_{I}\|_{0,\Omega}){\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathbf{e}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}\\\\[5.0pt]
&\lesssim(\|\mathbf{J}-\bm{\Pi}^{0}_{h}\mathbf{J}\|_{0,\Omega}+\|\mathbf{J}-\mathbf{J}_{I}\|_{0,\Omega}+\|\mathbf{J}-\bm{\Pi}^{0}_{h}\mathbf{J}_{I}\|_{0,\Omega}){\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathbf{e}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}\\\\[5.0pt]
&\lesssim(\|\mathbf{J}-\bm{\Pi}^{0}_{h}\mathbf{J}\|_{0,\Omega}+\|\mathbf{J}-\mathbf{J}_{I}\|_{0,\Omega}){\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathbf{e}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}\\\\[5.0pt]
&\lesssim
h^{s}(|\mathbf{J}|_{s,\Omega}+h^{1-s}\|\operatorname{\mathbf{curl}}\mathbf{J}\|_{0,\Omega}+h|\operatorname{\mathbf{curl}}\mathbf{J}|_{s,\Omega}){\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\mathbf{e}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}.\end{split}$
(59)

To estimate the term $T_{2}$, we introduce $\mathbf{c}$, the piecewise
constant average of $\mathbf{E}_{t}$ over $\Omega_{h}$, add and subtract
$(\widehat{\varepsilon}\mathbf{c},\mathbf{e}_{h})_{0,\Omega}=[\widehat{\varepsilon}\mathbf{c},\mathbf{e}_{h}]_{\operatorname{edge}}$,
note that $h\leq h^{s}$ for $s\leq 1$, and write

$\displaystyle\begin{split}T_{2}&\overset{\hskip
5.69054pt\eqref{consistency:edge}\hskip
8.53581pt}{=}[\widehat{\varepsilon}(\bm{\Pi}_{h}^{\text{curl}}\mathbf{E}_{t}-\mathbf{c}),\mathbf{e}_{h}]_{\operatorname{edge}}+(\widehat{\varepsilon}(\mathbf{c}-\mathbf{E}_{t}),\mathbf{e}_{h})_{0,\Omega}+((\widehat{\varepsilon}-\varepsilon)\mathbf{E}_{t},\mathbf{e}_{h})_{0,\Omega}\\\\[5.0pt]
&\overset{\eqref{assumption:coefficients},\eqref{stability:edge}}{\lesssim}\left(\|\mathbf{E}_{t}-\bm{\Pi}_{h}^{\text{curl}}\mathbf{E}_{t}\|_{0,\Omega}+\|\mathbf{E}_{t}-\mathbf{c}\|_{0,\Omega}+h\max_{K\in\Omega_{h}}|\varepsilon|_{W^{1,\infty}(K)}\|\mathbf{E}_{t}\|_{0,\Omega}\right){\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathbf{e}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}\\\\[5.0pt]
&\overset{\hskip 5.69054pt\eqref{approximation:Piboldh}\hskip
8.53581pt}{\lesssim}h^{s}\left(\|\mathbf{E}_{t}\|_{s,\Omega}+h^{1-s}\|\operatorname{\mathbf{curl}}\mathbf{E}_{t}\|_{0,\Omega}+h|\operatorname{\mathbf{curl}}\mathbf{E}_{t}|_{s,\Omega}\right){\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\mathbf{e}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}.\end{split}$

Recalling assumption (7) again, we treat the term $T_{3}$ analogously and
arrive at the bound

$T_{3}\lesssim
h^{s}\left(\|\mathbf{E}\|_{s,\Omega}+h^{1-s}\|\operatorname{\mathbf{curl}}\mathbf{E}\|_{0,\Omega}+h(|\operatorname{\mathbf{curl}}\mathbf{E}|_{s,\Omega}\right){\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\mathbf{e}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}.$
(60)

Introduce the regularity type term, which belong to $L^{1}(0,T)$ due to the
assumptions of the theorem,

$\displaystyle\begin{split}c_{\textbf{REG}}(t)=c_{\textbf{REG}}&=|\mathbf{J}|_{s,\Omega}+\|\mathbf{E}\|_{s,\Omega}+\|\mathbf{E}_{t}\|_{s,\Omega}+h^{1-s}\|\operatorname{\mathbf{curl}}\mathbf{J}\|_{0,\Omega}+h|\operatorname{\mathbf{curl}}\mathbf{J}|_{s,\Omega}\\\
&\quad+h^{1-s}(\|\operatorname{\mathbf{curl}}\mathbf{E}\|_{0,\Omega}+\|\operatorname{\mathbf{curl}}\mathbf{E}_{t}\|_{0,\Omega})+h(|\operatorname{\mathbf{curl}}\mathbf{E}|_{s,\Omega}).\end{split}$

Now, we collect the upper bounds on the terms $T_{1}$, $T_{2}$, and $T_{3}$ in
(58), recall (5), and deduce that

$\frac{1}{2}\partial_{t}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\mathbf{e}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}^{2}\leq\frac{1}{2}\left(\partial_{t}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\mathbf{e}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}^{2}+\partial_{t}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\mu}^{-1/2}\mathbf{b}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{face}}^{2}\right)\lesssim
c_{\textbf{REG}}h^{s}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\mathbf{e}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}.$
(61)

The following identity is valid:

$\displaystyle\frac{1}{2}\partial_{t}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\,\cdot\,\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}^{2}={\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\,\cdot\,\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}\,\partial_{t}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\,\cdot\,\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}.$

We use this identity in (61), so that, almost everywhere in time in $(0,T)$,

$\displaystyle\partial_{t}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\mathbf{e}_{h}(t)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}\lesssim
c_{\textbf{REG}}h^{s},$

and we integrate in time to obtain:

$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\mathbf{e}_{h}(t)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}\lesssim{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\mathbf{e}_{h}(0)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}+h^{s}\int_{0}^{t}c_{\textbf{REG}}(s){\rm
d}s\lesssim{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\mathbf{e}_{h}(0)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}+h^{s}.$
(62)

The error at the initial time $t=0$ is controlled as follows:

$\begin{split}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\mathbf{e}_{h}(0)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}&\overset{\hskip
8.53581pt\eqref{stability:edge},\eqref{eq:param:bound}\hskip
8.53581pt}{\lesssim}\|\mathbf{E}^{0}-\mathbf{E}_{\mathbf{I}}^{0}\|_{0,\Omega}+\|\mathbf{E}^{0}-\bm{\Pi}_{h}^{\text{curl}}\mathbf{E}^{0}\|_{0,\Omega}\\\\[5.0pt]
&\overset{\eqref{interpolation-
estimates:edge},\eqref{approximation:Piboldh}}{\lesssim}h^{s}\left(|\mathbf{E}^{0}|_{s,\Omega}+h^{1-s}\|\operatorname{\mathbf{curl}}(\mathbf{E}^{0})\|_{0,\Omega}+h|\operatorname{\mathbf{curl}}(\mathbf{E}^{0})|_{s,\Omega}\right)\lesssim
h^{s},\end{split}$ (63)

and using this inequality in (62) yields

$\|\mathbf{e}_{h}(t)\|_{0,\Omega}\overset{\eqref{stability:edge},\eqref{eq:param:bound}}{\lesssim}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\mathbf{e}_{h}(t)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}\lesssim
h^{s}.$ (64)

Thus, we can write

$\partial_{t}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\mu}^{-1/2}\mathbf{b}_{h}(t)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}_{\operatorname{face}}\overset{\eqref{bound:in-
time}}{\lesssim}c_{\textbf{REG}}h^{s}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\mathbf{e}_{h}(t)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}\overset{\eqref{star1}}{\lesssim}c_{\textbf{REG}}h^{2s}.$
(65)

Integrating in time (65) gives

$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\mu}^{-1/2}\mathbf{b}_{h}(t)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}_{\operatorname{face}}\lesssim{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\mu}^{-1/2}\mathbf{b}_{h}(0)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}_{\operatorname{face}}+h^{2s}.$

Observe that

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\mu}^{-1/2}\mathbf{b}_{h}(0)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{face}}\overset{\eqref{stability:face},\eqref{eq:param:bound}}{\lesssim}\|\mathbf{B}^{0}-\mathbf{P}_{h}\mathbf{B}^{0}\|_{0,\Omega}+\|\mathbf{B}^{0}-\mathbf{B}_{\mathbf{I}}^{0}\|_{0,\Omega}\overset{\eqref{interpolation-
estimates:face},\eqref{property-
Pboldh}}{\lesssim}h\|\mathbf{B}^{0}\|_{0,\Omega}+h^{s}|\mathbf{B}^{0}|_{s,\Omega}\lesssim
h^{s}.$ (66)

Then, we have

$\|\mathbf{b}_{h}(t)\|_{0,\Omega}\overset{\eqref{stability:face},\eqref{eq:param:bound}}{\lesssim}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\mu}^{-1/2}\mathbf{b}_{h}(t)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{face}}\lesssim
h^{s}.$ (67)

Finally, we add and subtract $\bm{\Pi}_{h}^{\text{curl}}(\mathbf{E}(t))$,
$\mathbf{P}_{h}(\mathbf{B}(t))$, and use the definitions of
$\mathbf{e}_{h}(t)$ and $\mathbf{b}_{h}(t)$ and the triangle inequality to
obtain

$\displaystyle\begin{split}&\|\mathbf{E}(t)-\mathbf{E}_{h}(t)\|_{0,\Omega}+\|\mathbf{B}(t)-\mathbf{B}_{h}(t)\|_{0,\Omega}\\\
&\lesssim\|\mathbf{E}_{h}(t)-\bm{\Pi}_{h}^{\text{curl}}(\mathbf{E}(t))\|_{0,\Omega}+\|\mathbf{B}_{h}(t)-\mathbf{P}_{h}(\mathbf{B}(t))\|_{0,\Omega}+\|\mathbf{e}_{h}(t)\|_{0,\Omega}+\|\mathbf{b}_{h}(t)\|_{0,\Omega}.\end{split}$

The assertion of the theorem follows on from using (44), (49), (64), and (67).
∎

##  5 Analysis of the fully-discrete scheme

In this section, we prove the convergence of the fully-discrete scheme (51).
Notably, we recall that we employ the implicit Euler scheme for the time
discretization and subdivide the time interval $[0,T]$ into $M$ sub-interval
of uniform length $\tau$. We can extend the result below to other, possibly
higher order, time discretization schemes.

###### Theorem 9

Let the geometric assumptions of Section 3.1 and assumption (7) be valid, and
$(\mathbf{E},\mathbf{B})$ be the solutions to Maxwell’s equations (4). We
assume that $\mathbf{E}(t)$ and $\mathbf{E}_{t}(t)$ belong to
$L^{\infty}((0,T),H^{s}(\operatorname{\mathbf{curl}},\Omega))$, and
$\partial_{tt}\mathbf{E}$ and $\partial_{tt}\mathbf{B}$ to
$L^{\infty}((0,T),[L^{2}(\Omega)]^{3})$, $1/2<s\leq 1$. Additionally, we
assume that the electric current density $\mathbf{J}(t)$ in the right-hand
side of (4) belongs to
$L^{\infty}((0,T),H^{s}(\operatorname{\mathbf{curl}},\Omega))$ for the same
value of $s$ as above. For all $m=0,\dots,M$, let
$(\mathbf{E}_{h}^{m},\mathbf{B}_{h}^{m})$ denote the solutions of the fully-
discrete scheme (51) at the time step $t^{m}$. Then, for sufficiently small
$\tau$ as required in (78), the following error estimate is valid:

$\begin{split}&\|\mathbf{E}(t^{m})-\mathbf{E}_{h}^{m}\|_{0,\Omega}+\|\mathbf{B}(t^{m})-\mathbf{B}_{h}^{m}\|_{0,\Omega}\lesssim
h^{s}+\tau\quad\quad\forall m=0,\dots,M.\end{split}$ (68)

Proof. Let $\bm{\Pi}_{h}^{\text{curl}}$ and $\mathbf{P}_{h}$ be the two
projectors introduced in (34) and (47), whose approximation properties are
detailed in (44) and (49), respectively, and introduce

$\mathbf{e}_{h}^{m}:=\bm{\Pi}_{h}^{\text{curl}}(\mathbf{E}(t^{m}))-\mathbf{E}_{h}^{m},\quad\mathbf{b}_{h}^{m}:=\mathbf{P}_{h}(\mathbf{B}(t^{m}))-\mathbf{B}_{h}^{m}\quad\quad\forall
m=0,\dots,M.$ (69)

As a first step, we show two error equations, which we can deduce from
definition (69) and the fully-discrete problem (51). The first error equation
reads as: for all
$\bm{\psi}_{h}\in\mathbf{\widetilde{V}}_{h}^{\operatorname{face}}$,

$\begin{split}&\frac{1}{\tau}[\widehat{\mu}^{-1}(\mathbf{b}_{h}^{m+1}-\mathbf{b}_{h}^{m}),\bm{\psi}_{h}]_{\operatorname{face}}+[\widehat{\mu}^{-1}\operatorname{\mathbf{curl}}(\mathbf{e}_{h}^{m+1}),\bm{\psi}_{h}]_{\operatorname{face}}\\\
&\overset{\eqref{VEM:Maxwell}}{=}\frac{1}{\tau}[\widehat{\mu}^{-1}\mathbf{P}_{h}(\mathbf{B}(t^{m+1})-\mathbf{B}(t^{m})),\bm{\psi}_{h}]_{\operatorname{face}}+[\widehat{\mu}^{-1}\operatorname{\mathbf{curl}}(\bm{\Pi}_{h}^{\text{curl}}\mathbf{E}(t^{m+1})),\bm{\psi}_{h}]_{\operatorname{face}}\\\
&\overset{\eqref{L2-projector-face},\eqref{commuting-
diagram}}{=}\frac{1}{\tau}(\mu^{-1}(\mathbf{B}(t^{m+1})-\mathbf{B}(t^{m})),\bm{\psi}_{h})_{0,\Omega}+[\widehat{\mu}^{-1}\mathbf{P}_{h}(\operatorname{\mathbf{curl}}(\mathbf{E}(t^{m+1}))),\bm{\psi}_{h}]_{\operatorname{face}}\\\
&\overset{\eqref{L2-projector-
face}}{=}(\mu^{-1}((\mathbf{B}(t^{m+1})-\mathbf{B}(t^{m}))/\tau+\operatorname{\mathbf{curl}}(\mathbf{E}(t^{m+1}))),\bm{\psi}_{h})_{0,\Omega}\\\
&\overset{\eqref{eq:Maxwell:strong}}{=}(\mu^{-1}((\mathbf{B}(t^{m+1})-\mathbf{B}(t^{m}))/\tau-\partial_{t}\mathbf{B}(t^{m+1})),\bm{\psi}_{h})_{0,\Omega}=:(\omega^{m},\bm{\psi}_{h})_{0,\Omega},\\\
\end{split}$ (70)

which intrinsically defines the last term $\omega^{m}$. The second error
equation reads

$\begin{split}&\frac{1}{\tau}[\widehat{\varepsilon}(\mathbf{e}_{h}^{m+1}-\mathbf{e}_{h}^{m}),\mathbf{e}_{h}^{m+1}]_{\operatorname{edge}}+[\widehat{\sigma}\mathbf{e}_{h}^{m+1},\mathbf{e}_{h}^{m+1}]_{\operatorname{edge}}-[\widehat{\mu}^{-1}\mathbf{b}_{h}^{m+1},\operatorname{\mathbf{curl}}(\mathbf{e}_{h}^{m+1})]_{\operatorname{face}}\\\
&\overset{\eqref{eq:Maxwell:weak},\eqref{L2-projector-
face}}{=}[-\mathbf{J}_{I},\mathbf{e}_{h}^{m+1}]_{\operatorname{edge}}+[\widehat{\varepsilon}(\bm{\Pi}_{h}^{\text{curl}}(\mathbf{E}(t^{m+1})-\mathbf{E}(t^{m})))/\tau+\widehat{\sigma}\bm{\Pi}_{h}^{\text{curl}}(\mathbf{E}(t^{m+1})),\mathbf{e}_{h}^{m+1}]_{\operatorname{edge}}\\\
&\quad\quad\quad-(\mu^{-1}\mathbf{B}(t^{m+1}),\operatorname{\mathbf{curl}}(\mathbf{e}_{h}^{m+1}))_{0,\Omega}.\end{split}$
(71)

We pick $\bm{\psi}_{h}=\mathbf{b}_{h}^{m+1}$ as a test function in (70), sum
the resulting equation with (71), and get

$\displaystyle\begin{split}&\frac{1}{\tau}[\widehat{\varepsilon}(\mathbf{e}_{h}^{m+1}-\mathbf{e}_{h}^{m}),\mathbf{e}_{h}^{m+1}]_{\operatorname{edge}}+[\widehat{\sigma}\mathbf{e}_{h}^{m+1},\mathbf{e}_{h}^{m+1}]_{\operatorname{edge}}+\frac{1}{\tau}[\widehat{\mu}^{-1}(\mathbf{b}_{h}^{m+1}-\mathbf{b}_{h}^{m}),\mathbf{b}_{h}^{m+1}]_{\operatorname{face}}\\\
&=[-\mathbf{J}_{I},\mathbf{e}_{h}^{m+1}]_{\operatorname{edge}}+[\widehat{\varepsilon}(\bm{\Pi}_{h}^{\text{curl}}(\mathbf{E}(t^{m+1})-\mathbf{E}(t^{m})))/\tau+\widehat{\sigma}\bm{\Pi}_{h}^{\text{curl}}(\mathbf{E}(t^{m+1})),\mathbf{e}_{h}^{m+1}]_{\operatorname{edge}}\\\
&\quad\quad\quad-(\mu^{-1}\mathbf{B}(t^{m+1}),\operatorname{\mathbf{curl}}(\mathbf{e}_{h}^{m+1}))_{0,\Omega}+(\omega^{m},\mathbf{b}_{h}^{m+1})_{0,\Omega}=:T_{1}+T_{2}+T_{3}+T_{4},\end{split}$

where, using the first equation of (4) with $\mathbf{e}_{h}^{m+1}$ as a test
function, we have set

$\displaystyle
T_{1}:=(\mathbf{J},\mathbf{e}_{h}^{m+1})_{0,\Omega}-[\mathbf{J}_{I},\mathbf{e}_{h}^{m+1}]_{\operatorname{edge}},$
(72) $\displaystyle
T_{2}:=[\widehat{\varepsilon}(\bm{\Pi}_{h}^{\text{curl}}(\mathbf{E}(t^{m+1})-\mathbf{E}(t^{m})))/\tau,\mathbf{e}_{h}^{m+1}]_{\operatorname{edge}}-(\varepsilon\partial_{t}\mathbf{E}(t^{m+1}),\mathbf{e}_{h}^{m+1})_{0,\Omega},$
(73) $\displaystyle
T_{3}:=[\widehat{\sigma}\bm{\Pi}_{h}^{\text{curl}}(\mathbf{E}(t^{m+1})),\mathbf{e}_{h}^{m+1}]_{\operatorname{edge}}-(\sigma\mathbf{E}(t^{m+1}),\mathbf{e}_{h}^{m+1})_{0,\Omega},$
(74) $\displaystyle
T_{4}:=(\mu^{-1}(\mathbf{B}(t^{m+1})-\mathbf{B}(t^{m}))/\tau-\partial_{t}\mathbf{B}(t^{m+1}),\mathbf{b}_{h}^{m+1})_{0,\Omega}.$
(75)

We easily deduce that

$\begin{split}&{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\mathbf{e}_{h}^{m+1}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}^{2}+{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\mu}^{1/2}\mathbf{b}_{h}^{m+1}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{face}}^{2}\\\
&\lesssim{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\mathbf{e}_{h}^{m}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\mathbf{e}_{h}^{m+1}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}+{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\mu}^{1/2}\mathbf{b}_{h}^{m}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{face}}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\mu}^{1/2}\mathbf{b}_{h}^{m+1}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{face}}+\tau(T_{1}+T_{2}+T_{3}+T_{4}).\end{split}$
(76)

For the time being, assume the following bound is valid:

$\tau(T_{1}+T_{2}+T_{3}+T_{4})\lesssim\tau(h^{s}+\tau){\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\mathbf{e}_{h}^{m+1}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}\lesssim\tau{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\mathbf{e}_{h}^{m+1}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}_{\operatorname{edge}}+\tau(h^{s}+\tau)^{2}.$
(77)

We shall show (77) at the end of the proof. Inserting (77) in (76) and some
standard manipulations yield, for a positive $c$ independent of $h$ and
$\tau$,

$\displaystyle\begin{split}&{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\mathbf{e}_{h}^{m+1}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}^{2}+{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\mu}^{1/2}\mathbf{b}_{h}^{m+1}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{face}}^{2}\\\
&\leq
c\left({\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\mathbf{e}_{h}^{m}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}^{2}+{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\mu}^{1/2}\mathbf{b}_{h}^{m}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{face}}^{2}+\tau{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\mathbf{e}_{h}^{m+1}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}_{\operatorname{edge}}+\tau(h^{s}+\tau)^{2}\right).\end{split}$

In other words, for $\tau\leq 1/c$, we have

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\mathbf{e}_{h}^{m+1}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}^{2}+{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\mu}^{1/2}\mathbf{b}_{h}^{m+1}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{face}}^{2}\leq\frac{1}{1-c\tau}({\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\mathbf{e}_{h}^{m}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}^{2}+{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\mu}^{1/2}\mathbf{b}_{h}^{m}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{face}}^{2})+\frac{c\tau}{1-c\tau}(h^{s}+\tau)^{2}.$
(78)

Bound (78) has the form

$\displaystyle
a_{m+1}\leq\frac{1}{1-c\tau}a_{m}+\frac{c\tau}{1-c\tau}(h^{s}+\tau)^{2}.$

Recalling that $M:=T/\tau$, we can iterate the above bound and write

$\displaystyle\begin{split}a_{m}&\leq\left(\frac{1}{1-c\tau}\right)^{M}a_{0}+c\tau\left(\sum_{i=1}^{M}\left(\frac{1}{1-c\tau}\right)^{i}\right)(h^{s}+\tau)^{2}\\\
&\leq\left(\frac{1}{1-c\tau}\right)^{M}a_{0}+cT\left(\frac{1}{1-c\tau}\right)^{M}(h^{s}+\tau)^{2}=\left(\frac{1}{1-c\tau}\right)^{M}\left(a_{0}+cT(h^{s}+\tau)^{2}\right)\quad\quad\forall
m=1,\dots,M.\end{split}$

By noting that

$\displaystyle\left(\frac{1}{1-c\tau}\right)^{M}=\left(\frac{1}{1-c\tau}\right)^{\frac{T}{\tau}}$

is uniformly bounded as $\tau\rightarrow 0$, we achieve

$\displaystyle a_{m}\lesssim a_{0}+(h^{s}+\tau)^{2}\quad\quad\forall
m=1,\dots,M.$

Thus, bound (78) yields

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\mathbf{e}_{h}^{m+1}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}^{2}+{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\mu}^{1/2}\mathbf{b}_{h}^{m+1}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{face}}^{2}\lesssim{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\mathbf{e}_{h}^{0}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}^{2}+{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\mu}^{1/2}\mathbf{b}_{h}^{0}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{face}}^{2}+(h^{s}+\tau)^{2}.$
(79)

The assertion of the theorem follows from the triangle inequality, the data
assumptions (5), stability properties (15) and (24), and bounds (63) and (66)
on the initial data approximation error.

Thus, we are left to show (77), i.e., the upper bounds on the terms $T_{i}$,
$i=1,\dots,4$, in (72)-(75). In the following bounds, we shall use the data
assumption (5) several times. Therefore, this will not be declared at every
instance. We can deal with the terms $T_{1}$ and $T_{3}$ as in the semi-
discrete analysis of Theorem 8. More precisely, proceeding as in (59) and
(60), we can write

$\displaystyle T_{1}\lesssim
h^{s}(|\mathbf{J}|_{s,\Omega}+h^{1-s}\|\operatorname{\mathbf{curl}}\mathbf{J}\|_{0,\Omega}+h|\operatorname{\mathbf{curl}}\mathbf{J}|_{s,\Omega}){\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathbf{e}_{h}^{m+1}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}\lesssim
h^{s}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\mathbf{e}_{h}^{m+1}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}$

and

$\displaystyle\begin{split}T_{3}&\overset{\eqref{assumption:coefficients}}{\lesssim}h^{s}(|\mathbf{E}(t^{m+1})|_{s,\Omega}+h^{1-s}\|\operatorname{\mathbf{curl}}(\mathbf{E}(t^{m+1}))\|_{0,\Omega}+h\max_{K\in\Omega_{h}}|\varepsilon|_{W^{1,\infty}(K)}|\operatorname{\mathbf{curl}}(\mathbf{E}(t^{m+1}))|_{s,\Omega}){\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathbf{e}_{h}^{m+1}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}\\\
&\lesssim
h^{s}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\mathbf{e}_{h}^{m+1}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}.\end{split}$

Next, we focus on the two remaining terms and start with $T_{4}$:

$\displaystyle\begin{split}T_{4}&:=(\mu^{-1}\partial_{t}\mathbf{B}(t^{m+1})-(\mathbf{B}(t^{m+1})-\mathbf{B}(t^{m}))/\tau,\mathbf{e}_{h}^{m+1})_{0,\Omega}\\\
&\overset{\eqref{stability:edge}}{\lesssim}\|\partial_{t}\mathbf{B}(t^{m+1})-(\mathbf{B}(t^{m+1})-\mathbf{B}(t^{m}))/\tau\|_{0,\Omega}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathbf{e}_{h}^{m+1}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}\\\
&=\|\frac{1}{\tau}\int_{t^{m}}^{t^{m+1}}(s-t^{m})\partial_{tt}\mathbf{B}(s)ds\|_{0,\Omega}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathbf{e}_{h}^{m+1}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}\\\
&\leq\tau\|\partial_{tt}\mathbf{B}\|_{L^{\infty}([t^{m},t^{m+1}],L^{2}(\Omega))}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathbf{e}_{h}^{m+1}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}\lesssim\tau{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\mathbf{e}_{h}^{m+1}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}.\\\
\end{split}$

As for the term $T_{2}$, we consider the splitting

$\displaystyle\begin{split}T_{2}&=[\widehat{\varepsilon}(\bm{\Pi}_{h}^{\text{curl}}(\mathbf{E}(t^{m+1})-\mathbf{E}(t^{m})))/\tau,\mathbf{e}_{h}^{m+1}]_{\operatorname{edge}}-(\varepsilon(\mathbf{E}(t^{m+1})-\mathbf{E}(t^{m}))/\tau,\mathbf{e}_{h}^{m+1})_{0,\Omega}\\\
&\quad+(\varepsilon(\mathbf{E}(t^{m+1})-\mathbf{E}(t^{m}))/\tau,\mathbf{e}_{h}^{m+1})_{0,\Omega}-(\varepsilon\partial_{t}\mathbf{E}(t^{m+1}),\mathbf{e}_{h}^{m+1})_{0,\Omega}=:T_{2,1}+T_{2,2}.\end{split}$

The term $T_{2,2}$ is dealt with as the term $T_{4}$:

$\displaystyle
T_{2,2}\leq\tau\|\partial_{tt}\mathbf{E}\|_{L^{\infty}([t^{m},t^{m+1}],L^{2}(\Omega))}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathbf{e}_{h}^{m+1}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}\lesssim\tau{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\mathbf{e}_{h}^{m+1}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}}.$

Finally, we show an upper bound on the term $T_{2,1}$ proceeding as for the
term $T_{3}$:

$\displaystyle T_{2,1}\lesssim
h^{s}\|(\mathbf{E}(t^{m+1})-\mathbf{E}(t^{m}))/\tau\|_{*,\Omega}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\varepsilon}^{1/2}\mathbf{e}_{h}^{m+1}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\operatorname{edge}},$

where we have used (7) again and set

$\displaystyle\|\cdot\|_{*,\Omega}:=|\cdot|_{s,\Omega}+h^{1-s}\|\operatorname{\mathbf{curl}}(\cdot)\|_{0,\Omega}+h|\operatorname{\mathbf{curl}}(\cdot)|_{s,\Omega}.$

We must prove that the *-norm of the difference quotient is finite. To this
aim, observe

$\displaystyle\begin{split}&\|(\mathbf{E}(t^{m+1})-\mathbf{E}(t^{m})/\tau)\|_{*,\Omega}=\left\|1/\tau\int_{t^{m}}^{t^{m+1}}\partial_{t}\mathbf{E}(s)ds\right\|_{*,\Omega}\\\
&\leq
1/\tau\int_{t^{m}}^{t^{m+1}}\|\partial_{t}\mathbf{E}(s)\|_{*,\Omega}ds\leq\|\partial_{t}\mathbf{E}\|_{L^{\infty}((0,T),H^{s}(\operatorname{\mathbf{curl}},\Omega))}.\end{split}$

Collecting the bounds on the terms $T_{1}$, $T_{2,1}$, $T_{2,2}$, $T_{3}$, and
$T_{4}$, we deduce (77), whence the assertion follows. ∎

##  6 Numerical results

In this section, we investigate the accuracy of the fully discrete scheme
(51). To this end, we consider three different mesh families:

  * •

cube: regular cubic meshes;

  * •

voro: Voronoi tessellations optimized by the Lloyd algorithm;

  * •

rand: Voronoi tessellations of a cloud of points that are randomly positioned
in the computational domain.

We selected these three types of meshes as they offer an increasing level of
geometric difficulty. Indeed, the meshes in cube are uniform; the meshes in
voro may have small edges and faces but the geometric shape of the mesh
elements is not distorted; finally meshes rand may have small edges and faces,
as well as stretched polyhedral elements. We refer to a specific partition of
$\Omega$ by the corresponding keyword (cube, voro, and rand) followed by the
number of elements. For example, “voro125” refers to a mesh made of $125$
Voronoi cells optimized by the Lloyd algorithm.

We numerically verify the optimal convergence rate in the $L^{2}$ norm of the
approximation to the electric field $\mathbf{E}$ and the magnetic flux field
$\mathbf{B}$ on a sequence of four refined meshes for each mesh family. These
four meshes have a decreasing mesh size. We show the third mesh of each family
in Figure 1.

cube1000 | voro1000 | rand1000
---|---|---
|  mm |  mm
Figure 1: The third mesh of the mesh families cube (left panel), voro (middle
panel), and rand (right panel).

The virtual element approximations $\mathbf{E}_{h}$ and $\mathbf{B}_{h}$ to
$\mathbf{E}$ and $\mathbf{B}$ are not available in closed form so we evaluate
the error in the $L^{2}$ norm at any time $T$ by using the polynomial
projections $\bm{\Pi}^{0}_{h}\mathbf{E}_{h}$ and
$\bm{\Pi}^{0}_{h}\mathbf{B}_{h}$.

We performed a sensitivity analysis on the stabilizations terms that appear in
the definition of edge and face discrete scalar products. More precisely, we
inserted two constant coefficients, $\eta_{\text{edge}}^{K}$ and
$\eta_{\text{face}}^{K}$ in front of
$S^{K}_{\operatorname{edge}}(\cdot,\cdot)$ and
$S^{K}_{\operatorname{face}}(\cdot,\cdot)$, respectively. We observe that the
best choice is given by $\eta_{\text{edge}}^{K}=0.01$ and
$\eta_{\text{face}}^{K}=0.5$ and use these values in all numerical tests.

###  6.1 Test Case 1: constant coefficients

We solve Maxwell’s equations on the computational domain $\Omega=(0,\,1)^{3}$
for $t\in[0,1]$ with constant coefficients $\varepsilon,\,\sigma$, and $\mu$
equal to $1$. The boundary condition and the current density vector are
computed by taking the exact solution fields

$\mathbf{E}(t,x,y,z)=t\operatorname{\mathbf{curl}}{\bm{\phi}}(x,y,z)+t^{2}{\bm{\psi}}(x,y,z)\quad\text{and}\quad\mathbf{B}(t,x,y,z)=\frac{t^{2}}{2}\operatorname{\mathbf{curl}}\operatorname{\mathbf{curl}}{\bm{\phi}}(x,y,z)\,,$

where the auxiliary vector-valued fields ${\bm{\phi}}$ and ${\bm{\psi}}$ are
defined as

${\bm{\phi}}=\left(\begin{array}[]{c}\sin^{2}(\pi
x)\,y^{2}(1-y)^{2}\,z^{2}(1-z)^{2}\\\\[5.0pt] x^{2}(1-x)^{2}\,\sin^{2}(\pi
y)\,z^{2}(1-z)^{2}\\\\[5.0pt] x^{2}(1-x)^{2}\,y^{2}(1-y)^{2}\,\sin^{2}(\pi
z)\end{array}\right),\quad\quad{\bm{\psi}}=\nabla\big{(}\sin(\pi x))\sin(\pi
y)\sin(\pi z)\big{)}.$

A straightforward calculation shows that $\operatorname{div}\mathbf{B}=0$,
i.e., the magnetic field $\mathbf{B}$ is solenoidal.

In Figure 2, we plot the $L^{2}$ errors at final time $T=1$ for simultaneous
refinements of $h$ and $\tau$ on the three mesh families and observe the
expected convergence rate, which we recall has to be proportional to $h+\tau$;
see Theorem 9.

|
---|---
Figure 2: Test Case 1: Error curves in the $L^{2}$ norm at final time $T=1$
for the virtual element approximation $\mathbf{E}_{h}$ (left panel) and
$\mathbf{B}_{h}$ (right panel) using mesh families cube, voro, and rand for
simultaneous refinements of $h$ and $\tau$.

In Tables 1 and 2 we report the approximation errors for rand meshes only
since we observe similar behavior of both cube and voro families.

Each column of the tables shows how the method converges with respect to space
discretization, i.e., by using a constant time step and refining the mesh.
Likewise, each row of the tables shows how the method converges with respect
to the time discretization, i.e., by halving the time step on a fixed mesh.
The error in space seems to be the dominant effect so it hides the convergence
in time. Indeed, the errors does not halve along rows while they do halve
along colums. However, the errors along the diagonal show how VEM behaves when
we simultaneously refine the numerical calculations in space and time.

$h/\tau$ | 1 | 1/2 | 1/4 | 1/8 | 1/16
---|---|---|---|---|---
rand27 | 7.59493e-01 | 6.87305e-01 | 6.61804e-01 | 6.53871e-01 | 6.51503e-01
rand125 | 6.75411e-01 | 5.28227e-01 | 4.65715e-01 | 4.42600e-01 | 4.34238e-01
rand1000 | 5.57938e-01 | 3.71547e-01 | 2.84666e-01 | 2.52453e-01 | 2.42061e-01
rand8000 | 5.15998e-01 | 3.05485e-01 | 1.93696e-01 | 1.45600e-01 | 1.29105e-01
Table 1: Test Case 1: $L^{2}$-norms at final time $T=1$ of the error for the virtual element approximation $\mathbf{E}_{h}$ using mesh family rand for various combinations of $h$ and $\tau$. $h/\tau$ | 1 | 1/2 | 1/4 | 1/8 | 1/16
---|---|---|---|---|---
rand27 | 1.44376e+00 | 1.41488e+00 | 1.40800e+00 | 1.41438e+00 | 1.42469e+00
rand125 | 8.04334e-01 | 7.98555e-01 | 7.96968e-01 | 7.96871e-01 | 7.96838e-01
rand1000 | 4.17143e-01 | 4.13694e-01 | 4.12693e-01 | 4.12505e-01 | 4.12469e-01
rand8000 | 2.15819e-01 | 2.12642e-01 | 2.11953e-01 | 2.11841e-01 | 2.11817e-01
Table 2: Test Case 1: $L^{2}$-norms of the error at final time $T=1$ for the
virtual element approximation $\mathbf{B}_{h}$ using mesh family rand for
various combinations of $h$ and $\tau$.

Finally, in Table 3 we report the $L^{2}$-norm of the divergence of
$\mathbf{B}_{h}$ for each combination of $h$ and $\tau$. This table confirms
that the numerical approximation to the magnetic field provided by the VEM is
divergence free. Indeed, all the values of the divergence are very small even
if a slight growth is visible during the refinement process for $h\to 0$,
which is very likely due to round-off effects related to the conditioning of
the final linear system.

Such interpretation is also supported by the results presented in [4]. Here,
it was noted that the $L^{2}$-norm of the $\operatorname{div}\mathbf{B}_{h}$
may be affected by the residual threshold at which the iterations of a
preconditioned Krilov method are arrested. More precisely, the authors of [4]
noted that the higher this threshold is, the bigger the $L^{2}$ norm of
$\operatorname{div}\mathbf{B}_{h}$ is. Consequently, we can infer that the
divergence-free property of $\mathbf{B}_{h}$ is related to _how well_ the
linear system is solved and we claim that this effect on the $L^{2}$-norm of
$\operatorname{div}\mathbf{B}_{h}$ is probably due a possible growth of the
condition number of the final linear system. We use the direct solver PARDISO
[3]. Thus, the divergence free condition is not affected by any parameters of
the solver; rather, it is related _only_ to the round-off error.

$h/\tau$ | 1 | 1/2 | 1/4 | 1/8 | 1/16
---|---|---|---|---|---
rand27 | 5.57005e-14 | 2.33412e-14 | 1.43632e-14 | 1.31831e-14 | 1.01836e-14
rand125 | 4.47399e-13 | 3.14928e-13 | 1.27115e-13 | 1.26818e-13 | 1.07065e-13
rand1000 | 5.41757e-12 | 2.49285e-12 | 1.74853e-12 | 1.17031e-12 | 9.19428e-13
rand8000 | 7.24382e-10 | 3.80017e-10 | 6.07471e-10 | 2.25296e-10 | 3.78838e-09
Table 3: Test Case 1: $L^{2}$-norm of $\operatorname{div}\mathbf{B}_{h}$ at
final time $T=1$ using mesh family rand for various combinations of $h$ and
$\tau$.

###  6.2 Test case 2: polarized fields with variable coefficients

We solve Maxwell’s equations on the computational domain $\Omega=(0,\,1)^{3}$
for $t\in[0,1]$ with the variable coefficients

$\mu(x,\,y,\,z):=\frac{1}{1+x^{2}+y^{2}+z^{2}},\quad\varepsilon(x,\,y,\,z):=2-x^{2}-z,\quad\sigma(x,\,y,\,z):=2-y^{2}+z.$
(80)

The boundary conditions and the current density vector $\mathbf{J}$ are
defined in accordance with (80) and the exact solution fields

$\mathbf{E}(\mathbf{x},\,t):=\left(\begin{array}[]{c}0\\\ 0\\\
\sin(\pi\,x)\,\sin(\pi\,y)\end{array}\right)\cos(2.2\,\pi\,t),\quad\quad\mathbf{B}(\mathbf{x},\,t):=\left(\begin{array}[]{c}-\cos(\pi\,y)\,\sin(\pi\,x)\\\
\phantom{-}\cos(\pi\,x)\,\sin(\pi\,y)\\\
0\end{array}\right)\sin(2.2\,\pi\,t)/2.2.$

The electromagnetic fields $\mathbf{E}$ and $\mathbf{B}$ are orthogonal at any
point in $\Omega$ and time in $[0,1]$. Consequently, this solution simulates a
polarized stationary electromagnetic wave with a polarization direction that
is parallel to $\mathbf{E}\times\mathbf{B}$. We underline that this second
test case is more complex than the previous one since the coefficients $\mu$,
$\epsilon$ and $\sigma$ are all variables in space.

In Figure 3, we plot the $L^{2}$ errors at final time $T=1$ for simultaneous
refinements of $h$ and $\tau$ on the three mesh families and observe the
expected convergence rate, which we recall is expected to be proportional to
$h+\tau$; see Theorem 9. Moreover, in Figure 4 we report the convergence rate
at time $T=0.25,0.50$ and 0.75. Also in this case the behaviour of the error
is the one predicted by Theorem 9. We show such convergence lines only for
rand meshes as the results for the other type of meshes are similar.

|
---|---
Figure 3: Test Case 2: error curves in the $L^{2}$-norm at final time $T=1$ for the virtual element approximation $\mathbf{E}_{h}$ (left panel) and $\mathbf{B}_{h}$ (right panel) using mesh families cube, voro, and rand for simultaneous refinement of $h$ and $\tau$. |
---|---
Figure 4: Test Case 2: error curves in the $L^{2}$-norm at times
$T=0.25,0.50,0.75$ and $1.00$ for the virtual element approximation
$\mathbf{E}_{h}$ (left panel) and $\mathbf{B}_{h}$ (right panel) using the
mesh family rand for simultaneous refinement of $h$ and $\tau$.

As in Test Case 1, we observe similar convergence behavior of the proposed VEM
scheme on each mesh families so in Tables 4 and 5, we report the approximation
errors measured in the $L^{2}$ norms only for voro meshes and we omit the
results for the other two mesh families.

Despite the increased complexity due to variable coefficients, we observe the
optimal convergence behavior of the error also in this example. Indeed, each
column shows the convergence with respect to the space discretization, each
row shows the convergence with respect to the time discretization and the
diagonal shows the convergence when we refine simultaneously in space and
time. As before, the error of the space discretization appears to dominate the
error of the time discretization. Thus, the convergence in time along the
rows, which should be proportional to $\tau$, is not clearly visible.

$h/\tau$ | 1/8 | 1/16 | 1/32 | 1/64 | 1/128 | 1/256 | 1/512
---|---|---|---|---|---|---|---
voro27 | 8.64460e-01 | 6.85496e-01 | 5.61864e-01 | 5.03712e-01 | 4.81700e-01 | 4.73853e-01 | 4.70936e-01
voro125 | 8.55032e-01 | 6.30865e-01 | 4.51896e-01 | 3.55114e-01 | 3.16074e-01 | 3.02671e-01 | 2.98186e-01
voro1000 | 8.44007e-01 | 5.92062e-01 | 3.76326e-01 | 2.43624e-01 | 1.81720e-01 | 1.59270e-01 | 1.52316e-01
voro8000 | 8.40933e-01 | 5.81424e-01 | 3.54326e-01 | 2.05935e-01 | 1.27080e-01 | 9.33219e-02 | 8.18718e-02
Table 4: Test Case 2: $L^{2}$-norms of the error at final time $T=1$ for the virtual element approximation $\mathbf{E}_{h}$ using mesh family voro for various combinations of $h$ and $\tau$. $h/\tau$ | 1/8 | 1/16 | 1/32 | 1/64 | 1/128 | 1/256 | 1/512
---|---|---|---|---|---|---|---
voro27 | 5.92004e-01 | 5.69835e-01 | 5.43713e-01 | 5.34765e-01 | 5.37512e-01 | 5.42274e-01 | 5.45764e-01
voro125 | 4.68226e-01 | 4.28029e-01 | 3.69947e-01 | 3.23317e-01 | 3.02412e-01 | 2.95933e-01 | 2.94413e-01
voro1000 | 4.06835e-01 | 3.58563e-01 | 2.82887e-01 | 2.10580e-01 | 1.70641e-01 | 1.55564e-01 | 1.51081e-01
voro8000 | 3.90261e-01 | 3.38876e-01 | 2.55327e-01 | 1.69054e-01 | 1.13551e-01 | 8.82538e-02 | 7.96198e-02
Table 5: Test Case 2: $L^{2}$-norms of the error at final time $T=1$ for the
virtual element approximation $\mathbf{B}_{h}$ using mesh family voro for
various combinations of $h$ and $\tau$.

Finally, Table 6 shows the values of the $L^{2}$-norm of
$\operatorname{div}\mathbf{B}_{h}$: the VEM does preserve the solenoidal
property of the magnetic induction, i.e., the discrete field $\mathbf{B}_{h}$
has a pointwise zero divergence up to machine precision. If we compare the
results in Tables 3 and Table 6, then we note that the latters are smaller.
This is a further numerical evidence of the fact that the divergence-free
property is affected by the condition number of the resulting linear system.
Indeed, voro meshes are more shape-regular with respect to rand ones so the
condition numbers of matrices associated with them are smaller than those
associated with rand meshes: the algebraic errors are smaller and we get a
smaller divergence.

$h/\tau$ | 1/8 | 1/16 | 1/32 | 1/64 | 1/128 | 1/256 | 1/512
---|---|---|---|---|---|---|---
voro27 | 1.18487e-13 | 6.65328e-14 | 7.74812e-14 | 9.42188e-14 | 1.02437e-13 | 9.95126e-14 | 1.04040e-13
voro125 | 1.79095e-15 | 4.39464e-15 | 2.87788e-15 | 4.01614e-15 | 6.12453e-15 | 8.84695e-15 | 1.02443e-14
voro1000 | 1.98678e-14 | 3.86910e-14 | 3.81884e-14 | 3.34432e-14 | 3.23360e-14 | 2.54928e-14 | 2.88468e-14
voro8000 | 1.75878e-13 | 5.95562e-13 | 2.50332e-13 | 2.40717e-13 | 1.63504e-13 | 1.48763e-13 | 9.27908e-14
Table 6: Test Case 2: $L^{2}$ norm of $\operatorname{div}\mathbf{B}_{h}$ at
final time $T=1$ using mesh family voro for various combinations of $h$ and
$\tau$.

##  7 Conclusions

In this paper, we have considered a low order virtual element approximation of
Maxwell’s equations based on a De Rahm sequence. After developing some
interpolation and stability properties of edge and face spaces, we showed
optimal a priori estimates for both the semi- and the fully- discrete schemes
and corroborated the theoretical predictions with numerical experiments.
Future works may cover the approximation of corner singularities and the
virtual element approximation of MHD problems. The extension to high order
methods requires high order interpolation estimates and stability properties
of edge and face VEM spaces, which is currently a work in progress.

## Acknowledgments

L. B. da V. and F. D. are partially supported by the European Research Council
through the H2020 Consolidator Grant (Grant No. 681162) CAVE - Challenges and
Advancements in Virtual Elements. L. B. da V. is also partially supported by
the MIUR through the PRIN grant n. 201744KLJL. G. M. has partially been
supported by the ERC Project CHANGE, which has received funding from the
European Research Council (ERC) under the European Union Horizon 2020 research
and innovation program (grant agreement no. 694515). L. M. acknowledges
support from the Austrian Science Fund (FWF) project P33477.

We further wish to thank Martin Costabel for an advice regarding a regularity
result.

## References

  * [1] R. A. Adams and J. J. F. Fournier.  Sobolev Spaces, volume 140.  Academic Press, 2003.
  * [2] B. Ahmad, A. Alsaedi, F. Brezzi, L.D. Marini, and A. Russo.  Equivalent projectors for virtual element methods.  Comput. Math. Appl., 66(3):376–391, 2013.
  * [3] Christie Alappat, Achim Basermann, Alan R. Bishop, Holger Fehske, Georg Hager, Olaf Schenk, Jonas Thies, and Gerhard Wellein.  A recursive algebraic coloring technique for hardware-efficient symmetric sparse matrix-vector multiplication.  ACM Trans. Parallel Comput., 7(3), June 2020.
  * [4] S. N. Alvarez, V. A. Bokil, V. Gyrya, and G. Manzini.  The virtual element method for resistive magnetohydrodynamics.  Comput. Methods Appl. Mech. Engrg., 381:113815, 2021.
  * [5] C. Amrouche, C. Bernardi, M. Dauge, and V. Girault.  Vector potentials in three-dimensional non-smooth domains.  Math. Methods Appl. Sci., 21(9):823–864, 1998.
  * [6] F. Assous, P. Degond, E. Heintze, P.-A. Raviart, and J. Segre.  On a finite-element method for solving the three-dimensional Maxwell equations.  J. Comput. Phys., 109(2):222–237, 1993.
  * [7] L. Beirão da Veiga and L. Mascotto.  Interpolation and stability properties of low order face and edge virtual element spaces.  https://arxiv.org/abs/2011.12834, 2020.
  * [8] L. Beirão da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L.D. Marini, and A. Russo.  Basic principles of virtual element methods.  Math. Models Methods Appl. Sci., 23(01):199–214, 2013.
  * [9] L. Beirão da Veiga, F. Brezzi, F. Dassi, L. D. Marini, and A. Russo.  Virtual element approximation of 2D magnetostatic problems.  Comput. Methods Appl. Mech. Engrg., 327:173–195, 2017.
  * [10] L. Beirão da Veiga, F. Brezzi, F. Dassi, L. D. Marini, and A. Russo.  A family of three-dimensional virtual elements with applications to magnetostatics.  SIAM J. Numer. Anal., 56(5):2940–2962, 2018.
  * [11] L. Beirão da Veiga, F. Brezzi, F. Dassi, L. D. Marini, and A. Russo.  Lowest order virtual element approximation of magnetostatic problems.  Comput. Methods Appl. Mech. Engrg., 332:343–362, 2018.
  * [12] L. Beirão Da Veiga, F. Brezzi, L. D. Marini, and A. Russo.  H(div) and H(curl)-conforming virtual element methods.  Numer. Math, 133(2):303–332, 2016.
  * [13] L. Beirao da Veiga, F. Brezzi, L. D. Marini, and A. Russo.  Mixed virtual element methods for general second order elliptic problems on polygonal meshes.  ESAIM Math. Model. Numer. Anal., 50(3):727–747, 2016.
  * [14] L. Beirão da Veiga, F. Dassi, and A. Russo.  High-order virtual element method on polyhedral meshes.  Comput. Math. Appl., 74(5):1110–1122, 2017.
  * [15] L. Beirão da Veiga, C. Lovadina, and A. Russo.  Stability analysis for the virtual element method.  Math. Models Methods Appl. Sci., 27(13):2557–2594, 2017.
  * [16] A. Bermúdez de Castro, D. Gómez, and P. Salgado.  Mathematical models and numerical simulation in electromagnetism, volume 74.  Springer, 2014.
  * [17] D. Boffi, F. Brezzi, and M. Fortin.  Mixed Finite Element Methods and Applications, volume 44.  Springer Series in Computational Mathematics, 2013.
  * [18] S. C. Brenner and L.-Y. Sung.  Virtual element methods on meshes with small edges or faces.  Math. Models Methods Appl. Sci., 268(07):1291–1336, 2018.
  * [19] F. Brezzi, R. S. Falk, and L.D. Marini.  Basic principles of mixed virtual element methods.  Math. Mod. Num. Anal., 48(4):1227–1240, 2014.
  * [20] F. Chave, D. A. Di Pietro, and S. Lemaire.  A three-dimensional hybrid high-order method for magnetostatics.  In International Conference on Finite Volumes for Complex Applications, pages 255–263. Springer, 2020.
  * [21] P. Ciarlet, Jr. and J. Zou.  Fully discrete finite element approaches for time-dependent Maxwell’s equations.  Numer. Math., 82(2):193–219, 1999.
  * [22] R. Coccioli, T. Itoh, G. Pelosi, and P. P. Silvester.  Finite-element methods in microwaves: A selected bibliography.  Antennas and Propagation Newsletter, IEEE Professional Group on, 38, Dec. 1996.
  * [23] F. Dassi, P. Di Barba, and A. Russo.  Virtual element method and permanent magnet simulations: potential and mixed formulations.  IET Science, Measurement & Technology, 14(10):1098–1104, 2021\.
  * [24] M. Dauge.  Elliptic boundary value problems on corner domains: smoothness and asymptotics of solutions, volume 1341.  Springer, 2006.
  * [25] D. A. Di Pietro, J. Droniou, and F. Rapetti.  Fully discrete polynomial de Rham sequences of arbitrary degree on polygons and polyhedra.  Math. Models Methods Appl. Sci., 30(9):1809–1855, 2020.
  * [26] T. Euler, R. Schuhmann, and T. Weiland.  Polygonal finite elements.  IEEE Trans. Magnetics, 42, 2006.
  * [27] L. C. Evans.  Partial Differential Equations.  American Mathematical Society, 2010.
  * [28] A. D. Greenwood and J.-M. Jin.  Finite-element analysis of complex axisymmetric radiating structures.  IEEE Trans. Antennas Propagation, 47(8):1260–1266, 1999.
  * [29] J.-M. Jin.  The finite element method in electromagnetics.  John Wiley & Sons, third edition edition, 2014.
  * [30] J.-M. Jin and D. J. Riley.  Finite Element Analysis of Antennas and Arrays.  John Wiley & Sons :, IEEE Press, 2009.
  * [31] A. Khebir, J. D’Angelo, and J. Joseph.  A new finite element formulation for RF scattering by complex bodies of revolution.  IEEE Trans. Antennas Propagation, 41(5):534–541, 1993.
  * [32] J. F. Lee, G.M. Wilkins, and R. Mitra.  Finite-element analysis of axisymmetric cavity resonator using a hybrid edge element technique.  IEEE Trans. Microwave Theory Techniques, 41:1981–1987, 1993.
  * [33] K. Lipnikov, G. Manzini, F. Brezzi, and A. Buffa.  The mimetic finite difference method for the 3D magnetostatic field problems on polyhedral meshes.  J. Comput. Phys., 230(2):305–328, 2011.
  * [34] Ch. G. Makridakis and P. Monk.  Time-discrete finite element schemes for Maxwell’s equations.  RAIRO Modél. Math. Anal. Numér., 29(2):171–197, 1995.
  * [35] L. Medgyesi-Mitschang and J. Putnam.  Electromagnetic scattering from axially inhomogeneous bodies of revolution.  IEEE Trans. Antennas Propagation, 32(8):797–806, 1984.
  * [36] P. Monk.  Finite Element Methods for Maxwell’s Equations.  Oxford University Press, 2003.
  * [37] P. B. Monk.  A mixed method for approximating Maxwell’s equations.  SIAM J. Numer. Anal., 28(6):1610–1634, 1991.
  * [38] D. Mora, G. Rivera, and R. Rodríguez.  A virtual element method for the Steklov eigenvalue problem.  Math. Models Methods Appl. Sci., 25(08):1421–1445, 2015.
  * [39] J.-C. Nédélec.  Mixed finite elements in ${\bf R}^{3}$.  Numer. Math., 35(3):315–341, 1980.
  * [40] X. Rui, J. Hu, and Q. H. Liu.  Higher order finite element method for inhomogeneous axisymmetric resonators.  Progress In Electromagnetics Research B, 21:189–201, 2010.
  * [41] F. L. Teixeira and J. R. Bergmann.  B-spline basis functions for moment-method analysis of axisymmetric reflector antennas.  Microwave and Optical Technology Letters, 14, 1997.
  * [42] F. L. Teixeira and J. R. Bergmann.  Moment-method analysis of circularly symmetric reflectors using bandlimited basis functions.  IEE Proc. Microw. Antennas Prop., 144(3):179–183, 1997.
  * [43] W. Tierens and D. De Zutter.  BOR-FDTD subgridding based on finite element principles.  J. Comput. Phys., 230, 2011.
  * [44] G.M. Wilkins, J.-F. Lee, and R. Mittra.  Numerical modeling of axisymmetric coaxial waveguide discontinuities.  IEEE Trans. Microwave Theory Techniques, 39, 1991.
  * [45] J. Zhao.  Analysis of finite element approximation for time-dependent Maxwell problems.  Math. Comp., 73(247):1089–1105, 2004.