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

HAL
open science

# hp-error analysis of mixed-order hybrid high-order methods for elliptic problems on simplicial meshes

Zhaonan Dong, Alexandre Ern

► To cite this version:

Zhaonan Dong, Alexandre Ern. *hp*-error analysis of mixed-order hybrid high-order methods for elliptic problems on simplicial meshes. 2025. hal-04720237v2

HAL Id: hal-04720237

https://inria.hal.science/hal-04720237v2

Preprint submitted on 24 Jul 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
hp-error analysis of mixed-order hybrid high-order methods for
elliptic problems on simplicial meshes

Zhaonan Dong*

Alexandre Ern†

July 24, 2025

Abstract

We present both hp-a priori and hp-a posteriori error analysis of a mixed-order hybrid high-order (HHO) method to approximate second-order elliptic problems on simplicial meshes. Our main result on the hp-a priori error analysis is a 1/2-order p-suboptimal error estimate. This result is, to our knowledge, the first of this kind for hybrid nonconforming methods and matches the state-of-the-art for other nonconforming methods (as discontinuous Galerkin methods) with general (mixed Dirichlet/Neumann) boundary conditions. Our second main result is a residual-based hp-a posteriori upper error bound, comprising residual, normal flux jump, tangential jump, and stabilization estimators (plus data oscillation terms). The first three terms are p-optimal and only the latter is 1/2-order p-suboptimal. This result is, to our knowledge, the first hp-a posteriori error estimate for HHO methods. A novel approach based on the partition-of-unity provided by hat basis functions and on local Helmholtz decompositions on vertex stars is devised to estimate the nonconforming error. Finally, we establish local lower error bounds. Remarkably, the normal flux jump estimator is only 1/2-order p-suboptimal, as it can be bounded by the stabilization owing to the local conservation property of HHO methods. Numerical examples illustrate the theory.

# 1 Introduction

Second-order elliptic PDEs are widely used in the modeling of diffusion phenomena. In the present work, we consider the following model problem:

$$ egin{cases} -\nabla \cdot (A \nabla u) = f & \text{in } \Omega, \\ u = g_D & \text{on } \Gamma_D, \\ (A \nabla u) \cdot n_{\Omega} = g_N & \text{on } \Gamma_N, \end{cases} \quad (1) $$

where the domain $\Omega$ is a Lipschitz polytopal domain (open, connected, bounded set) in $\mathbb{R}^d$, $d \in \{2, 3\}$, with boundary $\partial\Omega$ and unit outward normal $n_{\Omega}$. The boundary $\partial\Omega$ is split into two disjoint parts $\Gamma_D$ and $\Gamma_N$ with $|\Gamma_D| > 0$. In addition, the load $f \in L^2(\Omega)$, $g_D$ is a restriction to $\Gamma_D$ of a function in $H^{1/2}(\partial\Omega)$, $g_N \in L^2(\Gamma_N)$ and $A$ is a piecewise scalar-valued diffusion coefficient such that $0 < A_{\Omega}^b \le A(x) \le A_{\Omega}^s$ for a.e. $x \in \Omega$.

The hybrid high-order (HHO) method was introduced in [20] for linear diffusion and in [19] for locking-free linear elasticity. As shown in [16], the HHO method is closely related to hybridizable discontinuous Galerkin (HDG) and weak Galerkin (WG) methods. These links have been leveraged, e.g., in [22, 27] to devise a unified convergence analysis for the biharmonic problem and the acoustic wave equation. We also refer the reader to [41] and [15] for links to the nonconforming virtual element method (ncVEM), and to [13] for links to multiscale hybrid-mixed (MHM) methods. HHO methods are formulated in terms of broken cell and face polynomial spaces. The equal-order HHO method corresponds to cell and face unknowns having the same degree $k \ge 0$. Instead, the mixed-order HHO method corresponds to cell unknowns having degree $(k+1)$ and face unknowns having degree $k \ge 0$. Considering cell unknowns of

<footnote>*Inria, 48 rue Barrault, 75647 Paris, France, and CERMICS, ENPC, Institut Polytechnique de Paris, 6 & 8 avenue B. Pascal, 77455 Marne-la-Vallée, France zhaonan.dong@inria.fr.</footnote>

<footnote>†CERMICS, ENPC, Institut Polytechnique de Paris, 6 & 8 avenue B. Pascal, 77455 Marne-la-Vallée, France, and Inria, 48 Rue Barrault, 75647 Paris, France alexandre.ern@enpc.fr.</footnote>
degree ($k-1$) with $k \ge 1$ is also possible. One salient advantage of the mixed-order setting with cell unknowns of degree ($k+1$) is that $h$-optimal convergence can be achieved by using the simple Lehrenfeld–Schöberl (LS) HDG stabilization [40], in contrast with the more sophisticated HHO stabilization needed in the other settings.

The goal of the present work is to derive hp-a priori and hp-a posteriori error estimates for the mixed-order HHO method on simplicial meshes. Owing to the links highlighted above, the present results extend to HDG and WG methods in the same setting. Although HHO (and HDG, WG) methods can deal with polytopal meshes, we focus here on simplicial meshes because some of the hp-analysis tools we are going to invoke are only available on such meshes (and on tensor-product meshes as well). Our main result concerning the hp-a priori error analysis is Theorem 4.4, where we derive a $\frac{1}{2}$-order p-suboptimal and h-optimal error estimate. The only hp-a priori HHO error estimate we are aware of is derived in [1] for the equal-order HHO method and leads to a 1-order p-suboptimal error bound. Here, we use the same hp-scaling of the stabilization bilinear form, but the simpler form of the LS stabilization allows us to prove an error estimate with a tighter scaling in the polynomial degree. Notice in passing that $\frac{1}{2}$-order p-suboptimality corresponds to the state of the art for classical discontinuous Galerkin (dG) methods [35, 31, 8] under general Dirichlet/Neumann boundary conditions, and was also obtained in [23] for a hybrid dG method applied to Stokes flow. In the case of homogeneous Dirichlet conditions over the whole boundary, it is shown in [48] for dG methods applied to the Poisson problem and in [39] for HDG methods applied to the Stokes problem that an a priori error estimate with full p-optimality can be achieved. However, counterexamples in [31] confirm that the $\frac{1}{2}$-order p-suboptimality is sharp in the case of inhomogeneous Dirichlet boundary conditions.

The second main contribution of the present work is to derive a residual-based hp-a posteriori error estimate for the mixed-order HHO method in dimensions $d \in \{2, 3\}$. Our main result concerning the (global) upper error bound is Theorem 5.7, where all the terms in the upper bound are hp-optimal except one term which is $\frac{1}{2}$-order p-suboptimal. To the best of our knowledge, this is the first such estimate for HHO methods, whereas h-a posteriori error estimates for HHO methods were derived previously in [21, 5, 10], focusing on either the equal-order HHO setting or a stabilization-free variant of the method.

The main challenge in deriving an upper error bound for nonconforming methods is to estimate the nonconforming error, which essentially measures by how much the discrete solution departs from $H^1$. A first possibility is to invoke a nodal-averaging operator mapping to $H^1$, as done, e.g., in [36, 28, 2, 29] for dG methods. However, p-optimal approximation results for nodal-averaging operators are so far available only on tensor-product [6] and triangular [34] meshes, whereas the best bound available on tetrahedral meshes is p-suboptimal by one order [23, Lemma 7.6]. An alternative to using a nodal-averaging operator is to invoke a (global) Helmholtz decomposition on the nonconforming error, as in [18, 9, 11] for Crouzeix–Raviart finite elements and in [4, 7] for dG methods. Here, we adopt this technique, but we introduce a novel idea in that we additionally use the partition of unity provided by the hat basis functions to invoke a *local* Helmholtz decomposition on each vertex star (the subdomain covered by the mesh cells sharing the vertex). The benefit is that each vertex star is simply connected, and recent results on the stability constant in the local Helmholtz decomposition are available [32]. Instead, the stability constant for a global Helmholtz decomposition in a domain with $N$ holes grows unfavorably with $N$ [5]. Finally, we emphasize that our novel idea for bounding the nonconforming error can be directly applied to the a posteriori analysis of many other nonconforming discretization methods, such as dG, HDG, and WG.

We also address the efficiency of our a posteriori error estimate by establishing local lower error bounds. Our main result is Theorem 5.9 which leads to $\frac{3}{2}$-order p-suboptimality (only the tangential jump estimator leads to such suboptimality). Numerical experiments, though, indicate only $\frac{1}{2}$-order suboptimality in $p$. Our proof of the lower error bound uses bubble function techniques inspired from [44], but we introduce a novel argument in the proof in that we invoke the local conservation property of the HHO method to improve the efficiency result on the normal flux jump from $\frac{3}{2}$-order to $\frac{1}{2}$-order p-suboptimality. Another interesting numerical observation is that the normal flux jump is not the dominant component of the a posteriori error estimate for HHO methods, in contrast to the situation classically encountered with conforming finite elements [12].

The rest of this work is organised as follows. We present the weak formulation of the model problem together with the discrete setting in Section 2. In Section 3, we introduce the HHO method, and in Section 4, we derive the hp-a priori error estimate. In Section 5, we present the residual-based hp-a posteriori error analysis, leading to a (global) upper error bound and (local) lower error bounds. Numerical experiments are presented in Section 6 to illustrate the theory. Finally, in Section 7, we
collect several (technical) proofs related to the hp-a posteriori error analysis.

# 2 Weak form and discrete setting

In this section, we introduce some basic notation, the weak formulation of the model problem, and the discrete setting to formulate and analyze the HHO discretization.

## 2.1 Basic notation and weak formulation

We use standard notation for the Lebesgue and Sobolev spaces and, in particular, for the fractional-order Sobolev spaces, we consider the Sobolev-Slobodeckij seminorm based on the double integral. For an open, bounded, Lipschitz set $S \in \mathbb{R}^d$, $d \in \{1, 2, 3\}$, we denote by $(v, w)_S$ the $L^2(S)$-inner product, and we employ the same notation when $v$ and $w$ are vector- or matrix-valued fields. We denote by $\nabla w$ the (weak) gradient of $w$. We use boldface notation to denote vectors in $\mathbb{R}^d$, as well as $\mathbb{R}^d$-valued fields and functional spaces composed of such fields.

Setting $H_{g,D}^1(\Omega) := \{v \in H^1(\Omega) \mid v|_{\Gamma_D} = g_D\}$, the weak formulation of (1) is as follows: Find $u \in H_{g,D}^1(\Omega)$, such that

$$ (A\nabla u, \nabla v)_{\Omega} = (f, v)_{\Omega} + (g_N, v)_{\Gamma_N}, \quad (2) $$

for all $v \in H_{0,D}^1(\Omega)$. The well-posedness of (2) follows from the Lax-Milgram Lemma, see, e.g., [26, Proposition 31.21].

## 2.2 Mesh

Let $\mathcal{T}_h$ be a simplicial mesh that covers the domain $\Omega$ exactly and is compatible with the boundary partition as well as the domain partition on which $A$ is piecewise constant. A generic mesh cell is denoted by $K \in \mathcal{T}_h$, its diameter by $h_K$, and its unit outward normal by $n_K$. We set $A_K := A|_K$. We let $\text{st}(K)$ denote the collection of cells $\hat{K} \in \mathcal{T}_h$ sharing at least one vertex with the cell $K$ ($\text{st}(K)$ is often called cell star). Similarly, $\text{es}(K)$ denotes the collection of cells $\hat{K} \in \mathcal{T}_h$ sharing at least one vertex with $\text{st}(K)$ ($\text{es}(K)$ is often called extended cell star). For all $k \ge 0$, $\mathbb{P}^k(K)$ denotes the space of $d$-variate polynomials on $K$ of degree at most $k$, and $\Pi_K^k$ denotes the $L^2$-orthogonal projection onto $\mathbb{P}^k(K)$. Moreover, $\mathbb{P}^k(\mathcal{T}_h) := \{v_h \in L^2(\Omega) \mid v_h|_K \in \mathbb{P}^k(K)\}$ denotes the broken polynomial space of order $k$ on the mesh $\mathcal{T}_h$ (classically considered in dG methods).

The mesh faces are collected in the set $\mathcal{F}_h$, which is split as $\mathcal{F}_h = \mathcal{F}_h^i \cup \mathcal{F}_h^b$, where $\mathcal{F}_h^i$ is the collection of interfaces (shared by two distinct mesh cells) and $\mathcal{F}_h^b$ the collection of boundary faces. Moreover, we split $\mathcal{F}_h^b$ into the Dirichlet subset, $\mathcal{F}_h^D$, and the Neumann subset, $\mathcal{F}_h^N$. Let $n_F$ denote the unit normal vector orienting the mesh face $F \in \mathcal{F}_h$. For all $F \in \mathcal{F}_h^i$, the direction of $n_F$ is arbitrary, but fixed, whereas we set $n_F := n_\Omega|_F$ for all $F \in \mathcal{F}_h^b$. For every mesh cell $K \in \mathcal{T}_h$, the partition of its boundary $\partial K$ is defined as $\partial K = \partial K^i \cup \partial K^D \cup \partial K^N$ with obvious notation, and the mesh faces composing $\partial K$ are collected in the set $\mathcal{F}_{\partial K}$, which is partitioned as $\mathcal{F}_{\partial K} = \mathcal{F}_{\partial K^i} \cup \mathcal{F}_{\partial K^D} \cup \mathcal{F}_{\partial K^N}$ with obvious notation. For all $F \in \mathcal{F}_h$, $\Pi_F^k$ denotes the $L^2$-orthogonal projection onto $\mathbb{P}^k(F)$.

The set of mesh vertices is denoted by $\mathcal{V}_h$ and is decomposed into the subset of interior vertices, $\mathcal{V}_h^i$, and the subset of boundary vertices, $\mathcal{V}_h^b$. For all $a \in \mathcal{V}_h$, $\mathcal{T}_a$ denotes the collection of mesh cells which share $a$ and $\omega_a$ the corresponding open subdomain (often called vertex star). In addition, we define $\mathcal{F}_a$ as the collection of faces in $\mathcal{F}_h$ which share $a$.

For all $s > \frac{1}{2}$, we define the broken Sobolev spaces $H^s(\mathcal{T}_h; \mathbb{R}^q) := \{w \in L^2(\Omega; \mathbb{R}^q) \mid w_K := w|_K \in H^s(K; \mathbb{R}^q), \forall K \in \mathcal{T}_h\}$, $q \in \{1, d\}$. We define the jump $[[w]]_F$ of any function $w \in H^s(\mathcal{T}_h; \mathbb{R}^q)$ across any mesh interface $F = \partial K_1 \cap \partial K_2 \in \mathcal{F}_h^i$ as $[[w]]_F := w_{K_1|_F} - w_{K_2|_F}$, where $n_F$ points from $K_1$ to $K_2$. For any boundary face $F = \partial K \cap \partial \Omega \in \mathcal{F}_h^b$, we set $[[w]]_F := w_K|_F$. For all $K \in \mathcal{T}_h$, we define $[[w]]_{\partial K}|_F := [[w]]_F$ for all $F \in \mathcal{F}_{\partial K}$. Finally, we define the broken gradient $\nabla_{\mathcal{T}_h}$ as the gradient operator acting cellwise on $H^1(\mathcal{T}_h; \mathbb{R}^q)$.

## 2.3 hp-analysis tools

Let us briefly review the main hp-analysis tools used in this work. We use the symbol $C$ (sometimes with a subscript) to denote any positive generic constant whose value can change at each occurrence as
long as it is independent of the mesh size h and the polynomial degree k. The value of C can depend on the mesh shape-regularity and the space dimension d.

**Lemma 2.1** (Discrete trace inequality). The following holds for all $v \in \mathbb{P}^k(K)$, all $K \in \mathcal{T}_h$, and all $k \ge 0$,

$$ \|v\|_{\partial K} \le C \frac{(k+1)}{h_K^{\frac{1}{2}}} \|v\|_K, \qquad (3) $$

*Proof.* A proof can be found in [49, Theorem 5] (with explicit constant in terms of d). □

**Lemma 2.2** (Local $L^2$-orthogonal projection). The following holds for $v \in H^1(K)$, all $K \in \mathcal{T}_h$, and all $k \ge 0$,

$$ \|v - \Pi_K^k(v)\|_{\partial K} \le C \left( \frac{h_K}{k+1} \right)^{\frac{1}{2}} |v|_{H^1(K)}. \qquad (4) $$

*Proof.* A proof can be found in [14, 46]. □

**Lemma 2.3** (Local Babuška–Suri operator). There exists a positive constant $C_{BS}$ such that, for all $k \ge 1$ and all $K \in \mathcal{T}_h$, there exists an operator $\mathcal{I}_{BS,K}^k : L^2(K) \to \mathbb{P}^k(K)$, called Babuška–Suri approximation operator, such that, for all $r \in \{0, \dots, k\}$, all $m \in \{0, \dots, r\}$, and all $v \in H^r(K)$,

$$ |v - \mathcal{I}_{BS,K}^k(v)|_{H^m(K)} \le C_{BS} \left( \frac{h_K}{k} \right)^{r-m} \|v\|_{H^r(K)}. \qquad (5) $$

*Proof.* A proof can be found in [3]. □

**Lemma 2.4** (Global hp-Karkulik–Melenk operator). There exists a constant $C_{KM}$ such that, for all $k \ge 1$, there exists an operator $I_{KM}^k : H_{0,D}^1(\Omega) \to \mathbb{P}^k(\mathcal{T}_h) \cap H_{0,D}^1(\Omega)$, called Karkulik–Melenk interpolation operator, such that, for all $v \in H_{0,D}^1(\Omega)$ and all $K \in \mathcal{T}_h$,

$$ \left(\frac{k}{h_K}\right)^2 \|v - I_{KM}^k(v)\|_K^2 + \left(\frac{k}{h_K}\right) \|v - I_{KM}^k(v)\|_{\partial K}^2 + \|\nabla I_{KM}^k(v)\|_K^2 \le C_{KM} \|v\|_{H^1(\text{st}(K))}^2. \qquad (6) $$

*Proof.* A proof can be found in [43, 37]. □

**Corollary 2.5** (Modified hp-Karkulik–Melenk operator). There exists a constant $C_{mKM}$ such that, for all $k \ge 1$, there exists an operator $I_{mKM}^k : H_{0,D}^1(\Omega) \to \mathbb{P}^k(\mathcal{T}_h) \cap H_{0,D}^1(\Omega)$, called modified Karkulik–Melenk interpolation operator, such that, for all $v \in H_{0,D}^1(\Omega)$ and all $K \in \mathcal{T}_h$,

$$ \left(\frac{k}{h_K}\right)^2 \|v - I_{mKM}^k(v)\|_K^2 + \left(\frac{k}{h_K}\right) \|v - I_{mKM}^k(v)\|_{\partial K}^2 + \|\nabla I_{mKM}^k(v)\|_K^2 \le C_{mKM} \|\nabla v\|_{\text{es}(K)}^2. \qquad (7) $$

*Proof.* The idea is to set, for all $v \in H_{0,D}^1(\Omega)$,

$$ I_{mKM}^k(v) := I_{av,D}^1(v) + I_{KM}^k(v - I_{av,D}^1(v)), $$

where the (first-order) nodal-averaging operator $I_{av,D}^1 : H_{0,D}^1(\Omega) \to \mathbb{P}^1(\mathcal{T}_h) \cap H_{0,D}^1(\Omega)$ is devised, e.g., in [24] when $\Gamma_D = \partial\Omega$ and in [42] when $\Gamma_D$ is a proper subset of $\partial\Omega$. This operator satisfies, for all $v \in H_{0,D}^1(\Omega)$,

$$ \left(\frac{1}{h_K}\right)^2 \|v - I_{av,D}^1(v)\|_K^2 + \left(\frac{1}{h_K}\right) \|v - I_{av,D}^1(v)\|_{\partial K}^2 + \|\nabla I_{av,D}^1(v)\|_K^2 \le C \|\nabla v\|_{\text{st}(K)}^2. \qquad (8) $$

We can now prove (7). Using the approximation results (6) and (8) and the mesh shape-regularity, we infer that

$$ \left(\frac{k}{h_K}\right)^2 \|v - I_{mKM}^k(v)\|_K^2 = \left(\frac{k}{h_K}\right)^2 \left\| (v - I_{av,D}^1(v)) - I_{KM}^k(v - I_{av,D}^1(v)) \right\|_K^2 \\ \le C \|v - I_{av,D}^1(v)\|_{H^1(\text{st}(K))}^2 \le C \sum_{\widehat{K} \in \text{st}(K)} \|\nabla v\|_{\text{st}(\widehat{K})}^2 \le C \|\nabla v\|_{\text{es}(K)}^2. $$

This proves the bound on the first term on the left-hand side of (7), and the other two terms can be bounded in a similar way. □
**Remark 2.6 (Applications).** We use the Babuška–Suri approximation operator to establish our hp-a priori error estimate. Instead, we use the modified Karkulik–Melenk interpolation operator in the hp-a posteriori error analysis to establish the upper error bound. The advantage of the modified Karkulik–Melenk interpolation operator with respect to the original one is to invoke only the H¹-seminorm on the right-hand side of (7) (compare with (6)).

# 3 HHO method

Let $k \ge 0$ be the polynomial degree. We focus on the mixed-order HHO method where, for all $K \in \mathcal{T}_h$, the local HHO space is

$$ \hat{V}_K^k := \mathbb{P}^{k+1}(K) \times \mathbb{P}^k(\mathcal{F}_{\partial K}), \quad \mathbb{P}^k(\mathcal{F}_{\partial K}) := \bigtimes_{F \in \mathcal{F}_h} \mathbb{P}^k(F). \qquad (9) $$

A generic element in $\hat{V}_K^k$ is denoted by $\hat{v}_K := (v_K, v_{\partial K})$ with $v_K \in \mathbb{P}^{k+1}(K)$ and $v_{\partial K} \in \mathbb{P}^k(\mathcal{F}_{\partial K})$. The first component of the pair $\hat{v}_K$ aims at representing the solution inside the mesh cell and the second its trace at the cell boundary.

## 3.1 Reconstruction and stabilization

The HHO method is formulated locally by means of a reconstruction and a stabilization operator. The local reconstruction operator $R_K^{k+1} : \hat{V}_K^k \to \mathbb{P}^{k+1}(K)$ is such that, for all $\hat{v}_K := (v_K, v_{\partial K}) \in \hat{V}_K^k$, $R_K^{k+1}(\hat{v}_K) \in \mathbb{P}^{k+1}(K)$ is determined by solving the following well-posed problem:

$$ (\nabla R_K^{k+1}(\hat{v}_K), \nabla w)_K := (\nabla v_K, \nabla w)_K - (v_K - v_{\partial K}, \nabla w \cdot n_K)_{\partial K}, \qquad (10) $$

for all $w \in \mathbb{P}^{k+1}(K)/\mathbb{R}$ and $(R_K^{k+1}(\hat{v}_K), 1)_K = (v_K, 1)_K$. (Notice that (10) actually holds for all $w \in \mathbb{P}^{k+1}(K)$). Integration by parts gives

$$ (\nabla R_K^{k+1}(\hat{v}_K), \nabla w)_K = -(v_K, \Delta w)_K + (v_{\partial K}, \nabla w \cdot n_K)_{\partial K}. \qquad (11) $$

The local stabilization bilinear form $S_{\partial K}$ is defined such that, for all $(\hat{v}_K, \hat{w}_K) \in \hat{V}_K^k \times \hat{V}_K^k$,

$$ S_{\partial K}(\hat{v}_K, \hat{w}_K) := \frac{(k+1)^2}{h_K} (v_{\partial K} - \Pi_{\partial K}^k (v_K|_{\partial K}), w_{\partial K} - \Pi_{\partial K}^k (w_K|_{\partial K}))_{\partial K}, \qquad (12) $$

where $\Pi_{\partial K}^k$ denotes the $L^2$-orthogonal projection onto $\mathbb{P}^k(\mathcal{F}_{\partial K})$. The reconstruction and stabilization operators are combined together to build the local bilinear form $a_K$ on $\hat{V}_K^k \times \hat{V}_K^k$ such that, for all $K \in \mathcal{T}_h$,

$$ a_K(\hat{v}_K, \hat{w}_K) := A_K(\nabla R_K^{k+1}(\hat{v}_K), \nabla R_K^{k+1}(\hat{w}_K))_K + A_K S_{\partial K}(\hat{v}_K, \hat{w}_K). \qquad (13) $$

**Lemma 3.1 (Useful property).** The following holds for all $\hat{v}_K \in \hat{V}_K^k$ and all $K \in \mathcal{T}_h$:

$$ \|\nabla(R_K^{k+1}(\hat{v}_K) - v_K)\|_K^2 \le CS_{\partial K}(\hat{v}_K, \hat{v}_K). \qquad (14) $$

*Proof.* Using the definition (10) of the reconstruction operator with $w := R_K^{k+1}(\hat{v}_K) - v_K \in \mathbb{P}^{k+1}(K)$ gives

$$ \|\nabla(R_K^{k+1}(\hat{v}_K) - v_K)\|_K^2 := -(v_K - v_{\partial K}, \nabla(R_K^{k+1}(\hat{v}_K) - v_K) \cdot n_K)_{\partial K}. $$

Since $\nabla(R_K^{k+1}(\hat{v}_K) - v_K) \cdot n_K \in \mathbb{P}^k(\mathcal{F}_{\partial K})$, using the Cauchy–Schwarz inequality, the definition (12) of the stabilization, and the discrete trace inequality (3) implies that

$$ \begin{aligned} \|\nabla(R_K^{k+1}(\hat{v}_K) - v_K)\|_K^2 &= -(\Pi_{\partial K}^k(v_K|_{\partial K}) - v_{\partial K}, \nabla(R_K^{k+1}(\hat{v}_K) - v_K) \cdot n_K)_{\partial K} \\ &\le CS_{\partial K}(\hat{v}_K, \hat{v}_K)^{\frac{1}{2}} \|\nabla(R_K^{k+1}(\hat{v}_K) - v_K)\|_K. \end{aligned} $$

This concludes the proof of (14). □
## 3.2 Global discrete problem

We define the global HHO space as

$$ \hat{V}_h^k := \mathbb{P}^{k+1}(\mathcal{T}_h) \times \mathbb{P}^k(\mathcal{F}_h), \quad \mathbb{P}^{k+1}(\mathcal{T}_h) := \mathop{\times}_{K \in \mathcal{T}_h} \mathbb{P}^{k+1}(K), \quad \mathbb{P}^k(\mathcal{F}_h) := \mathop{\times}_{F \in \mathcal{F}_h} \mathbb{P}^k(F). \quad (15) $$

A generic element in $\hat{V}_h^k$ is denoted by $\hat{v}_h := (v_{\mathcal{T}_h}, v_{\mathcal{F}_h})$ with $v_{\mathcal{T}_h} := (v_K)_{K \in \mathcal{T}_h}$ and $v_{\mathcal{F}_h} := (v_F)_{F \in \mathcal{F}_h}$. For all $K \in \mathcal{T}_h$, the local components of $\hat{v}_h$ are collected in the pair $\hat{v}_K := (v_K, v_{\partial K}) \in \hat{V}_K^k$ with $v_{\partial K}|_F := v_F$ for all $F \in \mathcal{F}_{\partial K}$. Similarly, let $R_{\mathcal{T}_h}^{k+1}(\hat{v}_h) \in \mathbb{P}^{k+1}(\mathcal{T}_h)$ be such that $R_{\mathcal{T}_h}^{k+1}(\hat{v}_h)|_K = R_K^{k+1}(\hat{v}_K)$ for all $K \in \mathcal{T}_h$. To deal with the Dirichlet boundary condition, we define the (affine) subspaces

$$ \hat{V}_{h0,D}^k := \{\hat{v}_h \in \hat{V}_h^k \mid v_F = 0, \forall F \in \mathcal{F}_h^D\}, \quad (16a) $$

$$ \hat{V}_{hg,D}^k := \{\hat{v}_h \in \hat{V}_h^k \mid v_F = \Pi_F^k(g_D|_F), \forall F \in \mathcal{F}_h^D\}. \quad (16b) $$

The discrete HHO problem is as follows: Find $\hat{u}_h \in \hat{V}_{hg,D}^k$ such that

$$ a_h(\hat{u}_h, \hat{w}_h) = \ell_h(\hat{w}_h), \quad \forall \hat{w}_h \in \hat{V}_{h0,D}^k, \quad (17) $$

where the global discrete bilinear form $a_h$ and the global linear form $\ell_h$ are assembled cellwise as

$$ a_h(\hat{v}_h, \hat{w}_h) := \sum_{K \in \mathcal{T}_h} a_K(\hat{v}_K, \hat{w}_K), \quad \ell_h(\hat{w}_h) := \sum_{K \in \mathcal{T}_h} \{(f, w_K)_K + (g_N, w_{\partial K})_{\partial K^N}\}. \quad (18) $$

It is well-known that the discrete problem (17) is amenable to static condensation, i.e., the cell unknowns can be eliminated locally in every mesh cell, leading to a global problem where the only remaining unknowns are those attached to the mesh faces, i.e., those in $\mathbb{P}^k(\mathcal{F}_h \setminus \mathcal{F}_h^D)$.

An important property of the HHO method we exploit in the a posteriori error analysis is local conservation. For all $K \in \mathcal{T}_h$ and all $F \in \mathcal{F}_{\partial K}$, we define the flux

$$ \phi_{K,F}(\hat{u}_K) := -A_K \nabla R_K^{k+1}(\hat{u}_K) \cdot n_K|_F + A_K \frac{(k+1)^2}{h_K} \Pi_F^k(u_K|_F) - u_F \in \mathbb{P}^k(F). \quad (19) $$

Then, the following holds true [16]: (i) At every interface $F = \partial K_1 \cap \partial K_2 \in \mathcal{F}_h^i$, we have

$$ \phi_{K_1,F}(\hat{u}_{K_1}) + \phi_{K_2,F}(\hat{u}_{K_2}) = 0; \quad (20a) $$

(ii) At every Neumann boundary face $F = \partial K \cap \Gamma_N \in \mathcal{F}_h^N$, we have

$$ \phi_{K,F}(\hat{u}_K) + \Pi_F^k(g_N|_F) = 0. \quad (20b) $$

## 3.3 Stability and well-posedness

We equip the local HHO space $\hat{V}_K^k$ with the $H^1$-like seminorm such that, for all $\hat{v}_K \in \hat{V}_K^k$,

$$ |\hat{v}_K|_{\hat{V}_K^k}^2 := \|\nabla v_K\|_K^2 + \frac{(k+1)^2}{h_K} \|v_{\partial K} - \Pi_{\partial K}^k(v_K|_{\partial K})\|_{\partial K}^2. \quad (21) $$

**Lemma 3.2 (Local stability and boundedness).** There is a real number $\alpha > 0$, depending only on the mesh shape-regularity and the space dimension $d$, such that, for all $\hat{v}_K \in \hat{V}_K^k$ and all $K \in \mathcal{T}_h$,

$$ \alpha |\hat{v}_K|_{\hat{V}_K^k}^2 \le \|\nabla R_K^{k+1}(\hat{v}_K)\|_K^2 + S_{\partial K}(\hat{v}_K, \hat{v}_K) \le \alpha^{-1} |\hat{v}_K|_{\hat{V}_K^k}^2. \quad (22) $$

*Proof.* The proof proceeds as in [16] using the discrete trace inequality (3); this is the reason why the stabilization is scaled by $(k+1)^2$ and not just by $(k+1)$. $\square$

We equip the space $\hat{V}_{h0,D}^k$ with the norm

$$ \|\hat{v}_h\|_{\text{HHO}}^2 := \sum_{K \in \mathcal{T}_h} A_K |\hat{v}_K|_{\hat{V}_K^k}^2, \quad \forall \hat{v}_h \in \hat{V}_{h0,D}^k. \quad (23) $$

The fact that $\|\cdot\|_{\text{HHO}}$ defines indeed a norm on $\hat{V}_{h0,D}^k$ is shown, e.g., in [15, Lemma 1.6].
**Corollary 3.3** (Coercivity and well-posedness). The discrete bilinear form $a_h$ is coercive on $\hat{V}_{h0,D}^k$, and the discrete problem (17) is well-posed.

*Proof.* Summing (22) multiplied by $A_K$ over all the mesh cells shows the following coercivity and continuity properties:

$$ \alpha \|\hat{v}_h\|_{\text{H}0\text{H}}^2 \le a_h(\hat{v}_h, \hat{v}_h) \le \alpha^{-1} \|\hat{v}_h\|_{\text{H}0\text{H}}^2, \quad \forall \hat{v}_h \in \hat{V}_{h0,D}^k. \quad (24) $$

The well-posedness of (17) with homogeneous Dirichlet boundary conditions then follows from the Lax-Milgram lemma. For inhomogeneous Dirichlet boundary conditions, the proof follows by introducing a lifting, say $\hat{g}_h \in \hat{V}_{hg,D}^k$, such that $g_F = \Pi_F^k(g_D|_F)$, for all $F \in \mathcal{F}_h^D$. $\square$

# 4 hp-a priori error estimate

In this section, we establish our hp-a priori error estimate.

## 4.1 Approximation

For all $K \in \mathcal{T}_h$, we define the local reduction operator $\hat{\mathcal{I}}_K^k : H^1(K) \to \hat{V}_K^k$ such that, for all $v \in H^1(K)$,

$$ \hat{\mathcal{I}}_K^k(v) := (\Pi_K^{k+1}(v), \Pi_{\partial K}^k(v|_{\partial K})) \in \hat{V}_K^k. \quad (25) $$

Moreover, the elliptic projection $\mathcal{E}_K^{k+1} : H^1(K) \to \mathbb{P}^{k+1}(K)$ is defined such that

$$ \begin{gathered} (\nabla(\mathcal{E}_K^{k+1}(v) - v), \nabla w)_K = 0, \quad \forall w \in \mathbb{P}^{k+1}(K)/\mathbb{R}, \\ (\mathcal{E}_K^{k+1}(v) - v, 1)_K = 0. \end{gathered} \quad (26) $$

One readily verifies by proceeding as in [20, Lemma 3] that

$$ R_K^{k+1} \circ \hat{\mathcal{I}}_K^k = \mathcal{E}_K^{k+1}. \quad (27) $$

**Lemma 4.1** (Bound on stabilization). The following holds for all $K \in \mathcal{T}_h$ and all $v \in H^1(K)$:

$$ S_{\partial K}(\hat{\mathcal{I}}_K^k(v), \hat{\mathcal{I}}_K^k(v)) \le C(k+1) \|\nabla(v - \mathcal{I}_{BS,K}^{k+1}(v))\|_K^2. \quad (28) $$

*Proof.* Recalling the definition (12) of $S_{\partial K}$, the definition (25) of $\hat{\mathcal{I}}_K^k$ and since $\Pi_{\partial K}^k \circ \Pi_{\partial K}^k = \Pi_{\partial K}^k$, we have

$$ S_{\partial K}(\hat{\mathcal{I}}_K^k(v), \hat{\mathcal{I}}_K^k(v)) = \frac{(k+1)^2}{h_K} \|\Pi_{\partial K}^k((v - \Pi_K^{k+1}(v))|_{\partial K})\|_{\partial K}^2 \le \frac{(k+1)^2}{h_K} \|v - \Pi_K^{k+1}(v)\|_{\partial K}^2, $$

where we used the $L^2(\partial K)$-stability of $\Pi_{\partial K}^k$. Then, we invoke the approximation result on the $L^2$-orthogonal projection, see (4), and that $\Pi_K^{k+1} \circ \mathcal{I}_{BS,K}^{k+1} = \mathcal{I}_{BS,K}^{k+1}$, giving

$$ \begin{gathered} \frac{(k+1)^2}{h_K} \|v - \Pi_K^{k+1}(v)\|_{\partial K}^2 = \frac{(k+1)^2}{h_K} \|v - \mathcal{I}_{BS,K}^{k+1}(v) - \Pi_K^{k+1}(v - \mathcal{I}_{BS,K}^{k+1}(v))\|_{\partial K}^2 \\ \le C(k+1) \|\nabla(v - \mathcal{I}_{BS,K}^{k+1}(v))\|_{K}^2. \end{gathered} $$

Combining the above two bounds proves (28). $\square$

For all $K \in \mathcal{T}_h$ and all $v \in H^{1+s}(K)$, $s > \frac{1}{2}$, we consider the following norm:

$$ \|v\|_{\sharp,K}^2 := \|\nabla v\|_K^2 + \frac{h_K}{k+1} \|\nabla v\|_{\partial K}^2. \quad (29) $$

**Lemma 4.2** (Approximation). The following holds for all $K \in \mathcal{T}_h$ and all $v \in H^{1+s}(K)$, $s > \frac{1}{2}$:

$$ \|v - \mathcal{E}_K^{k+1}(v)\|_{\sharp,K} \le C(k+1)^{\frac{1}{2}} \|v - \mathcal{I}_{BS,K}^{k+1}(v)\|_{\sharp,K}. \quad (30) $$
Proof. Using the triangle inequality, we have

$$ \|v - \mathcal{E}_K^{k+1}(v)\|_{\sharp,K} \le \|v - \mathcal{I}_{\text{BS},K}^{k+1}(v)\|_{\sharp,K} + \|\mathcal{I}_{\text{BS},K}^{k+1}(v) - \mathcal{E}_K^{k+1}(v)\|_{\sharp,K}, $$

so that we only need to bound the second term on the right-hand side. Owing to the discrete trace inequality (3), we infer that

$$ \|\mathcal{E}_K^{k+1}(v) - \mathcal{I}_{\text{BS},K}^{k+1}(v)\|_{\sharp,K} \le C(k+1)^{\frac{1}{2}} \|\nabla(\mathcal{E}_K^{k+1}(v) - \mathcal{I}_{\text{BS},K}^{k+1}(v))\|_K, $$

and it remains to bound $\|\nabla(\mathcal{E}_K^{k+1}(v) - \mathcal{I}_{\text{BS},K}^{k+1}(v))\|_K$. Since $\mathcal{E}_K^{k+1} \circ \mathcal{I}_{\text{BS},K}^{k+1} = \mathcal{I}_{\text{BS},K}^{k+1}$, we infer that

$$ \|\nabla(\mathcal{E}_K^{k+1}(v) - \mathcal{I}_{\text{BS},K}^{k+1}(v))\|_K = \|\nabla(\mathcal{E}_K^{k+1}(v - \mathcal{I}_{\text{BS},K}^{k+1}(v)))\|_K \le \|\nabla(v - \mathcal{I}_{\text{BS},K}^{k+1}(v))\|_K, $$

where the last bound follows from the stability of the elliptic projection. Combining the above bounds completes the proof. □

## 4.2 Consistency

Let $u \in H_{g,D}^1(\Omega)$ be the exact solution to (2). We define the consistency error $\delta_h \in (\hat{V}_{h0,D}^k)'$ such that

$$ \langle \delta_h, \hat{w}_h \rangle := \ell_h(\hat{w}_h) - a_h(\hat{\mathcal{I}}_h^k(u), \hat{w}_h), \quad \forall \hat{w}_h \in \hat{V}_{h0,D}^k, \quad (31) $$

where $\langle \cdot, \cdot \rangle$ denotes the duality pairing between $(\hat{V}_{h0,D}^k)'$ and $\hat{V}_{h0,D}^k$, and where the global reduction operator $\hat{\mathcal{I}}_h^k: H^1(\Omega) \to \hat{V}_h^k$ is defined such that, for all $v \in H^1(\Omega)$,

$$ \hat{\mathcal{I}}_h^k(v) := (\Pi_K^{k+1}(v|_K))_{K \in \mathcal{T}_h}, (\Pi_F^k(v|_F))_{F \in \mathcal{F}_h} \in \hat{V}_h^k, \quad (32) $$

observing that $v$ is single-valued on every $F \in \mathcal{F}_h^i$. Notice that the local components of $\hat{\mathcal{I}}_h^k(v)$ attached to $K$ and its faces are $\hat{\mathcal{I}}_K^k(v|_K)$ for all $K \in \mathcal{T}_h$. The above definition of the consistency error is classical in the context of HHO methods. It is also rather natural since, as for most nonconforming discretization methods, it avoids defining any extension of the discrete bilinear form. Examples for other nonconforming methods can be found, e.g., in [26, Part VIII].

**Lemma 4.3 (Consistency).** Assume that $u \in H^{1+s}(\Omega)$ with $s > \frac{1}{2}$. The following holds true:

$$ \|\delta_h\|_{\text{HHO}'} := \sup_{\hat{w}_h \in \hat{V}_{h0,D}^k} \frac{|\langle \delta_h, \hat{w}_h \rangle|}{\|\hat{w}_h\|_{\text{HHO}}} \le C \left( \sum_{K \in \mathcal{T}_h} A_K(k+1) \|u - \mathcal{I}_{\text{BS},K}^{k+1}(u)\|_{\sharp,K}^2 \right)^{\frac{1}{2}}. \quad (33) $$

Proof. Let $\hat{w}_h \in \hat{V}_{h0,D}^k$. Using the definition of $\ell_h$ in (18), the PDE and the boundary conditions satisfied by the exact solution $u$, and integrating by parts cellwise, we infer that

$$ \ell_h(\hat{w}_h) = \sum_{K \in \mathcal{T}_h} A_K \left\{ (\nabla u, \nabla w_K)_K - (\nabla u \cdot \mathbf{n}_K, w_K)_{\partial K} + (\nabla u \cdot \mathbf{n}_K, w_{\partial K})_{\partial K^\mathrm{N}} \right\}. $$

The assumption $u \in H^{1+s}(\Omega)$ with $s > \frac{1}{2}$ implies that $(\nabla u \cdot \mathbf{n}_K)|_{\partial K}$ is meaningful in $L^2(\partial K)$ and single-valued at every mesh interface. Moreover, since $w_{\partial K}$ is single-valued on $\partial K^i$ and vanishes on $\partial K^D$, we infer that

$$ \ell_h(\hat{w}_h) = \sum_{K \in \mathcal{T}_h} A_K \left\{ (\nabla u, \nabla w_K)_K - (\nabla u \cdot \mathbf{n}_K, w_K - w_{\partial K})_{\partial K} \right\}. $$

Since $a_h$ is assembled cellwise and the local components of $\hat{\mathcal{I}}_h^k(u)$ are $\hat{\mathcal{I}}_K^k(u|_K)$ for all $K \in \mathcal{T}_h$, we infer that $a_h(\hat{\mathcal{I}}_h^k(u), \hat{w}_h) = \sum_{K \in \mathcal{T}_h} a_K(\hat{\mathcal{I}}_K^k(u|_K), \hat{w}_K)$. Using the definition (13) of $a_K$, the definition (10) of $R_K^{k+1}(\hat{w}_K)$ and the identity (27) leads to

$$ a_h(\hat{\mathcal{I}}_h^k(u), \hat{w}_h) = \sum_{K \in \mathcal{T}_h} A_K \left\{ (\nabla \mathcal{E}_K^{k+1}(u|_K), \nabla w_K)_K - (\nabla \mathcal{E}_K^{k+1}(u|_K) \cdot \mathbf{n}_K, w_K - w_{\partial K})_{\partial K} + S_{\partial K}(\hat{\mathcal{I}}_K^k(u|_K), \hat{w}_K) \right\}. $$
Defining the function $\eta$ cellwise as $\eta|_K := u|_K - \mathcal{E}_K^{k+1}(u|_K)$ for all $K \in \mathcal{T}_h$, we infer that

$$ \langle \delta_h, \hat{w}_h \rangle = \sum_{K \in \mathcal{T}_h} A_K \left\{ - (\nabla \eta \cdot n_K, w_K - w_{\partial K})_{\partial K} - S_{\partial K}(\hat{I}_K^k(u|_K), \hat{w}_K) \right\}, \quad (34) $$

where we used that $(\nabla\eta, \nabla w_K)_K = 0$. Let us denote by $\mathcal{T}_{1,K}$ the first term and by $\mathcal{T}_{2,K}$ the second term inside braces on the right-hand side of (34). We bound $\mathcal{T}_{1,K}$ by the Cauchy–Schwarz inequality, the triangle inequality, the trivial bound $\|w_K - \Pi_{\partial K}^k(w_K|_{\partial K})\|_{\partial K} \le \|w_K - \Pi_K^k(w_K)\|_{\partial K}$, and the approximation result (4) on the $L^2$-projection. Recalling the definition (29) of the $\|\cdot\|_{\sharp,K}$-norm, this yields

$$ \begin{aligned} |\mathcal{T}_{1,K}| &\le \|\nabla\eta\|_{\partial K} (\|\Pi_{\partial K}^k(w_K|_{\partial K}) - w_{\partial K}\|_{\partial K} + \|w_K - \Pi_{\partial K}^k(w_K|_{\partial K})\|_{\partial K}) \\ &\le \left(\frac{h_K}{k+1}\right)^{\frac{1}{2}} \|\nabla\eta\|_{\partial K} \left(\frac{k+1}{h_K}\right)^{\frac{1}{2}} \left(\|\Pi_{\partial K}^k(w_K|_{\partial K}) - w_{\partial K}\|_{\partial K} + \|w_K - \Pi_K^k(w_K)\|_{\partial K}\right) \\ &\le \|\eta\|_{\sharp,K} \left(\left(\frac{k+1}{h_K}\right)^{\frac{1}{2}} \|\Pi_{\partial K}^k(w_K|_{\partial K}) - w_{\partial K}\|_{\partial K} + C \|\nabla w_K\|_K\right) \le C \|\eta\|_{\sharp,K} |\hat{w}_K|_{\hat{V}_K^k}. \end{aligned} $$

Moreover, owing to (28) and the upper bound in (22), we have

$$ |\mathcal{T}_{2,K}| \le S_{\partial K}(\hat{\mathcal{I}}_K^k(u|_K), \hat{\mathcal{I}}_K^k(u|_K))^{\frac{1}{2}} S_{\partial K}(\hat{w}_K, \hat{w}_K)^{\frac{1}{2}} \le C(k+1)^{\frac{1}{2}} \|\nabla(u - \mathcal{I}_{\text{BS},K}^{k+1}(u))\|_K |\hat{w}_K|_{\hat{V}_K^k}. $$

Altogether, this implies that

$$ |\langle \delta_h, \hat{w}_h \rangle| \le C \left( \sum_{K \in \mathcal{T}_h} A_K \left\{ \|\eta\|_{\sharp,K}^2 + (k+1) \|\nabla(u - \mathcal{I}_{\text{BS},K}^{k+1}(u))\|_K^2 \right\} \right)^{\frac{1}{2}} \|\hat{w}_h\|_{\text{HHO}}. $$

Invoking Lemma 4.2 to bound $\|\eta\|_{\sharp,K}$ completes the proof. $\square$

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

## 4.3 Error estimate

We are now ready to establish our main result concerning the hp-a priori error analysis. The estimate is $\frac{1}{2}$-order p-suboptimal.

**Theorem 4.4 (hp-a priori error estimate).** Let $u$ be the weak solution to (2), and let $\hat{u}_h$ be the discrete solution to (17). Assume that $u \in H^{1+s}(\Omega)$ with $s > \frac{1}{2}$. The following holds:

$$ \sum_{K \in \mathcal{T}_h} A_K \left\{ \|\nabla(u - R_K^{k+1}(\hat{u}_K))\|_K^2 + \|\nabla(u - u_K)\|_K^2 + S_{\partial K}(\hat{u}_K, \hat{u}_K) \right\} \le C \sum_{K \in \mathcal{T}_h} A_K(k+1) \|u - \mathcal{I}_{\text{BS},K}^{k+1}(u)\|_{\sharp,K}^2. \quad (35a) $$

Moreover, if $u \in H^{l+1}(\mathcal{T}_h; \mathbb{R})$, with $l \in \{1, \dots, k+1\}$, we have

$$ \sum_{K \in \mathcal{T}_h} A_K \left\{ \|\nabla(u - R_K^{k+1}(\hat{u}_K))\|_K^2 + \|\nabla(u - u_K)\|_K^2 + S_{\partial K}(\hat{u}_K, \hat{u}_K) \right\} \le C \sum_{K \in \mathcal{T}_h} A_K(k+1) \left(\frac{h_K}{k+1}\right)^{2l} \|u\|_{H^{l+1}(K)}^2. \quad (35b) $$

**Proof.** (i) Proof of (35a). We set $\hat{e}_h := \hat{\mathcal{I}}_h^k(u) - \hat{u}_h$ and observe that $\hat{e}_h \in \hat{V}_{h0,D}^k$. Moreover, since $a_h(\hat{e}_h, \hat{e}_h) = -\langle \delta_h, \hat{e}_h \rangle$, the coercivity property (24) implies that

$$ \alpha \|\hat{e}_h\|_{\text{HHO}}^2 \le a_h(\hat{e}_h, \hat{e}_h) = -\langle \delta_h, \hat{e}_h \rangle \le \|\delta_h\|_{\text{HHO}'} \|\hat{e}_h\|_{\text{HHO}}, $$

so that $\|\hat{e}_h\|_{\text{HHO}} \le \frac{1}{\alpha} \|\delta_h\|_{\text{HHO}'}$. Since $\sum_{K \in \mathcal{T}_h} A_K \{\|\nabla R_K^{k+1}(\hat{e}_K)\|_K^2 + S_{\partial K}(\hat{e}_K, \hat{e}_K)\} \le \alpha^{-1} \|\hat{e}_h\|_{\text{HHO}}^2$ follows from (24), we infer from Lemma 4.3 that

$$ \sum_{K \in \mathcal{T}_h} A_K \left\{ \|\nabla R_K^{k+1}(\hat{e}_K)\|_K^2 + S_{\partial K}(\hat{e}_K, \hat{e}_K) \right\} \le C \sum_{K \in \mathcal{T}_h} A_K(k+1) \|u - \mathcal{I}_{\text{BS},K}^{k+1}(u)\|_{\sharp,K}^2. $$

Since $u - R_K^{k+1}(\hat{u}_K) = (u - \mathcal{E}_K^{k+1}(u)) + R_K^{k+1}(\hat{e}_K)$ because $\hat{u}_K = \hat{\mathcal{I}}_K^k(u) - \hat{e}_K$, the triangle inequality combined with Lemma 4.1, Lemma 4.2, and the above bound proves that

$$ \sum_{K \in \mathcal{T}_h} A_K \left\{ \|\nabla(u - R_K^{k+1}(\hat{u}_K))\|_K^2 + S_{\partial K}(\hat{u}_K, \hat{u}_K) \right\} \le C \sum_{K \in \mathcal{T}_h} A_K(k+1) \|u - \mathcal{I}_{\text{BS},K}^{k+1}(u)\|_{\sharp,K}^2. $$
Moreover, using $\|\nabla(u - u_K)\|_K^2 \le 2(\|\nabla(u - R_K^{k+1}(\hat{u}_K))\|_K^2 + \|\nabla(u_K - R_K^{k+1}(\hat{u}_K))\|_K^2)$ and (14), we infer that

$$ \sum_{K \in \mathcal{T}_h} A_K \|\nabla(u - u_K)\|_K^2 \le 2 \sum_{K \in \mathcal{T}_h} A_K \left\{ \|\nabla(u - R_K^{k+1}(\hat{u}_K))\|_K^2 + C S_{\partial K}(\hat{u}_K, \hat{u}_K) \right\} \\ \le C \sum_{K \in \mathcal{T}_h} A_K (k+1) \|u - \mathcal{I}_{\text{BS},K}^{k+1}(u)\|_{\sharp,K}^2. $$

This completes the proof of (35a).

(ii) The proof of (35b) follows by invoking on all $K \in \mathcal{T}_h$ the approximation properties of $\mathcal{I}_{\text{BS},K}^{k+1}$ (see Lemma 2.3), together with a multiplicative trace inequality (see, e.g., [25, Sec. 12.3.2]) to estimate the term involving the normal derivative in $\|u - \mathcal{I}_{\text{BS},K}^{k+1}(u)\|_{\sharp,K}$. $\square$

# 5 hp-a posteriori error analysis

In this section, we perform the hp-a posteriori error analysis of the above HHO discretization. We establish both (global) upper and (local) lower error bounds. Since we are interested in the hp-a posteriori error analysis, we only consider in what follows the case $k \ge 1$ so that we can invoke the (modified) Karkulik-Melenk interpolation operator. We shall also assume that the Dirichlet datum satisfies $g_D \in H^{\frac{1}{2}}(\partial\Omega)$. For all $K \in \mathcal{T}_h$, we define the following data oscillation terms:

$$ O_K(f) := A_K^{-\frac{1}{2}} \left( \frac{h_K}{k+1} \right) \|f - \Pi_K^{k+1}(f)\|_K, \quad (36a) $$

$$ O_K(g_N) := A_K^{-\frac{1}{2}} \left( \frac{h_K}{k+1} \right)^{\frac{1}{2}} \|g_N - \Pi_{\partial K}^k (g_N|_{\partial K})\|_{\partial K^N}, \quad (36b) $$

$$ O_K(g_D) := A_K^{\frac{1}{2}} \left( \frac{h_K}{k+1} \right)^{\frac{1}{2}} \| \nabla (g_D - \Pi_{\partial K}^{k+1} (g_D|_{\partial K})) \times \mathbf{n}_{\Omega} \|_{\partial K^{\mathrm{D}}}. \quad (36c) $$

It is useful to define some contrast factors related to the diffusion coefficient $A$. For all $K \in \mathcal{T}_h$ and all $\mathbf{a} \in \mathcal{V}_h$, we set

$$ \chi_K(A) := A_K \max_{\widehat{K} \in \text{es}(K)} A_{\widehat{K}}^{-1}, \quad \chi_a(A) := A_a^{\sharp} (A_a^{\flat})^{-1}, \quad (37) $$

with $A_a^{\sharp} := \max_{K \in \mathcal{T}_h} A_K$ and $A_a^{\flat} := \min_{K \in \mathcal{T}_h} A_K$. Moreover, it is convenient to set

$$ A_F^{\flat} := \min(A_{K_1}, A_{K_2}), \quad \forall F = \partial K_1 \cap \partial K_2 \in \mathcal{F}_h^i, \quad (38a) $$

$$ A_F^{\flat} := A_K, \quad \forall F = \partial K \cap \partial \Omega \in \mathcal{F}_h^{\flat}. \quad (38b) $$

For all $K \in \mathcal{T}_h$, we also define $A_{\partial K}^{\flat}|_F := A_F^{\flat}$ for all $F \in \mathcal{F}_{\partial K}$.

## 5.1 Global upper error bound

We decompose the error into two components as follows:

$$ e := u - u_{\mathcal{T}_h} = (u - u_c) + (u_c - u_{\mathcal{T}_h}) =: e_c + e_d, \quad u_c \in H_{g,D}^{1}(\Omega), \quad (39) $$

where $u_{\mathcal{T}_h}$ denotes the cell component of the HHO solution such that $u_{\mathcal{T}_h}|_K := u_K$ for all $K \in \mathcal{T}_h$. The function $u_c$ is constructed from $u_{\mathcal{T}_h}$ as detailed in Section 5.1.2. The precise definition of $u_c$ is irrelevant for bounding $e_c$ (we only use that $e_c \in H_{0,D}^{1}(\Omega)$), and is only relevant for bounding $e_d$. We call $e_c$ the conforming error and $e_d$ the nonconforming error.

### 5.1.1 Bound on conforming error $e_c$

In this section, we derive two bounds on the conforming error $e_c$. The first bound avoids the normal flux jump (classically considered in the context of finite elements) and is $\frac{1}{2}$-order p-suboptimal. The second bound includes the normal flux jump and is p-optimal. The first bound is, however, interesting in its own right. We will also see that, in the context of the lower error bound, the normal flux jump leads to $\frac{1}{2}$-order p-suboptimality anyway. The proofs of the following two results are postponed to Section 7.
**Lemma 5.1** (Conforming error (p-suboptimal bound)). The following holds:

$$ \|A^{\frac{1}{2}} \nabla e_c\|_{\Omega} \le C_{c,1} \left\{ \sum_{K \in \mathcal{T}_h} \left\{ A_K^{-1} \left( \frac{h_K}{k+1} \right)^2 \|\Pi_K^{k+1}(f) + A_K \Delta R_K^{k+1}(\hat{u}_K)\|_K^2 + A_K(k+1) S_{\partial K}(\hat{u}_K, \hat{u}_K) \right. \right. \\ \left. \left. + O_K(f)^2 + O_K(g_N)^2 \right\} \right\}^{\frac{1}{2}} + \|A^{\frac{1}{2}} \nabla_{\mathcal{T}_h} e_d\|_{\Omega}, \quad (40) $$

where the constant $C_{c,1}$ depends on the mesh shape-regularity and on $\max_{K \in \mathcal{T}_h} \chi_K(A)^{\frac{1}{2}}$, but is independent of $h$ and $k$.

**Lemma 5.2** (Conforming error (p-optimal bound)). The following holds:

$$ \|A^{\frac{1}{2}} \nabla e_c\|_{\Omega} \le C_{c,2} \left\{ \sum_{K \in \mathcal{T}_h} \left\{ A_K^{-1} \left( \frac{h_K}{k+1} \right)^2 \|\Pi_K^{k+1}(f) + A_K \Delta R_K^{k+1}(\hat{u}_K)\|_K^2 + A_K S_{\partial K}(\hat{u}_K, \hat{u}_K) \right. \right. \\ \left. + A_K^{-1} \left( \frac{h_K}{k+1} \right) \left( \|\llbracket A \nabla R_{\mathcal{T}_h}^{k+1}(\hat{u}_h) \rrbracket_{\partial K} \cdot \mathbf{n}_K \|_{\partial K^i}^2 + \|A_K \nabla R_K^{k+1}(\hat{u}_K) \cdot \mathbf{n}_K - \Pi_{\partial K}^k(g_N)\|_{\partial K^N}^2 \right) \right. \\ \left. + O_K(f)^2 + O_K(g_N)^2 \right\} \right\}^{\frac{1}{2}} + \|A^{\frac{1}{2}} \nabla_{\mathcal{T}_h} e_d\|_{\Omega}, \quad (41) $$

where the constant $C_{c,2}$ depends on the mesh shape-regularity and on $\max_{K \in \mathcal{T}_h} \chi_K(A)^{\frac{1}{2}}$, but is independent of $h$ and $k$.

### 5.1.2 Bound on nonconforming error $e_d$

We start by defining the function $u_c \in H_{g,D}^1(\Omega)$ introduced in (39). To this purpose, in the spirit of [30], we solve local minimization problems in $H^1(\omega_a)$ with suitable boundary conditions for every mesh vertex $a \in V_h$ (recall that $\omega_a$ denotes the star associated with the vertex $a$). Let $\psi_a$ be the hat basis function equal to 1 at $a$ and having a support in the vertex star $\omega_a$. Recall that the hat basis functions satisfy the following partition-of-unity property:

$$ \sum_{a \in \mathcal{V}_h} \psi_a = 1. \quad (42) $$

**Definition 5.3** (Patchwise and global potential reconstruction). For all $a \in V_h$, let $u_c^a \in H_{g,D}^1(\omega_a)$ be the solution of the following well-posed problem:

$$ (A \nabla u_c^a, \nabla v_a)_{\omega_a} = (A \nabla_{\mathcal{T}_h}(\psi_a u_{\mathcal{T}_h}), \nabla v_a)_{\omega_a}, \quad \forall v_a \in H_{0,D}^1(\omega_a), \quad (43) $$

with

$$ H_{g,D}^1(\omega_a) := \{v \in H^1(\omega_a) \mid v|_{\partial \omega_a \cap \Gamma_D} = \psi_a g_D \text{ and } v|_{\partial \omega_a \cap \Omega} = 0\}. \quad (44) $$

An equivalent definition is

$$ u_c^a := \arg \min_{\rho_a \in H_{g,D}^1(\omega_a)} \|A^{\frac{1}{2}}(\nabla \rho_a - \nabla_{\mathcal{T}_h}(\psi_a u_{\mathcal{T}_h}))\|_{\omega_a}. \quad (45) $$

Then, extending $u_c^a$ by zero to $\Omega$, we set

$$ u_c := \sum_{a \in \mathcal{V}_h} u_c^a. \quad (46) $$

Notice that we indeed have $u_c \in H_{g,D}^1(\Omega)$ as required in (39); this follows from the partition-of-unity property (42) and the definition (44) of $H_{g,D}^1(\omega_a)$. The proof of the following result is postponed to Section 7.
**Lemma 5.4 (Local nonconforming error).** Let $u_c^a$ solve (43) and set $e_d^a := u_c^a - \psi_a u_{\mathcal{T}_h}$ for all $a \in \mathcal{V}_h$. The following holds:

$$ \|A^{\frac{1}{2}} \nabla_{\mathcal{T}_h} e_d^a\|_{\omega_a} \le C_d^a \left\{ \sum_{K \in \mathcal{T}_a} \left\{ A_K S_{\partial K}(\hat{u}_K, \hat{u}_K) + O_K(g_D)^2 \right\} + \sum_{F \in \mathcal{F}_a \cap \mathcal{F}_h^i} A_F^b \left( \frac{h_F}{k+1} \right) \|\|\nabla u_{\mathcal{T}_h}\|\|_F \times \mathbf{n}_F \|\|_F^2 \right. \\ \left. + \sum_{F \in \mathcal{F}_a \cap \mathcal{F}_h^D} A_F^b \left( \frac{h_F}{k+1} \right) \|\nabla(u_{\mathcal{T}_h} - \Pi_{\partial K}^{k+1}(g_D|_{\partial K})) \times \mathbf{n}_{\Omega} \|\|_F^2 \right\}^{\frac{1}{2}}, \quad (47) $$

where the constant $C_d^a$ depends on the mesh shape-regularity and on $\chi_a(A)^{\frac{1}{2}}$, but is independent of $h$ and $k$.

**Corollary 5.5 (Global bound on nonconforming error).** The following holds:

$$ \|A^{\frac{1}{2}} \nabla_{\mathcal{T}_h} e_d\|_{\Omega} \le C_d \left\{ \sum_{K \in \mathcal{T}_h} \left\{ A_K S_{\partial K}(\hat{u}_K, \hat{u}_K) + O_K(g_D)^2 + A_{\partial K}^b \left( \frac{h_K}{k+1} \right) \|\|\nabla u_{\mathcal{T}_h}\|\|_{\partial K} \times \mathbf{n}_K \|\|_{\partial K}^2 \right. \right. \\ \left. \left. + A_K \left( \frac{h_K}{k+1} \right) \|\nabla(u_K - \Pi_{\partial K}^{k+1}(g_D|_{\partial K})) \times \mathbf{n}_{\Omega} \|\|_{\partial K}^2 \right\} \right\}^{\frac{1}{2}}, \quad (48) $$

where the constant $C_d$ depends on the mesh shape-regularity and on $\max_{a \in \mathcal{V}_h} \chi_a(A)^{\frac{1}{2}}$, but is independent of $h$ and $k$.

*Proof.* Using that $u_c|_K = \sum_{a \in \mathcal{V}_K} u_c^a|_K$ for all $K \in \mathcal{T}_h$, where $\mathcal{V}_K$ is the set of vertices of $K$, the partition-of-unity property (42), and the triangle inequality, we infer that

$$ \|A^{\frac{1}{2}} \nabla_{\mathcal{T}_h} e_d\|_{\Omega}^2 = \sum_{K \in \mathcal{T}_h} A_K \|\nabla(u_c - u_K)\|_K^2 = \sum_{K \in \mathcal{T}_h} A_K \left\| \sum_{a \in \mathcal{V}_K} \nabla(u_c^a - \psi_a u_K) \right\|_K^2 \\ \le \sum_{K \in \mathcal{T}_h} \sum_{a \in \mathcal{V}_K} (d+1) A_K \|\nabla(u_c^a - \psi_a u_K)\|_K^2 = \sum_{a \in \mathcal{V}_h} (d+1) \|A^{\frac{1}{2}} \nabla_{\mathcal{T}_h} e_d^a\|_{\omega_a}^2. $$

Invoking Lemma 5.4 and using that $h_F \le h_K$ for all $F \in \partial K$ and all $K \in \mathcal{T}_h$ completes the proof. $\square$

**Remark 5.6 (Constant $C_d$).** The constant $C_d$ in (48) does not depend on the topology of $\Omega$. The reason is that we do not invoke a global Helmholtz decomposition in $\Omega$, but instead invoke a local decomposition in each vertex star $\omega_a$ (see Section 7.3 for further insight).

### 5.1.3 Main result

For all $K \in \mathcal{T}_h$, we define the following error indicators:

$$ \eta_{K,\mathrm{res}} := A_K^{-\frac{1}{2}} \left( \frac{h_K}{k+1} \right) \| \Pi_K^{k+1}(f) + A_K \Delta R_K^{k+1}(\hat{u}_K) \|_K, \quad (49a) $$

$$ \eta_{K,\mathrm{sta}} := A_K^{\frac{1}{2}} S_{\partial K}(\hat{u}_K, \hat{u}_K)^{\frac{1}{2}}, \quad (49b) $$

$$ \eta_{K,\mathrm{tan}} := (A_{\partial K}^b)^{\frac{1}{2}} \left( \frac{h_K}{k+1} \right)^{\frac{1}{2}} \left\{ \|\|\nabla u_{\mathcal{T}_h}\|\|_{\partial K} \times \mathbf{n}_K \|\|_{\partial K}^2 + \|\nabla(u_K - \Pi_{\partial K}^{k+1}(g_D|_{\partial K})) \times \mathbf{n}_{\Omega} \|\|_{\partial K^D} \right\}, \quad (49c) $$

$$ \eta_{K,\mathrm{nor}} := A_K^{-\frac{1}{2}} \left( \frac{h_K}{k+1} \right)^{\frac{1}{2}} \left\{ \|\|\nabla R_{\mathcal{T}_h}^{k+1}(\hat{u}_h)\|\|_{\partial K} \cdot \mathbf{n}_K \|\|_{\partial K}^2 + \|A_K \nabla R_K^{k+1}(\hat{u}_K) \cdot \mathbf{n}_{\Omega} - \Pi_{\partial K}^k (g_N|_{\partial K}) \|\|_{\partial K}^N \right\}, \quad (49d) $$

$$ O_{K,\mathrm{dat}} := O_K(f) + O_K(g_N) + O_K(g_D), \quad (49e) $$

where the three data oscillation terms on the right-hand side of (49e) are defined in (36).
**Theorem 5.7 (hp-upper error bound).** Let $u$ be the weak solution to (2), and let $\hat{u}_h$ be the discrete solution to (17). The following holds:

$$ \sum_{K \in \mathcal{T}_h} \{ \|A^{\frac{1}{2}} \nabla e\|_K^2 + A_K S_{\partial K}(\hat{u}_K, \hat{u}_K) \} \le C_u \left\{ \sum_{K \in \mathcal{T}_h} \left\{ \eta_{K, \text{res}}^2 + \eta_{K, \text{tan}}^2 + \eta_{K, \text{sta}}^2 + O_K^2 \right\} \right. \\ \left. + \min \left( \sum_{K \in \mathcal{T}_h} k \eta_{K, \text{sta}}^2, \sum_{K \in \mathcal{T}_h} \eta_{K, \text{nor}}^2 \right) \right\}, \quad (50) $$

where $C_u$ depends on the mesh shape-regularity and on $\max_{K \in \mathcal{T}_h} \chi_K(A)^{\frac{1}{2}}$ but is independent of $h$ and $k$. Proof. Combining the bounds from Lemma 5.1 and Lemma 5.2 gives

$$ \|A^{\frac{1}{2}} \nabla e_c\|_{\Omega} \le C \left\{ \sum_{K \in \mathcal{T}_h} \left\{ \eta_{K, \text{res}}^2 + \eta_{K, \text{sta}}^2 + O_K(f)^2 + O_K(g_N)^2 \right\} + \min \left( \sum_{K \in \mathcal{T}_h} k \eta_{K, \text{sta}}^2, \sum_{K \in \mathcal{T}_h} \eta_{K, \text{nor}}^2 \right) \right\}^{\frac{1}{2}}, $$

where the constant $C$ depends on $\max_{K \in \mathcal{T}_h} \chi_K(A)^{\frac{1}{2}}$. Moreover, we can rewrite the bound of Corollary 5.5 as follows:

$$ \|A^{\frac{1}{2}} \nabla_{\mathcal{T}_h} e_d\|_{\Omega} \le C \left\{ \sum_{K \in \mathcal{T}_h} \left\{ \eta_{K, \text{sta}}^2 + \eta_{K, \text{tan}}^2 + O_K(g_D)^2 \right\} \right\}^{\frac{1}{2}}, $$

where the constant $C$ depends on $\max_{a \in \mathcal{V}_h} \chi_a(A)^{\frac{1}{2}}$. Combining the above two bounds and since $\max_{a \in \mathcal{V}_h} \chi_a(A) \le \max_{K \in \mathcal{T}_h} \chi_K(A)$ proves that $\sum_{K \in \mathcal{T}_h} \|A^{\frac{1}{2}} \nabla e\|_K^2$ is bounded by the right-hand side of (50). Since $\sum_{K \in \mathcal{T}_h} A_K S_{\partial K}(\hat{u}_K, \hat{u}_K) = \sum_{K \in \mathcal{T}_h} \eta_{K, \text{sta}}^2$, the proof is complete. $\square$

**Remark 5.8 (Estimator without normal flux jump).** Notice that (50) implies that

$$ \sum_{K \in \mathcal{T}_h} \{ \|A^{\frac{1}{2}} \nabla e\|_K^2 + A_K S_{\partial K}(\hat{u}_K, \hat{u}_K) \} \le C_u \sum_{K \in \mathcal{T}_h} \{ \eta_{K, \text{res}}^2 + \eta_{K, \text{tan}}^2 + (k+1)\eta_{K, \text{sta}}^2 + O_K^2 \}. $$

This upper bound does not contain the normal flux jump, which is often the dominant component of the error estimator for $H^1$-conforming FEM [12]. The price to pay is a $\frac{1}{2}$-order p-suboptimality for the stabilization term. Our numerical experiments confirm that the normal flux jump term does not dominate the total error estimator.

## 5.2 Local lower error bound

In this section, we establish a local lower error bound. Specifically, we bound the local error indicators $\eta_{K, \text{res}}$, $\eta_{K, \text{nor}}$, and $\eta_{K, \text{tan}}$, for all $K \in \mathcal{T}_h$, in terms of the error $e = u - u_{\mathcal{T}_h}$ in (a neighborhood of) $K$ and the data oscillation indicators defined in (36). We do not bound the local error indicator $\eta_{K, \text{sta}}$ since it is present on both sides of the upper error bound (50). This is classical in a posteriori error estimates for nonconforming methods. The proof of the following result is postponed to Section 7.

**Theorem 5.9 (hp-local lower error bound).** The following holds for all $K \in \mathcal{T}_h$:

$$ \eta_{K, \text{res}} \le C_1(k+1) \left( \|A^{\frac{1}{2}} \nabla e\|_K + (A_K S_{\partial K}(\hat{u}_K, \hat{u}_K))^{\frac{1}{2}} + O_K(f) \right), \quad (51a) $$

$$ \eta_{K, \text{nor}} \le C_1(k+1)^{\frac{1}{2}} \left\{ \sum_{K \in \omega_K} A_K S_{\partial K}(\hat{u}_K, \hat{u}_K) \right\}^{\frac{1}{2}}, \quad (51b) $$

$$ \eta_{K, \text{tan}} \le C_1(k+1)^{\frac{3}{2}} \left\{ \sum_{K \in \omega_K} \|A^{\frac{1}{2}} \nabla e\|_K^2 \right\}^{\frac{1}{2}}, \quad (51c) $$

where $\omega_K$ collects all the mesh cells (including $K$) sharing at least an interface with $K$, and the constant $C_1$ depends on the mesh shape-regularity and on the diffusion contrast factor $\chi'_K(A) := A_K^{-1} \max_{K' \in \omega_K} A_{K'}$, but is independent of $h$ and $k$.

**Remark 5.10 (Bound on normal flux jump).** We observe that the upper bound on $\eta_{K, \text{nor}}$ has only a $\frac{1}{2}$-order p-suboptimality. This rather sharp result is achieved by exploiting the local conservation property of the HHO method, and is in contrast with the upper bound that could be obtained using bubble function techniques and which would feature a $\frac{3}{2}$-order p-suboptimality (details not shown for brevity).
<figure><img src="image_4.png" /><figcaption>Figure 1: Example 1. Energy error and a posteriori error estimator for $k \in \{0, 1, 2, 3\}$ as a function of DoFs (left) and effectivity index as a function of DoFs (right).</figcaption></figure>

# 6 Numerical examples

In this section, we present numerical examples to illustrate our theoretical results. Our goal here is to illustrate the main findings of the above theoretical analysis and also to illustrate how the present a posteriori error estimators can be used to drive an h-adaptive procedure that behaves in a robust way with respect to the polynomial degree. A further step forward, which we leave to future work, is to devise and test a full hp-adaptive procedure. This entails substantial algorithmic developments, in particular to set up a criterion to select between h- and p-refinement locally. We notice that hp-adaptive procedures are already available in various contexts in the literature, as reflected, e.g., in [45, 34, 17].

## 6.1 Example 1: Convergence rates for smooth solution

We select $f, A := I_{2 \times 2}$ and Dirichlet boundary conditions on the unit square $\Omega := (-1, 1)^2$, so that the exact solution is

$$u(x, y) := \sin(\pi x) \sin(\pi y). \qquad (52)$$

We employ the polynomial degrees $k \in \{0, \dots, 3\}$ and a sequence of successively refined triangular meshes consisting of $\{32, 128, 512, 2048, 8192\}$ cells.

Let us first verify the convergence rates obtained with the HHO methods with $k \in \{0, \dots, 3\}$. We measure the error in the energy norm defined as $\{\sum_{K \in \mathcal{T}_h} A_K \|\nabla(u - u_K)\|_K^2 + S_{\partial K}(\hat{u}_K, \hat{u}_K)\}^{\frac{2}{3}}$. The energy error and the a posteriori estimator on the right-hand side of (50) (with constant $C_u$ set to one) are reported in Figure 1. The rates are computed as a function of DoFs, which denotes the total number of globally coupled discrete unknowns (that is, the total number of face unknowns except those located on the boundary faces). We observe that the energy error and the a posteriori estimator both converge at the optimal rate $O(\text{DoFs}^{-\frac{k+1}{2}})$. The convergence rate of the energy error is optimal in view of the result of Theorem 4.4. Moreover, the effectivity index, defined as the ratio of the a posteriori estimator to the energy error, remains well behaved as a function of DoFs. The effectivity index takes values between 2 and 2.8 for $k \ge 1$, whereas the effectivity index is almost 3 for $k = 0$.

As the results in Theorem 5.7 (upper error bound) and Theorem 5.9 (lower error bound) differ by an algebraic rate in the polynomial degree $k$, we investigate the dependence of the effectivity index on $k$. In Figure 2, we report the effectivity index as a function of the polynomial degree $k \in \{1, \dots, 9\}$ on a mesh consisting of 128 triangular cells. We observe an algebraic rate of $p^{\frac{1}{2}}$, which matches the statement in Theorem 5.7.
<figure><img src="image_5.png" /><figcaption>Figure 2: Example 1. Effectivity index with $k \in \{1, \dots, 9\}$ on a mesh composed of 128 cells.</figcaption></figure>

Finally, we compare the relative contributions of the various terms composing the a posteriori error estimator (all in percentage). Setting $\eta_X := \{\sum_{K \in \mathcal{T}} \eta_{K,X}^2\}^{\frac{1}{2}}$ for $X \in \{\text{res}, \text{sta}, \text{nor}, \text{tan}\}$, we report in Table 1 the relative contribution of the residual estimator $\eta_{\text{res}}$, the stabilization estimator $\eta_{\text{sta}}$, the normal flux jump estimator $\eta_{\text{nor}}$, and the tangential flux jump estimator $\eta_{\text{tan}}$, for polynomial degrees $k \in \{0, 1, 2\}$. For $k \in \{1, 2\}$, the residual estimator dominates the total estimator (by more than 60%), followed by $\eta_{\text{nor}}$ (about 20%), while $\eta_{\text{tan}}$ and $\eta_{\text{sta}}$ are both below 10%. For $k = 0$, the residual dominates the total estimator (by more than 50%), followed by $\eta_{\text{tan}}$ (25%), $\eta_{\text{nor}}$ (15%), and $\eta_{\text{sta}}$ (6%).

<table><caption>Table 1: Example 1. Relative contribution (in %) of the various terms composing the a posteriori error estimator for $k \in \{0, 1, 2\}$.</caption><thead><tr><th rowspan="2"># cell</th><th colspan="4">k = 0</th><th colspan="4">k = 1</th><th colspan="4">k = 2</th></tr><tr><th>η<sub>res</sub></th><th>η<sub>sta</sub></th><th>η<sub>nor</sub></th><th>η<sub>tan</sub></th><th>η<sub>res</sub></th><th>η<sub>sta</sub></th><th>η<sub>nor</sub></th><th>η<sub>tan</sub></th><th>η<sub>res</sub></th><th>η<sub>sta</sub></th><th>η<sub>nor</sub></th><th>η<sub>tan</sub></th></tr></thead><tbody><tr><td>128</td><td>54</td><td>6</td><td>15</td><td>25</td><td>64</td><td>9</td><td>19</td><td>8</td><td>66</td><td>8</td><td>19</td><td>7</td></tr><tr><td>512</td><td>54</td><td>6</td><td>15</td><td>25</td><td>62</td><td>9</td><td>22</td><td>7</td><td>65</td><td>8</td><td>20</td><td>7</td></tr><tr><td>2048</td><td>53</td><td>6</td><td>15</td><td>26</td><td>62</td><td>9</td><td>24</td><td>5</td><td>65</td><td>8</td><td>20</td><td>7</td></tr><tr><td>8192</td><td>53</td><td>6</td><td>15</td><td>26</td><td>61</td><td>9</td><td>25</td><td>5</td><td>65</td><td>8</td><td>20</td><td>7</td></tr></tbody></table>

## 6.2 Example 2: Adaptive algorithm for singular solution.

We select $f, A := I_{2 \times 2}$ and Dirichlet boundary conditions on the L-shaped domain $\Omega := (-1, 1)^2 \setminus \{(0, 1) \times (-1, 0)\}$, so that the exact solution is in polar coordinates

$$u = r^{\frac{2}{3}} \sin(2\theta/3). \qquad (53)$$

We test an h-adaptive algorithm driven by the a posteriori error estimator from Section 5. The adaptive algorithm starts from a coarse mesh and uses the estimator on the right-hand side of (50) to mark mesh cells for refinement through a bulk-chasing criterion (also known as Dörfler's marking). The adaptive algorithm can be classically described as

SOLVE → ESTIMATE → MARK → REFINE.

We first test the convergence rate of the above adaptive algorithm with $k \in \{0, 1, 2, 3\}$ and a bulk-chasing criterion set to 40%. The energy error and the a posteriori error estimator are reported in
<figure><img src="image_6.png" /><figcaption>Figure 3: Example 2. Energy error and a posteriori error estimator as a function of DoFs for $k \in \{0, 1, 2, 3\}$.</figcaption></figure>

Figure 3 with convergence rates computed as a function of DoFs. We observe that the energy error and the a posteriori error estimator converge at the optimal rate $O(DoFs^{-\frac{k+1}{2}})$. In contrast, both the energy error and the a posteriori error estimator converge at the suboptimal rate $O(DoFs^{-\frac{1}{3}})$ under uniform refinement, independently of $k$. Moreover, we observe in the left panel of Figure 4 that the effectivity index remains well behaved as a function of DoFs and that it slightly increases with the polynomial degree $k \ge 1$, with values between 2 and 2.8. For $k = 0$, the effectivity index is almost 3. To gain further insight, we report in the right panel of Figure 4 the effectivity index as a function of the polynomial degree $k \in \{1, \dots, 10\}$ on a mesh consisting of 96 triangular cells. We observe an algebraic rate of $p^{\frac{1}{2}}$, again in agreement with Theorem 5.7.

Finally, we report in Table 2 the relative contributions of the various terms composing the a posteriori error estimator (all in percentage) with polynomial degrees $k \in \{0, 1, 2\}$. For $k = 1$, the tangential jump estimator dominates the total estimator (by more than 58%), followed by $\eta_{\text{nor}}$ (25%), $\eta_{\text{sta}}$ (10%) and $\eta_{\text{res}}$ (0%). For $k = 2$, $\eta_{\text{tan}}$, $\eta_{\text{nor}}$ and $\eta_{\text{res}}$ are fairly well equilibrated (all at about 30%), whereas $\eta_{\text{sta}}$ is below 10%. For $k = 0$, only the tangential jump estimator is nonzero. Clearly, $\eta_{\text{res}} = 0$ since $f = 0$ and $u_K$ is affine for $k = 0$. Moreover, one can show that for $f = 0$ and $k = 0$, the Crouzeix-Raviart FEM solution
<figure><img src="image_7.png" /><figcaption>Figure 4: Example 2. Effectivity index as a function of DoFs for $k \in \{0, 1, 2, 3\}$ (left). Effectivity index as a function of $k \in \{1, \dots, 10\}$ on a mesh composed of 96 cells (right).</figcaption></figure>

is the cellwise component of the HHO solution, while the facewise component is the mean-value of the trace of the cell components. This explains why $\eta_{\text{sta}} = 0$ in Table 2, and consequently $\eta_{\text{nor}} = 0$ owing to the local conservation property (20).

<table><caption>Table 2: Example 2. Relative contribution (in %) of the various terms composing the a posteriori error estimator for $k \in \{0, 1, 2\}$.</caption><thead><tr><th colspan="5">k = 0</th><th colspan="5">k = 1</th><th colspan="5">k = 2</th></tr><tr><th># cell</th><th>η<sub>res</sub></th><th>η<sub>sta</sub></th><th>η<sub>nor</sub></th><th>η<sub>tan</sub></th><th># cell</th><th>η<sub>res</sub></th><th>η<sub>sta</sub></th><th>η<sub>nor</sub></th><th>η<sub>tan</sub></th><th># cell</th><th>η<sub>res</sub></th><th>η<sub>sta</sub></th><th>η<sub>nor</sub></th><th>η<sub>tan</sub></th></tr></thead><tbody><tr><td>116</td><td>0</td><td>0</td><td>0</td><td>100</td><td>118</td><td>0</td><td>10</td><td>32</td><td>58</td><td>172</td><td>28</td><td>9</td><td>29</td><td>34</td></tr><tr><td>1118</td><td>0</td><td>0</td><td>0</td><td>100</td><td>1207</td><td>0</td><td>10</td><td>26</td><td>64</td><td>1348</td><td>28</td><td>8</td><td>30</td><td>34</td></tr><tr><td>5948</td><td>0</td><td>0</td><td>0</td><td>100</td><td>6098</td><td>0</td><td>10</td><td>26</td><td>64</td><td>5856</td><td>27</td><td>8</td><td>32</td><td>33</td></tr><tr><td>22306</td><td>0</td><td>0</td><td>0</td><td>100</td><td>21762</td><td>0</td><td>10</td><td>25</td><td>65</td><td>21574</td><td>26</td><td>8</td><td>33</td><td>33</td></tr></tbody></table>

## 6.3 Example 3: Adaptive algorithm for Kellogg’s test case

Our last example is Kellogg’s test case [38], i.e., a diffusion problem on the square domain $\Omega := (-1, 1)^2$ with a checkerboard pattern for the diffusion coefficient $A$, namely $A := b$ for $xy \ge 0$ and $A := 1$ otherwise. The exact solution (for zero right-hand side and suitable Dirichlet boundary conditions) takes the form $u := r^\alpha \phi(\theta)$ in polar coordinates, where the explicit expression for the function $\phi$ can be found in [47, Section 5.3]. We select $\alpha = 0.1$ and $b = 161.4476387975881$, so that $u \in H^{1.1-\epsilon}(\Omega)$ with $\epsilon > 0$ arbitrarily small.

We first test the convergence rate of the h-adaptive algorithm described in the previous section with $k \in \{0, 1, 2, 3\}$ and a bulk-chasing criterion set to 10%. The energy error and the a posteriori error estimator are reported in Figure 5 with convergence rates reported as a function of DoFs. We observe that the energy error and the a posteriori error estimator converge at the optimal rate $O(\text{DoFs}^{-\frac{k+1}{2}})$. In contrast to Figure 3, we do not include here the convergence plots under uniform mesh refinement because the poor regularity of the exact solution leads to an extremely slow decay; for instance, the energy error is still above $0.1$ for $10^5$ DoFs, independently of $k$. Moreover, we observe in Figure 6 that the effectivity index remains well behaved as a function of DoFs and that it slightly increases with the polynomial degree $k \ge 1$, with values between 1.9 and 2.6. For $k = 0$, the effectivity index is almost 2.8.

Finally, we report in Table 3 the relative contributions of the various terms composing the a posteriori error estimator (all in percentage) with polynomial degrees $k \in \{0, 1, 2\}$. For $k = 1$, the tangential jump
<figure><img src="image_8.png" /><figcaption>Figure 5: Example 3. Energy error and a posteriori error estimator as a function of DoFs for $k \in \{0, 1, 2, 3\}$.</figcaption></figure>

estimator dominates the total estimator (by more than 59%), followed by $\eta_{nor}$ (27%), while $\eta_{sta}$ is below 10% and $\eta_{res}$ is zero. For $k = 2$, the tangential jump estimator dominates the total estimator on the coarsest mesh, whereas the contributions of $\eta_{tan}$, $\eta_{nor}$ and $\eta_{res}$ are fairly well equilibrated on the other meshes. For $k = 0$, only the tangential jump residual is nonzero, for the same reasons as discussed in the previous example.

<table><caption>Table 3: Example 3. Relative contribution (in %) of the various terms composing the a posteriori error estimator for $k \in \{0, 1, 2\}$.</caption><thead><tr><th colspan="5">k = 0</th><th colspan="5">k = 1</th><th colspan="5">k = 2</th></tr><tr><th># cell</th><th>η<sub>res</sub></th><th>η<sub>sta</sub></th><th>η<sub>nor</sub></th><th>η<sub>tan</sub></th><th># cell</th><th>η<sub>res</sub></th><th>η<sub>sta</sub></th><th>η<sub>nor</sub></th><th>η<sub>tan</sub></th><th># cell</th><th>η<sub>res</sub></th><th>η<sub>sta</sub></th><th>η<sub>nor</sub></th><th>η<sub>tan</sub></th></tr></thead><tbody><tr><td>122</td><td>0</td><td>0</td><td>0</td><td>100</td><td>124</td><td>0</td><td>7</td><td>8</td><td>85</td><td>142</td><td>22</td><td>9</td><td>11</td><td>58</td></tr><tr><td>1386</td><td>0</td><td>0</td><td>0</td><td>100</td><td>1376</td><td>0</td><td>8</td><td>33</td><td>59</td><td>1326</td><td>29</td><td>9</td><td>30</td><td>32</td></tr><tr><td>5962</td><td>0</td><td>0</td><td>0</td><td>100</td><td>5864</td><td>0</td><td>9</td><td>28</td><td>63</td><td>5486</td><td>28</td><td>8</td><td>32</td><td>32</td></tr><tr><td>15642</td><td>0</td><td>0</td><td>0</td><td>100</td><td>15342</td><td>0</td><td>10</td><td>27</td><td>63</td><td>14894</td><td>28</td><td>8</td><td>32</td><td>32</td></tr></tbody></table>
<figure><img src="image_9.png" /><figcaption>Figure 6: Example 3. Effectivity index as a function of DoFs with $k \in \{0, 1, 2, 3\}$.</figcaption></figure>

# 7 Proofs

This section collects the proofs of all the preparatory results stated in the previous section, namely Lemma 5.1, Lemma 5.2, Lemma 5.4, and Theorem 5.9.

## 7.1 Proof of Lemma 5.1

We will prove the following result:

$$ \|A^{\frac{1}{2}} \nabla e_c\|_{\Omega} \le C_{c,1} \left\{ \sum_{K \in \mathcal{T}_h} \left\{ A_K^{-1} \left( \frac{h_K}{k+1} \right)^2 \| \Pi_K^{k+1}(f) + A_K \Delta R_K^{k+1}(\hat{u}_K) \|_{K}^2 + A_K(k+1) S_{\partial K}(\hat{u}_K, \hat{u}_K) \right. \right. \\ \left. \left. + O_K(f)^2 + O_K(g_N)^2 \right\} \right\}^{\frac{1}{2}} + \|A^{\frac{1}{2}} \nabla_{\tau_h} e_d\|_{\Omega}. \quad (54) $$

The proof is split into two steps.

(i) Since $e_c \in H_{0,D}^1(\Omega)$, using (39) leads to

$$ \|A^{\frac{1}{2}} \nabla e_c\|_{\Omega}^2 = (A \nabla_{\tau_h} e, \nabla e_c)_{\Omega} - (A \nabla_{\tau_h} e_d, \nabla e_c)_{\Omega}. \quad (55) $$

For the second term, the Cauchy–Schwarz inequality gives

$$ |(A \nabla_{\tau_h} e_d, \nabla e_c)|_{\Omega} \le \|A^{\frac{1}{2}} \nabla e_c\|_{\Omega} \|A^{\frac{1}{2}} \nabla_{\tau_h} e_d\|_{\Omega}. \quad (56) $$

Next, we focus on bounding $(A \nabla_{\tau_h} e, \nabla e_c)_{\Omega}$. Using the weak form of the PDE (2), adding and subtracting the term $(A_K \nabla R_K^{k+1}(\hat{u}_K), \nabla e_c)_K$ for all $K \in \mathcal{T}_h$, and using the discrete problem (18) with some test function $\hat{w}_h \in \hat{V}_{h0}^k$, we infer that

$$ (A \nabla_{\tau_h} e, \nabla e_c)_{\Omega} = \sum_{K \in \mathcal{T}_h} \left\{ (f, e_c)_K + (g_N, e_c)_{\partial K^{\mathbb{N}}} - (A_K \nabla R_K^{k+1}(\hat{u}_K), \nabla e_c)_K \right. \\ \left. + (A_K \nabla (R_K^{k+1}(\hat{u}_K) - u_K), \nabla e_c)_K + (A_K \nabla R_K^{k+1}(\hat{u}_K), \nabla R_K^{k+1}(\hat{w}_K))_K \right. \\ \left. + S_{\partial K}(\hat{u}_K, \hat{w}_K) - (f, w_K)_K - (g_N, w_{\partial K})_{\partial K^{\mathbb{N}}} \right\}. $$

Let us set $\hat{w}_h := ((I_{mKM}^k(e_c)|_K)_{K \in \mathcal{T}_h}, (I_{mKM}^k(e_c)|_F)_{F \in \mathcal{F}_h})$. This definition is meaningful since $I_{mKM}^k(e_c)$ is single-valued at the mesh interfaces; moreover, since $I_{mKM}^k(e_c)$ vanishes on the boundary faces in $\mathcal{F}_h^D$,
we have $\hat{w}_h \in \hat{V}_{h0}^k$. We observe that, for all $K \in \mathcal{T}_h$, the definition (10) of the reconstruction operator gives $R_K^{k+1}(\hat{w}_K) = w_K$, and we also have $S_{\partial K}(\hat{w}_K, \hat{w}_K) = 0$. Hence, defining $\hat{\eta}_K = (\eta_K, \eta_{\partial K}) := (e_c|_K - w_K, e_c|_{\partial K} - w_{\partial K})$ for all $K \in \mathcal{T}_h$, we infer that

$$ (A \nabla_{\mathcal{T}_h} e, \nabla e_c)_{\Omega} = \sum_{K \in \mathcal{T}_h} \left\{ (f, \eta_K)_K + (g_N, \eta_{\partial K})_{\partial K^N} - (A_K \nabla R_K^{k+1}(\hat{u}_K), \nabla \eta_K)_K \right. \\ \left. + (A_K \nabla (R_K^{k+1}(\hat{u}_K) - u_K), \nabla e_c)_K \right\}. $$

Integrating by parts and using that $\eta_{\partial K^D} = 0$ and that $A_K$ is constant on $K$ gives

$$ (A \nabla_{\mathcal{T}_h} e, \nabla e_c)_{\Omega} = \sum_{K \in \mathcal{T}_h} \left\{ (f + A_K \Delta R_K^{k+1}(\hat{u}_K), \eta_K)_K - (A_K \nabla R_K^{k+1}(\hat{u}_K) \cdot n_K, \eta_{\partial K})_{\partial K^i} \right. \\ \left. + (g_N - A_K \nabla R_K^{k+1}(\hat{u}_K) \cdot n_K, \eta_{\partial K})_{\partial K^N} + (A_K \nabla (R_K^{k+1}(\hat{u}_K) - u_K), \nabla e_c)_K \right\}. \quad (57) $$

(ii) Since $\eta_{\partial K}|_F = (e_c - I_{mKM}(e_c))|_F$ is single-valued on each interface $F \in \mathcal{F}_{\partial K^i}$, using the local conservation property (20) of the HHO method gives

$$ (A \nabla_{\mathcal{T}_h} e, \nabla e_c)_{\Omega} = \sum_{K \in \mathcal{T}_h} \left\{ (f + A_K \Delta R_K^{k+1}(\hat{u}_K)), \eta_K)_K - A_K \frac{(k+1)^2}{h_K} (\Pi_{\partial K}^k(u_K|_{\partial K}) - u_{\partial K}, \eta_{\partial K})_{\partial K^i} \right. \\ \left. + (g_N - \Pi_{\partial K}^k(g_N|_{\partial K}), \eta_{\partial K})_{\partial K^N} - A_K \frac{(k+1)^2}{h_K} (\Pi_{\partial K}^k(u_K|_{\partial K}) - u_{\partial K}, \eta_{\partial K})_{\partial K^N} \right. \\ \left. + (A_K \nabla (R_K^{k+1}(\hat{u}_K) - u_K), \nabla e_c)_K \right\}. $$

Invoking the Cauchy-Schwarz inequality leads to

$$ \begin{aligned} |(A \nabla_{\mathcal{T}_h} e, \nabla e_c)_{\Omega}| & \le \sum_{K \in \mathcal{T}_h} \left\{ \left(A_K^{-\frac{1}{2}} \frac{h_K}{k+1}\right) \|f + A_K \Delta R_K^{k+1}(\hat{u}_K)\|_K \left(A_K^{\frac{1}{2}} \frac{k+1}{h_K}\right) \|\eta_K\|_K \right. \\ & \qquad + \left(A_K \frac{(k+1)^3}{h_K}\right)^{\frac{1}{2}} \|\Pi_{\partial K}^k(u_K|_{\partial K}) - u_{\partial K}\|_{\partial K^i} \left(A_K \frac{k+1}{h_K}\right)^{\frac{1}{2}} \|\eta_{\partial K}\|_{\partial K^i} \\ & \qquad + \left(A_K^{-1} \frac{h_K}{k+1}\right)^{\frac{1}{2}} \|g_N - \Pi_{\partial K}^k(g_N|_{\partial K})\|_{\partial K^N} \left(A_K \frac{k+1}{h_K}\right)^{\frac{1}{2}} \|\eta_{\partial K}\|_{\partial K^N} \\ & \qquad + \left(A_K \frac{(k+1)^3}{h_K}\right)^{\frac{1}{2}} \|\Pi_{\partial K}^k(u_K|_{\partial K}) - u_{\partial K}\|_{\partial K^N} \left(A_K \frac{k+1}{h_K}\right)^{\frac{1}{2}} \|\eta_{\partial K}\|_{\partial K^N} \\ & \qquad \left. + A_K^{\frac{1}{2}} \|\nabla(R_K^{k+1}(\hat{u}_K) - u_K)\|_K A_K^{\frac{1}{2}} \|\nabla e_c\|_K \right\}. \end{aligned} $$

Using the approximation result (7), we infer that, for all $K \in \mathcal{T}_h$,

$$ A_K^{\frac{1}{2}} \left\{ \left( \frac{k+1}{h_K} \right) \|\eta_K\|_K + \left( \frac{k+1}{h_K} \right)^{\frac{1}{2}} \|\eta_{\partial K}\|_{\partial K} \right\} \le C A_K^{\frac{1}{2}} \|\nabla e_c\|_{\text{es}(K)} \le C \chi_K(A) \|A^{\frac{1}{2}} \nabla e_c\|_{\text{es}(K)}. \quad (58) $$

Using the above bounds, the bound (14) on $\|\nabla(R_K^{k+1}(\hat{u}_K) - u_K)\|_K$, the mesh shape-regularity and the triangle inequality for the residual term gives

$$ \left| (A \nabla_{\mathcal{T}_h} e, \nabla e_c)_{\Omega} \right| \le C \max_{K \in \mathcal{T}_h} \chi_K(A) \left\{ \sum_{K \in \mathcal{T}_h} \left\{ A_K^{-1} \left( \frac{h_K}{k+1} \right)^2 \| \Pi_K^{k+1}(f) + A_K \Delta R_K^{k+1}(\hat{u}_K) \|_K^2 \right. \right. \\ \left. + A_K(k+1) S_{\partial K}(\hat{u}_K, \hat{u}_K) + A_K S_{\partial K}(\hat{u}_K, \hat{u}_K) + O_K(f)^2 + O_K(g_N)^2 \right\} \right\}^{\frac{1}{2}} \| A^{\frac{1}{2}} \nabla e_c \|_{\Omega}. $$

Combining (55) and (56) with the above bound completes the proof.
## 7.2 Proof of Lemma 5.2

We will prove the following result:

$$ \|A^{\frac{1}{2}} \nabla e_c\|_{\Omega} \le C_{c,2} \left\{ \sum_{K \in \mathcal{T}_h} \left\{ A_K^{-1} \left( \frac{h_K}{k+1} \right)^2 \| \Pi_K^{k+1}(f) + A_K \Delta R_K^{k+1}(\hat{u}_K) \|_K^2 + A_K S_{\partial K}(\hat{u}_K, \hat{u}_K) \right. \right. \\ \left. \left. + A_K^{-1} \left( \frac{h_K}{k+1} \right) \left( \| \llbracket A \nabla R_{\mathcal{T}_h}^{k+1}(\hat{u}_h) \rrbracket_{\partial K} \cdot \mathbf{n}_K \|_{\partial K^i}^2 + \| A_K \nabla R_K^{k+1}(\hat{u}_K) \cdot \mathbf{n}_K - \Pi_{\partial K}^k (g_N|_{\partial K}) \|_{\partial K^N}^2 \right) \right. \right. \\ \left. \left. + O_K(f)^2 + O_K(g_N)^2 \right) \right\}^{\frac{1}{2}} + \| A^{\frac{1}{2}} \nabla_{\mathcal{T}_h} e_d \|_{\Omega}. \tag{59} $$

The starting point is the bound (57) obtained at the end of the first step of the proof of Lemma 5.1. However, we no longer invoke the local conservation property (20), but consider the jump of $A \nabla R_{\mathcal{T}_h}^{k+1}(\hat{u}_h)$ across all the mesh interfaces and its value at all the Neumann boundary faces. Recalling that $\eta_{\partial K}|_F = (e_c - I_{mKM}^k(e_c))|_F$ is single-valued on every interface $F \in \mathcal{F}_{\partial K^i}$, we infer from (57) that

$$ (A \nabla_{\mathcal{T}_h} e, \nabla e_c)_{\Omega} = \sum_{K \in \mathcal{T}_h} \left\{ \begin{aligned} & (f + A_K \Delta R_K^{k+1}(\hat{u}_K), \eta_K)_K + (A_K \nabla(R_K^{k+1}(\hat{u}_K) - u_K), \nabla e_c)_K \\ & - \sum_{F \in \mathcal{F}_{\partial K^i}} \frac{1}{2} (\| \llbracket A \nabla R_{\mathcal{T}_h}^{k+1}(\hat{u}_h) \rrbracket_F \cdot \mathbf{n}_F, \eta_{\partial K})_F \\ & - (A_K \nabla R_K^{k+1}(\hat{u}_K) \cdot \mathbf{n}_K - \Pi_{\partial K}^k (g_N|_{\partial K}), \eta_{\partial K})_{\partial K^N} + (g_N - \Pi_{\partial K}^k (g_N|_{\partial K}), \eta_{\partial K})_{\partial K^N} \end{aligned} \right\} $$

Invoking the Cauchy-Schwarz inequality leads to

$$ \begin{aligned} |(A \nabla_{\mathcal{T}_h} e, \nabla e_c)_{\Omega}| \le & \sum_{K \in \mathcal{T}_h} \left\{ \begin{aligned} & (A_K^{-\frac{1}{2}} \frac{h_K}{k+1}) \|f + A_K \Delta R_K^{k+1}(\hat{u}_K)\|_K (A_K^{\frac{1}{2}} \frac{k+1}{h_K}) \|\eta_K\|_K \\ & + A_K^{\frac{1}{2}} \|\nabla(R_K^{k+1}(\hat{u}_K) - u_K)\|_K A_K^{\frac{1}{2}} \nabla e_c \|_K \\ & + \sum_{F \in \mathcal{F}_{\partial K^i}} \frac{1}{2} (A_K^{-1} \frac{h_K}{k+1})^{\frac{1}{2}} \| \llbracket A \nabla R_{\mathcal{T}_h}^{k+1}(\hat{u}_h) \rrbracket_F \cdot \mathbf{n}_F \|_F (A_K \frac{k+1}{h_K})^{\frac{1}{2}} \|\eta_{\partial K}\|_F \\ & + (A_K^{-1} \frac{h_K}{k+1})^{\frac{1}{2}} \| A_K \nabla R_K^{k+1}(\hat{u}_K) \cdot \mathbf{n}_\Omega - \Pi_{\partial K}^k (g_N|_{\partial K}) \|_{\partial K^N} (A_K \frac{k+1}{h_K})^{\frac{1}{2}} \|\eta_{\partial K}\|_{\partial K^N} \\ & + (A_K^{-1} \frac{h_K}{k+1})^{\frac{1}{2}} \| g_N - \Pi_{\partial K}^k (g_N|_{\partial K}) \|_{\partial K^N} (A_K \frac{k+1}{h_K})^{\frac{1}{2}} \|\eta_{\partial K}\|_{\partial K^N} \end{aligned} \right\}. \end{aligned} $$

Using the Cauchy-Schwarz inequality, the bound (14) on $\|\nabla(R_K^{k+1}(\hat{u}_K) - u_K)\|_K$, the approximation result (58) on $\eta$, and the triangle inequality for the residual term, we obtain

$$ |(A \nabla_{\mathcal{T}_h} e, \nabla e_c)_{\Omega}| \le C \max_{K \in \mathcal{T}_h} \chi_K(A) \left\{ \sum_{K \in \mathcal{T}_h} \left\{ A_K^{-1} \left( \frac{h_K}{k+1} \right)^2 \| \Pi_K^{k+1}(f) + A_K \Delta R_K^{k+1}(\hat{u}_K) \|_K^2 + A_K S_{\partial K}(\hat{u}_K, \hat{u}_K) \right. \right. \\ \left. + A_K^{-1} \left( \frac{h_K}{k+1} \right) \left( \| \llbracket A \nabla R_{\mathcal{T}_h}^{k+1}(\hat{u}_h) \rrbracket_{\partial K} \cdot \mathbf{n}_K \|_{\partial K^i}^2 + \| A_K \nabla R_K^{k+1}(\hat{u}_K) \cdot \mathbf{n}_K - \Pi_{\partial K}^k (g_N|_{\partial K}) \|_{\partial K^N}^2 \right) \right. \\ \left. + O_K(f)^2 + O_K(g_N)^2 \right\}^{\frac{1}{2}} \| A^{\frac{1}{2}} \nabla e_c \|_{\Omega}. \right\} $$

Combining (55) and (56) with the above bound completes the proof.
## 7.3 Proof of Lemma 5.4

We will prove the following result: For all $\boldsymbol{a} \in \mathcal{V}_h$,

$$ \|A^{\frac{1}{2}} \nabla_{\mathcal{T}_h} e_d^a\|_{\omega_a} \le C_d^a \left\{ \sum_{K \in \mathcal{T}_a} \left\{ A_K S_{\partial K}(\hat{u}_K, \hat{u}_K) + O_K(g_D)^2 \right\} + \sum_{F \in \mathcal{F}_a \cap \mathcal{F}_h^1} A_F^b \left( \frac{h_F}{k+1} \right) \|\| \nabla u_{\mathcal{T}_h} \|\|_F \times \mathbf{n}_F \|\|_F^2 \right. \\ \left. + \sum_{F \in \mathcal{F}_a \cap \mathcal{F}_h^D} A_F^b \left( \frac{h_F}{k+1} \right) \|\nabla (u_{\mathcal{T}_h} - \Pi_{\partial K}^{k+1}(g_D)) \times \mathbf{n}_{\Omega} \|\|_F^2 \right\}^{\frac{1}{2}}. \quad (60) $$

To fix the ideas, we assume that $d=3$; all the results also hold for $d=2$.

(1) First, we recall a result on the Helmholtz decomposition on simply connected domains [33, Theorem 3]. Let $\omega$ be a bounded, simply connected, Lipschitz domain and consider a partition of its boundary $\partial\omega$ into two disjoint Lipschitz parts $\gamma_D$ and $\gamma_N$. Then, for all $\boldsymbol{w} \in \mathcal{L}^2(\omega)$, there exist

$$ \xi \in H_{0,D}^1(\omega) := \{\zeta \in H^1(\omega) \mid \zeta|_{\gamma_D} = 0\}, \quad (61a) $$

$$ \phi \in \mathbf{H}_{0,N}^1(\omega) := \{\psi \in \mathbf{H}^1(\omega) \mid \psi|_{\gamma_N} = \mathbf{0}\}, \quad (61b) $$

such that

$$ A\boldsymbol{w} = A\nabla\xi + \mathbf{curl}\phi \quad \text{in } \omega. \quad (62) $$

Moreover, the following holds with a constant $C_\omega$ only depending on $\omega$ and its boundary partition:

$$ \|A^{\frac{1}{2}} \boldsymbol{w}\|_{\omega}^{2} = \|A^{\frac{1}{2}} \nabla \xi\|_{\omega}^{2} + \|A^{-\frac{1}{2}} \mathbf{curl} \phi\|_{\omega}^{2}, \quad (63a) $$

$$ \|\nabla \phi\|_{\omega} \le C_{\omega} \|\mathbf{curl} \phi\|_{\omega}. \quad (63b) $$

For $d=2$, one has $C_\omega = 1$. For $d=3$, the constant $C_\omega$ depends on $\omega$ and is more delicate to estimate. In particular, the constant $C_\omega$ grows with the number of holes in $\omega$. An explicit bound on $C_\omega$ is derived in [32] when $\omega$ is the union of overlapping star-shaped domains. This is the case, in particular, when $\omega$ is a vertex star, which is the application setting we consider in what follows.

(2) We apply the Helmholtz decomposition (62) to $\boldsymbol{w} := \nabla_{\mathcal{T}_h} e_d^a$ on the vertex star $\omega := \omega_a$ with the boundary partition such that $\gamma_D := \partial\omega_a \cap (\Omega \cup \Gamma_D)$ and $\gamma_N := \partial\omega_a \cap \Gamma_N$. This gives $\xi \in H_{0,D}^1(\omega_a)$ and $\phi \in \mathbf{H}_{0,N}^1(\omega_a)$ such that $A\nabla_{\mathcal{T}_h} e_d^a = A\nabla\xi + \mathbf{curl}\phi$ in $\omega_a$. Notice that $H_{0,D}^1(\omega_a)$ is indeed the functional space defined in (44) for $g=0$. Notice also that the stability constant $C_\omega$ in (63) only depends on the mesh shape-regularity, since, in particular, there is only a finite number of possible partitions of the boundary of $\omega_a$ and this number is bounded in terms of the mesh shape-regularity. Taking the $L^2(\omega_a)$-inner product with $\nabla_{\mathcal{T}_h} e_d^a$ and observing that $(A\nabla_{\mathcal{T}_h} e_d^a, \nabla\xi)_{\omega_a} = 0$ owing to the definition (43) of $u_c^a$, we infer that

$$ \|A^{\frac{1}{2}} \nabla_{\mathcal{T}_h} e_d^a\|_{\omega_a}^2 = (\nabla_{\mathcal{T}_h} e_d^a, \mathbf{curl} \phi)_{\omega_a}. $$

(3) We consider a modified hp-Karkulik-Melenk interpolation operator on the vertex star $\omega_a$ preserving the homogeneous boundary condition on $\gamma_N = \partial\omega_a \cap \Gamma_N$, say $I_{mKM}^{k,a}$ (recall that $k \ge 1$ by assumption). This operator is constructed as in the proof of Corollary 2.5 and leads to the following local approximation result: For all $\boldsymbol{a} \in \mathcal{V}_h$, all $\phi \in \mathbf{H}_{0,N}^1(\omega_a)$, and all $K \in \mathcal{T}_a$,

$$ \left(\frac{k}{h_K}\right)^2 \|\phi - I_{mKM}^{k,a}(\phi)\|_K^2 + \left(\frac{k}{h_K}\right) \|\phi - I_{mKM}^{k,a}(\phi)\|_{\partial K}^2 + \|\nabla I_{mKM}^{k,a}(\phi)\|_K^2 \le C_{mKM}^{\mathcal{V}_h} \|\nabla\phi\|_{\omega_a}^2, \quad (64) $$

where the constant $C_{mKM}^{\mathcal{V}_h}$ is bounded in terms of the mesh shape-regularity. Notice that the bound on $C_{mKM}^{\mathcal{V}_h}$ is uniform over the mesh vertices; this follows by observing that the Karkulik-Melenk interpolation operator is itself constructed vertex-wise and that the nodal averaging operator is uniformly stable over the mesh vertices. Notice also that $\mathrm{es}(K) = \omega_a$ since $I_{mKM}^{k,a}$ is locally constructed in the vertex star $\omega_a$. Considering the vector-valued version of the above operator (simply acting componentwise), we obtain

$$ \|A^{\frac{1}{2}} \nabla_{\mathcal{T}_h} e_d^a\|_{\omega_a}^2 = (\nabla_{\mathcal{T}_h} e_d^a, \mathbf{curl}(\phi - I_{mKM}^{k,a}(\phi)))_{\omega_a} + (\nabla_{\mathcal{T}_h} e_d^a, \mathbf{curl} I_{mKM}^{k,a}(\phi))_{\omega_a} =: T_1 + T_2, $$

and it remains to bound $T_1$ and $T_2$.
(4) Bound on $T_1$. Since $\mathrm{curl} \nabla u_c^a = 0$, $(\phi - I_{\mathrm{mKM}}^{k,a}(\phi))|_{\partial\omega_a \cap \Gamma_N} = 0$, and $u_c^a \in H_{g,D}^1(\omega_a)$ (see (44)), integrating by parts the curl operator gives

$$ (\nabla u_c^a, \mathrm{curl}(\phi - I_{\mathrm{mKM}}^{k,a}(\phi)))_{\omega_a} = (\phi - I_{\mathrm{mKM}}^{k,a}(\phi), \nabla(\psi_a g_D) \times n_\Omega)|_{\partial\omega_a \cap \Gamma_D}. $$

Similar arguments, together with the fact that $\phi - I_{\mathrm{mKM}}^{k,a}(\phi)$ is single-valued across every interface $F \in \mathcal{F}_a \cap \mathcal{F}_h^i$, give

$$ (\nabla_{\mathcal{T}_h}(\psi_a u_{\mathcal{T}_h}), \mathrm{curl}(\phi - I_{\mathrm{mKM}}^{k,a}(\phi)))_{\omega_a} = \sum_{F \in \mathcal{F}_a \cap \mathcal{F}_h^i} (\phi - I_{\mathrm{mKM}}^{k,a}(\phi), [[\nabla_{\mathcal{T}_h}(\psi_a u_{\mathcal{T}_h})]]_F \times n_F)_F + \sum_{F \in \mathcal{F}_a \cap \mathcal{F}_h^D} (\phi - I_{\mathrm{mKM}}^{k,a}(\phi), \nabla(\psi_a u_{\mathcal{T}_h}) \times n_\Omega)_F. $$

Putting the above two identities together and using that

$$ \nabla_{\mathcal{T}_h} e_d^a = \nabla u_c^a - \nabla_{\mathcal{T}_h}(\psi_a u_{\mathcal{T}_h}), \quad \nabla_{\mathcal{T}_h}(\psi_a u_{\mathcal{T}_h}) = \psi_a \nabla_{\mathcal{T}_h} u_{\mathcal{T}_h} + u_{\mathcal{T}_h} \nabla \psi_a, $$

and that $(\nabla \psi_a \times n_F)_F$ is single-valued across every interface $F \in \mathcal{F}_a \cap \mathcal{F}_h^i$, we infer that

$$ T_1 = - \sum_{F \in \mathcal{F}_a \cap \mathcal{F}_h^i} (\phi - I_{\mathrm{mKM}}^{k,a}(\phi), \psi_a [[\nabla u_{\mathcal{T}_h}]]_F \times n_F)_F - \sum_{F \in \mathcal{F}_a \cap \mathcal{F}_h^D} (\phi - I_{\mathrm{mKM}}^{k,a}(\phi), \psi_a \nabla(u_{\mathcal{T}_h} - g_D) \times n_\Omega)_F \\ - \sum_{F \in \mathcal{F}_a \cap \mathcal{F}_h^i} (\phi - I_{\mathrm{mKM}}^{k,a}(\phi), (\nabla \psi_a \times n_F)[[u_{\mathcal{T}_h}]]_F)_F - \sum_{F \in \mathcal{F}_a \cap \mathcal{F}_h^D} (\phi - I_{\mathrm{mKM}}^{k,a}(\phi), (\nabla \psi_a \times n_\Omega)(u_{\mathcal{T}_h} - g_D))_F. $$

Let $T_{11}$ and $T_{12}$ denote, respectively, the terms on the first and second lines on the above right-hand side. Using the Cauchy–Schwarz inequality and $\|\psi_a\|_{L^\infty(F)} = 1$, we infer that

$$ |T_{11}| \le \left\{ \sum_{F \in \mathcal{F}_a \cap \mathcal{F}_h^i} \left(\frac{h_F}{k+1}\right) \|[[\nabla u_{\mathcal{T}_h}]]_F \times n_F\|_F^2 + \sum_{F \in \mathcal{F}_a \cap \mathcal{F}_h^D} \left(\frac{h_F}{k+1}\right) \|\nabla(u_{\mathcal{T}_h} - g_D) \times n_\Omega\|_F^2 \right\}^{\frac{1}{2}} \\ \times \left\{ \sum_{F \in \mathcal{F}_a \cap (\mathcal{F}_h^i \cup \mathcal{F}_h^D)} \left(\frac{k+1}{h_F}\right) \|\phi - I_{\mathrm{mKM}}^{k,a}(\phi)\|_F^2 \right\}^{\frac{1}{2}}. $$

For every $F \in \mathcal{F}_a \cap (\mathcal{F}_h^i \cup \mathcal{F}_h^D)$, we can pick a mesh cell $K \in \mathcal{T}_a$ of which $F$ is a face and obtain from (64) that

$$ \left(\frac{k+1}{h_F}\right)^{\frac{1}{2}} \|\phi - I_{\mathrm{mKM}}^{k,a}(\phi)\|_F \le C \|\nabla \phi\|_{\omega_a}. $$

Putting the above two bounds together gives

$$ |T_{11}| \le C \left\{ \sum_{F \in \mathcal{F}_a \cap \mathcal{F}_h^i} \left(\frac{h_F}{k+1}\right) \|[[\nabla u_{\mathcal{T}_h}]]_F \times n_F\|_F^2 + \sum_{F \in \mathcal{F}_a \cap \mathcal{F}_h^D} \left(\frac{h_F}{k+1}\right) \|\nabla(u_{\mathcal{T}_h} - g_D) \times n_\Omega\|_F^2 \right\}^{\frac{1}{2}} \|\nabla \phi\|_{\omega_a}. $$

Owing to (63), we infer that

$$ \|\nabla \phi\|_{\omega_a} \le C_{\omega_a} \|\mathrm{curl} \phi\|_{\omega_a} \le C_{\omega_a} (A_a^\sharp)^{\frac{1}{2}} \|A^{-\frac{1}{2}} \mathrm{curl} \phi\|_{\omega_a} \le C_{\omega_a} (A_a^\sharp)^{\frac{1}{2}} \|A^{\frac{1}{2}} \nabla_{\mathcal{T}_h} e_d\|_{\omega_a}. \quad (65) $$

Combining the above bounds and recalling the definition (37) of $\chi_a(A)$ and the definition (38) of $A_F^b$, we obtain

$$ |T_{11}| \le C \chi_a(A)^{\frac{1}{2}} \left\{ \sum_{F \in \mathcal{F}_a \cap \mathcal{F}_h^i} A_F^b \left(\frac{h_F}{k+1}\right) \|[[\nabla u_{\mathcal{T}_h}]]_F \times n_F\|_F^2 + \sum_{F \in \mathcal{F}_a \cap \mathcal{F}_h^D} A_F^b \left(\frac{h_F}{k+1}\right) \|\nabla(u_{\mathcal{T}_h} - g_D) \times n_\Omega\|_F^2 \right\}^{\frac{1}{2}} \\ \times \|A^{\frac{1}{2}} \nabla_{\mathcal{T}_h} e_d\|_{\omega_a}. $$
Let us now bound $T_{12}$. Letting $\mathbb{I}$ be the identity operator, we write $T_{12} = T_{12a} + T_{12b}$ with

$$ T_{12a} = - \sum_{F \in \mathcal{F}_a \cap \mathcal{F}_h^i} (\phi - I_{mKM}^{k,a}(\phi), (\nabla \psi_a \times n_F)(\mathbb{I} - \Pi_F^k)[[u_{\tau_h}]]_F)_F \\ - \sum_{F \in \mathcal{F}_a \cap \mathcal{F}_h^D} (\phi - I_{mKM}^{k,a}(\phi), (\nabla \psi_a \times n_\Omega)(\mathbb{I} - \Pi_F^k)(u_{\tau_h} - g_D)|_F)_F $$

$$ T_{12b} = - \sum_{F \in \mathcal{F}_a \cap \mathcal{F}_h^i} (\phi - I_{mKM}^{k,a}(\phi), (\nabla \psi_a \times n_F) \Pi_F^k [[u_{\tau_h}]]_F)_F \\ - \sum_{F \in \mathcal{F}_a \cap \mathcal{F}_h^D} (\phi - I_{mKM}^{k,a}(\phi), (\nabla \psi_a \times n_\Omega) (\Pi_F^k (u_{\tau_h} - g_D)|_F)_F. $$

Using the Cauchy–Schwarz inequality and $\|\nabla \psi_a \times n_F\|_{L^\infty(F)} \le Ch_F^{-1}$ gives

$$ |T_{12a}| \le C \left\{ \sum_{F \in \mathcal{F}_a \cap \mathcal{F}_h^i} \left( \frac{h_F}{k+1} \right) h_F^{-2} \|(\mathbb{I} - \Pi_F^k)[[u_{\tau_h}]]_F\|_F^2 + \sum_{F \in \mathcal{F}_a \cap \mathcal{F}_h^D} \left( \frac{h_F}{k+1} \right) h_F^{-2} \|(\mathbb{I} - \Pi_F^k)(u_{\tau_h} - g_D)|_F\|_F^2 \right\}^{\frac{1}{2}} \\ \times \left\{ \sum_{F \in \mathcal{F}_a \cap (\mathcal{F}_h^i \cup \mathcal{F}_h^D)} \left( \frac{k+1}{h_F} \right) \|\phi - I_{mKM}^{k,a}(\phi)\|_F^2 \right\}^{\frac{1}{2}}. $$

Invoking the Poincaré inequalities $\|(\mathbb{I} - \Pi_F^k)[[u_{\tau_h}]]_F\|_F \le Ch_F \|[\nabla u_{\tau_h}]_F \times n_F\|_F$ for all $F \in \mathcal{F}_a \cap \mathcal{F}_h^i$ and $\|(\mathbb{I} - \Pi_F^k)(u_{\tau_h} - g_D)|_F\|_F \le Ch_F \|\nabla(u_{\tau_h} - g_D) \times n_\Omega\|_F$ for all $F \in \mathcal{F}_a \cap \mathcal{F}_h^D$, we infer that

$$ |T_{12a}| \le C \left\{ \sum_{F \in \mathcal{F}_a \cap \mathcal{F}_h^i} \left( \frac{h_F}{k+1} \right) \|[\nabla u_{\tau_h}]_F \times n_F\|_F^2 + \sum_{F \in \mathcal{F}_a \cap \mathcal{F}_h^D} \left( \frac{h_F}{k+1} \right) \|\nabla(u_{\tau_h} - g_D) \times n_\Omega\|_F^2 \right\}^{\frac{1}{2}} \\ \times \left\{ \sum_{F \in \mathcal{F}_a \cap (\mathcal{F}_h^i \cup \mathcal{F}_h^D)} \left( \frac{k+1}{h_F} \right) \|\phi - I_{mKM}^{k,a}(\phi)\|_F^2 \right\}^{\frac{1}{2}}. $$

Invoking the same arguments as above to bound the term involving $\phi$, we conclude that

$$ |T_{12a}| \le C \chi_a(A)^{\frac{1}{2}} \left\{ \sum_{F \in \mathcal{F}_a \cap \mathcal{F}_h^i} A_F^b \left( \frac{h_F}{k+1} \right) \|[\nabla u_{\tau_h}]_F \times n_F\|_F^2 + \sum_{F \in \mathcal{F}_a \cap \mathcal{F}_h^D} A_F^b \left( \frac{h_F}{k+1} \right) \|\nabla(u_{\tau_h} - g_D) \times n_\Omega\|_F^2 \right\}^{\frac{1}{2}} \\ \times \|A^{\frac{1}{2}} \nabla_{\tau_h} e_d\|_{\omega_a}. $$

Turning our attention to $T_{12b}$, we observe that $u_F$ is single-valued at every interface $F \in \mathcal{F}_a \cap \mathcal{F}_h^i$, and that $u_F = \Pi_F^k(g_D|_F)$ on every Dirichlet boundary face $F \in \mathcal{F}_a \cap \mathcal{F}_h^D$. Therefore, we have

$$ T_{12b} = - \sum_{F \in \mathcal{F}_a \cap \mathcal{F}_h^i} (\phi - I_{mKM}^{k,a}(\phi), (\nabla \psi_a \times n_F) \Pi_F^k [[u_{\tau_h} - u_F]]_F)_F \\ - \sum_{F \in \mathcal{F}_a \cap \mathcal{F}_h^D} (\phi - I_{mKM}^{k,a}(\phi), (\nabla \psi_a \times n_\Omega) (\Pi_F^k (u_{\tau_h}|_F) - u_F)_F. $$

Using the Cauchy–Schwarz inequality and $\|\nabla \psi_a \times n_F\|_{L^\infty(F)} \le Ch_F^{-1}$ gives

$$ |T_{12b}| \le C \left\{ \sum_{K \in \mathcal{T}_a} \sum_{F \in \mathcal{F}_{\partial K}} \left( \frac{h_F}{k+1} \right) h_F^{-2} \| \Pi_{\partial K}^k (u_K|_{\partial K}) - u_F \|_F^2 \right\}^{\frac{1}{2}} \\ \times \left\{ \sum_{F \in \mathcal{F}_a \cap (\mathcal{F}_h^i \cup \mathcal{F}_h^D)} \left( \frac{k+1}{h_F} \right) \| \phi - I_{mKM}^{k,a}(\phi) \|_F^2 \right\}^{\frac{1}{2}}, $$

where, in the first factor, we re-organized the sum over the mesh faces in $\mathcal{F}_a$ as a sum over the mesh cells in $\mathcal{T}_a$. Invoking the same arguments as above to bound the factor involving $\phi$ and recalling the
definition (12) of the stabilization bilinear form $S_{\partial K}$, we obtain

$$
\begin{aligned}
|T_{12b}| & \le C\chi_a(A)^{\frac{1}{2}}(k+1)^{-\frac{3}{2}} \left\{ \sum_{K \in \mathcal{T}_a} A_K S_{\partial K}(\hat{u}_K, \hat{u}_K) \right\}^{\frac{1}{2}} \|A^{\frac{1}{2}} \nabla_{\tau_h} e_d\|_{\omega_a} \\
& \le C\chi_a(A)^{\frac{1}{2}} \left\{ \sum_{K \in \mathcal{T}_a} A_K S_{\partial K}(\hat{u}_K, \hat{u}_K) \right\}^{\frac{1}{2}} \|A^{\frac{1}{2}} \nabla_{\tau_h} e_d\|_{\omega_a}.
\end{aligned}
$$

Notice that we dropped the (favorable) factor $(k+1)^{-3/2}$ since the bound on $T_2$ derived below does not involve this factor. Putting everything together, we conclude that

$$
|T_1| \le C\chi_a(A)^{\frac{1}{2}} \left\{ \sum_{F \in \mathcal{F}_a \cap \mathcal{F}_h^i} A_F^b \left(\frac{h_F}{k+1}\right) \|\llbracket \nabla u_{\tau_h} \rrbracket_F \times \mathbf{n}_F \|_F^2 + \sum_{F \in \mathcal{F}_a \cap \mathcal{F}_h^D} A_F^b \left(\frac{h_F}{k+1}\right) \|\nabla (u_{\tau_h} - g_D) \times \mathbf{n}_\Omega \|_F^2 \right. \\
\left. \qquad + \sum_{K \in \mathcal{T}_a} A_K S_{\partial K}(\hat{u}_K, \hat{u}_K) \right\}^{\frac{1}{2}} \times \|A^{\frac{1}{2}} \nabla_{\tau_h} e_d\|_{\omega_a}.
$$

(5) Bound on $T_2$. Since $I_{\text{mKM}}^{k,a}(\phi) \in H_{0,N}^1(\omega_a)$, $[\text{curl } I_{\text{mKM}}^{k,a}(\phi) \cdot \mathbf{n}]_F = 0$ across every interface $F \in \mathcal{F}_a \cap \mathcal{F}_h^i$, $\text{curl } I_{\text{mKM}}^{k,a}(\phi) \cdot \mathbf{n}_F = 0$ on every boundary face $F \in \partial\omega_a \cap \Gamma_N$, and $\nabla \cdot \text{curl } I_{\text{mKM}}^{k,a}(\phi) = 0$ on $\omega_a$, integrating by parts the gradient operator gives

$$ (\nabla u_c^a, \text{curl } I_{\text{mKM}}^{k,a}(\phi))_{\omega_a} = (\psi_a g_D, \text{curl } I_{\text{mKM}}^{k,a}(\phi) \cdot \mathbf{n}_\Omega)_{\partial\omega_a \cap \Gamma_D} = \sum_{K \in \mathcal{T}_a} (\psi_a g_D, \text{curl } I_{\text{mKM}}^{k,a}(\phi) \cdot \mathbf{n}_\Omega)_{\partial K^D}. $$

Similarly, we have

$$ (\nabla_{\tau_h} (\psi_a u_{\tau_h}), \text{curl } I_{\text{mKM}}^{k,a}(\phi))_{\omega_a} = \sum_{K \in \mathcal{T}_a} (\psi_a u_K, \text{curl } I_{\text{mKM}}^{k,a}(\phi) \cdot \mathbf{n}_K)_{\partial K}. $$

Recalling that $T_2 = (\nabla_{\tau_h} e_d^a, \text{curl } I_{\text{mKM}}^{k,a}(\phi))_{\omega_a}$ and $e_d^a = u_c^a - \psi_a u_{\tau_h}$, we infer that

$$
\begin{aligned}
T_2 &= \sum_{K \in \mathcal{T}_a} \left\{ (g_D - u_K, \psi_a \text{curl } I_{\text{mKM}}^{k,a}(\phi) \cdot \mathbf{n}_\Omega)_{\partial K^D} - (u_K, \psi_a \text{curl } I_{\text{mKM}}^{k,a}(\phi) \cdot \mathbf{n}_K)_{\partial K^i} \right\} \\
&= \sum_{K \in \mathcal{T}_a} \left\{ (\Pi_{\partial K}^k (g_D|_{\partial K} - u_K|_{\partial K}), \psi_a \text{curl } I_{\text{mKM}}^{k,a}(\phi) \cdot \mathbf{n}_\Omega)_{\partial K^D} - (\Pi_{\partial K}^k (u_K|_{\partial K}), \psi_a \text{curl } I_{\text{mKM}}^{k,a}(\phi) \cdot \mathbf{n}_K)_{\partial K^i} \right\}.
\end{aligned}
$$

Notice that we used here that $\psi_a \text{curl } I_{\text{mKM}}^{k,a}(\phi) \cdot \mathbf{n}_K$ is in $\mathbb{P}^k(\mathcal{F}_{\partial K})$ to introduce the $L^2$-orthogonal projection $\Pi_{\partial K}^k$ in both terms on the right-hand side. Since $u_{\partial K}$ is single-valued on every interface $F \in \mathcal{F}_{\partial K^i}$ and $u_{\partial K}|_F = \Pi_F^k(g_D|_F)$ on every boundary face $F \in \mathcal{F}_{\partial K^D}$, and since $\|\psi_a\|_{L^\infty(F)} = 1$, we obtain

$$ T_2 = \sum_{K \in \mathcal{T}_a} - (\Pi_{\partial K}^k (u_K|_{\partial K}) - u_{\partial K}, \psi_a \text{curl } I_{\text{mKM}}^{k,a}(\phi) \cdot \mathbf{n}_K)_{\partial K^i \cup \partial K^D}. $$

Invoking the Cauchy–Schwarz inequality, the discrete trace inverse inequality (3), and the definition of the stabilization bilinear form $S_{\partial K}$, we infer that

$$
\begin{aligned}
|T_2| & \le (A_a^b)^{-\frac{1}{2}} \left\{ \sum_{K \in \mathcal{T}_a} A_K \frac{(k+1)^2}{h_K} \| \Pi_{\partial K}^k (u_K|_{\partial K}) - u_{\partial K} \|_{\partial K^i \cup \partial K^D}^2 \right\}^{\frac{1}{2}} \left\{ \sum_{K \in \mathcal{T}_a} \frac{h_K}{(k+1)^2} \| \text{curl } I_{\text{mKM}}^{k,a}(\phi) \|_{\partial K^i \cup \partial K^D}^2 \right\}^{\frac{1}{2}} \\
& \le C(A_a^b)^{-\frac{1}{2}} \left\{ \sum_{K \in \mathcal{T}_a} A_K S_{\partial K}(\hat{u}_K, \hat{u}_K) \right\}^{\frac{1}{2}} \| \text{curl } I_{\text{mKM}}^{k,a}(\phi) \|_{\omega_a} \\
& \le C'(A_a^b)^{-\frac{1}{2}} \left\{ \sum_{K \in \mathcal{T}_a} A_K S_{\partial K}(\hat{u}_K, \hat{u}_K) \right\}^{\frac{1}{2}} \| \nabla \phi \|_{\omega_a},
\end{aligned}
$$
where the last bound follows from $\|\text{curl } I_{\text{mKM}}^{k,a}(\phi)\|_{\omega_a} \le (d-1)\|\nabla I_{\text{mKM}}^{k,a}(\phi)\|_{\omega_a}$ and the $H^1$-stability of $I_{\text{mKM}}^{k,a}$ (see (64)). Proceeding as above to bound $\|\nabla\phi\|_{\omega_a}$, we conclude that

$$|T_2| \le C\chi_a(A)^{\frac{1}{2}} \left\{ \sum_{K \in \mathcal{T}_a} A_K S_{\partial K}(\hat{u}_K, \hat{u}_K) \right\}^{\frac{1}{2}} \times \|A^{\frac{1}{2}} \nabla_{\tau_h} e_d\|_{\omega_a}.$$

(3) Putting the bounds on $T_1$ and $T_2$ together and invoking the triangle inequality to introduce the oscillation term on $g_D$ proves the assertion.

**Remark 7.1** ($k=0$). Recall that the assumption $k \ge 1$ is needed to invoke the (modified) Karkulik-Melenk interpolation operator. In the case $k=0$, the nonconforming error can be bounded by using a (piecewise affine) nodal-averaging operator.

## 7.4 Proof of Theorem 5.9

We will prove the following result: For all $K \in \mathcal{T}_h$,

$$\eta_{K,\text{res}} \le C_1(k+1) \left( \|A^{\frac{1}{2}} \nabla e\|_K + (A_K S_{\partial K}(\hat{u}_K, \hat{u}_K))^{\frac{1}{2}} + O_K(f) \right), \quad (66a)$$

$$\eta_{K,\text{nor}} \le C_1(k+1)^{\frac{1}{2}} \left\{ \sum_{K \in \omega_K} A_K S_{\partial K}(\hat{u}_K, \hat{u}_K) \right\}^{\frac{1}{2}}, \quad (66b)$$

$$\eta_{K,\text{tan}} \le C_1(k+1)^{\frac{3}{2}} \left\{ \sum_{K \in \omega_K} \|A^{\frac{1}{2}} \nabla e\|_K^2 \right\}^{\frac{1}{2}}. \quad (66c)$$

We consider the following bubble functions from [44]: For all $K \in \mathcal{T}_h$ and all $F \in \mathcal{F}_h$, we set

$$b_K(x) := h_K^{-1} \text{dist}(x, \partial K), \quad \forall x \in K, \qquad b_F(x) := h_F^{-1} \text{dist}(x, \partial F), \quad \forall x \in F. \quad (67)$$

We recall the following hp-inverse inequalities from [44, Theorem 2.5]: For all $v \in \mathbb{P}^{k+1}(K)$, all $K \in \mathcal{T}_h$, and all $F \in \mathcal{F}_{\partial K}$, we have

$$\|v\|_K \le C(k+1) \|b_K^{\frac{1}{2}} v\|_K, \quad (68a)$$

$$\|\nabla(vb_K)\|_K \le C\left(\frac{k+1}{h_K}\right) \|b_K^{\frac{1}{2}}v\|_K, \quad (68b)$$

$$\|v\|_F \le C(k+1) \|b_F^{\frac{1}{2}}v\|_F. \quad (68c)$$

We also recall the following hp-polynomial extension result from [44, Lemma 2.6] (with $\varepsilon := (k+1)^{-2}$ and $\alpha := 1$): For all $F \in \mathcal{F}_{\partial K}$ and all $K \in \mathcal{T}_h$, there exists an extension operator $E_{K,F}: H_0^1(F) \to H^1(K)$ such that, for all $v \in \mathbb{P}^k(F)$,

$$E_{K,F}(b_F v)|_F = b_F v, \quad E_{K,F}(b_F v)|_{\partial K \setminus F} = 0, \quad \|\nabla E_{K,F}(b_F v)\|_K \le C\left(\frac{(k+1)^2}{h_K}\right)^{\frac{1}{2}} \|b_F^{\frac{1}{2}}v\|_F. \quad (69)$$

(1) Proof of (66a). Set $v_K := \Pi_K^{k+1}(f) + A_K \Delta R_K^{k+1}(\hat{u}_K) \in \mathbb{P}^{k+1}(K)$, $w_K := b_K v_K$, observe that $w_K$ vanishes on $\partial K$, and let $w$ be the zero-extension $w_K$ to $\Omega$. Since $(f, w_K)_K - (A_K \nabla u, \nabla w_K)_K = (f, w)_K - (A_K \nabla u, \nabla w)_K = 0$, an integration by parts for $\Delta R_K^{k+1}(\hat{u}_K)$ gives

$$\begin{align*}
\|b_K^{\frac{1}{2}} v_K\|_K^2 &= \|b_K^{\frac{1}{2}} (\Pi_K^{k+1}(f) + A_K \Delta R_K^{k+1}(\hat{u}_K))\|_K^2 \\
&= (\Pi_K^{k+1}(f) + A_K \Delta R_K^{k+1}(\hat{u}_K), w_K)_K \\
&= ((\Pi_K^{k+1}(f) - f), w_K)_K + (A_K \nabla e, \nabla w_K)_K + (A_K \nabla (u_K - R_K^{k+1}(\hat{u}_K)), \nabla w_K)_K \\
&\le \|\Pi_K^{k+1}(f) - f\|_K \|b_K v_K\|_K + A_K^{\frac{3}{2}} (\|A^{\frac{1}{2}} \nabla e\|_K + A_K^{\frac{3}{2}} \|\nabla (u_K - R_K^{k+1}(\hat{u}_K))\|_K) \|\nabla (b_K v_K)\|_K \\
&\le \left( \|\Pi_K^{k+1}(f) - f\|_K + C A_K^{\frac{1}{2}} \left( \frac{k+1}{h_K} \right) \left( \|A^{\frac{1}{2}} \nabla e\|_K + (A_K S_{\partial K}(\hat{u}_K, \hat{u}_K))^{\frac{1}{2}} \right) \right) \|b_K^{\frac{1}{2}} v_K\|_K,
\end{align*}$$
where the last bound follows from $b_K \le 1$, (68b) and (14). Hence,

$$
\begin{aligned}
\|b_K^{\frac{1}{2}} v_K\|_K &\le \|\Pi_K^{k+1}(f) - f\|_K + C A_K^{\frac{1}{2}} \left(\frac{k+1}{h_K}\right) \left(\|A^{\frac{1}{2}} \nabla e\|_K + (A_K S_{\partial K}(\hat{u}_K, \hat{u}_K))^{\frac{1}{2}}\right) \\
&= A_K^{\frac{1}{2}} \left(\frac{k+1}{h_K}\right) \left(O_K(f) + C \left(\|A^{\frac{1}{2}} \nabla e\|_K + A_K^{\frac{1}{2}} S_{\partial K}(\hat{u}_K, \hat{u}_K)^{\frac{1}{2}}\right)\right),
\end{aligned}
$$

where we used the definition (36a) of $O_K(f)$. Invoking (68a), we infer that

$$
\eta_{K,\text{res}} = A_K^{-\frac{1}{2}} \left(\frac{h_K}{k+1}\right) \|v_K\|_K \le C(k+1) A_K^{-\frac{1}{2}} \left(\frac{h_K}{k+1}\right) \|b_K^{\frac{1}{2}} v_K\|_K.
$$

Combining the above two bounds proves (66a).

(2) Proof of (66b). For every interface $F \in \mathcal{F}_{\partial K^i}$, using the local conservation property (20) and the triangle inequality gives

$$
\| \llbracket A \nabla R_{\tau_h}^{k+1}(\hat{u}_h) \rrbracket_F \cdot \mathbf{n}_F \|_F^2 \le \sum_{\mathcal{K} \in \{K, K'\}} 2 \left( \frac{A_{\mathcal{K}}(k+1)^2}{h_{\mathcal{K}}} \right)^2 \| \Pi_{\partial \mathcal{K}}^k (u_{\mathcal{K}}|_{\partial \mathcal{K}}) - u_{\partial \mathcal{K}} \|_F^2,
$$

where $K'$ denotes the mesh cell sharing $F$ with $K$. Using the mesh shape-regularity, we infer that

$$
\begin{aligned}
A_K^{-1} \left(\frac{h_K}{k+1}\right) \| \llbracket A \nabla R_{\tau_h}^{k+1}(\hat{u}_h) \rrbracket_F \cdot \mathbf{n}_F \|_F^2 &\le \sum_{\mathcal{K} \in \{K, K'\}} 2 \left(\frac{A_{\mathcal{K}}}{A_K}\right) \left(\frac{(k+1)h_K}{h_{\mathcal{K}}}\right) \left(\frac{A_{\mathcal{K}}(k+1)^2}{h_{\mathcal{K}}}\right) \|\Pi_{\partial \mathcal{K}}^k (u_{\mathcal{K}}|_{\partial \mathcal{K}}) - u_{\partial \mathcal{K}}\|_F^2 \\
&\le C \chi_K'(A) (k+1) \sum_{\mathcal{K} \in \{K, K'\}} A_{\mathcal{K}} S_{\partial \mathcal{K}}(\hat{u}_{\mathcal{K}}, \hat{u}_{\mathcal{K}}).
\end{aligned}
$$

Moreover, for every boundary face $F \in \mathcal{F}_{\partial K^N}$, we have

$$
\begin{aligned}
A_K^{-1} \left( \frac{h_K}{k+1} \right) \| A \nabla R_K^{k+1} (\hat{u}_K) \cdot \mathbf{n}_F - \Pi_{\partial K}^k (g_N|_{\partial K}) \|_F^2 &= (k+1) A_K \frac{(k+1)^2}{h_K} \| \Pi_{\partial K}^k (u_K|_{\partial K}) - u_{\partial K} \|_F^2 \\
&\le (k+1) A_K S_{\partial K} (\hat{u}_K, \hat{u}_K).
\end{aligned}
$$

Summing over all the faces $F \in \mathcal{F}_{\partial K^i} \cup \mathcal{F}_{\partial K^N}$ completes the proof of (66b).

(3) Proof of (66c). For every interface $F \in \mathcal{F}_{\partial K^i}$, we set $v_F := \llbracket \nabla u_{\tau_h} \rrbracket_F \times \mathbf{n}_F$ and $w_{\mathcal{K},F} := E_{\mathcal{K},F}(b_F v_F)$ for all $\mathcal{K} \in \{K, K'\}$, where $K'$ denotes, as above, the mesh cell sharing $F$ with $K$. Since $\llbracket \nabla u_{\tau_h} \rrbracket_F \times \mathbf{n}_F = \sum_{\mathcal{K} \in \{K, K'\}} \nabla u_{\mathcal{K}} \times \mathbf{n}_{\mathcal{K}}$ and $w_{\mathcal{K},F}|_F = b_F v_F$ owing to (69), we infer that

$$
\|b_F^{\frac{1}{2}} v_F\|_F^2 = \sum_{\mathcal{K} \in \{K, K'\}} (w_{\mathcal{K},F}, \nabla u_{\mathcal{K}} \times \mathbf{n}_{\mathcal{K}})_F = \sum_{\mathcal{K} \in \{K, K'\}} (\text{curl } w_{\mathcal{K},F}, \nabla u_{\mathcal{K}})_\mathcal{K} = \sum_{\mathcal{K} \in \{K, K'\}} (\text{curl } w_{\mathcal{K},F}, \nabla (u_{\mathcal{K}} - u))_\mathcal{K},
$$

where the second equality follows by integration by parts and the third equality additionally uses that $\llbracket \nabla u \rrbracket_F \times \mathbf{n}_F = 0$. The Cauchy–Schwarz inequality gives

$$
\|b_F^{\frac{1}{2}} v_F\|_F^2 \le (A_F^b)^{-\frac{1}{2}} \sum_{\mathcal{K} \in \{K, K'\}} \|A^{\frac{1}{2}} \nabla e\|_{\mathcal{K}} \|\mathrm{curl} \, w_{\mathcal{K},F}\|_{\mathcal{K}}.
$$

Invoking (69) gives

$$
\| \mathbf{curl} w_{\mathcal{K},F} \|_{\mathcal{K}} \le C \| \nabla E_{\mathcal{K},F} (b_F v_F) \|_{\mathcal{K}} \le C (k+1) h_K^{-\frac{1}{2}} \| b_F^{\frac{1}{2}} v_F \|_F.
$$

Hence,

$$
\|b_F^{\frac{1}{2}} v_F\|_F \le C(A_F^b)^{-\frac{1}{2}} (k+1) h_K^{-\frac{1}{2}} \sum_{\mathcal{K} \in \{K, K'\}} \|A^{\frac{1}{2}} \nabla e\|_{\mathcal{K}}.
$$

Since $\|v_F\|_F \le C(k+1)\|b_F^{\frac{1}{2}}v_F\|_F$ owing to (68c), combining the above bounds gives

$$
\|v_F\|_F \le C(A_F^b)^{-\frac{1}{2}} (k+1)^2 h_K^{-\frac{1}{2}} \sum_{\mathcal{K} \in \{K, K'\}} \|A^{\frac{1}{2}} \nabla e\|_{\mathcal{K}}.
$$
This shows that

$$ (A_F^b)^{\frac{1}{2}} \left( \frac{h_K}{k+1} \right)^{\frac{1}{2}} \| [\nabla u_{\tau_h}] \|_F \times \mathbf{n}_F \|_F = (A_F^b)^{\frac{1}{2}} \left( \frac{h_K}{k+1} \right)^{\frac{1}{2}} \| v_F \|_F \le C(k+1)^{\frac{3}{2}} \sum_{\mathcal{K} \in \{K, K'\}} \| A^{\frac{1}{2}} \nabla e \|_{\mathcal{K}}. $$

A similar bound can be established for all $F \in \mathcal{F}_{\partial K^D}$. Summing over all $F \in \mathcal{F}_{\partial K^i} \cup \mathcal{F}_{\partial K^D}$ completes the proof of (66c).

## Acknowledgment

The use of the Cleps computing platform at INRIA Paris is gratefully acknowledged.

## References

[1] J. AGHILI, D. A. DI PIETRO, AND B. RUFFINI, An hp-hybrid high-order method for variable diffusion on general meshes, Comput. Methods Appl. Math., 17 (2017), pp. 359-376.

[2] M. AINSWORTH, A posteriori error estimation for discontinuous Galerkin finite element approximation, SIAM J. Numer. Anal., 45 (2007), pp. 1777-1798.

[3] I. BABUŠKA AND M. SURI, The optimal convergence rate of the p-version of the finite element method, SIAM J. Numer. Anal., 24 (1987), pp. 750-776.

[4] R. BECKER, P. HANSBO, AND M. G. LARSON, Energy norm a posteriori error estimation for discontinuous Galerkin methods, Comput. Methods Appl. Mech. Engrg., 192 (2003), pp. 723-733.

[5] F. BERTRAND, C. CARSTENSEN, B. GRÄSSLE, AND N. TRAN, Stabilization-free HHO a posteriori error control, Numer. Math., 154 (2023), pp. 369-408.

[6] E. BURMAN AND A. ERN, Continuous interior penalty hp-finite element methods for advection and advection-diffusion equations, Math. Comp., 76 (2007), pp. 1119-1140.

[7] A. CANGIANI, Z. DONG, AND E. H. GEORGOULIS, A posteriori error estimates for discontinuous Galerkin methods on polygonal and polyhedral meshes, SIAM J. Numer. Anal., 61 (2023), pp. 2352-2380.

[8] A. CANGIANI, Z. DONG, E. H. GEORGOULIS, AND P. HOUSTON, hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes, SpringerBriefs in Mathematics, Springer Nature, Cham, 2017.

[9] C. CARSTENSEN, S. BARTELS, AND S. JANSCHE, A posteriori error estimates for nonconforming finite element methods, Numer. Math., 92 (2002), pp. 233-256.

[10] C. CARSTENSEN, B. GRÄSSLE, AND N. TRAN, Adaptive hybrid high-order method for guaranteed lower eigenvalue bounds, Numer. Math., (2024), pp. 1-39.

[11] C. CARSTENSEN AND J. HU, A unifying theory of a posteriori error control for nonconforming finite element methods, Numer. Math., 107 (2007), pp. 473-502.

[12] C. CARSTENSEN AND R. VERFÜRTH, Edge residuals dominate a posteriori error estimates for low order finite element methods, SIAM J. Numer. Anal., 36 (1999), pp. 1571-1587.

[13] T. CHAUMONT-FRELET, A. ERN, S. LEMAIRE, AND F. VALENTIN, Bridging the multiscale hybrid-mixed and multiscale hybrid high-order methods, ESAIM Math. Model. Numer. Anal., 56 (2022), pp. 261-285.

[14] A. CHERNOV, Optimal convergence estimates for the trace of the polynomial $L^2$-projection operator on a simplex, Math. Comp., 81 (2012), pp. 765-787.

[15] M. CICUTTIN, A. ERN, AND N. PIGNET, Hybrid high-order methods. A primer with application to solid mechanics, SpringerBriefs in Mathematics, Springer Nature, Cham, 2021.
[16] B. COCKBURN, D. A. DI PIETRO, AND A. ERN, *Bridging the Hybrid High-Order and hybridizable discontinuous Galerkin methods*, ESAIM Math. Model. Numer. Anal., 50 (2016), pp. 635-650.

[17] P. DANIEL, A. ERN, I. SMEARS, AND M. VOHRALÍK, *An adaptive hp-refinement strategy with com- putable guaranteed bound on the error reduction factor*, Comput. Math. Appl., 76 (2018), pp. 967–983.

[18] E. DARI, R. DURAN, C. PADRA, AND V. VAMPA, *A posteriori error estimators for nonconforming finite element methods*, ESAIM Math. Model. Numer. Anal., 30 (1996), pp. 385-400.

[19] 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 (2015), pp. 1–21.

[20] 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 (2014), pp. 461–472.

[21] D. A. DI PIETRO AND R. SPECOGNA, *An a posteriori-driven adaptive mixed high-order method with application to electrostatics*, J. Comput. Phys., 326 (2016), pp. 35–55.

[22] Z. DONG AND A. ERN, *Hybrid high-order and weak Galerkin methods for the biharmonic problem*, SIAM J. Numer. Anal., 60 (2022), pp. 2626–2656.

[23] H. EGGER AND C. WALUGA, *hp analysis of a hybrid DG method for Stokes flow*, IMA J. Numer. Anal., 33 (2013), pp. 687–721.

[24] A. ERN AND J.-L. GUERMOND, *Finite element quasi-interpolation and best approximation*, ESAIM Math. Model. Numer. Anal. (M2AN), 51 (2017), pp. 1367-1385.

[25] ——, *Finite Elements I: Approximation and Interpolation*, vol. 72 of Texts in Applied Mathematics, Springer Nature, Cham, Switzerland, 2021.

[26] ——, *Finite Elements II: Galerkin Approximation, Elliptic and Mixed PDEs*, vol. 73 of Texts in Applied Mathematics, Springer Nature, Cham, Switzerland, 2021.

[27] A. ERN AND M. STEINS, *Convergence analysis for the wave equation discretized with hybrid methods in space (HHO, HDG and WG) and the leapfrog scheme in time*, J. Sci. Comput., 101 (2024), p. 7.

[28] A. ERN AND A. F. STEPHANSEN, *A posteriori energy-norm error estimates for advection-diffusion equations approximated by weighted interior penalty methods*, J. Comput. Math., 26 (2008), pp. 488–510.

[29] A. ERN AND M. VOHRALÍK, *Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations*, SIAM J. Numer. Anal., 53 (2015), pp. 1058–1081.

[30] A. ERN AND M. VOHRALÍK, *Stable broken H¹ and H(div) polynomial extensions for polynomial-degree-robust potential and flux reconstruction in three space dimensions*, Math. Comp., 89 (2020), pp. 551–594.

[31] E. H. GEORGOULIS, E. HALL, AND J. M. MELENK, *On the suboptimality of the p-version interior penalty discontinuous Galerkin method*, J. Sci. Comput., 42 (2010), pp. 54–67.

[32] J. GUZMÁN AND A. J. SALGADO, *Estimation of the continuity constants for Bogovskií and regularized Poincaré integral operators*, J. Math. Anal. Appl., 502 (2021), pp. Paper No. 125246, 36.

[33] R. HIPTMAIR AND C. PECHSTEIN, *Regular decompositions of vector fields - continuous, discrete, and structure-preserving*, Tech. Rep. 2019-18, Seminar for Applied Mathematics, ETH Zürich, Switzerland, 2019.

[34] P. HOUSTON, D. SCHÖTZAU, AND T. WIHLER, *Energy norm a posteriori error estimation of hp-adaptive discontinuous Galerkin methods for elliptic problems*, Math. Models Methods Appl. Sci., 17 (2007), pp. 33–62.
[35] P. HOUSTON, C. SCHWAB, AND E. SÜLI, *Discontinuous hp-finite element methods for advection-diffusion-reaction problems*, SIAM J. Numer. Anal., 39 (2002), pp. 2133–2163.

[36] O. A. KARAKASHIAN AND F. PASCAL, *A posteriori error estimates for a discontinuous Galerkin approximation of second-order elliptic problems*, SIAM J. Numer. Anal., 41 (2003), pp. 2374–2399.

[37] M. KARKULIK AND J. M. MELENK, *Local high-order regularization and applications to hp-methods*, Comput. Math. Appl., 70 (2015), pp. 1606–1639.

[38] B. R. KELLOGG, *On the Poisson equation with intersecting interfaces*, Appl. Anal., 4 (1974), pp. 101–129.

[39] P. L. LEDERER, C. LEHRENFELD, AND J. SCHÖBERL, *Hybrid discontinuous Galerkin methods with relaxed $H(\operatorname{div})$-conformity for incompressible flows. Part II*, ESAIM Math. Model. Numer. Anal., 53 (2019), pp. 503–522.

[40] C. LEHRENFELD AND J. SCHÖBERL, *High order exactly divergence-free hybrid discontinuous Galerkin methods for unsteady incompressible flows*, Comput. Methods Appl. Mech. Engrg., 307 (2016), pp. 339–361.

[41] S. LEMAIRE, *Bridging the hybrid high-order and virtual element methods*, IMA J. Numer. Anal., 41 (2021), pp. 549–593.

[42] M. W. LICHT, *Smoothed Projections and Mixed Boundary Conditions*, Math. Comp., 88 (2017), pp. 607–635.

[43] J. M. MELENK, *hp-interpolation of nonsmooth functions and an application to hp-a posteriori error estimation*, SIAM J. Numer. Anal., 43 (2005), pp. 127–155.

[44] J. M. MELENK AND B. I. WOHLMUTH, *On residual-based a posteriori error estimation in hp-FEM*, Adv. Comput. Math., 15 (2001), pp. 311–331.

[45] ——, *On residual-based a posteriori error estimation in hp-FEM*, vol. 15, 2001, pp. 311–331.

[46] J. M. MELENK AND T. WURZER, *On the stability of the boundary trace of the polynomial $L^2$-projection on triangles and tetrahedra*, Comput. Math. Appl., 67 (2014), pp. 944–965.

[47] P. MORIN, R. H. NOCHETTO, AND K. G. SIEBERT, *Convergence of adaptive finite element methods*, SIAM Rev., 44 (2002), pp. 631–658.

[48] B. STAMM AND T. P. WIHLER, *hp-optimal discontinuous Galerkin methods for linear elliptic problems*, Math. Comp., 79 (2010), pp. 2117–2133.

[49] T. WARBURTON AND J. S. HESTHAVEN, *On the constants in hp-finite element trace inverse inequalities*, Comput. Methods Appl. Mech. Engrg., 192 (2003), pp. 2765–2773.