# Homogeneous multigrid for hybrid discretizations: application to HHO methods

Daniele A. Di Pietro  IMAG, Univ. Montpellier, CNRS, Montpellier, France
<EMAIL_ADDRESS>,  Zhaonan Dong  Inria, 2 rue Simone Iff,
75589 Paris, France, 2) CERMICS, Ecole des Ponts, 77455 Marne-la- Vallée,
France<EMAIL_ADDRESS>,  Guido Kanschat  Interdisciplinary Center for
Scientific Computing (IWR), Heidelberg University, Mathematikon, Im
Neuenheimer Feld 205, 69120 Heidelberg, Germany<EMAIL_ADDRESS>,
Pierre Matalon  CMAP, CNRS, École Polytechnique, Institut Polytechnique de
Paris, 91120 Palaiseau, France<EMAIL_ADDRESS>and  Andreas
Rupp  School of Engineering Sciences, Lappeenranta–Lahti University of
Technology, P.O. Box 20, 53851 Lappeenranta, Finland<EMAIL_ADDRESS>
(Date: September 3, 2024)

###### Abstract.

We prove the uniform convergence of the geometric multigrid V-cycle for hybrid
high-order (HHO) and other discontinuous skeletal methods. Our results
generalize previously established results for HDG methods, and our multigrid
method uses standard smoothers and local solvers that are bounded, convergent,
and consistent. We use a weak version of elliptic regularity in our proofs.
Numerical experiments confirm our theoretical results.
Keywords. Geometric multigrid, homogeneous multigrid, hybrid high-order
methods, skeleton methods.

###### 2010 Mathematics Subject Classification:

65F10, 65N30, 65N50

##  1\. Introduction

In the context of fast solvers for linear systems arising from hybrid
discretizations of second-order elliptic equations, we propose in this work a
generalized framework for the convergence analysis of geometric multigrid
methods [BVEM00, TOS01].

Hybrid discretization methods have been part of the numerical analyst’s
toolbox to solve partial differential equations since the seventies, starting
with hybridized versions of mixed methods such as the Raviart-Thomas (RT-H)
[RT77] and the Brezzi-Douglas-Marini (BDM-H) [BDM85] methods. They have gained
growing interest in recent years with the outbreak of modern schemes such as
the unifying framework of Hybridizable Discontinuous Galerkin (HDG) methods
[CGL09] and Hybrid High-Order (HHO) methods [DPEL14, DPE15, DPD20]. The list
also includes, but is not limited to, Compatible Discrete Operators (CDO)
[BE14], Mixed and Hybrid Finite Volumes (MHFV) [DE06, EGH10, DEGH10], Mimetic
Finite Differences (MFD) [BdVLM14], weak Galerkin [WY13], local discontinuous
Galerkin-hybridizable (LDG-H) [CDG08] and discontinuous Petrov–Galerkin (DPG)
[DG11] methods. Hybrid discretization methods are characterized by the
location of their degrees of freedom (DOFs), placed both within the mesh cells
and on the faces. In this configuration, the discrete scheme is built so that
the cell DOFs are only locally coupled, leaving the face DOFs in charge of the
global coupling. Algebraically, this feature enables the local elimination of
the cell unknowns from the arising linear system, resulting in a Schur
complement of reduced size, where only the face unknowns remain. The
mechanical engineering terminology refers to this elimination process as
_static condensation_ , and the resulting system is often called the
(statically) _condensed_ system. The terminology _trace_ or _skeleton_ system
is also employed. For an extended introduction to hybrid methods and
hybridization, we refer to the preface of [DPD20] and the first pages of
[CGL09].

This paper addresses the fast solution of the condensed systems arising from
hybrid discretizations of second-order elliptic equations, focusing on
multigrid methods. Various multigrid solvers and preconditioners have been
designed over the last decade, focusing primarily on HDG [CDGT14, SA17,
WMBT19, MBTS20, LRK21, LRK22a] and HHO [DPHM+21a, DPHM+21b, BDP21, DPMMR23],
but also covering specifically DPG [RC17, PD21] and RT-H/BDM-H [GT09]. The
above references focus on geometric multigrid methods, while a fully algebraic
multigrid method has been considered in [DPHM+23].

The main difficulty in designing a geometric multigrid algorithm for a trace
system resides in the fact that the DOFs are supported by the mesh skeleton,
which makes classical intergrid transfer operators designed for element-based
DOFs unsuitable. While earlier approaches [CDGT14, GT09, KW18] recast trace
functions into bulk functions in order to make use of a known efficient solver
(typically, a standard piecewise linear continuous Finite Element multigrid
solver), the more recent developments seem to converge toward the so-called
notion of _homogeneous_ multigrid, where the hybrid discretization is
conserved at every level of the mesh hierarchy and “interskeleton” transfer
operators are designed to communicate trace functions from one mesh skeleton
to the other. To this end, multiple skeletal injection operators have already
been proposed [LRK22a, DPHM+21a].

In this paper, we build upon the work of [LRK21, LRK22a, LRK22b, LWKR23] on
HDG methods to propose a generalized demonstration framework for the uniform
convergence of _homogeneous_ multigrid methods with V-cycle. The main
motivation behind this generalization is to include the HHO methods in its
application scope. Indeed, if recent multigrid solvers for HHO methods
[DPHM+21a, DPHM+21b, DPMMR23] have experimentally shown their optimal
behavior, the supporting theory is still missing. In the second part of this
paper, we then prove that the classical HHO method fulfills the assumptions of
this new framework.

The present theory does not assume full elliptic regularity of the problem,
thus allowing complex domains with re-entrant corners. However, although
modern hybrid methods natively handle polyhedral elements, it is restricted to
simplicial meshes. The multigrid method is built in a standard fashion from an
abstract injection operator: standard smoothers are used, and the restriction
operator is chosen as the adjoint of the injection operator. A classical,
symmetric V-cycle is used if the problem exhibits full elliptic regularity. If
not, then a variable V-cycle is employed, in which the number of (symmetric)
smoothing steps increases as the level decreases in the hierarchy.

The paper is organized as follows. In Section 2, we introduce notations and
the model problem. Section 3 describes the abstract framework: (i) the hybrid
method is described abstractly, in the form of approximation spaces, local
linear operators, stabilization term and bilinear form; (ii) the injection
operators used in the multigrid method and its analysis are introduced; (iii)
the properties assumed from the hybrid method ((HM1)-(HM8)) and the injection
operator ((IA1)-(IA2)) are listed. Section 4 describes the multigrid method
and asserts the associated convergence results, while Section 5 carries out
the convergence analysis. In Section 6, we verify our framework’s assumptions
for the standard HHO method, and analyze various injection operators in
Section 7. Finally, Section 8 presents the numerical experiments supporting
our theoretical results, realized with the HHO method on two and three-
dimensional test cases.

##  2\. Problem formulation and notation

Let $\Omega\subset\mathbb{R}^{d}$ be a polygonally bounded Lipschitz domain
with boundary $\partial\Omega$. We approximate the unknown $u\in
H^{1}_{0}(\Omega)$ satisfying

$-\Delta u=f\qquad\text{ in }\Omega$ (2.1)

for some right-hand side $f$. We assume the problem to be regular in the sense
that

$\|u\|_{1+\alpha}\leq c\|f\|_{\alpha-1}\qquad\text{ for some
}\alpha\in(\tfrac{1}{2},1].$ (2.2)

Here, $\|\cdot\|_{\alpha}$ for $\alpha\in\mathbb{R}_{\geq 0}$ denotes the
induced norm of the Sobolev–Slobodeckij space $H^{\alpha}(\Omega)$. By
$\|\cdot\|_{-\alpha}$ and $H^{-\alpha}(\Omega)$ we denote the dual norm and
the dual space of $H^{\alpha}_{0}(\Omega)$ with respect to the extension of
$L^{2}$ duality, respectively. We denote by $\|\cdot\|_{X}$ and
$(\cdot,\cdot)_{X}$ the $L^{2}$-norm and $L^{2}$-scalar product with respect
to the set $X\subset\Omega$, respectively. The same notation is also used for
the inner product of $[L^{2}(X)]^{d}$ (the exact meaning can be inferred from
the context) such that for all
$\boldsymbol{p},\boldsymbol{q}\in[L^{2}(X)]^{d}$,
$(p,q)_{X}\coloneq\int_{X}\boldsymbol{p}\cdot\boldsymbol{q}$, where $\cdot$
denotes the dot product. Without loss of generality, we assume that $\Omega$
has diameter 1.

######  Remark 2.1 (More general boundary conditions).

The following approaches directly transfer to non-homogeneous Dirichlet
boundary conditions since those only influence the right-hand sides of the
numerical schemes. However, this is not true (without some further arguments)
for locally changing the type of the boundary condition.

##  3\. Discontinuous skeletal methods

We will render a general framework for all discontinuous skeletal methods
covering HDG and HHO methods. To this end, we start with a family
$(\mathcal{T}_{\ell})_{\ell=0,\dots,L}$ of successively refined simplicial
meshes and their corresponding face sets
$(\mathcal{F}_{\ell})_{\ell=0,\dots,L}$. For all level $\ell$, we define
$h_{\ell}=\max_{T\in\mathcal{T}_{\ell}}h_{T}$, where $h_{T}$ denotes the
diameter of $T$. The family $(\mathcal{T}_{\ell})_{\ell=0,\dots,L}$ is regular
(see [EG21, Def. 11.2] for a precise definition), which implies that mesh
elements do not deteriorate, and all of its elements are geometrically
conforming, which excludes hanging nodes. Additionally, we assume that either
$F\subset\partial\Omega$, or $F\cap\partial\Omega$ comprises at most one
point, and that refinement does not progress too fast, i.e., there is
$c_{\text{ref}}>0$ such that

$h_{\ell}\geq c_{\text{ref}}h_{\ell-1}\qquad\forall\ell\in\\{1,\dots,L\\}.$
(3.1)

Finally, we assume that our mesh is quasi-uniform implying that there is a
constant $c_{\textup{uni}}$ such that $h_{\ell}\leq
c_{\textup{uni}}\min_{T\in\mathcal{T}_{\ell}}h_{T}$.

Next, we choose a finite-dimensional space $M_{F}\subset L^{2}(F)$ of
functions living on face $F$ for all $F\in\mathcal{F}_{\ell}$. In the
numerical scheme, this space will be the test and trial spaces for the
skeletal variable $m_{\ell}$ approximating the trace of the unknown $u$ on the
skeleton, which we denote by $\mathcal{F}_{\ell}$. Moreover, finite-
dimensional approximation spaces $V_{T}\subset L^{2}(T)$ and
$\boldsymbol{W}_{T}\subset L^{2}(T)^{d}$ are defined element-wise. The space
$V_{T}$ is the local test and trial space for the primary unknown $u$, while
the space $\boldsymbol{W}_{T}$ is the local test and trial space for the dual
unknown $\boldsymbol{q}=-\nabla u$. Global spaces are defined by concatenation
via

$\displaystyle M_{\ell}$ $\displaystyle\coloneq\left\\{m\in
L^{2}(\mathcal{F}_{\ell})\;\middle|\;\begin{array}[]{r@{\,}c@{\,}ll}m_{|F}&\in&\mathcal{P}_{p}&\forall
F\in\mathcal{F}_{\ell}\\\ m_{|F}&=&0&\forall
F\subset\partial\Omega\end{array}\right\\},$ (3.4) $\displaystyle V_{\ell}$
$\displaystyle\coloneq\bigl{\\{}v\in L^{2}(\Omega)\big{|}\;v_{|T}\in
V_{T},\forall T\in\mathcal{T}_{\ell}\bigr{\\}},$ (3.5)
$\displaystyle\boldsymbol{W}_{\ell}$
$\displaystyle\coloneq\bigl{\\{}\boldsymbol{q}\in
L^{2}(\Omega;\mathbb{R}^{d})\big{|}\;\boldsymbol{q}_{|T}\in\boldsymbol{W}_{T},\forall
T\in\mathcal{T}_{\ell}\bigr{\\}}.$ (3.6)

We define element-wise, abstract linear operators

$\displaystyle\mathcal{U}_{T}\colon M_{\ell}|_{\partial T}\to V_{T},$
$\displaystyle\mathcal{V}_{T}\colon L^{2}(T)\to V_{T},$
$\displaystyle\boldsymbol{\mathcal{Q}}_{T}\colon M_{\ell}|_{\partial
T}\to\boldsymbol{W}_{T}.$

The choice of these operators, called local solvers, influences the numerical
schemes. Global linear operators are constructed by concatenation of their
element-local analogous:

$\displaystyle\mathcal{U}_{\ell}\colon M_{\ell}\to V_{\ell},$
$\displaystyle\mathcal{V}_{\ell}\colon L^{2}(\Omega)\to V_{\ell},$
$\displaystyle\boldsymbol{\mathcal{Q}}_{\ell}\colon
M_{\ell}\to\boldsymbol{W}_{\ell}.$

One can show that $m_{\ell}$ approximates the trace of $u$ on the skeleton and
that

$u\approx\mathcal{U}_{\ell}m_{\ell}+\mathcal{V}_{\ell}f$

in the bulk of $\Omega$ if $m_{\ell}$ satisfies

$a_{\ell}(m_{\ell},\mu)=\int_{\Omega}f\,\mathcal{U}_{\ell}\mu\qquad\text{ for
all }\mu\in M_{\ell},$ (3.7)

where the elliptic and continuous bilinear form $a_{\ell}$ has the form

$a_{\ell}(m,\mu)=\int_{\Omega}\boldsymbol{\mathcal{Q}}_{\ell}m\cdot\boldsymbol{\mathcal{Q}}_{\ell}\mu+s_{\ell}(m,\mu).$
(3.8)

The symmetric positive semi-definite $s_{\ell}$ is usually referred to as
(condensed) _penalty_ or _stabilizing term_.

The choices of the spaces $M_{F}$, $V_{T}$, $\boldsymbol{W}_{T}$, the local
solvers $\mathcal{U}_{\ell}$, $\mathcal{V}_{\ell}$,
$\boldsymbol{\mathcal{Q}}_{\ell}$ and the stabilizing term $s_{\ell}$
completely define a discontinuous skeletal method.

######  Remark 3.1 (Possible choices of hybrid methods).

  * •

For the LDG-H method, we set $M_{F}=\mathcal{P}_{p}(F)$,
$V_{T}=\mathcal{P}_{p}(T)$, and $\boldsymbol{W}_{T}=\mathcal{P}^{d}_{p}(T)$.
The operators $\mathcal{U}_{T}$ and $\boldsymbol{\mathcal{Q}}_{T}$ map
$m_{\partial T}$ to the element-wise solutions $u_{T}\in V_{T}$ and
$\boldsymbol{q}_{T}\in\boldsymbol{W}_{T}$ of

$\displaystyle\int_{T}\boldsymbol{q}_{T}\cdot\boldsymbol{p}_{T}-\int_{T}u_{T}\nabla\\!\cdot\\!\boldsymbol{p}_{T}$
$\displaystyle=-\int_{\partial T}m\boldsymbol{p}_{T}\cdot\boldsymbol{\nu}$
(3.9a) $\displaystyle\int_{\partial
T}(\boldsymbol{q}_{T}\cdot\boldsymbol{\nu}+\tau_{\ell}u_{T})v_{T}-\int_{T}\boldsymbol{q}_{T}\cdot\nabla
v_{T}$ $\displaystyle=\tau_{\ell}\int_{\partial T}mv_{T}$ (3.9b)

with test functions $v_{T}\in V_{T}$ and
$\boldsymbol{p}_{T}\in\boldsymbol{W}_{T}$, and $\boldsymbol{\nu}$ the outward
normal vector to $\partial T$. In this sense $\mathcal{U}_{T}\colon
m_{\partial T}\mapsto u_{T}$, $\boldsymbol{\mathcal{Q}}_{T}:m_{\partial
T}\mapsto\boldsymbol{q}_{T}$ for all $T\in\mathcal{T}_{\ell}$.

Finally, the bilinear form $a_{\ell}$ is defined as

$a_{\ell}(m_{\ell},\mu)=\int_{\Omega}\boldsymbol{\mathcal{Q}}_{\ell}m_{\ell}\cdot\boldsymbol{\mathcal{Q}}_{\ell}\mu+\underbrace{\sum_{T\in\mathcal{T}_{\ell}}\tau_{\ell}\int_{\partial
T}(\mathcal{U}_{\ell}m_{\ell}-m_{\ell})(\mathcal{U}_{\ell}\mu-\mu)}_{=s_{\ell}(m_{\ell},\mu)}$

for some parameter $\tau_{\ell}>0$.

  * •

For the RT-H and BDM-H methods, one uses the framework of the LDG-H method
replacing the local bulk spaces by the RT-H and BDM-H spaces with
$\tau_{\ell}=0$.

  * •

For the HHO method, we set $M_{F}=\mathcal{P}_{p}(F)$,
$V_{T}=\mathcal{P}_{p}(T)$, and
$\boldsymbol{W}_{T}=\nabla\mathcal{P}_{p+1}(T)$. For $(u_{T},m_{\partial
T})\in V_{T}\times M_{\partial T}$, we define
$\boldsymbol{q}_{T}(u_{T},m_{\partial T})$ as the element-wise solution of

$\int_{T}\boldsymbol{q}_{T}(u_{T},m_{\partial
T})\cdot\boldsymbol{p}_{T}-\int_{T}u_{T}\nabla\\!\cdot\\!\boldsymbol{p}_{T}=-\int_{\partial
T}m_{\partial T}\boldsymbol{p}_{T}\cdot\boldsymbol{\nu}$ (3.10)

for all $\boldsymbol{p}_{T}\in\boldsymbol{W}_{T}$. Note that this relation is
the same as (3.9a). In HHO terminology, $\boldsymbol{q}_{T}(u_{T},m_{\partial
T})$ corresponds to $-\nabla\theta_{T}^{p+1}(u_{T},m_{\partial T})$, where
$\theta_{T}^{p+1}$ denotes the so-called local higher-order reconstruction
operator, defined by (3.10) and the closure condition

$\int_{T}\theta_{T}^{p+1}(u_{T},m_{\partial T})=\int_{\partial T}u_{T}.$
(3.11)

In a second step, given a local bilinear stabilizer
$\underline{s}_{T}((u_{T},m_{\partial T}),(v_{T},\mu))$, we introduce the
local bilinear form

$\underline{a}_{T}((u_{T},m_{\partial T}),(v_{T},\mu))\\\
=\int_{T}\boldsymbol{q}_{T}(u_{T},m_{\partial
T})\cdot\boldsymbol{q}_{T}(v_{T},\mu)+\underline{s}_{T}((u_{T},m_{\partial
T}),(v_{T},\mu)).$ (3.12)

Consider the following local problems:
(i) For $m_{\partial T}\in M_{\partial T}$, find $u_{T}^{1}\in V_{T}$ such
that

$\underline{a}_{T}((u_{T}^{1},0),(v_{T},0))=-\underline{a}_{T}((0,m_{\partial
T}),(v_{T},0))\qquad\forall v_{T}\in V_{T}.$ (3.13)

(ii) For $f\in L^{2}(\Omega)$, find $u_{T}^{2}\in V_{T}$ such that

$\underline{a}_{T}((u_{T}^{2},0),(v_{T},0))=\int_{T}fv_{T}\qquad\forall
v_{T}\in V_{T}.$ (3.14)

We define $\mathcal{U}_{T}\colon m_{\partial T}\mapsto u_{T}^{1}$ solution of
(3.13), $\mathcal{V}_{T}\colon f\mapsto u_{T}^{2}$ solution of (3.14), and we
set element-by-element

$\boldsymbol{\mathcal{Q}}_{T}\colon m_{\partial
T}\mapsto\boldsymbol{q}_{T}(\mathcal{U}_{T}m_{\partial T},m_{\partial T})$
(3.15)

Finally, the global bilinear form is defined as

$\displaystyle a_{\ell}(m_{\ell},\mu)=:$
$\displaystyle\sum_{T\in\mathcal{T}_{\ell}}\underline{a}_{T}((\mathcal{U}_{T}m_{\partial
T},m_{\partial T}),(\mathcal{U}_{T}\mu,\mu))$ (3.16) $\displaystyle=$
$\displaystyle\int_{\Omega}\boldsymbol{\mathcal{Q}}_{\ell}m_{\ell}\cdot\boldsymbol{\mathcal{Q}}_{\ell}\mu+\underbrace{\sum_{T\in\mathcal{T}_{\ell}}\underline{s}_{T}((\mathcal{U}_{\ell}m_{\ell},m_{\ell}),(\mathcal{U}_{\ell}\mu,\mu))}_{=s_{\ell}(m_{\ell},\mu)}.$

###  3.1. Operators for the multigrid method and analysis

For functions $\rho,\mu\in L^{2}(\mathcal{F}_{\ell})$, we set

$\langle\rho,\mu\rangle_{\ell}=\sum_{T\in\mathcal{T}}\frac{|T|}{|\partial
T|}\int_{\partial T}\rho\mu.$

Importantly, $\langle\cdot,\cdot\rangle_{\ell}$ defines a scalar product on
$M_{\ell}$, and the norm of this scalar product

$\|\cdot\|_{\ell}\coloneq\langle\cdot,\cdot\rangle_{\ell}^{\frac{1}{2}}$

scales like the $L^{2}$-norm in the bulk of the domain. Moreover, the scalar
product can readily be applied to a combination of bulk and skeleton
functions: For example, $\langle\mu,v\rangle_{\ell}$ can readily be evaluated
for $\mu\in M_{\ell}$ and $v\in V_{\ell}$ with the understanding that the
trace of $v$ on $\partial T$ is used inside the integral. We relate the
bilinear forms $a_{\ell}$ and $\langle\cdot,\cdot\rangle_{\ell}$ using the
operator $A_{\ell}$, which is defined via

$\langle
A_{\ell}\rho,\mu\rangle_{\ell}=a_{\ell}(\rho,\mu)\qquad\forall\rho,\mu\in
M_{\ell}.$ (3.17)

The injection operator $I_{\ell}\colon M_{\ell-1}\to M_{\ell}$ remains
abstract at this level. Its properties securing our analytical findings are
listed in Section 3.3. Possible realizations of injection operators can be
found in [LRK22a]. The multigrid operator $B_{\ell}\colon M_{\ell}\to
M_{\ell}$ for preconditioning $A_{\ell}$ will be defined in Section 4.2.

The operators $P^{\langle\cdot\rangle}_{\ell-1}$ and $P^{a}_{\ell-1}$, which
are defined via

$\displaystyle P^{\langle\cdot\rangle}_{\ell-1}$ $\displaystyle\colon
M_{\ell}\to M_{\ell-1},$ $\displaystyle\langle
P^{\langle\cdot\rangle}_{\ell-1}\rho,\mu\rangle_{\ell-1}$
$\displaystyle=\langle\rho,I_{\ell}\mu\rangle_{\ell}$
$\displaystyle\forall\mu\in M_{\ell-1},$ (3.18) $\displaystyle P^{a}_{\ell-1}$
$\displaystyle\colon M_{\ell}\to M_{\ell-1},$ $\displaystyle
a_{\ell-1}(P^{a}_{\ell-1}\rho,\mu)$ $\displaystyle=a_{\ell}(\rho,I_{\ell}\mu)$
$\displaystyle\forall\mu\in M_{\ell-1},$ (3.19)

replace the $L^{2}$ and the Ritz projections of conforming methods,
respectively. While the former (or a discrete variation of it) is relevant for
implementing multigrid methods, the latter is key to the analysis. We also
introduce the $L^{2}$-projections

$\displaystyle\Pi^{\partial}_{\ell}\colon$ $\displaystyle H^{1}_{0}(\Omega)\to
M_{\ell},$
$\displaystyle\langle\Pi^{\partial}_{\ell}v,\mu\rangle_{\ell}=\langle
v,\mu\rangle_{\ell}$ $\displaystyle\forall\mu\in M_{\ell},$
$\displaystyle\Pi^{\textup{d}}_{\ell}\colon$ $\displaystyle H^{1}(\Omega)\to
V_{\ell},$ $\displaystyle(\Pi^{\textup{d}}_{\ell}v,w)_{\Omega}=(v,w)_{\Omega}$
$\displaystyle\forall w\in V_{\ell}.$

###  3.2. Assumptions on discontinuous skeletal methods

Here and in the following, $\lesssim$ means smaller than or equal to, up to a
constant independent of the mesh size $h_{\ell}$ and the multigrid level
$\ell$. We also write $A\simeq B$ as a shortcut for “$A\lesssim B$ and
$B\lesssim A$”. We assume that the numerical hybrid method, characterized in
Section 3, satisfies the following conditions for any $\mu\in M_{\ell}$:

  * •

The trace of the bulk unknown approximates the skeletal unknown:

$\|\mathcal{U}_{\ell}\mu-\mu\|_{\ell}\lesssim h_{\ell}\|\mu\|_{a_{\ell}},$
(HM1)

where $\|\cdot\|_{a_{\ell}}$ denotes the norm induced by $a_{\ell}$ on
$M_{\ell}$, i.e.,

$\|\cdot\|_{a_{\ell}}\coloneq a_{\ell}(\cdot,\cdot)^{\frac{1}{2}}.$ (3.20)

  * •

The operators $\boldsymbol{\mathcal{Q}}_{\ell}\mu$ and $\mathcal{U}_{\ell}\mu$
are continuous:

$\|\boldsymbol{\mathcal{Q}}_{\ell}\mu\|_{\Omega}\lesssim
h^{-1}_{\ell}\|\mu\|_{\ell}\quad\text{ and
}\quad\|\mathcal{U}_{\ell}\mu\|_{\Omega}\lesssim\|\mu\|_{\ell}.$ (HM2)

  * •

The dual approximation
$\boldsymbol{\mathcal{Q}}_{\ell}\mu\sim-\nabla_{\ell}\mathcal{U}_{\ell}\mu$,
where $\nabla_{\ell}$ denotes the broken gradient. That is,

$\|\boldsymbol{\mathcal{Q}}_{\ell}\mu+\nabla_{\ell}\mathcal{U}_{\ell}\mu\|_{\Omega}\lesssim
h_{\ell}^{-1}\|\mathcal{U}_{\ell}\mu-\mu\|_{\ell}.$ (HM3)

  * •

Consistency with the standard linear finite element method: If
$w\in\overline{V}^{\textup{c}}_{\ell}$, we have

$\mathcal{U}_{\ell}\gamma_{\ell}w=w\qquad\text{ and
}\qquad\boldsymbol{\mathcal{Q}}_{\ell}\gamma_{\ell}w=-\nabla w,$ (HM4)

where $\gamma_{\ell}$ is the trace operator to the skeleton
$\mathcal{F}_{\ell}$ and

$\overline{V}^{\textup{c}}_{\ell}\coloneq\bigl{\\{}v\in
H^{1}_{0}(\Omega)\;\big{|}\;v_{|T}\in\mathcal{P}_{1}(T)\;\;\forall
T\in\mathcal{T}_{\ell}\bigr{\\}}.$

  * •

Convergence of the skeletal unknown to the traces of the analytical solution.
That is, if $m_{\ell}$ is the skeletal function of the hybrid approximation to
$u\in H^{1+\alpha}(\Omega)$, we have

$\|\boldsymbol{\mathcal{Q}}_{\ell}m_{\ell}-\boldsymbol{\mathcal{Q}}_{\ell}\Pi^{\partial}_{\ell}u\|_{\Omega}\lesssim
h_{\ell}^{\alpha}\|u\|_{\alpha+1}.$ (HM5)

  * •

The usual bounds on the eigenvalues of the condensed discretization matrix
hold:

$\|\mu\|^{2}_{\ell}\lesssim a_{\ell}(\mu,\mu)\lesssim
h^{-2}_{\ell}\|\mu\|^{2}_{\ell}.$ (HM6)

  * •

If $\rho=\gamma_{\ell}w$ for $w\in\overline{V}^{\textup{c}}_{\ell}$, the
stabilization satisfies

$s_{\ell}(\rho,\mu)=0.$ (HM7)

  * •

The dual variable bounds the stabilization term

$s_{\ell}(\mu,\mu)\lesssim\|\boldsymbol{\mathcal{Q}}\mu\|^{2}_{\Omega}.$ (HM8)

######  Remark 3.2 (HDG methods).

For the LDG-H (with $\tau_{\ell}h_{\ell}\lesssim 1$), the RT-H, and the BDM-H
methods, (HM7) follows directly from (HM4), and (HM8) is an obvious
consequence of their analog to (HM1). They use almost identical assumptions to
prove multigrid convergence for the HDG, BDM-H, and RT-H methods. (HM2)–(HM6)
have been used to prove the convergence of multigrid methods for HDG, while
(HM1), (HM7), and (HM8) are novel assumptions that allow to generalize the
convergence theory to HHO.

###  3.3. Assumptions on injection operators

Our analytical findings will rely on the following conditions to hold on all
mesh levels $\ell$:

  1. (1)

Boundedness:

$\|I_{\ell}\rho\|_{\ell}\lesssim\|\rho\|_{\ell-1}\qquad\forall\rho\in
M_{\ell-1}.$ (IA1)

  2. (2)

Conformity with linear finite elements:

$I_{\ell}\gamma_{\ell-1}w=\gamma_{\ell}w\qquad\forall
w\in\overline{V}^{\textup{c}}_{\ell-1}.$ (IA2)

This way, we do not restrict ourselves to one specific injection operator. All
injection operators in [LRK22a] satisfy the desired properties. An immediate
consequence of (HM4) and (IA2) is

######  Lemma 3.3 (Quasi-orthogonality).

Let (HM4), (HM7), and (IA2) hold, then we have for all $\mu\in M_{\ell}$ and
$w\in\overline{V}^{\textup{c}}_{\ell-1}$ that

$(\boldsymbol{\mathcal{Q}}_{\ell}\mu-\boldsymbol{\mathcal{Q}}_{\ell-1}P^{a}_{\ell-1}\mu,\nabla
w)_{\Omega}=0.$

###### Proof.

For $\rho=\gamma_{\ell-1}w$, assumptions (IA2) and (HM4), we obtain using the
embedding
$\overline{V}^{\textup{c}}_{\ell-1}\subset\overline{V}^{\textup{c}}_{\ell}$

$\displaystyle\boldsymbol{\mathcal{Q}}_{\ell-1}\rho=\boldsymbol{\mathcal{Q}}_{\ell-1}\gamma_{\ell-1}w=-\nabla
w=\boldsymbol{\mathcal{Q}}_{\ell}\gamma_{\ell}w=\boldsymbol{\mathcal{Q}}_{\ell}I_{\ell}\gamma_{\ell-1}w=\boldsymbol{\mathcal{Q}}_{\ell}I_{\ell}\rho.$

The definitions of $a_{\ell}$, $a_{\ell-1}$, and (HM7) yield

$\displaystyle a_{\ell-1}(P^{a}_{\ell-1}\mu,\rho)$
$\displaystyle=(\boldsymbol{\mathcal{Q}}_{\ell-1}P^{a}_{\ell-1}\mu,\boldsymbol{\mathcal{Q}}_{\ell-1}\rho)_{\Omega},$
$\displaystyle a_{\ell}(\mu,I_{\ell}\rho)$
$\displaystyle=(\boldsymbol{\mathcal{Q}}_{\ell}\mu,\boldsymbol{\mathcal{Q}}_{\ell}I_{\ell}\rho)_{\Omega}.$

The two lines are equal by definition (3.19) of the projection operator
$P^{a}_{\ell-1}$. Thus, taking the difference gives the result. ∎

##  4\. Multigrid algorithm and abstract convergence results

If (2.1) has full elliptic regularity, which means that (2.2) holds for
$\alpha=1$, we use a standard (symmetric) V-cycle multigrid method. Otherwise,
we use a variable V-cycle multigrid method to solve the system of linear
equations arising from (3.7).

We follow the lines of [LRK22a] and avoid a discussion about good smoothing
operators. That is, we allow our smoother in smoothing step $i$

$R^{i}_{\ell}\colon M_{\ell}\to M_{\ell}$

to fit the criteria of [BP92], which allow pointwise Jacobi and Gauss-Seidel
methods. Next, we present the multigrid method as in [BPX91]. Afterwards, we
present the abstract convergence results of [DGTZ07].

###  4.1. Multigrid algorithm

We recursively define the multigrid operator of the refinement level $\ell$

$B_{\ell}\colon M_{\ell}\ni\mu\mapsto B_{\ell}\mu\in M_{\ell}$

with $n_{\ell}\in\mathbb{N}\setminus\\{0\\}$ smoothing steps on level $\ell$:
Let $B_{0}=A^{-1}_{0}$ be the exact inverse of $A_{\ell}$ on the coarsest
level. For $\ell>0$, set $x^{0}=0\in M_{\ell}$.

  1. (1)

Perform $n_{\ell}$ smoothing steps

$x^{i}=x^{i-1}+R_{\ell}^{i}(\mu-A_{\ell}x^{i-1}).$

  2. (2)

Perform recursive multigrid step

$q=B_{\ell-1}P^{\langle\cdot\rangle}_{\ell-1}(\mu-A_{\ell}x^{n_{\ell}}),$

and set $y^{0}=x^{n_{\ell}}+I_{\ell}q$.

  3. (3)

Perform $n_{\ell}$ smoothing steps

$y^{i}=y^{i-1}+R^{i+n_{\ell}}_{\ell}(\mu-A_{\ell}y^{i-1}).$

  4. (4)

Set $B_{\ell}\mu=y^{n_{\ell}}$.

If $n_{\ell}=n$ independent of the level, we obtain the standard V-cycle. The
variable V-cycle is characterized by

$\displaystyle\rho_{1}n_{\ell}\leq n_{\ell-1}\leq\rho_{2}n_{\ell},$ (4.1)

which needs to hold for all $\ell\leq L$ and uniform constants
$1<\rho_{1}\leq\rho_{2}$.

###  4.2. Main convergence results

We use the results obtained in [BPX91, DGTZ07]. These publications show
convergence under three abstract assumptions, which are assumed to hold for
all level $\ell$. In order to illustrate these assumptions, let
$\lambda^{A}_{\ell}$ be the largest eigenvalue of $A_{\ell}$, and

$\displaystyle
K_{\ell}\coloneq\bigl{(}\mathds{1}-(\mathds{1}-R_{\ell}A_{\ell})(\mathds{1}-R^{\prime}_{\ell}A_{\ell})\bigr{)}A^{-1}_{\ell},$

where $\mathds{1}$ is the identity matrix and $R^{\prime}_{\ell}$ is the
transpose of $R_{\ell}$. The aforementioned, abstract assumptions claim
existence of constants $C_{1},C_{2},C_{3}>0$ independent of the mesh level
$\ell$ and of the function $\mu\in M_{\ell}$ satisfying:

  * •

Regularity of approximation:

$|a_{\ell}(\mu-I_{\ell}P^{a}_{\ell-1}\mu,\mu)|\leq
C_{1}\left(\frac{\|A_{\ell}\mu\|^{2}_{\ell}}{\lambda^{A}_{\ell}}\right)^{\alpha}a_{\ell}(\mu,\mu)^{1-\alpha}.$
(A1)

If $\alpha=1$ in (2.2), (A1) simplifies to

$|a_{\ell}(\mu-I_{\ell}P^{a}_{\ell-1}\mu,\mu)|\leq
C_{1}\frac{\|A_{\ell}\mu\|^{2}_{\ell}}{\lambda^{A}_{\ell}}.$ (A1’)

  * •

Boundedness of the composition $I_{\ell}\circ P^{a}_{\ell-1}\colon M_{\ell}\to
M_{\ell}$ of the injection and Ritz quasi-projection operators:

$\|\mu-I_{\ell}P^{a}_{\ell-1}\mu\|_{a_{\ell}}\leq C_{2}\|\mu\|_{a_{\ell}}.$
(A2)

  * •

Smoothing hypothesis:

$\frac{\|\mu\|^{2}_{\ell}}{\lambda^{A}_{\ell}}\leq C_{3}\langle
K_{\ell}\mu,\mu\rangle_{\ell}.$ (A3)

######  Theorem 4.1.

Let (2.2) hold with $\alpha=1$, as well as (A1’), (A2), and (A3). Then, for
the standard V-cycle, for all $\ell\geq 0$, and for all $\mu\in M_{\ell}$

$|a_{\ell}(\mu-B_{\ell}A_{\ell}\mu,\mu)|\leq\delta a_{\ell}(\mu,\mu),$

where

$\delta=\frac{C_{1}C_{3}}{n-C_{1}C_{3}}\qquad\text{with}\qquad n>2C_{1}C_{3}.$

###### Proof.

This result is an immediate consequence of [DGTZ07, Thm. 3.1]. ∎

######  Theorem 4.2.

Let (2.2) hold for $\alpha\in(\tfrac{1}{2},1]$. Further, assumptions (A1),
(A2), and (A3) shall be true. If, additionally, (4.1) holds with
$\rho_{1},\rho_{2}>1$, then we have, for all $\ell\geq 0$ and all $\mu\in
M_{\ell}$, that

$\eta_{0}a_{\ell}(\mu,\mu)\leq
a_{\ell}(B_{\ell}A_{\ell}\mu,\mu)\leq\eta_{1}a_{\ell}(\mu,\mu)$

holds with

$\eta_{0}\geq\frac{n_{\ell}^{\alpha}}{M+n_{\ell}^{\alpha}},\qquad\eta_{1}\leq\frac{M+n_{\ell}^{\alpha}}{n_{\ell}^{\alpha}},$

where the constant $M>0$ is independent of $\ell$. Hence, the condition number
of $B_{\ell}A_{\ell}$ does not depend on $\ell$.

###### Proof.

This result is an immediate consequence of [BPX91, Thm. 6]. ∎

##  5\. Convergence analysis

We follow the lines of [LRK22a] and extend their results to our more general
framework. Thus, the theorems and lemmas in the following sections will be
very similar to the ones in [LRK22a]. If Lu et al. [LRK22a] have not used any
HDG-specific arguments in their proofs, we also accept the results to hold in
our framework and cite them. Otherwise, we will redo the respective proofs.

###  5.1. Energy boundedness of the injection and proof of (A2)

We provide a boundedness result for $I_{\ell}$ concerning several measures. In
particular, we investigate the injection operator’s energy boundedness and the
“Ritz projection” to show (A2). To this end, we use the averaging linear
interpolation

$I^{\textup{avg}}_{\ell}\colon
V_{\ell}\to\overline{V}^{\textup{c}}_{\ell},\qquad
I^{\textup{avg}}_{\ell}u({\boldsymbol{x}}_{\boldsymbol{a}})\coloneq\\{\\!\\{u\\}\\!\\}_{\boldsymbol{a}}\text{
for any internal mesh vertex $\boldsymbol{a}$}.$ (5.1)

Here, $\boldsymbol{x}_{\boldsymbol{a}}$ is the point that corresponds to
vertex $\boldsymbol{a}$, and $\\{\\!\\{u\\}\\!\\}_{\boldsymbol{a}}$ describes
the arithmetic mean of all values that $u$ attains in $\boldsymbol{a}$. For
$\boldsymbol{x}\in\partial\Omega$, we set
$I^{\textup{avg}}_{\ell}u(\boldsymbol{x})\coloneq 0$.

######  Lemma 5.1.

Assuming (HM1) and (HM3), we have for all $\ell$ that

$\displaystyle\|\nabla
I^{\textup{avg}}_{\ell}\mathcal{U}_{\ell}\mu\|_{\Omega}\leavevmode\nobreak\
\lesssim\leavevmode\nobreak\ $ $\displaystyle\|\mu\|_{a_{\ell}},$
$\displaystyle\forall\mu\in M_{\ell},$ (5.2)
$\displaystyle\|\mathcal{U}_{\ell}\lambda-I^{\textup{avg}}_{\ell}\mathcal{U}_{\ell}\lambda\|_{\Omega}\leavevmode\nobreak\
\lesssim\leavevmode\nobreak\ $ $\displaystyle h_{\ell}\|\mu\|_{a_{\ell}}$
$\displaystyle\forall\lambda\in M_{\ell}.$ (5.3)

###### Proof.

This can be shown analogously to [LRK22a, Lem. 5.2]. However, their (LS1)
needs to be replaced by our version of (HM1), which changes the right-hand
sides in (5.2) and (5.3). ∎

######  Lemma 5.2.

Assuming (HM1) and (HM3), we have

$\|\mu-\gamma_{\ell}I^{\textup{avg}}_{\ell}\mathcal{U}_{\ell}\mu\|_{\ell}\lesssim
h_{\ell}\|\mu\|_{a_{\ell}}\qquad\forall\mu\in M_{\ell}.$

###### Proof.

This is [LRK22a, Lem. 5.3], but with our version of (HM1). ∎

######  Lemma 5.3.

Let (IA1), (IA2), (HM1)–(HM4) hold. We have

$\displaystyle\|\boldsymbol{\mathcal{Q}}_{\ell}I_{\ell}\mu\|_{\Omega}\leavevmode\nobreak\
\lesssim\leavevmode\nobreak\ $ $\displaystyle\|\mu\|_{a_{\ell-1}},$ (5.4)
$\displaystyle\|\mathcal{U}_{\ell-1}\mu-\mathcal{U}_{\ell}I_{\ell}\mu\|_{\Omega}\leavevmode\nobreak\
\lesssim\leavevmode\nobreak\ $ $\displaystyle
h_{\ell-1}\|\mu\|_{a_{\ell-1}}\lesssim h_{\ell}\|\mu\|_{a_{\ell-1}}$ (5.5)

for all $\mu\in M_{\ell-1}$. If, additionally, (HM8) holds, we have

$\displaystyle\|I_{\ell}\mu\|_{a_{\ell}}$
$\displaystyle\lesssim\|\mu\|_{a_{\ell-1}}$ $\displaystyle\forall\mu\in
M_{\ell-1},$ (5.6) $\displaystyle\|P^{a}_{\ell-1}\mu\|_{a_{\ell-1}}$
$\displaystyle\lesssim\|\mu\|_{a_{\ell}}$ $\displaystyle\forall\mu\in
M_{\ell}.$ (5.7)

###### Proof.

Inequalities (5.4) and (5.5) can be obtained as in [LRK22a, Lem. 5.1].

Inequality (5.6) is obtained as follows:

$\|I_{\ell}\mu\|_{a_{\ell}}^{2}=a_{\ell}(I_{\ell}\mu,I_{\ell}\mu)\overset{\eqref{EQ:bilinear}}{=}\|\boldsymbol{\mathcal{Q}}_{\ell}I_{\ell}\mu\|_{\Omega}^{2}+s_{\ell}(I_{\ell}\mu,I_{\ell}\mu)\overset{\eqref{EQ:LS8}}{\lesssim}\|\boldsymbol{\mathcal{Q}}_{\ell}I_{\ell}\mu\|_{\Omega}^{2}\overset{\eqref{EQ:Q_mehses}}{\lesssim}\|\mu\|_{a_{\ell-1}}^{2},$

while (5.7) results from

$\|P^{a}_{\ell-1}\mu\|^{2}_{a_{\ell-1}}=a_{\ell-1}(P^{a}_{\ell-1}\mu,P^{a}_{\ell-1}\mu)\overset{\eqref{EQ:projection_definition}}{=}a_{\ell}(\mu,I_{\ell}P^{a}_{\ell-1}\mu)\\\
\leq\|\mu\|_{a_{\ell}}\|I_{\ell}P^{a}_{\ell-1}\mu\|_{a_{\ell}}\overset{\eqref{EQ:i_stable}}{\lesssim}\|\mu\|_{a_{\ell}}\|P^{a}_{\ell-1}\mu\|_{a_{\ell-1}},$

where we have used the Bunyakovsky–Cauchy–Schwarz inequality to pass to the
second line. ∎

######  Lemma 5.4.

Under the assumptions of Lemma 5.3, (A2) holds.

###### Proof.

Using first the linearity and symmetry of $a_{\ell}$ and then the definition
(3.19) of the Ritz quasi-projector $P^{a}_{\ell-1}$, we obtain

$\displaystyle a_{\ell}(\mu-I_{\ell}P^{a}_{\ell-1}\mu,\mu-
I_{\ell}P^{a}_{\ell-1}\mu)$
$\displaystyle\quad=a_{\ell}(\mu,\mu)-2a_{\ell}(\mu,I_{\ell}P^{a}_{\ell-1}\mu)+a_{\ell}(I_{\ell}P^{a}_{\ell-1}\mu,I_{\ell}P^{a}_{\ell-1}\mu)$
$\displaystyle\quad\leq\|\mu\|_{a_{\ell}}^{2}\underbrace{-2a_{\ell-1}(P^{a}_{\ell-1}\mu,P^{a}_{\ell-1}\mu)}_{\leq
0}+\|I_{\ell}P^{a}_{\ell-1}\mu\|_{a_{\ell}}^{2}.$

To conclude, we write for the rightmost term

$\|I_{\ell}P^{a}_{\ell-1}\mu\|_{a_{\ell}}^{2}\overset{\eqref{EQ:i_stable}}{\lesssim}\|P^{a}_{\ell-1}\mu\|_{a_{\ell-1}}^{2}\overset{\eqref{EQ:ap_stable}}{\lesssim}\|\mu\|^{2}_{a_{\ell}}.\qed$

###  5.2. Proofs of (A1) and (A1’)

A classical approach for multigrid proofs consists is considering
$A_{\ell}\mu$ (cf. (3.17)) as a function in $L^{2}(\Omega)$ and using it as
right-hand side in an auxiliary problem. However, this is not possible for
discontinuous skeletal methods, since $A_{\ell}\mu\in M_{\ell}$ is defined
only on the mesh skeleton and not on the whole domain $\Omega$. We address
this issue and construct a suitable right-hand side $f_{\mu}$ from $\mu$ with
the help of the surjective lifting operator $S_{\ell}\colon M_{\ell}\to
V_{\textup{disc}}$ (with $V_{\textup{disc}}\subset C(\Omega)\cap
L^{2}(\Omega)$ being a suitable discrete space), which is inspired by [Tan09].
It is also used in [GR86, Lem. A.3] for the two dimensional case and in
[Mon03, Def. 5.46] for three dimensional settings:

$\displaystyle(S_{\ell}\mu,v)_{T}$
$\displaystyle=(\mathcal{U}_{\ell}\mu,v)_{T}$ $\displaystyle\forall
v\in\mathcal{P}_{p}(T),\,\forall T\in\mathcal{T}_{\ell},$ (5.8a)
$\displaystyle\langle S_{\ell}\mu,\eta\rangle_{F}$
$\displaystyle=\langle\mu,\eta\rangle_{F}$
$\displaystyle\forall\eta\in\mathcal{P}_{p+1}(F),\,\forall
F\in\mathcal{F}_{\ell},$ (5.8b) $\displaystyle S_{\ell}\mu(\boldsymbol{a})$
$\displaystyle=\\{\\!\\{\mu\\}\\!\\}_{\boldsymbol{a}}$
$\displaystyle\forall\boldsymbol{a}\text{ is a vertex in
}\mathcal{T}_{\ell},\,\boldsymbol{a}\not\in\partial\Omega,$ (5.8c)
$\displaystyle S_{\ell}\mu(\boldsymbol{a})$ $\displaystyle=0$
$\displaystyle\forall\boldsymbol{a}\text{ is a vertex in
}\mathcal{T}_{\ell},\,\boldsymbol{a}\in\partial\Omega$ (5.8d) in two
dimensions. In three spatial dimensions, we add the constraints
$\displaystyle\langle
S_{\ell}\mu,\eta\rangle_{\Gamma}=\langle\\{\\!\\{\mu\\}\\!\\}_{\Gamma},\eta\rangle_{\Gamma}\qquad\forall\eta\in\mathcal{P}_{p+2}(\Gamma),$
(5.8e)

Here, $\\{\\!\\{m\\}\\!\\}_{\Gamma}$ is the average taken over all cells
adjacent to $\Gamma$. This construction extends to dimensions higher than
three. For a precise definition of $V_{\textup{disc}}$ and a proof of the fact
that $S_{\ell}$ is well-defined, please refer to [GR86, Lem. A.3] and [Mon03,
Def. 5.46]. Moreover, (HM4) implies that

$S_{\ell}\gamma_{\ell}w=w\qquad\forall w\in\overline{V}^{\textup{c}}_{\ell}.$
(5.9)

Eventually, we define the auxiliary right-hand side $f_{m}\in
V_{\textup{disc}}$ as the unique solution of

$(f_{\mu},S_{\ell}\eta)=\langle
A_{\ell}\mu,\eta\rangle_{\ell}=a_{\ell}(\mu,\eta)\qquad\forall\eta\in
M_{\ell}.$ (5.10)

######  Lemma 5.5.

Let (HM1)–(HM4) hold. Then we have, for all $\mu\in M_{\ell}$,

$\|f_{\mu}\|_{\Omega}\lesssim\|A_{\ell}\mu\|_{\ell}.$ (5.11)

###### Proof.

This is part of [LRK22a, Lem. 5.6]. The change in the meaning of (HM1) does
not influence this result. ∎

######  Lemma 5.6 (Reconstruction approximation).

Let (IA1), (IA2), and (HM1)–(HM5) hold. Then, if the model problem has
regularity (2.2), for all $\mu\in M_{\ell}$ there exists an auxiliary function
$\overline{u}\in\overline{V}^{\textup{c}}_{\ell-1}$ such that

$\|\boldsymbol{\mathcal{Q}}_{\ell}\mu+\nabla\overline{u}\|_{\Omega}+\|\boldsymbol{\mathcal{Q}}_{\ell-1}P^{a}_{\ell-1}\mu+\nabla\overline{u}\|_{\Omega}\lesssim
h_{\ell}^{\alpha}\|f_{\mu}\|_{\alpha-1}.$ (5.12)

###### Proof.

This is [LRK22a, Lem. 5.8]. ∎

######  Lemma 5.7.

Let (HM1)–(HM4), (HM6), and (HM7) hold. Then,

$\displaystyle\|f_{\mu}\|_{-1}\lesssim\|\mu\|_{a_{\ell}}.$

###### Proof.

By the definition of negative norms and properties of $\sup$,

$\|f_{\mu}\|_{-1}\leq\sup_{\psi\in H^{1}_{0}(\Omega)}\frac{(f_{\mu},\psi-
S_{\ell}\gamma_{\ell}\overline{\Pi}^{\textup{c}}_{\ell}\psi)}{|\psi|_{1}}+\sup_{\psi\in
H^{1}_{0}(\Omega)}\frac{(f_{\mu},S_{\ell}\gamma_{\ell}\overline{\Pi}^{\textup{c}}_{\ell}\psi)}{|\psi|_{1}}.$
(5.13)

Here, $\overline{\Pi}^{\textup{c}}_{\ell}$ is an (quasi-)interpolator onto
$\overline{V}^{\textup{c}}_{\ell}$ that satisfies

$\displaystyle|\overline{\Pi}^{\textup{c}}_{\ell}v|_{1}$
$\displaystyle\lesssim|v|_{1}$ $\displaystyle\forall v\in H^{1}(\Omega),$
(5.14) $\displaystyle\|v-\overline{\Pi}^{\textup{c}}_{\ell}v\|_{\Omega}$
$\displaystyle\lesssim h_{\ell}^{1-k+\alpha}|v|_{\alpha+1},$
$\displaystyle\forall v\in H^{\alpha+1}(\Omega),\;k=0,1.$ (5.15)

An example is given in [SZ90]. We continue with bounding

$\|\psi-
S_{\ell}\gamma_{\ell}\overline{\Pi}^{\textup{c}}_{\ell}\psi\|_{\Omega}\overset{\eqref{EQ:lift_iden}}{=}\|\psi-\overline{\Pi}^{\textup{c}}_{\ell}\psi\|_{\Omega}\overset{\eqref{EQ:H1_stab}}{\lesssim}h_{\ell}|\psi|_{1},$
(5.16)

and observe that

$(f_{\mu},S_{\ell}\gamma_{\ell}\overline{\Pi}^{\textup{c}}_{\ell}\psi)_{\Omega}\overset{\eqref{EQ:scp_def}}{=}a_{\ell}(\mu,\gamma_{\ell}\overline{\Pi}^{\textup{c}}_{\ell}\psi),$

which immediately yields, using a Cauchy–Schwarz inequality,

$|(f_{\mu},S_{\ell}\gamma_{\ell}\overline{\Pi}^{\textup{c}}_{\ell}\psi)_{\Omega}|\leq\|\mu\|_{a_{\ell}}\|\gamma_{\ell}\overline{\Pi}^{\textup{c}}_{\ell}\psi\|_{a_{\ell}}=\|\mu\|_{a_{\ell}}|\overline{\Pi}^{\textup{c}}_{\ell}\psi|_{1}\lesssim\|\mu\|_{a_{\ell}}|\psi|_{1},$
(5.17)

where $|\cdot|_{1}$ denotes the $H^{1}(\Omega)$-seminorm. Applying a
Cauchy–Schwarz inequality to the numerator of the first supremum in (5.13)
followed by (5.16) and using (5.17) to estimate the numerator of the second
supremum, we get, after simplification,

$\|f_{\mu}\|_{-1}\leq
h_{\ell}\|f_{\mu}\|_{\Omega}+\|\mu\|_{a_{\ell}}\overset{\eqref{EQ:bound_flambda}}{\lesssim}h_{\ell}\|A_{\ell}\mu\|_{\ell}+\|\mu\|_{a_{\ell}}\lesssim\|\mu\|_{a_{\ell}},$

where the last inequality is the rightmost inequality of (HM6). ∎

With these preliminary results at hand, we can formulate a generalization of
[LRK21, Theo. 4.1] and [LRK22a, Theo. 5.10]:

######  Theorem 5.8.

If (2.2) holds with $\alpha\in(1/2,1]$, and if the assumptions (IA1), (IA2),
and (HM1)–(HM8) hold, then (A1) is satisfied.

######  Remark 5.9.

Theorem 5.8 implies (A1’) if $\alpha=1$.

###### Proof.

We use the definitions (3.8) of $a_{\ell}$ and (3.19) of $P^{a}_{\ell}$ to
obtain

$\displaystyle a_{\ell}(\mu-I_{\ell}P^{a}_{\ell-1}\mu,\mu)$
$\displaystyle=a_{\ell}(\mu,\mu)-a_{\ell-1}(P^{a}_{\ell-1}\mu,P^{a}_{\ell-1}\mu)$
$\displaystyle=\underbrace{(\boldsymbol{\mathcal{Q}}_{\ell}\mu,\boldsymbol{\mathcal{Q}}_{\ell}\mu)_{\Omega}-(\boldsymbol{\mathcal{Q}}_{\ell-1}P^{a}_{\ell-1}\mu,\boldsymbol{\mathcal{Q}}_{\ell-1}P^{a}_{\ell-1}\mu)_{\Omega}}_{\mathfrak{T}_{1}}$
$\displaystyle\quad+\underbrace{s_{\ell}(\mu,\mu)-s_{\ell-1}(P^{a}_{\ell-1}\mu,P^{a}_{\ell-1}\mu),}_{\mathfrak{T}_{2}}$

and analyze the respective contributions $\mathfrak{T}_{1}$ and
$\mathfrak{T}_{2}$ separately.

First, binomial factorization yields

$\displaystyle\mathfrak{T}_{1}=(\boldsymbol{\mathcal{Q}}_{\ell}\mu+\boldsymbol{\mathcal{Q}}_{\ell-1}P^{a}_{\ell-1}\mu,\boldsymbol{\mathcal{Q}}_{\ell}\mu-\boldsymbol{\mathcal{Q}}_{\ell-1}P^{a}_{\ell-1}\mu)_{\Omega}.$

Let now $w\in\overline{V}^{\textup{c}}_{\ell-1}$. Invoking the orthogonality
property stated in Lemma 3.3 to insert $2\nabla w$ into the first argument of
the $L^{2}$-product, and adding and subtracting $\nabla w$ to the second
argument, we obtain

$\displaystyle\mathfrak{T}_{1}$
$\displaystyle=(\boldsymbol{\mathcal{Q}}_{\ell}\mu+2\nabla
w+\boldsymbol{\mathcal{Q}}_{\ell-1}P^{a}_{\ell-1}\mu,\boldsymbol{\mathcal{Q}}_{\ell}\mu+\nabla
w-\nabla w-\boldsymbol{\mathcal{Q}}_{\ell-1}P^{a}_{\ell-1}\mu)_{\Omega}$
(5.18) $\displaystyle\leq\left(\|\boldsymbol{\mathcal{Q}}_{\ell}+\nabla
w\|_{\Omega}+\|\boldsymbol{\mathcal{Q}}_{\ell-1}P^{a}_{\ell-1}\mu+\nabla
w\|_{\Omega}\right)^{2}\overset{\eqref{EQ:reconst_approx}}{\lesssim}h_{\ell}^{2\alpha}\|f_{\mu}\|_{\alpha-1}^{2},$

where we have used Bunyakovsky–Cauchy–Schwarz and triangle inequalities to
pass to the second line.

Second, we obtain a similar estimate for $\mathfrak{T}_{2}$. To this end, we
detail the treatment of the first summand of $\mathfrak{T}_{2}$, the second
summand can be treated analogously. For any
$w\in\overline{V}^{\textup{c}}_{\ell}$, it holds

$s_{\ell}(\mu,\mu)\overset{\eqref{EQ:LS7}}{=}s_{\ell}(\mu-\gamma_{\ell}w,\mu-\gamma_{\ell}w)\overset{\eqref{EQ:LS8}}{\lesssim}\|\boldsymbol{\mathcal{Q}}_{\ell}(\mu-\gamma_{\ell}w)\|^{2}_{\Omega}\overset{\eqref{EQ:reconst_approx}}{\lesssim}h_{\ell}^{2\alpha}\|f_{\mu}\|_{\alpha-1}^{2}.$

Proceeding similarly for the second summand in $\mathfrak{T}_{2}$, we infer

$\mathfrak{T}_{2}\lesssim h_{\ell}^{2\alpha}\|f_{\mu}\|_{\alpha-1}^{2}.$
(5.19)

Using (5.18) and (5.19) to estimate the right-hand side of (5.2) and invoking
Sobolev interpolation provides us with

$|a_{\ell}(\mu-I_{\ell}P^{a}_{\ell-1}\mu,\mu)|\lesssim
h_{\ell}^{2\alpha}\|f_{\mu}\|_{\alpha-1}^{2}\lesssim
h_{\ell}^{2\alpha}\|f_{\mu}\|^{2(1-\alpha)}_{-1}\|f_{\mu}\|_{\Omega}^{2\alpha},$

and Lemmas 5.5 and 5.7 yield

$|a_{\ell}(\mu-I_{\ell}P^{a}_{\ell-1}\mu,\mu)|\lesssim
h_{\ell}^{2\alpha}\|A_{\ell}\mu\|^{2\alpha}_{\ell}\|\mu\|^{2(1-\alpha)}_{a_{\ell}}.$

This is the result, since (HM6) implies that
$\underline{\lambda}^{A}_{\ell}\lesssim h_{\ell}^{-2}$. ∎

##  6\. Verification of assumptions for HHO

###  6.1. Preliminaries

For all $T\in\mathcal{T}_{\ell}$ (resp. $F\in\mathcal{F}_{\ell}$), we denote
by $\|\cdot\|_{T}$ (resp. $\|\cdot\|_{F}$) the standard $L^{2}$-norm over $T$
(resp. $F$). For all $T\in\mathcal{T}_{\ell}$, we introduce the set
$\mathcal{F}_{T}$ collecting the faces of $T$.

####  6.1.1. General properties of used norms and function spaces

Let $T\in\mathcal{T}_{\ell}$, $F\in\mathcal{F}_{T}$.

  * •

The _discrete trace inequality_ (see, e.g., [DPD20, Lem. 1.32]), stipulates
that

$\|v\|_{F}\lesssim h_{T}^{-1/2}\|v\|_{T}\qquad\forall v\in\mathcal{P}(T).$
(6.1)

  * •

The _discrete inverse inequality_ (see, e.g., [DPD20, Lem. 1.28]) stipulates
that

$\|\nabla v\|_{T}\lesssim h_{T}^{-1}\|v\|_{T}\qquad\forall
v\in\mathcal{P}(T).$ (6.2)

  * •

The _Poincaré-Friedrichs inequality_ (see, e.g., [EG21, Lem. 3.30]) stipulates
that

$\|\rho\|_{T}\lesssim h_{T}\|\nabla\rho\|_{T}+h_{T}^{1/2}\|\rho\|_{\partial
T}\qquad\forall\rho\in H^{1}(T).$ (6.3)

  * •

Finally, we have the useful relation

$\|\mu\|_{\ell}^{2}\simeq\sum_{F\in\mathcal{F}_{\ell}}h_{F}\|\mu\|_{F}^{2}\qquad\forall\mu\in
M_{\ell}$ (6.4)

with $\mathcal{F}_{\ell}$ containing all the faces of $\mathcal{T}_{\ell}$.
Thus, the square root of the expression in the right-hand side induces a
global norm on $M_{\ell}$.

####  6.1.2. Properties of HHO

In this part, we enlist some properties of HHO that will be used to show
(HM1)–(HM8).

Recalling the hybrid bilinear form $\underline{a}_{\ell}$ introduced in
(3.12), we consider the hybrid, discrete problem: find $(u_{\ell},m_{\ell})\in
V_{\ell}\times M_{\ell}$ such that

$\underline{a}_{\ell}((u_{\ell},m_{\ell}),(v,\mu))=\int_{\Omega}fv\qquad\forall(v,\mu)\in
V_{\ell}\times M_{\ell}.$ (6.5)

The pair
$\underline{u}_{\ell}\coloneq(\mathcal{U}_{\ell}m_{\ell}+\mathcal{V}_{\ell}f,m_{\ell})\in
V_{\ell}\times M_{\ell}$ is solution of (6.5) if and only if $m_{\ell}\in
M_{\ell}$ is solution of the condensed problem (3.7) (see [CDPE16, Prop. 4]).
In particular, $\mathcal{U}_{\ell}m_{\ell}+\mathcal{V}_{\ell}f$ approximates
the exact solution $u$ of the continuous problem (2.1).

As in [DPD20, 2.1.2 (2.7), 2.2.2 (2.35)], we define

$\displaystyle|(v,\mu)|_{\underline{1,\ell}}^{2}\coloneq\sum_{T\in\mathcal{T}_{\ell}}|(v,\mu)|_{\underline{1,T}}^{2}=\sum_{T\in\mathcal{T}_{\ell}}\sum_{F\in\mathcal{F}_{T}}h_{F}^{-1}\|v-\mu\|_{F}^{2},$
(6.6a)
$\displaystyle\|(v,\mu)\|_{\underline{1,\ell}}^{2}\coloneq\sum_{T\in\mathcal{T}_{\ell}}\|(v,\mu)\|_{\underline{1,T}}^{2}=\sum_{T\in\mathcal{T}_{\ell}}\left(\|\nabla
v\|_{T}^{2}+|(v,\mu)|_{\underline{1,T}}^{2}\right),$ (6.6b)

as the $H^{1}$-like seminorm on the hybrid space $V_{\ell}\times M_{\ell}$.
Moreover, we assume that the properties of [DPD20, Assumption 2.4] for the
stabilization bilinear form $\underline{s}_{T}$ hold true. Then, we have:

  * •

(Boundedness and stability of $\underline{a}_{T}$ [CEP21, Lem. 2.6])

$\|(v,\mu)\|_{\underline{1,T}}^{2}\lesssim\underline{a}_{T}((v,\mu),(v,\mu))\lesssim\|(v,\mu)\|_{\underline{1,T}}^{2}.$
(6.7)

  * •

(Energy error estimate [CEP21, Lem. 2.8, Lem. 2.9])

$\|\boldsymbol{q}\underline{u}_{\ell}-\boldsymbol{q}(\Pi^{\textup{d}}_{\ell}u,\Pi^{\partial}_{\ell}u)\|_{\underline{a}_{\ell}}\lesssim
h_{\ell}^{\alpha}\|u\|_{\alpha+1},$ (6.8)

where $\alpha>\frac{1}{2}$ and $\|\cdot\|_{\underline{a}_{\ell}}$ is the norm
induced by the bilinear form $\underline{a}_{\ell}$.

  * •

(Cell unknown $L^{2}$-error estimate [CEP21, Lem. 2.11])
If $p\geq 1$, $\alpha>\frac{1}{2}$, then

$\|\mathcal{U}_{\ell}m_{\ell}+\mathcal{V}_{\ell}f-\Pi^{\textup{d}}_{\ell}u\|_{\Omega}\lesssim
h_{\ell}^{\alpha+\delta}\|u\|_{\alpha+1},$ (6.9)

where $\delta\coloneq\min\\{\alpha,1\\}$.

  * •

(Face unknown $L^{2}$-error estimate [CEP21, Lem. 2.9])
If $p\geq 1$, then

$\|m_{\ell}-\Pi^{\partial}_{\ell}u\|_{\ell}\lesssim
h_{\ell}^{\alpha+1}\|u\|_{\alpha+1}.$ (6.10)

Now, let us introduce the seminorms $\|\cdot\|_{1,\ell}$ and
$|\cdot|_{1,\ell}$ on the skeletal space $M_{\ell}$, which are based on the
seminorms $\|\cdot\|_{\underline{1,\ell}}$ and $|\cdot|_{\underline{1,\ell}}$
defined in (6.6):

$\|\mu\|_{1,\ell}\coloneq\|(\mathcal{U}_{\ell}\mu,\mu)\|_{\underline{1,\ell}}\qquad\text{
and
}\qquad|\mu|_{1,\ell}\coloneq|(\mathcal{U}_{\ell}\mu,\mu)|_{\underline{1,\ell}}.$
(6.11)

Similarly, the _condensed_ bilinear form $a_{\ell}$ introduced in (3.16) is
built from the _hybrid_ bilinear form $\underline{a}_{\ell}$ (3.12), in which
the generic cell unknown variable is also recovered from the skeletal variable
through $\mathcal{U}_{\ell}$, i.e.

$a_{\ell}(m,\mu)=\underline{a}_{\ell}((\mathcal{U}_{\ell}m,m),(\mathcal{U}_{\ell}\mu,\mu)).$

In the same fashion, we also have

$s_{\ell}(m,\mu)=\underline{s}_{\ell}((\mathcal{U}_{\ell}m,m),(\mathcal{U}_{\ell}\mu,\mu)).$

Thus, useful properties of $a_{\ell}$ and $s_{\ell}$ derive in a natural way
from those of $\underline{a}_{\ell}$ and $\underline{s}_{\ell}$. In
particular, we shall use in this work the following results:

  * •

(Boundedness and stability of $a_{\ell}$)

$\|\mu\|_{1,\ell}^{2}\lesssim a_{\ell}(\mu,\mu)\lesssim\|\mu\|_{1,\ell}^{2}.$
(6.12)

  * •

Since $\underline{s}_{\ell}$ is symmetric positive semi-definite [DPD20,
Assumption 2.4], so is $s_{\ell}$.

####  6.1.3. The norm $\|\cdot\|_{\ell}$ is weaker than $\|\cdot\|_{1,\ell}$

In this section, we prove that

$\|\mu\|_{\ell}\lesssim\|\mu\|_{1,\ell}$ (6.13)

Starting from (6.4), the triangle inequality is used via the insertion of
$0=\mathcal{U}\mu-\mathcal{U}\mu$ into the norms in the sum, yielding

$\displaystyle\|\mu\|^{2}_{\ell}$
$\displaystyle\lesssim\sum_{F\in\mathcal{F}}h_{F}\left(\|\mathcal{U}_{\ell}\mu-\mu\|^{2}_{F}+\|\mathcal{U}_{\ell}\mu\|^{2}_{F}\right)$
(6.14) $\displaystyle\overset{\eqref{EQ:discr_trace_inequality}}{}$
$\displaystyle\lesssim\sum_{F\in\mathcal{F}}h_{F}\|\mathcal{U}_{\ell}\mu-\mu\|^{2}_{F}+\sum_{T\in\mathcal{T}_{\ell}}\|\mathcal{U}_{\ell}\mu\|^{2}_{T}$
$\displaystyle\lesssim\sum_{F\in\mathcal{F}}h_{F}\|\mathcal{U}_{\ell}\mu-\mu\|^{2}_{F}+\|(\mathcal{U}_{\ell}\mu,\mu)\|_{\underline{1,\ell}}^{2},$

where the last passage follows from the discrete Poincaré inequality [DPD20,
Lem. 2.15] applied to $(\mathcal{U}_{\ell}\mu,\mu)$ (notice that $\mu$
vanishes on $\partial\Omega$). Hence, recalling (6.11),

$\|\mu\|_{\ell}^{2}\lesssim\sum_{F\in\mathcal{F}}h_{F}\|\mathcal{U}_{\ell}\mu-\mu\|_{F}^{2}+\|\nabla_{\ell}\mathcal{U}_{\ell}\mu\|_{\Omega}^{2}+|\mu|_{1,\ell}$
(6.15)

The fact that $h_{F}\lesssim h_{F}^{-1}$ is a consequence of the fact that we
assumed that $\Omega$ has diameter 1 without loss of generality so that
$h_{F}\lesssim 1$.

$\|\mu\|_{\ell}^{2}\lesssim\|\nabla_{\ell}\mathcal{U}_{\ell}\mu\|_{\Omega}^{2}+|\mu|_{1,\ell}=\|\mu\|_{1,\ell}^{2},$

and the result follows by taking the square root.

###  6.2. Verification of (HM1)

We write

$\displaystyle\|\mathcal{U}_{\ell}\mu-\mu\|_{\ell}^{2}$
$\displaystyle\lesssim\sum_{T\in\mathcal{T}_{h}}\sum_{F\in\mathcal{F}_{T}}h_{F}\|\mathcal{U}_{\ell}\mu-\mu\|_{L^{2}(F)}^{2}$
$\displaystyle\leq
h_{\ell}^{2}\sum_{T\in\mathcal{T}_{h}}\sum_{F\in\mathcal{F}_{T}}h_{F}^{-1}\|\mathcal{U}_{\ell}\mu-\mu\|_{L^{2}(F)}^{2}$
$\displaystyle\overset{\eqref{eq:norm.u1l}}{}$ $\displaystyle\leq
h_{\ell}^{2}\|(\mathcal{U}_{\ell}\mu,\mu)\|_{\underline{1,h}}^{2}\overset{\eqref{EQ:hho_hybrid_local_stability}}{\lesssim}h_{\ell}^{2}\|(\mathcal{U}_{\ell}\mu,\mu)\|_{\underline{a}_{\ell}}^{2}\overset{\eqref{eq:norm.al}}{=}h_{\ell}^{2}\|\mu\|_{a_{\ell}}^{2},$

where the first inequality follows from mesh regularity. Taking the square
root of the above inequality proves (HM1).

###  6.3. Verification of (HM2)

Let us start with the second inequality in (HM2) on an arbitrary element
$T\in\mathcal{T}_{\ell}$. Recalling the definition (3.13) of
$\mathcal{U}_{\ell}$ and choosing $v_{T}=u_{T}^{1}=\mathcal{U}_{T}\mu$ and
$m_{\partial T}=\mu$, we have, denoting by $\|\cdot\|_{\underline{a}_{T}}$ the
seminorm induced by $\underline{a}_{T}$,

$\|(\mathcal{U}_{T}\mu,0)\|_{\underline{a}_{T}}^{2}=\underline{a}_{T}((\mathcal{U}_{T}\mu,0),(\mathcal{U}_{T}\mu,0))\\\
=-\underline{a}_{T}((0,\mu),(\mathcal{U}_{T}\mu,0))\leq\|(0,\mu)\|_{\underline{a}_{T}}\,\|(\mathcal{U}_{T}\mu,0)\|_{\underline{a}_{T}},$

where the conclusion follows from the Cauchy–Schwarz inequality. Hence,

$\|(\mathcal{U}_{T}\mu,0)\|_{\underline{a}_{T}}\leq\|(0,\mu)\|_{\underline{a}_{T}}$
(6.16)

We write

$\|(\mathcal{U}_{T}\mu,0)\|_{\underline{1,T}}^{2}\overset{\eqref{EQ:hho_hybrid_local_stability}}{\lesssim}\|(\mathcal{U}_{T}\mu,0)\|_{\underline{a}_{T}}^{2}\overset{\eqref{eq:norm.Umu.le.mu}}{\leq}\|(0,\mu)\|_{\underline{a}_{T}}^{2}\overset{\eqref{EQ:hho_hybrid_local_stability}}{\lesssim}\|(0,\mu)\|_{\underline{1,T}}^{2}\lesssim
h_{T}^{-1}\|\mu\|_{\partial T}^{2},$ (6.17)

where the last inequality follows recalling the definition (6.6b) of
$\|\cdot\|_{\underline{1,T}}$ and noticing that $h_{F}^{-1}\lesssim
h_{T}^{-1}$ for all $F\in\mathcal{F}_{T}$ by mesh regularity. Next, using the
Poincaré-Friedrichs inequality (6.3) and definition (6.6b) of the
$\|(\cdot,\cdot)\|_{\underline{1,T}}$-norm along with the fact that
$h_{T}^{-1}\lesssim h_{F}^{-1}$ for all $F\in\mathcal{F}_{T}$ by mesh
regularity, we have

$\|\mathcal{U}_{T}\mu\|_{T}^{2}\lesssim
h_{T}^{2}\|\nabla(\mathcal{U}_{T}\mu)\|_{T}^{2}+h_{T}\|\mathcal{U}_{T}\mu\|_{\partial
T}^{2}\lesssim
h_{T}^{2}\|(\mathcal{U}_{T}\mu,0)\|_{\underline{1,T}}^{2}\overset{\eqref{eq:norm.U.1T.le.norm.mu.1T}}{\lesssim}h_{T}\|\mu\|_{\partial
T}^{2}.$ (6.18)

The second inequality in (HM2) is derived by using (6.18) and summing over all
elements.

Second, consider the first inequality on an arbitrary element
$T\in\mathcal{T}_{\ell}$. Recalling the definition (3.10) of
$\boldsymbol{q}_{T}(\cdot,\cdot)$ with $(u_{T},m_{\partial
T})=(\mathcal{U}_{T}\mu,\mu)$ and choosing
$\boldsymbol{p}_{T}=\boldsymbol{q}_{T}(\mathcal{U}_{T}\mu,\mu)\overset{\eqref{eq:localQ}}{=}\boldsymbol{\mathcal{Q}}_{T}\mu$,
$u_{T}=\mathcal{U}_{T}\mu$, and $m_{\partial T}=\mu$, we have

$\int_{T}\boldsymbol{\mathcal{Q}}_{T}\mu\cdot\boldsymbol{\mathcal{Q}}_{T}\mu=\int_{T}\mathcal{U}_{T}\mu\,(\nabla\\!\cdot\\!\boldsymbol{\mathcal{Q}}_{T}\mu)-\int_{\partial
T}\mu\,(\boldsymbol{\mathcal{Q}}_{T}\mu\cdot\boldsymbol{\nu}).$

Estimating the right-hand side of the above expression with Cauchy–Schwarz
inequalities followed by the discrete inverse inequality (6.2) for the first
term, the discrete trace inequality (6.1) for the second term, we have

$\|\boldsymbol{\mathcal{Q}}_{T}\mu\|_{T}^{2}\lesssim
h^{-2}_{T}\|\mathcal{U}_{T}\mu\|_{T}^{2}+h_{T}^{-1}\|\mu\|_{\partial
T}^{2}\overset{\eqref{Proof:LS2 second}}{\leq}h_{T}^{-1}\|\mu\|_{\partial
T}^{2}.$ (6.19)

Finally, the first inequality in (HM2) is derived by using (6.19), summing
over all elements, and using the mesh quasi-uniformity assumption to write
$\|\boldsymbol{\mathcal{Q}}_{\ell}\mu\|_{\Omega}\lesssim
h^{-1}_{\ell}\|\mu\|_{\ell}.$

###  6.4. Verification of (HM3)

Plugging the definition (3.15) of $\boldsymbol{\mathcal{Q}}_{T}$ into (3.10)
we have,

$\int_{T}\boldsymbol{\boldsymbol{\mathcal{Q}}}_{T}\mu\cdot\boldsymbol{p}_{T}-\int_{T}\mathcal{U}_{T}\mu\,(\nabla\\!\cdot\\!\boldsymbol{p}_{T})=-\int_{\partial
T}\mu\,(\boldsymbol{p}_{T}\cdot\boldsymbol{\nu})\quad\forall\boldsymbol{p}_{T}\in\boldsymbol{W}_{T}.$

Then, we integrate by parts the second term of the left-hand side and
rearrange to infer

$\int_{T}(\boldsymbol{\boldsymbol{\mathcal{Q}}}_{T}\mu+\nabla\mathcal{U}_{T}\mu)\cdot\boldsymbol{p}_{T}=\int_{\partial
T}(\mathcal{U}_{T}\mu-\mu)\boldsymbol{p}_{T}\cdot\boldsymbol{\nu}\quad\forall\boldsymbol{p}_{T}\in\boldsymbol{W}_{T}.$

Now, after noticing that
$\nabla\mathcal{U}_{T}\mu\in\nabla\mathcal{P}_{p}(T)\subset\nabla\mathcal{P}_{p+1}(T)=\boldsymbol{W}_{T}$,
we can specify this relation for
$\boldsymbol{p}_{T}=\boldsymbol{\mathcal{Q}}_{T}\mu+\nabla\mathcal{U}_{T}\mu$,
which gives

$\|\boldsymbol{\mathcal{Q}}_{T}\mu+\nabla\mathcal{U}_{T}\mu\|_{T}^{2}=\int_{\partial
T}(\mathcal{U}_{T}\mu-\mu)\,(\boldsymbol{\mathcal{Q}}_{T}\mu+\nabla\mathcal{U}_{T}\mu)\cdot\boldsymbol{\nu}.$

Using the Cauchy–Schwarz inequality on the right-hand side, we have

$\displaystyle\|\boldsymbol{\mathcal{Q}}_{T}\mu+\nabla\mathcal{U}_{T}\mu\|_{T}\leq\|\mathcal{U}_{T}\mu-\mu\|_{\partial
T}\|\boldsymbol{\mathcal{Q}}_{T}\mu+\nabla\mathcal{U}_{T}\mu\|_{\partial T}.$

At this point, we refer to the discrete trace inequality (6.1), which we use
component by component to bound the last term and obtain

$\|\boldsymbol{\mathcal{Q}}_{T}\mu+\nabla\mathcal{U}_{T}\mu\|_{T}\lesssim
h_{T}^{-1/2}\|\mathcal{U}_{T}\mu-\mu\|_{\partial
T}\|\boldsymbol{\mathcal{Q}}_{T}\mu+\nabla\mathcal{U}_{T}\mu\|_{T}.$

Simplifying and squaring, we get

$\|\boldsymbol{\mathcal{Q}}_{T}\mu+\nabla\mathcal{U}_{T}\mu\|_{T}^{2}\lesssim
h_{T}^{-1}\|\mathcal{U}_{T}\mu-\mu\|_{\partial T}^{2}.$

Summing over all elements and using mesh quasi-uniformity yields

$\displaystyle\|\boldsymbol{\mathcal{Q}}_{T}\mu+\nabla\mathcal{U}_{T}\mu\|_{\Omega}^{2}$
$\displaystyle\lesssim\sum_{T\in\mathcal{T}_{\ell}}h_{T}^{-1}\|\mathcal{U}_{T}\mu-\mu\|_{\partial
T}^{2}\lesssim\sum_{T\in\mathcal{T}_{\ell}}h_{T}^{-2}\sum_{F\in\mathcal{F}_{T}}h_{F}\|\mathcal{U}_{T}\mu-\mu\|_{\partial
T}^{2}$ $\displaystyle\lesssim
h_{\ell}^{-2}\sum_{T\in\mathcal{T}_{\ell}}\sum_{F\in\mathcal{F}_{T}}h_{F}\|\mathcal{U}_{T}\mu-\mu\|_{\partial
T}^{2}\lesssim h_{\ell}^{-1}\|\mathcal{U}_{T}\mu-\mu\|_{\ell}^{2},$

where we have used the fact that $1\lesssim\frac{h_{F}}{h_{T}}$ for all
$T\in\mathcal{T}_{\ell}$ and all $F\in\mathcal{F}_{T}$ by mesh regularity in
the second inequality, the mesh quasi-uniformity assumption to write
$h_{T}^{-2}\lesssim h_{\ell}^{-2}$ in the third inequality, and (6.4) to
conclude

###  6.5. Verification of (HM4)

We show (HM4) for a generic element $T\in\mathcal{T}_{\ell}$. If the
identities hold there, they will hold on all elements. We need that
$\gamma_{\ell}\overline{V}^{\textup{c}}_{\ell}$ is a subspace of $M_{\ell}$
and $V_{\ell}$. That is, we need that $p\geq 1$. We must show that HHO
reproduces the $w\in\overline{V}^{\textup{c}}_{\ell}$.

We take $w_{\partial T}\coloneq\gamma_{\ell}w$.

$\int_{T}\boldsymbol{q}_{T}(w,w_{\partial
T})\cdot\boldsymbol{p}_{T}=\int_{T}w\nabla\\!\cdot\\!\boldsymbol{p}_{T}-\int_{\partial
T}w_{\partial
T}\boldsymbol{p}_{T}\cdot\boldsymbol{\nu}\quad\forall\boldsymbol{p}_{T}\in\boldsymbol{W}_{T}.$

Integrating by parts the first term on the right-hand side, it holds that

$\displaystyle\int_{T}w\nabla\\!\cdot\\!\boldsymbol{p}_{T}=-\int_{T}\nabla
w\cdot p_{T}+\int_{\partial T}w_{\partial
T}\boldsymbol{p}_{T}\cdot\boldsymbol{\nu}$

Plugging that yields

$\int_{T}\boldsymbol{q}_{T}(w,w_{\partial
T})\cdot\boldsymbol{p}_{T}=-\int_{T}\nabla w\cdot
p_{T}\quad\forall\boldsymbol{p}_{T}\in\boldsymbol{W}_{T},$

which shows that

$\boldsymbol{q}_{T}(w,w_{\partial T})=-\nabla w.$ (6.20)

Given the definition (3.13) of $\mathcal{U}_{T}$, we need to show that

$\underline{a}_{T}((w,w_{\partial T}),(v_{T},0))=0\qquad\forall v_{T}\in
V_{T}.$

We have

$\displaystyle\underline{a}_{T}(w,w_{\partial T},(v_{T},0))$
$\displaystyle=\int_{T}\boldsymbol{q}_{T}(w,w_{\partial
T})\cdot\boldsymbol{q}_{T}(v_{T},0)+\underbrace{\underline{s}_{T}((w,w_{\partial
T}),(v_{T},0))}_{0}$ $\displaystyle=-\int_{T}\nabla
w\cdot\boldsymbol{q}_{T}(v_{T},0)$

where the stabilization term vanishes due to the polynomial consistency of
$\underline{s}_{T}$ [DPD20, Assumption 2.4 (S3)], and we have used (6.20) for
the last equality. For the remaining term, (3.10) gives

$\int_{T}\boldsymbol{q}_{T}(v_{T},0)\cdot\boldsymbol{p}_{T}-\int_{T}v_{T}\nabla\\!\cdot\\!\boldsymbol{p}_{T}=0\quad\forall\boldsymbol{p}_{T}\in\boldsymbol{W}_{T},$

which we specialize to $\boldsymbol{p}_{T}=-\nabla w$ to infer

$-\int_{T}\boldsymbol{q}_{T}(v_{T},0)\cdot\nabla w+\int_{T}v_{T}\Delta w=0.$

As $w\in\overline{V}^{\textup{c}}_{\ell}$, $\Delta w=0$, which concludes the
proof of $\mathcal{U}\gamma_{\ell}w=w$.

###  6.6. Verification of (HM5)

Remark that $\boldsymbol{q}_{\ell}$ is linear, i.e., for all
$(v_{1},\mu_{1})\in V_{\ell}\times M_{\ell}$ and $(v_{2},\mu_{2})\in
V_{\ell}\times M_{\ell}$,

$\displaystyle\boldsymbol{q}_{\ell}(v_{1}+v_{2},\mu_{1}+\mu_{2})=\boldsymbol{q}_{\ell}(v_{1},\mu_{1})+\boldsymbol{q}_{\ell}(v_{2},\mu_{2}).$

Recalling that
$\underline{u}_{\ell}\coloneq(\mathcal{U}_{\ell}m_{\ell}+\mathcal{V}_{\ell}f,m_{\ell})$,
we then have

$\boldsymbol{\mathcal{Q}}_{\ell}m_{\ell}\overset{\eqref{eq:localQ}}{=}\boldsymbol{q}_{\ell}(\mathcal{U}_{\ell}m_{\ell},m_{\ell})=\boldsymbol{q}_{\ell}(\mathcal{U}_{\ell}m_{\ell}+\mathcal{V}_{\ell}f,m_{\ell})-\boldsymbol{q}_{\ell}(\mathcal{V}_{\ell}f,0)=\boldsymbol{q}_{\ell}\underline{u}_{\ell}-\boldsymbol{q}_{\ell}(\mathcal{V}_{\ell}f,0).$

Consequently,

$\displaystyle\|\boldsymbol{\mathcal{Q}}_{\ell}m_{\ell}-\boldsymbol{\mathcal{Q}}_{\ell}\Pi^{\partial}_{\ell}u\|_{\Omega}=\|\boldsymbol{q}_{\ell}\underline{u}_{\ell}-\boldsymbol{q}_{\ell}(\mathcal{V}_{\ell}f,0)-\boldsymbol{\mathcal{Q}}_{\ell}\Pi^{\partial}_{\ell}u\|_{\Omega}.$

By inserting
$0=\boldsymbol{q}_{\ell}(\Pi^{\textup{d}}_{\ell}u,\Pi^{\partial}_{\ell}u)-\boldsymbol{q}_{\ell}(\Pi^{\textup{d}}_{\ell}u,\Pi^{\partial}_{\ell}u)$
into the right-hand side and using the triangle inequality, we have

$\|\boldsymbol{\mathcal{Q}}_{\ell}m_{\ell}-\boldsymbol{\mathcal{Q}}_{\ell}\Pi^{\partial}_{\ell}u\|_{\Omega}\\\
\leq\|\underbrace{\boldsymbol{q}_{\ell}\underline{u}_{\ell}-\boldsymbol{q}_{\ell}(\Pi^{\textup{d}}_{\ell}u,\Pi^{\partial}_{\ell}u)}_{\text{\normalsize$=:\Xi_{1}$}}\|_{\Omega}+\|\underbrace{\boldsymbol{q}_{\ell}(\Pi^{\textup{d}}_{\ell}u,\Pi^{\partial}_{\ell}u)-\boldsymbol{q}_{\ell}(\mathcal{V}_{\ell}f,0)-\boldsymbol{\mathcal{Q}}_{\ell}\Pi^{\partial}_{\ell}u}_{\text{\normalsize$=:\Xi_{2}$}}\|_{\Omega}.$

The definition of $\underline{a}_{\ell}$ and the energy error estimate (6.8)
readily provides

$\|\Xi_{1}\|_{\Omega}\leq\|\boldsymbol{q}_{\ell}\underline{u}_{\ell}-\boldsymbol{q}_{\ell}(\Pi^{\textup{d}}_{\ell}u,\Pi^{\partial}_{\ell}u)\|_{\underline{a}_{\ell}}\lesssim
h_{\ell}^{\alpha}\|u\|_{\alpha+1}.$

For $\|\Xi_{2}\|_{\Omega}$, we start by combining all terms as follows after
recalling that
$\boldsymbol{\mathcal{Q}}_{\ell}\Pi^{\partial}_{\ell}u\overset{\eqref{eq:localQ}}{=}\boldsymbol{q}_{\ell}(\mathcal{U}_{\ell}\Pi^{\partial}_{\ell}u,\Pi^{\partial}_{\ell}u)$:

$\displaystyle\Xi_{2}$
$\displaystyle=\boldsymbol{q}_{\ell}(\Pi^{\textup{d}}_{\ell}u,\Pi^{\partial}_{\ell}u)-\boldsymbol{q}_{\ell}(\mathcal{V}_{\ell}f,0)-\boldsymbol{q}_{\ell}(\mathcal{U}_{\ell}\Pi^{\partial}_{\ell}u,\Pi^{\partial}_{\ell}u)$
(6.21)
$\displaystyle=\boldsymbol{q}_{\ell}(\underbrace{\Pi^{\textup{d}}_{\ell}u-(\mathcal{U}_{\ell}\Pi^{\partial}_{\ell}u+\mathcal{V}_{\ell}f}_{\text{\normalsize$=:\mathfrak{A}$}}),0).$

Next, we make use of the defining equation (3.10) of $\boldsymbol{q}_{T}$ to
infer

$\int_{T}\boldsymbol{q}_{T}(\mathfrak{A},0)\cdot\boldsymbol{p}_{T}=\int_{T}\mathfrak{A}\nabla\\!\cdot\\!\boldsymbol{p}_{T}\qquad\forall\boldsymbol{p}_{T}\in\boldsymbol{W}_{T}.$

Specializing for $\boldsymbol{p}_{T}=\boldsymbol{q}_{T}(\mathfrak{A},0)$, it
holds that

$\|\Xi_{2}\|_{\Omega}^{2}=\int_{T}\mathfrak{A}\,(\nabla\\!\cdot\\!\boldsymbol{q}_{\ell}(\mathfrak{A},0))\leq\|\mathfrak{A}\|_{\Omega}\|\nabla\\!\cdot\\!\boldsymbol{q}_{\ell}(\mathfrak{A},0)\|_{\Omega},$

where we have used the Cauchy–Schwarz inequality. The inverse inequality (6.2)
used on the last term along with the mesh quasi-uniformity assumption
$h_{T}^{-1}\lesssim h_{\ell}^{-1}$ yields

$\|\Xi_{2}\|_{\Omega}^{2}\leq
h_{\ell}^{-1}\|\mathfrak{A}\|_{\Omega}\|\boldsymbol{q}_{\ell}(\mathfrak{A},0)\|_{\Omega}.$

simplifying, we obtain

$\|\Xi_{2}\|_{\Omega}\leq h_{\ell}^{-1}\|\mathfrak{A}\|_{\Omega}.$

By inserting $0=\mathcal{U}_{\ell}m_{\ell}-\mathcal{U}_{\ell}m_{\ell}$ into
the norm in the right-hand side, recalling the definition (6.21) of
$\mathfrak{A}$, and using triangle inequalities, we get

$\displaystyle\|\Xi_{2}\|_{\Omega}$ $\displaystyle\leq
h_{\ell}^{-1}\left(\|\Pi^{\textup{d}}_{\ell}u-(\mathcal{U}_{\ell}m_{\ell}+\mathcal{V}_{\ell}f)\|_{\Omega}+\|\Pi^{\partial}_{\ell}u-\mathcal{U}_{\ell}m_{\ell}\|_{\Omega}\right)$
$\displaystyle\lesssim
h_{\ell}^{-1}\left(h_{\ell}^{\alpha+1}\|u\|_{\alpha+1}+h_{\ell}^{\alpha+1}\|u\|_{\alpha+1}\right)\lesssim
h_{\ell}^{\alpha}\|u\|_{\alpha+1},$

where we have used the convergence property of the cell unknowns (6.9) for the
first term and the convergence property of the face unknowns (6.10) for the
second term.

###  6.7. Verification of (HM6)

Combining some preliminary results, we have that

$\|\mu\|_{\ell}^{2}\overset{\eqref{EQ:norm_equiv1}}{\lesssim}\|\mu\|_{1,\ell}^{2}\overset{\eqref{EQ:hho_stability}}{\lesssim}a_{\ell}(\mu,\mu)\overset{\eqref{EQ:hho_stability}}{\lesssim}\|\mu\|_{1,\ell}^{2}.$

Thus, we only need to prove that

$\|\mu\|_{1,\ell}\lesssim h_{\ell}^{-1}\|\mu\|_{\ell}.$ (6.22)

Inserting
$0=\boldsymbol{\mathcal{Q}}_{\ell}\mu-\boldsymbol{\mathcal{Q}}_{\ell}\mu$ into
the gradient term of (6.11) and using the triangle inequality gives

$\displaystyle\|\mu\|_{1,\ell}^{2}$
$\displaystyle\lesssim\|\nabla_{\ell}\mathcal{U}_{\ell}\mu+\boldsymbol{\mathcal{Q}}_{\ell}\mu\|_{\Omega}^{2}+\|\boldsymbol{\mathcal{Q}}_{\ell}\mu\|_{\Omega}^{2}+|\mu|_{1,\ell}$
$\displaystyle\lesssim
h_{\ell}^{-2}\|\mathcal{U}_{\ell}\mu-\mu\|_{\ell}^{2}+\|\boldsymbol{\mathcal{Q}}_{\ell}\mu\|_{\Omega}^{2}+|\mu|_{1,\ell}.$
(6.23)

In order to make $\|\cdot\|_{\ell}$ appear in $|\mu|_{1,\ell}$, we observe
that, by mesh quasi-uniformity,

$\displaystyle|\mu|_{1,\ell}$ $\displaystyle\simeq
h_{\ell}^{-2}\sum_{T\in\mathcal{T}_{\ell}}\sum_{F\in\mathcal{F}_{T}}h_{\ell}\|\mathcal{U}_{\ell}\mu-\mu\|_{F}^{2}\overset{\eqref{EQ:norm_l_rewrite}}{\simeq}h_{\ell}^{-2}\|\mathcal{U}_{\ell}\mu-\mu\|_{\ell}^{2}.$
(6.24)

Plugging (6.24) into (6.23) gives

$\displaystyle\|\mu\|_{1,\ell}^{2}$ $\displaystyle\lesssim
h_{\ell}^{-2}\|\mathcal{U}_{\ell}\mu-\mu\|_{\ell}^{2}+\|\boldsymbol{\mathcal{Q}}_{\ell}\mu\|_{\Omega}^{2}\overset{\eqref{EQ:LS1}}{\lesssim}\|\boldsymbol{\mathcal{Q}}_{\ell}\mu\|_{\Omega}^{2}.$
(6.25)

Finally, (6.22) follows from (HM2).

###  6.8. Verification of (HM7)

This equality corresponds to [DPD20, Assumption 2.4 (S3)].

###  6.9. Verification of (HM8)

By the definition of $a_{\ell}$, we have

$\displaystyle s_{\ell}(\mu,\mu)\leq
a_{\ell}(\mu,\mu)\overset{\eqref{EQ:hho_stability}}{\lesssim}\|\mu\|_{1,\ell}^{2}\lesssim\|\boldsymbol{\mathcal{Q}}_{\ell}\mu\|_{\Omega}^{2},$

where the last inequality is (6.25), which uses (HM1).

##  7\. Injection operators for HHO

Consider two successive levels $\ell$ (fine) and $\ell-1$ (coarse). Given the
coarse faces $\mathcal{F}_{\ell-1}$, the mesh nestedness allows us to
decompose $\mathcal{F}_{\ell}$ as the disjoint union
$\widehat{\mathcal{F}}_{\ell}\cup\mathring{\mathcal{F}}_{\ell}$, where

$\widehat{\mathcal{F}}_{\ell}\coloneq\\{F\in\mathcal{F}_{\ell}\mid\exists
F_{\ell-1}\in\mathcal{F}_{\ell-1}\textrm{ s.t. }F\subset
F_{\ell-1}\\},\qquad\mathring{\mathcal{F}}_{\ell}\coloneq\mathcal{F}_{\ell}\setminus\widehat{\mathcal{F}}_{\ell}.$

In the following, we introduce three injection operators used in [DPHM+21a],
denoted by $I_{\ell}^{i}$, $i\in\\{1,2,3\\}$. The first one is defined as

$(I_{\ell}^{1}\mu)_{|F}\coloneq\begin{cases}\mu_{|F}&\text{ if
}F\in\widehat{\mathcal{F}}_{\ell},\\\ (\mathcal{U}_{\ell-1}\mu)_{|F}&\text{
otherwise.}\end{cases}$ (7.1)

This first operator exploits the nestedness of the meshes to straightforwardly
transfer values from coarse faces to their embedded fine faces. Regarding the
fine faces that are _not_ geometrically included in the coarse skeleton (i.e.
$\mathring{\mathcal{F}}_{\ell}$), we make use of the local solver:
$\mathcal{U}_{\ell-1}$ builds a bulk function in the coarse cells, whose
traces provide admissible approximations on the fine faces.

Instead of the straight injection for the embedded faces, an alternative is to
also use the bulk functions reconstructed from the local solver. As
$\mathcal{U}_{\ell-1}$ yields a discontinuous polynomial, and the faces of
$\widehat{\mathcal{F}}_{\ell}$ are located at the interface of two coarse
cells, we propose to take the average of the respective traces, i.e.

$(I_{\ell}^{2}\mu)_{|F}\coloneq\\{\\!\\{\mathcal{U}_{\ell-1}\mu\\}\\!\\}_{F}.$
(7.2)

Remark that this formula also holds for $F\in\mathring{\mathcal{F}}_{\ell}$,
as the average trace of a continuous bulk function on both sides of a face
reduces to its regular trace. This injection operator comes with advantages.
First, it allows the mesh nestedness condition to be relaxed, paving the way
to generalized multigrid methods on non-nested meshes, as developed in
[DPHM+21b]. Additionally, a more complex formula can handle difficulties
occurring at cell interfaces. Typically, using an adequate weighted average
formula can yield robust convergence in the presence of large jumps in the
diffusion coefficient; see [DPHM+21a].

Finally, the third injection operator leverages the higher-order
reconstruction operator, a salient feature of the HHO methods:

$(I_{\ell}^{3}\mu)_{|F}\coloneq\pi_{\ell}^{p}\\{\\!\\{\theta_{\ell-1}^{p+1}(\mathcal{U}_{\ell-1}\mu,\mu)\\}\\!\\}_{F},$
(7.3)

where $\pi_{\ell}^{p}\colon L^{2}(\mathcal{F}_{\ell})\to M_{\ell}$ denotes the
$L^{2}$-orthogonal projector onto the skeletal polynomial space of degree $p$.
$I_{\ell}^{3}$ is based on the same principle of the averaged trace as
$I_{\ell}^{2}$, except that the higher-order reconstruction operator
$\theta_{\ell-1}^{p+1}$ enables the gain of one extra polynomial degree in the
approximation of the coarse error. Then, after computing the average trace,
the polynomial degree is lowered back to its original value by applying the
$L^{2}$-orthogonal projector onto the lower-order space.

To fit these injection operators into the framework, we demonstrate

######  Lemma 7.1.

Under the assumptions (HM2) and (HM4), the injection operators $I_{\ell}^{i}$,
$i\in\\{1,2,3\\}$ verify (IA1)–(IA2).

###### Proof.

Regarding $I_{\ell}^{1}$, we refer to [LRK22a, Lem. 3.2], as they investigate
the same injection operator: $I_{\ell}^{1}$ is their third operator.

Let us consider $I_{\ell}^{2}$: condition (IA2) follows directly from (HM4)
recalling that we assume here $p\geq 1$ (cf. Remark 8.1 below). The discrete
trace inequality (6.1) along with mesh regularity, i.e., $\frac{|T|}{|\partial
T|}\|\mu\|_{\partial T}^{2}\lesssim\|\mu\|_{T}^{2}$ for all
$T\in\mathcal{T}_{\ell}$, gives

$\|\mathcal{U}_{\ell-1}\mu\|_{\ell}\lesssim\|\mathcal{U}_{\ell-1}\mu\|_{\Omega}\lesssim\|\mu\|_{\ell-1},$

where the second inequality is (HM2). (IA1) now follows from the mean value
property $\|\\{\\!\\{\mathcal{U}_{\ell-1}\mu\\}\\!\\}_{F}\|_{\ell}\leq
2\|\mathcal{U}_{\ell-1}\mu\|_{\ell}$.

Let us consider (IA2) for $I_{\ell}^{3}$. On one hand, (HM4) says that

$-\boldsymbol{\mathcal{Q}}_{\ell-1}\gamma_{\ell-1}w=\nabla_{\ell}\theta_{\ell-1}^{p+1}(\mathcal{U}_{\ell-1}\gamma_{\ell-1}w,\gamma_{\ell-1}w)=\nabla
w$

for $w\in\overline{V}^{\textup{c}}$. This, in turn, implies that

$\theta_{\ell-1}^{p+1}(\mathcal{U}_{\ell-1}\gamma_{\ell-1}w,\gamma_{\ell-1}w)=w+c$

for some constant $c\in\mathbb{R}$. On the other hand, the high-order
reconstruction guarantees (cf. (3.11)) that the mean value of
$\theta_{\ell-1}^{p+1}(\mathcal{U}_{\ell-1}\gamma_{\ell-1}w,\gamma_{\ell-1}w)$
equals the mean values of $\mathcal{U}_{\ell-1}\gamma_{\ell-1}w=w$ (through
(HM4)), which implies that $c=0$.

Next, we prove (IA1) for $I_{\ell}^{3}$: Let us consider a coarse element $T$.
We have that

$\displaystyle\sum_{T\in\mathcal{T}_{\ell-1}}\sum_{\bar{}T\in\mathcal{T}_{\ell}}^{\bar{}T\subset
T}$
$\displaystyle\sqrt{\frac{|\bar{}T|}{|\partial\bar{}T|}}\|\pi_{\ell}^{p}\theta_{\ell-1}^{p+1}(\mathcal{U}_{\ell-1}\mu|_{T},\mu|_{\partial\bar{}T})\|_{\partial
T}\lesssim\sum_{T\in\mathcal{T}_{\ell-1}}\|\theta_{\ell-1}^{p+1}(\mathcal{U}_{\ell-1}\mu|_{T},\mu|_{\partial
T})\|_{T}$
$\displaystyle\leq\sum_{T\in\mathcal{T}_{\ell-1}}\|\theta_{\ell-1}^{p+1}(\mathcal{U}_{\ell-1}\mu|_{T},\mu|_{\partial
T})-\mathcal{U}_{\ell-1}\mu|_{T}\|_{T}+\|\mathcal{U}_{\ell-1}\mu|_{T}\|_{\Omega}$
$\displaystyle\lesssim\sum_{T\in\mathcal{T}_{\ell-1}}h_{T}\|\underbrace{\nabla\theta_{\ell-1}^{p+1}(\mathcal{U}_{\ell-1}\mu|_{T},\mu|_{\partial
T})}_{=-\boldsymbol{\mathcal{Q}}\mu}-\nabla\mathcal{U}_{\ell-1}\mu|_{T}\|_{T}+\|\mathcal{U}_{\ell-1}\mu|_{T}\|_{\Omega}$
$\displaystyle\lesssim
h_{T}\|\mu\|_{a_{\ell-1}}+\|\mathcal{U}_{\ell-1}\mu|_{T}\|_{\Omega}\lesssim\|\mu\|_{\ell-1}.$

Here, the first inequality uses the boundedness of the $L^{2}$-orthogonal
projector along with discrete trace inequalities. The second passage follows
inserting $\pm\mathcal{U}_{\ell-1}\mu$ into the norm and using a triangle
inequality. The third inequality is obtained using a local Poincaré–Wirtinger
inequality. The fourth passage is obtained using (HM3) and (HM1). The last
inequality uses (HM6) and (HM2). ∎

##  8\. Numerical experiments

###  8.1. Experimental setup

The numerical tests reported in this section are performed on the Poisson
problem (2.1), on two- and three-dimensional domains. Namely, the fully
elliptic problem shall be tested on the unit square and cube, and the low-
regularity problem shall be tested on the L-shape domain. The problem is
discretized using the standard HHO method as described in Section 3, with the
classical stabilizing term [DPD20, Ex. 2.7 (2.22)]. Tests will span the
polynomial degrees $p=1,2,3$, leaving the special case of $p=0$ out of the
scope, as it does not verify the hypotheses of our framework; see Remark 8.1.
Modal polynomial bases are used in cells and on faces to assemble the method
(specifically, we use $L^{2}$-orthogonal Legendre bases).

######  Remark 8.1.

(The lowest order case) The HHO method with $p=0$ does not verify the error
estimates (6.9) and (6.10). Consequently, (HM5) does not hold and the
multigrid method might not be uniformly convergent. It is worth noting that
this limitation is not just theoretical. It has been observed in practice in
[DPHM+21a, DPHM+21b], where suboptimal convergence is reported for $p=0$ while
the method exhibits optimal convergence for the higher orders, with no
variation in any other parameter.

The multigrid method is constructed as described in Section 4.1. Section 7
describes the injection operators we consider. The V-cycle is carried out
using pointwise Gauss–Seidel smoothing iterations, arranged so that both the
smoothing step and the multigrid iteration remain symmetric. Namely, for
V(1,1), the pre-smoothing iteration is performed in the forward order and the
post-smoothing iteration in the backward order. For V(2,2), the pre- and post-
smoothing procedures involve a single forward iteration followed by a single
backward iteration. The grid hierarchy is built by successive refinements of
an initial simplicial mesh, and the condensed problem is assembled at every
level. The coarsest system is solved with the Cholesky factorization. The
stopping criterion relies on the backward error
$\|\mathbf{r}\|_{2}/\|\mathbf{b}\|_{2}$, where $\mathbf{r}$ denotes the
residual of the algebraic system, $\mathbf{b}$ the right-hand side, and
$\|\cdot\|_{2}$ the standard Euclidean norm applied to the vector space of
coordinates. In all experiments, convergence is considered to be reached once
$\|\mathbf{r}\|_{2}/\|\mathbf{b}\|_{2}<10^{-6}$ is satisfied.

###  8.2. Numerical tests

####  8.2.1. Full regularity: the unit square

The Poisson problem is set up for the manufactured solution
$u\colon(x,y)\mapsto\sin(4\pi x)\sin(4\pi y)$. The unit square is discretized
by a hierarchy of 7 structured triangular meshes. For each problem solved on
this hierarchy, the number of face unknowns is given in Table 1. The largest
problem size we consider involves three million unknowns. Considering the
injection operator $I_{\ell}^{1}$, we first state that the V(1,1) cycle
diverges regardless of the value of $p$. One must raise the number of
smoothing steps to achieve convergence. With the V(2,2) cycle, the method
converges for $p\in\\{1,2\\}$ and still diverges for $p=3$. The numbers of
iterations are reported in Table 2 (the symbol $\infty$ indicating divergence
or a number of iterations $>100$). Their mild increase as the number of levels
grows indicates asymptotic optimality. Table 3 now considers the injection
operator $I_{\ell}^{2}$. In V(1,1), contrary to the same cycle with
$I_{\ell}^{1}$, the method now converges for $p\in\\{1,2\\}$, indicating
better robustness of $I_{\ell}^{2}$ compared to $I_{\ell}^{1}$. The number of
iterations of the V(2,2) cycle illustrates the uniform convergence of the
method for all values of $p$. Finally, using the injection operator
$I_{\ell}^{3}$, we see from the results of Table 4 that the method converges
uniformly with both cycles and for all values of $p$, making this injection
operator the most robust amongst those considered here.

Levels | 3 | 4 | 5 | 6 | 7
---|---|---|---|---|---
$p=1$ | $6016$ | $24320$ | $97792$ | $392192$ | $1570816$
$p=2$ | $9024$ | $36480$ | $146688$ | $588288$ | $2356224$
$p=3$ | $12032$ | $48640$ | $195584$ | $784384$ | $3141632$
Table 1. For the square domain, number of face unknowns at each level and for each value of $p$. Levels | 3 | 4 | 5 | 6 | 7
---|---|---|---|---|---
$p=1$ | 13 | 14 | 14 | 15 | 15
$p=2$ | 36 | 38 | 38 | 39 | 39
$p=3$ | $\infty$ | $\infty$ | $\infty$ | $\infty$ | $\infty$
Table 2. In the square domain, number of V(2,2) iterations with the injection operator $I_{\ell}^{1}$. | V(1,1) | V(2,2)
---|---|---
Levels | 3 | 4 | 5 | 6 | 7 | 3 | 4 | 5 | 6 | 7
$p=1$ | 24 | 25 | 25 | 26 | 26 | 13 | 13 | 14 | 14 | 14
$p=2$ | 20 | 22 | 25 | 26 | 27 | 10 | 10 | 11 | 11 | 11
$p=3$ | $\infty$ | $\infty$ | $\infty$ | $\infty$ | $\infty$ | 13 | 13 | 13 | 14 | 14
Table 3. In the square domain, number of iterations with the injection operator $I_{\ell}^{2}$. | V(1,1) | V(2,2)
---|---|---
Levels | 3 | 4 | 5 | 6 | 7 | 3 | 4 | 5 | 6 | 7
$p=1$ | 18 | 18 | 19 | 19 | 20 | 10 | 10 | 11 | 11 | 11
$p=2$ | 17 | 17 | 17 | 17 | 18 | 9 | 10 | 10 | 10 | 10
$p=3$ | 20 | 21 | 21 | 21 | 21 | 11 | 11 | 11 | 11 | 11
Table 4. In the square domain, number of iterations with the injection
operator $I_{\ell}^{3}$.

####  8.2.2. Low regularity: the L-shape domain

The computational domain is now the L-shape domain
$\Omega=(-1,1)^{2}\setminus([0,1]\times[0,-1])$. The source function of the
Poisson problem is set to zero. The Dirichlet boundary condition follows the
manufactured exact solution $u\colon(r,\varphi)\mapsto
r^{2/3}\sin(\frac{2}{3}\varphi)$, where $(r,\varphi)$ represent polar
coordinates of domain points. An unstructured Delaunay triangulation builds
the coarsest mesh, and a classical refinement technique by edge bisection is
employed to build finer meshes successively. Although the theory requires the
variable cycle to conclude, we notice that fixed cycles exhibit uniform
convergence in practice. Therefore, the test results presented here are
obtained with classical V(1,1) and V(2,2) cycles. Similarly to the preceding
test case, the problem sizes are indicated in Table 5. The results, presented
in Tables 6, 7, and 8 are also qualitatively similar. With the least robust
injection operator, $I_{\ell}^{1}$, the V(1,1) cycle diverges for all values
of $p$, and the V(2,2) cycle is uniformly convergent only for $p=1$ (cf. Table
6). With $I_{\ell}^{2}$, V(1,1) provides a uniformly convergent solver up to
$p=2$, while V(2,2) converges uniformly for all $p\in\\{1,2,3\\}$ (cf. Table
7). Finally, using $I_{\ell}^{3}$, Table 8 shows that the convergence is
optimal for both cycles and all values of $p$. Based on these experiments, one
can conclude that the convergence of the multigrid solver does not seem to
suffer from the lower regularity of the solution. However, if we compare the
number of iterations obtained with $I_{\ell}^{2}$ to those obtained with
$I_{\ell}^{3}$, we remark that robustness aside, $I_{\ell}^{3}$ does not
generally yield a faster convergence (see especially V(2,2)). Recall that
$I_{\ell}^{3}$ is obtained from $I_{\ell}^{2}$ by inserting the higher-order
reconstruction operator. Since it is known that sufficient regularity is
required for higher orders to have beneficial effects, the low regularity of
the solution might explain this lack of noticeable improvement. Looking back
at the fully regular problem on the square, one sees that the higher-order
reconstruction did improve convergence in that case (Table 3 vs. Table 4).

Levels | 3 | 4 | 5 | 6 | 7
---|---|---|---|---|---
$p=1$ | $4480$ | $18176$ | $73216$ | $293888$ | $1177600$
$p=2$ | $6720$ | $27264$ | $109824$ | $440832$ | $1766400$
$p=3$ | $8960$ | $36352$ | $146432$ | $587776$ | $2355200$
Table 5. For the L-shape domain, number of face unknowns at each level and for each value of $p$. Levels | 3 | 4 | 5 | 6 | 7
---|---|---|---|---|---
$p=1$ | 15 | 12 | 14 | 15 | 15
$p=2$ | $\infty$ | $\infty$ | $\infty$ | $\infty$ | $\infty$
$p=3$ | $\infty$ | $\infty$ | $\infty$ | $\infty$ | $\infty$
Table 6. In the L-shape domain, number of V(2,2) iterations with the injection operator $I_{\ell}^{1}$. | V(1,1) | V(2,2)
---|---|---
Levels | 3 | 4 | 5 | 6 | 7 | 3 | 4 | 5 | 6 | 7
$p=1$ | 20 | 20 | 20 | 20 | 20 | 11 | 11 | 11 | 11 | 11
$p=2$ | 17 | 17 | 17 | 17 | 17 | 9 | 9 | 9 | 9 | 9
$p=3$ | $\infty$ | $\infty$ | $\infty$ | $\infty$ | $\infty$ | 11 | 11 | 11 | 11 | 11
Table 7. In the L-shape domain, number of iterations with the injection operator $I_{\ell}^{2}$. | V(1,1) | V(2,2)
---|---|---
Levels | 3 | 4 | 5 | 6 | 7 | 3 | 4 | 5 | 6 | 7
$p=1$ | 16 | 16 | 16 | 16 | 16 | 9 | 9 | 10 | 11 | 11
$p=2$ | 17 | 17 | 17 | 17 | 17 | 8 | 9 | 9 | 9 | 9
$p=3$ | 20 | 20 | 20 | 20 | 20 | 10 | 10 | 10 | 10 | 10
Table 8. In the L-shape domain, number of iterations with the injection
operator $I_{\ell}^{3}$.

####  8.2.3. 3D test case: the cubic domain

The unit cube is discretized by a Cartesian grid, where each element is
decomposed into six geometrically similar tetrahedra, following [Bey95, Fig.
9]. Repeating the same procedure on successive levels of embedded Cartesian
grids ensures that the subsequent tetrahedral grids are also embedded. Five
levels are built, and the considered problem sizes are given in Table 9. The
manufactured solution is $u\colon(x,y)\mapsto\sin(4\pi x)\sin(4\pi y)\sin(4\pi
z)$. This problem is more challenging, so we restrict ourselves to the most
efficient injection operator $I_{\ell}^{3}$ and the V(2,2) cycle (V(1,1)
diverges). The convergence results presented in Table 10, while still
exhibiting asymptotically optimal convergence, also show higher numbers of
iterations than in 2D. More numerical tests are available in [DPHM+21a], where
it is shown that better convergence can be achieved for comparable cost
through simple parameter tunings such as the use of cycles with post-smoothing
only, non-alternating directions in the Gauss-Seidel sweeps, or blockwise
Gauss-Seidel smoothers.

Levels | 3 | 4 | 5
---|---|---|---
$p=1$ | $17280$ | $142848$ | $1161216$
$p=2$ | $34560$ | $285696$ | $2322432$
$p=3$ | $57600$ | $476160$ | $3870720$
Table 9. For the cubic domain, number of face unknowns at each level and for each value of $p$. Levels | 3 | 4 | 5
---|---|---|---
$p=1$ | 18 | 23 | 25
$p=2$ | 23 | 24 | 22
$p=3$ | 22 | 23 | 23
Table 10. In the cubic domain, number of V(2,2) iterations with the injection
operator $I_{\ell}^{3}$.

###  8.3. Conclusions

We have constructed and rigorously analyzed a homogeneous multigrid method for
HHO methods, which does not internally change the discretization scheme. To
this end, we generalized the homogeneous multigrid framework for HDG methods
and verified the resulting, more general assumptions for HHO methods. We
verified our analytical findings numerically using several injection
operators.

## Acknowledgements

The work of D. Di Pietro was partially funded by the European Union (ERC
Synergy, NEMESIS, project number 101115663). Views and opinions expressed are
however those of the authors only and do not necessarily reflect those of the
European Union or the European Research Council Executive Agency. Neither the
European Union nor the granting authority can be held responsible for them.

G. Kanschat has been supported by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) under Germany’s Excellence Strategy EXC 2181/1 -
390900948 (the Heidelberg STRUCTURES Excellence Cluster).

A. Rupp has been supported by the Academy of Finland’s grant number 350101
_Mathematical models and numerical methods for water management in soils_ ,
grant number 354489 _Uncertainty quantification for PDEs on hypergraphs_ ,
grant number 359633 _Localized orthogonal decomposition for high-order, hybrid
finite elements_ , Business Finland’s project number 539/31/2023 _3D-Cure: 3D
printing for personalized medicine and customized drug delivery_ , and the
Finnish _Flagship of advanced mathematics for sensing, imaging and modelling_
, decision number 358944.

The authors thank Ulrich Rüde (FAU Erlangen-Nuremberg) for his support.

## References

  * [BDM85] F. Brezzi, J. Douglas, and L.D. Marini.  Two families of mixed finite elements for second order elliptic problems.  Numerische Mathematik, 47(2):217–235, 1985. 10.1007/BF01389710.
  * [BDP21] L. Botti and D. A. Di Pietro.  $p$-multilevel preconditioners for HHO discretizations of the Stokes equations with static condensation.  Communications on Applied Mathematics and Computation, 2021. 10.1007/s42967-021-00142-5.
  * [BdVLM14] L. Beirão da Veiga, K. Lipnikov, and G. Manzini.  The Mimetic Finite Difference Method for Elliptic Problems.  Springer International Publishing, Cham, 2014. 10.1007/978-3-319-02663-3.
  * [BE14] J. Bonelle and A. Ern.  Analysis of compatible discrete operator schemes for elliptic problems on polyhedral meshes.  ESAIM: M2AN, 48(2):553–581, 2014. 10.1051/m2an/2013104.
  * [Bey95] J. Bey.  Tetrahedral grid refinement.  Computing, 55(4):355–378, 1995. 10.1007/BF02238487.
  * [BP92] J.H. Bramble and J.E. Pasciak.  The analysis of smoothers for multigrid algorithms.  Mathematics of Computation, 58(198):467–488, 1992. 10.1090/S0025-5718-1992-1122058-0.
  * [BPX91] J.H. Bramble, J.E. Pasciak, and J. Xu.  The analysis of multigrid algorithms with nonnested spaces or noninherited quadratic forms.  Mathematics of Computation, 56(193):1–34, 1991.
  * [BVEM00] W.L. Briggs, H. Van Emden, and S.F. McCormick.  A multigrid tutorial.  Society for Industrial and Applied Mathematics, 2000.
  * [CDG08] B. Cockburn, B. Dong, and J. Guzmán.  A superconvergent LDG-hybridizable Galerkin method for second-order elliptic problems.  Mathematics of Computation, 77(264):1887–1916, 2008.
  * [CDGT14] B. Cockburn, O. Dubois, J. Gopalakrishnan, and S. Tan.  Multigrid for an HDG method.  IMA Journal of Numerical Analysis, 34(4):1386–1425, 2014.
  * [CDPE16] B. Cockburn, D.A. Di Pietro, and A. Ern.  Bridging the hybrid high-order and hybridizable discontinuous Galerkin methods.  ESAIM: Mathematical Modelling and Numerical Analysis, 50(3):635–650, May 2016. 10.1051/m2an/2015051.
  * [CEP21] M. Cicuttin, A. Ern, and N. Pignet.  Hybrid high-order methods—a primer with applications to solid mechanics.  SpringerBriefs in Mathematics. Springer, Cham, 2021. 10.1007/978-3-030-81477-9.
  * [CGL09] B. Cockburn, J. Gopalakrishnan, and R. Lazarov.  Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems.  SIAM Journal on Numerical Analysis, 47(2):1319–1365, 2009. 10.1137/070706616.
  * [DE06] J. Droniou and R. Eymard.  A mixed finite volume scheme for anisotropic diffusion problems on any grid.  Numerische Mathematik, 105(1):35–71, 2006. 10.1007/s00211-006-0034-1.
  * [DEGH10] J. Droniou, R. Eymard, T. Gallouët, and R. Herbin.  A unified approach to mimetic finite difference, hybrid finite volume and mixed finite volume methods.  Mathematical Models and Methods in Applied Sciences, 20(02):265–295, 2010. 10.1142/S0218202510004222.
  * [DG11] L. Demkowicz and J. Gopalakrishnan.  Analysis of the DPG method for the Poisson equation.  SIAM Journal on Numerical Analysis, 49(5/6):1788–1809, 2011.
  * [DGTZ07] H.Y. Duan, S.Q. Gao, R.C.E. Tan, and S. Zhang.  A generalized BPX multigrid framework covering nonnested V-cycle methods.  Mathematics of Computation, 76(257):137–152, 2007.
  * [DPD20] D.A. Di Pietro and J. Droniou.  The Hybrid High-Order method for polytopal meshes.  Number 19 in Modeling, Simulation and Application. Springer International Publishing, 2020. 10.1007/978-3-030-37203-3.
  * [DPE15] D.A. Di Pietro and A. Ern.  A hybrid high-order locking-free method for linear elasticity on general meshes.  Comput. Meth. Appl. Mech. Engrg., 283:1–21, 2015. 10.1016/j.cma.2014.09.009.
  * [DPEL14] D.A. Di Pietro, A. Ern, and S. Lemaire.  An arbitrary-order and compact-stencil discretization of diffusion on general meshes based on local reconstruction operators.  Comput. Meth. Appl. Math., 14(4):461–472, 2014. 10.1515/cmam-2014-0018.
  * [DPHM+21a] D.A. Di Pietro, F. Hülsemann, P. Matalon, P. Mycek, U. Rüde, and D. Ruiz.  An $h$-multigrid method for hybrid high-order discretizations.  SIAM Journal on Scientific Computing, 43(5):S839–S861, 2021. 10.1137/20M1342471.
  * [DPHM+21b] D.A. Di Pietro, F. Hülsemann, P. Matalon, P. Mycek, U. Rüde, and D. Ruiz.  Towards robust, fast solutions of elliptic equations on complex domains through hybrid high‐order discretizations and non‐nested multigrid methods.  International Journal for Numerical Methods in Engineering, 122(22):6576–6595, 2021. 10.1002/nme.6803.
  * [DPHM+23] D.A. Di Pietro, F. Hülsemann, P. Matalon, P. Mycek, and U. Rüde.  Algebraic multigrid preconditioner for statically condensed systems arising from lowest-order hybrid discretizations.  SIAM Journal on Scientific Computing, 45(4):S329–S350, 2023. 10.1137/21M1429849.
  * [DPMMR23] D.A. Di Pietro, P. Matalon, P. Mycek, and U. Rüde.  High-order multigrid strategies for HHO discretizations of elliptic equations.  Numerical Linear Algebra with Applications, 30(1):e2456, 2023. 10.1002/nla.2456.
  * [EG21] A. Ern and J.L. Guermond.  Finite elements I—Approximation and interpolation, volume 72 of Texts in Applied Mathematics.  Springer, Cham, [2021] ©2021. 10.1007/978-3-030-56341-7.
  * [EGH10] R. Eymard, T. Gallouët, and R. Herbin.  Discretization of heterogeneous and anisotropic diffusion problems on general nonconforming meshes SUSHI: a scheme using stabilization and hybrid interfaces.  IMA Journal of Numerical Analysis, 30(4):1009–1043, 2010. 10.1093/imanum/drn084.
  * [GR86] V. Girault and P.A. Raviart.  Finite Element Methods for Navier-Stokes Equations.  Springer-Verlag, Berlin Heidelberg, 1986. 10.1007/978-3-642-61623-5.
  * [GT09] J. Gopalakrishnan and S. Tan.  A convergent multigrid cycle for the hybridized mixed method.  Numerical Linear Algebra with Applications, 16(9):689–714, 2009\. 10.1002/nla.636.
  * [KW18] M. Kronbichler and W. Wall.  A performance comparison of continuous and discontinuous Galerkin methods with fast multigrid solvers.  SIAM Journal on Scientific Computing, 40(5):A3423–A3448, 2018. 10.1137/16M110455X.
  * [LRK21] P. Lu, A. Rupp, and G. Kanschat.  Homogeneous multigrid for HDG.  IMA Journal of Numerical Analysis, 42:3135–3153, 2021. 10.1093/imanum/drab055.
  * [LRK22a] P. Lu, A. Rupp, and G. Kanschat.  Analysis of injection operators in geometric multigrid solvers for HDG methods.  SIAM Journal on Numerical Analysis, 60(4):2293–2317, 2022. 10.1137/21M1400110.
  * [LRK22b] P. Lu, A. Rupp, and G. Kanschat.  Homogeneous multigrid for embedded discontinuous Galerkin methods.  BIT Numerical Mathematics, 62:1029–1048, 2022. 10.1007/s10543-021-00902-y.
  * [LWKR23] P. Lu, W. Wang, G. Kanschat, and A. Rupp.  Homogeneous multigrid for HDG applied to the Stokes equation.  IMA Journal of Numerical Analysis, 2023. 10.1093/imanum/drad079.
  * [MBTS20] S. Muralikrishnan, T. Bui-Thanh, and J. N. Shadid.  A multilevel approach for trace system in HDG discretizations.  Journal of Computational Physics, 407:109240, 2020. 10.1016/j.jcp.2020.109240.
  * [Mon03] P. Monk.  Finite Element Methods for Maxwell’s Equations.  Oxford University Press, New York, 2003.
  * [PD21] S. Petrides and L. Demkowicz.  An adaptive multigrid solver for DPG methods with applications in linear acoustics and electromagnetics.  Computers & Mathematics with Applications, 87:12–26, 2021. 10.1016/j.camwa.2021.01.017.
  * [RC17] N. V. Roberts and J. Chan.  A geometric multigrid preconditioning strategy for DPG system matrices.  Computers & Mathematics with Applications, 74(8):2018–2043, 2017\.
  * [RT77] P.A. Raviart and J.M. Thomas.  A mixed finite element method for 2-nd order elliptic problems.  In Mathematical Aspects of Finite Element Methods, pages 292–315, Berlin, Heidelberg, 1977. Springer Berlin Heidelberg.
  * [SA17] J. Schütz and V. Aizinger.  A hierarchical scale separation approach for the hybridized discontinuous Galerkin method.  Journal of Computational and Applied Mathematics, 317:500–509, 2017\. 10.1016/j.cam.2016.12.018.
  * [SZ90] L. Scott and S. Zhang.  Finite element interpolation of nonsmooth functions satisfying boundary conditions.  Math. Comp., 54(190):483–493, 1990.
  * [Tan09] S. Tan.  Iterative solvers for hybridized finite element methods.  PhD thesis, University of Florida, 2009.
  * [TOS01] U. Trottenberg, C.W. Oosterlee, and A. Schuller.  Multigrid.  Academic Press, 2001.
  * [WMBT19] T. Wildey, S. Muralikrishnan, and T. Bui-Thanh.  Unified geometric multigrid algorithm for hybridized high-order finite element methods.  SIAM Journal on Scientific Computing, 41(5):S172–S195, 2019. 10.1137/18M1193505.
  * [WY13] J. Wang and X. Ye.  A weak Galerkin finite element method for second-order elliptic problems.  Journal of Computational and Applied Mathematics, 241:103–115, 2013\. 10.1016/j.cam.2012.10.003.