# Homogeneous multigrid for hybrid discretizations: application to HHO methods

Daniele A. Di Pietro, Zhaonan Dong, Guido Kanschat, Pierre Matalon, Andreas Rupp

► To cite this version:

Daniele A. Di Pietro, Zhaonan Dong, Guido Kanschat, Pierre Matalon, Andreas Rupp. Homogeneous multigrid for hybrid discretizations: application to HHO methods. Numerical Methods for Partial Differential Equations, 2025, 41 (5), 10.1002/num.70023. hal-04518103v3

HAL Id: hal-04518103

https://inria.hal.science/hal-04518103v3

Submitted on 24 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
HOMOGENEOUS MULTIGRID FOR HYBRID
DISCRETIZATIONS: APPLICATION TO HHO METHODS

DANIELE A. DI PIETRO, ZHAONAN DONG, GUIDO KANSCHAT, PIERRE MATALON, AND ANDREAS RUPP

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.

# 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].

Date: May 24, 2025.

2010 Mathematics Subject Classification. 65F10, 65N30, 65N50.
This paper addresses the fast solution of the condensed systems arising from hy-
brid discretizations of second-order elliptic equations, focusing on multigrid meth-
ods. 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<sup>+</sup>21a, DPHM<sup>+</sup>21b, BDP21, DPMMR23], but also cov-
ering specifically DPG [RC17, PD21] and RT-H/BDM-H [GT09]. The above ref-
erences focus on geometric multigrid methods, while a fully algebraic multigrid
method has been considered in [DPHM<sup>+</sup>23].

The main difficulty in designing a geometric multigrid algorithm for a trace sys-
tem resides in the fact that the DOFs are supported by the mesh skeleton, which
makes classical intergrid transfer operators designed for element-based DOFs un-
suitable. 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 injec-
tion operators have already been proposed [LRK22a, DPHM<sup>+</sup>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, recent multigrid solvers for HHO methods [DPHM<sup>+</sup>21a, DPHM<sup>+</sup>21b, DPMMR23] have experimentally shown their optimal behavior, but the supporting theory is still missing.

The main theoretical challenges of devising a multigrid theory for HHO methods over a multigrid theory for HDG methods lie in HHO's more abstract stabilization term, which requires a more abstract theory. That is, we need to address two main generalizations to cope with issues that did not come up in [LRK22a]. These are (i) the reduced smoothing property of HHO (assumption (LS1) in [LRK22a] does not hold and needs to be replaced by (HM1) with larger right-hand side) and (ii) the explicit expression of the HDG stabilizer cannot be used in the proofs resulting in more technical proofs while several intermediate results to prove multigrid convergence remain true (or very similar results hold). These issues will be addressed in the first theoretical part. In the second theory 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 since it relies on a compatibility condition with linear finite elements, cf.
(HM4) in Section 3. 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.
Despite our mesh restrictions, we expect this theory to lay the foundation for considering polyhedral meshes that can be decomposed into simplicial meshes for a vast group of hybrid methods. This goal has already been achieved for many other methods, such as mimetic finite differences [AVZ15], (non-hybridized) discontinuous Galerkin [AP19], which do not even require nested mesh sequences. Schwarz-based preconditioners have been derived under similar assumptions for the virtual finite elements in [Cal19]. A related future research prospect lies in addressing multi-grid methods on agglomerated meshes as it has been done for other finite element methods with the ability to cope with polyhedral meshes, such as discontinuous Galerkin [PP22] and virtual finite elements [ABBV23]. As this setting presents further challenges in our case, we leave it to future work. Indeed, since condensed systems arising from hybrid methods rely on face DOFs, it is required that faces be coarsened between levels [DPHM+21a, Sec. 4.4.3], thus making suitable coarsening strategies more difficult to design.

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 multi-grid 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$, $d \ge 1$, be a polytopal bounded Lipschitz domain with boundary $\partial\Omega$. We approximate the unknown $u \in H_0^1(\Omega)$ satisfying

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

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

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

Here, $\| \cdot \|_{\alpha}$ for $\alpha \in \mathbb{R}_{\ge 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_0^{\alpha}(\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 $\mathbf{p}, \mathbf{q} \in [L^2(X)]^d$, $(p, q)_X := \int_X \mathbf{p} \cdot \mathbf{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. Adaptations are needed for other types of boundary conditions (e.g., Neumann or mixed).
### 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,...,L}$ of successively refined simplicial meshes and their corresponding face sets $(\mathcal{F}_\ell)_{\ell=0,...,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,...,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, for face $F$, 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} \quad \forall \ell \in \{1, \dots, L\}. \qquad (3.1)$$

Finally, we assume that our mesh is quasi-uniform implying that there is a constant $c_{\text{uni}}$ such that $h_{\ell} \le c_{\text{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 space for the skeletal variable $m_\ell$ approximating the trace of the unknown $u$ on the skeleton. Moreover, finite-dimensional approximation spaces $V_T \subset L^2(T)$ and $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 $W_T$ is the local test and trial space for the dual unknown $q = -\nabla u$. Global spaces are defined by concatenation via

$$M_{\ell} := \left\{ m \in L^2(\mathcal{F}_{\ell}) \mid \begin{array}{ll} m_{|F} \in \mathcal{P}_p(F) & \forall F \in \mathcal{F}_{\ell} \\ m_{|F} = 0 & \forall F \subset \partial\Omega \end{array} \right\}, \qquad (3.2)$$

$$V_{\ell} := \{v \in L^2(\Omega) \mid v|_{T} \in V_{T} \quad \forall T \in \mathcal{T}_{\ell}\}, \qquad (3.3)$$

$$W_{\ell} := \{q \in L^2(\Omega; \mathbb{R}^d) \mid q|_T \in W_T \quad \forall T \in \mathcal{T}_{\ell}\}, \qquad (3.4)$$

where $\mathcal{P}_p(F)$ comprises polynomials of degree at most $p$ over $F$. We define element-wise, abstract linear operators

$$\begin{align*} \mathcal{U}_T: M_\ell|_{\partial T} &\to V_T, & m_\ell|_{\partial T} &\mapsto u_T, \\ \mathcal{V}_T: L^2(T) &\to V_T, & m_\ell|_{\partial T} &\mapsto u_T, \\ \mathcal{Q}_T: M_\ell|_{\partial T} &\to W_T, & f &\mapsto q_T. \end{align*}$$

The choice of these operators, called local solvers, influences the numerical schemes. Examples for $\mathcal{U}$ and $\mathcal{Q}$, which are relevant to the multigrid method, are provided in Remark 3.1. Global linear operators are constructed by concatenation of their element-local analogous:

$$\mathcal{U}_{\ell}: M_{\ell} \to V_{\ell}, \quad \mathcal{V}_{\ell}: L^2(\Omega) \to V_{\ell}, \quad \mathcal{Q}_{\ell}: M_{\ell} \to W_{\ell}.$$

One can show that $m_\ell \in 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 \quad \text{for all } \mu \in M_{\ell}, \qquad (3.5)$$

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

$$a_{\ell}(m, \mu) = \int_{\Omega} \mathcal{Q}_{\ell} m \cdot \mathcal{Q}_{\ell} \mu + s_{\ell}(m, \mu). \qquad (3.6)$$
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$, $W_T$, the local solvers $U_\ell$, $V_\ell$, $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)$, $W_T = \mathcal{P}_p^d(T)$, and $\tau_\ell > 0$. The operators $U_T$ and $Q_T$ map $m_{\partial T}$ to the element-wise solutions $u_T \in V_T$ and $q_T \in W_T$ of

$$ \int_T \boldsymbol{q}_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} \quad (3.7a) $$

$$ \int_{\partial T} (\boldsymbol{q}_T \cdot \boldsymbol{\nu} + \tau_\ell u_T) v_T - \int_T \boldsymbol{q}_T \cdot \nabla v_T = \tau_\ell \int_{\partial T} m_{\partial T} v_T \quad (3.7b) $$

with test functions $v_T \in V_T$ and $p_T \in W_T$, and $\nu$ the outward normal vector to $\partial T$. In this sense $U_T: m_{\partial T} \mapsto u_T$, $Q_T: m_{\partial T} \mapsto 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{Q}_\ell m_\ell \cdot \boldsymbol{Q}_\ell \mu + \underbrace{\sum_{T \in \mathcal{T}_\ell} \tau_\ell \int_{\partial T} (U_\ell m_\ell - m_\ell)(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 $W_T = \nabla \mathcal{P}_{p+1}(T)$. For $(u_T, m_{\partial T}) \in V_T \times M_{\partial T}$, we define $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} \quad (3.8) $$

for all $p_T \in W_T$. Note that this relation is the same as (3.7a) except for the different choice of $W_T$. In HHO terminology, $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.8) and the closure condition

$$ \int_T \theta_T^{p+1}(u_T, m_{\partial T}) = \int_{\partial T} u_T. \quad (3.9) $$

In a second step, given a local bilinear stabilizer $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)). \quad (3.10) $$

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)) \quad \forall v_T \in V_T. \quad (3.11) $$
(ii) For $f \in L^2(\Omega)$, find $u_T^2 \in V_T$ such that

$$a_T((u_T^2, 0), (v_T, 0)) = \int_T f v_T \quad \forall v_T \in V_T. \qquad (3.12)$$

We define $\mathcal{U}_T: m_{\partial T} \mapsto u_T^1$ solution of (3.11), $\mathcal{V}_T: f \mapsto u_T^2$ solution of (3.12), and we set element-by-element

$$\mathcal{Q}_T: m_{\partial T} \mapsto q_T(\mathcal{U}_T m_{\partial T}, m_{\partial T}) \qquad (3.13)$$

Finally, the global bilinear form is defined as

$$\begin{aligned} a_{\ell}(m_{\ell}, \mu) &=: \sum_{T \in \mathcal{T}_{\ell}} a_{T}((\mathcal{U}_{T} m_{\partial T}, m_{\partial T})), (\mathcal{U}_{T} \mu, \mu)) \\ &= \int_{\Omega} \mathcal{Q}_{\ell} m_{\ell} \cdot \mathcal{Q}_{\ell} \mu + \underbrace{\sum_{T \in \mathcal{T}_{\ell}} s_{T}((\mathcal{U}_{\ell} m_{\ell}, m_{\ell}), (\mathcal{U}_{\ell} \mu, \mu))}_{=s_{\ell}(m_{\ell}, \mu)}. \end{aligned} \qquad (3.14)$$

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} := \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 over $\partial T$. 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) \quad \forall \rho, \mu \in M_{\ell}. \qquad (3.15)$$

The injection operator $I_\ell: 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: M_\ell \to M_\ell$ for preconditioning $A_\ell$ will be defined in Section 4.2.

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

$$P_{\ell-1}^{\langle\rangle}: M_{\ell} \to M_{\ell-1}, \quad \langle P_{\ell-1}^{\langle\rangle} \rho, \mu \rangle_{\ell-1} = \langle \rho, I_{\ell} \mu \rangle_{\ell} \quad \forall \mu \in M_{\ell-1}, \qquad (3.16)$$

$$P_{\ell-1}^{a}: M_{\ell} \to M_{\ell-1}, \quad a_{\ell-1}(P_{\ell-1}^{a} \rho, \mu) = a_{\ell}(\rho, I_{\ell} \mu) \quad \forall \mu \in M_{\ell-1}, \qquad (3.17)$$

replace the $L^2$-orthogonal (for short, $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

$$\Pi_{\ell}^{\partial}: H_{0}^{1}(\Omega) \to M_{\ell}, \quad \langle \Pi_{\ell}^{\partial} v, \mu \rangle_{\ell} = \langle v, \mu \rangle_{\ell} \quad \forall \mu \in M_{\ell},$$

$$\Pi_{\ell}^{d}: H^{1}(\Omega) \to V_{\ell}, \quad (\Pi_{\ell}^{d} v, w)_{\Omega} = (v, w)_{\Omega} \quad \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}, \qquad (\text{HM1}) $$

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

$$ \| \cdot \|_{a_\ell} := a_\ell(\cdot, \cdot)^{\frac{1}{2}}. \qquad (3.18) $$

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

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

* The quantity $\mathcal{Q}_\ell\mu$ approximates $-\nabla_\ell \mathcal{U}_\ell\mu$, where $\nabla_\ell$ denotes the broken gradient:

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

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

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

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

$$ \bar{V}_\ell^c := \{v \in H_0^1(\Omega) \mid v|_{T} \in \mathcal{P}_1(T) \ \forall T \in \mathcal{T}_\ell\}. $$

* 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

$$ \|m_\ell - \Pi_\ell^\partial u\|_{a_\ell} \lesssim h_\ell^\alpha \|u\|_{\alpha+1}. \qquad (\text{HM5}) $$

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

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

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

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

* The global bilinear form is bounded by an HDG-type norm:

$$ a_{\ell}(\mu, \mu) \lesssim \|\mathcal{Q}\mu\|_{\Omega}^{2} + \sum_{T \in \mathcal{T}_{\ell}} \frac{1}{h_{T}^{2}} \|\mathcal{U}\mu - \mu\|_{\partial T}^{2}. \qquad (\text{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) follows directly from the definition of the bilinear form $a_\ell$. Thus, one can use almost identical (but slightly simpler) 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.
Notably, assumption (HM6) typically requires mesh regularity, which might limit the method’s flexibility. Moreover, we emphasize that (HM1) has been weakened compared to our previous work, complicating the analysis below. Moreover, the more general stabilization (as compared to HDG methods) needs to be controlled in our multigrid analysis, which is realized in (HM7) and (HM8).

### 3.3. Assumptions on injection operators.
Our analytical findings will rely on the following conditions to hold on all mesh levels $ℓ$:

(1) Boundedness:

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

(2) Conformity with linear finite elements:

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

(3) Locality: There is a parameter $\bar{\sigma} > 0$ such that for any $T \in \mathcal{T}_{\ell-1}$

$$ (I_{\ell}\mu)|_{\partial T} = (I_{\ell}\rho)|_{\partial T} \qquad (\text{IA3}) $$

for all $\mu, \rho \in M_{\ell-1}$ with $\mu|_{B(x_T, \bar{\sigma}h_{\ell-1})} = \rho|_{B(x_T, \bar{\sigma}h_{\ell-1})}$ and $B(x_T, \bar{\sigma}h_{\ell-1})$ denoting the ball with radius $\bar{\sigma}h_{\ell-1}$ around the barycenter of $T$ called $x_T$.

(The second condition can be interpreted as conformity with overall continuous, piecewise linear finite element spaces, which are nested.) 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 all $w \in \bar{V}_{\ell-1}^{c}$,

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

*Proof.* Let $w \in \bar{V}_{\ell-1}^{c}$ and set $\rho := \gamma_{\ell-1}w$. We have

$$ \mathcal{Q}_{\ell-1}\rho = \mathcal{Q}_{\ell-1}\gamma_{\ell-1}w \stackrel{\text{(HM4)}}{=} -\nabla w \stackrel{\text{(HM4)}}{=} \mathcal{Q}_{\ell}\gamma_{\ell}w \stackrel{\text{(IA2)}}{=} \mathcal{Q}_{\ell}I_{\ell}\gamma_{\ell-1}w = \mathcal{Q}_{\ell}I_{\ell}\rho, \quad (3.19) $$

where we have additionally used the embedding $\bar{V}_{\ell-1}^{c} \subset \bar{V}_{\ell}^{c}$ in the third equality. The definitions of $a_{\ell}$, $a_{\ell-1}$, and (HM7) yield

$$ a_{\ell-1}(P_{\ell-1}^{a}\mu, \rho) = (\mathcal{Q}_{\ell-1}P_{\ell-1}^{a}\mu, \mathcal{Q}_{\ell-1}\rho)_{\Omega} \stackrel{(3.19)}{=} -(\mathcal{Q}_{\ell-1}P_{\ell-1}^{a}\mu, \nabla w)_{\Omega}, $$

$$ a_{\ell}(\mu, I_{\ell}\rho) = (\mathcal{Q}_{\ell}\mu, \mathcal{Q}_{\ell}I_{\ell}\rho)_{\Omega} \stackrel{(3.19)}{=} -(\mathcal{Q}_{\ell}\mu, \nabla w)_{\Omega} $$

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

## 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.5).

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_{\ell}^{i}: 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}: 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_0^{-1}$ be the exact inverse of $A_{\ell}$ on the coarsest level. For $\ell > 0$, set $x^0 = 0 \in M_{\ell}$.

(1) Perform $n_{\ell}$ smoothing steps

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

(2) Perform recursive multigrid step

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

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

(3) Perform $n_{\ell}$ smoothing steps

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

(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

$$\rho_1 n_{\ell} \le n_{\ell-1} \le \rho_2 n_{\ell}, \quad (4.1)$$

which needs to hold for all $\ell \le L$ and uniform constants $1 < \rho_1 \le \rho_2$, see [BPX91].

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_{\ell}^A$ be the largest eigenvalue of $A_{\ell}$, and set

$$K_{\ell} := (\mathbf{1} - (\mathbf{1} - R_{\ell}A_{\ell})(\mathbf{1} - R'_{\ell}A_{\ell}))A_{\ell}^{-1},$$

where $\mathbf{1}$ is the identity matrix and $R'_{\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_{\ell-1}^{a}\mu, \mu)| \le C_1 \left( \frac{\|A_{\ell}\mu\|_{\ell}^2}{\lambda_{\ell}^A} \right)^{\alpha} a_{\ell}(\mu, \mu)^{1-\alpha}. \quad (A1)$$

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

$$|a_{\ell}(\mu - I_{\ell}P_{\ell-1}^{a}\mu, \mu)| \le C_1 \frac{\|A_{\ell}\mu\|_{\ell}^2}{\lambda_{\ell}^A}. \quad (A1')$$

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

$$\|\mu - I_{\ell}P_{\ell-1}^{a}\mu\|_{a_{\ell}} \le C_2 \|\mu\|_{a_{\ell}}. \quad (A2)$$

* Smoothing hypothesis:

$$\frac{\|\mu\|_{\ell}^2}{\lambda_{\ell}^A} \le C_3 \langle K_{\ell}\mu, \mu \rangle_{\ell}. \quad (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 \ge 0$, and for all $\mu \in M_\ell$

$$|a_\ell(\mu - B_\ell A_\ell \mu, \mu)| \le \delta a_\ell(\mu, \mu),$$

where

$$\delta = \frac{C_1 C_3}{n - C_1 C_3} \quad \text{with} \quad n > 2C_1 C_3.$$

*Proof.* This result is [DGTZ07, Thm. 3.1]. □

<figure><img src="image_2.png" /></figure>

**Theorem 4.2.** Let (2.2) hold for $\alpha \in (\frac{1}{2}, 1]$, and further assume (A1), (A2), (A3), and that (4.1) holds with $\rho_1, \rho_2 > 1$. Then, for all $\ell \ge 0$ and all $\mu \in M_\ell$, it holds

$$\eta_0 a_\ell(\mu, \mu) \le a_\ell(B_\ell A_\ell \mu, \mu) \le \eta_1 a_\ell(\mu, \mu)$$

with

$$\eta_0 \ge \frac{n_\ell^\alpha}{M + n_\ell^\alpha}, \quad \eta_1 \le \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 obtained combining [BPX91, Thm. 6] with the relation [BPX91, (A.4) ⇐⇒ (3.4)]. □

<figure><img src="image_3.png" /></figure>

## 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).** In the first half of this section, we prove the energy bounds

$$\|I_\ell P_{\ell-1}^a \mu\|_{a_\ell}^2 \lesssim \|P_{\ell-1}^a \mu\|_{a_{\ell-1}}^2 \lesssim \|\mu\|_{a_\ell}^2 \quad \forall \mu \in M_\ell,$$

which are the key properties needed to prove (A2) (see (5.14) in the final proof of this section). The key step is proving that $\|I_\ell \mu\|_{a_\ell} \lesssim \|\mu\|_{a_{\ell-1}}$ holds for all $\mu \in M_{\ell-1}$ (i.e., the first inequality holds). Proving the second inequality is then a rather straightforward task and can be done in a few lines, cf. (5.13). Thus, the main target of this section can be reduced to proving the first inequality, which is rigorously stated in (5.7).

To this end, we use the averaging linear interpolation

$$I_\ell^{\text{avg}} : V_\ell \to \bar{V}_\ell^c, \quad I_\ell^{\text{avg}} u(x_\alpha) := \{\{u\}\}_\alpha \quad \text{for any internal mesh vertex } \alpha. \quad (5.1)$$

Here, $x_\alpha$ is the point that corresponds to vertex $\alpha$, and $\{\{u\}\}_\alpha$ describes the arithmetic mean of all values that $u$ attains in $\alpha$. For $x \in \partial\Omega$, we set $I_\ell^{\text{avg}} u(x) := 0$.

This averaging operator allows us to go from the coarse mesh (level $\ell - 1$) to the fine mesh (level $\ell$) in way that increases the $\| \cdot \|_{a_\ell}$ norm in a controllable way. Thus, this (auxiliary) averaging operator can be used as a stepping stone to prove the boundedness of the injection operator $I_\ell$. This idea is solidified by the properties in the following lemma.
**Lemma 5.1.** *Assuming (HM1) and (HM3), we have, for all $\ell$, that*

$$ \|\nabla I_{\ell}^{\text{avg}} U_{\ell} \mu\|_{\Omega} \lesssim \|\mu\|_{a_{\ell}}, \quad \forall \mu \in M_{\ell}, \qquad (5.2) $$

$$ \|U_{\ell} \lambda - I_{\ell}^{\text{avg}} U_{\ell} \lambda\|_{\Omega} \lesssim h_{\ell} \|\mu\|_{a_{\ell}}, \quad \forall \lambda \in M_{\ell}, \qquad (5.3) $$

$$ \|\mu - \gamma_{\ell} I_{\ell}^{\text{avg}} U_{\ell} \mu\|_{\ell} \lesssim h_{\ell} \|\mu\|_{a_{\ell}}, \quad \forall \mu \in M_{\ell}. \qquad (5.4) $$

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

Now, we can transport these boundedness arguments to the injection operator in the below lemma. Notice the similarities between (5.2) and (5.5), and between (5.3) and (5.6).

**Lemma 5.2.** Let (IA1), (IA2), (HM1)-(HM4) hold. We have

$$ \|Q_{\ell} I_{\ell} \mu\|_{\Omega} \lesssim \|\mu\|_{a_{\ell-1}}, \qquad (5.5) $$

$$ \|U_{\ell-1} \mu - U_{\ell} I_{\ell} \mu\|_{\Omega} \lesssim h_{\ell-1} \|\mu\|_{a_{\ell-1}} \lesssim h_{\ell} \|\mu\|_{a_{\ell-1}} \qquad (5.6) $$

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

$$ \|I_{\ell} \mu\|_{a_{\ell}} \lesssim \|\mu\|_{a_{\ell-1}}, \quad \forall \mu \in M_{\ell-1}, \qquad (5.7) $$

$$ \|P_{\ell-1}^{a} \mu\|_{a_{\ell-1}} \lesssim \|\mu\|_{a_{\ell}}, \quad \forall \mu \in M_{\ell}. \qquad (5.8) $$

To prove this lemma, we use another operator to compare with: the lifting operator $S_{\ell}: M_{\ell} \to V_{\text{disc}}$ (with $V_{\text{disc}} \subset C(\Omega) \cap L^2(\Omega)$ being a suitable discrete space). This operator serves as an intermediate: We transport the boundedness properties of the averaging operator to the lifting operator, and use the lifting operator to deduce boundedness of the injection operator (one might, loosely speaking, think of a triangle inequality-like argument). It 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. For any $\mu \in M_{\ell}$

$$ (S_{\ell} \mu, v)_{T} = (U_{\ell} \mu, v)_{T} \quad \forall v \in \mathcal{P}_{p}(T), \forall T \in \mathcal{T}_{\ell}, \qquad (5.9a) $$

$$ \langle S_{\ell} \mu, \eta \rangle_{F} = \langle \mu, \eta \rangle_{F} \quad \forall \eta \in \mathcal{P}_{p+1}(F), \forall F \in \mathcal{F}_{\ell}, \qquad (5.9b) $$

$$ S_{\ell} \mu(\mathbf{a}) = \{\{\mu\}\}_{\mathbf{a}} \quad \forall \mathbf{a} \text{ vertex in } \mathcal{T}_{\ell}, \mathbf{a} \notin \partial\Omega, \qquad (5.9c) $$

$$ S_{\ell} \mu(\mathbf{a}) = 0 \quad \forall \mathbf{a} \text{ vertex in } \mathcal{T}_{\ell}, \mathbf{a} \in \partial\Omega \qquad (5.9d) $$

in two dimensions. In three spatial dimensions, we add the constraints

$$ \langle S_{\ell} \mu, \eta \rangle_{\Gamma} = \langle \{\{\mu\}\}_{\Gamma}, \eta \rangle_{\Gamma} \quad \forall \eta \in \mathcal{P}_{p+2}(\Gamma), \qquad (5.9e) $$

Here, $\{\{m\}\}_{\Gamma}$ is the average taken over all cells adjacent to $\Gamma$. For a precise definition of $V_{\text{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]. This construction extends to dimensions higher than three. Moreover, we have

**Lemma 5.3** (Properties of $S_{\ell} \lambda$). Under assumptions (HM1)-(HM4), we have

$$ \|S_{\ell} \mu\|_{\Omega} \cong \|\mu\|_{\ell} \quad \forall \lambda \in M_{\ell}, \qquad (\text{norm equivalence}) \qquad (5.10) $$

$$ S_{\ell} \gamma_{\ell} w = w \quad \forall w \in \bar{V}_{\ell}^{c}, \qquad (\text{lifting identity}) \qquad (5.11) $$

$$ |S_{\ell} \mu|_{1, \Omega} \lesssim \|\mu\|_{a_{\ell}} \quad \forall \lambda \in M_{\ell}, \qquad (\text{lifting bound}) \qquad (5.12) $$

where $| \cdot |_{1, \Omega}$ denotes the $H^1(\Omega)$-seminorm.
*Proof.* The proof can be performed similarly to the proof of [LRK22a, Lem. 5.5] with our modified version of (HM1) as compared to their (LS1). □

With all preparations ready, we can now prove Lemma 5.2 with the key property (5.7):

*Proof of Lemma 5.2.* Inequalities (5.5) and (5.6) can be obtained as in [LRK22a, Lem. 5.1]. Let us obtain inequality (5.7). To this end, we observe that

$$ \|I_{\ell}\mu\|_{a_{\ell}}^{2} \stackrel{(\text{HM8})}{\lesssim} \|\mathcal{Q}_{\ell}I_{\ell}\mu\|_{\Omega}^{2} + \sum_{T \in \mathcal{T}_{\ell}} \frac{1}{h_{T}} \|\mathcal{U}_{\ell}I_{\ell}\mu - \mathcal{U}_{\ell-1}\mu\|_{\partial T}^{2} + \sum_{T \in \mathcal{T}_{\ell}} \frac{1}{h_{T}} \|\mathcal{U}_{\ell-1}\mu - I_{\ell}\mu\|_{\partial T}^{2}, $$

where (5.5) allows us to bound $\|\mathcal{Q}_{\ell}I_{\ell}\mu\|_{\Omega}^{2}$, (5.6) allows us to bound

$$ \sum_{T \in \mathcal{T}_{\ell}} h_{T}^{-1} \|\mathcal{U}_{\ell}I_{\ell}\mu - \mathcal{U}_{\ell-1}\mu\|_{\partial T}^{2} \lesssim h_{\ell}^{-2} \|\mathcal{U}_{\ell}I_{\ell}\mu - \mathcal{U}_{\ell-1}\mu\|_{\Omega}^{2}. $$

Assumption (IA3) implies that $\|\mathcal{U}_{\ell-1}\mu - I_{\ell}\mu\|_{\partial T}$ depends only on $\mu|_{B(x_{T}, \bar{\sigma}h)}$. Thus, we use the face-wise $L^2$ orthogonal projection $\pi: L^2(\tilde{B}(x_{T}, \bar{\sigma}h) \cap \mathcal{F}) \to M_{\tilde{B}(x_{T}, \bar{\sigma}h) \cap \mathcal{F}}$ to define the operator

$$ \mathcal{G}: H^1(\tilde{B}(x_T, \bar{\sigma}h)) \ni v \mapsto \mathcal{U}_{\ell-1}\pi v - I_{\ell}\pi v \in M_{\ell}|_{\partial T} \subset L^2(\partial T), $$

where, $\tilde{B}(x_{T}, \bar{\sigma}h) = B(x_{T}, \bar{\sigma}h) \cap \Omega$. Operator $\mathcal{G}$ is continuous (because of (IA1) and (HM2)) and vanishes if $v$ is an overall continuous, element-wise linear polynomial, c.f. (IA2). Thus, the standard scaling argument, where $B(x_{T}, h)$ is mapped to the unity ball $B(0, 1)$ can be applied. This allows us to deduce that

$$ \sum_{T \in \mathcal{T}_{\ell}} \frac{1}{h_{T}} \|\mathcal{U}_{\ell-1}\pi v - I_{\ell}\pi v\|_{\partial T}^{2} \lesssim |v|_{H^1(\Omega)}^{2}. $$

Finally, setting $v = S_{\ell-1}\mu$ with the lifting operator $S_{\ell-1}$ as defined in (5.9) yields the result (in conjunction with (5.12)).

Relation (5.8) results from

$$ \begin{aligned} \|P_{\ell-1}^a \mu\|_{a_{\ell-1}}^2 &\stackrel{(3.18)}{=} a_{\ell-1}(P_{\ell-1}^a \mu, P_{\ell-1}^a \mu) \stackrel{(3.17)}{=} a_{\ell}(\mu, I_{\ell}P_{\ell-1}^a \mu) \\ &\le \|\mu\|_{a_{\ell}} \|I_{\ell}P_{\ell-1}^a \mu\|_{a_{\ell}} \stackrel{(5.7)}{\lesssim} \|\mu\|_{a_{\ell}} \|P_{\ell-1}^a \mu\|_{a_{\ell-1}}, \end{aligned} \quad (5.13) $$

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

**Lemma 5.4.** *Under the assumptions of Lemma 5.2, (A2) holds.*

*Proof.* Using first the linearity and symmetry of $a_{\ell}$ and then the definition (3.17) of the Ritz quasi-projector $P_{\ell-1}^a$, we obtain, for all $\mu \in M_{\ell}$,

$$ \begin{aligned} & a_{\ell}(\mu - I_{\ell}P_{\ell-1}^{a}\mu, \mu - I_{\ell}P_{\ell-1}^{a}\mu) \\ &= a_{\ell}(\mu, \mu) - 2a_{\ell}(\mu, I_{\ell}P_{\ell-1}^{a}\mu) + a_{\ell}(I_{\ell}P_{\ell-1}^{a}\mu, I_{\ell}P_{\ell-1}^{a}\mu) \\ &\le \|\mu\|_{a_{\ell}}^{2} \underbrace{-2a_{\ell-1}(P_{\ell-1}^{a}\mu, P_{\ell-1}^{a}\mu)}_{\le 0} + \|I_{\ell}P_{\ell-1}^{a}\mu\|_{a_{\ell}}^{2}. \end{aligned} $$

To conclude, we write for the rightmost term

$$ \|I_{\ell}P_{\ell-1}^{a}\mu\|_{a_{\ell}}^{2} \stackrel{(5.7)}{\lesssim} \|P_{\ell-1}^{a}\mu\|_{a_{\ell-1}}^{2} \stackrel{(5.8)}{\lesssim} \|\mu\|_{a_{\ell}}^{2}. \quad (5.14) $$
5.2. **Proofs of (A1) and (A1')**. A classical approach for multigrid proofs consists is considering $A_\ell\mu$ (cf. (3.15)) 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$. To this end, we use the lifting operator $S_\ell$. That is, we define the auxiliary right-hand side $f_\mu \in V_{\text{disc}}$ as the unique solution of

$$ (f_\mu, S_\ell\eta) = \langle A_\ell\mu, \eta \rangle_\ell = a_\ell(\mu, \eta) \quad \forall \eta \in M_\ell, \qquad (5.15) $$

and its perturbed skeletal approximation $\tilde{\mu} \in M_\ell$ with

$$ a_\ell(\tilde{\mu}, \eta) = (f_\mu, U_\ell\eta) \quad \forall \eta \in M_\ell. \qquad (5.16) $$

Notably, both $\mu$ and $\tilde{\mu}$ approximate the solution of the continuous problem

$$ (\nabla \tilde{u}, \nabla v) = (f_\mu, v) \quad \forall v \in H_0^1(\Omega). \qquad (5.17) $$

Thus, they are close to one another, which is formalized in the lemma below:

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

$$ \|\mu - \tilde{\mu}\|_{a_\ell} \lesssim h_\ell \|f_\mu\|_{\Omega}, \quad \text{and} \quad \|f_\mu\|_{\Omega} \lesssim \|A_\ell\mu\|_{\ell}. \qquad (5.18) $$

The key trick in the proof of (A1') and (A1) is the observation that, on one hand, functions from $\bar{V}_{\ell-1}^c$ are not changed by the injection/projection operators or the local solvers (as implied by the quasi orthogonality in Lemma 3.3), while, on the other hand, this function space has approximation properties. Thus, we can introduce an arbitrary $\bar{V}_{\ell-1}^c$ in a scalar product resembling that in Lemma 3.3 and hope to exploit the approximation properties described in Lemma 5.6, whose right-hand side can be controlled by the means of Lemma 5.7.

**Lemma 5.6** (Reconstruction approximation). Let (IA1), (IA2), and (HM1)–(HM5) hold. Assume further that the model problem admits the regularity estimate (2.2). Then, for all $\mu \in M_\ell$, there exists an auxiliary function $\bar{u} \in \bar{V}_{\ell-1}^c$ such that

$$ \| \mathcal{Q}_\ell \mu + \nabla \bar{u} \|_{\Omega} + \| \mathcal{Q}_{\ell-1} P_{\ell-1}^a \mu + \nabla \bar{u} \|_{\Omega} \lesssim h_\ell^\alpha \| f_\mu \|_{\alpha-1}. \qquad (5.19) $$

*Proof.* This is [LRK22a, Lem. 5.8]. □

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

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

*Proof.* By the definition of negative norms and properties of sup,

$$ \|f_\mu\|_{-1} \le \sup_{\psi \in H_0^1(\Omega)} \frac{(f_\mu, \psi - S_\ell \gamma_\ell \bar{\Pi}_\ell^c \psi)}{|\psi|_{1,\Omega}} + \sup_{\psi \in H_0^1(\Omega)} \frac{(f_\mu, S_\ell \gamma_\ell \bar{\Pi}_\ell^c \psi)}{|\psi|_{1,\Omega}}. \qquad (5.20) $$

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

$$ |\bar{\Pi}_\ell^c v|_{1,\Omega} \lesssim |v|_{1,\Omega} \quad \forall v \in H^1(\Omega), \qquad (5.21) $$

$$ \|v - \bar{\Pi}_\ell^c v\|_{\Omega} \lesssim h_\ell^{1-k+\alpha} |v|_{\alpha+1}, \quad \forall v \in H^{\alpha+1}(\Omega), \quad k = 0, 1. \qquad (5.22) $$

An example is given in [SZ90]. We continue writing

$$ \|\psi - S_\ell \gamma_\ell \bar{\Pi}_\ell^c \psi\|_{\Omega} \stackrel{(5.11)}{=} \|\psi - \bar{\Pi}_\ell^c \psi\|_{\Omega} \lesssim h_\ell |\psi|_{1,\Omega}, \qquad (5.23) $$

and observe that

$$ (f_\mu, S_\ell \gamma_\ell \bar{\Pi}_\ell^c \psi)_{\Omega} \stackrel{(5.15)}{=} a_\ell(\mu, \gamma_\ell \bar{\Pi}_\ell^c \psi), $$
which immediately yields, using a Cauchy–Schwarz inequality,

$$|(f_\mu, S_\ell \gamma_\ell \bar{\Pi}_\ell^c \psi)_\Omega| \le \|\mu\|_{a_\ell} \|\gamma_\ell \bar{\Pi}_\ell^c \psi\|_{a_\ell} = \|\mu\|_{a_\ell} |\bar{\Pi}_\ell^c \psi|_{1,\Omega} \lesssim \|\mu\|_{a_\ell} |\psi|_{1,\Omega}. \quad (5.24)$$

Applying a Cauchy–Schwarz inequality to the numerator of the first supremum in (5.20) followed by (5.23) and using (5.24) to estimate the numerator of the second supremum, we get, after simplification,

$$\|f_\mu\|_{-1} \le h_\ell \|f_\mu\|_{\Omega} + \|\mu\|_{a_\ell} \lesssim h_\ell \|A_\ell \mu\|_{\ell} + \|\mu\|_{a_\ell} \lesssim \|\mu\|_{a_\ell}, \quad (5.18)$$

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 (\frac{1}{2}, 1]$, and if the assumptions (IA1), (IA2), and (HM1)–(HM8) hold, then (A1) is satisfied.

Remark 5.9 (Full elliptic regularity). Theorem 5.8 implies (A1') if $\alpha = 1$.

*Proof.* We use the definitions (3.6) of $a_\ell$ and (3.17) of $P_{\ell-1}^a$ to obtain

$$\begin{aligned} a_\ell(\mu - I_\ell P_{\ell-1}^a \mu, \mu) &= a_\ell(\mu, \mu) - a_{\ell-1}(P_{\ell-1}^a \mu, P_{\ell-1}^a \mu) \\ &= \underbrace{(\mathcal{Q}_\ell \mu, \mathcal{Q}_\ell \mu)_\Omega - (\mathcal{Q}_{\ell-1} P_{\ell-1}^a \mu, \mathcal{Q}_{\ell-1} P_{\ell-1}^a \mu)_\Omega}_{\mathcal{T}_1} \\ &\quad + \underbrace{s_\ell(\mu, \mu) - s_{\ell-1}(P_{\ell-1}^a \mu, P_{\ell-1}^a \mu)}_{\mathcal{T}_2}, \end{aligned}$$

and analyze the respective contributions $\mathcal{T}_1$ and $\mathcal{T}_2$ separately.

First, binomial factorization yields

$$\mathcal{T}_1 = (\mathcal{Q}_\ell \mu + \mathcal{Q}_{\ell-1} P_{\ell-1}^a \mu, \mathcal{Q}_\ell \mu - \mathcal{Q}_{\ell-1} P_{\ell-1}^a \mu)_\Omega.$$

Let now $w \in \bar{V}_{\ell-1}^c$. Invoking the quasi-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

$$\begin{aligned} \mathcal{T}_1 &= ((\mathcal{Q}_\ell \mu + \nabla w) + (\mathcal{Q}_{\ell-1} P_{\ell-1}^a \mu + \nabla w), (\mathcal{Q}_\ell \mu + \nabla w) - (\nabla w + \mathcal{Q}_{\ell-1} P_{\ell-1}^a \mu))_\Omega \\ &\le (\|\mathcal{Q}_\ell + \nabla w\|_\Omega + \|\mathcal{Q}_{\ell-1} P_{\ell-1}^a \mu + \nabla w\|_\Omega)^2 \stackrel{(5.19)}{\lesssim} h_\ell^{2\alpha} \|f_\mu\|_{\alpha-1}^2, \end{aligned} \quad (5.25)$$

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

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

$$s_\ell(\mu, \mu) \stackrel{(\text{HM}7)}{=} s_\ell(\mu - \gamma_\ell w, \mu - \gamma_\ell w) \le a_\ell(\mu - \gamma_\ell w, \mu - \gamma_\ell w) = \|\mu - \gamma_\ell w\|_{a_\ell}^2.$$

Using the triangle inequality and Lemma 5.5, we have that

$$\|\mu - \gamma_\ell w\|_{a_\ell} \le \|\mu - \tilde{\mu}\|_{a_\ell} + \|\tilde{\mu} - \gamma_\ell w\|_{a_\ell} \le h_\ell \|f_\mu\|_{\Omega} + \|\tilde{\mu} - \gamma_\ell w\|_{a_\ell}$$

Defining $w$ as the $L^2$ orthogonal projection of $\tilde{u}$ into $\bar{V}_\ell^c$ yields the result with (HM5) as we can write

$$\|\tilde{\mu} - \gamma_\ell w\|_{a_\ell} \le \|\tilde{\mu} - \Pi_\ell^\partial \tilde{u}\|_{a_\ell} + \|\Pi_\ell^\partial \tilde{u} - \gamma_\ell w\|_{a_\ell},$$
where the first summand is bounded via (HM5) and the second summand can be bounded using (HM6) and standard approximation properties. 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. \qquad (5.26) $$

Using (5.25) and (5.26) to estimate the right-hand side of (5.25) and invoking Sobolev interpolation provides us with

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

and Lemmas 5.5 and 5.7 yield

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

This is the result, since (HM6) implies that $\lambda_{\ell}^{A} \lesssim h_{\ell}^{-2}$. $\square$

<figure><img src="image_4.png" /></figure>

## 6. VERIFICATION OF ASSUMPTIONS FOR HHO

6.1. Preliminaries. For all $T \in \mathcal{T}_{\ell}$, we introduce the set $\mathcal{F}_T$ collecting the faces of $T$ and denote a polynomial space on $T$ by $\mathcal{P}(T)$.

6.1.1. General properties of used norms and function spaces. Let $T \in \mathcal{T}_{\ell}$, $F \in \mathcal{F}_T$. We recall the following results:

* The discrete trace inequality (see, e.g., [DPD20, Lem. 1.32]):

$$ \|v\|_F \lesssim h_T^{-\frac{1}{2}} \|v\|_T \quad \forall v \in \mathcal{P}(T). \qquad (6.1) $$

* The discrete inverse inequality (see, e.g., [DPD20, Lem. 1.28]):

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

* The Poincaré-Friedrichs inequality (see, e.g., [EG21, Lem. 3.30]):

$$ \|\rho\|_T \lesssim h_T \|\nabla \rho\|_T + h_T^{\frac{1}{2}} \|\rho\|_{\partial T} \quad \forall \rho \in H^1(T). \qquad (6.3) $$

* The norm equivalence

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

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 $a_{\ell}$ introduced in (3.10), we consider the hybrid, discrete problem: Find $(u_{\ell}, m_{\ell}) \in V_{\ell} \times M_{\ell}$ such that

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

The pair $a_{\ell} := (U_{\ell}m_{\ell} + 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.5) (see [CDPE16, Prop. 4]). In particular, $U_{\ell}m_{\ell} + 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

$$|(v, \mu)|_{1, \ell}^2 := \sum_{T \in \mathcal{T}_\ell} |(v, \mu)|_{1, T}^2 = \sum_{T \in \mathcal{T}_\ell} \sum_{F \in \mathcal{F}_T} h_F^{-1} \|v - \mu\|_F^2, \quad (6.6a)$$

$$\|(v, \mu)\|_{1, \ell}^2 := \sum_{T \in \mathcal{T}_\ell} \|(v, \mu)\|_{1, T}^2 = \sum_{T \in \mathcal{T}_\ell} \left( \|\nabla v\|_T^2 + |(v, \mu)|_{1, T}^2 \right), \quad (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 $s_T$ hold true. Then, we have:

* (Boundedness and stability of $a_T$ [CEP21, Lem. 2.6])

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

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

$$\|\boldsymbol{q}u_{\ell} - \boldsymbol{q}(\Pi_{\ell}^{\mathrm{d}}u, \Pi_{\ell}^{\mathrm{\theta}}u)\|_{a_{\ell}} \lesssim h_{\ell}^{\alpha} \|u\|_{\alpha+1}, \quad (6.8)$$

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

* (Cell unknown $L^2$-error estimate [CEP21, Lem. 2.11])

If $p \ge 1$, $\alpha > \frac{1}{2}$, then

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

where $\delta := \min\{\alpha, 1\}$.

* (Face unknown $L^2$-error estimate [CEP21, Lem. 2.9])

If $p \ge 1$, then

$$\|m_{\ell} - \Pi_{\ell}^{\partial} u\|_{\ell} \lesssim h_{\ell}^{\alpha+1} \|u\|_{\alpha+1}. \quad (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 \|_{1, \ell}$ and $| \cdot |_{1, \ell}$ defined in (6.6):

$$\|\mu\|_{1, \ell} := \|(\mathcal{U}_\ell \mu, \mu)\|_{\underline{1}, \ell} \quad \text{and} \quad |\mu|_{1, \ell} := |(\mathcal{U}_\ell \mu, \mu)|_{\underline{1}, \ell}. \quad (6.11)$$

Similarly, the condensed bilinear form $a_\ell$ introduced in (3.14) is built from the hybrid bilinear form $a_\ell$ (3.10), in which the generic cell unknown variable is also recovered from the skeletal variable through $\mathcal{U}_\ell$, i.e.

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

In the same fashion, we also have

$$s_{\ell}(m, \mu) = 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 $a_\ell$ and $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. \quad (6.12)$$

* Since $s_\ell$ is symmetric positive semi-definite [DPD20, Assumption 2.4], so is $s_\ell$.

**Lemma 6.1.** The norm $\| \cdot \|_\ell$ is weaker than $\| \cdot \|_{1, \ell}$, i.e.,

$$\|\mu\|_{\ell} \lesssim \|\mu\|_{1, \ell} \quad \forall \mu \in M_{\ell}. \quad (6.13)$$
*Proof.* Starting from (6.4), the triangle inequality is used via the insertion of 0 = $U\mu - U\mu$ into the norms in the sum, yielding

$$
\begin{align*}
\|\mu\|_{\ell}^{2} &\lesssim \sum_{F \in \mathcal{F}} h_{F} (\|U_{\ell}\mu - \mu\|_{F}^{2} + \|U_{\ell}\mu\|_{F}^{2}) \\
&\lesssim \sum_{F \in \mathcal{F}} h_{F} \|U_{\ell}\mu - \mu\|_{F}^{2} + \sum_{T \in \mathcal{T}_{\ell}} \|U_{\ell}\mu\|_{T}^{2} \tag{6.14} \\
&\lesssim \sum_{F \in \mathcal{F}} h_{F} \|U_{\ell}\mu - \mu\|_{F}^{2} + \|(U_{\ell}\mu, \mu)\|_{1,\ell}^{2},
\end{align*}
$$

where the last passage follows from the discrete Poincaré inequality [DPD20, Lem. 2.15] applied to $(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} \|U_{\ell}\mu - \mu\|_{F}^{2} + \|\nabla_{\ell} U_{\ell}\mu\|_{\Omega}^{2} + |\mu|_{1,\ell}^{2}. \quad (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$. Hence, we can go on writing

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

and the result follows by taking the square root.

<figure><img src="image_5.png" /></figure>

6.2. Verification of (HM1). We write

$$
\begin{align*}
\|U_{\ell}\mu - \mu\|_{\ell}^{2} &\lesssim \sum_{T \in \mathcal{T}_{\ell}} \sum_{F \in \mathcal{F}_{T}} h_{F} \|U_{\ell}\mu - \mu\|_{F}^{2} \\
&\leq h_{\ell}^{2} \sum_{T \in \mathcal{T}_{\ell}} \sum_{F \in \mathcal{F}_{T}} h_{F}^{-1} \|U_{\ell}\mu - \mu\|_{F}^{2} \\
&\stackrel{(6.6b)}{\lesssim} h_{\ell}^{2} \|(U_{\ell}\mu, \mu)\|_{\underline{1},\underline{h}}^{2} \stackrel{(6.7)}{\lesssim} h_{\ell}^{2} \|(U_{\ell}\mu, \mu)\|_{\underline{a}_{\ell}}^{2} \stackrel{(3.18)}{=} h_{\ell}^{2} \|\mu\|_{\underline{a}_{\ell}}^{2},
\end{align*}
$$

where the first inequality follows from mesh regularity implying $h_T \lesssim h_F$. 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.11) of $U_\ell$ and choosing $v_T = u_T^1 = U_T\mu$ and $m_{\partial T} = \mu$, we have, denoting by $\| \cdot \|_{\underline{a}_T}$ the seminorm induced by $\underline{a}_T$,

$$
\begin{align*}
\|(U_T \mu, 0)\|_{\underline{a}_T}^2 &= \underline{a}_T((U_T \mu, 0), (U_T \mu, 0)) \\
&\stackrel{(3.11)}{=} -\underline{a}_T((0, \mu), (U_T \mu, 0)) \le \|(0, \mu)\|_{\underline{a}_T} \|(U_T \mu, 0)\|_{\underline{a}_T},
\end{align*}
$$

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

$$
\|(U_T \mu, 0)\|_{\underline{a}_T} \le \|(0, \mu)\|_{\underline{a}_T}. \quad (6.16)
$$

We write

$$
\| (U_T \mu, 0) \|_{1,T}^{2} \stackrel{(6.7)}{\lesssim} \| (U_T \mu, 0) \|_{\underline{a}_T}^{2} \stackrel{(6.16)}{\lesssim} \| (0, \mu) \|_{\underline{a}_T}^{2} \stackrel{(6.7)}{\lesssim} \| (0, \mu) \|_{1,T}^{2} \lesssim h_T^{-1} \| \mu \|_{\partial T}^{2}, \quad (6.17)
$$
where the last inequality follows recalling the definition (6.6b) of $\|\cdot\|_{1,T}$ and noticing that $h_T^{-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)\|_{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)\|_{\frac{1}{T}, T}^2 \stackrel{(6.17)}{\lesssim} h_T \|\mu\|_{\partial T}^2. \quad (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$. Re-
calling the definition (3.8) of $\mathbf{q}_T(\cdot, \cdot)$ with $(u_T, m_{\partial T}) = (\mathcal{U}_T \mu, \mu)$ and choosing
$\mathbf{p}_T = \mathbf{q}_T(\mathcal{U}_T \mu, \mu) \stackrel{(3.13)}{=} \mathcal{Q}_T \mu$, $u_T = \mathcal{U}_T \mu$, and $m_{\partial T} = \mu$, we have

$$ \int_T \mathcal{Q}_T \mu \cdot \mathcal{Q}_T \mu = \int_T \mathcal{U}_T \mu (\nabla \cdot \mathcal{Q}_T \mu) - \int_{\partial T} \mu (\mathcal{Q}_T \mu \cdot \mathbf{\nu}). $$

Estimating the right-hand side of the above expression with Cauchy–Schwarz in-
equalities followed by the discrete inverse inequality (6.2) for the first term, the
discrete trace inequality (6.1) for the second term, we have after simplifying and
raising to the square,

$$ \|\mathcal{Q}_T \mu\|_T^2 \lesssim h_T^{-2} \|\mathcal{U}_T \mu\|_T^2 + h_T^{-1} \|\mu\|_{\partial T}^2 \stackrel{(6.18)}{\le} h_T^{-1} \|\mu\|_{\partial T}^2. \quad (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 $\|\mathcal{Q}_\ell \mu\|_{\Omega} \lesssim h_\ell^{-1} \|\mu\|_{\ell}$.

6.4. Verification of (HM3). Plugging the definition (3.13) of $\mathcal{Q}_T$ into (3.8) we have,

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

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

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

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

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

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

$$ \|\mathcal{Q}_T \mu + \nabla \mathcal{U}_T \mu\|_T \le \|\mathcal{U}_T \mu - \mu\|_{\partial T} \|\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

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

Simplifying and squaring, we get

$$ \|\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

$$
\begin{aligned}
\| \mathbf{Q}_T \mu + \nabla \mathcal{U}_T \mu \|_{\Omega}^2 &\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 \|_F^2 \\
&\lesssim h_{\ell}^{-2} \sum_{T \in \mathcal{T}_{\ell}} \sum_{F \in \mathcal{F}_T} h_F \| \mathcal{U}_T \mu - \mu \|_F^2 \\
&\lesssim h_{\ell}^{-1} \| \mathcal{U}_T \mu - \mu \|_{\ell}^2,
\end{aligned}
$$

where we have used the fact that $h_T^{-1}h_F \lesssim 1$ 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 must show that HHO reproduces the $w \in \bar{V}_\ell^c$. We take $w_{\partial T} = \gamma_\ell w$.

$$ \int_T \mathbf{q}_T(w, w_{\partial T}) \cdot \mathbf{p}_T = \int_T w \nabla \cdot \mathbf{p}_T - \sum_{F \in \mathcal{F}_T} \int_F w_{\partial T} \mathbf{p}_T \cdot \mathbf{v}, $$

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

$$ \int_T w \nabla \cdot \mathbf{p}_T = - \int_T \nabla w \cdot \mathbf{p}_T + \int_{\partial T} w_{\partial T} \mathbf{p}_T \cdot \mathbf{v} $$

Plugging that yields

$$ \int_T \mathbf{q}_T(w, w_{\partial T}) \cdot \mathbf{p}_T = - \int_T \nabla w \cdot \mathbf{p}_T \quad \forall \mathbf{p}_T \in \mathbf{W}_T, $$

which shows that

$$ \mathbf{q}_T(w, w_{\partial T}) = -\nabla w. \qquad (6.20) $$

Given the definition (3.11) of $\mathcal{U}_T$, we need to show that

$$ \underline{a}_T((w, w_{\partial T}), (v_T, 0)) = 0 \quad \forall v_T \in V_T. $$

We have

$$
\begin{aligned}
\underline{a}_T(w, w_{\partial T}, (v_T, 0)) &= \int_T \mathbf{q}_T(w, w_{\partial T}) \cdot \mathbf{q}_T(v_T, 0) + \underbrace{s_T((w, w_{\partial T}), (v_T, 0))}_{0} \\
&= -\int_T \nabla w \cdot \mathbf{q}_T(v_T, 0),
\end{aligned}
$$

where the stabilization term vanishes due to the polynomial consistency of $s_T$ [DPD20, Assumption 2.4 (S3)], and we have used (6.20) for the last equality. For the remaining term, (3.8) gives

$$ \int_T \mathbf{q}_T(v_T, 0) \cdot \mathbf{p}_T - \int_T v_T \nabla \cdot \mathbf{p}_T = 0 \quad \forall \mathbf{p}_T \in \mathbf{W}_T, $$

which we specialize to $\mathbf{p}_T = -\nabla w$ to infer

$$ - \int_T \mathbf{q}_T(v_T, 0) \cdot \nabla w + \int_T v_T \Delta w = 0. $$

As $w \in \bar{V}_\ell^c$, $\Delta w = 0$, which concludes the proof of $\mathcal{U}\gamma_\ell w = w$.
Remark 6.2 (Generalization to $p=0$). Equation (HM4) is one reason why we assume that $p \ge 1$ throughout this manuscript. However, this assumption can be relaxed by inserting the $L^2$-orthogonal projector at several locations in the above argument. This observation motivates that our multigrid framework can be extended to the case $p=0$ if one can find ways to make the following criterion work, which might work only under additional constraints, cf. Remark 8.1.

6.6. Verification of (HM5). This is an extension of Theorem 2.27 in [CEP21] by Sobolev interpolation. Following Lemma 2.9, the left-hand side of (HM5) is bounded by the consistency error estimates. The (HM5) is derived using Theorem 2.10 in [CEP21].

6.7. Verification of (HM6). Combining some preliminary results, we have that

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

Thus, we only need to prove that

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

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

$$ \begin{aligned} \|\mu\|_{1,\ell}^2 &\lesssim \|\nabla_{\ell} U_{\ell}\mu + Q_{\ell}\mu\|_{\Omega}^2 + \|Q_{\ell}\mu\|_{\Omega}^2 + |\mu|_{1,\ell} \\ &\lesssim h_{\ell}^{-2} \|U_{\ell}\mu - \mu\|_{\ell}^2 + \|Q_{\ell}\mu\|_{\Omega}^2 + |\mu|_{1,\ell}. \end{aligned} \qquad (6.22) $$

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

$$ |\mu|_{1,\ell} \approx h_{\ell}^{-2} \sum_{T \in \mathcal{T}_{\ell}} \sum_{F \in \mathcal{F}_{T}} h_{\ell} \|U_{\ell}\mu - \mu\|_{F}^{2} \stackrel{(6.4)}{\approx} h_{\ell}^{-2} \|U_{\ell}\mu - \mu\|_{\ell}^{2}. \qquad (6.23) $$

Plugging (6.23) into (6.22) gives

$$ \|\mu\|_{1,\ell}^2 \lesssim h_{\ell}^{-2} \|U_{\ell}\mu - \mu\|_{\ell}^2 + \|Q_{\ell}\mu\|_{\Omega}^2 \stackrel{(\text{HM1})}{\lesssim} \|\mu\|_{a_{\ell}}^2. \qquad (6.24) $$

Finally, (6.21) follows from (HM2).

6.8. Verification of (HM7). This equality corresponds to [DPD20, Assumption 2.4 (S3)].

6.9. Verification of (HM8). This is relation (6.12).

## 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 $\hat{\mathcal{F}}_{\ell} \cup \hat{\mathcal{F}}_{\ell}$, where

$$ \hat{\mathcal{F}}_{\ell} := \{F \in \mathcal{F}_{\ell} \mid \exists F_{\ell-1} \in \mathcal{F}_{\ell-1} \text{ s.t. } F \subset F_{\ell-1}\}, \quad \hat{\mathcal{F}}_{\ell} := \mathcal{F}_{\ell} \setminus \hat{\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} := \begin{cases} \mu|_{F} & \text{if } F \in \hat{\mathcal{F}}_{\ell}, \\ (U_{\ell-1}\mu)|_{F} & \text{otherwise.} \end{cases} \qquad (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. $\hat{\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 $\hat{\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} := \{\{\mathcal{U}_{\ell-1}\mu\}\}_{F}. \qquad (7.2) $$

Remark that this formula also holds for $F \in \hat{\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} := \pi_{\ell}^{p}\{\{\theta_{\ell-1}^{p+1}(\mathcal{U}_{\ell-1}\mu, \mu)\}\}_{F}, \qquad (7.3) $$

where $\pi_{\ell}^{p}: 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.

While the first injection operator is local with respect to each coarse element, the later injection operators use information from their neighboring elements, increasing their domain of dependence. Thus, the first injection operator can be implemented more efficiently than the latter ones. On the contrary, the evaluation of the HHO reconstruction produces negligible cost (solving a local system of equations per cell) if the degree of parallelization is high.

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}^{i}$ is their third operator.

Let us consider $I_{\ell}^{2}$: condition (IA2) follows directly from (HM4) recalling that we assume here $p \ge 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 triangle inequality $\|\{\mathcal{U}_{\ell-1}\mu\}_{F}\|_{\ell} \le \|\mathcal{U}_{\ell-1}\mu\|_{\ell}$.
Let us consider (IA2) for $I_{\ell}^3$. On one hand, (HM4) says that

$$-\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 \bar{V}^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.9)) 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

$$\begin{align*}
\sum_{T \in \mathcal{T}_{\ell-1}} \sum_{\bar{T} \subset T} \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} \\
&\le \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} \\
&\lesssim \sum_{T \in \mathcal{T}_{\ell-1}} h_T \|\nabla \theta_{\ell-1}^{p+1}(\mathcal{U}_{\ell-1}\mu|_{T}, \mu|_{\partial T}) - \nabla \mathcal{U}_{\ell-1}\mu|_{T}\|_{T} + \|\mathcal{U}_{\ell-1}\mu|_{T}\|_{\Omega} \\
&\lesssim h_T \|\mu\|_{a_{\ell-1}} + \|\mathcal{U}_{\ell-1}\mu|_{T}\|_{\Omega} \lesssim \|\mu\|_{\ell-1}.
\end{align*}$$

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). $\square$

## 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 an L-shaped 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$. 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$ is not considered in the numerical tests as it 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. This limitation does not appear to be 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. Specifically, for V(1,1), the pre-smoothing iteration is performed in the forward order and the post-smoothing
<table><caption>TABLE 1. For the square domain, number of face unknowns at each level and for each value of p.</caption><thead><tr><th>Levels</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th></tr></thead><tbody><tr><td>p = 1</td><td>6016</td><td>24,320</td><td>97,792</td><td>392,192</td><td>1,570,816</td></tr><tr><td>p = 2</td><td>9024</td><td>36,480</td><td>146,688</td><td>588,288</td><td>2,356,224</td></tr><tr><td>p = 3</td><td>12,032</td><td>48,640</td><td>195,584</td><td>784,384</td><td>3,141,632</td></tr></tbody></table>

<table><caption>TABLE 2. In the square domain, number of V(2,2) iterations with the injection operator $I_{\ell}^{1}$.</caption><thead><tr><th>Levels</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th></tr></thead><tbody><tr><td>p = 1</td><td>13</td><td>14</td><td>14</td><td>15</td><td>15</td></tr><tr><td>p = 2</td><td>36</td><td>38</td><td>38</td><td>39</td><td>39</td></tr><tr><td>p = 3</td><td>&infin;</td><td>&infin;</td><td>&infin;</td><td>&infin;</td><td>&infin;</td></tr></tbody></table>

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 $||r||_2/||b||_2$, where r denotes the residual of the algebraic system, 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 $||r||_2/||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: (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 notice that the V(1,1) cycle diverges regardless of the value of p. Increasing the number of smoothing steps restores 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.

**Remark 8.2.** (Preconditioning) For the cases where the multigrid solver diverges, good performance may nonetheless be achieved if it is used as a preconditioner. For
<table><caption>TABLE 3. In the square domain, number of iterations with the injection operator $I_{\ell}^2$.</caption><thead><tr><th rowspan="2">Levels</th><th colspan="5">V(1,1)</th><th colspan="5">V(2,2)</th></tr><tr><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th></tr></thead><tbody><tr><th>p = 1</th><td>24</td><td>25</td><td>25</td><td>26</td><td>26</td><td>13</td><td>13</td><td>14</td><td>14</td><td>14</td></tr><tr><th>p = 2</th><td>20</td><td>22</td><td>25</td><td>26</td><td>27</td><td>10</td><td>10</td><td>11</td><td>11</td><td>11</td></tr><tr><th>p = 3</th><td>&infin;</td><td>&infin;</td><td>&infin;</td><td>&infin;</td><td>&infin;</td><td>13</td><td>13</td><td>13</td><td>14</td><td>14</td></tr></tbody></table>

<table><caption>TABLE 4. In the square domain, number of iterations with the injection operator $I_{\ell}^3$.</caption><thead><tr><th rowspan="2">Levels</th><th colspan="5">V(1,1)</th><th colspan="5">V(2,2)</th></tr><tr><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th></tr></thead><tbody><tr><th>p = 1</th><td>18</td><td>18</td><td>19</td><td>19</td><td>20</td><td>10</td><td>10</td><td>11</td><td>11</td><td>11</td></tr><tr><th>p = 2</th><td>17</td><td>17</td><td>17</td><td>17</td><td>18</td><td>9</td><td>10</td><td>10</td><td>10</td><td>10</td></tr><tr><th>p = 3</th><td>20</td><td>21</td><td>21</td><td>21</td><td>21</td><td>11</td><td>11</td><td>11</td><td>11</td><td>11</td></tr></tbody></table>

<table><caption>TABLE 5. In the square domain, number conjugate gradient iterations, preconditioned with one V(1,1) multigrid cycle, using the injection operator $I_{\ell}^1$.</caption><thead><tr><th>Levels</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th></tr></thead><tbody><tr><th>p = 1</th><td>18</td><td>19</td><td>20</td><td>21</td><td>21</td></tr><tr><th>p = 2</th><td>20</td><td>22</td><td>24</td><td>25</td><td>26</td></tr><tr><th>p = 3</th><td>21</td><td>23</td><td>25</td><td>26</td><td>27</td></tr></tbody></table>

instance, consider the multigrid solver with V(1,1) cycle and the injection operator $I_{\ell}^1$. When used as a solver, the method diverges for all values of $p$. On the other hand, used as a preconditioner for the conjugate gradient, the method only exhibits a mild dependency on the number of unknowns, as shown in Table 5.

8.2.2. *Low regularity: the L-shaped domain.* The computational domain is now the L-shaped 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: (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 would require the variable cycle, 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 6. The results, presented in Tables 7, 8, and 9 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 7). With $I_{\ell}^2$, V(1,1) provides a uniformly convergent solver up to $p = 2$, while V(2,2) converges uniformly for all
<table><caption>TABLE 6. For the L-shaped domain, number of face unknowns at each level and for each value of $p$.</caption><thead><tr><th>Levels</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th></tr></thead><tbody><tr><td>p = 1</td><td>4480</td><td>18,176</td><td>73,216</td><td>293,888</td><td>1,177,600</td></tr><tr><td>p = 2</td><td>6720</td><td>27,264</td><td>109,824</td><td>440,832</td><td>1,766,400</td></tr><tr><td>p = 3</td><td>8960</td><td>36,352</td><td>146,432</td><td>587,776</td><td>2,355,200</td></tr></tbody></table>

<table><caption>TABLE 7. In the L-shaped domain, number of V(2,2) iterations with the injection operator $I_{\ell}^{1}$.</caption><thead><tr><th>Levels</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th></tr></thead><tbody><tr><td>p = 1</td><td>15</td><td>12</td><td>14</td><td>15</td><td>15</td></tr><tr><td>p = 2</td><td>&infin;</td><td>&infin;</td><td>&infin;</td><td>&infin;</td><td>&infin;</td></tr><tr><td>p = 3</td><td>&infin;</td><td>&infin;</td><td>&infin;</td><td>&infin;</td><td>&infin;</td></tr></tbody></table>

<table><caption>TABLE 8. In the L-shaped domain, number of iterations with the injection operator $I_{\ell}^{2}$.</caption><thead><tr><th rowspan="2">Levels</th><th colspan="5">V(1,1)</th><th colspan="5">V(2,2)</th></tr><tr><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th></tr></thead><tbody><tr><td>p = 1</td><td>20</td><td>20</td><td>20</td><td>20</td><td>20</td><td>11</td><td>11</td><td>11</td><td>11</td><td>11</td></tr><tr><td>p = 2</td><td>17</td><td>17</td><td>17</td><td>17</td><td>17</td><td>9</td><td>9</td><td>9</td><td>9</td><td>9</td></tr><tr><td>p = 3</td><td>&infin;</td><td>&infin;</td><td>&infin;</td><td>&infin;</td><td>&infin;</td><td>11</td><td>11</td><td>11</td><td>11</td><td>11</td></tr></tbody></table>

$p \in \{1, 2, 3\}$ (cf. Table 8). Finally, using $I_{\ell}^{3}$, Table 9 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).

### 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 10. The manufactured solution is $u: (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
<table><caption>TABLE 9. In the L-shaped domain, number of iterations with the injection operator $I_{\ell}^{3}$.</caption><thead><tr><th rowspan="2">Levels</th><th colspan="5">V(1,1)</th><th colspan="5">V(2,2)</th></tr><tr><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th></tr></thead><tbody><tr><td>p = 1</td><td>16</td><td>16</td><td>16</td><td>16</td><td>16</td><td>9</td><td>9</td><td>10</td><td>11</td><td>11</td></tr><tr><td>p = 2</td><td>17</td><td>17</td><td>17</td><td>17</td><td>17</td><td>8</td><td>9</td><td>9</td><td>9</td><td>9</td></tr><tr><td>p = 3</td><td>20</td><td>20</td><td>20</td><td>20</td><td>20</td><td>10</td><td>10</td><td>10</td><td>10</td><td>10</td></tr></tbody></table>

<table><caption>TABLE 10. For the cubic domain, number of face unknowns at each level and for each value of $p$.</caption><thead><tr><th>Levels</th><th>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>p = 1</td><td>17,280</td><td>142,848</td><td>1,161,216</td></tr><tr><td>p = 2</td><td>34,560</td><td>285,696</td><td>2,322,432</td></tr><tr><td>p = 3</td><td>57,600</td><td>476,160</td><td>3,870,720</td></tr></tbody></table>

<table><caption>TABLE 11. In the cubic domain, number of V(2,2) iterations with the injection operator $I_{\ell}^{3}$.</caption><thead><tr><th>Levels</th><th>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>p = 1</td><td>18</td><td>23</td><td>25</td></tr><tr><td>p = 2</td><td>23</td><td>24</td><td>22</td></tr><tr><td>p = 3</td><td>22</td><td>23</td><td>23</td></tr></tbody></table>

presented in Table 11, 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.

### 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.

There are three important generalizations to our approach that we would like to address in future work: First, we would like to extend our theory to polyhedral meshes. Second, we would like to extend our multigrid framework to work as a preconditioner. Third, we want to extend this work to piece-wise constant approximations as detailed in Remark 6.2.
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

[ABBV23] P. F. Antonietti, S. Berrone, M. Busetto, and M. Verani. Agglomeration-based geometric multigrid schemes for the virtual element method. *SIAM Journal on Numerical Analysis*, 61(1):223–249, 2023. 10.1137/21M1466864.

[AP19] P. F. Antonietti and G. Pennesi. *Journal of Scientific Computing*, 78(1):625–652, 2019. 10.1007/s10915-018-0783-x.

[AVZ15] P. F. Antonietti, M. Verani, and L. Zikatanov. A two-level method for mimetic finite difference discretizations of elliptic problems. *Computers & Mathematics with Applications*, 70(11):2674–2687, 2015. Numerical Methods for Scientific Computations and Advanced Applications. 10.1016/j.camwa.2015.06.010.

[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.

[Cal19] J. G. Calvo. An overlapping schwarz method for virtual element discretizations in two dimensions. *Computers & Mathematics with Applications*, 77(4):1163–1177, 2019. 10.1016/j.camwa.2018.10.043.
[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.

[PP22] Y. Pan and P.-O. Persson. Agglomeration-based geometric multigrid solvers for compact discontinuous galerkin discretizations on unstructured meshes. *Journal of Computational Physics*, 449:110775, 2022. 10.1016/j.jcp.2021.110775.

[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.
IMAG, UNIV. MONTPELLIER, CNRS, MONTPELLIER, FRANCE

*Email address:* daniele.di-pietro@umontpellier.fr

INRIA, 48 RUE BARRAULT, 75647 PARIS, FRANCE, 2) CERMICS, ECOLE DES PONTS, 77455 MARNE-LA- VALLÉE, FRANCE

*Email address:* zhaonan.dong@inria.fr

INTERDISCIPLINARY CENTER FOR SCIENTIFIC COMPUTING (IWR), HEIDELBERG UNIVERSITY, MATHEMATIKON, IM NEUENHEIMER FELD 205, 69120 HEIDELBERG, GERMANY

*Email address:* kanschat@uni-heidelberg.de

CMAP, CNRS, ÉCOLE POLYTECHNIQUE, INSTITUT POLYTECHNIQUE DE PARIS, 91120 PALAISEAU, FRANCE

*Email address:* pierre.matalon@polytechnique.edu

DEPARTMENT OF MATHEMATICS, SAARLAND UNIVERSITY, DE-66123 SAARBRÜCKEN, GERMANY

*Email address:* andreas.rupp@uni-saarland.de