# Adaptive virtual elements based on hybridized, reliable, and efficient flux
reconstructions

Franco Dassi and Joscha Gedicke and Lorenzo Mascotto  Dip. di Matematica e
Applicazioni, Università degli Studi di Milano-Bicocca, Italy
<EMAIL_ADDRESS>Institut für Numerische Simulation, Universität Bonn,
53115 Bonn<EMAIL_ADDRESS>Fakultät für Mathematik, Universität Wien,
1090 Vienna, Austria<EMAIL_ADDRESS>
###### Abstract.

We present two a posteriori error estimators for the virtual element method
(VEM) based on global and local flux reconstruction in the spirit of [6]. The
proposed error estimators are reliable and efficient for the $h$-, $p$-, and
$hp$-versions of the VEM. This solves a partial limitation of our former
approach in [7], which was based on solving local nonhybridized mixed
problems. Differently from the finite element setting, the proof of the
efficiency turns out to be simpler, as the flux reconstruction in the VEM does
not require the existence of polynomial, stable, divergence right-inverse
operators. Rather, we only need to construct right-inverse operators in
virtual element spaces, exploiting only the implicit definition of virtual
element functions. The theoretical results are validated by some numerical
experiments on a benchmark problem.

AMS subject classification: 65N12, 65N30, 65N5

Keywords: virtual element method, hybridized methods, equilibrated fluxes,
$hp$-adaptivity, polygonal meshes

##  1\. Introduction

In [7], we analyzed certain error estimators for the virtual element method
(VEM) [3], based on global and local flux reconstruction via mixed problems.
To the aim, we mainly followed the finite element (FE) and discontinuous
Galerkin (dG) approach of [9] and realized that the VEM contains too many
variational crimes to allow for a localized, efficient error estimator with
such an approach. In this paper, we undertake a slightly different path, i.e.,
we design error estimators by means of flux reconstruction via hybridized
mixed methods. Notably, we adapt the finite element framework [6] to the
virtual element one, which is based on polygonal meshes that are more suited
for mesh adaptation as they naturally handle hanging nodes.

We shall design two error estimators based on global and local flux
reconstructions, and prove their reliability and efficiency for the $h$-,
$p$-, and $hp$ versions of the VEM. In both cases, the proof of the efficiency
turns out to be more straightforward than the finite element one. Indeed, when
using virtual elements, the existence of a polynomial, divergence right-
inverse operator, see [6, Theorem $5$ and Conjecture $6$], is not required.

As a pivot result, we also prove the existence of a stabilization term for
mixed virtual elements, with stability constants having explicit dependence on
the degree of accuracy of the method. Furthermore, we discuss in details the
flux reconstruction, well posedness of the hybridized mixed VEM; see also [8]
for the hybridization of the mixed VEM in linear elasticity.

##### Notation.

We consider standard notation for Sobolev spaces. Given a domain $D$,
$H^{s}(D)$ denotes the standard Sobolev space of order $s$, which we endow
with the standard inner product $(\cdot,\cdot)_{s,D}$, norm $\|\cdot\|_{s,D}$,
and seminorm $|\cdot|_{s,D}$. Fractional and negative Sobolev spaces are
defined via interpolation and duality, respectively. Recall the definition of
standard differential equations in two dimensions: given
$v:\mathbb{R}^{2}\to\mathbb{R}$ and
$\boldsymbol{\tau}=(\boldsymbol{\tau}_{1},\boldsymbol{\tau}_{2}):\mathbb{R}^{2}\to\mathbb{R}^{2}$,

$\begin{split}\operatorname{div}\boldsymbol{\tau}:=\partial_{x}\boldsymbol{\tau}_{1}+\partial_{y}\boldsymbol{\tau}_{2},\qquad\operatorname{rot}\boldsymbol{\tau}:=\partial_{y}\boldsymbol{\tau}_{1}-\partial_{x}\boldsymbol{\tau}_{2},\qquad\operatorname{\bf{curl}}v:=(\partial_{y}v;-\partial_{x}v).\end{split}$

In light of the above definitions, we also introduce the spaces

$H(\operatorname{div},D):=\\{\boldsymbol{\tau}\in[L^{2}(D)]^{2}\mid\operatorname{div}\boldsymbol{\tau}\in
L^{2}(D)\\},\quad
H(\operatorname{rot},D):=\\{\boldsymbol{\tau}\in[L^{2}(D)]^{2}\mid\operatorname{rot}\boldsymbol{\tau}\in
L^{2}(D)\\}.$

Given two positive quantities $a$ and $b$, we shall write “$a\lesssim b$”
instead of “there exists a positive constant $C$, independent of the
discretization parameters, such that $a\leq C\,b$”.

##### The model problem.

Given a polygonal domain $\Omega\subset\mathbb{R}^{2}$ and $f\in
L^{2}(\Omega)$, we consider the Poisson problem: find $\widetilde{u}$ such
that

$\begin{cases}-\Delta\widetilde{u}=f&\text{in }\Omega\\\\[3.0pt]
\widetilde{u}=0&\text{on }\partial\Omega,\end{cases}$ (1)

where $\partial\Omega$ denotes the boundary of $\Omega$.

The weak formulation of (1) reads

$\begin{cases}\text{find }\widetilde{u}\in\widetilde{V}\text{ such that}\\\
\widetilde{a}(\widetilde{u},\widetilde{v})=(f,\widetilde{v})_{0,\Omega}\qquad\forall\widetilde{v}\in\widetilde{V},\end{cases}$
(2)

where we have set

$\begin{split}&\widetilde{V}:=H^{1}_{0}(\Omega),\qquad\widetilde{a}(\widetilde{u},\widetilde{v}):=(\nabla\widetilde{u},\nabla\widetilde{v})_{0,\Omega}\qquad\forall\widetilde{u},\widetilde{v}\in\widetilde{V}.\end{split}$

Problem (2) is well posed thanks to the Lax-Milgram lemma.

##### Mesh assumptions.

We employ standard notation for polygonal meshes. More precisely, we consider
sequences $\\{\mathcal{T}_{n}\\}_{n\in\mathbb{N}}$ of polygonal meshes with
straight edges. Hanging nodes are naturally incorporated in such meshes. Given
a mesh $\mathcal{T}_{n}$, we denote its set of edges and vertices by
$\mathcal{E}_{n}$ and $\mathcal{V}_{n}$. We split $\mathcal{E}_{n}$ into the
set of boundary and internal edges as follows:

$\mathcal{E}_{n}^{D}:=\\{e\in\mathcal{E}_{n}\mid
e\subset\partial\Omega\\},\qquad\mathcal{E}_{n}^{I}:=\mathcal{E}_{n}\setminus\mathcal{E}_{n}^{D}.$

With each element $K\in\mathcal{T}_{n}$, we associate its diameter $h_{K}$,
the outward pointing unit normal vector $\mathbf{n}_{K}$, its centroid
$\mathbf{x}_{K}$, and its set of edges $\mathcal{E}^{K}$ and vertices
$\mathcal{V}^{K}$. Moreover, we associate with each edge $e\in\mathcal{E}_{n}$
its length $h_{e}$ and any of the two unit normal vectors $\mathbf{n}_{e}$.

We consider standard properties for the sequence of meshes [3]: there exist
two positive constants $\gamma$ and $\widetilde{\gamma}$, such that, for all
$n\in\mathbb{N}$,

  * (A1)

every $K\in\mathcal{T}_{n}$ is star-shaped with respect to a ball $B$ with
radius $r_{B}\geq\gamma h_{K}$;

  * (A2)

given $K\in\mathcal{T}_{n}$, for every of its edges $e\in\mathcal{E}^{K}$,
$h_{K}\leq\widetilde{\gamma}h_{e}$.

Associated with each mesh $\mathcal{T}_{n}$ and $p\in\mathbb{N}$, we define
standard broken polynomial and Sobolev spaces:

$\begin{split}&\mathbb{P}_{p}(\mathcal{T}_{n}):=\\{q_{p}\in L^{2}(\Omega)\mid
q_{p}{}_{|K}\in\mathbb{P}_{p}(K)\;\forall K\in\mathcal{T}_{n}\\};\\\\[3.0pt]
&H^{1}(\mathcal{T}_{n}):=\\{v\in L^{2}(\Omega)\mid v{}_{|K}\in
H^{1}(K)\;\forall K\in\mathcal{T}_{n}\\};\\\\[3.0pt]
&H(\operatorname{div},\mathcal{T}_{n}):=\\{\boldsymbol{\tau}\in[L^{2}(\Omega)]^{2}\mid\boldsymbol{\tau}_{|K}\in
H(\operatorname{div},\Omega)\\}.\end{split}$

The first and second spaces can be generalized to the corresponding vector
spaces.

Henceforth, for ease of presentation, we assume that
$f\in\mathbb{P}_{p-1}(\mathcal{T}_{n})$. Moreover, throughout the paper, we
assume a uniform distribution of the polynomial degree $p$ for the sake of
presentation. The general case involving nonuniform distributions of
polynomial degrees is a trivial generalization.

##### Structure of the paper.

In Section 2, we review the hybridized global and local flux reconstruction in
the finite element setting. After describing virtual elements for the primal
and the hybridized mixed formulations in Section 3, we design and analyze the
hybridized global and local flux reconstructions in the virtual element
setting in Section 4. Notably, we introduce a reliable and efficient error
estimator. Section 5 devotes to presenting numerical experiments validating
the theoretical predictions. Finally, we state conclusions in Section 6

##  2\. Hybridized global/local flux reconstruction in the finite element
setting

Here, we briefly recall the approach of [6] providing $p$ robust, reliable and
efficient error estimators on simplicial and Cartesian meshes for the finite
element method. This will pave the way to design reliable and efficient error
estimators for the virtual element method.

Let $\mathcal{T}_{n}$ be a decomposition of $\Omega$ into triangles. We
consider the standard finite element approximation of (2):

$\begin{cases}\text{find }\widetilde{u}_{p}\in\widetilde{V}_{p}\text{ such
that}\\\
\widetilde{a}(\widetilde{u}_{p},\widetilde{v}_{p})=(f,\widetilde{v}_{p})_{0,\Omega}\qquad\forall\widetilde{v}_{p}\in\widetilde{V}_{p},\end{cases}$
(3)

where we have set

$\widetilde{V}_{p}:=\mathbb{P}_{p}(\mathcal{T}_{n})\cap H^{1}_{0}(\Omega).$

We define a continuous residual $r\in[H^{1}_{0}(\Omega)]^{*}$ as

$\langle
r,\widetilde{v}\rangle:=(f,\widetilde{v})_{0,\Omega}-\widetilde{a}(\widetilde{u}_{p},\widetilde{v})\qquad\forall\widetilde{v}\in
H^{1}_{0}(\Omega).$ (4)

An integration by parts yields

$\langle
r,\widetilde{v}\rangle:=\sum_{K\in\mathcal{T}_{n}}(r^{K},\widetilde{v})_{0,K}+\sum_{e\in\mathcal{E}_{n}}(r^{e},\widetilde{v})_{0,e}\qquad\forall\widetilde{v}\in\widetilde{V},$

where we have introduced the bulk and internal edge residuals

$r^{K}:=f+\Delta\widetilde{u}_{p}\in\mathbb{P}_{p-1}(K)\quad\forall
K\in\mathcal{T}_{n},\qquad
r^{e}:=\llbracket\nabla\widetilde{u}_{p}\rrbracket\in\mathbb{P}_{p-1}(e)\quad\forall
e\in\mathcal{E}_{n}^{I},$ (5)

for the standard definition of the jump operator on an internal edge $e$:

$\llbracket\boldsymbol{\tau}\rrbracket_{e}:=\llbracket\boldsymbol{\tau}\rrbracket:=\mathbf{n}_{K_{1}}\cdot\boldsymbol{\tau}_{K_{1}}+\mathbf{n}_{K_{2}}\cdot\boldsymbol{\tau}_{K_{2}},\quad
e\subset\overline{K_{1}}\cap\overline{K_{2}},\qquad\forall\boldsymbol{\tau}_{p}\in
H(\operatorname{div},\mathcal{T}_{n}).$

The initial step towards the design of the hybridized flux reconstruction
consists in letting the bulk and edge residual belong to different polynomial
spaces. Notably, we assume that $r^{K}\in\widetilde{V}_{p}$, whereas we let
$r^{e}\in\Lambda_{p}^{I}$, where we have set
$\Lambda_{p}^{I}:=\mathbb{P}_{p}(\mathcal{E}_{n}^{I})$.

Define a flux $\boldsymbol{\sigma}_{p-1}^{\Delta}$ in a distributional sense
as follows:

$\operatorname{div}\boldsymbol{\sigma}_{p-1}^{\Delta}=r\in(H^{1}_{0}(\Omega))^{*}.$
(6)

We can interpret $\boldsymbol{\sigma}_{p-1}^{\Delta}$ as follows:

$\boldsymbol{\sigma}_{p}^{\Delta}\approx\nabla\widetilde{u}_{p}-\nabla\widetilde{u}\in(H^{1}_{0}(\Omega))^{*}.$

Using (4) and (6) implies that $\boldsymbol{\sigma}_{p-1}^{\Delta}$ satisfies

$\begin{cases}\operatorname{div}\boldsymbol{\sigma}_{p-1}^{\Delta}=r^{K}\in\mathbb{P}_{p-1}(K)&\text{in
all }K\in\mathcal{T}_{n}\\\\[3.0pt]
\llbracket\boldsymbol{\sigma}_{p-1}^{\Delta}\rrbracket=r^{e}\in\mathbb{P}_{p-1}(e)&\text{on
all }e\in\mathcal{E}_{n}^{I}.\end{cases}$ (7)

The information about the homogeneous Dirichlet boundary conditions is
implicitly contained in the two residuals through $\widetilde{u}_{p}$. We
rewrite (7) as a hybridized mixed method with flux and primal variable given
by $\boldsymbol{\sigma}_{p}^{\Delta}$ and $u_{p-1}$, respectively. More
precisely, we write

$\begin{cases}\text{find
}(\boldsymbol{\sigma}_{p-1}^{\Delta},u_{p-1},\lambda_{p-1}^{I})\in\boldsymbol{\Sigma}_{p-1}\times
V_{p-1}\times\Lambda_{p-1}^{I}\text{ such that}\\\\[3.0pt]
a(\boldsymbol{\sigma}_{p-1}^{\Delta},\boldsymbol{\tau}_{p-1})+b(\boldsymbol{\tau}_{p-1},u_{p-1})+c(\boldsymbol{\tau}_{p-1},\lambda_{p-1}^{I})=0&\forall\boldsymbol{\tau}_{p-1}\in\boldsymbol{\Sigma}_{p-1}\\\\[3.0pt]
b(\boldsymbol{\sigma}_{p-1}^{\Delta},v_{p-1})=(r^{K},v_{p-1})_{0,\Omega}&\forall
v_{p-1}\in V_{p-1}\\\\[3.0pt]
c(\boldsymbol{\sigma}_{p-1}^{\Delta},\mu_{p-1}^{I})=(r^{e},\mu_{p-1}^{I})_{0,\mathcal{E}_{n}^{I}}&\forall\mu_{p-1}^{I}\in\Lambda_{p-1}^{I},\end{cases}$
(8)

where we picked the spaces

$\begin{split}\boldsymbol{\Sigma}_{p-1}&:=\left\\{\boldsymbol{\tau}_{p-1}\in[L^{2}(\Omega)]^{2}\mid\boldsymbol{\tau}_{p-1}{}_{|K}\in\mathbb{RT}_{p-1}(K)\;\forall
K\in\mathcal{T}_{n}\right\\},\\\
V_{p-1}&:=\mathbb{P}_{p-1}(\mathcal{T}_{n}),\\\
\Lambda_{p-1}&:=\mathbb{P}_{p-1}(\mathcal{E}_{n}^{I}),\end{split}$ (9)

being $\mathbb{RT}_{p}(K)$ the local Raviart-Thomas space of order
$p\in\mathbb{N}_{0}$ on element $K$, and the bilinear forms

$\begin{split}a(\boldsymbol{\sigma}_{p-1},\boldsymbol{\tau}_{p-1})&:=(\boldsymbol{\sigma}_{p-1},\boldsymbol{\tau}_{p-1})_{0,\Omega},\\\
b(\boldsymbol{\tau}_{p-1},u_{p-1})&:=-(\operatorname{div}\boldsymbol{\tau}_{p-1},u_{p-1})_{0,\Omega},\\\
c(\boldsymbol{\tau}_{p-1},\mu_{p-1}^{I})&:=(\llbracket\boldsymbol{\tau}_{p-1}\rrbracket,\mu_{p-1}^{I})_{0,\mathcal{E}_{n}^{I}}.\end{split}$
(10)

In the literature, the space $\boldsymbol{\Sigma}_{p-1}$ is known as the
broken, hybridized Raviart-Thomas space of order $p-1$. No degrees of freedom
coupling takes place at the interface between two elements. Rather, the
communication between local Raviart-Thomas spaces takes place in a mortar
fashion through the Lagrange multipliers in the space $\Lambda_{p-1}$.

Next, we discuss the localized version of the above problem. With each vertex
$\nu\in\mathcal{V}_{n}$, we associate the partition of unity function
$\widetilde{\varphi}_{\nu}\in\widetilde{V}_{p}$, i.e., the standard hat
function with value $1$ at $\nu$ and $0$ at all other vertices, and localize
the bulk and edge residuals in (5) as follows: for all
$\nu\in\mathcal{V}_{n}$,

$r^{K}_{\nu}:=\widetilde{\varphi}_{\nu}(f+\Delta\widetilde{u}_{p})\in\mathbb{P}_{p}(K)\quad\forall
K\in\mathcal{T}_{n},\qquad
r^{e}_{\nu}:=\widetilde{\varphi}_{\nu}\llbracket\nabla\widetilde{u}_{p}\rrbracket\in\mathbb{P}_{p}(e)\qquad\forall
e\in\mathcal{E}_{n}^{I}.$ (11)

Proceeding as above, on each vertex patch

$\omega_{\nu}:=\bigcup\\{K\in\mathcal{T}_{n}\mid\nu\text{ is a vertex of
}K\\},$ (12)

we introduce the following local problems: if $\nu$ is an interior vertex,

$\begin{cases}\text{find
}(\boldsymbol{\sigma}_{p}^{\nu},u_{p}^{\nu},\lambda_{p}^{\nu,I})\in\boldsymbol{\Sigma}_{p}^{0}(\omega_{\nu})\times
V_{p}^{*}(\omega_{\nu})\times\Lambda_{p}^{I}(\omega_{\nu})\text{ such
that}\\\\[3.0pt]
a^{\omega_{\nu}}(\boldsymbol{\sigma}_{p}^{\nu},\boldsymbol{\tau}_{p}^{\nu})+b^{\omega_{\nu}}(\boldsymbol{\tau}_{p}^{\nu},u_{p}^{\nu})+c^{\omega_{\nu}}(\boldsymbol{\tau}_{p}^{\nu},\lambda_{p}^{\nu,I})=0&\forall\boldsymbol{\tau}_{p}^{\nu}\in\boldsymbol{\Sigma}_{p}^{0}(\omega_{\nu})\\\\[3.0pt]
b^{\omega_{\nu}}(\boldsymbol{\sigma}_{p}^{\nu},v_{p}^{\nu})=(r^{K}_{\nu},v_{p}^{\nu})_{0,\omega_{\nu}}&\forall
v_{p}^{\nu}\in V_{p}^{*}(\omega_{\nu})\\\\[3.0pt]
c^{\omega_{\nu}}(\boldsymbol{\sigma}_{p}^{\nu},\mu_{p}^{\nu,I})=(r^{e}_{\nu},\mu_{p}^{\nu,I})_{0,\mathcal{E}_{n}^{I}(\omega_{\nu})}&\forall\mu_{p}^{\nu,I}\in\Lambda_{p}^{I}(\mathcal{E}_{n}^{I}(\omega_{\nu})),\end{cases}$

whereas, if $\nu$ is a boundary vertex,

$\begin{cases}\text{find
}(\boldsymbol{\sigma}_{p}^{\nu},u_{p}^{\nu},\lambda_{p}^{\nu,I})\in\boldsymbol{\Sigma}_{p}^{0,\partial\Omega}(\omega_{\nu})\times
V_{p}(\omega_{\nu})\times\Lambda_{p}^{I}(\omega_{\nu})\text{ such
that}\\\\[3.0pt]
a^{\omega_{\nu}}(\boldsymbol{\sigma}_{p}^{\nu},\boldsymbol{\tau}_{p}^{\nu})+b^{\omega_{\nu}}(\boldsymbol{\tau}_{p}^{\nu},u_{p}^{\nu})+c^{\omega_{\nu}}(\boldsymbol{\tau}_{p}^{\nu},\lambda_{p}^{\nu,I})=0&\forall\boldsymbol{\tau}_{p}^{\nu}\in\boldsymbol{\Sigma}_{p}^{0,\partial\Omega}(\omega_{\nu})\\\\[3.0pt]
b^{\omega_{\nu}}(\boldsymbol{\sigma}_{p}^{\nu},v_{p})=(r^{K}_{\nu},v_{p}^{\nu})_{0,\omega_{\nu}}&\forall
v_{p}^{\nu}\in V_{p}(\omega_{\nu})\\\\[3.0pt]
c^{\omega_{\nu}}(\boldsymbol{\sigma}_{p}^{\nu},\mu_{p}^{\nu,I})=(r^{e}_{\nu},\mu_{p}^{\nu,I})_{0,\mathcal{E}_{n}^{I}(\omega_{\nu})}&\forall\mu_{p}^{\nu,I}\in\Lambda_{p}^{I}(\mathcal{E}_{n}^{I}(\omega_{\nu})),\end{cases}$

where the involved spaces and bilinear forms are the obvious counterparts of
those in (9) and (10), and

$\begin{split}\boldsymbol{\Sigma}_{p}^{0}(\omega_{\nu})&:=\\{\boldsymbol{\tau}_{p}\in\boldsymbol{\Sigma}_{p}(\omega_{\nu})\mid\mathbf{n}_{e}\cdot\boldsymbol{\tau}_{p}{}_{|e}=0\quad\forall
e\subset\partial\omega_{\nu}\\},\\\\[3.0pt]
\boldsymbol{\Sigma}_{p}^{0,\partial\Omega}(\omega_{\nu})&:=\\{\boldsymbol{\tau}_{p}\in\boldsymbol{\Sigma}_{p}(\omega_{\nu})\mid\mathbf{n}_{e}\cdot\boldsymbol{\tau}_{p}{}_{|e}=0\quad\forall
e\subset\partial\omega_{\nu},\,\nu\nsubseteq e\\},\\\\[3.0pt]
V_{p}^{*}(\omega_{\nu})&:=\left\\{v_{p}\in
V_{p}(\omega_{\nu})\Big{|}\int_{\omega_{\nu}}v_{p}=0\right\\}.\end{split}$

In [6], it was proven that the error estimator

$\boldsymbol{\sigma}_{p}:=\sum_{\nu\in\mathcal{V}_{n}}\boldsymbol{\sigma}_{p}^{\nu}$

is reliable and efficient with constants independent of $h$ and $p$.

In Sections 3 and 4 below, we shall generalize the finite element setting to
the virtual element one.

##  3\. The virtual element setting

Here, we introduce several virtual element spaces, which we shall use in the
design of the hybridized local flux reconstruction in Section 4 below. Primal
virtual elements are detailed in Section 3.1, whereas we design couples of
broken, mixed virtual elements in Section 3.2. We also introduce skeletal,
broken polynomial spaces in Section 3.3, that will play the role of the
hybridization space. As a matter of notation, we shall denote the primal,
scalar variable with a tilde on top when considered as solution to the primal
formulation, and without tilde when appearing in the mixed formulation.

###  3.1. Primal virtual elements

Here, we follow [3] and briefly recall the primal virtual element method for
the Poisson problem (2).

####  3.1.1. Primal virtual element spaces.

Given $p\in\mathbb{N}$, for all $K\in\mathcal{T}_{n}$, we introduce the primal
local space

$\widetilde{V}_{p}(K):=\\{\widetilde{v}_{p}\in\mathcal{C}^{0}(\overline{K})\mid\widetilde{v}_{p}{}_{|e}\in\mathbb{P}_{p}(e)\,\forall
e\in\mathcal{E}_{n},\;\Delta\widetilde{v}_{p}\in\mathbb{P}_{p-2}(K)\\}.$

We have that $\mathbb{P}_{p}(K)\subseteq\widetilde{V}_{p}(K)$. Next, consider
$\\{m^{K}_{\boldsymbol{\alpha}}\\}_{|\boldsymbol{\alpha}|=0}^{p-2}$ any basis
of $\mathbb{P}_{p-2}(K)$, whose elements are invariant with respect to
translations and dilations. The following set of linear functionals of
$\widetilde{V}_{p}(K)$ is a set of unisolvent degrees of freedom [3]: given
$\widetilde{v}_{p}\in\widetilde{V}_{p}(K)$,

  * •

the point values of $\widetilde{v}_{p}$ at the vertices of $K$;

  * •

for each $e\in\mathcal{E}^{K}$, the point values of $\widetilde{v}_{p}$ at the
$p$ internal Gauß-Lobatto nodes of $e$;

  * •

the scaled moments

$\frac{1}{|K|}\int_{K}\widetilde{v}_{p}m^{K}_{\boldsymbol{\alpha}}\qquad\forall|\boldsymbol{\alpha}|=0,\dots,p-2.$

The global virtual element space together with its set of degrees of freedom
is constructed via a standard $H^{1}$-conforming coupling of their local
counterparts and strongly imposing the Dirichlet boundary conditions:

$\widetilde{V}_{p}:=\\{\widetilde{v}_{p}\in\mathcal{C}^{0}(\overline{\Omega})\cap
H^{1}_{0}(\Omega)\mid\widetilde{v}_{p}{}_{|K}\in\widetilde{V}_{p}(K)\;\forall
K\in\mathcal{T}_{n}\\}.$

####  3.1.2. Projectors

Through the degrees of freedom detailed in Section 3.1.1, it is possible to
compute two local orthogonal projectors onto polynomial space. The first one
is the energy projector
$\widetilde{\Pi}^{\nabla}_{p}:\widetilde{V}_{p}(K)\rightarrow\mathbb{P}_{p}(K)$
defined as follows: given $\widetilde{v}_{p}$ in $\widetilde{V}_{p}(K)$,

$\int_{K}\nabla
q^{K}_{p}\cdot\nabla(\widetilde{v}_{p}-\widetilde{\Pi}^{\nabla}_{p}\widetilde{v}_{p})=0\quad\forall
q^{K}_{p}\in\mathbb{P}_{p}(K),\qquad\int_{\partial
K}(\widetilde{v}_{p}-\widetilde{\Pi}^{\nabla}_{p}\widetilde{v}_{p})=0.$

The second one is the $L^{2}$ projector
$\widetilde{\Pi}^{0}_{p-2}:\widetilde{V}_{p}(K)\rightarrow\mathbb{P}_{p-2}(K)$
defined as follows: for all $q^{K}_{p-2}$ in $\mathbb{P}_{p}(K)$ and
$\widetilde{v}_{p}$ in $\widetilde{V}_{p}(K)$,

$\int_{K}q^{K}_{p-2}(\widetilde{v}_{p}-\widetilde{\Pi}^{0}_{p-2}\widetilde{v}_{p})=0.$

####  3.1.3. The primal virtual element method.

Following [3], we split the bilinear form $\widetilde{a}(\cdot,\cdot)$ into
local contributions:

$\widetilde{a}(u,v)=\sum_{K\in\mathcal{T}_{n}}\widetilde{a}^{K}(u_{|K},v_{|K}):=\sum_{K\in\mathcal{T}_{n}}(\nabla
u,\nabla v)_{0,K}\qquad\forall u,\,v\in H^{1}(\Omega).$

We approximate each local bilinear form $\widetilde{a}^{K}(\cdot,\cdot)$ with
the following discrete counterpart: for all $\widetilde{u}_{p}$ and
$\widetilde{v}_{p}$ in $\widetilde{V}_{p}(K)$,

$\widetilde{a}_{p}^{K}(\widetilde{u}_{p},\widetilde{v}_{p}):=\widetilde{a}^{K}(\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p},\widetilde{\Pi}^{\nabla}_{p}\widetilde{v}_{p})+\widetilde{S}^{K}((I-\widetilde{\Pi}^{\nabla}_{p})\widetilde{u}_{p},(I-\widetilde{\Pi}^{\nabla}_{p})\widetilde{v}_{p}),$

where the stabilizing term $\widetilde{S}^{K}(\cdot,\cdot)$ is any bilinear
form computable via the degrees of freedom and satisfying the following
bounds: there exist positive stability constants $\widetilde{\alpha}_{*}$ and
$\widetilde{\alpha}^{*}$ such that

$\widetilde{\alpha}_{*}|\widetilde{v}_{p}|^{2}_{1,K}\leq\widetilde{S}^{K}(\widetilde{v}_{p},\widetilde{v}_{p})\leq\widetilde{\alpha}^{*}|\widetilde{v}_{p}|^{2}_{1,K}\qquad\forall\widetilde{v}_{p}\in\ker(\widetilde{\Pi}^{\nabla}_{p}).$
(13)

Explicit choices of $\widetilde{S}^{K}(\cdot,\cdot)$ will be discussed in
Section 5 below. In the literature [5], it is possible to find explicit
choices of the stabilizing term such that the stability constants
$\widetilde{\alpha}_{*}$ and $\widetilde{\alpha}^{*}$ depend on $p$
algebraically only.

The global discrete bilinear form for the primal formulation is given by the
sum of the above local contributions:

$\widetilde{a}_{p}(\widetilde{u}_{p},\widetilde{v}_{p}):=\sum_{K\in\mathcal{T}_{n}}\widetilde{a}_{p}^{K}(\widetilde{u}_{p}{}_{|K},\widetilde{v}_{p}{}_{|K}).$

As for the discretization of the right-hand side, we pick

$(f_{p},\widetilde{v}_{p})_{0,\Omega}:=\sum_{K\in\mathcal{T}_{n}}(f,\widetilde{\Pi}^{0}_{p-2}\widetilde{v}_{p})_{0,K}.$

The primal virtual element method reads as follows:

$\begin{cases}\text{find }\widetilde{u}_{p}\in\widetilde{V}_{p}\text{ such
that}\\\
\widetilde{a}_{p}(\widetilde{u}_{p},\widetilde{v}_{p})=(f_{p},\widetilde{v}_{p})_{0,\Omega}\qquad\forall\widetilde{v}_{p}\in\widetilde{V}_{p}.\end{cases}$
(14)

Method (14) is well posed; see [3].

###  3.2. Mixed virtual elements

Here, we introduce stable couples of broken, mixed virtual element spaces. As
for the flux space, the local spaces are those introduced in [4], yet the
global space consists of piecewise discontinuous functions as in [8] for the
elasticity case. Furthermore, we design suitable bilinear forms appearing in
the hybridized mixed formulations of Section 4 below.

####  3.2.1. Mixed virtual element spaces.

Given $p\in\mathbb{N}_{0}$, for all $K\in\mathcal{T}_{n}$, we introduce the
local mixed space [4]

$\begin{split}\boldsymbol{\Sigma}_{p}(K):=\\{\boldsymbol{\tau}_{p}\in
H(\operatorname{div},K)\cap
H(\operatorname{rot},\Omega)&\mid\mathbf{n}_{K}\cdot\boldsymbol{\tau}_{p}{}_{|e}\in\mathbb{P}_{p}(e)\,\forall
e\in\mathcal{E}^{K},\\\\[3.0pt]
&\quad\operatorname{div}\boldsymbol{\tau}_{p}\in\mathbb{P}_{p}(K),\;\operatorname{rot}\boldsymbol{\tau}_{p}\in\mathbb{P}_{p-1}(K)\\}.\end{split}$

We have that $[\mathbb{P}_{p}(K)]^{2}\subseteq\boldsymbol{\Sigma}_{p}(K)$.
Next, consider the splitting of $[\mathbb{P}_{p}(K)]^{2}$ into

$[\mathbb{P}_{p}(K)]^{2}=\mathcal{G}_{p}(K)\oplus\mathcal{G}^{\perp}_{p}(K),\qquad\mathcal{G}_{p}(K):=\nabla\mathbb{P}_{p+1}(K),$

where $\mathcal{G}^{\perp}(K)$ is any completion of $\mathcal{G}_{p}(K)$ in
$[\mathbb{P}_{p}(K)]^{2}$. Denote the dimension of $\mathcal{G}_{p}(K)$ and
$\mathcal{G}^{\perp}_{p}(K)$ by $\mathfrak{g}_{p}$ and
$\mathfrak{g}^{\perp}_{p}$, respectively. Moreover, consider
$\\{\mathbf{g}_{\alpha}\\}_{\alpha=1}^{\mathfrak{g}_{p}}$ and
$\\{\mathbf{g}^{\perp}_{\alpha}\\}_{\alpha=1}^{\mathfrak{g}^{\perp}_{p}}$ any
bases of $\mathcal{G}_{p}(K)$ and $\mathcal{G}^{\perp}_{p}(K)$, consisting of
elements that are invariant with respect to translations and dilations; we
refer to [bricksVEMmixed] for explicit choices of such bases. The following
set of linear functionals of $\boldsymbol{\Sigma}_{p}(K)$ is a set of
unisolvent degrees of freedom [4]: given
$\boldsymbol{\tau}_{p}\in\boldsymbol{\Sigma}_{p}(K)$,

  * •

for each $e\in\mathcal{E}^{K}$, the point values of $\boldsymbol{\tau}_{p}$ at
the $p+1$ Gauß quadrature nodes of $e$;

  * •

the gradient-type moments up to order $p-1$

$\frac{1}{|K|}\int_{K}\boldsymbol{\tau}_{p}\cdot\mathbf{g}_{\alpha}\qquad\forall\alpha=1,\dots,\mathfrak{g}_{p-1};$

  * •

the “completion”-type moments up to order $p$

$\frac{1}{|K|}\int_{K}\boldsymbol{\tau}_{p}\cdot\mathbf{g}^{\perp}_{\alpha}\qquad\forall\alpha=1,\dots,\mathfrak{g}^{\perp}_{p}.$

Differently from the standard Raviart-Thomas approach of [4], we do not
construct the corresponding global space by imposing $H(\operatorname{div})$
conformity. Rather, we hybridize the space, see [8, 1], and therefore consider
the following global discontinuous mixed space:

$\boldsymbol{\Sigma}_{p}:=\\{\boldsymbol{\tau}_{p}\in[L^{2}(\Omega)]^{2}\mid\boldsymbol{\tau}_{p}{}_{|K}\in\boldsymbol{\Sigma}_{p}(K)\\}.$

Since we do not match the edge degrees of freedom, we do not impose any
boundary condition in the space. Rather, the boundary conditions will be
imposed in a mortar way through the hybridized variable.

We further introduce the space taking care of the scalar variable. To the aim,
we simply set

$V_{p}:=\mathbb{P}_{p}(\mathcal{T}_{n}).$

####  3.2.2. Projectors

Through the degrees of freedom detailed in Section 3.2.1, it is possible to
compute [4] a vector, piecewise $L^{2}$ orthogonal projector
$\boldsymbol{\Pi}^{0}_{p}:\boldsymbol{\Sigma}_{p}(K)\to[\mathbb{P}_{p}(K)]^{2}$,
which is defined as follows: for all $\mathbf{q}^{K}_{p}$ and
$\boldsymbol{\tau}_{p}\in\boldsymbol{\Sigma}_{p}(K)$,

$(\mathbf{q}^{K}_{p},\boldsymbol{\tau}_{p}-\boldsymbol{\Pi}^{0}_{p}\boldsymbol{\tau}_{p})_{0,K}=0.$

####  3.2.3. Mixed virtual element bilinear forms.

Following [4], we split the bilinear form
$a(\cdot,\cdot):=(\cdot,\cdot)_{0,\Omega}$ into local contributions:

$a(\boldsymbol{\sigma},\boldsymbol{\tau})=\sum_{K\in\mathcal{T}_{n}}a^{K}(\boldsymbol{\sigma}{}_{|K},\boldsymbol{\tau}{}_{|K}):=\sum_{K\in\mathcal{T}_{n}}(\boldsymbol{\sigma},\boldsymbol{\tau})_{0,K}\qquad\forall\boldsymbol{\sigma},\,\boldsymbol{\tau}\in[L^{2}(\Omega)]^{2}.$

We approximate each local bilinear form $a_{p}^{K}(\cdot,\cdot)$ with the
following discrete counterpart: for all $\boldsymbol{\sigma}_{p}$ and
$\boldsymbol{\tau}_{p}$ in $\boldsymbol{\Sigma}_{p}(K)$,

$a_{p}^{K}(\boldsymbol{\sigma}_{p},\boldsymbol{\tau}_{p})=a^{K}(\boldsymbol{\Pi}^{0}_{p}\boldsymbol{\sigma}_{p},\boldsymbol{\Pi}^{0}_{p}\boldsymbol{\tau}_{p})+S^{K}((\mathbf{I}-\boldsymbol{\Pi}^{0}_{p})\boldsymbol{\sigma}_{p},(\mathbf{I}-\boldsymbol{\Pi}^{0}_{p})\boldsymbol{\tau}_{p}),$
(15)

where the stabilizing term $S^{K}(\cdot,\cdot)$ is any bilinear form
computable via the degrees of freedom and satisfying the following bounds:
there exist positive stability constants $\alpha_{*}$ and $\alpha^{*}$ such
that

$\alpha_{*}\|\boldsymbol{\tau}_{p}\|^{2}_{0,K}\leq
S^{K}(\boldsymbol{\tau}_{p},\boldsymbol{\tau}_{p})\leq\alpha^{*}\|\boldsymbol{\tau}_{p}\|^{2}_{0,K}\qquad\forall\boldsymbol{\tau}_{p}\in\ker(\boldsymbol{\Pi}^{0}_{p}).$
(16)

An explicit choice of $S^{K}(\cdot,\cdot)$ will be derived and analyzed in
Appendix A below. Importantly, we shall show that the stability constants
depend on $p$ algebraically only.

The global discrete bilinear form for mixed virtual elements is given by the
sum of the above local contributions:

$a_{p}(\boldsymbol{\sigma}_{p},\boldsymbol{\tau}_{p}):=\sum_{K\in\mathcal{T}_{n}}a_{p}^{K}(\boldsymbol{\sigma}_{p}{}_{|K},\boldsymbol{\tau}_{p}{}_{|K}).$

The choice of the degrees of freedom in Section 3.2.1 allows for the explicit
computation of the local divergence of any function
$\boldsymbol{\tau}_{p}\in\boldsymbol{\Sigma}_{p}(K)$ on each element
$K\in\mathcal{T}_{n}$; see [4]. Therefore, we also define the exactly
computable bilinear form

$b(\boldsymbol{\tau}_{p},v_{p}):=\sum_{K\in\mathcal{T}_{n}}(\operatorname{div}\boldsymbol{\tau}_{p},v_{p})_{0,K}\qquad\forall\boldsymbol{\tau}_{p}\in\boldsymbol{\Sigma}_{p},\;\forall
v_{p}\in V_{p}.$

###  3.3. Hybridized virtual elements

Here, we introduce piecewise discontinuous polynomial spaces, and the related
bilinear forms and their properties, which will be instrumental in the design
of the hybridized flux reconstruction of Section 4 below.

####  3.3.1. Hybrid virtual element spaces.

Given $p\in\mathcal{T}_{n}$, we introduce the hybridization space

$\Lambda_{p}^{I}:=\mathbb{P}_{p}(\mathcal{E}_{n}^{I}).$

The following set of linear functionals is a global set of degrees of freedom
for $\Lambda_{p}^{I}$: given $\lambda_{p}^{I}\in\Lambda_{p}^{I}$,

  * •

for each $e\in\mathcal{E}_{n}^{I}$, the point values of $\lambda_{p}^{I}$ at
the $p+1$ Gauß nodes of $e$.

####  3.3.2. Hybrid virtual element bilinear forms.

We introduce a bilinear form that will be instrumental in the design of the
hybridized virtual element method. To the aim, we first define the jump
operator across an internal edge $e\in\mathcal{E}_{n}^{I}$: for any
$\boldsymbol{\tau}\in H(\operatorname{div},\mathcal{T}_{n})$,

$\llbracket\boldsymbol{\tau}\rrbracket_{e}:=\llbracket\boldsymbol{\tau}\rrbracket=\boldsymbol{\tau}_{|K^{+}}\cdot\mathbf{n}_{K^{+}}+\boldsymbol{\tau}_{|K^{-}}\cdot\mathbf{n}_{K^{-}}\quad\forall
e\in\mathcal{E}_{n}^{I},\;e\in\mathcal{E}^{K^{+}}\cap\mathcal{E}^{K^{-}}.$

We define the bilinear form
$c(\cdot,\cdot):\boldsymbol{\Sigma}_{p}\times\Lambda_{p}^{I}\to\mathbb{R}$ as

$c(\boldsymbol{\tau}_{p},\mu_{p}^{I}):=\sum_{e\in\mathcal{E}_{n}^{I}}(\llbracket\boldsymbol{\tau}_{p}\rrbracket,\mu_{p}^{I})_{0,e}\qquad\forall\boldsymbol{\tau}_{p}\in\boldsymbol{\Sigma}_{p},\,\forall\mu_{p}^{I}\in\Lambda_{p}^{I}.$

Introduce the weighted, broken divergence norm:

$\|\boldsymbol{\tau}\|_{\operatorname{div},\mathcal{T}_{n}}^{2}=\sum_{K\in\mathcal{T}_{n}}(h_{K}^{-1}\|\boldsymbol{\tau}\|^{2}_{0,K}+h_{K}\|\operatorname{div}\boldsymbol{\tau}\|^{2}_{0,K})\qquad\forall\boldsymbol{\tau}\in
H(\operatorname{div},\mathcal{T}_{n}).$

To analyze hybridized virtual elements, it is convenient to prove the
existence of a stable, broken $H(\operatorname{div})$ lifting operator
$\mathcal{L}:L^{2}(\mathcal{E}_{n}^{I})\to
H(\operatorname{div},\mathcal{T}_{n})$ such that

$\|\mathcal{L}(\mu_{p}^{I})\|_{\operatorname{div},\mathcal{T}_{n}}\leq\beta_{c}\|\mu_{p}^{I}\|_{0,\mathcal{E}_{n}^{I}},\qquad\llbracket\mathcal{L}(\mu_{p}^{I})\rrbracket_{e}=\mu_{p}^{I}\quad\forall
e\in\mathcal{E}_{n}^{I},$ (17)

where $\beta_{c}$ is a positive constant, independent of the discretization
parameters.

######  Proposition 3.1.

There exists a stable lifting operator $\mathcal{L}$ as in (17).

###### Proof.

For all $\mu_{p}^{I}\in\Lambda_{p}^{I}$, we have to construct a flux
$\mathcal{L}(\mu_{p}^{I})=\boldsymbol{\tau}_{p}\in\boldsymbol{\Sigma}_{p}$
satisfying (17). To the aim, we define $\boldsymbol{\tau}_{p}$ in each element
$K$ as $\boldsymbol{\tau}_{p}=\nabla\varphi$, where $\varphi$ is the solution
to the following local elliptic problems on each $K\in\mathcal{T}_{n}$:

$\begin{cases}\Delta\varphi=c_{K}&\text{in }K\\\\[3.0pt]
\mathbf{n}_{K}\cdot\nabla\varphi=\frac{1}{2}\mu_{p}^{I}&\text{on each
$e\in\mathcal{E}^{K}\cap\mathcal{E}_{n}^{I}$}\\\\[3.0pt]
\mathbf{n}_{K}\cdot\nabla\varphi=0&\text{on each
$e\in\mathcal{E}^{K}\cap\mathcal{E}_{n}^{D}$}\\\\[3.0pt]
\int_{K}\varphi=0,\end{cases}$

where $c_{K}\in\mathbb{R}$ is fixed so that the following compatibility
condition is valid:

$c_{K}=\frac{1}{2}|K|^{-1}\sum_{e\in\mathcal{E}^{K}\cap\mathcal{E}_{n}^{I}}(\mu_{p}^{I},1)_{0,e}.$
(18)

The above Neumann problem is well posed since the compatibility condition (18)
is valid.

The flux $\boldsymbol{\tau}_{p}=\nabla\varphi$ satisfies the following
properties:

$\begin{split}&\operatorname{div}\boldsymbol{\tau}_{p}{}_{|K}=c_{K},\qquad\|\boldsymbol{\tau}_{p}\|_{0,K}=|\varphi|_{1,K},\qquad\|\operatorname{div}\boldsymbol{\tau}_{p}\|_{0,K}=|K|^{\frac{1}{2}}|c_{K}|\qquad\forall
K\in\mathcal{T}_{n},\\\\[3.0pt]
&\llbracket\boldsymbol{\tau}_{p}\rrbracket_{e}=\mu_{p}^{I}{}_{|e}\qquad\forall
e\in\mathcal{E}_{n}^{I}.\end{split}$ (19)

Additionally, for all $K\in\mathcal{T}_{n}$, we observe that

$\|\boldsymbol{\tau}_{p}\|_{0,K}^{2}=\int_{K}\nabla\varphi\cdot\nabla\varphi=-\int_{K}(\Delta\varphi)\varphi+\int_{\partial
K}(\mathbf{n}_{K}\cdot\nabla\varphi)\varphi.$

The first term on the right-hand side is zero, since $\Delta\varphi$ is
constant and $\varphi$ has zero average over $K$. Thus, the trace and the
Poincaré inequalities yield

$\|\boldsymbol{\tau}_{p}\|_{0,K}^{2}\leq\|\mathbf{n}_{K}\cdot\nabla\varphi\|_{0,\partial
K\setminus\partial\Omega}\|\varphi\|_{0,\partial
K}\lesssim\|\mu_{p}^{I}\|_{0,\partial
K\setminus\partial\Omega}h_{K}^{\frac{1}{2}}|\varphi|_{1,K}\overset{\eqref{properties:tauboldn:inf-
sup}}{=}\|\mu_{p}^{I}\|_{0,\partial
K\setminus\partial\Omega}h_{K}^{\frac{1}{2}}\|\boldsymbol{\tau}_{p}\|_{0,K}.$

Furthermore, we have

$\|\operatorname{div}\boldsymbol{\tau}_{p}\|_{0,K}\overset{\eqref{compatibility:condition}}{=}|K|^{\frac{1}{2}}|c_{K}|\overset{\eqref{properties:tauboldn:inf-
sup}}{\lesssim}|K|^{-\frac{1}{2}}\|\mu_{p}^{I}\|_{0,\partial
K\setminus\partial\Omega}|e|^{\frac{1}{2}}\lesssim
h_{K}^{-\frac{1}{2}}\|\mu_{p}^{I}\|_{0,\partial K\setminus\partial\Omega}.$

Collecting the last two bounds and summing over all the elements yield the
assertion. ∎

##  4\. The hybridized virtual element hypercircle and local flux
reconstruction

Here, we introduce the global and local flux reconstruction in hybridized
virtual elements. The two formulations allow us to propose computable error
estimators that are reliable and efficient. The global and local versions are
detailed in Sections 4.1 and 4.2, respectively.

###  4.1. The hybridized global flux reconstruction with the virtual element
method

Here, we introduce a hybridized mixed VEM mimicking its finite element
counterpart (8).

Set

$r^{K}:=f+\Delta\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p}\in\mathbb{P}_{p-1}(K)\quad\forall
K\in\mathcal{T}_{n},\qquad
r^{e}:=\llbracket\nabla\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p}\rrbracket\in\mathbb{P}_{p-1}(e)\quad\forall
e\in\mathcal{E}_{n}^{I},$ (20)

where $\widetilde{u}_{p}$ solves the primal VEM (14).

We consider the method

$\begin{cases}\text{find
}(\boldsymbol{\sigma}_{p-1}^{\Delta},u_{p-1},\lambda_{p-1}^{I})\in\boldsymbol{\Sigma}_{p-1}\times
V_{p-1}\times\Lambda_{p-1}^{I}\text{ such that}\\\\[3.0pt]
a_{p-1}(\boldsymbol{\sigma}_{p-1}^{\Delta},\boldsymbol{\tau}_{p-1})+b(\boldsymbol{\tau}_{p-1},u_{p-1})+c(\boldsymbol{\tau}_{p-1},\lambda_{p-1}^{I})=0&\forall\boldsymbol{\tau}_{p-1}\in\boldsymbol{\Sigma}_{p-1}\\\\[3.0pt]
b(\boldsymbol{\sigma}_{p-1}^{\Delta},v_{p-1})=(r^{K},v_{p-1})_{0,\Omega}&\forall
v_{p-1}\in V_{p-1}\\\\[3.0pt]
c(\boldsymbol{\sigma}_{p-1}^{\Delta},\mu_{p-1}^{I})=(r^{e},\mu_{p-1}^{I})_{0,\mathcal{E}_{n}^{I}}&\forall\mu_{p-1}^{I}\in\Lambda_{p-1}^{I},\end{cases}$
(21)

where the discrete spaces and bilinear form are those introduced in Section 3.

Problem (21) is well posed. To see this, we observe that Proposition 3.1
guarantees the existence of a stable lifting operator, which allows us to
rewrite problem (21) with zero right-hand side in the third equation.

This implies that the flux and primal solutions are in fact solutions to a
standard mixed problem, which is well posed. The bilinear form
$a_{p}(\cdot,\cdot)$ is coercive on $\boldsymbol{\Sigma}_{p}$. Moreover, the
bilinear form $b(\cdot,\cdot)$ satisfies the following discrete inf-sup
condition:

$\inf_{v_{p}\in
V_{p}}\sup_{\boldsymbol{\tau}_{p}\in\boldsymbol{\Sigma}_{p}\cap
H(\operatorname{div},\Omega)}\frac{b(\boldsymbol{\tau}_{p},v_{p})}{(\sum_{K\in\mathcal{T}_{n}}\|\boldsymbol{\tau}_{p}\|_{\operatorname{div},K}^{2})^{\frac{1}{2}}\|v_{p}\|_{0,\Omega}}\geq\beta_{b}\qquad\forall
p\in\mathbb{N}_{0},$ (22)

where $\beta_{b}$ is a positive constant independent of $h$ and $p$. The
$p$-independence of $\beta_{b}$ is proven in [7, Theorem $2.8$].

Therefore, there exist unique $\boldsymbol{\sigma}_{p-1}^{\Delta}$ and
$u_{p-1}$ solutions to (21), with continuous dependence on the data. The
existence, uniqueness, and continuous dependence on the data of the hybridized
solution $\lambda_{p-1}^{I}$ are a consequence of the first line of the
method, picking a test function $\boldsymbol{\tau}_{p-1}$ such that

$c(\boldsymbol{\tau}_{p-1},\lambda_{p-1}^{I})=c(\lambda_{p-1}^{I},\lambda_{p-1}^{I})$

the $L^{2}$ skeleton coercivity of the bilinear form $c(\cdot,\cdot)$, and the
stability property of the lifting operator $\mathcal{L}$ in (17).

Eventually, we introduce a global error estimator based on the global flux
reconstruction:

$\eta^{2}:=\|\boldsymbol{\sigma}_{p-1}^{\Delta}\|^{2}_{a_{p-1}}+\sum_{K\in\mathcal{T}_{n}}\widetilde{S}^{K}((I-\widetilde{\Pi}^{\nabla}_{p})\widetilde{u}_{p},(I-\widetilde{\Pi}^{\nabla}_{p})\widetilde{u}_{p}),$
(23)

which we shall analyze in Section 4.2 below.

###  4.2. The hybridized global flux reconstruction: reliability and
efficiency

Here, we prove the reliability and efficiency of the error estimator $\eta$
defined in (23). We begin by proving its reliability.

######  Theorem 4.1.

The error estimator $\eta$ defined in (23) is reliable, i.e., given
$\widetilde{u}$ and $\widetilde{u}_{p}$ the primal solutions to (2) and (14),
and $\boldsymbol{\sigma}_{p-1}^{\Delta}$ the solution to (21), the following
bound is valid:

$|\widetilde{u}-\widetilde{u}_{p}|_{1,K}^{2}\leq(1+\alpha_{*}^{-1})\|\boldsymbol{\sigma}_{p-1}^{\Delta}\|^{2}_{a_{p-1}}+\widetilde{\alpha}_{*}^{-1}\sum_{K\in\mathcal{T}_{n}}\widetilde{S}^{K}((I-\widetilde{\Pi}^{\nabla}_{p})\widetilde{u}_{p},(I-\widetilde{\Pi}^{\nabla}_{p})\widetilde{u}_{p}),$

where $\widetilde{\alpha}_{*}$ and $\alpha_{*}$ are stability constants
appearing in (13) and (16).

###### Proof.

For all $\widetilde{v}\in\widetilde{V}$, we have

$\begin{split}&(\nabla(\widetilde{u}-\widetilde{u}_{p}),\nabla\widetilde{v})_{0,\Omega}=-(\Delta\widetilde{u},\widetilde{v})_{0,\Omega}-(\nabla\widetilde{u}_{p},\nabla\widetilde{v})_{0,\Omega}\\\\[3.0pt]
&\quad=\sum_{K\in\mathcal{T}_{n}}[(f+\Delta\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p},\widetilde{v})_{0,K}-(\Delta\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p},\widetilde{v})_{0,K}]-(\nabla\widetilde{u}_{p},\nabla\widetilde{v})_{0,\Omega}\\\\[3.0pt]
&\quad=\sum_{K\in\mathcal{T}_{n}}(r^{K},\widetilde{v})_{0,K}-\sum_{e\in\mathcal{E}_{n}^{I}}(r^{e},\widetilde{v})_{0,e}+\sum_{K\in\mathcal{T}_{n}}(\nabla(\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p}-\widetilde{u}_{p}),\nabla\widetilde{v})_{0,K}.\end{split}$

From the definition of the virtual element spaces and the global hybridized
problem (21), we have that

$\operatorname{div}\boldsymbol{\sigma}_{p-1}^{\Delta}{}_{|K}=r^{K}\quad\forall
K\in\mathcal{T}_{n},\qquad\llbracket\boldsymbol{\sigma}_{p-1}^{\Delta}\rrbracket{}_{e}=r^{e}\quad\forall
e\in\mathcal{E}_{n}^{I}.$

Therefore, we can write

$\begin{split}&(\nabla(\widetilde{u}-\widetilde{u}_{p}),\nabla\widetilde{v})_{0,\Omega}=\sum_{K\in\mathcal{T}_{n}}\left[-(\boldsymbol{\sigma}_{p-1}^{\Delta},\nabla\widetilde{v})_{0,K}+(\nabla(\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p}-\widetilde{u}_{p}),\nabla\widetilde{v})_{0,K}\right].\end{split}$

We focus on the local contributions of the second term on the right-hand side:

$\begin{split}(\nabla(\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p}-\widetilde{u}_{p}),\nabla\widetilde{v})_{0,K}&\leq|\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p}-\widetilde{u}_{p}|_{1,K}|\widetilde{v}|_{1,K}\\\\[3.0pt]
&\overset{\eqref{bounds:mixed-
stab:p}}{\leq}\widetilde{\alpha}_{*}^{-1}\widetilde{S}^{K}((I-\widetilde{\Pi}^{\nabla}_{p})\widetilde{u}_{p},(I-\widetilde{\Pi}^{\nabla}_{p})\widetilde{u}_{p})^{\frac{1}{2}}|\widetilde{v}|_{1,K}.\end{split}$

Summing over all the elements and collecting the above bounds yield the
assertion picking $\widetilde{v}=\widetilde{u}-\widetilde{u}_{p}$. ∎

Next, we investigate the efficiency of the error estimator $\eta$ defined in
(23), which is proven up to the presence of a higher order converging term.

######  Theorem 4.2.

The error estimator $\eta$ defined in (23) is efficient, i.e., given
$\widetilde{u}$ and $\widetilde{u}_{p}$ the primal solutions to (2) and (14),
and $\boldsymbol{\sigma}_{p-1}^{\Delta}$ the solution to (21), the following
bound is valid:

$\|\boldsymbol{\sigma}_{p-1}^{\Delta}\|_{a_{p-1}}^{2}\leq(1+\alpha^{*})\left[|\widetilde{u}-\widetilde{u}_{p}|_{1,\mathcal{T}_{n}}^{2}+|\widetilde{u}-\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}|_{1,\mathcal{T}_{n}}^{2}\right],$
(24)

where $\alpha^{*}$ is a stability constant appearing in (16).

###### Proof.

Define the discrete kernel

$\begin{split}\mathbf{Z}_{p-1}:=\\{\boldsymbol{\tau}_{p-1}\in\boldsymbol{\Sigma}_{p-1}\mid&b(\boldsymbol{\tau}_{p-1},v_{p-1})=(r^{K},v_{p-1})_{0,\Omega}\quad\forall
v_{p-1}\in V_{p-1},\\\
&c(\boldsymbol{\tau}_{p-1},\mu_{p-1}^{I})=(r^{e},\mu_{p-1}^{I})_{0,\mathcal{E}_{n}^{I}}\quad\forall\mu_{p-1}^{I}\in\Lambda_{p-1}^{I}\\}\end{split}$
(25)

and

$\begin{split}\mathbf{Z}:=\\{\boldsymbol{\tau}\in
H(\operatorname{div},\mathcal{T}_{n})\cap
H(\operatorname{rot},\mathcal{T}_{n})\mid&\operatorname{div}\boldsymbol{\tau}{}_{|K}=r^{K}\quad\forall
K\in\mathcal{E}^{K},\\\\[3.0pt]
&c(\boldsymbol{\tau},\mu_{p-1}^{I})=(r^{e},\mu_{p-1}^{I})_{0,\mathcal{E}_{n}^{I}}\quad\forall\mu_{p-1}^{I}\in\Lambda_{p-1}^{I}\\}.\end{split}$

A consequence of standard saddle point problem theory is that the solution
$\boldsymbol{\sigma}_{p-1}^{\Delta}$ to (21) satisfies

$\|\boldsymbol{\sigma}_{p-1}^{\Delta}\|_{a_{p-1}}=\inf_{\boldsymbol{\tau}_{p-1}\in\mathbf{Z}_{p-1}}\|\boldsymbol{\tau}_{p-1}\|_{a_{p-1}}\leq\|\boldsymbol{\tau}^{\Delta}_{p-1}\|_{a_{p-1}}\qquad\forall\boldsymbol{\tau}^{\Delta}_{p-1}\in\mathbf{Z}_{p-1}.$

Therefore, it suffices to show that there exists
$\boldsymbol{\tau}^{\Delta}_{p-1}\in\mathbf{Z}_{p-1}$ satisfying (24). We
define $\boldsymbol{\tau}^{\Delta}_{p-1}$ as the interpolant of a function
$\boldsymbol{\tau}^{\Delta}\in\mathbf{Z}$ in $\mathbf{Z}_{p-1}$ that we now
design.

On each $K\in\mathcal{T}_{n}$, we set $\boldsymbol{\tau}^{\Delta}=\nabla w$,
where $w$ solves

$\begin{cases}\Delta
w=r^{K}(=f+\Delta\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p})&\text{in
}K\\\\[3.0pt] \mathbf{n}_{K}\cdot\nabla
w=\mathbf{n}_{K}\cdot(\nabla\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p}-\nabla\widetilde{u})&\forall
e\in\mathcal{E}^{K}\\\\[3.0pt] \int_{K}w=0.\end{cases}$

The above local problems are well posed: the compatibility condition between
the right-hand side and the Neumann data follows from an integration by parts.
It is immediate to check that $\boldsymbol{\tau}^{\Delta}\in\mathbf{Z}$:
$\widetilde{u}$ has continuous Neumann traces across internal edges.

Recalling that the discrete bilinear form $a_{p-1}(\cdot,\cdot)$ is computed
only through the degrees of freedom of functions in mixed virtual element
spaces, we also have

$\|\boldsymbol{\tau}^{\Delta}_{p-1}\|_{a_{p-1}}=\|\boldsymbol{\tau}^{\Delta}\|_{a_{p-1}}.$

This, Remark 1 below, the fact that
$\operatorname{div}\boldsymbol{\tau}^{\Delta}{}_{|K}$ belongs to
$\mathbb{P}_{p-1}(K)$ for all $K\in\mathcal{T}_{n}$, and (16), imply

$\|\boldsymbol{\tau}^{\Delta}_{p-1}\|_{a_{p-1}}^{2}=\|\boldsymbol{\tau}^{\Delta}\|_{a_{p-1}}^{2}\leq(1+\alpha^{*})\|\boldsymbol{\tau}^{\Delta}\|^{2}_{0,\Omega}.$

Furthermore, for all $K\in\mathcal{T}_{n}$, we have

$\begin{split}\|\boldsymbol{\tau}^{\Delta}\|^{2}_{0,K}&=(\boldsymbol{\tau}^{\Delta},\nabla
w)_{0,K}=-(\operatorname{div}\boldsymbol{\tau}^{\Delta},w)_{0,K}+(\mathbf{n}_{K}\cdot\boldsymbol{\tau}^{\Delta},w)_{0,\partial
K}\\\\[3.0pt]
&=-(f+\Delta\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p},w)_{0,K}+(\mathbf{n}_{K}\cdot(\nabla\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p}-\nabla\widetilde{u}),w)_{0,\partial
K}\\\\[3.0pt]
&=(\nabla(\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p}-\widetilde{u}),\nabla
w)_{0,K}\leq|\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p}-\widetilde{u}|_{1,K}\|\boldsymbol{\tau}^{\Delta}\|_{0,K}.\end{split}$

After collecting the above estimates and summing over all the elements, we can
write

$\|\boldsymbol{\sigma}_{p-1}^{\Delta}\|_{a_{p-1}}^{2}\leq(1+\alpha^{*})|\widetilde{u}-\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p}|_{1,\mathcal{T}_{n}}^{2}\leq(1+\alpha^{*})(|\widetilde{u}-\widetilde{u}_{p}|_{1,\mathcal{T}_{n}}^{2}+|\widetilde{u}-\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}|_{1,\mathcal{T}_{n}}^{2}).$

The assertion follows upon using (16). ∎

###  4.3. The hybridized localized virtual element flux reconstruction

Here, we introduce an error estimator based on local hybridized flux
reconstruction.

To this aim, with each vertex $\nu\in\mathcal{V}_{n}$, we associate the
partition of unity function $\widetilde{\varphi}_{\nu}\in\widetilde{V}_{1}$.
Each $\widetilde{\varphi}_{\nu}$ belongs to the lowest order primal virtual
element space $\widetilde{V}_{1}$. Notably, we require that
$\widetilde{\varphi}_{\nu}$ is harmonic in each element, piecewise linear on
the skeleton, and equal to $1$ at vertex $\nu$ and $0$ at all the others.

We localize the bulk and edge residuals in (5) as follows: for all
$\nu\in\mathcal{V}_{n}$,

$\begin{split}r^{K}_{\nu}&:=\widetilde{\Pi}^{\nabla}_{1}\widetilde{\varphi}_{\nu}(f+\Delta\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p})\in\mathbb{P}_{p}(K)\quad\forall
K\in\mathcal{T}_{n}(\omega_{\nu}),\\\
r^{e}_{\nu}&:=\widetilde{\varphi}_{\nu}\llbracket\nabla\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p}\rrbracket\in\mathbb{P}_{p}(e)\qquad\qquad\quad\forall
e\in\mathcal{E}_{n}^{I}(\omega_{\nu}).\end{split}$ (26)

Proceeding as above, on each vertex patch $\omega_{\nu}$ defined in (12), we
introduce the following local problems: if $\nu$ is an interior vertex,

$\begin{cases}\text{find
}(\boldsymbol{\sigma}_{p}^{\nu},u_{p}^{\nu},\lambda_{p}^{\nu,I})\in\boldsymbol{\Sigma}_{p}^{0}(\omega_{\nu})\times
V_{p}^{*}(\omega_{\nu})\times\Lambda_{p}^{I}(\omega_{\nu})\text{ such
that}\\\\[3.0pt]
a^{\omega_{\nu}}(\boldsymbol{\sigma}_{p}^{\nu},\boldsymbol{\tau}_{p}^{\nu})+b^{\omega_{\nu}}(\boldsymbol{\tau}_{p}^{\nu},u_{p}^{\nu})+c^{\omega_{\nu}}(\boldsymbol{\tau}_{p}^{\nu},\lambda_{p}^{\nu,I})=0&\forall\boldsymbol{\tau}_{p}^{\nu}\in\boldsymbol{\Sigma}_{p}^{0}(\omega_{\nu})\\\\[3.0pt]
b^{\omega_{\nu}}(\boldsymbol{\sigma}_{p}^{\nu},v_{p}^{\nu})=(r^{K}_{\nu},v_{p}^{\nu})_{0,\omega_{\nu}}&\forall
v_{p}^{\nu}\in V_{p}^{*}(\omega_{\nu})\\\\[3.0pt]
c^{\omega_{\nu}}(\boldsymbol{\sigma}_{p}^{\nu},\mu_{p}^{\nu,I})=(r^{e}_{\nu},\mu_{p}^{\nu,I})_{0,\mathcal{E}_{n}^{I}(\omega_{\nu})}&\forall\mu_{p}^{\nu,I}\in\Lambda_{p}^{I}(\omega_{\nu}),\end{cases}$
(27)

whereas, if $\nu$ is a boundary vertex,

$\begin{cases}\text{find
}(\boldsymbol{\sigma}_{p}^{\nu},u_{p}^{\nu},\lambda_{p}^{\nu,I})\in\boldsymbol{\Sigma}_{p}^{0,\partial\Omega}(\omega_{\nu})\times
V_{p}(\omega_{\nu})\times\Lambda_{p}^{I}(\omega_{\nu})\text{ such
that}\\\\[3.0pt]
a^{\omega_{\nu}}(\boldsymbol{\sigma}_{p}^{\nu},\boldsymbol{\tau}_{p}^{\nu})+b^{\omega_{\nu}}(\boldsymbol{\tau}_{p}^{\nu},u_{p}^{\nu})+c^{\omega_{\nu}}(\boldsymbol{\tau}_{p}^{\nu},\lambda_{p}^{\nu,I})=0&\forall\boldsymbol{\tau}_{p}^{\nu}\in\boldsymbol{\Sigma}_{p}^{0,\partial\Omega}(\omega_{\nu})\\\\[3.0pt]
b^{\omega_{\nu}}(\boldsymbol{\sigma}_{p}^{\nu},v_{p})=(r^{K}_{\nu},v_{p}^{\nu})_{0,\omega_{\nu}}&\forall
v_{p}^{\nu}\in V_{p}(\omega_{\nu})\\\\[3.0pt]
c^{\omega_{\nu}}(\boldsymbol{\sigma}_{p}^{\nu},\mu_{p}^{\nu,I})=(r^{e}_{\nu},\mu_{p}^{\nu,I})_{0,\mathcal{E}_{n}^{I}(\omega_{\nu})}&\forall\mu_{p}^{\nu,I}\in\Lambda_{p}^{I}(\omega_{\nu}),\end{cases}$
(28)

where the involved spaces and bilinear forms are the obvious counterparts of
those in (9) and (10), and where

$\begin{split}\boldsymbol{\Sigma}_{p}^{0}(\omega_{\nu})&:=\\{\boldsymbol{\tau}_{p}\in\boldsymbol{\Sigma}_{p}(\omega_{\nu})\mid\mathbf{n}_{e}\cdot\boldsymbol{\tau}_{p}{}_{|e}=0\quad\forall
e\subset\partial\omega_{\nu}\\},\\\\[3.0pt]
\boldsymbol{\Sigma}_{p}^{0,\partial\Omega}(\omega_{\nu})&:=\\{\boldsymbol{\tau}_{p}\in\boldsymbol{\Sigma}_{p}(\omega_{\nu})\mid\mathbf{n}_{e}\cdot\boldsymbol{\tau}_{p}{}_{|e}=0\quad\forall
e\subset\partial\omega_{\nu},\,\nu\nsubseteq e\\},\\\\[3.0pt]
V_{p}^{*}(\omega_{\nu})&:=\left\\{v_{p}\in
V_{p}(\omega_{\nu})\Big{|}\int_{\omega_{\nu}}v_{p}=0\right\\}.\end{split}$

The well posedness of the above local problems follows as for the global
problem (21).

We introduce the following error estimator:

$\eta_{\operatorname{loc}}^{2}:=\sum_{K\in\mathcal{T}_{n}}\|\boldsymbol{\sigma}_{p}^{\Delta}\|^{2}_{a_{p}^{K}}+\sum_{K\in\mathcal{T}_{n}}\widetilde{S}^{K}((I-\widetilde{\Pi}^{\nabla}_{p})\widetilde{u}_{p},(I-\widetilde{\Pi}^{\nabla}_{p})\widetilde{u}_{p}),$
(29)

where

$\boldsymbol{\sigma}_{p}^{\Delta}:=\sum_{\nu\in\mathcal{V}_{n}}\boldsymbol{\sigma}_{p}^{\nu}.$
(30)

###  4.4. The hybridized localized virtual element flux reconstruction:
reliability and efficiency

Here, we prove the reliability and efficiency of the error estimator
$\eta_{\operatorname{loc}}$ in (29). We begin by proving the reliability.

######  Theorem 4.3.

The error estimator $\eta_{\operatorname{loc}}$ in (29) is reliable, i.e.,
given $\widetilde{u}$ and $\widetilde{u}_{p}$ the primal solutions to (2) and
(14), and $\boldsymbol{\sigma}_{p}^{\Delta}$ as in (30) the following bound is
valid:

$|\widetilde{u}-\widetilde{u}_{p}|_{1,\Omega}^{2}\leq(1+\alpha_{*}^{-1})\|\boldsymbol{\sigma}_{p}^{\Delta}\|^{2}_{a_{p}}+\widetilde{\alpha}_{*}^{-1}\sum_{K\in\mathcal{T}_{n}}\widetilde{S}^{K}((I-\widetilde{\Pi}^{\nabla}_{p})\widetilde{u}_{p},(I-\widetilde{\Pi}^{\nabla}_{p})\widetilde{u}_{p}),$

where $\widetilde{\alpha}_{*}$ and $\alpha_{*}$ are stability constants
appearing in (13) and (16).

###### Proof.

Proceeding as in the proof of Theorem 4.1, we further use the property of the
partition of unity functions, i.e,

$\sum_{\nu\in\mathcal{V}^{K}}\widetilde{\Pi}^{\nabla}_{1}\widetilde{\varphi}_{\nu}=1.$

For all $\widetilde{v}\in\widetilde{V}$, we write

$\begin{split}&(\nabla(\widetilde{u}-\widetilde{u}_{p}),\nabla\widetilde{v})_{0,\Omega}=\sum_{K\in\mathcal{T}_{n}}[(f+\Delta\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p},\widetilde{v})_{0,K}-(\Delta\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p},\widetilde{v})_{0,K}]-(\nabla\widetilde{u}_{p},\nabla\widetilde{v})_{0,\Omega}\\\\[3.0pt]
&=\sum_{\nu\in\mathcal{V}_{n}}\left[\sum_{K\in\mathcal{T}_{n}}(r^{K}_{\nu},\widetilde{v})_{0,K}-\sum_{e\in\mathcal{E}_{n}^{I}}(r^{e}_{\nu},\widetilde{v})_{0,e}\right]+\sum_{K}(\nabla(\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p}-\widetilde{u}_{p}),\nabla\widetilde{v})_{0,\Omega}.\end{split}$

Recalling the local problems (27) and (28), and using again the properties of
the partition of unity, we arrive at

$\begin{split}(\nabla(\widetilde{u}-\widetilde{u}_{p}),\nabla\widetilde{v})_{0,\Omega}&=\sum_{K\in\mathcal{T}_{n}}\sum_{\nu\in\mathcal{V}^{K}}(\boldsymbol{\sigma}_{p}^{\nu},\nabla\widetilde{v})_{0,K}+\sum_{K\in\mathcal{T}_{n}}(\nabla(\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p}-\widetilde{u}_{p}),\nabla\widetilde{v})_{0,K}\\\\[3.0pt]
&=\sum_{K\in\mathcal{T}_{n}}(\boldsymbol{\sigma}_{p}^{\Delta},\nabla\widetilde{v})_{0,K}+\sum_{K\in\mathcal{T}_{n}}(\nabla(\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p}-\widetilde{u}_{p}),\nabla\widetilde{v})_{0,K}.\end{split}$

The assertion follows picking $\widetilde{v}=\widetilde{u}-\widetilde{u}_{p}$,
and using (13) and (16). ∎

Next, we investigate the efficiency of the error estimator
$\eta_{\operatorname{loc}}$ in (29). We point out that the efficiency is
proven up to the presence of a higher order converging term.

######  Theorem 4.4.

The error estimator $\eta_{\operatorname{loc}}$ in (29) is efficient, i.e.,
given $\widetilde{u}$ and $\widetilde{u}_{p}$ the primal solutions to (2) and
(14), and $\boldsymbol{\sigma}_{p}^{\Delta}$ as in (30), the following bound
is valid:

$\sum_{\nu\in\mathcal{V}_{n}}\|\boldsymbol{\sigma}_{p}^{\Delta}\|_{a_{p}^{\omega_{\nu}}}^{2}\leq(1+\alpha^{*})\left[|\widetilde{u}-\widetilde{u}_{p}|_{1,\mathcal{T}_{n}}^{2}+|\widetilde{u}-\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}|_{1,\mathcal{T}_{n}}^{2}\right],$
(31)

where $\widetilde{\alpha}^{*}$ and $\alpha^{*}$ are stability constants
appearing in (13) and (16).

###### Proof.

The proof follows along the same lines as of Theorem 4.2. Moreover, we only
focus on the case of $\nu$ internal vertex. Boundary vertices yield to a
similar analysis. In overall, each local flux $\boldsymbol{\sigma}_{p}^{\nu}$
minimizes the $a_{p}^{\omega_{\nu}}$ norm amongst all functions in a discrete
kernel. Therefore, the proof boils down to find a discrete flux
$\boldsymbol{\tau}^{\Delta,\nu}_{p}$ that we can estimate from above by means
of the primal error. Such $\boldsymbol{\tau}^{\Delta,\nu}_{p}$ is the virtual
element interpolant of $\boldsymbol{\tau}^{\Delta,\nu}$, which is defined as
$\nabla w^{\nu}$, where, for all $K\in\mathcal{T}_{n}$

$\begin{cases}\Delta
w^{\nu}=r^{K}_{\nu}(=\widetilde{\Pi}^{\nabla}_{1}\widetilde{\varphi}_{\nu}(f+\Delta\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p}))&\text{in
}K\\\\[3.0pt] \mathbf{n}_{K}\cdot\nabla
w^{\nu}=\widetilde{\varphi}_{\nu}(\mathbf{n}_{K}\cdot(\nabla\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p}-\nabla\widetilde{u}))-c_{\nu}^{e,K}&\forall
e\in\mathcal{E}^{K}\\\\[3.0pt] \int_{K}w^{\nu}=0,\end{cases}$

where $c_{\nu}^{e,K}$ is zero on all edges of $K$ except the two having
endpoint $\nu$. In particular, $c_{\nu}^{e,K}$ is always zero on all the
boundary edges of $\omega_{\nu}$. The values of $c_{\nu}^{e,K}$ are to be
sought such that (_i_) the jump of the Neumann traces across $e$ are zero;
(_ii_) the above problem admits compatibility conditions and thence is well
posed.

Flicking through all the $N_{\nu}$ elements abutting $\nu$, we have 2
$N_{\nu}$ (each internal edge of $\omega_{\nu}$ is shared by two elements)
constants to fix by means of $N_{\nu}$ conditions given by the zero jumps and
$N_{\nu}$ conditions given by the fact that we want to have compatability
conditions for the above Laplace-Neumann problem.

Having a closer look, simple computations give that the compatibility
conditions are as follows: for all $K\in\mathcal{T}_{n}(\omega_{\nu})$ and
given $e_{1}$ and $e_{2}$ the two edges in $\mathcal{E}^{K}$ with endpoint
$\nu$, we have that

$\begin{split}\sum_{j=1}^{2}|e_{j}|c_{\nu}^{e_{j},K}&=(\widetilde{\varphi}_{\nu},\mathbf{n}_{K}\cdot(\nabla\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p}-\nabla\widetilde{u}))_{0,\partial
K}-(\widetilde{\Pi}^{\nabla}_{1}\widetilde{\varphi}_{\nu},-\Delta\widetilde{u}+\Delta\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p})_{0,K}\\\
&=(\widetilde{\varphi}_{\nu},\mathbf{n}_{K}\cdot(\nabla\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p}-\nabla\widetilde{u}))_{0,\partial
K}-(\widetilde{\Pi}^{\nabla}_{1}\widetilde{\varphi}_{\nu},\mathbf{n}_{K}\cdot(\nabla\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p}-\nabla\widetilde{u}))_{0,\partial
K}\\\\[3.0pt]
&\quad-(\nabla\widetilde{\Pi}^{\nabla}_{1}\widetilde{\varphi}_{\nu},\nabla(\widetilde{u}-\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p}))_{0,K}.\end{split}$

Importantly, the properties of the partition of unity entails that

$\sum_{\nu\in\mathcal{V}^{K}}\sum_{j=1}^{2}|e_{j}|c_{\nu}^{e_{j}}=0.$ (32)

To fully fix the $c_{\nu}^{e,K}$ constants, we simply impose zero jump
conditions across edges of the normal components of the gradients.

Further, observe that $\boldsymbol{\tau}^{\Delta,\nu}_{p}$ belongs to
$\mathbf{Z}_{p}(\omega_{\nu})$, where $\mathbf{Z}_{p}(\omega_{\nu})$ is the
obvious local counterpart of (25).

We are left to show the upper bound on the sum of the
$\|\boldsymbol{\tau}^{\Delta,\nu}\|_{a_{p}}$ on an element
$K\in\mathcal{T}_{n}$. For each $\nu\in\mathcal{V}^{K}$, we have

$\begin{split}\|\boldsymbol{\tau}^{\Delta,\nu}\|^{2}_{0,K}&=(\boldsymbol{\tau}^{\Delta},\nabla
w^{\nu})_{0,K}=-(\operatorname{div}\boldsymbol{\tau}^{\Delta,\nu},w^{\nu})_{0,K}+(\mathbf{n}_{K}\cdot\boldsymbol{\tau}^{\Delta,\nu},w^{\nu})_{0,\partial
K}\\\\[3.0pt]
&=-(\widetilde{\Pi}^{\nabla}_{1}\widetilde{\varphi}_{\nu}(f+\Delta\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p}),w^{\nu})_{0,K}+(\widetilde{\varphi}_{\nu}(\mathbf{n}_{K}\cdot(\nabla\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p}-\nabla\widetilde{u}))-c_{\nu},w^{\nu})_{0,\partial
K}.\\\\[3.0pt] \end{split}$

Summing up over all the vertices of $K$ and using (32), we deduce

$\sum_{\nu\in\mathcal{V}^{K}}\|\boldsymbol{\tau}^{\Delta,\nu}\|^{2}_{0,K}=-(f+\Delta\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p},w^{\nu})_{0,K}+(\mathbf{n}_{K}\cdot(\nabla\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p}-\nabla\widetilde{u}),w^{\nu})_{0,\partial
K}.$

To conclude, it suffices to proceed exactly as in the proof of Theorem 4.2. ∎

##  5\. Numerical results

Here, we numerically validate the reliability and efficiency of the error
estimators $\eta$ and $\eta_{\operatorname{loc}}$ in (23) and (29).

Since the primal formulation error
$|\widetilde{u}-\widetilde{u}_{p}|_{1,\Omega}$ is not computable, we rather
consider the projected error

$|\widetilde{u}-\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p}|_{1,\mathcal{T}_{n}}.$
(33)

Introduce the L-shaped domain

$\Omega:=(-1,1)^{2}\setminus(-1,0]\times(-1,0].$

Given $(r,\theta)$ the polar coordinates at $(0,0)$, we consider the test case

$u(r,\theta)=r^{\frac{2}{3}}\sin\left(\frac{2}{3}\theta\right)\quad\text{in
}\Omega.$ (34)

The remainder of this section is organized as follows. In Section 5.1, we
numerically investigate the efficiency and $p$-robustness of the proposed
error estimators by observing the performance of the efficiency index. After
presenting the adaptive algorithm and the concept of $hp$-adaptive mesh
refinements in Section 5.1, we compare the performance of the $h$\- and
$hp$-adaptive methods in Section 5.3.

###  5.1. Efficiency index

Introduce the global and local efficiency indices

$I^{2}:=\frac{\eta^{2}}{|\widetilde{u}-\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p}|_{1,\mathcal{T}_{n}}},\qquad\qquad
I_{\operatorname{loc}}^{2}:=\frac{\eta_{\operatorname{loc}}^{2}}{|\widetilde{u}-\widetilde{\Pi}^{\nabla}_{p}\widetilde{u}_{p}|_{1,\mathcal{T}_{n}}}.$
(35)

We run a uniform $p$-version of the virtual element method on a Cartesian mesh
of $48$ elements for the test case with exact solution $u$, compute the
corresponding global and localized efficiency indices in (35), and display
them in Figure 1.

Figure 1. Behaviour of the efficiency indices in (35) for the $p$-version of
the method. _Left panel:_ $I$. _Right panel:_ $I_{\operatorname{loc}}$.

From Figure 1, we observe that both the efficiency indices remain close to fix
values: the global one approximately to $2.3$; the local one approximately to
$1.3$. This confirms the $p$-robustness of the proposed error estimator. Even
more interestingly, the stability constants appearing in the efficiency
estimates (24) and (31) seem to be $p$-independent in practice.

###  5.2. The adaptive algorithm and $hp$-adaptive mesh refinements

A standard adaptive algorithm reads

$\textbf{SOLVE}\quad\Longrightarrow\quad\textbf{ESTIMATE}\quad\Longrightarrow\quad\textbf{MARK}\quad\Longrightarrow\quad\textbf{REFINE.}$

Since the solving and estimation parts are clear, we provide few details for
the marking and refining steps.

As for the marking strategy, we follow [11]: we compute the average of the
local error estimators

$\eta_{\text{average}}:=\frac{1}{\operatorname{card}(\mathcal{T}_{n})}\sum_{K\in\mathcal{T}_{n}}\eta_{K},$

where $\eta_{K}$ denotes any local error estimator on the element $K$, and
mark the elements whose error estimator is larger than
$\theta\eta_{\text{average}}$ for a given positive parameter $\theta$. In the
numerical experiments below, we pick $\theta=0.75$.

Further, we follow the $hp$-refinement strategy of Melenk and Wohlmuth [11] to
choose whether performing $h$\- or $p$-refinements. We refer to [7, Section
4], for more details, notably on the choice of refining parameters, as well as
on how the elements are refined.

###  5.3. The $h$\- and $hp$-adaptive algorithm

Here, based on the adaptive algorithm discussed in Section 5.2, we investigate
the performance of the $h$\- and $hp$-versions of the VEM upon using the
corresponding mesh refinements. We adapt meshes and spaces using both the
global (23) and localized (29) error estimators.

In Figure 2, we plot the computable error (33) for the $h$-version with $p=1$,
$2$, and $3$, as well as the $hp$-version of the method for the test case $u$.
As a starting mesh, we employ a coarse uniform Cartesian mesh consisting of
$12$ squares.

Figure 2. The computable error (33) for the $h$-version with $p=1$, $2$, and
$3$, as well as the $hp$-version of the method for the test case with exact
solution in (34). The starting mesh is a uniform Cartesian mesh made of $12$
squares. $h$\- and $hp$-refinements using (_left panel:_) the global error
estimator (23); (_right panel:_) the local error estimator (29).

For both versions, we observe the expected rate of convergence: algebraic
convergence for the $h$-version; exponential convergence for the $hp$-version.

##  6\. Conclusions

We designed two reliable and efficient error estimators for the virtual
element method based on global and local flux reconstruction. The
corresponding mixed problem are hybridized. This allows us to overcome the
limitation of our previous approach hinging upon employing a standard mixed
formulation [7]. The theoretical results are validated by numerical
experiments. As _per se_ interesting, pivot results, we highlight the
existence of a stabilization term for mixed virtual elements, with stability
constants having explicit dependence on the degree of accuracy of the method,
and the well posedness of the hybridized mixed VEM for diffusion problems.

##  7\. Acknowledgements

L. M. acknowledges support from the Austrian Science Fund (FWF) project
P33477.

##  Appendix A Stabilizations for mixed virtual elements with explicit bounds
in terms of the degree of accuracy

Here, we introduce an explicit stabilization $S^{K}(\cdot,\cdot)$ as in (15)
and prove the two bounds in (16) with explicit dependence in terms of the
degree of accuracy $p$. To the aim, we extend the strategies undertaken in
[2].

We define the stabilizing bilinear form

$S^{K}(\boldsymbol{\sigma}_{p},\boldsymbol{\tau}_{p}):=(\mathbf{n}_{K}\\!\cdot\\!\boldsymbol{\sigma}_{p},\mathbf{n}_{K}\\!\cdot\\!\boldsymbol{\tau}_{p})_{-\frac{1}{2},\partial
K}+h_{K}^{2}(\operatorname{div}\boldsymbol{\sigma}_{p},\operatorname{div}\boldsymbol{\tau}_{p})_{0,K}\quad\forall\boldsymbol{\sigma}_{p},\boldsymbol{\tau}_{p}\in\ker(\boldsymbol{\Pi}^{0}_{p})\cap\boldsymbol{\Sigma}_{p}(K).$
(36)

######  Proposition A.1.

For all $K\in\mathcal{T}_{n}$, the stabilization $S^{K}(\cdot,\cdot)$ defined
in (36) satisfies the two bounds in (16) with

$\alpha_{*}\gtrsim p^{-8},\qquad\alpha^{*}\lesssim p^{4}.$ (37)

###### Proof.

Without loss of generality, we assume that $h_{K}=1$. The general statement
follows from a scaling argument.

Given
$\boldsymbol{\tau}_{p}\in\ker(\boldsymbol{\Pi}^{0}_{p})\cap\boldsymbol{\Sigma}_{p}(K)$,
we consider the following Helmholtz-type decomposition of
$\boldsymbol{\tau}_{p}$:

$\boldsymbol{\tau}_{p}:=\nabla\psi+\operatorname{\bf{curl}}\rho.$ (38)

We require the two functions $\psi$ and $\rho$ to satisfy

$\begin{cases}\Delta\psi=\operatorname{div}\boldsymbol{\tau}_{p}&\text{in
}K\\\
\mathbf{n}_{K}\cdot\nabla\psi=\mathbf{n}_{K}\cdot\boldsymbol{\tau}_{p}&\text{on
}\partial K\\\ \int_{K}\psi=0\end{cases}$ (39)

and

$\begin{cases}(\Delta\rho=)\operatorname{rot}\operatorname{\bf{curl}}\rho=\operatorname{rot}\boldsymbol{\tau}_{p}&\text{in
}K\\\ \rho=0&\text{on }\partial K.\end{cases}$ (40)

It is apparent the $L^{2}(K)$ orthogonality of $\nabla\psi$ and
$\operatorname{\bf{curl}}\rho$, whence the following identity follows:

$\|\boldsymbol{\tau}_{p}\|^{2}_{0,K}=\|\operatorname{\bf{curl}}\rho\|^{2}_{0,K}+|\psi|^{2}_{1,K}.$
(41)

In order to derive the estimates on $\alpha_{*}$, we have to show upper bounds
for the two terms on the right-hand side of (41). We begin with that involving
$\psi$ using the trace and Poincaré inequalities:

$\begin{split}|\psi|^{2}_{1,K}&=-\int_{K}(\Delta\psi)\psi+\int_{\partial
K}(\mathbf{n}_{K}\cdot\nabla\psi)\,\psi\overset{\eqref{definition:psi}}{=}-\int_{K}\operatorname{div}\boldsymbol{\tau}_{p}\psi+\int_{\partial
K}(\mathbf{n}_{K}\cdot\nabla\psi)\,\psi\\\
&\lesssim\|\operatorname{div}\boldsymbol{\tau}_{p}\|_{0,K}|\psi|_{1,K}+\|\mathbf{n}_{K}\cdot\boldsymbol{\tau}_{p}\|_{-\frac{1}{2},\partial
K}|\psi|_{1,K}.\end{split}$ (42)

Thus, we are left with proving the bound on the term involving $\rho$. To the
aim, we write

$\|\operatorname{\bf{curl}}\rho\|^{2}_{0,K}=\int_{K}\operatorname{\bf{curl}}\rho\cdot\operatorname{\bf{curl}}\rho=\int_{K}\rho(\operatorname{rot}\operatorname{\bf{curl}}\rho)\overset{\eqref{definition:rho}}{=}\int_{K}\rho(\operatorname{rot}\boldsymbol{\tau}_{p}).$

We recall the 2D version of [10, Lemma $2.1$]: there exists a stable
polynomial inverse operator $\mathcal{K}$ of the $\operatorname{rot}$
operator. More precisely, for any $q_{p-1}\in\mathbb{P}_{p-1}(K)$ there exists
$\mathbf{q}_{p}\in[\mathbb{P}_{p}(K)]^{2}$ such that

$\mathbf{q}_{p}=\mathcal{K}q_{p-1},\qquad\operatorname{rot}\mathbf{q}_{p}=q_{p-1},\qquad\|\mathcal{K}q_{p-1}\|_{0,K}=\|\mathbf{q}_{p}\|_{0,K}\lesssim\|\operatorname{rot}\mathbf{q}_{p}\|_{0,K}=\|q_{p-1}\|_{0,K}.$
(43)

The hidden constant depends only on the shape and size of $K$. Since
$\operatorname{rot}\boldsymbol{\tau}_{p}\in\mathbb{P}_{p-1}(K)$, there exists
$\mathbf{q}_{p}\in[\mathbb{P}_{p}(K)]^{2}$ satisfying (43) with
$q_{p-1}=\operatorname{rot}\boldsymbol{\tau}_{p}$.

In light of (43), we can write

$\|\operatorname{\bf{curl}}\rho\|^{2}_{0,K}=\int_{K}\rho(\operatorname{rot}\mathbf{q}_{p}).$

We deduce

$\begin{split}\|\operatorname{\bf{curl}}\rho\|^{2}_{0,K}&\overset{\eqref{definition:rho}}{=}\int_{K}(\operatorname{\bf{curl}}\rho)\cdot\mathbf{q}_{p}\overset{\eqref{Helmholtz-
decomposition:stab}}{=}\int_{K}(\boldsymbol{\tau}_{p}-\nabla\psi)\cdot\mathbf{q}_{p}.\end{split}$

Recalling that $\boldsymbol{\tau}_{p}$ belongs to
$\ker(\boldsymbol{\Pi}^{0}_{p})$, we can write

$\begin{split}\|\operatorname{\bf{curl}}\rho\|^{2}_{0,K}&=-\int_{K}\nabla\psi\cdot\mathbf{q}_{p}\leq|\psi|_{1,K}\|\mathbf{q}_{p}\|_{0,K}\overset{\eqref{properties:Gopalakrishnan}}{\lesssim}|\psi|_{1,K}\|\operatorname{rot}\mathbf{q}_{p}\|_{0,K}.\end{split}$
(44)

We focus on the last term on the right-hand side. To this aim, recall the
following inverse estimate on polygons, which trivially follows from [13,
Propositions $3.37$ and $3.46$]: given $b_{K}$ the standard piecewise cubic
bubble over a shape-regular triangulation of $K$ satisfying
$\|b_{K}\|_{\infty,K}=1$,

$\|q_{p}\|_{0,K}\lesssim p\|b_{K}^{\frac{1}{2}}q_{p}\|_{0,K}\qquad\forall
q_{p}\in\mathbb{P}_{p}(K).$ (45)

We have

$\begin{split}\|\operatorname{rot}\mathbf{q}_{p}\|_{0,K}^{2}&\overset{\eqref{definition:rho},\eqref{properties:Gopalakrishnan}}{=}\|\operatorname{rot}\operatorname{\bf{curl}}\rho\|^{2}_{0,K}\overset{\eqref{inverse-
estimate:bubble}}{\lesssim}p^{2}\|b_{K}^{\frac{1}{2}}\operatorname{rot}\operatorname{\bf{curl}}\rho\|^{2}_{0,K}=p^{2}\int_{K}(b_{K}\operatorname{rot}\operatorname{\bf{curl}}\rho)\operatorname{rot}\operatorname{\bf{curl}}\rho\\\
&=p^{2}\int_{K}\operatorname{\bf{curl}}(b_{K}\operatorname{rot}\operatorname{\bf{curl}}\rho)\operatorname{\bf{curl}}\rho\leq
p^{2}\|\operatorname{\bf{curl}}(b_{K}\operatorname{rot}\operatorname{\bf{curl}}\rho)\|_{0,K}\|\operatorname{\bf{curl}}\rho\|_{0,K}.\end{split}$
(46)

Recall the standard polynomial inverse inequality on polygons, see, e.g., [5,
Theorem $5$]:

$\|q_{p}\|_{0,K}\lesssim p^{2}\|q_{p}\|_{-1,K}.$ (47)

The inverse inequality (47) remains valid for vector valued polynomials. We
deduce that

$\begin{split}\|\operatorname{\bf{curl}}(b_{K}\operatorname{rot}\operatorname{\bf{curl}}\rho)\|_{0,K}&\lesssim
p^{2}\sup_{\mathbf{v}\in[H^{1}_{0}(K)]^{2}}\frac{(\operatorname{\bf{curl}}(b_{K}\operatorname{rot}\operatorname{\bf{curl}}\rho),\mathbf{v})_{0,K}}{|\mathbf{v}|_{1,K}}\\\
&=p^{2}\sup_{\mathbf{v}\in[H^{1}_{0}(K)]^{2}}\frac{(b_{K}\operatorname{rot}\operatorname{\bf{curl}}\rho,\operatorname{rot}\mathbf{v})_{0,K}}{|\mathbf{v}|_{1,K}}\lesssim
p^{2}\|\operatorname{rot}\operatorname{\bf{curl}}\rho\|_{0,K}.\end{split}$
(48)

Using (46) and the above inverse estimate, we arrive at

$\|\operatorname{rot}\mathbf{q}_{p}\|_{0,K}^{2}\lesssim
p^{4}\|\operatorname{rot}\operatorname{\bf{curl}}\rho\|_{0,K}\|\operatorname{\bf{curl}}\rho\|_{0,K}=p^{4}\|\operatorname{rot}\mathbf{q}_{p}\|_{0,K}\|\operatorname{\bf{curl}}\rho\|_{0,K}.$

Inserting this bound in (44) entails

$\|\operatorname{\bf{curl}}\rho\|_{0,K}\lesssim p^{4}|\psi|_{1,K}.$

Combining the above estimate with (42) yields the desired bound on
$\alpha_{*}$.

As for the bound on $\alpha^{*}$, we refer to [7, Proposition $2.7$] for full
details. Here, we only sketch the proof. We need to estimate from above the
two terms on the right-hand side of (36). As for the boundary one, we apply a
$\operatorname{div}$ trace inequality, see, e.g., [12]:

$\|\mathbf{n}_{K}\cdot\boldsymbol{\tau}_{p}\|_{-\frac{1}{2},\partial
K}\lesssim\|\boldsymbol{\tau}_{p}\|_{0,K}+\|\operatorname{div}\boldsymbol{\tau}_{p}\|_{0,K}.$

Therefore, it suffices to estimate the divergence term. This can be done via
an inverse inequality as detailed in the proof of [7, Proposition $2.7$],
which can be proven as in (48):

$\|\operatorname{div}\boldsymbol{\tau}_{p}\|_{0,K}\lesssim
p^{2}\|\boldsymbol{\tau}_{p}\|_{0,K}.$

The assertion follows. ∎

######  Remark 1.

From the proof of Proposition A.1, it is immediate to check the bound on
$\alpha^{*}$ in (37) is valid not only for functions in
$\ker(\boldsymbol{\Pi}^{0})\cap\boldsymbol{\Sigma}_{p}(K)$, but for all
functions $\boldsymbol{\tau}\in H(\operatorname{div},K)\cap
H(\operatorname{rot},K)$ such that $\operatorname{div}\boldsymbol{\tau}$
belongs to $\mathbb{P}_{p}(K)$.

In practical computations, the boundary contribution in (36) can be
substituted by

$h_{K}(\mathbf{n}_{K}\cdot\boldsymbol{\sigma}_{p},\mathbf{n}_{K}\cdot\boldsymbol{\tau}_{p})_{0,\partial
K}.$

In this case, it is possible to check that the bounds on the stability
constants become

$\alpha_{*}\gtrsim p^{-8},\qquad\alpha^{*}\lesssim p^{6}.$

## References

  * [1] D. N. Arnold and F. Brezzi.  Mixed and nonconforming finite element methods: implementation, postprocessing and error estimates.  ESAIM Math. Model. Numer. Anal., 19(1):7–32, 1985.
  * [2] 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.
  * [3] 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.
  * [4] L. Beirão 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.
  * [5] L. Beirão da Veiga, A. Chernov, L. Mascotto, and A. Russo.  Exponential convergence of the $hp$ virtual element method with corner singularity.  Numer. Math., 138(3):581–613, 2018.
  * [6] D. Braess, V. Pillwein, and J. Schöberl.  Equilibrated residual error estimates are $p$-robust.  Comput. Methods Appl. Mech. Engrg., 198(13-14):1189–1197, 2009\.
  * [7] F. Dassi, J. Gedicke, and L. Mascotto.  Adaptive virtual elements with equilibrated fluxes.  https://arxiv.org/abs/2004.11220, 2020.
  * [8] F. Dassi, C. Lovadina, and M. Visinoni.  Hybridization of the virtual element method for linear elasticity problems.  http://arxiv.org/abs/2103.01164, 2021.
  * [9] A. Ern and M. Vohralík.  Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations.  SIAM J. Numer. Anal., 53(2):1058–1081, 2015.
  * [10] J. Gopalakrishnan and L. F. Demkowicz.  Quasioptimality of some spectral mixed methods.  J. Comput. Appl. Math., 167(1):163–182, 2004.
  * [11] J. M. Melenk and B. I. Wohlmuth.  On residual-based a posteriori error estimation in $hp$-FEM.  Adv. Comput. Math., 15(1-4):311–331, 2001.
  * [12] P. Monk.  Finite Element Methods for Maxwell’s Equations.  Oxford University Press, 2003.
  * [13] R. Verfürth.  A posteriori error estimation techniques for finite element methods.  OUP Oxford, 2013.