Stabilization-free HHO a posteriori error control

Abstract

The known a posteriori error analysis of hybrid high-order methods treats the stabilization contribution as part of the error and as part of the error estimator for an efficient and reliable error control. This paper circumvents the stabilization contribution on simplicial meshes and arrives at a stabilization-free error analysis with an explicit residual-based a posteriori error estimator for adaptive mesh-refining as well as an equilibrium-based guaranteed upper error bound (GUB). Numerical evidence in a Poisson model problem supports that the GUB leads to realistic upper bounds for the displacement error in the piecewise energy norm. The adaptive mesh-refining algorithm associated to the explicit residual-based a posteriori error estimator recovers the optimal convergence rates in computational benchmarks.

1 Introduction

Hybrid high-order methods (HHO) were introduced in [26, 27] and are examined in the textbooks [25, 29] as a promising class of flexible nonconforming discretization methods for partial differential equations that involve a parameter-free stabilization term for the link between the volume and skeletal variables.

1.1 Known a posteriori error estimator

The a priori error analysis of HHO involves the stability terms in extended norms as part of the methodology and motivated a first explicit residual-based a posteriori error estimator in [25] with a reformulation of the stabilization in the upper bound. Let $s_h(u_h,u_h)$ denote the stabilization at the discrete solution $u_h\in V_h$ and let the (elliptic) reconstruction $Ru_h$ of $u_h$ denote a piecewise polynomial of degree at most $k+1$ that approximates $u\in H^1(\Omega )$, cf. (2) and Sect. 3 below for further details. Then a possible error term reads

$$\begin{aligned} \text {total error}^2:=\Vert \nabla _{\textrm{pw}}(u-Ru_h)\Vert _{L^2(\Omega )}^2 +s_h(u_h,u_h). \end{aligned}$$

(1)

It is disputable if $s_h(u_h,u_h)\ge 0$ is an error contribution, but if the total error includes $s_h(u_h,u_h)$ (or an equivalent form), then the error estimator may also include this term (or a computable equivalent) for a reliable and efficient a posteriori error control. Amongst the many skeletal schemes like (nonconforming) virtual elements, hybridized (weak) discontinuous Galerkin schemes et al., the HHO methodology has a clear and efficacious stabilization

$$\begin{aligned} s_h(v_h, w_h) :=\sum _{T\in \mathcal {T}}\sum _{F\in \mathcal {F}(T)} h_F^{-1}\langle S_{TF}v_h, S_{TF}w_h\rangle _{L^2(F)} \end{aligned}$$

(2)

with the abbreviation $S_{TF}v_h :=\Pi _{F,k} \left( v_{\mathcal {T}} + (1-\Pi _{T,k})R v_{h} \right) |_{T}-v_{\mathcal {F}}|_{F}$ for $v_h=(v_\mathcal {T}, v_\mathcal {F})\in V_h$ in terms of the $L^2$ projections $\Pi _{K,k}$ onto polynomials of degree at most k on a facet or simplex $K\in \mathcal {F}\cup \mathcal {T}$ of diameter $h_K = \textrm{diam}(K)$; cf. Sect. 1.4 for further details. The original residual-based estimator $\eta _{\textrm{HHO}}$ from the textbook [25] for the Poisson model problem $-\Delta u = f$ includes (2) and an interpolation ${\mathcal {A}} Ru_h\in V$ of $Ru_h$ by nodal averaging in

$$\begin{aligned} \eta _{\textrm{HHO}}^2&= \Vert h_{\mathcal {T}}(1 - \Pi _0)(f+{\Delta }_{\textrm{pw}} Ru_h)\Vert _{L^2(\Omega )}^2 + \Vert \nabla _{\textrm{pw}}(1-{\mathcal {A}})Ru_h\Vert _{L^2(\Omega )}^2\\&+ s_h(u_h, u_h). \end{aligned}$$

(Multiplicative constants are undisplayed in this introduction for simplicity.) The results from Theorem 4.3 and 4.7 in [25] show reliability and efficiency for the total error (1) and piecewise polynomial source terms $f\in P_{k+1}(\mathcal {T})$,

$$\begin{aligned} \text {total error}^2 \approx \eta _{\textrm{HHO}}^2. \end{aligned}$$

1.2 Stabilization-free a posteriori error control

There are objections against the double role of $s_h(u_h,u_h)$ on both sides of the efficiency and reliability estimate. First, the term $s_h(u,u_h)$ may dominate both sides of the error estimate. In other words, the total error might be equivalent to $s_h(u_h,u_h)$, but the quantity of interest may exclusively be

$$\begin{aligned} \text {error}^2:=\Vert \nabla _{\textrm{pw}}(u-Ru_h)\Vert _{L^2(\Omega )}^2. \end{aligned}$$

Second, since the stabilization (2) incorporates a negative power of the mesh-size, a reduction property for local refinements remains unclear but is inevitable in the proofs of optimal convergence of an adaptive algorithm [5, 16]. This paper, therefore, asks a different question about the control of the error without the stabilization term (2) in the upper bound and introduces two stabilization-free error estimators (multiplicative constants are undisplayed)

$$\begin{aligned} \eta _{\textrm{res}}^2=\ {}&\Vert h_\mathcal {T}( f+{\Delta }_{\textrm{pw}} Ru_h)\Vert _{L^2(\Omega )}^2 +\sum _{F\in {\mathcal {F}}} h_F\Vert [ \nabla _{\textrm{pw}}Ru_h ]_F \Vert _{L^2(F)}^2,\\ \eta _{\textrm{eq},p}^{2}=\ {}&\textrm{osc}_{k+p}^{2}(f, \mathcal {T}) + \Vert Q_p - \nabla _\textrm{pw}R u_h \Vert _{L^2(\Omega )}^2 + \Vert \nabla _{\textrm{pw}}(1-{\mathcal {A}})Ru_h\Vert _{L^2(\Omega )}^2 \end{aligned}$$

for some parameter $p \in {\mathbb {N}}_0$. The explicit residual-based a posteriori error estimator $\eta _{\textrm{res}}$ follows from the a posteriori methodology in the spirit of [14, 18, 20, 21] with a piecewise volume residual $f+\Delta _\textrm{pw}Ru_h$ and the jumps $[\nabla _{\textrm{pw}}Ru_h]_F$ across a facet F (on the boundary this is only the tangential component of $\nabla _{\textrm{pw}}Ru_h$). The equilibrated error estimator $\eta _{\textrm{eq},p}$ includes the post-processed quantity $Q_p \in RT_{k+p}(\mathcal {T})$ in the space $RT_{k+p}(\mathcal {T})$ of Raviart-Thomas functions of degree $k+p$ for $p \in {\mathbb {N}}_0$ and the nodal average ${\mathcal {A}} R u_h \in S^{k+1}_0({\mathcal {T}})$ of $R u_h$. The main results establish reliability and efficiency

$$\begin{aligned} \text {error}^2 + \textrm{osc}_{k-1}^2(f, \mathcal {T}) \lesssim \eta _{\textrm{res}}^2 \approx \eta _{\textrm{eq},p}^2\lesssim \text {error}^2 + \textrm{osc}_{q}^2(f, \mathcal {T}) \end{aligned}$$

for any $p, q\in {\mathbb {N}}_0$ up to data-oscillations $\textrm{osc}_{q}^2(f, \mathcal {T}):=\Vert h_\mathcal {T}(1-\Pi _{q})f\Vert ^2_{L^2(\Omega )}$ of the source term $f\in L^2(\Omega )$ and without any stabilization terms. Computational benchmarks with adaptive mesh-refinement driven by any of these estimators provide numerical evidence for optimal convergence rates.

1.3 Further contributions and outline

The higher-order Crouzeix-Raviart finite element schemes are complicated at least in 3D [23] and then the HHO methodology is an attractive alternative even for simplicial triangulations with partly unexpected advantages like the computation of higher-order guaranteed eigenvalue bounds [15]. Higher convergence rates rely on an appropriate adaptive mesh-refining algorithm and hence stabilization-free a posteriori error estimators are of particular interest. The recent paper [36] establishes the latter for virtual elements with an over-penalization strategy as an extension of [8] for the discontinuous Galerkin schemes. A disadvantage is the quantification of the restriction on the stabilization parameter in practise and poor condition for larger parameters. The stabilization-free a posteriori error control in this paper is based on two observations for the HHO schemes on simplicial triangulations. First, the $P_1$-conforming finite element functions let the stabilization vanish and, second, the divergence-free lowest-order Raviart-Thomas functions are $L^2$ perpendicular to the piecewise gradients $\nabla _{\textrm{pw}}R u_h$. In fact, those two fairly general properties lead in Sect. 2 to a reliable explicit residual-based a posteriori error estimator. In contrast to the simplified introduction above, the paper also focuses on multiplicative constants that lead to the GUB

$$\begin{aligned} \text {error}\le \eta _{\textrm{res}}\quad \text {and error}\le \eta _{\textrm{eq},p}; \end{aligned}$$

cf. Table 1 for explicit quantities and Theorem 1 and Theorem 3 for further details.

Table 1 Explicit constants $C_1, C_\textrm{st}, C_2$ for right-isosceles triangles with respect to the maximum interior angle $\omega _{\textrm{max}}$ of the polygonal domain $\Omega$

Numerical comparisons of $\eta _{\textrm{HHO}}$ with $\eta _{\textrm{res}}$ and $\eta _{\textrm{eq},p}$ favour the latter. Section 2 identifies general building blocks of the a posteriori error analysis for discontinuous schemes with emphasis on explicit constants. An application to HHO leads to the new stabilization-free residual-based estimator $\eta _{\textrm{res}}$ in Sect. 3. The alternative stabilization-free error estimator $\eta _{\textrm{eq},p}$ follows from an equilibration strategy plus post-processing in Sect. 4. This paper also contributes to the HHO literature a local equivalence of two stabilizations and the efficiency of the stabilization terms up to data-oscillations in extension of [32]. Numerical comparisons of the different error estimators and an error estimator competition for guaranteed error control of the piecewise energy norm in 2D conclude this paper in Sect. 5. Three computational benchmarks provide striking numerical evidences for the optimality of the associated adaptive algorithms. The appendix provides algorithmic details on the computation of the post-processed contribution $\Vert Q_p-\nabla _{\textrm{pw}}Ru_h\Vert _{L^2(\Omega )}$ in $\eta _{\textrm{eq}, p}$.

1.4 Overall notation

Standard notation for Sobolev and Lebesgue spaces and norms apply with $\Vert \bullet \Vert :=\Vert \bullet \Vert _{L^2(\Omega )}$ and $|\!|\!|\bullet |\!|\!|:=\Vert \nabla \bullet \Vert _{L^2(\Omega )}$. In particular, $H(\mathop {\textrm{div}}\limits , \Omega )$ is the space of Sobolev functions with weak divergence in $L^2(\Omega )$ and $H(\mathop {\textrm{div}}\limits =0, \Omega )$ contains only divergence-free functions in $H(\mathop {\textrm{div}}\limits , \Omega )$. Throughout this paper, $\mathcal {T}$ denotes a shape-regular triangulation of the polyhedral bounded Lipschitz domain $\Omega \subset \mathbb R^n$ into n-simplices with facets $\mathcal {F}$ (edges for $n=2$ and faces for $n=3$) and vertices ${\mathcal {V}}$. Let $\mathcal {F}(\Omega )$ (resp. ${\mathcal {V}}(\Omega )$) denote the set of interior facets (resp. vertices) and (resp. ). Given $v \in \Omega \rightarrow {\mathbb {R}}^n$ and $w: \Omega \rightarrow {\mathbb {R}}^{2n-3}$, let if $n = 2$ and and $\mathop {\textrm{Curl}}\limits w :=\mathop {\textrm{curl}}\limits w$ if $n = 3$. For $s\in \mathbb R$, let $H^s(\mathcal {T}), H(\mathop {\textrm{div}}\limits ,\mathcal {T})$, and $H(\mathop {\textrm{curl}}\limits ,\mathcal {T})$ denote the space of piecewise Sobolev functions with restriction to $T \in \mathcal {T}$ in $H^s(T), H(\mathop {\textrm{div}}\limits ,T)$, and $H(\mathop {\textrm{curl}}\limits ,T)$. To simplify notation, $H^s(K)$ abbreviates $H^s(\textrm{int}(K))$ for the open interior $\textrm{int}(K)$ of a compact set K. The $L^2$-scalar product reads $(\bullet , \bullet )_{L^{2}(\omega )}$ for volumes $\omega \subseteq \Omega$ and $\langle \bullet , \bullet \rangle _{L^2(\gamma )}$ for surfaces $\gamma \subset {\overline{\Omega }}$ of co-dimension one; the same symbol applies to scalars and to vectors. For and , let $\langle \bullet , \bullet \rangle$ denote the duality-brackets in $V^* \times V$ for the dual space $V^*$ of V equipped with the operator norm .

Define the energy scalar product $a(v,w) :=(\nabla v, \nabla w)_{L^2(\Omega )}$ for $v, w \in H^{1}(\Omega )$ and its piecewise version $a_{\textrm{pw}}(v, w)=(\nabla _{\textrm{pw}}v,\nabla _{\textrm{pw}}w)_{L^2(\Omega )}$ for $v, w\in H^1(\mathcal {T})$. The latter induces the seminorm $|\!|\!|\bullet |\!|\!|_\textrm{pw}:=a_\textrm{pw}(\bullet ,\bullet )^{1/2}$ in $H^1({\mathcal {T}})$. Here and throughout the paper, $\nabla _{\textrm{pw}}, \textrm{div}_{\textrm{pw}}, \textrm{curl}_{\textrm{pw}}, {\Delta }_{\textrm{pw}}$, denote the piecewise evaluation of the differential operators $\nabla$, div, $\textrm{curl}_{\textrm{pw}}, \Delta$ without explicit reference to the underlying shape-regular triangulation $\mathcal {T}$.

The vector space $P_k(K)$ of polynomials of degrees at most $k\in {\mathbb N}_0$ over a facet or simplex $K\in \mathcal {F}\cup \mathcal {T}$ defines the piecewise polynomial spaces

$$\begin{aligned} P_k(\mathcal {T})&:=\{ p \in L^2(\Omega ) \ :\ p_{|T} \in P_{k}(T) \text { for all } T \in \mathcal {T}\}, \\ P_k(\mathcal {F})&:=\{ p \in L^2(\mathcal F) \ :\ p_{|F} \in P_k(F) \text { for all } F \in \mathcal {F}\} \ \end{aligned}$$

and the space of piecewise Raviart-Thomas functions

$$\begin{aligned} RT_{k}^\textrm{pw}(\mathcal {T}) :=P_k(\mathcal {T};\mathbb R^n) + xP_k(\mathcal {T}). \end{aligned}$$

The associated $L^2$ projections read $\Pi _{K,k}:L^2(\Omega )\rightarrow P_k(K), \Pi _k: L^2(\Omega ) \rightarrow P_k(\mathcal {T}),$ and $\Pi _{\mathcal F, k}: L^2(\Omega ) \rightarrow P_k(\mathcal {F})$ with the convention $\Pi _{-1}:=0$. Abbreviate $S^{k+1}_0(\mathcal {T}) :=P_{k+1}(\mathcal {T}) \cap V$ and $RT_k(\mathcal {T}) :=RT^\textrm{pw}_k(\mathcal {T}) \cap H(\mathop {\textrm{div}}\limits ,\Omega )$ for all $k \in {\mathbb {N}}_0$. The piecewise constant mesh-size function $h_\mathcal {T}\in P_0(\mathcal {T})$ satisfies $h_{\mathcal {T}|T}:=h_T$ for $T\in \mathcal {T}$ with the diameter $h_K :=\textrm{diam}(K)\in P_0(K)$ of $K\in \mathcal {F}\cup \mathcal {T}$.

If not explicitly stated otherwise, constants are independent of the mesh-size in the triangulation but may depend on the shape-regularity and on the polynomial degree k. The abbreviation $A\lesssim B$ hides a generic constant C (independent of the mesh-size) in $A\le C\; B; A\approx B$ abbreviates $A\lesssim B\lesssim A$.

2 Foundations of the a posteriori error analysis

This section investigates general building blocks of the a posteriori error analysis and revisits arguments from [14, 18, 20, 21] with emphasis on multiply connected domains $\Omega \subset \mathbb R^n$ for $n=2,3$. The general setting of this section results in reliability for an error estimator that is applicable beyond the HHO methodology. Consider the weak solution $u\in V = H^1_0(\Omega )$ to the Poisson model problem $-\Delta u=f$ a.e. in $\Omega$ and $u = 0$ on $\partial \Omega$ for a given source $f\in L^2(\Omega )$; i.e., $u\in V$ satisfies

$$\begin{aligned} a(u,v)= (f,v)_{L^2(\Omega )} \quad \text {for all } v\in V . \end{aligned}$$

(3)

An approximation of the gradient gives rise to the residual seen as a linear functional on V, i.e.,

Let $\nu _T$ denote the unit outer normal along the boundary $\partial T$ of each simplex $T\in \mathcal {T}$ and fix the orientation of a unit normal $\nu _F=\pm \nu _T$ for each facet $F\in \mathcal {F}(T)$ of T such that it matches the outer unit normal $\nu$ of $\partial \Omega$ at the boundary. The jump $[{G}]_F$ of a piecewise function in $m\in {\mathbb N}$ components ${G}\in H^1(\mathcal {T};\mathbb R^m)$ reads $[{G}]_F:={G}_{|T_+}-{G}_{|T_-}$ on interior facets $F=T_+\cap T_-\in \mathcal {F}(\Omega )$ (with $T_\pm$ labelled such that $\nu _{T_+|F} = \nu _F = -\nu _{T_-|F})$ and $[{G}]_F:={G}$ on the boundary $F\in \mathcal {F}(\partial \Omega )$. The main result of this section establishes the residual-based error estimator

$$\begin{aligned} \eta ^2(\mathcal {T}, {G}) :=&\left( C_{1}\Vert h_\mathcal {T}( f+{\mathop {\textrm{div}}\limits }_{\textrm{pw}} {G})\Vert +C_{2}\sqrt{\sum _{F\in {\mathcal {F}}(\Omega )} \ell (F)\Vert [ {G}]_F\cdot \nu _F \Vert _{L^2(F)}^2}\right) ^2 \nonumber \\&+C_{\textrm{H}}^2\left( C_{1}\Vert h_\mathcal {T}{\mathop {\textrm{curl}}\limits }_{\textrm{pw}} {G}\Vert +C_{2}\sqrt{\sum _{F\in {\mathcal {F}}} \ell (F)\Vert [ {G}]_F\times \nu _F \Vert _{L^2(F)}^2}\right) ^2 \end{aligned}$$

(4)

as a GUB $\Vert \nabla u - {G}\Vert \le \eta (\mathcal {T}, {G})$ under minimal assumptions on the approximation ${G}\in H^1(\mathcal {T};\mathbb R^n)\subset L^2(\Omega ; \mathbb R^n)$. The constants $C_1, C_2$, and $C_H$ (or upper bounds thereof) are computable; cf. Table 1 for an example in 2D with details in Example 1 at the end of Sect. 2. The first assumption is a weakened discrete solution property

$$\begin{aligned} ({G}, \nabla w_C)_{L^2(\Omega )}=(f,w_C)_{L^2(\Omega )}\quad \text {for all }w_C\in S_0^1(\mathcal T). \end{aligned}$$

(5)

The second assumption is the orthogonality to the lowest-order divergence-free Raviart-Thomas functions

$$\begin{aligned} ({G}, r)_{L^2(\Omega )}=0\quad \text {for all } r\in RT_0(\mathcal {T})\cap H(\mathop {\textrm{div}}\limits =0, \Omega ). \end{aligned}$$

(6)

Theorem 1

(residual-based GUB) Suppose that ${G}\in H^1(\mathcal T;\mathbb R^n)$ and $f\in L^2(\Omega )$ satisfy (5)–(6). Then the error estimator $\eta (\mathcal {T}, {G})$ from (4) is a GUB

$$\begin{aligned} \Vert \nabla u-{G}\Vert \le \eta (\mathcal {T}, {G}) \end{aligned}$$

of the error $\Vert \nabla u-{G}\Vert$ for the solution $u\in V$ to (3). The constants $C_\textrm{1}, C_{2}, C_{\textrm{H}}$ exclusively depend on $\Omega$ and the shape-regularity of $\mathcal T$.

The remaining parts of this section are devoted to the proof of Theorem 1 and the computation of (upper bounds of) the constants $C_\textrm{1}, C_\textrm{2}$, and $C_{\textrm{H}}$ in (4). The point of departure is the subsequent decomposition that appears necessary in the nonconforming and mixed finite element a posteriori error analysis. It leads to a split of the error $\Vert \nabla u-{G}\Vert$ into some divergence part and some consistency part.

Lemma 1

(decomposition) Any $v\in V$ and ${G}\in L^2(\Omega ; \mathbb R^n)$ satisfy the decomposition

(7)

with the (unique) minimizer $w\in V$ of the distance

of ${G}$ to the gradients $\nabla V$ of Sobolev functions. The solution $u\in V$ to (3) satisfies

$$\begin{aligned} \mu :=&\;|\!|\!| f+\mathop {\textrm{div}}\limits {G}|\!|\!|_{*}=|\!|\!|u-w|\!|\!|\qquad \text {and} \nonumber \\ \Vert \nabla u-{G}\Vert ^2 =&\;|\!|\!| f+\mathop {\textrm{div}}\limits {G}|\!|\!|_{*}^2+\Vert {G}- \nabla w\Vert ^2=\mu ^2+\delta ^2. \end{aligned}$$

(8)

Proof

The minimizer $w \in V$ of $\Vert {G}- \nabla \varphi \Vert$ among $\varphi \in V$ satisfies the variational formulation $a(w, \varphi )_{L^2(\Omega )} = ({G},\nabla \varphi )_{L^2(\Omega )}$ for all $\varphi \in V$. (Notice that w is the unique weak solution to the Poisson model problem $-\Delta w = -\mathop {\textrm{div}}\limits {G}\in V^*$.) In particular, ${G}- \nabla w$ is $L^2$ orthogonal onto $\nabla V$ and the Pythagoras theorem proves (7). Given $\varphi \in V$ with $|\!|\!|\varphi |\!|\!| =1$, the orthogonality of ${G}-\nabla w$ to $\nabla \varphi$ and (3) show

$$\begin{aligned} a(u-w,\varphi )_{L^2(\Omega )}= (\nabla u- {G},\nabla \varphi )_{L^2(\Omega )}=\langle f+\mathop {\textrm{div}}\limits {G},\varphi \rangle \end{aligned}$$

(9)

with the duality brackets $\langle \bullet , \bullet \rangle$ in $V^* \times V$. Since the supremum of (9) over all $\varphi \in V$ with $|\!|\!|\varphi |\!|\!| = 1$ is equal to $|\!|\!|u-w|\!|\!|= |\!|\!| f+\mathop {\textrm{div}}\limits {G}|\!|\!|_{*}$, this and (7) conclude the proof of (8). $\square$

The split (7) of the error $\Vert \nabla u-{G}\Vert$ allows for and enforces a separate estimation of the equilibrium and consistency contribution in residual-based a posteriori error estimators.

In order to derive explicit constants, two lemmas are recalled. The first has a long tradition in the a posteriori error control in form of a Helmholtz decomposition on simply connected domains [4, 13] and introduces the constant $C_{\textrm{H}}$ from Theorem 1. The following version includes the general case of multiply connected domains as in [33] for $n=2$ or $n=3$ dimensions and weak assumptions on a divergence-free function $\varrho \in H(\mathop {\textrm{div}}\limits =0,\Omega )$.

Lemma 2

(Helmholtz-decomposition) Suppose the divergence-free function $\varrho \in H(\mathop {\textrm{div}}\limits =0,\Omega )$ is $L^2$ orthogonal onto $RT_0(\mathcal T)\cap H(\mathop {\textrm{div}}\limits =0, \Omega )$. Then there exists $\beta \in H^1(\Omega ;\mathbb R^N), N=2n-3$, such that any $\beta _C\in S^1(\mathcal T)^N$ satisfies

$$\begin{aligned} \Vert \varrho \Vert ^2=\int _\Omega \varrho \cdot \textrm{Curl}(\beta -\beta _C)\ \textrm{d} x \quad \text {and}\quad |\!|\!|\beta |\!|\!|\le C_{\textrm{H}}\, \Vert \varrho \Vert . \end{aligned}$$

(10)

The constant $C_{\textrm{H}}>0$ exclusively depends on $\Omega$.

Proof

The compact polyhedral boundary $\partial \Omega$ of the bounded Lipschitz domain $\Omega$ has $J+1$ connectivity components $\Gamma _0,\dots ,\Gamma _J$ for some finite $J\in \mathbb N_0$. Those connectivity components have a positive surface measure $|\Gamma _j|$ and a positive distance of each other. So the integral mean

$$\begin{aligned} \gamma _j:= \int _{\Gamma _j} \varrho \cdot \nu \, \ \textrm{d} s /|\Gamma _j| \end{aligned}$$

is well defined and depends continuously on $\varrho \in H(\mathop {\textrm{div}}\limits =0,\Omega )$ in the sense that $|\gamma _j|\le c_1\Vert \varrho \Vert$ (recall $\mathop {\textrm{div}}\limits \varrho =0$) for each $j=0,\dots , J$ and $c_1>0$. This constant $c_1$ and the constants $c_2,c_3,c_4$ below exclusively depend on the domain $\Omega$. The finite real numbers $\gamma _0,\dots ,\gamma _J$ define the Neumann data for the harmonic function $z\in H^1(\Omega )/\mathbb R$ with

$$\begin{aligned} \Delta z=0 \text { in }\Omega \quad \text {and}\quad \partial z/\partial \nu = \gamma _j\text { on }\Gamma _j \text { for all }j=0,\dots , J. \end{aligned}$$

The elliptic regularity theory for polyhedral domains lead to $z\in H^{1+\alpha }(\Omega )$ for some $\alpha >1/2$ and $c_2>0$ with $\Vert z \Vert _{H^{1+\alpha }(\Omega )}\le c_2\, (|\gamma _0|+\dots +|\gamma _J|)$. The Raviart-Thomas interpolation operator defines a bounded linear operator on $H(\mathop {\textrm{div}}\limits ,\Omega )\cap L^p(\Omega ;\mathbb R^n)$ for $p>2$. It is generally accepted that, for $\alpha >0$ and $\nabla z\in H(\mathop {\textrm{div}}\limits =0,\Omega )\cap H^\alpha (\Omega ;\mathbb R^n)$, the Fortin interpolation $I_{\textrm{F}} \nabla z \in RT_0(\mathcal T)\cap H(\mathop {\textrm{div}}\limits =0,\Omega )$ is well defined and $\Vert I_{\textrm{F}} \nabla z\Vert \le c_3 \Vert \nabla z\Vert _{H^\alpha (\Omega )}$ follows for some $c_3>0$. The additional property $\nabla z\in L^p(\Omega ;\mathbb R^n)$ for some $p>2$ allows the definition of $\int _F \nabla z\cdot \nu _F \ \textrm{d} s$ as a Lebesgue integral over a facet $F\in \mathcal {F}$. One consequence for the boundary facets is the vanishing integral

$$\begin{aligned} \int _{\Gamma _j} (\varrho -I_{\textrm{F}} \nabla z)\cdot \nu \, \textrm{d} s =0\quad \text {for all }j=0,\dots , J. \end{aligned}$$

Since $\varrho -I_{\textrm{F}} \nabla z\in H(\mathop {\textrm{div}}\limits =0,\Omega )$ is divergence-free, Theorems 3.1 and 3.4 in [33] prove the existence of $c_4>0$ and $\beta \in H^1(\Omega ;\mathbb R^N)$ with

$$\begin{aligned} \varrho =I_{\textrm{F}} \nabla z+ \mathop {\textrm{Curl}}\limits \beta \quad \text {and}\quad |\!|\!|\beta |\!|\!|\le c_4 \Vert \varrho -I_{\textrm{F}} \nabla z\Vert . \end{aligned}$$

Recall that $\varrho \perp I_{\textrm{F}}\nabla z$ and $\varrho \perp \mathop {\textrm{curl}}\limits \beta _C\in RT_0(\mathcal {T})\cap H(\mathop {\textrm{div}}\limits =0, \Omega )$. This concludes the proof of (10) with $C_{\textrm{H}}:=\sqrt{1 + (c_1c_2c_3(1+J))^2}c_4$. $\square$

The subsequent version of the trace inequality on the facets $\mathcal {F}$ leads to the piecewise constant $\ell \in P_0(\mathcal {F})$ defined by

$$\begin{aligned} \ell (F):={\left\{ \begin{array}{ll}{} (n+1)h_{T}^{2}|F|/|T|&{}\text {for }F\in \mathcal {F}(\partial \Omega )\cap \mathcal {F}(T),\\ (n+1) |F|/(h_{T_+}^{-2}|T_+| + h_{T_-}^{-2}|T_-|)&{}\text {for }F=\partial T_+\cap \partial T_-\in \mathcal {F}(\Omega ). \end{array}\right. } \end{aligned}$$

Lemma 3

(trace inequality) Any $f\in H^1(\Omega )$ satisfies

$$\begin{aligned} \sum ^{}_{F\in \mathcal {F}} \ell (F)^{-1}\left\| f\right\| _{L^2(F)}^{2} \le \left\| h_\mathcal {T}^{-1}f\right\| ^{2} + \frac{2C_{\textrm{tr}}}{n}\left\| h_\mathcal {T}^{-1}f\right\| {\left| \hspace{-1.0625pt}\left| \hspace{-1.0625pt}\left| f \right| \hspace{-1.0625pt}\right| \hspace{-1.0625pt}\right| } \end{aligned}$$

with the constant $C_{\textrm{tr}}:=\max _{T\in \mathcal {T}}\max _{x\in T}|x-\textrm{mid}(T)|/h_T<n/(n+1)$.

Proof

The center of inertia $\textrm{mid}(T)=\sum _{j=0}^n P_j/(n+1)$ of the n-simplex $T=\textrm{conv}\{P_0,..., P_n\}\in \mathcal {T}$ and the $n+1$ faces $F_j=\textrm{conv}\{P_0,..., P_{j-1},P_{j+1},..., P_n\}\in {\mathcal {F}}(T)$ for $j=0,..., n$ give rise to the decomposition of T into $n+1$ sub-simplices $T_j' = \textrm{conv}(F_j, \textrm{mid}(T))$ with volume

$|T_j'|=|T|/(n+1)$.

Standard arguments like the trace identity on $T_j'\subset T$ [17, Lemma 2.1] for $|f|^2\in W^{1,1}(T')$ and a Cauchy inequality show

$$\begin{aligned} \frac{1}{|F_j|}\left\| f\right\| _{L^2(F_j)}^{2}\le \frac{1}{|T_j'|}\left\| f\right\| _{L^2(T_j')}^{2}+ \frac{2}{n|T_j'|}\max _{x\in T_j'}|x-\textrm{mid}(T)|\left\| f\right\| _{L^2({T_j'})}^{}\Vert \nabla f\Vert _{L^2(T')}. \end{aligned}$$

The distance $\max _{x\in T_j'}|x-\textrm{mid}(T)|=|P_k - \textrm{mid}(T)|\le C_{\textrm{tr}}h_T$ is attained at a vertex $P_k$ for $k\in \{0,..., j-1, j+1,..., n\}$. Since the centroid $\textrm{mid}(T)$ divides each median of T in the ratio n to 1 and the length of each median is strictly bounded by $h_T$,

the bound $C_{\textrm{tr}}<n/(n+1)$ follows and cannot be improved in the absence of further assumptions on the shape of the simplex T. Since $|T_j'| = |T|/(n+1)$, the previously displayed estimate leads to

$$\begin{aligned} \frac{|T|}{(n+1)h_{T}^2|F_j|}\left\| f\right\| _{L^2(F_j)}^{2}&\le \left\| h_{T}^{-1}f\right\| _{L^2(T_j')}^{2}+ \frac{2C_{\textrm{tr}}}{n}\left\| h_T^{-1}f\right\| _{L^2({T_j'})}^{}\Vert \nabla f\Vert _{L^2(T_j')}. \end{aligned}$$

Let $\mathcal {T}'$ be the refinement of $\mathcal {T}$, obtained by replacing $T\in \mathcal {T}$ with $T_0',..., T_d'$ from above. The triangulation $\mathcal {T}'$ allows for the facet based decomposition $\{\omega '(F)\}_{F\in \mathcal {F}}$ of $\Omega$, where $\omega '(F)$ is either the patch $\omega '(F)=\textrm{int}(T'_+\cup T'_-)$ for an interior facet $F=T'_+\cap T'_-$ or $\omega '(F)=\textrm{int}(T')$ for $F\in \mathcal {F}(\partial \Omega )\cap \mathcal {F}(T')$. This establishes, for any $F\in \mathcal {F}$, the estimate

$$\begin{aligned} \ell (F)^{-1}\left\| f\right\| _{L^2(F)}^{2}&\le \left\| h_{\mathcal {T}}^{-1}f\right\| _{L^2(\omega '(F))}^{2} + \frac{2C_{\textrm{tr}}}{n}\left\| h_{\mathcal {T}}^{-1}f\right\| _{L^2(\omega '(F))}\Vert \nabla f\Vert _{L^2(\omega '(F))}. \end{aligned}$$

Since the family $\{\omega '(F): F\in \mathcal {F}\}$ has no overlap, the sum of the last displayed inequality over all $F \in \mathcal {F}$ and a Cauchy inequality conclude the proof of Lemma 3. $\square$

The next lemma utilizes a quasi-interpolation operator $J: H^1(\Omega ) \rightarrow S^1(\mathcal {T})$ with the restriction $J(V) \subset S^1_0(\mathcal {T})$, e.g., $J = J_1 \circ I_{NC}$ from [19, Section 5] with explicit constants for $n=2$, and the approximation and stability properties

$$\begin{aligned} \left\| h_\mathcal {T}^{-1}(\varphi -J\varphi )\right\|&\le C_{1} {\left| \hspace{-1.0625pt}\left| \hspace{-1.0625pt}\left| \varphi \right| \hspace{-1.0625pt}\right| \hspace{-1.0625pt}\right| }\quad \text {and}\quad {\left| \hspace{-1.0625pt}\left| \hspace{-1.0625pt}\left| \varphi -J\varphi \right| \hspace{-1.0625pt}\right| \hspace{-1.0625pt}\right| }\le C_\textrm{st}{\left| \hspace{-1.0625pt}\left| \hspace{-1.0625pt}\left| \varphi \right| \hspace{-1.0625pt}\right| \hspace{-1.0625pt}\right| } \end{aligned}$$

(11)

for constants $C_{1}$ and $C_\textrm{st}$ exclusively depending on the shape-regularity of $\mathcal {T}$. For the precise definition of $J_1$ and $I_{NC}$, we refer to [19, eq. (47) and Section 5]. Recall the constant $C_{\textrm{tr}}$ from Lemma 3 and set $C_{2}:=(C_{1}(C_{1}+2C_{\textrm{tr}}C_\textrm{st}\,/n))^{1/2}$.

Lemma 4

(equilibrium) Suppose that ${G}\in H(\mathop {\textrm{div}}\limits ,\mathcal T)$ and $f\in L^2(\Omega )$ satisfy (5) and suppose $({G}|_T)|_F\cdot \nu _F\in L^2(F)$ for all $F\in {\mathcal {F}}(T)$ and $T\in \mathcal T$. Then

$$\begin{aligned} |\!|\!|f+\mathop {\textrm{div}}\limits {G}|\!|\!|_* \le {C_{1}}\Vert h_\mathcal T( f+{\mathop {\textrm{div}}\limits }_{\textrm{pw}} {G}) \Vert + C_{2}\sqrt{\sum _{F\in {\mathcal {F}}(\Omega )} \ell (F) \Vert [ {G}]_F\cdot \nu _F \Vert _{L^2(F)}^2}. \end{aligned}$$

Proof

Given $\varphi \in V$ with $|\!|\!|\varphi |\!|\!| =1$, set $\psi :=\varphi -\varphi _C$ for some quasi-interpolation $\varphi _C:=J\varphi \in S^1_0(\mathcal T)$ with (11). Since (5) implies $\langle f+\mathop {\textrm{div}}\limits {G},\varphi \rangle =\langle f+\mathop {\textrm{div}}\limits {G},\psi \rangle$, a piecewise integration by parts and the collection of jump contributions show

$$\begin{aligned} \langle f+\mathop {\textrm{div}}\limits {G},\varphi \rangle = (f + {\mathop {\textrm{div}}\limits }_{\textrm{pw}}{G}, \psi )_{L^2(\Omega )} -\sum _{F\in {\mathcal {F}}(\Omega )} \langle [{G}]_F\cdot \nu _F, \psi \rangle _{L^2(F)}. \end{aligned}$$

(12)

The first bound follows from a Cauchy inequality and (11),

$$\begin{aligned} (f + {\mathop {\textrm{div}}\limits }_{\textrm{pw}}{G}, \psi )_{L^2(\Omega )}&\le \Vert h_\mathcal {T}(f+{\mathop {\textrm{div}}\limits }_{\textrm{pw}}{G})\Vert \, \Vert h_\mathcal {T}^{-1}\psi \Vert \nonumber \\&\le {C_{1}} \Vert h_\mathcal {T}(f+{\mathop {\textrm{div}}\limits }_{\textrm{pw}}{G})\Vert \, |\!|\!|\varphi |\!|\!|. \end{aligned}$$

(13)

The second bound additionally exploits the trace inequality of Lemma 3,

$$\begin{aligned} \sum _{F\in {\mathcal {F}}(\Omega )}&\langle [{G}]_F\cdot \nu _F, \psi \rangle _{L^2(F)} \nonumber \\&\le \sqrt{\sum _{F\in {\mathcal {F}}(\Omega )} \ell (F)\Vert [{G}]_F\cdot \nu _F\Vert _{L^2(F)}^2} \sqrt{\sum _{F\in {\mathcal {F}}(\Omega )}\ell (F)^{-1}\Vert \psi \Vert _{L^2(F)}^2}\nonumber \\&\le \sqrt{{C_{1}}^2 + \frac{2C_{\textrm{tr}}}{n}C_\textrm{st}{C_{1}}}\sqrt{\sum _{F\in {\mathcal {F}}(\Omega )} \ell (F)\Vert [{G}]_F\cdot \nu _F\Vert _{L^2(F)}^2} \end{aligned}$$

(14)

with ${\left| \hspace{-1.0625pt}\left| \hspace{-1.0625pt}\left| \varphi \right| \hspace{-1.0625pt}\right| \hspace{-1.0625pt}\right| }=1$ in the last step. Since (13)–(14) hold for all $\varphi \in V$ with $|\!|\!|\varphi |\!|\!|=1$, the supremum in (12) over all such $\varphi$ concludes the proof.$\square$

The final ingredient for the proof of Theorem 1 controls the second term $\delta$ in the decomposition of Lemma 1 for $\varrho :={G}-\nabla w$. Recall $C_{\textrm{H}}$ from Lemma 2 and ${C_{1}}$ from (11), and $C_{2}$ from page 9.

Lemma 5

(conformity) Suppose the divergence-free function $\varrho \in H(\mathop {\textrm{div}}\limits =0,\Omega )\cap H(\mathop {\textrm{curl}}\limits ,\mathcal T)$ is $L^2$ orthogonal onto $RT_0(\mathcal T) \cap H(\mathop {\textrm{div}}\limits =0,\Omega )$ and satisfies $(\varrho |_T)|_F\times \nu _F\in L^2(F)$ for all $F\in {\mathcal {F}}(T)$ and $T\in T$. Then

$$\begin{aligned} C_{\textrm{H}}^{-1}\Vert \varrho \Vert \le {C_{1}}\Vert h_\mathcal T{\mathop {\textrm{curl}}\limits }_{\textrm{pw}} \varrho \Vert + C_{2}\sqrt{\sum _{F\in {\mathcal {F}}} \ell (F) \Vert [ \varrho ]_F \times \nu _F \Vert _{L^2(F)}^2}.\end{aligned}$$

Proof

Lemma 2 provides $\beta \in H^1(\Omega ; \mathbb R^N)$ with (10) for a (component-wise) quasi-interpolation $\beta _C\in S^1(\mathcal {T})^N$ with (11) as in the proof of Lemma 4; set $\psi :=\beta -\beta _C$. A piecewise integration by parts and the collection of jump contributions shows

$$\begin{aligned} \Vert \varrho \Vert ^2 = \int _\Omega \varrho \cdot \mathop {\textrm{Curl}}\limits (\beta -\beta _C)dx = \int _\Omega \psi \cdot {\mathop {\textrm{curl}}\limits }_{\textrm{pw}} \varrho \, dx +\sum _{F\in {\mathcal {F}}} \int _F \psi [ \varrho ]_F\times \nu _F \, ds. \end{aligned}$$

Stability and approximation properties of the quasi-interpolation (11) and the trace inequality of Lemma 3 eventually lead to

$$\begin{aligned} \left\| \varrho \right\| ^{2}&\le {C_{1}}\left\| h_\mathcal {T}{\mathop {\textrm{curl}}\limits }_{\textrm{pw}} \varrho \right\| _{}^{}|\!|\!|\beta |\!|\!| + C_{2} \sqrt{\sum _{F\in {\mathcal {F}}} \ell (F)\Vert [\varrho ]_F\times \nu _F\Vert _{L^2(F)}^2}|\!|\!|\beta |\!|\!|. \end{aligned}$$

In fact, the routine estimation with element and jump terms is completely analogous to the proof of Lemma 4 and leads to the same constants ${C_{1}}, C_{2}$. This and $|\!|\!|\beta |\!|\!|\le C_{\textrm{H}}\Vert \varrho \Vert$ conclude the proof. $\square$

Proof

(Theorem 1) The trace of ${G}|_T\in H^1(T; \mathbb R^n)$ is well defined on any facet $F\in \mathcal {F}(T)$ of the simplex $T\in \mathcal {T}$. Lemma 1 provides $w\in V=H^1_0(\Omega )$ with $\Vert \nabla u-{G}\Vert ^2 = |\!|\!| f + \mathop {\textrm{div}}\limits {G}|\!|\!|_{*}^2+\Vert {G}- \nabla w\Vert ^2$. Since ${G}$ satisfies (5), Lemma 4 establishes

$$\begin{aligned} |\!|\!| f + \mathop {\textrm{div}}\limits {G}|\!|\!|_{*} \le C_\textrm{1}\Vert h_\mathcal T( f+{\mathop {\textrm{div}}\limits }_{\textrm{pw}} {G}) \Vert + C_{2}\sqrt{\sum _{F\in {\mathcal {F}}(\Omega )} \ell (F) \Vert [ {G}]_F\cdot \nu _F \Vert _{L^2(F)}^2}. \end{aligned}$$

The assumption (6) on ${G}$ and an integration by parts prove that Lemma 5 is applicable to $\varrho :={G}-\nabla w\in H(\mathop {\textrm{div}}\limits =0, \Omega )\cap H(\mathop {\textrm{curl}}\limits , \mathcal {T}; \mathbb R^n)$. Since $\mathop {\textrm{curl}}\limits \nabla w = 0$ and $\nabla w\times \nu = 0$, this reveals

$$\begin{aligned} C_{\textrm{H}}^{-1}\Vert {G}- \nabla w\Vert \le C_\textrm{1}\Vert h_\mathcal T{\mathop {\textrm{curl}}\limits }_{\textrm{pw}} {G}\Vert + C_{2}\sqrt{\sum _{F\in {\mathcal {F}}} \ell (F) \Vert [{G}]_F \times \nu _F \Vert _{L^2(F)}^2}. \end{aligned}$$

The above estimates together with the decomposition of Lemma 1 establish $\eta (\mathcal {T}, {G})$ as a GUB for the error $\Vert \nabla u - {G}\Vert$. $\square$

Example 1

(constants for right-isosceles triangles) In two space dimensions, $\Vert \mathop {\textrm{Curl}}\limits \bullet \Vert = |\!|\!|\bullet |\!|\!|$ and so $C_{\textrm{H}}\le 1$ for a simply connected domain $\Omega$ in Lemma 2. The choice $J:=J_1\circ I_{\textrm{NC}}$ from [19, Section 5] of the quasi-interpolation operator J in the proof of Lemma 4 allows for the explicit estimates

$$\begin{aligned} {C_{1}}\le \sqrt{48^{-1} + j_{1,1}^{-2} + c_{\textrm{apx}}^2} \text { and } C_\textrm{st}\le 1+\sqrt{72}c_{\textrm{apx}}, \end{aligned}$$

where $j_{1,1} = 3.8317$ denotes the first positive root of the first Bessel function. For triangulations into right-isosceles triangles, the constant $c_{\textrm{apx}}\le \sqrt{3}/(2 - 2\cos (\pi /\max \{4,M_{\textrm{bd}}\}))$ from [19, Lemma 4.8] depends on the domain by the maximal number $M_{\textrm{bd}}\le 4\max \{\pi ,\omega _{\textrm{max}}\}/\pi \le 8$ of triangles sharing a boundary vertex. Given the maximal interior angle $\omega _{\textrm{max}}$ of $\Omega$, Table 1 displays those constants for the maximal possible value $M_{\textrm{bd}}= 4\max \{\pi ,\omega _{\textrm{max}}\}/\pi$. The geometric quantity $\max _{x\in T}|x-\textrm{mid}(T)|$ equals two-thirds of the maximum median of T. Thus, $C_{\textrm{tr}}=\sqrt{5}/(3\sqrt{2})\le 0.5271$ and $\ell (F) = 6h_F$ for interior edges $F\in \mathcal {F}(\Omega )$ and $\ell (F)=12h_F$ for boundary edges $F\in \mathcal {F}(\partial \Omega )$ of triangulations into right-isosceles triangles. Consequently,

$$\begin{aligned} {C_{1}}&\le \sqrt{48^{-1} + j_{1,1}^{-2} + c_{\textrm{apx}}^2}=:C_\mathcal {T}, \end{aligned}$$

(15)

$$\begin{aligned} {C_{2}}&\le \sqrt{C_\mathcal {T}(C_\mathcal {T}+0.5271(1+\sqrt{72}c_{\textrm{apx}}))}=:C_\mathcal {F}. \end{aligned}$$

(16)

3 Explicit residual-based a posteriori HHO error estimator

The arguments from Sect. 2 apply to the HHO method and result in a stabilization-free reliable a posteriori error control. In combination with the efficiency estimate from this section, this leads to a new explicit residual-based a posteriori error estimator for the HHO method that is equivalent to the error up to data oscillations.

3.1 Hybrid high-order methodology

The HHO ansatz space reads $V_h:=P_k(\mathcal {T}) \times P_k(\mathcal {F}(\Omega ))$ for $k\in {\mathbb N}_0$ with the subspace $P_k(\mathcal {F}(\Omega )) \subset P_k(\mathcal {F})$ of piecewise polynomials $p\in P_k(\mathcal {F})$ under the convention $p_{|\partial \Omega }=0$. The interpolation $\textrm{I}: V \rightarrow V_h$ maps $v \in V$ onto $\textrm{I} v :=(\Pi _k v, \Pi _{{\mathcal {F}},k} v) \in V_h$. Given any $v_h = (v_{\mathcal {T}},v_{\mathcal {F}})\in V_h$, the reconstruction operator $R: V_h\rightarrow P_{k+1}(\mathcal {T})$ defines the unique piecewise polynomial $Rv_h \in P_{k+1}(\mathcal {T})$ with $\Pi _{0}(Rv_h - v_{\mathcal {T}}) = 0$ such that, for all $w_{k+1} \in P_{k+1}(\mathcal {T})$,

$$\begin{aligned}&a_{\textrm{pw}}(Rv_h,w_{k+1})\nonumber \\&\qquad = a_{\textrm{pw}}(v_{\mathcal {T}}, w_{k+1}) -\sum _{T\in \mathcal {T}}\langle v_{\mathcal {T}|T}-v_{\mathcal {F}}, \nabla w_{k+1|T} \cdot \nu _T \rangle _{L^2(\partial T)}. \end{aligned}$$

(17)

Let $u_h\in V_h$ solve the HHO discrete formulation of (3) with

$$\begin{aligned} a_h(u_h,v_h)= (f,v_{\mathcal {T}})_{L^2(\Omega )} \qquad \text {for all } v_h=(v_{\mathcal {T}},v_{\mathcal {F}})\in V_h\end{aligned}$$

(18)

for the HHO bilinear form

$$\begin{aligned} a_h(u_h,v_h)&:=a_{\textrm{pw}}( Ru_h, Rv_h)+ s_h(u_h, v_h) \end{aligned}$$

(19)

and the stabilization term $s_h(u_h, v_h)$ from (2). Given any $w_C \in S^1_0(\mathcal {T})= P_1(\mathcal {T})\cap H^1_0(\Omega )$, the definition of the reconstruction operator R in (17) verifies $R \textrm{I} w_C=w_C$ with the interpolation $\textrm{I}$ onto $V_h$. Hence, $S_{TF}\textrm{I} w_C=0$ vanishes for all $F \in \mathcal {F}(T)$ and $T \in \mathcal {T}$. This and (18) show, for all $w_C\in S^1_0(\mathcal {T})$, that

$$\begin{aligned} a_{\textrm{pw}}(Ru_h, w_C) = (f, \Pi _k w_C)_{L^2(\Omega )} = (\Pi _k f, w_C)_{L^2(\Omega )}. \end{aligned}$$

(20)

3.2 Explicit a posteriori error estimator

As a result of (20), $\nabla _{\textrm{pw}}Ru_h$ satisfies the solution property (5) if $k\ge 1$ and, in the lowest order case $k=0$, (5) holds with f replaced by $\Pi _0 f$. This allows the application of the theory from Sect. 2 to the HHO method with minor modifications for the case $k=0$. Define the error estimator contributions

$$\begin{aligned} \begin{aligned} \eta ^2_{\textrm{res}, 1}(\mathcal {T})&:={\left\{ \begin{array}{ll}\Vert h_\mathcal {T}( f+\Delta _{\textrm{pw}} Ru_h)\Vert ^2&{}\text {for } k\ge 1,\\ \Vert h_\mathcal {T}\Pi _0 f\Vert ^2&{}\text {for } k=0, \end{array}\right. }\\ \eta ^2_{\textrm{res}, 2}(\mathcal {T})&:={\left\{ \begin{array}{ll}0&{}\text {for } k\ge 1,\\ \textrm{osc}_{0}^2(f, \mathcal {T})&{}\text {for } k=0, \end{array}\right. }\\ \eta ^2_{\textrm{res}, 3}(\mathcal {T})&:=\sum _{F\in {\mathcal {F}}(\Omega )} \ell (F)\Vert [ \nabla _{\textrm{pw}}Ru_h ]_F\cdot \nu _F \Vert _{L^2(F)}^2,\\ \eta ^2_{\textrm{res}, 4}(\mathcal {T})&:=\sum _{F\in {\mathcal {F}}} \ell (F)\Vert [ \nabla _{\textrm{pw}}Ru_h ]_F\times \nu _F \Vert _{L^2(F)}^2. \end{aligned} \end{aligned}$$

(21)

Since $\nabla _{\textrm{pw}}Ru_h$ is a piecewise gradient, its piecewise $\mathop {\textrm{curl}}\limits$ vanishes. This leads to the explicit residual-based a posteriori error estimator

$$\begin{aligned} \eta _{\textrm{res}}^2(\mathcal {T}):=\left( {C_{1}}\eta _{\textrm{res}, 1}(\mathcal {T}) + C_P\eta _{\textrm{res}, 2}(\mathcal {T}) +C_{2}\eta _{\textrm{res}, 3}(\mathcal {T})\right) ^2+ C_{\textrm{H}}^2C_{2}^2\eta _{\textrm{res}, 4}^2(\mathcal {T}). \end{aligned}$$

(22)

(Recall ${C_{1}}, C_{2}$ from Lemma 4 and $C_{\textrm{H}}$ from Lemma 2 as well as the Poincaré constant $C_P\le \pi ^{-1}$.) The main result of this section verifies the assumptions in Theorem 1 and proves reliability and efficiency of $\eta _{\textrm{res}}(\mathcal {T})$.

Theorem 2

(residual-based GUB for HHO) Let $u\in V$ solve the Poisson equation (3) and let $u_h \in V_h$ solve the discrete formulation (18). Then

$$\begin{aligned} |\!|\!| u- Ru_h |\!|\!|_{\textrm{pw}} \le \eta _{\textrm{res}}(\mathcal {T})\le C_3\big ( |\!|\!| u- Ru_h |\!|\!|_{\textrm{pw}}+\textrm{osc}_{q}(f, \mathcal T)\big ) \end{aligned}$$

and $\textrm{osc}_{k-1}(f, \mathcal T)\le C_{4} \eta _{\textrm{res}}(\mathcal {T})$ hold for any $q\in {\mathbb N}_0$. The constants $C_{3}$ and $C_{4}$ exclusively depend on k, q and on the shape-regularity of the triangulation $\mathcal T$.

3.3 Proof of Theorem 2

The orthogonality of $\nabla _\textrm{pw}R u_h$ to the divergence-free Raviart-Thomas space of lowest degree is an assumption in Theorem 1 and verified below.

Lemma 6

(orthogonality) The piecewise gradients $\nabla _{\textrm{pw}}R V_h$ are $L^2$ orthogonal to the space  $RT_0(\mathcal {T})\cap H(\textrm{div} = 0,\Omega )$, i.e., any $v_h\in V_h$ and $q_{RT}\in RT_0(\mathcal {T})\cap H(\textrm{div} = 0,\Omega )$ satisfy

$$\begin{aligned} (\nabla _{\textrm{pw}}Rv_h, q_{RT})_{L^2(\Omega )} = 0. \end{aligned}$$

(23)

Proof

Given any $q_{RT}\in RT_0(\mathcal {T})\cap H(\textrm{div} = 0,\Omega )$, $\mathop {\textrm{div}}\limits q_{RT}=0$ shows $q_{RT}\in P_0(\mathcal {T};\mathbb R^n)$ [28, Lemma 14.9]. Since $P_0(\mathcal {T};\mathbb R^n) = \nabla _{\textrm{pw}}P_1(\mathcal {T})$, there exists a piecewise affine function $\phi _1 \in P_1(\mathcal {T})$ with $q_{RT}= \nabla _\textrm{pw}\phi _1$ a.e. in $\Omega$. This and the definition of $Rv_h$ from (17) imply, for any $v_h = (v_{\mathcal {T}},v_{\mathcal {F}})\in V_h$, that

$$\begin{aligned} (\nabla _{\textrm{pw}}Rv_h,q_{RT})_{L^2(\Omega )}&= a_\textrm{pw}(Rv_h, \phi _1)\\&= a_\textrm{pw}(v_\mathcal {T},\phi _1) - \sum \limits _{T\in {\mathcal {T}}} \langle v_{\mathcal {T}}|_T-v_{\mathcal {F}}, \nabla _\textrm{pw}\phi _1 \cdot \nu _T \rangle _{L^2(\partial T)}. \end{aligned}$$

This, a piecewise integration by parts, and $\Delta _\textrm{pw}\phi _1 \equiv 0$ lead to

$$\begin{aligned} (\nabla _{\textrm{pw}} Rv_h,q_{RT})_{L^2(\Omega )}&= \sum \limits _{T\in {\mathcal {T}}}\langle v_{\mathcal {F}},\nabla _\textrm{pw}\phi _1 \cdot \nu _T \rangle _{L^2(\partial T)}\nonumber \\&= \sum \limits _{F\in {\mathcal {F}}}\langle v_{\mathcal {F}}, [q_{RT}\cdot \nu _F]_F \rangle _{L^2(F)}. \end{aligned}$$

(24)

Since $q_{RT}\in RT_0(\mathcal {T})$ has continuous normal components, the jump term $\langle v_{\mathcal {F}}, [q_{RT}]_F\cdot \nu _F\rangle _{L^2(F)} =0$ vanishes for all $F\in \mathcal F(\Omega )$. This, $v_{\mathcal {F}}\equiv 0$ on $\partial \Omega$, and (24) conclude $(\nabla _{\textrm{pw}} Rv_h,q_{RT})_{L^2(\Omega )} = 0$. $\square$

The following lemma concerns the efficiency of the jump contributions. Each facet $F\in \mathcal {F}$ has at most two adjacent simplices that define a triangulation $\mathcal {T}(F):=\{T\in \mathcal {T}: F\in \mathcal {F}(T)\}$ of the facet-patch $\omega (F):=\textrm{int}(\bigcup _{T\in \mathcal {T}(F)}T)$.

Lemma 7

(efficiency of jumps) The solution $u \in V$ to (3) and the discrete solution $u_h \in V_h$ to (18) satisfy (a) for all $F \in \mathcal {F}$ and (b) for all $F \in \mathcal {F}(\Omega )$.

• (a) $h_F^{1/2}\Vert [\nabla _{\textrm{pw}} Ru_h]_F\times \nu _F \Vert _{L^2(F)} \lesssim \min _{v \in V} \Vert \nabla v - \nabla _{\textrm{pw}} Ru_h\Vert _{L^2(\omega (F))}$,

• (b) $h_F^{1/2}\Vert [\nabla _{\textrm{pw}} Ru_h]_F\cdot \nu _F \Vert _{L^2(F)} \lesssim \Vert \nabla u - \nabla _{\textrm{pw}} Ru_h\Vert _{L^2(\omega (F))} + \textrm{osc}_k(f, \mathcal {T}(F))$.

Proof

The proof is based on the following extension argument. Given a polynomial $p\in P_k(F)$ of degree at most k along the side $F\in \mathcal {F}$, the coefficients determine a polynomial (also denoted by p) along the hyperplane H that enlarges F. The intersection ${\widehat{F}}:= H\cap \textrm{conv}({ \omega (F)})$ of the hyperplane H with the convex hull of the facet-patch $\omega (F)$ may be strictly larger than F. The shape-regularity of $\mathcal {T}$ bounds the size of ${\widehat{F}}$ in terms of F and an inverse estimate leads to a bound $\Vert p \Vert _{L^\infty ({\widehat{F}})} \le C(k)\Vert p \Vert _{L^\infty (F)}$ with a constant C(k) that depends on the shape-regularity of $\mathcal {T}$ and on k. The extension of p from H to $\mathbb R^n$ by constant values along the side normal $\nu _F$ leads to a polynomial ${\widehat{p}}\in P_k(\mathbb R^ n)$ with

$$\begin{aligned} \Vert {\widehat{p}} \Vert _{L^\infty (\omega (F))} \le \Vert p \Vert _{L^\infty ({\widehat{F}})} \le C(k)\Vert p \Vert _{L^\infty (F)}. \end{aligned}$$

(25)

Proof of (a). The tangential jump $\varrho _F:= [\nabla _{\textrm{pw}} Rv_h]_F\times \nu _F \in P_{k}(F;{\mathbb {R}}^N)$ is a polynomial in $N=2n-3$ components on $F\in \mathcal {F}$ for $n=2,3$. Let $p=\varrho _F(j)$ be one of the components of $\varrho _F\in P_k(F)^N$, for $j=1,...,N$, and extend it as explained above to ${\hat{p}}\in P_k(\mathbb R^n)$ and call this ${\widehat{\varrho }}(j)$ in the vector $\widehat{\varrho _F}\in P_k(\mathbb R^n;{\mathbb {R}}^N)$. The proof involves the piecewise polynomial facet-bubble function $b_F:=n^n \Pi _{j=1}^n \varphi _j$ for the n nodal basis function $\varphi _1, \dots , \varphi _n \in S^1(\mathcal T)$ associated with the vertices of F. An inverse estimate [38, Proposition 3.37] shows

$$\begin{aligned} \Vert \varrho _F \Vert _{L^2(F)}^2 \lesssim \Vert b_F^{1/2}\varrho _F \Vert _{L^2(F)}^2=\langle b_F\varrho _F, \varrho _F\rangle _{L^2(F)}. \end{aligned}$$

(26)

Since $\varrho := b_F\widehat{\varrho _F} \in S^{k+n}_0(\mathcal {T}(F);{\mathbb {R}}^N)$ vanishes on $\partial \omega (F){\setminus }\textrm{int}(F)$, (26) and a piecewise integration by parts show

$$\begin{aligned} \Vert \varrho _F \Vert _{L^2(F)}^2 \lesssim \langle \varrho , [\nabla _{\textrm{pw}} Rv_h]_F\times \nu _F \rangle _{L^2(F)} = (\mathop {\textrm{curl}}\limits \varrho , \nabla _{\textrm{pw}}R v_h )_{L^2(\omega (F))}. \end{aligned}$$

(27)

This and $(\mathop {\textrm{curl}}\limits \varrho , \nabla v)_{L^2(\Omega )} = 0$ for any $v\in V$ imply

$$\begin{aligned} \Vert \varrho _F \Vert _{L^2(F)}^2&\lesssim (\mathop {\textrm{curl}}\limits \varrho ,\nabla _{\textrm{pw}}(R v_h-v))_{L^2(\omega (F))}\nonumber \\&\le \Vert \mathop {\textrm{curl}}\limits \varrho \Vert _{L^2(\omega (F))} \Vert \nabla _{\textrm{pw}}(R v_h-v) \Vert _{L^2(\omega (F))}. \end{aligned}$$

(28)

An inverse estimate, $\Vert b_F\Vert _{L^\infty (\omega (F))} = 1$, and (25) imply

$$\begin{aligned} \Vert \mathop {\textrm{curl}}\limits \varrho \Vert _{L^2(\omega (F))}&\lesssim \Vert \nabla \varrho \Vert _{L^2(\omega (F))} \lesssim h_F^{-1+n/2} \Vert \varrho \Vert _{L^\infty (\omega (F))}\nonumber \\&\lesssim h_F^{-1+n/2} h_F^{-(n-1)/2}\Vert \varrho _F \Vert _{L^2(F)} = h_F^{-1/2}\Vert \varrho _F \Vert _{L^2(F)}. \end{aligned}$$

(29)

In combination with (28), this concludes the proof of (a).

Proof of (b). The efficiency of normal jumps (b) follows from the arguments for conforming FEMs, cf. [38, Section 1.4.5]; further details are omitted. $\square$

The following lemma reveals that the order $k \ge {\mathbb {N}}_0$ of the oscillations $\textrm{osc}_{k}(f, T)$ in Lemma 7 (b) can be any natural number. It is certainly known to the experts but hard to find in the literature. Recall the convention $\Pi _{-1}:=0$.

Lemma 8

(efficiency of lower-order oscillations) Given any simplex $T \in {\mathcal {T}}$ and parameters $k, q \in {\mathbb {N}}_0$, the solution $u \in V$ to (3) satisfies

$$\begin{aligned}&{C_{5}}^{-1}\textrm{osc}_{k-1}^2( f,T) \le \min _{v_{k+1} \in P_{k+1}(T)}\Vert \nabla _\textrm{pw}(u - v_{k+1})\Vert _{L^2(T)}^2+\textrm{osc}_{q}^2(f,T). \end{aligned}$$

(30)

The constant $C_5$ exclusively depends on q and the shape of T.

Proof

The assertion (30) is trivial for $q \le k-1$, so suppose $k \le q$. Any $v_{k+1}\in P_{k+1}(T)$ and $\varrho _T:=\Pi _q f + \Delta v_{k+1}\in P_q(T)$ satisfy

$$\begin{aligned} \textrm{osc}_{k-1}^2(f, T) \le h_T^2\Vert f + \Delta v_{k+1}\Vert _{L^2(T)}^2 = \textrm{osc}_{q}^2(f, T) + h_T^2\Vert \varrho _T\Vert _{L^2(T)}^2. \end{aligned}$$

(31)

Let $b_T\in S^{n+1}_0(T)$ with $0\le b_T\le 1=\max b_T$ denote the volume bubble-function on $T\in \mathcal {T}$. The equivalence of norms in the finite-dimensional space $P_q(T)$ provides

$$\begin{aligned} \Vert b_T^{1/2}\varrho _T\Vert _{L^2(T)} \le \Vert \varrho _T\Vert _{L^2(T)}\le C_6 \Vert b_T^{1/2}\varrho _T\Vert _{L^2(T)}. \end{aligned}$$

(32)

A more detailed analysis of the mass matrices reveals that the constant $C_{6}$ exclusively depends on the polynomial degree q. An integration by parts with $b_T\varrho _T\in S^{q+n+1}_0(T)\subset V$ and the weak formulation (3) result in

$$\begin{aligned} \Vert b_T^{1/2}\varrho _T\Vert _{L^2(T)}^2&= (\Pi _q f + \Delta v_{k+1}, b_T\varrho _T)_{L^2(T)}\\&= (\nabla (u-v_{k+1}), \nabla (b_T\varrho _T))_{L^2(T)} - (f - \Pi _q f, b_T\varrho _T)_{L^2(T)}. \end{aligned}$$

A Cauchy inequality, the inverse estimate $h_T\Vert \nabla (b_T\varrho _T)\Vert _{L^2(T)}\le C_7\Vert \varrho _T\Vert _{L^2(T)}$ with a constant $C_7$ that exclusively depends on $q+n+1$ and the shape of T, and (32) lead to

$$\begin{aligned} C_{6}^{-2}h_T\Vert \varrho _T\Vert _{L^2(T)}\le C_{7}\Vert \nabla (u-v_{k+1})\Vert _{L^2(T)} + \textrm{osc}_{q}(f, T)\big . \end{aligned}$$

(33)

The combination of (31) with (33) and a Cauchy inequality conclude the proof of (30), e.g., with $C_{5}= 1 + {C_{6}}^4(1 + {C_{7}}^2).$ $\square$

Proof

(of Theorem 2) Recall the definition of $\eta _{\textrm{res}}(\mathcal {T})$ for $k\ge 1$ and $k=0$ in (21). Since $\textrm{osc}_{k-1}^2(f, \mathcal {T})\le \eta _{\textrm{res},1}^2(\mathcal {T}) + \eta _{\textrm{res},2}^2(\mathcal {T})\lesssim \eta _{\textrm{res}}^2(\mathcal {T})$, the remaining parts of this proof discuss the reliability and efficiency of $\eta _{\textrm{res}}(\mathcal {T})$.

Lemma 6 provides the orthogonality of $\nabla _{\textrm{pw}}Ru_h\in H^1(\mathcal {T}; \mathbb R^n)$ to the divergence-free Raviart-Thomas function $RT_0(\mathcal {T})\cap H(\textrm{div} = 0,\Omega )$. This and (20) show that the assumptions in Theorem 1 hold for ${G}:=\nabla _{\textrm{pw}}Ru_h$ and $k\ge 1$, whence the reliability of $\eta _{\textrm{res}}(\mathcal {T})$ follows with a reliability constant 1.

Minor modifications to the proof of Theorem 1 lead to reliability in the case $k=0$. In fact, the only modifications required concern the upper bound of $|\!|\!|f + \textrm{div}\,\nabla _{\textrm{pw}}R u_h|\!|\!|_* \le |\!|\!|(1 - \Pi _0) f|\!|\!|_* + |\!|\!|\Pi _0 f + \textrm{div}\,\nabla _{\textrm{pw}}R u_h|\!|\!|_*$. A piecewise Poincaré inequality shows $|\!|\!|(1 - \Pi _0) f|\!|\!|_* \le C_P\textrm{osc}_{0}(\mathcal {T},f)$ with the Poincaré constant $C_P$($\le 1/\pi$ for simplices). Lemma 4 proves $|\!|\!|\Pi _0 f + \textrm{div}\,\nabla _{\textrm{pw}}R u_h|\!|\!|_* \le C_\textrm{1}\eta _{\textrm{res},1}(\mathcal {T}) + C_\textrm{2}\eta _{\textrm{res},3}(\mathcal {T})$. Hence, the decomposition of Lemma 1 and Lemma 5 result in $|\!|\!|u-Ru_h|\!|\!|_\textrm{pw}\le \eta _{\textrm{res}}(\mathcal {T})$.

This provides the reliability and it remains to verify the efficiency $\eta _{\textrm{res}}(\mathcal {T})\lesssim |\!|\!|u-Ru_h|\!|\!|_\textrm{pw}+ \textrm{osc}_{q}(f,\mathcal {T})$ for any $q\in {\mathbb {N}}_0$. The Pythagoras theorem and (33) with $\varrho _T:=\Pi _kf + \Delta R u_h\in P_k(T)$ and $v_{k+1}:=Ru_h\in P_{k+1}(T)$ lead to the local efficiency of the volume contributions

$$\begin{aligned} \Vert h_T(f+\Delta _\textrm{pw}Ru_h)\Vert _{L^2(T)}^2&=\textrm{osc}_k^2(f, T) + \Vert h_T\varrho _T\Vert _{L^2(T)}^2\\&\lesssim \Vert \nabla (u- Ru_h)\Vert _{L^2(T)}^2 + \textrm{osc}_k^2(f, T). \end{aligned}$$

Lemma 7 considers the remaining terms in the error estimator and establishes their efficiency namely,

$$\begin{aligned} \sum _{F\in \mathcal {F}}h_F\Vert [ \nabla _\textrm{pw}R u_h ]_F\Vert _{L^2(F)}\lesssim |\!|\!|u-Ru_h|\!|\!|_\textrm{pw}+ \textrm{osc}_k(f,\mathcal {T}) \end{aligned}$$

with the modified jump $[\nabla _\textrm{pw}R u_h]_F = \nabla _\textrm{pw}R u_h \times \nu _F$ on boundary facets $F \in {\mathcal {F}}(\partial \Omega )$. This and Lemma 8 establish the existence of some mesh-independent constant $C_3>0$ with $C_3^{-1}\eta _{\textrm{res}}(\mathcal {T})\le |\!|\!|u-Ru_h|\!|\!|_\textrm{pw}+ \textrm{osc}_{q}(f,\mathcal {T})$ for arbitrary $q\in {\mathbb {N}}_0$. This concludes the proof.$\square$

While the focus of this paper is on the HHO methodology, the framework of Sect. 2 also applies to other skeletal methods as well. The following example covers a hybridized discontinuous Galerkin (HDG) FEM from [35] with the Lehrenfeld-Schöberl stabilization.

Example 2

(HDG FEM) Let $V_h :=P_{k+1}({\mathcal {T}}) \times P_k({\mathcal {F}}(\Omega ))$ for $k \in {\mathbb {N}}_0$. An equivalent formulation to the HDG FEM from [35] seeks $u_h \in V_h$ with

$$\begin{aligned} a_\textrm{pw}(R u_h, R v_h) + s_h(u_h,v_h) = (f, v_{\mathcal {T}}) \quad \textrm{for any } v_h = (v_{\mathcal {T}},v_{\mathcal {F}}) \in V_h. \end{aligned}$$

Here, $R: V_h \rightarrow P_{k+1}({\mathcal {T}})$ is defined as in (17) and

$$\begin{aligned} s_h(v_h,w_h) :=\sum _{T\in \mathcal {T}}\sum _{F\in \mathcal {F}(T)} h_F^{-1}\langle \Pi _{F,k}(v_{\mathcal {T}}|_T - v_{\mathcal {F}}|_F), w_{\mathcal {T}}|_T - w_{\mathcal {F}}|_F \rangle _{L^2(F)} \end{aligned}$$

for any $v_h = (v_{\mathcal {T}}, v_{\mathcal {F}}), w_h = (w_{\mathcal {T}}, w_{\mathcal {F}}) \in V_h$. This method is also known under the label of weak Galerkin FEM [39]. It is straightforward to verify that $\nabla _\textrm{pw}R u_h$ satisfies (5)–(6). Notice that (5) also holds for $k = 0$ without any modification. Therefore, Theorem 1 leads to the reliable a posteriori estimate

$$\begin{aligned}&|\!|\!|u-R u_h|\!|\!|_\textrm{pw}^2 \le \eta _{\textrm{res}}^2({\mathcal {T}}) :=\Vert h_\mathcal {T}(f + \Delta _{\textrm{pw}} Ru_h)\Vert ^2\\&~ + \sum _{F\in {\mathcal {F}}(\Omega )} \ell (F)\Vert [\nabla _{\textrm{pw}}Ru_h ]_F\cdot \nu _F \Vert _{L^2(F)}^2 + \sum _{F\in {\mathcal {F}}} \ell (F)\Vert [\nabla _{\textrm{pw}}Ru_h ]_F\times \nu _F \Vert _{L^2(F)}^2. \end{aligned}$$

The efficiency $\eta _{\textrm{res}}({\mathcal {T}}) \lesssim |\!|\!|u-R u_h|\!|\!|_\textrm{pw}$ follows from the arguments in the proofs of Lemma 7–8.

4 Equilibrium-based a posteriori HHO error analysis

The residual-based guaranteed upper bound (GUB) of the error $|\!|\!|u - R u_h|\!|\!|_\textrm{pw}$ from Sect. 3.2 employs explicit constants that may lead to overestimation in higher dimensions and for different triangular shapes. This section utilizes equilibrated flux reconstructions [1, 2, 5, 6, 30] to establish, up to the well-known Poincaré constant $C_P \le 1/\pi$, a constant-free guaranteed upper bounds for a tight error control.

4.1 Guaranteed error control

The guaranteed upper bounds of this section involves two post-processings of the potential reconstruction $R u_h \in P_{k+1}(\mathcal {T})$ of the discrete solution $u_h$ to (18). First, the patch-wise design of a flux reconstruction $Q_p \in RT_{k+p}(\mathcal {T})$ with $p \in {\mathbb {N}}_0$ from [3, 10, 31] provides an $H(\textrm{div},\Omega )$-conforming approximation to $\nabla _{\textrm{pw}}R u_h$ with the equilibrium $\Pi _{r} f + {\text { div }}Q_p = 0$ in $\Omega$ and r from (35) below. Second, the nodal average ${\mathcal {A}} R u_h \in S^{k+1}_0(\mathcal {T})\subset V$ results in an V-conforming approximation of $R u_h$ by averaging all values of the discontinuous function $R u_h$ at each Lagrange point of $S^{k+1}_0(\mathcal {T})$. This, the split (7), and the solution property (20) give rise to the guaranteed upper bound (GUB)

$$\begin{aligned} \eta _{\textrm{eq},p}^2(\mathcal {T}) := \left( C_P\textrm{osc}_{r}(f,\mathcal {T}) + \Vert Q_p- \nabla _{\textrm{pw}}R u_h\Vert \right) ^2 + |\!|\!|(1 - {\mathcal {A}})R u_h|\!|\!|_\textrm{pw}^2 \end{aligned}$$

(34)

with $r \in {\mathbb {N}}_0$ defined by

$$\begin{aligned} r :=0 \text { if } k = 0 \quad \text {and}\quad r :=k + p \text { if } k \ge 1. \end{aligned}$$

(35)

The main result of this section states the reliability and efficiency (up to data oscillations) of $\eta _{\textrm{eq},p}$ for all parameters $p \in {\mathbb {N}}_0$.

Theorem 3

(equilibrium-based GUB for HHO) Let $u \in V$ resp. $u_h \in V_h$ solve (3) resp. (18). Given a parameter $p \in {\mathbb {N}}_0$, there exists $Q_p\in RT_{k+p}({\mathcal {T}})$ such that the error estimator $\eta _{\textrm{eq},p}(\mathcal {T})$ from (34) is an efficient GUB

$$\begin{aligned} |\!|\!|u-R u_h|\!|\!|_\textrm{pw}\le \eta _{\textrm{eq},p}(\mathcal {T}) \le {C_{8}}\big (|\!|\!|u-R u_h|\!|\!|_\textrm{pw}+ \textrm{osc}_{q}(f, {\mathcal {T}})\big ). \end{aligned}$$

(36)

for any $q\in {\mathbb {N}}_0$ and $\textrm{osc}_{k-1}(f, {\mathcal {T}}) \le {C_{9}}\eta _{\textrm{eq},p}(\mathcal {T})$.

The constants $C_8$ and $C_9$ exclusively depend on the polynomial degree $k \in {\mathbb {N}}_0$, the parameter $q \in {\mathbb {N}}_0$, and the shape-regularity of $\mathcal {T}$.

At least two technical contributions for the proof of Theorem 3 are of broader interest. A first contribution to the HHO literature is the local equivalence of the original HHO stabilization $s_h$ from (2) and the alternative stabilization ${\tilde{s}}_h(v_h,v_h) :=\sum _{T \in \mathcal {T}} {\tilde{s}}_T(v_h,v_h)$ from [25] defined, for $v_h=(v_\mathcal {T}, v_\mathcal {F})\in V_h$, by

$$\begin{aligned} {\tilde{s}}_T(v_h,v_h)&:=h_T^{-2} \Vert \Pi _{T,k} (v_\mathcal {T}- R v_h)\Vert _{L^2(T)}^2 \nonumber \\&\quad + \sum _{F \in \mathcal {F}(T)} h_F^{-1} \Vert \Pi _{F,k} (v_\mathcal {F}- R v_{h}|_{T})\Vert _{L^2(F)}^2. \end{aligned}$$

(37)

A second result of separate interest in the HHO literature (cf. [25] where the efficiency in (38) is left open) is the efficiency of the stabilizations from Theorems 4–5 below,

$$\begin{aligned} {\tilde{s}}_h(v_h,v_h)^{1/2}\approx s_h(u_h,u_h)^{1/2} \lesssim |\!|\!|u - R u_h|\!|\!|_\textrm{pw}+ \textrm{osc}_{k+p}(f,\mathcal T). \end{aligned}$$

(38)

The subsequent subsection continues with some explanations on the flux reconstruction $Q_p\in RT_{k+p}({\mathcal {T}})$ that is defined by local minimization problems on each vertex patch. Appendix A complements the discussion with an algorithmic two-step procedure for the computation of $Q_{p}-\nabla _{\textrm{pw}}R u_h$ in 2D. The efficiency of the averaging $|\!|\!|(1 - {\mathcal {A}})R u_h|\!|\!|_\textrm{pw}$ up to data oscillations follows in Sects. 4.3–4.4. Section 4.6 concludes with the proof of Theorem 3.

4.2 Construction of equilibrated flux

This subsection defines the post-processed $H(\textrm{div},\Omega )$-conforming equilibrated flux $Q_p\in RT_{k+p}({\mathcal {T}})$ that enters the GUB $\eta _{\textrm{eq},p}$ from (34) based on local patch-wise minimization problems in the spirit of [10, 30, 31].

Consider the shape-regular vertex-patch $\omega (z):=\textrm{int}(\bigcup \{T\in \mathcal {T}(z)\})$ covered by the neighbouring simplices $\mathcal {T}(z):=\{T\in \mathcal {T}: z\in T\}$ sharing a given vertex $z \in {\mathcal {V}}$ with the facet spider $\mathcal {F}(z):=\{F\in \mathcal {F}\,\ z\in F\}$. Recall the space of piecewise Raviart-Thomas functions $RT_{k}^{\textrm{pw}}(\mathcal {T})$ from Sect. 1.4 and define

$$\begin{aligned} L^2_0(\omega (z))&:=\{f\in L^2(\omega (z)) : (f,1)_{L^2(\omega (z))}= 0\},\\ L^2_*(\omega (z))&:={\left\{ \begin{array}{ll} L^2_0(\omega (z))&{} \text { if } z \in {\mathcal {V}}(\Omega ),\\ L^2(\omega (z)) &{}\text { else}, \end{array}\right. }\\ H^1_*(\omega (z))&:={\left\{ \begin{array}{ll} H^1(\omega (z)) \cap L^2_0(\omega (z))\quad \text { if } z \in {\mathcal {V}}(\Omega ),\\ \{v \in H^1(\omega (z)) : v = 0 \text { on } \partial \Omega \cap \bigcup \mathcal {F}(z)\}\quad \text { else}, \end{array}\right. }\\ H_0(\mathop {\textrm{div}}\limits ,\omega (z))&:={\left\{ \begin{array}{ll} \{r \in H(\textrm{div},\omega (z)) : r \cdot \nu = 0 \text { on } \partial \omega (z)\} \quad \text { if } z \in {\mathcal {V}}(\Omega ),\\ \{r \in H(\textrm{div},\omega (z)) : r \cdot \nu = 0 \text { on } \partial \omega (z) \setminus \bigcup \mathcal {F}(z)\}\quad \text { else}, \end{array}\right. }\\ RT_k^0(\mathcal {T}(z))&:=\;RT_{k}^{\textrm{pw}}(\mathcal {T}(z)) \cap H_0(\mathop {\textrm{div}}\limits ,\omega (z)). \end{aligned}$$

Throughout the remaining parts of this section, abbreviate $G_h:=\nabla _\textrm{pw}R u_h \in P_k({\mathcal {T}};\mathbb R^n)$. Given a vertex $z\in {\mathcal {V}}$ with the $P_1$-conforming nodal basis function $\varphi _z\in S^1(\mathcal {T})$, the property (20) provides compatible data

$$\begin{aligned} f_z :=\, {\left\{ \begin{array}{ll} \Pi _{p}(\varphi _z \Pi _0 f- G_h\cdot \nabla \varphi _z) &{}\text{ if } k = 0,\\ \Pi _{k+p}(\varphi _z f- G_h\cdot \nabla \varphi _z) &{}\text{ if } k \ge 1 \end{array}\right. }\quad \in L^2_*(\omega (z)) \end{aligned}$$

(39)

such that the discrete affine space

$$\begin{aligned} {\mathcal {Q}}_h(z) :=\{\tau _z \in RT^0_{k+p}(\mathcal {T}(z))\ :\ \textrm{div}\, \tau _z + f_z = 0 \text { in } \Omega \} \ne \emptyset \end{aligned}$$

(40)

is not empty. Consequently,

$$\begin{aligned} Q_{z,h}:=\mathop {\textrm{arg}\,\textrm{min}}\limits _{\tau _z \in {\mathcal {Q}}_h(z)} \Vert \tau _z - {\mathcal {I}}_{RT} (\varphi _z G_h)\Vert _{L^2(\omega (z))}= \Pi _{{\mathcal {Q}}_h(z)} {\mathcal {I}}_{RT} (\varphi _z G_h) \end{aligned}$$

(41)

is well defined as the $L^2$ projection $\Pi _{{\mathcal {Q}}_h(z)} {\mathcal {I}}_{RT} (\varphi _z G_h)$ of ${\mathcal {I}}_{RT} (\varphi _z G_h)$ onto ${\mathcal {Q}}_h(z)$ with the piecewise Raviart-Thomas interpolation ${\mathcal {I}}_{RT}: H^1(\mathcal {T};\mathbb R^n)\rightarrow RT^\textrm{pw}_{k+p}(\mathcal {T})$ [7, Section III.3.1]. In the case $p \ge 1$, $\varphi _z G_h\in P_{k+1}({\mathcal {T}}(z)) \subset RT_{k+p}^\textrm{pw}({\mathcal {T}}(z))$ is a piecewise Raviart-Thomas function of degree $k+p$. Hence ${\mathcal {I}}_{RT}(\varphi _z G_h) = \varphi _z G_h$ and ${\mathcal {I}}_{RT}$ could be omitted in the formula (41). The partition of unity $\sum _{z\in {\mathcal {V}}}\varphi _z\equiv 1$ and $G_h= {\mathcal {I}}_{RT} G_h= \sum _{z \in {\mathcal {V}}} {\mathcal {I}}_{RT} (\varphi _z G_h)$ show that the sum $Q_p = \sum _{z \in {\mathcal {V}}} Q_{z,h}\in H(\mathop {\textrm{div}}\limits , \Omega )$ of the patch-wise contributions satisfies

$$\begin{aligned}&\mathop {\textrm{div}}\limits Q_p = {\left\{ \begin{array}{ll} \sum \limits _{z \in {\mathcal {V}}} \Pi _{p} (G_h\cdot \nabla \varphi _z - \varphi _z \Pi _0 f) = -\Pi _{0} f &{}\text{ if } k = 0,\\ \sum \limits _{z \in {\mathcal {V}}} \Pi _{k+p} (G_h\cdot \nabla \varphi _z - \varphi _z f) = -\Pi _{k+p} f &{}\text{ if } k \ge 1, \end{array}\right. } \end{aligned}$$

(42)

$$\begin{aligned}&\Vert Q_{p}-G_h\Vert _{L^2(\Omega )}^2 \lesssim \sum _{z\in {\mathcal {V}}} \Vert Q_{z,h}- {\mathcal {I}}_{RT}(\varphi _z G_h) \Vert _{L^2(\omega (z))}^2. \end{aligned}$$

(43)

This establishes the flux reconstruction $Q_p$. The efficiency of the flux reconstruction will be based on the following general equivalence.

Lemma 9

(control of $H(\mathop {\textrm{div}}\limits )$ minimization by residual [10, 31]) Given any vertex $z\in {\mathcal {V}}$, a piecewise Raviart-Thomas function $\sigma _z\in RT_{q}^{\textrm{pw}}(\mathcal {T}(z))$ and a piecewise polynomial $r_z \in {P}_{q}\left( {\mathcal {T}}(z)\right)$ of degree $q \in {\mathbb {N}}_0$, define the residual

$$\begin{aligned} \textrm{Res}_z(v):=\sum _{T \in {\mathcal {T}}(z)}\Big ( \left( r_{z}, v\right) _{L^2(T)}+\left\langle \sigma _z\cdot \nu _T, v\right\rangle _{L^2(\partial T)}\Big ) \end{aligned}$$

(44)

for all $v\in H^1(\Omega )$. If $z\in {\mathcal {V}}(\Omega )$ is an interior vertex, then suppose additionally that $\textrm{Res}_z(1) = 0$. Then

$$\begin{aligned} \min \limits _{\begin{array}{c} \tau _z \in RT^0_q({\mathcal {T}}(z))\\ \textrm{div} \tau _z = r_z + \textrm{div}_\textrm{pw}\sigma _z \end{array}} \left\| \tau _z-\sigma _z\right\| _{L^2(\omega (z))} \le C_{\textrm{s}}\max \limits _{ \begin{array}{c} v \in H^{1}_*(\omega (z)) \\ \Vert \nabla v\Vert _{L^2(\omega (z))}=1 \end{array} } \textrm{Res}_z(v) \end{aligned}$$

(45)

holds for a constant $C_{\textrm{s}}$ that exclusively depends on the shape-regularity (and is in particular independent of the polynomial degree q).

Proof

The assertion follows from [10, Theorem 7] in $n=2$ dimensions and [31, Corollaries 3.3, 3.6, and 3.8] in $n=3$ dimensions.$\square$

Remark 1

The patch-wise construction of $Q_p$ in Sect. 4.2 typically generates local data oscillation $\textrm{osc}_{k+p}(\varphi _z f,{\mathcal {T}}(z))$ in the error analysis as in the proof of Theorem 3 in Sect. 4.6 below or, e.g., [30, Theorem 3.17]. A straightforward computation $\textrm{osc}_{k+p}(\varphi _z f,{\mathcal {T}}(z)) \le \textrm{osc}_{k+p-1}(f,{\mathcal {T}}(z))$ apparently leads to a loss of one degree in the data oscillation but Lemma 8 verifies

$$\begin{aligned} \textrm{osc}_{k+p}^2(\varphi _z f, \mathcal T(z))\lesssim \min _{v_{h}\in P_{k+1}(\mathcal {T}(z))} \Vert \nabla _{\textrm{pw}}(u - v_{h})\Vert _{L^2(\omega (z)}^2 + \textrm{osc}_{q}^2(f,\mathcal T(z)) \end{aligned}$$

(46)

for any $p, q\in {\mathbb N}_0$. This allows for efficiency of the data oscillations on the right-hand side of the efficiency estimate [30, Formula (3.42)] and leads to a corresponding refinement in [30, Theorem 3.17].

4.3 Local equivalence of stabilizations

The first improvement to the current HHO literature is the local equivalence of the two stabilizations ${\tilde{s}}_h$ from (37) and $s_h$ from (2). The authors of this paper could not find any motivation for the alternative stabilization ${\tilde{s}}_h$ in the error analysis of [25, Section 4] and suggest to apply Theorem 4 below to [25, Theorem 4.7] to recover the results therein for the original HHO stabilization $s_h$. Recall the local stabilization ${\tilde{s}}_T$ in ${\tilde{s}}_h(v_h,v_h) :=\sum _{T \in \mathcal {T}} {\tilde{s}}_T(v_h,v_h)$ from (37) and $S_{TF}v_h = \Pi _{F,k} \left( v_{\mathcal {T}} + (1-\Pi _{T,k})R v_{h} \right) |_{T}-v_{\mathcal {F}}|_{F}$ in the definition of $s_h$ from (2).

Theorem 4

(local equivalence of stabilizations) Any $v_h = (v_{\mathcal T},v_{\mathcal {F}}) \in V_h$ and $T \in \mathcal T$ satisfy

$$\begin{aligned} {C_{10}}^{-1}{\tilde{s}}_T(v_h,v_h) \le \sum _{F \in {\mathcal {F}}(T)} h_F^{-1} \Vert S_{TF}v_h\Vert _{L^2(F)}^2 \le {C_{11}}{\tilde{s}}_T(v_h,v_h). \end{aligned}$$

(47)

The constants $C_10$ and $C_11$ exclusively depend on the polynomial degree k and the shape regularity of ${\mathcal {T}}$.

Proof

The second inequality in (47) follows directly from a triangle inequality and an inverse estimate. Therefore, the proof focuses on the first inequality in (47). Given $v_h = (v_\mathcal T,v_{\mathcal {F}}) \in V_h$ and $T \in \mathcal T$, let $\varphi _k :=(\Pi _k R v_h - v_\mathcal {T})|_{T} \in P_k(T)$. Since $S_{TF}v_h = \Pi _{F,k}(R v_{h|T} - v_\mathcal {F}|_F - \varphi _k)$, the triangle inequality $\Vert \Pi _{F,k}(R v_{h|T} - v_{\mathcal {F}})\Vert _{L^2(F)} \le \Vert S_{TF}v_h\Vert _{L^2(F)} + \Vert \varphi _k\Vert _{L^2(F)}$, the discrete trace inequality $\Vert \varphi _k\Vert _{L^2(F)} \lesssim h_F^{-1/2}\Vert \varphi _k\Vert _{L^2(T)}$, and the shape-regularity $h_F \approx h_T$ for all $F \in {\mathcal {F}}(T)$ reveal

$$\begin{aligned}&\sum _{F \in {\mathcal {F}}(T)} h_F^{-1}\Vert \Pi _{F,k}(R v_{h|T} - v_{\mathcal {F}})\Vert _{L^2(F)}^2\nonumber \\&\qquad \lesssim \sum _{F \in {\mathcal {F}}(T)} h_F^{-1} \Vert S_{TF}v_h\Vert _{L^2(F)}^2 + h_T^{-2}\Vert \varphi _k\Vert _{L^2(T)}^2. \end{aligned}$$

(48)

Since $\Pi _0 \varphi _k = 0$ (from the design of $Rv_h$), a Poincaré inequality shows

$$\begin{aligned} h_T^{-2}\Vert \varphi _k\Vert _{L^2(T)}^2 \le C_P^2\Vert \nabla \varphi _k\Vert _{L^2(T)}^2. \end{aligned}$$

(49)

On the one hand, an integration by parts provides

$$\begin{aligned} \Vert \nabla \varphi _k\Vert _{L^2(T)}^2 = - (\Pi _k R v_h - v_\mathcal {T}, \Delta \varphi _k)_{L^2(T)} + \langle \varphi _k, \nabla \varphi _k \cdot \nu _T\rangle _{L^2(\partial T)}. \end{aligned}$$

(50)

On the other hand, an integration by parts and the definition of R imply

$$\begin{aligned}&- (R v_h, \Delta \varphi _k)_{L^2(T)} = (\nabla R v_h, \nabla \varphi _k)_{L^2(T)} - \langle R v_{h|T}, \nabla \varphi _k \cdot \nu _T\rangle _{L^2(\partial T)}\nonumber \\&\qquad = - (v_{\mathcal {T}}, \Delta \varphi _k)_{L^2(T)} + \sum _{F \in {\mathcal {F}}(T)} \langle v_{\mathcal {F}} - R v_{h|T}, \nabla \varphi _k \cdot \nu _T\rangle _{L^2(F)}. \end{aligned}$$

(51)

Since $\Delta \varphi _k \in P_k(T)$, the $L^2$ projection $\Pi _k$ on the right-hand side of (50) is redundant. Hence, the combination of (50)–(51) with $\nabla \varphi _k \cdot \nu _{T|F} \in P_k(F)$ for all $F \in {\mathcal {F}}(T)$ results in

$$\begin{aligned} \Vert \nabla \varphi _k\Vert _{L^2(T)}^2 = \sum _{F \in {\mathcal {F}}(T)} \langle \Pi _{F,k}(v_{\mathcal {F}} - R v_{h|T} + \varphi _k), \nabla \varphi _k \cdot \nu _T \rangle _{L^2(F)}. \end{aligned}$$

(52)

A Cauchy inequality on the right-hand side of (52), a discrete trace inequality, and $S_{TF}v_h = \Pi _{F,k}(R v_{h|T} - v_{\mathcal {F}|F} - \varphi _k)$ for all $F \in {\mathcal {F}}(T)$ lead to

$$\begin{aligned} \Vert \nabla \varphi _k\Vert _{L^2(T)}^2 \lesssim \sum _{F \in {\mathcal {F}}(T)} h_F^{-1} \Vert S_{TF}v_h\Vert _{L^2(F)}^2. \end{aligned}$$

(53)

Since ${\tilde{s}}_T(v_h,v_h) = \sum _{F \in {\mathcal {F}}(T)} h_F^{-1}\Vert \Pi _{F,k}(R v_{h|T} - v_{\mathcal {F}})\Vert _{L^2(F)}^2 + h_T^{-2}\Vert \varphi _k\Vert _{L^2(T)}^2$, the combination of (48)–(49) with (53) concludes the proof of (47).$\square$

4.4 Efficiency of the stabilization

The second improvement to the HHO literature is a quasi-best approximation estimate along the lines of the seminal paper [32, Theorem 4.10]. In combination with Theorem 4, this, in particular, provides the efficiency (38) of the stabilization up to data oscillation.

Theorem 5

(quasi-best approximation up to data oscillation) For any $p\in {\mathbb {N}}_0$, the solution u to (3) and the discrete solution $u_h$ to (18) satisfy

$$\begin{aligned}&|\!|\!|u - R u_h|\!|\!|_\textrm{pw}+ s_h(u_h,u_h)^{1/2}\\&\qquad \le {C_{12}}\Big (\min _{v_{k+1} \in P_{k+1}(\mathcal T)} |\!|\!|u - v_{k+1}|\!|\!|_\textrm{pw}+ \textrm{osc}_{k+p}(f,\mathcal T)\Big ). \end{aligned}$$

The constant $C_12$ exclusively depends on k, p, and the shape regularity of ${\mathcal {T}}$.

Proof

Given $k,p \in {\mathbb {N}}_0$, let ${\widetilde{u}} \in V$ solve the Poisson model problem $-\Delta {\widetilde{u}} = \Pi _{k+p} f$ with the right-hand side $\Pi _{k+p} f$. Section 4.3 in [32] constructs a stable enriching operator ${\mathcal {J}}: V_h \rightarrow V$ with local bubble-functions from [38]. A modification, where the polynomial degree $k-1$ in [32, Eq. (4.16)] is replaced by $k+p$, leads to a right-inverse ${\mathcal {J}}: V_h \rightarrow V$ of the interpolation $\textrm{I}: V_h \rightarrow V$ with the stability $|\!|\!|{\mathcal {J}} v_h|\!|\!|^2 \lesssim a_h(v_h,v_h)$ and the additional $L^2$ orthogonality ${\mathcal {J}} v_h - v_{\mathcal {T}} \perp P_{k+p}({\mathcal {T}})$ for all $v_h = (v_{\mathcal {T}},v_{\mathcal {F}}) \in V_h$. The extra orthogonality allows for

$$\begin{aligned} (\Pi _{k+p} f, {\mathcal {J}} v_h)_{L^2(\Omega )} = (\Pi _{k+p} f, v_{\mathcal {T}})_{L^2(\Omega )} = (f,v_{\mathcal {T}})_{L^2(\Omega )}\quad \text {for all } (v_\mathcal {T}, v_\mathcal {F})\in V_h. \end{aligned}$$

Consequently, the smoother ${\mathcal {J}}$ leads to a discrete solution $u_h=(u_\mathcal {T},u_\mathcal {F})\in V_h$ in the modified HHO discretization of [32] as a quasi-best approximation of the above solution ${{\tilde{u}}}$. The point is that $u_h\in V_h$ coincides with the original HHO solution $u_h$ from (18). The arguments from the proof of [32, Theorem 4.10] reveal the quasi-best approximation

$$\begin{aligned} |\!|\!|{\tilde{u}} - R u_h|\!|\!|_\textrm{pw}+ s_h(u_h,u_h)^{1/2} \lesssim \min _{v_{k+1} \in P_{k+1}(\mathcal {T})} |\!|\!|{\tilde{u}} - v_{k+1}|\!|\!|_\textrm{pw}\end{aligned}$$

also for the above modified smoother ${\mathcal {J}}$. This, the triangle inequalities $|\!|\!|u - R u_h|\!|\!|_\textrm{pw}\le |\!|\!|u - {\tilde{u}}|\!|\!| + |\!|\!|{\tilde{u}} - R u_h|\!|\!|_\textrm{pw}$ and $|\!|\!|{\tilde{u}} - v_{k+1}|\!|\!|_\textrm{pw}\le |\!|\!|u - {\tilde{u}}|\!|\!| + |\!|\!|u - v_{k+1}|\!|\!|_\textrm{pw}$ for any $v_{k+1} \in P_{k+1}({\mathcal {T}})$, and the standard perturbation bound $|\!|\!|u - {\tilde{u}}|\!|\!| \le C_P\textrm{osc}_{k+p}(f,\mathcal {T})$ conclude the proof. $\square$

4.5 Stabilization-free efficiency of averaging

The main result of this subsection establishes the stabilization-free efficiency of the nodal averaging technique.

Theorem 6

(averaging is efficient) Let $u \in V$ resp. $u_h \in V_h$ solve (3) resp. (18). Then $R u_h$ and ${\mathcal {A}} R u_h$ satisfy, for any $p \in {\mathbb {N}}_0$, that

$$\begin{aligned} {C_{13}}^{-1}|\!|\!|(1 - {\mathcal {A}}) R u_h|\!|\!|_\textrm{pw}&\le |\!|\!|u - R u_h|\!|\!|_\textrm{pw}+ \textrm{osc}_{k+p}(f,\mathcal T). \end{aligned}$$

(54)

The constant $C_12$ exclusively depends on k, p, and the shape regularity of the triangulation ${\mathcal {T}}$.

The proof can follow the proof of [25, Theorem 4.7] but additionally utilizes the two significant improvements from Sects. 4.3–4.4 that allow a stabilization-free efficiency in (54).

Proof

Theorem 4.7 in [25] shows $|\!|\!|(1 - {\mathcal {A}})R u_h|\!|\!|_\textrm{pw}^2 \lesssim |\!|\!|u - R u_h|\!|\!|_\textrm{pw}^2 + {\tilde{s}}_h(u_h,u_h)$. This, the equivalence ${\tilde{s}}_h(u_h,u_h) \approx s_h(u_h,u_h)$ of stabilizations (from Theorem 4), and the efficiency $s_h(u_h,u_h) \lesssim |\!|\!|u - R u_h|\!|\!|_\textrm{pw}^2 + \textrm{osc}_{k+p}^2(f,{\mathcal {T}})$ (from Theorem 5) imply the assertion. $\square$

Remark 2

(p-robustness) The $H^1(\Omega )$-conforming post-processing of $R u_h$ from [31] provides an efficiency constant independent of the polynomial degree k. The right-hand side of [31, Corollary 4.2] involves the stabilization-related term $\sum _{F \in \mathcal {F}}h_F^{-1}\Vert \Pi _{F,0}[R u_h]_F\Vert _{L^2(F)}^2$. It remains an open question whether this term can be controlled p-robustly by $|\!|\!|u - R u_h|\!|\!|_\textrm{pw}+ \textrm{osc}_{k}(f,\mathcal {T})$ (with a multiplicative constant independent of the polynomial degree k).

4.6 Proof of Theorem 3

Let $p \in {\mathbb {N}}_0$ and $r = 0$ if $k = 0$ and $r = k+p$ if $k \ge 1$ as in (35) be given. Recall the abbreviation $G_h= \nabla _\textrm{pw}R u_h \in P_k({\mathcal {T}})$ with the discrete solution $u_h \in V_h$ to (18) and let $Q_p = \sum _{z \in {\mathcal {V}}} Q_{z,h}\in H(\mathop {\textrm{div}}\limits , \Omega )$ denote the flux reconstruction from Sect. 4.2. The proof establishes (36) in five steps.

Step 1 provides the GUB $|\!|\!|u - R u_h|\!|\!|_\textrm{pw}\le \eta _{\textrm{eq},p}({\mathcal {T}})$. This can follow from the paradigm of [1, 2, 30] as outlined below. The choice ${G}:=G_h$ and $w :={\mathcal {A}} R u_h$ in (8) and a triangle inequality lead to

$$\begin{aligned} |\!|\!|u - R u_h|\!|\!|_\textrm{pw}^2 \le (|\!|\!|f + \mathop {\textrm{div}}\limits Q_p|\!|\!|_* + |\!|\!|\textrm{div}(Q_p - G_h)|\!|\!|_*)^2 + |\!|\!|(1 - {\mathcal {A}}) R u_h|\!|\!|_\textrm{pw}^2. \end{aligned}$$

Since $\mathop {\textrm{div}}\limits Q_p + \Pi _{r} f = 0$ vanishes in $\Omega$ by (42), a piecewise Poincaré inequality shows $|\!|\!|f + \mathop {\textrm{div}}\limits Q_p|\!|\!|_* \le C_P\textrm{osc}_{r}(f,\mathcal {T})$. This, the bound $|\!|\!|\textrm{div}(Q_p - G_h)|\!|\!|_* \le \Vert Q_p - G_h\Vert$ from the definition of $|\!|\!|\bullet |\!|\!|_*$, and the previously displayed formula result in the reliability $|\!|\!|u - R u_h|\!|\!|_\textrm{pw}\le \eta _{\textrm{eq,p}}({\mathcal {T}})$.

Step 2 establishes $\textrm{osc}_{k-1}(f,{\mathcal {T}}) \lesssim \eta _{\textrm{eq},p}({\mathcal {T}})$. Lemma 8 provides

$$\begin{aligned} \textrm{osc}_{k-1}(f,{\mathcal {T}})\lesssim |\!|\!|u - R u_h|\!|\!|_\textrm{pw}+\textrm{osc}_{r}(f, \mathcal {T}). \end{aligned}$$

(55)

This, Step 1 and $\textrm{osc}_{r}(f, \mathcal {T})\lesssim \eta _{\textrm{eq}, p}(\mathcal {T})$ from (34) conclude the proof of Step 2.

Step 3 reveals, for any $q\in {\mathbb N}_0$, the efficiency of the flux reconstruction

$$\begin{aligned} \Vert Q_p - G_h\Vert \lesssim |\!|\!|u-R u_h|\!|\!|_\textrm{pw}+ \textrm{osc}_{q}(\mathcal T, f) \end{aligned}$$

(56)

for any polynomial degree $k \ge 1$. The case $k=0$ follows in Step 4 below. Recall $f_z = \Pi _{k+p}(\varphi _z f - G_h\cdot \nabla \varphi _z)$ from (39) for any vertex $z\in {\mathcal {V}}$ in the construction of $Q_p$ from Sect. 4.2 and set $\sigma _z :={\mathcal {I}}_{\textrm{RT}}(\varphi _z G_h)\in RT_{k+p}^{\textrm{pw}}(\mathcal {T}(z))$. The piecewise Raviart-Thomas interpolation ${\mathcal {I}}_{\textrm{RT}}:H^1(\mathcal {T};\mathbb R^n)\rightarrow RT_{k+p}^\textrm{pw}(\mathcal {T})$ satisfies the well-known commuting diagram properties [7, Section 2.5.1]

$$\begin{aligned} {\mathop {\textrm{div}}\limits }_{\textrm{pw}}\circ {\mathcal {I}}_{\textrm{RT}}= \Pi _{k+p}\circ {\mathop {\textrm{div}}\limits }_{\textrm{pw}}{} & {} \text {and}{} & {} \gamma _T\circ {\mathcal {I}}_{\textrm{RT}}|_F=\Pi _{F,k+p}\circ \gamma _T \end{aligned}$$

(57)

for any facet $F\in \mathcal {F}(T)$ of a simplex $T\in \mathcal {T}$ and the normal trace $\gamma _T:H^1(T;\mathbb R^n)\rightarrow L^2(\partial T)$ with $\gamma _T \sigma :=\sigma \cdot \nu _T$ for $\sigma \in H^1(T;\mathbb R^n)$. This and elementary algebra with the product rule ${\mathop {\textrm{div}}\limits }_{\textrm{pw}}(\varphi _z G_h)=\varphi _z{\mathop {\textrm{div}}\limits }_{\textrm{pw}}G_h+ \nabla \varphi _z\cdot G_h\in P_k({\mathcal {T}}(z))$ imply

$$\begin{aligned} r_z :=-\textrm{div}_\textrm{pw}\sigma _z - f_z = - \varphi _z \textrm{div}_\textrm{pw}G_h- \Pi _{k+p}(\varphi _zf)\in P_{k+p}(\mathcal {T}(z)). \end{aligned}$$

(58)

Recall the residual $\textrm{Res}_z(v)$ with $v\in H^1(\Omega )$ for the vertex $z\in {\mathcal {V}}$ from (44). The commuting diagram property (57)–(58) establish the identity

$$\begin{aligned} \textrm{Res}_z(1) = -(\textrm{div}_\textrm{pw}G_h+ f, \varphi _z)_{L^2(\omega (z))} + \sum _{T \in {\mathcal {T}}(z)} \langle G_h\cdot \nu _T, \varphi _z\rangle _{L^2(\partial T)}. \end{aligned}$$

This, a piecewise integration by parts, and the property (20) verify $\textrm{Res}(1) = (G_h,\nabla \varphi _z)_{L^2(\omega (z))} - (f,\varphi _z)_{L^2(\omega (z))}=0$ for any interior vertex $z \in {\mathcal {V}}(\Omega )$. Hence, Lemma 9 applies for any vertex $z\in {\mathcal {V}}$ and provides

$$\begin{aligned} \Vert Q_{z,h}- {\mathcal {I}}_{RT}(\varphi _z G_h) \Vert _{L^2(\omega (z))} \le C_{\textrm{s}}\sup \limits _{ \begin{array}{c} v \in H^{1}_*(\omega (z))\\ \Vert \nabla v\Vert _{L^2(\omega (z))}=1 \end{array}} \textrm{Res}_z(v) \end{aligned}$$

(59)

for the local contributions $Q_{z,h}$ of $Q_p = \sum _{z \in {\mathcal {V}}} Q_{z,h}$ from (41). The identity ${\mathcal {I}}_{\textrm{RT}}(\varphi _z G_h) = \varphi _z G_h$ for $p\ge 1$ from $\varphi _zG_h\in P_{k+1}(\mathcal {T}(z))\subset RT_{k+p}^{\textrm{pw}}(\mathcal {T}(z))$ allows for a k- and p-robust efficiency control of the flux reconstruction error

$$\begin{aligned} {C_{14}}^{-1} \Vert Q_p - G_h\Vert ^2 \le \Vert \nabla u-G_h\Vert ^2 + \sum _{z \in {\mathcal {V}}} \textrm{osc}_{r}^2(\varphi _z f,\mathcal T(z)) \end{aligned}$$

(60)

with a constant $C_14$ that solely depends on the shape regularity of ${\mathcal {T}}$. This is deemed noteworthy and motivates two different approaches for the bound of the residual $\textrm{Res}_z(v)$ on the right-hand side of (59) for $p=0$ and $p \ge 1$ below.

Step 3.1 provides (56) for $p=0$. Given any normalized $v \in H^1_*(\omega (z))$ with $\Vert \nabla v\Vert _{L^2(\omega (z))}=1$, the product ${\mathcal {I}}_{RT}(\varphi _z G_h) \cdot \nu _F\,v$ vanishes on any boundary facet $F \in {\mathcal {F}}(\partial \omega (z))$ of the patch $\omega (z)$. This, the commuting diagram property (57), and (58) with $\varphi _z{\mathop {\textrm{div}}\limits }_{\textrm{pw}}G_h\in P_k({\mathcal {T}}(z))$ in the definition of the residual $\textrm{Res}_z(v)$ from (44) verify

$$\begin{aligned} \textrm{Res}_z(v)&= -(\varphi _z f + \varphi _z{\mathop {\textrm{div}}\limits }_{\textrm{pw}}G_h,\Pi _k v)_{L^2(\omega (z))}\nonumber \\&\quad + \sum _{F\in \mathcal {F}(z) \cap \mathcal {F}(\Omega )} \langle \varphi _z [G_h]_F\cdot \nu _F, \Pi _{F,k} v\rangle _{L^2(F)}. \end{aligned}$$

(61)

The shape regularity of ${\mathcal {T}}$, a Poincaré inequality for interior vertices $z \in {\mathcal {V}}(\Omega )$, and a Friechrichs inequality for boundary vertices $z \in {\mathcal {V}}(\partial \Omega )$ provide

$$\begin{aligned} \Vert h_{\mathcal {T}}^{-1} v\Vert _{L^2(\omega (z))} \approx \textrm{diam}(\omega (z))^{-1} \Vert v\Vert _{L^2(\omega (z))} \lesssim \Vert \nabla v\Vert _{L^2(\omega (z))} = 1. \end{aligned}$$

(62)

Given any facet $F \in {\mathcal {F}}(z) \cap {\mathcal {F}}(\Omega )$ in the facet spider $F\in \mathcal {F}(z)$ with facet patch $\omega (F)$, a trace inequality thus shows

$$\begin{aligned} h_F^{-1/2}\Vert v\Vert _{L^2(F)} \lesssim \Vert h_{\mathcal {T}}^{-1} v\Vert _{L^2(\omega (F))}+\Vert \nabla v\Vert _{L^2(\omega (F))}\lesssim 1. \end{aligned}$$

(63)

Cauchy inequalities on the right-hand side of (61), the stability of the $L^2$ projection, $\Vert \varphi _z\Vert _{L^\infty (\omega _z)} = 1$, and (62)–(63) prove that ${\textrm{Res}}_z(v)$ is controlled by

$$\begin{aligned} \Vert h_{\mathcal {T}}(f + {\mathop {\textrm{div}}\limits }_{\textrm{pw}}G_h)\Vert _{L^2(\omega (z))}+ \left( \sum _{F \in {\mathcal {F}}(z) \cap {\mathcal {F}}(\Omega )} h_F\Vert [G_h]_F \cdot \nu _F\Vert _{L^2(F)}^2\right) ^{1/2}. \end{aligned}$$

The efficiency of those residual terms follows from the proof of Theorem 2 in Sect. 3. In combination with (43) and (59), this results in the global efficiency (56) for $p = 0$ and concludes the proof of Step 3.1.

Step 3.2 provides (56) for $p\ge 1$. Given any normalized $v \in H^1_*(\omega (z))$ with $\Vert \nabla v\Vert _{L^2(\omega (z))}=1$, (58) and the identity ${\mathcal {I}}_{\textrm{RT}}(\varphi _z G_h) = \varphi _z G_h$ in the definition of the residual $\textrm{Res}_z(v)$ from (44) verify that $\textrm{Res}_z(v)$ is equal to

$$\begin{aligned} - (\varphi _z \textrm{div}_\textrm{pw}G_h+ \Pi _{k+p}(\varphi _z f) , v)_{L^2(\omega (z))} + \sum _{T \in {\mathcal {T}}(z)} \langle G_h\cdot \nu _T, \varphi _z v \rangle _{L^2(\partial T)}. \end{aligned}$$

This, a piecewise integration by parts, and the product rule $\nabla (\varphi _z v)=\varphi _z\nabla v+ v\nabla \varphi _z$ reveal

$$\begin{aligned} \textrm{Res}_z(v)&= - (\Pi _{k+p}(\varphi _z f) , v)_{L^2(\omega (z))}+ ( G_h, \nabla (\varphi _z v))_{L^2(\omega (z))}. \end{aligned}$$

(64)

The weak formulation (3) with the test function $\varphi _z v \in H^1_0(\omega (z))\subset V$ shows $(\nabla u, \nabla (\varphi _z v))_{L^2(\omega (z))} = (f, \varphi _z v)_{L^2(\omega (z))}$. Consequently, (64) implies

$$\begin{aligned} \textrm{Res}_z(v) =&\, ((1 - \Pi _{k+p})(\varphi _z f), v)_{L^2(\omega (z))}- (\nabla u - G_h, \nabla (\varphi _z v))_{L^2(\omega (z))}. \end{aligned}$$

This, a Cauchy-Schwarz, and a piecewise Poincaré inequality with the normalization $\Vert \nabla v\Vert _{L^2(\omega (z))} = 1$ provide

$$\begin{aligned} \textrm{Res}_z(v)&\le \Vert \nabla u - G_h\Vert _{L^2(\omega (z))}\Vert \nabla (\varphi _zv)\Vert _{L^2(\omega (z))} + C_P\textrm{osc}_{k+p}(\varphi _z f,\mathcal {T}(z)). \end{aligned}$$

(65)

Since $\Vert \nabla \varphi _z\Vert _{L^\infty (\omega (z))}\approx h^{-1}$, the Leibniz rule, a triangle inequality, and (62) show that

$$\begin{aligned} \Vert \nabla (\varphi _zv)\Vert _{L^2(\omega (z))}&\le \Vert \varphi _z\Vert _{L^\infty (\omega (z))}\Vert \nabla v\Vert _{L^2(\omega (z))}+\Vert \nabla \varphi _z\Vert _{L^\infty (\omega (z))}\Vert v\Vert _{L^2(\omega (z))}\\ {}&\lesssim \Vert \nabla v\Vert _{L^2(\omega (z))} = 1 \end{aligned}$$

can be bounded by a constant independent of h. Thus, the combination of (43) with (59) and (65) results in the efficiency of $\Vert Q_p - G_h\Vert$ in (60) with a k- and p-robust constant $C_{14}$. Since $\varphi _z \Pi _{k-1} f \in P_k(T) \subset P_{k+p}(T)$, the Pythagoras theorem and $\Vert \varphi _z\Vert _{L^\infty (\omega (z))} = 1$ show $\Vert (1 - \Pi _{k+p})(\varphi _z f)\Vert _{L^2(T)} \le \Vert \varphi _z(1 - \Pi _{k-1}) f\Vert _{L^2(T)} \le \Vert (1-\Pi _{k-1})f\Vert _{L^2(T)}$ for all $T \in {\mathcal {T}}(z)$, whence

$$\begin{aligned} \textrm{osc}_{k+p}(\varphi _z f,{\mathcal {T}}(z)) \le \textrm{osc}_{k-1}(f,{\mathcal {T}}(z)). \end{aligned}$$

This and Lemma 8 implies the efficiency

$$\begin{aligned} \textrm{osc}_{k+p}(\varphi _z f, \mathcal T(z)) \lesssim \Vert \nabla _{\textrm{pw}}(u - R u_h)\Vert _{L^2(\omega (z))} + \textrm{osc}_{q}(f,\mathcal T(z)) \end{aligned}$$

of the local oscillations from (65) for any $q\in {\mathbb {N}}_0$ as in Remark 1. This and (60) prove (56) for $p \ge 1$ and conclude the proof of Step 3.2.

Step 4 affirms the efficiency of the flux reconstruction

$$\begin{aligned} {C_{15}}^{-1}\Vert Q_p - G_h\Vert \le |\!|\!|u - R u_h|\!|\!|_\textrm{pw}+ \textrm{osc}_{q}(f,{\mathcal {T}}) \end{aligned}$$

(66)

for the polynomial degree $k = 0$ and any $q\in {\mathbb {N}}_0$. Let ${\widetilde{u}} \in V$ solve the Poisson model problem $-\Delta {\widetilde{u}} = \Pi _0 f$ with piecewise constant right-hand side $\Pi _0 f \in P_0({\mathcal {T}})$. A careful inspection reveals that all arguments from Step 3 apply to the case $k = 0$ for f replaced by $\Pi _0 f$. This leads to $\Vert Q_p - G_h\Vert \le C_{15}|\!|\!|{\widetilde{u}} - R u_h|\!|\!|_\textrm{pw}$ with a constant $C_{15}$ that solely depends on the shape of ${\mathcal {T}}$ (because the data oscillations on the right-hand side of (60) vanish). Therefore, the triangle inequality $|\!|\!|{\widetilde{u}} - R u_h|\!|\!|_\textrm{pw}\le |\!|\!|u - {\widetilde{u}}|\!|\!| + |\!|\!|u - R u_h|\!|\!|_\textrm{pw}$, the standard bound $|\!|\!|u - {\widetilde{u}}|\!|\!| \le C_P\textrm{osc}_{0}(f,{\mathcal {T}})$, and $C_P = 1/\pi < 1$ on simplicial domains lead to (66). This and the efficiency of the data oscillations from Lemma 8 conclude the proof of Step 4. Notice that the constant ${C_{15}}$ is independent of the parameter p.

Step 5 finishes the proof. On the one hand, the reliability $|\!|\!|u - R u_h|\!|\!|_\textrm{pw}+ \textrm{osc}_{k-1}(f,{\mathcal {T}}) \lesssim \eta _{\textrm{eq},p}$ is established in Step 1–2. On the other hand, the efficiency of the averaging operator from Theorem 6, Lemma 8, and the efficiency of the flux reconstruction in (56) for $k \ge 1$ and in (66) for $k = 0$ imply the efficiency $\eta _{\textrm{eq,p}} \lesssim |\!|\!|u - R u_h|\!|\!|_\textrm{pw}+ \textrm{osc}_{q}(f,{\mathcal {T}})$ for any $q \in {\mathbb {N}}_0$. This concludes the proof of Theorem 3.

5 Numerical experiments

This section provides numerical evidence for optimal convergence and a comparison of the stabilization-free GUBs $\eta _{\textrm{res}}$ and $\eta _{\textrm{eq,p}}$ from the Sects. 3 and 4 with the original error estimator $\eta _{\textrm{HHO}}$ from [25, Theorem 4.3] for the HHO method in three 2D benchmarks.

5.1 A posteriori error estimation with explicit constants

All triangulations in this section consist of right-isosceles triangles with the Poincaré constant $C_P = (\sqrt{2}\pi )^{-1}$ [34]. With the estimates $C_\textrm{1}\le C_\mathcal {T}$, $C_{\textrm{H}}\le 1$, and $C_\textrm{2}\le C_\mathcal {F}$ from Example 1, the residual-based error estimator from (22) reads

$$\begin{aligned} \eta ^2_{\textrm{res}}(\mathcal {T}):=&\left( C_\mathcal {T}\eta _{\textrm{res}, 1}(\mathcal {T}) + C_P\eta _{\textrm{res}, 2}(\mathcal {T}) +C_\mathcal {F}\eta _{\textrm{res}, 3}(\mathcal {T})\right) ^2+ C_\mathcal {F}^2\eta ^2_{\textrm{res}, 4}(\mathcal {T}). \end{aligned}$$

The equilibrated GUB from Theorem 3,

$$\begin{aligned} \eta _{\textrm{eq},p}^2(\mathcal {T}):=\left( C_P\textrm{osc}_{k+p}(f, \mathcal {T}) + \Vert Q_p^\Delta \Vert _{L^2(\Omega )}\right) ^2 + |\!|\!|(1-{\mathcal {A}}) Ru_h|\!|\!|_\textrm{pw}^2, \end{aligned}$$

depends on the post-processed quantity $Q_p^\Delta :=Q_p-\nabla _{\textrm{pw}}Ru_h$ with the Raviart-Thomas function $Q_p\in RT_{k+p}(\mathcal {T})$ of degree $k+p$ for $p\in {\mathbb {N}}_0$ from Sect. 4.2. An algorithmic description of the computation of $Q_p^\Delta$ for arbitrary p follows in Appendix A. Theorem 3 shows that $\eta _{\textrm{eq}, p}$ is efficient and reliable for all $p\in {\mathbb {N}}_0$. The residual-based error estimator $\eta _{\textrm{HHO}}$ from [25, Theorem 4.3] reads

$$\begin{aligned} \eta _{\textrm{HHO}}^2(\mathcal {T}) :=&\sum _{T\in \mathcal {T}} \bigg ( C_Ph_{T}\Vert (I - \Pi _{T,0})(f+{\Delta }_{\textrm{pw}} Ru_h)\Vert _{L^2(T)}\\&\quad + \sqrt{\sum _{F\in \mathcal {F}(T)} C_{\partial T} h_T \Vert R_{T, F}^k u_h\Vert _{L^2(F)}^2}\bigg )^2 + |\!|\!|(I-{\mathcal {A}}) Ru_h|\!|\!|_\textrm{pw}^2 \end{aligned}$$

with the operator $R_{T, F}^k$ from [25, Eq. (2.59)] for the original HHO stabilization [25, Eq. (2.22)] that induces the global stabilization $s_h$. The constant $C_{\partial T}=12C_P(C_P+C_{\textrm{tr}})\le 2.0315$ for right-isosceles triangles bounds the trace of $f-\Pi _0 f$ for $f\in H^1(T)$ and improves on the estimate $C_{\partial T}\le (C_P(h_T|\partial T|/|T|)(1+C_P))^2=2.524$ from [25, Subs. 4.1.1].

Lemma 10

(Poincaré-type inequality on trace) Given a simplex $T \subset {\mathbb {R}}^n$, any $f\in H^1(T)$ and $C_{\partial T}:=(n+1)h_T^2|T|^{-1} C_P(C_P + 2C_{\textrm{tr}}/n)$ satisfy

$$\begin{aligned} \Vert f-\Pi _{0} f\Vert _{L^2(\partial T)}^2 \le C_{\partial T}h_T\Vert \nabla f\Vert _{L^2(T)}^2. \end{aligned}$$

(67)

Proof

Abbreviate ${\tilde{\ell }}(F):=(n+1)h_Fh_T^2|T|^{-1}$ for the facet $F\in \mathcal {F}(T)$ of T and apply Lemma 3 to the singleton triangulation $\{T\}$ and $f-\Pi _0f\in H^1(T)$. This and the Poincaré inequality $\Vert f-\Pi _{0, T}f\Vert _{L^2(T)}\le C_P h_T\Vert \nabla f\Vert _{L^2(T)}$ show

$$\begin{aligned} \sum _{F\in \mathcal {F}(T)} {\tilde{\ell }}(F)^{-1}\Vert f-\Pi _{0, T} f\Vert _{L^2(F)}^2&\le (C_P+2C_{\textrm{tr}}/n)C_P \Vert \nabla f\Vert _{L^2(T)}^2. \end{aligned}$$

The assertion (67) follows from the observation that ${\tilde{\ell }}(F)$ is maximized on the facet $F\in \mathcal {F}(T)$ with $h_F=h_T$.$\square$

5.2 Implementation and adaptive algorithm

Our implementation of the HHO method in MATLAB uses nodal bases for the spaces $P_k(\mathcal {F}), P_k(\mathcal {T})$, and $P_{k+1}(\mathcal {T})$ and the direct solver mldivide (behind the \-operator) for the discrete system of equations representing (18). For implementation details on the HHO method itself we refer to [25, Appendix B]. The integration of polynomial expressions is carried out exactly. The errors in approximating non-polynomial expressions, such as exact solutions u and source terms f, by polynomials of sufficiently high degree are expected to be very small and are neglected for simplicity.

Algorithm 1 displays the standard adaptive algorithm (AFEM) [16, 22] driven by the refinement indicators, for any triangle $T\in \mathcal {T}$,

$$\begin{aligned} \eta ^2_{\textrm{res}}(T)&:=|T|\Vert f+{\Delta }_{\textrm{pw}} Ru_h\Vert ^2_{L^2(T)} +|T|^{1/2}\sum _{F\in {\mathcal {F}}(T)}\Vert [ {\nabla }_{\textrm{pw}} Ru_h ]_F \Vert _{L^2(F)}^2 \end{aligned}$$

(68)

with the modified jump $[\bullet ]_F:=\bullet \times n_F$ along a boundary side $F\in \mathcal {F}(\partial \Omega )$ and the newest-vertex-bisection (NVB). The sum of this over all triangles is, up to some multiplicative constants, equivalent to $\eta _\textrm{res}^2({\mathcal {T}})$.

Algorithm 1

AFEM algorithm

5.3 High oscillations on the unit square

This benchmark on the unit square $\Omega :=(0,1)^2$ considers the Laplace equation $-\Delta u=f$ with source term f matching the smooth exact solution

$$\begin{aligned} u(x, y) = x(x-1)y(y-1) e^{-100((x-1/2)^2 + (y-117/1000)^2)}\in C^\infty (\Omega ). \end{aligned}$$

Figure 1 displays the energy norm $|\!|\!|e|\!|\!|_\textrm{pw}$ of the error $e:=u-Ru_h$ for uniform and adaptive refinement by Algorithm 1 on the left. The smooth solution allows for optimal convergence rates $(k+1)/2$ in the number ndof of degrees of freedom, while the adaptive mesh sequence leads to a lower energy error with respect to ndof. The GUBs $\eta _{\textrm{res}}, \eta _{\textrm{eq}, p}$, and $\eta _{\textrm{HHO}}$ are efficient and therefore equivalent to the energy error $|\!|\!|e|\!|\!|_\textrm{pw}$, see Fig. 1 on the right for $k=0, 2$.

Figure 2 shows the efficiency indices $EF(\eta ):=\eta /|\!|\!|e|\!|\!|_\textrm{pw}$ for the residual-based GUBs $\eta =\eta _{\textrm{res}}, \eta _{\textrm{HHO}}$ and the equilibration-based GUB $\eta =\eta _{\textrm{eq}, p}$ for $p=0,1$. Higher values of p for a more expensive postprocessing in $\eta _{\textrm{eq}, p}$ do not significantly improve on $\eta _{\textrm{eq}, 1}$.

Fig. 1

Convergence history of the energy error $|\!|\!|e|\!|\!|_{\textrm{pw}}$ (left) and the GUBs $\eta _{\textrm{res}}, \eta _{\textrm{HHO}}, \eta _{\textrm{eq}, 0}$ (right) on the square domain

Fig. 2

History of the overestimation factor $EF(\eta )=\eta /|\!|\!|e|\!|\!|_\textrm{pw}$ for the residual-based error estimators $\eta _{\textrm{res}}, \eta _{\textrm{HHO}}$ (left) and $\eta _{\textrm{eq}, 0}, \eta _{\textrm{eq}, 1}$ (right) on the unit square

5.4 Analytical solution for the slit domain

The source term $f\in L^2(\Omega )$ in the second benchmark on the slit domain $\Omega :=(0,1)^2\!\setminus \!([0,1)\times \{0\})$ matches the singular solution (in polar coordinates)

$$\begin{aligned} u(r, \varphi ) = r^{1/2}(r^2\sin (\varphi )^2 - 1)(r^2\cos (\varphi )^2 - 1)\sin (\varphi /2). \end{aligned}$$

The singularity of u at the origin (0, 0) leads to reduced convergence rates 1/4 under uniform refinement, regardless of the polynomial degree k. Figure 3 shows that the adaptive algorithm recovers optimal rates and verifies the equivalence of the GUBs $\eta _{\textrm{res}}, \eta _{\textrm{eq}, p}$, and $\eta _{\textrm{HHO}}$ to the energy error $|\!|\!|e|\!|\!|_\textrm{pw}$.

The efficiency indices in Fig. 4 show a strong overestimaton by $\eta _{\textrm{res}}$ in the preasymptotic regime (undisplayed) with values $EF(\eta _{\textrm{res}})>60$. However, asymptotically the quotients $EF(\eta )=\eta /|\!|\!|e|\!|\!|_\textrm{pw}$ for the two residual-based GUBs $\eta _{\textrm{res}}, \eta _{\textrm{HHO}}$ differ only by a factor 10, while the equilibrated GUBs $\eta _{\textrm{eq}, p}$ provide the closest values to 1.

Fig. 3

Convergence history of the energy error $|\!|\!|e|\!|\!|_{\textrm{pw}}$ (left) and the GUBs $\eta _{\textrm{res}}, \eta _{\textrm{HHO}}, \eta _{\textrm{eq}, 0}$ (right) on the slit domain

Fig. 4

History of the overestimation factor $EF(\eta )=\eta /|\!|\!|e|\!|\!|_\textrm{pw}$ for the residual-based error estimators $\eta _{\textrm{res}}, \eta _{\textrm{HHO}}$ (left) and $\eta _{\textrm{eq}, 0}, \eta _{\textrm{eq}, 1}$ (right) on the slit domain

5.5 Corner singularity in the L-shaped domain

The third benchmark problem is set in the L-shaped domain $\Omega = (-1, 1)^2\setminus [0, 1)^2$ with constant right-hand side $f\equiv 1$ with an exact solution $u \in H^{1+s}(\Omega )$ for all $0 \le s < 2/3$ [24, Theorem 14.6]. Figure 5 displays the convergence history of the error $e:=u-Ru_h$ and compares the adaptive scheme, Algorithm 1, driven by the refinement indicators $\eta _{\textrm{res}}(T)$ from (68) and

$$\begin{aligned} \eta _{\textrm{HHO}}^2(T)&:=|T|\Vert (I - \Pi _0)(f+{\Delta }_{\textrm{pw}} Ru_h)\Vert _{L^2(T)}^2 + \Vert \nabla (1-{\mathcal {A}}) Ru_h\Vert _{L^2(T)}^2 \\&\quad + |T|^{1/2}\sum _{F\in \mathcal {F}(T)} \Vert R_{T, F}^k u_h\Vert _{L^2(F)}^2,\\ \eta _{\textrm{eq}, 0}^2(T)&:=\textrm{osc}_{k}^2(f, T) + \Vert Q_0^\Delta \Vert _{L^2(T)}^2 + \Vert \nabla (1-{\mathcal {A}}) Ru_h\Vert _{L^2(T)}^2 \end{aligned}$$

for $T\in \mathcal {T}$ that are induced from the GUB $\eta _{\textrm{res}}, \eta _{\textrm{HHO}}$, and $\eta _{\textrm{eq}, 0}$. Here, the norm $|\!|\!|e|\!|\!|_\textrm{pw}$ of the distance e from the discrete solution $Ru_h\in P_{k+1}(\mathcal {T})$ over $\mathcal {T}$ to the unknown solution $u\in H^1(\Omega )$ is approximated by $|\!|\!|\widehat{u}-Ru_h|\!|\!|_\textrm{pw}$, where ${{\widehat{u}}}\in P_{k+1}({\widehat{\mathcal {T}}})$ is the HHO approximation of u on an adaptive refinement ${\widehat{\mathcal {T}}}$ of $\mathcal {T}$ with at least $2|\mathcal {T}|\le |{\widehat{\mathcal {T}}}|$ elements. The three adaptive schemes (Algorithm 1, driven by $\eta _{\textrm{res}}(T), \eta _{\textrm{HHO}}(T)$, or $\eta _{\textrm{eq}, 0}(T)$) recover optimal rates of convergence and lead to similar local refinement of the adaptive mesh sequences as in Fig. 6.

Fig. 5

Convergence history plot of the energy error $|\!|\!|e|\!|\!|_\textrm{pw}$ on the L-shaped domain with uniform and adaptive refinement with AFEM, driven by $\eta _{\textrm{res}}(T)$, (left) and for AFEM, driven by $\eta _{\textrm{res}}(T), \eta _{\textrm{HHO}}(T)$, and $\eta _{\textrm{eq}, 0}(T)$, (right)

Fig. 6

Adaptive triangulations on the L-shaped domain for $k=3$ from AFEM, driven by $\eta _{\textrm{res}}(T)$ (left, $|\mathcal {T}|=882$), driven by $\eta _{\textrm{HHO}}(T)$ (middle, $|\mathcal {T}|=907$), and driven by $\eta _{\textrm{eq}, 0}(T)$ (right, $|\mathcal {T}|=919$)

5.6 Conclusion

The adaptive mesh-refining algorithm recovers optimal convergence rates in all three benchmarks. This holds for AFEM driven by any of the three refinement indicators derived from the GUB $\eta _{\textrm{res}}, \eta _{\textrm{HHO}}$, and $\eta _{\textrm{eq}, p}$. The generated mesh sequences from the adaptive schemes, driven by the different estimators, display a very similar concentration of the local mesh-refinement as in Fig. 6. All three benchmarks verify that the considered error estimators are GUB with reliability constant 1, while the post-processing in the equilibrated GUB $\eta _{\textrm{eq}, p}$ produces minimal overestimation with significant additional computational costs.

Appendix A: Equilibration algorithm for higher order

The post-processed quantity $Q_p\in RT_{k+p}(\mathcal {T})$ from Sect. 4.2 enters the equilibrated error estimator $\eta _{\textrm{eq}, p}(\mathcal {T})$ in Theorem 3 and could be computed by a minimization problem on the vertex patches. The solution property (20) gives rise to the two cases $r = 0$ if $k = 0$ and $r = k+p$ if $k \ge 1$ for the polynomial degree r in the equilibrium $\mathop {\textrm{div}}\limits Q_p + \Pi _{r} f = 0$ in $\Omega$ from (42). This appendix follows [6, 9, 12, 37] to compute the quantity of interest $Q_p - \nabla _\textrm{pw}R u_h$ directly in an efficient two-step procedure in 2D. Throughout this appendix, fix $k,p\in {\mathbb {N}}_0$ and abbreviate ${q}:=k+p$ and $G_h:=\nabla _\textrm{pw}R u_h\in P_{k}(\mathcal {T};\mathbb R^2)$. Let the data $f\in L^2(\Omega )$ be given and assume, for the sake of brevity, that $f \in P_0({\mathcal {T}})$ if $k = 0$.

1.1 A.1 Overview

Recall the definition (41) of the summand $Q_{z,h}$ in $Q_p :=\sum _{z\in {\mathcal {V}}}Q_{z,h}$ from Sect. 4.2 with the piecewise Raviart-Thomas interpolation ${\mathcal {I}}_{\textrm{RT}}:H^1(\mathcal {T};\mathbb R^2)\rightarrow RT_{{q}}^{\textrm{pw}}(\mathcal {T})$ [7, Section III.3.1]. The focus is on one vertex $z\in \mathcal V$ with vertex-patch $\omega (z)$ and its triangulation $\mathcal {T}(z)=\{T\in \mathcal {T}\ |\ z\in T\}$. Consider set of edges ${\mathcal {F}}$ and the facet-spider $\mathcal {F}(z) =\{E\in \mathcal {F}\ |\ z\in E\}$ as in Fig. 7. The nodal basis function $\varphi _z\in S^1(\mathcal {T}(z))$ gives rise to the discrete spaces

$$\begin{aligned}&RT_{{q}}^{\textrm{pw}, 0}({\mathcal {T}}(z)):=\{\sigma _z\in RT_{{q}}^{\text {pw}}(\mathcal {T}(z))\ :\ \sigma _{z}\cdot \nu _E=0 \text { for all }E\in \mathcal {F}\!\setminus \!\mathcal {F}(z)\},\\&{\mathcal {S}}(z):=\left\{ \sigma _z\in RT_{{q}}^{\textrm{pw}, 0}(z)\ :\ \begin{array}{llr} \mathop {\textrm{div}}\limits \sigma _z&{}=-\Pi _{T,{{q}}}(\varphi _z(f+\mathop {\textrm{div}}\limits G_h))&{}\text { for all } T\in \mathcal {T}(z)\\ {[}\sigma _z\cdot \nu _E{]_E} &{}= -\Pi _{E, {{q}}}(\varphi _z[G_h\cdot \nu _E]_E)&{}\text { for all } E\in \mathcal {F}(\omega (z)) \end{array}\right\} . \end{aligned}$$

Proposition 1

(alternative minimization) It holds

$$\begin{aligned} Q_{z,h}^\Delta :=Q_{z,h}- {\mathcal {I}}_{\textrm{RT}}(\varphi _z G_h) = \mathop {\textrm{arg}\,\textrm{min}}\limits _{ \sigma _z\in {\mathcal {S}}(z)} \Vert \sigma _z\Vert _{L^2(\omega (z))}. \end{aligned}$$

(69)

Proof

Recall $f_z = \Pi _{{{\tilde{p}}}}(\varphi _z f - G_h\cdot \nabla \varphi _z)$ from (39). (Notice that this formula coincide with the definition (39) for $k = 0$ because $f \in P_0({\mathcal {T}})$.) Given any $\sigma _z \in {\mathcal {S}}(z)$, the commuting diagram property ${\mathop {\textrm{div}}\limits }_{\textrm{pw}}\circ {\mathcal {I}}_{\textrm{RT}}= \Pi _{{q}}\circ {\mathop {\textrm{div}}\limits }_{\textrm{pw}}$ [7, Proposition 2.5.2] shows

$$\begin{aligned} \textrm{div}_\textrm{pw}(\sigma _z + {\mathcal {I}}_{RT} (\varphi _z G_h)) = \textrm{div}_\textrm{pw}\sigma _z + \Pi _{{{\tilde{p}}}}\,\textrm{div}_\textrm{pw}(\varphi _z G_h) = - f_z. \end{aligned}$$

(70)

By design of the interpolation ${\mathcal {I}}_{RT}$, $({\mathcal {I}}_{RT} (\varphi _z G_h)_{|T}\cdot \nu _E)_{|E} = \Pi _{E,{q}}(\varphi _z G_h)_{|T}\cdot \nu _E$ holds and so $[{\mathcal {I}}_{\textrm{RT}}(\varphi _z G_h)\cdot \nu _E]_E = \Pi _{E, {{q}}}(\varphi _z[G_h\cdot \nu _E]_E)$ follows for any $E \in {\mathcal {F}}(T)$ and $T\in \mathcal {T}$. Therefore, the jump $[\sigma _z + {\mathcal {I}}_{\textrm{RT}}(\varphi _z G_h)]_E\cdot \nu _E \equiv 0$ vanishes on $E \in {\mathcal {F}}(\omega (z))$, whence $\sigma _z + {\mathcal {I}}_{\textrm{RT}}(\varphi _z G_h) \in RT_{{\tilde{p}}}({\mathcal {T}}(z))$. Since $RT_{{\tilde{p}}}^{\textrm{pw},0}({\mathcal {T}}(z)) \cap H(\textrm{div},\omega (z)) = RT_{{\tilde{p}}}^0({\mathcal {T}}(z))$, this and (70) imply $\sigma _z + {\mathcal {I}}_{\textrm{RT}}(\varphi _z G_h) \in {\mathcal {Q}}_h(z)$ for any $\sigma _z \in {\mathcal {S}}(z)$. In particular, ${\mathcal {S}}(z)+ {\mathcal {I}}_{\textrm{RT}}(\varphi _z G_h) \subseteq {\mathcal {Q}}_h(z)$. On the other hand, similar arguments verify the reverse inclusion ${\mathcal {Q}}_h(z) \subseteq {\mathcal {S}}(z)+ {\mathcal {I}}_{\textrm{RT}}(\varphi _z G_h)$. The substitution ${\mathcal {S}}(z)+ {\mathcal {I}}_{\textrm{RT}}(\varphi _z G_h) = {\mathcal {Q}}_h(z)$ in (41) concludes the proof.$\square$

This establishes that the norm $\Vert Q^\Delta _p\Vert$ of $Q^\Delta _p:=\sum _{z\in {\mathcal {V}}}Q_{z,h}^\Delta = Q_p - G_h$ contributes to the equilibrated error estimator and the remaining parts of this appendix compute the minimizer $Q_{z,h}^\Delta$ of (69) in a two-step procedure.

Fig. 7

Triangulation $\mathcal {T}(z)$ and enumeration of the edges $\mathcal {F}(z)$ of the vertex-patch $\omega (z)$ for an internal vertex $z\in \mathcal V(\Omega )$ with $N=|\mathcal {T}(z)|=5$ (left) or boundary vertex $z\in {\mathcal {V}}(\partial \Omega )$ with $N=|\mathcal {T}(z)|=4$ (right)

First, Algorithm 2 generates the coefficients of a particular solution ${\tilde{\sigma }}_z^\Delta \in {\mathcal {S}}(z)$ in terms of the finite element basis ${\mathcal {B}}_{RT}$ of $RT_{{q}}^{\textrm{pw}}({\mathcal {T}}(z))$ from Subsection A.2. The second step computes the correction

$$\begin{aligned} \sigma _z^\Delta :=Q_{z,h}^\Delta - {\tilde{\sigma }}_z^\Delta = \mathop {\textrm{arg}\,\textrm{min}}\limits _{\sigma _z\in V(z)}\Vert \sigma _z+{\tilde{\sigma }}_z^\Delta \Vert _{L^2(\omega (z))} \end{aligned}$$

(71)

in terms of the low-dimensional unconstrained minimization problem over the linear space $V(z):={\mathcal {S}}(z)-{\tilde{\sigma }}_z^\Delta =RT_{{q}}^{\textrm{pw}}({\mathcal {T}}(z))\cap H(\mathop {\textrm{div}}\limits =0, \omega (z))$ characterized in Lemma 12. Because [12, Lemma 3.1] is wrong (take, e.g., $\tau _K=\textrm{curl}\,b_K\ne 0$ for the element bubble function $b_K$ in their notation to see that uniqueness for general polynomial degrees ${q}$ cannot hold) and [12] omits algorithmic details, this appendix focuses on the explicit characterization of the degrees of freedom for the minimization problem (71) over $V(z)$.

1.2 A.2 Degrees of freedom for $RT_{{q}}^{\textrm{pw}}(T)$

This subsection introduces a basis for the Raviart-Thomas finite element on $T\in \mathcal {T}(z)$ and starts with the definition of some linear functionals on $H(\mathop {\textrm{div}}\limits , \mathcal {T})$. For any $\sigma \in H(\mathop {\textrm{div}}\limits , \mathcal {T})$, set

$$\begin{aligned} \lambda _{T, \alpha }^{\ell , m}(\sigma ) :=&\int _T \sigma \cdot e_\alpha \ x_1^\ell x_2^m \textrm{dx}{} & {} 0\le \ell +m\le {q}-1,\alpha =1, 2,\\ \lambda _{T, \mathop {\textrm{div}}\limits }^{\ell , m}(\sigma ) :=&\int _T \mathop {\textrm{div}}\limits \sigma \ x_1^\ell x_2^m \textrm{dx}{} & {} 0\le \ell +m\le {q},\\ \lambda _{T, E}^j(\sigma ) :=&\int _E \sigma _{|T} \cdot \nu _T\ s^j \textrm{ds}{} & {} 0\le j\le {q}, E\in \mathcal {F}(T). \end{aligned}$$

Here and throughout, $e_1=(1,0)$ and $e_2=(0,1)$ denote the canonical unit vectors in $\mathbb R^2$. Note that the (classical) degrees of freedom for the Raviart-Thomas finite element $RT_{q}(T)$ of degree ${q}\in {\mathbb {N}}_0$ from [11] read

$$\begin{aligned} \widetilde{\Lambda _T}:=\{\lambda _{T, E}^j, \lambda _{T, \alpha }^{\ell , m}\,\ \text {for }0\le j\le {q}, 0\le \ell +m\le {q}-1, \alpha =1,2, E\in \mathcal {F}(T)\}. \end{aligned}$$

This appendix requires, for the construction of ${\tilde{\sigma }}_z^\Delta \in {\mathcal {S}}(z)$, a different set of (unisolvent) degrees of freedom $\Lambda _T$ for $RT_{q}(T)$ that includes the edge and divergence moments

$$\begin{aligned} \Lambda _T^0:=\{\lambda _{T, E}^j, \lambda _{T, \mathop {\textrm{div}}\limits }^{\ell , m}\ :\ \text {for }0\le j\le {q}, 1\le \ell +m\le {q}, E\in \mathcal {F}(T)\}\subseteq \Lambda _T. \end{aligned}$$

(72)

(The set $\Lambda _T^0$ itself is linear independent [37, Lemma 3.1].) Given any $\Lambda _T$ with (72), denote the remaining $N_{q}={q}({q}-1)/2$ degrees of freedom $\Lambda _T\setminus \Lambda _T^0$ by $\lambda _T^{r}$ for $r=1,\dots ,N_{q}$. Let ${\mathcal {B}}_{RT,T}=\{\varphi _{T, E}^j, \varphi _{T, \mathop {\textrm{div}}\limits }^{\ell , m}, \varphi _{T}^r\}$ be the unique basis of $RT_{q}(T)$ dual to $\Lambda _T$ with inferred indices from $\Lambda _T$. Then, the collection ${\mathcal {B}}_{RT}:=\bigcup _{T\in \mathcal {T}(z)} {\mathcal {B}}_{RT,T}$ is a basis of $RT_{{q}}^{\textrm{pw}}({\mathcal {T}}(z))$ and any function $\sigma _z\in RT_{{q}}^{\textrm{pw}}({\mathcal {T}}(z))$ has the representation

$$\begin{aligned} \sigma _z:=\sum _{T\in \mathcal {T}(z)}\left( \sum _{E\in \mathcal {F}(T)}\sum _{j=0}^{q}c_{T, E}^j\varphi _{T, E}^j+ \sum _{1\le \ell +m\le {q}}c_{T, \mathop {\textrm{div}}\limits }^{\ell , m}\varphi _{T, \mathop {\textrm{div}}\limits }^{\ell , m}+ \sum _{r=1}^{N_{q}}c_{T}^r\varphi _{T}^r\right) \end{aligned}$$

(73)

with coefficients $c_{T, E}^j= \lambda _{T, E}^j(\sigma _z), c_{T, \mathop {\textrm{div}}\limits }^{\ell , m}= \lambda _{T, \mathop {\textrm{div}}\limits }^{\ell , m}(\sigma _z)$, and $c_{T}^r= \lambda _T^{r}(\sigma _z)$ for all $T\in \mathcal {T}(z), E\in \mathcal {F}(T)$, and $0\le j\le {q}, 1\le \ell +m\le {q}, 1\le r\le N_{q}$. By duality, the coefficients $c_{T, E}^j$ with $0\le j\le {q}$ uniquely determine the normal trace $(\sigma _{z|T})_{|E}\cdot \nu _E\in P_{q}(E)$ on the edge $E\in \mathcal {F}(T)$ of $T\in \mathcal {T}(z)$. Any set of degrees of freedom $\Lambda _T$ with (72) works with the equilibration algorithm in A.3.

Example 3

(Construction of $\Lambda _T$) This example presents a generic procedure to obtain such a set from $\widetilde{\Lambda _T}$. The integration by parts formula shows that the lowest-order divergence moment $\lambda _{T, {\text { div }}}^{0,0} = \sum _{E\in \mathcal {F}(T)} \lambda _{T, E}^0$ depends linearly on the lowest-order edge moments and, similarly, the sums

$$\begin{aligned} \lambda _{T, \mathop {\textrm{div}}\limits }^{\ell , m}+ \ell \lambda _{T, 1}^{\ell -1, m} + m\lambda _{T, 2}^{\ell , m-1}\in (P_{\ell + m}(\mathcal {F}(T)))^* \end{aligned}$$

(74)

are functionals on $P_{\ell + m}(\mathcal {F}(T))$ (summands with negative indices are understood as zero). This relation allows for the substitution of volume moments in $\widetilde{\Lambda _T}$ for divergence moments $\lambda _{T, \mathop {\textrm{div}}\limits }^{\ell , m}$, $1\le \ell +m\le {q}$, and leads to $\Lambda _T$ with (72). The remaining degrees of freedom $\Lambda _T\setminus \Lambda _T^0$ are volume moments of the form $\lambda _T^{r}=\lambda _{T, \alpha }^{\ell , m}$ for a fixed $\alpha \in \{1,2\}$, e.g.,

$$\begin{aligned} \Lambda _T\setminus \Lambda _T^0=\{\lambda _{T, 2}^{\ell , m}\,\ 1\le \ell \le {q}-1, 0\le m\le {q}-1-\ell \}. \end{aligned}$$

1.3 A.3 Equilibration algorithm

Algorithm 2

Particular solution in ${\mathcal {S}}(z)$

This subsection presents the equilibration procedure, starting with Algorithm 2, that computes an admissible function ${\tilde{\sigma }}_z^\Delta \in {\mathcal {S}}(z)$ in terms of the representation (73). Enumerate the $N:=|\mathcal {T}(z)|$ triangles $T\in \mathcal {T}(z)$ from 1 to N as in Fig. 7. Any two neighbouring triangles $T_a, T_{a+1}$ share an edge $E_a:=T_a\cap T_{a+1}$ for $a=1,...,N-1$. If $z\in {\mathcal {V}}(\Omega )$ is an interior vertex, $T_1$ and $T_{N}$ share an additional edge $E_0:=E_{N}:=T_1\cap T_{N}$. For a boundary vertex $z\in {\mathcal {V}}(\partial \Omega )$, $T_1, T_{N}$ have the distinct boundary edges $E_0, E_{N}\in \mathcal {F}(z)\cap \mathcal {F}(\partial \Omega )$. The following lemma shows correctness of Algorithm 2 under the compatibility condition (5) and represents step one of the equilibration algorithm. The final step is the local minimization problem in Lemma 12 that provides $Q_{z,h}^\Delta$ from (69). Both proofs are provided in A.4.

Lemma 11

Given $z\in {\mathcal {V}}$, let $\{\varphi _{T, E}^j, \varphi _{T, \mathop {\textrm{div}}\limits }^{\ell , m}, \varphi _{T}^r\}$ be the basis of $RT_{q}(T)$ dual to $\Lambda _T$ with (72) for all $T\in \mathcal {T}(z)$. Suppose $f \in L^2(\omega (z))$ and $G_h\in P_{q}(\mathcal {T}(z);{\mathbb {R}}^2)$ satisfy

$$\begin{aligned} (G_h, \nabla \varphi _z)_{L^2(\omega (z))}=(f,\varphi _z)_{L^2(\omega (z))}\quad \text {if }z\in {\mathcal {V}}(\Omega ). \end{aligned}$$

(75)

Then the output of Algorithm 2 with input f and $G_h$ defines a function ${\tilde{\sigma }}_z^\Delta \in {\mathcal {S}}(z)$.

Note that (75) is a local version of (5) and therefore holds for the HHO method with the choice $G_h:=\nabla _\textrm{pw}R u_h$ as proven in (20). This allows for the computation of $Q^\Delta _p :=\sum _{z\in {\mathcal {V}}}Q_{z,h}^\Delta = Q_p - G_h$ in terms of local and unconstrained minimization problems on the vertex-patches $\omega (z)$.

Lemma 12

Given $z\in {\mathcal {V}}$, let $\{\varphi _{T, E}^j, \varphi _{T, \mathop {\textrm{div}}\limits }^{\ell , m}, \varphi _{T}^r\}$ be as in Lemma 11 for all $T\in \mathcal {T}(z)$ and let ${\tilde{\sigma }}_z^\Delta \in {\mathcal {S}}(z)$ be arbitrary. Then $V(z):={\mathcal {S}}(z)- {\tilde{\sigma }}_z^\Delta$ is a linear vector space and consists of all functions of the form

$$\begin{aligned} \sum _{a=1}^{N}\left( d_0 (\varphi _{T_a, E_{a-1}}^0-\varphi _{T_{a}, E_{a}}^0) + \sum _{\ell =1}^{q}(d_{E_{a-1}}^\ell \varphi _{T_a, E_{a-1}}^\ell -d_{E_{a}}^\ell \varphi _{T_a, F_a}^\ell ) +\sum _{r=1}^{N_{q}}d_{T_a}^r\varphi _{T_a}^r\right) \end{aligned}$$

(76)

for arbitrary $d_0, d_{E_a}^\ell , d_{T_a}^r\in \mathbb R$ with $\ell =1,...,{q}, r=1,...,N_{q}, a=1,...,N$ (and $d_{E_0}^\ell =d_{E_{N}}^\ell$ for $z\in {\mathcal {V}}(\Omega )$) and the enumeration of $\mathcal {T}(z)$ as in Fig. 7. Furthermore, $Q_{z,h}^\Delta ={\tilde{\sigma }}_z^\Delta + \sigma _z^\Delta$ holds for the solution $\sigma _z^\Delta \in V(z)$ to the $1+{q}|\mathcal {F}(z)|+{q}({q}-1)/2N$-dimensional minimization problem (71).

1.4 A. 4 Proofs

The remaining parts of this appendix are devoted to the verification of Lemmas 11–12.

Proof

(of Lemma 11) Enumerate $\mathcal {T}(z)$ as in A.3 and recall the definition of the jump $[G_h]_E = G_h|_{T_+} - G_h|_{T_-}$ on the interior edge $E=T_+\cap T_-$ shared by $T_+, T_-\in \mathcal {T}$, and $[G_h]_E = G_h|_{T_+}$ for the unique triangle $T_+\in \mathcal {T}$ with $E\subset T_+$ for the boundary edge $E\in \mathcal {F}(\partial \Omega )$. First, observe that $\sigma _z\in RT_{{q}}^{\textrm{pw}}({\mathcal {T}}(z))$ lies in $RT_{{q}}^{\textrm{pw}, 0}(\mathcal {T}(z))$ if and only if the coefficients $c_{T, E}^j=0$ in the representation (73) are zero for $0\le j\le {q}$ at the edge $E\in \mathcal {F}(T){\setminus }\mathcal {F}(z)$ in $T\in \mathcal {T}$ opposing z. By definition, $\sigma _z\in RT_{{q}}^{\textrm{pw}, 0}(\mathcal {T}(z))$ belongs to ${\mathcal {S}}(z)$ if and only if

$$\begin{aligned} \lambda _{T, \mathop {\textrm{div}}\limits }^{\ell , m}(\sigma _z)&= (f+\mathop {\textrm{div}}\limits G_h, \varphi _zx_1^\ell x_2^m)_{L^2(T)}{} & {} \text {for all }0\le \ell +m\le {q},T\in \mathcal {T}(z), \end{aligned}$$

(77)

$$\begin{aligned} (\lambda _{T_+, E}^{j} + \lambda _{T_-, E}^{j})(\sigma _z)&= ([G_h]_E\cdot \nu _E, \varphi _zs^j)_{L^2(E)}{} & {} \text {for all }0\le j\le {q}, E\in \mathcal {F}(\omega (z)). \end{aligned}$$

(78)

This translates into equivalent conditions on the coefficients of $\sigma _z$ in the representation (73), namely, for all $a=1,..., N$,

$$\begin{aligned} c_{T_a, E_{a}}^{0}&= (f+\mathop {\textrm{div}}\limits G_h, \varphi _z)_{L^2(T_a)} - c_{T_a, E_{a-1}}^{0},{} & {} \end{aligned}$$

(79)

$$\begin{aligned} c_{T_a, \mathop {\textrm{div}}\limits }^{\ell , m}&= (f+\mathop {\textrm{div}}\limits G_h, \varphi _zx_1^\ell x_2^m)_{L^2(T_a)}{} & {} \text {for all }1\le \ell +m\le {q}, \end{aligned}$$

(80)

$$\begin{aligned} c_{T_{a}, E_{a-1}}^j&= d_{E_{a-1}}^j{} & {} \text {for all }0\le j\le {q}, \end{aligned}$$

(81)

$$\begin{aligned} c_{T_a, E_a}^j&= ([G_h]_E\cdot \nu _E, \varphi _zs^j)_{L^2(E_a)} - d_{E_a}^j{} & {} \text {for all }0\le j\le {q}, \end{aligned}$$

(82)

where $d_{E_a}^\ell \in \mathbb R$. Since $\sigma _z\in RT_{{q}}^{\textrm{pw}, 0}(\mathcal {T}(z))$ vanishes at the other edges $E\in \mathcal {F}\setminus \mathcal {F}(z)$, $\lambda _{T_a, \mathop {\textrm{div}}\limits }^{0,0}(\sigma _z) = c_{T_a, E_{a-1}}^{0} + c_{T_a, E_{a}}^{0}$ and (79)–(80) are equivalent to (77). The identification $d_{E_0}^\ell =d_{E_{N}}^\ell$ for an interior vertex $z\in \mathcal V(\Omega )$ with $E_0=E_N\in \mathcal {F}(z)$ shows that (81)–(82) are equivalent to (78). This identification is well defined. Note that, whereas there is no condition on $d_{E_a}^\ell$ for $1\le \ell \le {q}$, the combination of (79) and (81) with (82) shows the implicit extra condition

$$\begin{aligned} d_{E_a}^0 = d_{E_0}^0 + \sum _{\alpha =1}^a\left( ([G_h]_E\cdot \nu _E, \varphi _z)_{L^2(E_\alpha )} - (f+\mathop {\textrm{div}}\limits G_h, \varphi _z)_{L^2(T_\alpha )}\right) \quad \text {for }a=0,...,N. \end{aligned}$$

For an interior vertex $z\in {\mathcal {V}}(\Omega )$, an integration by parts and (75) show that the sum on the right-hand side above vanishes for $a=N$, whence $d_{E_N}^0 = d_{E_0}^0$ is indeed well defined. Furthermore, there is no condition on the coefficients $c_{T_a}^{r}$ for all $T_a\in \mathcal {T}(z)$ and $r=1,...,N_{q}$ and $c_{T_a}^{r}=d_{T_a}^r$ is a further degree of freedom.

Algorithm 2 finds coefficients that satisfy (79)–(82) in a loop over $a=1,...,N$ and therefore defines ${\tilde{\sigma }}_z^\Delta \in {\mathcal {S}}(z)$ by (73). $\square$

Proof

(of Lemma 12) This follows immediately after revisiting the proof of Lemma 11 for an arbitrary function $\sigma _z\in RT_{{q}}^{\textrm{pw}, 0}(\mathcal {T}(z))$. Since $\sigma _z\in {\mathcal {S}}(z)$ is equivalent to (79)–(82) for the representation (73) of $\sigma _z$ in the given basis, all functions $\sigma _z\in {\mathcal {S}}(z)-{\tilde{\sigma }}_z^\Delta$ are of the form (76) for arbitrary $d_0, d_{E_a}^\ell , d_{T_a}^r\in \mathbb R$ with $\ell =1,...,{q}$, $r=1,...,N_{q}$, and $a=1,...,N$ (and $d_{E_0}^\ell =d_{E_{N}}^\ell$ for $z\in {\mathcal {V}}(\Omega )$). Hence, the dimension of the linear space $V(z)={\mathcal {S}}(z)-{\tilde{\sigma }}_z^\Delta$ is $1+{q}|\mathcal {F}(z)|+{q}({q}-1)/2N$. The claim follows from Proposition 1 by observing

$$\begin{aligned} Q_{z,h}^\Delta :=&\mathop {\textrm{arg}\,\textrm{min}}\limits _{ \sigma _z\in {\mathcal {S}}(z)} \Vert \sigma _z\Vert _{L^2(\omega (z))}= {\tilde{\sigma }}_z^\Delta + \mathop {\textrm{arg}\,\textrm{min}}\limits _{ \sigma _z\in {\mathcal {S}}(z)- {\tilde{\sigma }}_z^\Delta } \Vert {\tilde{\sigma }}_z^\Delta + \sigma _z\Vert _{L^2(\omega (z))}.\quad \square \end{aligned}$$

$\square$