Stabilization-free HHO a posteriori error control

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.


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 parameterfree stabilization term for the link between the volume and skeletal variables.

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 ∈ 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 ∈ H 1 (Ω), cf.(2) and Section 3 below for further details.Then a possible error term reads total error 2 := ∇ pw (u It is disputable if s h (u h , u h ) ≥ 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 with the abbreviation terms of the L 2 projections Π K,k onto polynomials of degree at most k on a facet or simplex K ∈ F ∪ T of diameter h K = diam(K); cf.Subsection 1.4 for further details.The original residual-based estimator η HHO from the textbook [25] for the Poisson model problem −∆u = f includes (2) and an interpolation ARu h ∈ V of Ru h by nodal averaging in + s h (u h , u h ).
(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 ∈ P k+1 (T ), total error 2 ≈ η 2 HHO .

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 Second, since the stabilization (2) incorporates a negative power of the meshsize, 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) for some parameter p ∈ N 0 .The explicit residual-based a posteriori error estimator η res follows from the a posteriori methodology in the spirit of [14,20,18,21] with a piecewise volume residual f +∆ pw Ru h and the jumps [∇ pw Ru h ] F across a facet F (on the boundary this is only the tangential component of ∇ pw Ru h ).The equilibrated error estimator η eq,p includes the post-processed quantity Q p ∈ RT k+p (T ) in the space RT k+p (T ) of Raviart-Thomas functions of degree k + p for p ∈ N 0 and the nodal average ARu h ∈ S k+1 0 (T ) of Ru h .The main results establish reliability and efficiency error 2 + osc 2  k−1 (f, T ) η 2 res ≈ η 2 eq,p error 2 + osc 2 q (f, T ) for any p, q ∈ N 0 up to data-oscillations osc 2 q (f, T ) := h T (1 of the source term f ∈ L 2 (Ω) and without any stabilization terms.Computational benchmarks with adaptive mesh-refinement driven by any of these estimators provide numerical evidence for optimal convergence rates.

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 1conforming 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 ∇ pw Ru h .In fact, those two fairly general properties lead in Section 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 error ≤ η res and error ≤ η eq,p ; cf. Table 1 for explicit quantities and Theorem 1 and Theorem 3 for further details.Numerical comparisons of η HHO with η res and η 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 η res in Section 3. The alternative stabilization-free error estimator η eq,p follows from an equilibration strategy plus post-processing in Section 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 Section 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 Q p − ∇ pw Ru h L 2 (Ω) in η eq,p .

Overall notation
Standard notation for Sobolev and Lebesgue spaces and norms apply with In particular, H(div, Ω) is the space of Sobolev functions with weak divergence in L 2 (Ω) and H(div = 0, Ω) contains only divergence-free functions in H(div, Ω).Throughout this paper, T denotes a shape-regular triangulation of the polyhedral bounded Lipschitz domain Ω ⊂ R n into n-simplices with facets F (edges for n = 2 and faces for n = 3) and vertices V. Let F(Ω) (resp.V(Ω)) denote the set of interior facets (resp.vertices) and F(∂Ω t and Curl w := curl w if n = 3.For s ∈ R, let H s (T ), H(div, T ), and H(curl, T ) denote the space of piecewise Sobolev functions with restriction to T ∈ T in H s (T ), H(div, T ), and H(curl, T ).To simplify notation, H s (K) abbreviates H s (int(K)) for the open interior int(K) of a compact set K. The L 2 -scalar product reads (•, •) L 2 (ω) for volumes ω ⊆ Ω and •, • L 2 (γ) for surfaces γ ⊂ Ω of co-dimension one; the same symbol applies to scalars and to vectors.For Define the energy scalar product a(v, w) := (∇v, ∇w) L 2 (Ω) for v, w ∈ H 1 (Ω) and its piecewise version Here and throughout the paper, ∇ pw , div pw , curl pw , ∆ pw , denote the piecewise evaluation of the differential operators ∇, div, curl pw , ∆ without explicit reference to the underlying shape-regular triangulation T .
The vector space P k (K) of polynomials of degrees at most k ∈ N 0 over a facet or simplex K ∈ F ∪ T defines the piecewise polynomial spaces and the space of piecewise Raviart-Thomas functions The associated If not explicitly stated otherwise, constants are independent of the meshsize in the triangulation but may depend on the shape-regularity and on the polynomial degree k.The abbreviation A B hides a generic constant C (independent of the mesh-size) in A ≤ C B; A ≈ B abbreviates A B A.

Foundations of the a posteriori error analysis
This section investigates general building blocks of the a posteriori error analysis and revisits arguments from [14,20,18,21] with emphasis on multiply connected domains Ω ⊂ 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 ∈ V = H 1 0 (Ω) to the Poisson model problem −∆u = f a.e. in Ω and u = 0 on ∂Ω for a given source f ∈ L 2 (Ω); i.e., u ∈ V satisfies An approximation G ∈ L 2 (Ω; R n ) of the gradient ∇u ∈ H(div, Ω) gives rise to the residual f + div G ∈ V * = H −1 (Ω) seen as a linear functional on V , i.e., Let ν T denote the unit outer normal along the boundary ∂T of each simplex T ∈ T and fix the orientation of a unit normal ν F = ±ν T for each facet F ∈ F(T ) of T such that it matches the outer unit normal ν of ∂Ω at the boundary.
The jump The main result of this section establishes the residual-based error estimator The second assumption is the orthogonality to the lowest-order divergence-free Raviart-Thomas functions (G, r) L 2 (Ω) = 0 for all r ∈ RT 0 (T ) ∩ H(div = 0, Ω).
Theorem 1 (residual-based GUB) Suppose that G ∈ H 1 (T ; R n ) and f ∈ L 2 (Ω) satisfy (5)- (6).Then the error estimator η(T , G) from (4) is a GUB of the error ∇u − G for the solution u ∈ V to (3).The constants C 1 , C 2 , C H exclusively depend on Ω and the shape-regularity of 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 1 , C 2 , and C 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 ∇u − G into some divergence part and some consistency part.
The split (7) of the error ∇u−G 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 [13,4] and introduces the constant C 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 ∈ H(div = 0, Ω).
Lemma 2 (Helmholtz-decomposition) Suppose the divergence-free function The constant C H > 0 exclusively depends on Ω.
Proof The compact polyhedral boundary ∂Ω of the bounded Lipschitz domain Ω has J + 1 connectivity components Γ 0 , . . ., Γ J for some finite J ∈ N 0 .Those connectivity components have a positive surface measure |Γ j | and a positive distance of each other.So the integral mean is well defined and depends continuously on ∈ H(div = 0, Ω) in the sense that |γ j | ≤ c 1 (recall div = 0) for each j = 0, . . ., J and c 1 > 0. This constant c 1 and the constants c 2 , c 3 , c 4 below exclusively depend on the domain Ω.The finite real numbers γ 0 , . . ., γ J define the Neumann data for the harmonic function z ∈ H 1 (Ω)/R with ∆z = 0 in Ω and ∂z/∂ν = γ j on Γ j for all j = 0, . . ., J.
The elliptic regularity theory for polyhedral domains lead to z ∈ H 1+α (Ω) for some α > 1/2 and c 2 > 0 with z The Raviart-Thomas interpolation operator defines a bounded linear operator on H(div, Ω) ∩ L p (Ω; R n ) for p > 2. It is generally accepted that, for α > 0 and ∇z ∈ H(div = 0, Ω) ∩ H α (Ω; R n ), the Fortin interpolation I F ∇z ∈ RT 0 (T ) ∩ H(div = 0, Ω) is well defined and I F ∇z ≤ c 3 ∇z H α (Ω) follows for some c 3 > 0. The additional property ∇z ∈ L p (Ω; R n ) for some p > 2 allows the definition of F ∇z • ν F ds as a Lebesgue integral over a facet F ∈ F. One consequence for the boundary facets is the vanishing integral Γj ( − I F ∇z) • ν ds = 0 for all j = 0, . . ., J.
Standard arguments like the trace identity on T j ⊂ T [17, Lemma 2.1] for |f | 2 ∈ W 1,1 (T ) and a Cauchy inequality show The distance max Since the centroid 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 tr < n/(n + 1) follows and cannot be improved in the absence of further assumptions on the shape of the simplex T .Since , the previously displayed estimate leads to Let T be the refinement of T , obtained by replacing T ∈ T with T 0 , ..., T d from above.The triangulation T allows for the facet based decomposition Since the family {ω (F ) : F ∈ F} has no overlap, the sum of the last displayed inequality over all F ∈ F and a Cauchy inequality conclude the proof of Lemma 3.
The next lemma utilizes a quasi-interpolation operator J : H 1 (Ω) → S 1 (T ) with the restriction J(V ) ⊂ S 1 0 (T ), e.g., J = J 1 •I N C from [19, Section 5] with explicit constants for n = 2, and the approximation and stability properties for constants C 1 and C st exclusively depending on the shape-regularity of T .
For the precise definition of J 1 and I N C , we refer to [19, eq. (47) and Section 5].
The final ingredient for the proof of Theorem 1 controls the second term δ in the decomposition of Lemma 1 for := G − ∇w.Recall C H from Lemma 2 and C 1 from ( 11), and C 2 from page 9.
Lemma 5 (conformity) Suppose the divergence-free function Proof Lemma 2 provides β ∈ H 1 (Ω; R N ) with (10) for a (component-wise) quasi-interpolation β C ∈ S 1 (T ) N with (11) as in the proof of Lemma 4; set ψ := β − β C .A piecewise integration by parts and the collection of jump contributions shows Stability and approximation properties of the quasi-interpolation (11) and the trace inequality of Lemma 3 eventually lead to 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 |||β||| ≤ C H conclude the proof.
Proof (Theorem 1) The trace of The assumption ( 6) on G and an integration by parts prove that Lemma 5 is applicable to := G−∇w ∈ H(div = 0, Ω)∩H(curl, T ; R n ).Since curl ∇w = 0 and ∇w × ν = 0, this reveals The above estimates together with the decomposition of Lemma 1 establish η(T , G) as a GUB for the error ∇u − G .
Example 1 (constants for right-isosceles triangles) Section 5] of the quasi-interpolation operator J in the proof of Lemma 4 allows for the explicit estimates where j 1,1 = 3.8317 denotes the first positive root of the first Bessel function.
For triangulations into right-isosceles triangles, the constant Lemma 4.8] depends on the domain by the maximal number M bd ≤ 4 max{π, ω max }/π ≤ 8 of triangles sharing a boundary vertex.Given the maximal interior angle ω max of Ω, Table 1 displays those constants for the maximal possible value M bd = 4 max{π, ω max }/π.The geometric quantity max x∈T |x − mid(T )| equals two-thirds of the maximum median of T .Thus, 3 Explicit residual-based a posteriori HHO error estimator The arguments from Section 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.

Hybrid high-order methodology
The HHO ansatz space reads Let u h ∈ V h solve the HHO discrete formulation of (3) with for the HHO bilinear form and the stabilization term s h (u h , v h ) from (2).Given any w C ∈ S 1 0 (T ) = P 1 (T ) ∩ H 1 0 (Ω), the definition of the reconstruction operator R in (17) verifies RIw C = w C with the interpolation I onto V h .Hence, S T F Iw C = 0 vanishes for all F ∈ F(T ) and T ∈ T .This and (18) show, for all w C ∈ S 1 0 (T ), that

Explicit a posteriori error estimator
As a result of (20), ∇ pw Ru h satisfies the solution property (5) if k ≥ 1 and, in the lowest order case k = 0, (5) holds with f replaced by Π 0 f .This allows the application of the theory from Section 2 to the HHO method with minor modifications for the case k = 0. Define the error estimator contributions Since ∇ pw Ru h is a piecewise gradient, its piecewise curl vanishes.This leads to the explicit residual-based a posteriori error estimator (Recall C 1 , C 2 from Lemma 4 and C H from Lemma 2 as well as the Poincaré constant C P ≤ π −1 .)The main result of this section verifies the assumptions in Theorem 1 and proves reliability and efficiency of η res (T ).
Theorem 2 (residual-based GUB for HHO) Let u ∈ V solve the Poisson equation (3) and let u h ∈ V h solve the discrete formulation (18).Then and osc k−1 (f, T ) ≤ C 4 η res (T ) hold for any q ∈ N 0 .The constants C 3 and C 4 exclusively depend on k, q and on the shape-regularity of the triangulation T .

Proof of Theorem 2
The orthogonality of ∇ pw Ru h to the divergence-free Raviart-Thomas space of lowest degree is an assumption in Theorem 1 and verified below.
The following lemma concerns the efficiency of the jump contributions.Each facet F ∈ F has at most two adjacent simplices that define a triangulation Lemma 7 (efficiency of jumps) The solution u ∈ V to (3) and the discrete solution u h ∈ V h to (18) satisfy (a) for all F ∈ F and (b) for all F ∈ F(Ω).
(a) h Proof The proof is based on the following extension argument.Given a polynomial p ∈ P k (F ) of degree at most k along the side F ∈ F, the coefficients determine a polynomial (also denoted by p) along the hyperplane H that enlarges F .The intersection F := H ∩ conv(ω(F )) of the hyperplane H with the convex hull of the facet-patch ω(F ) may be strictly larger than F .The shaperegularity of T bounds the size of F in terms of F and an inverse estimate leads to a bound p L ∞ ( F ) ≤ C(k) p L ∞ (F ) with a constant C(k) that depends on the shape-regularity of T and on k.The extension of p from H to R n by constant values along the side normal ν F leads to a polynomial p ∈ P k (R n ) with Proof of (a).The tangential jump be one of the components of F ∈ P k (F ) N , for j = 1, ..., N , and extend it as explained above to p ∈ P k (R n ) and call this (j) in the vector F ∈ P k (R n ; R N ).The proof involves the piecewise polynomial facet-bubble function b F := n n Π n j=1 ϕ j for the n nodal basis function ϕ 1 , . . ., ϕ n ∈ S 1 (T ) associated with the vertices of F .An inverse estimate [38,Proposition 3.37] shows 26) and a piecewise integration by parts show This and (curl , ∇v) L 2 (Ω) = 0 for any v ∈ V imply An inverse estimate, b F L ∞ (ω(F )) = 1, and (25) imply 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.
The following lemma reveals that the order k ≥ N 0 of the oscillations 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 Π −1 := 0.
Lemma 8 (efficiency of lower-order oscillations) Given any simplex T ∈ T and parameters k, q ∈ N 0 , the solution u ∈ V to (3) satisfies The constant C 5 exclusively depends on q and the shape of T .
Proof The assertion ( 30) is trivial for q ≤ k − 1, so suppose k ≤ q.Any v k+1 ∈ P k+1 (T ) and T := Π q f + ∆v k+1 ∈ P q (T ) satisfy 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 T ∈ S q+n+1 0 (T ) ⊂ V and the weak formulation (3) result in A Cauchy inequality, the inverse estimate with a constant C 7 that exclusively depends on q + n + 1 and the shape of T , and (32) lead to The combination of ( 31) with ( 33) and a Cauchy inequality conclude the proof of ( 30), e.g., with Proof (of Theorem 2) Recall the definition of η res (T ) for k ≥ 1 and k = 0 in (21).Since osc 2 k−1 (f, T ) ≤ η 2 res,1 (T ) + η 2 res,2 (T ) η 2 res (T ), the remaining parts of this proof discuss the reliability and efficiency of η res (T ).
Lemma 6 provides the orthogonality of ∇ pw Ru h ∈ H 1 (T ; R n ) to the divergence-free Raviart-Thomas function RT 0 (T ) ∩ H(div = 0, Ω).This and (20) show that the assumptions in Theorem 1 hold for G := ∇ pw Ru h and k ≥ 1, whence the reliability of η res (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 This provides the reliability and it remains to verify the efficiency η res (T ) |||u − Ru h ||| pw + osc q (f, T ) for any q ∈ N 0 .The Pythagoras theorem and (33) with T := Π k f + ∆Ru h ∈ P k (T ) and v k+1 := Ru h ∈ P k+1 (T ) lead to the local efficiency of the volume contributions Lemma 7 considers the remaining terms in the error estimator and establishes their efficiency namely, with the modified jump [∇ pw Ru h ] F = ∇ pw Ru h × ν F on boundary facets F ∈ F(∂Ω).This and Lemma 8 establish the existence of some mesh-independent constant C 3 > 0 with C −1 3 η res (T ) ≤ |||u − Ru h ||| pw + osc q (f, T ) for arbitrary q ∈ N 0 .This concludes the proof.
While the focus of this paper is on the HHO methodology, the framework of Section 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 (T ) × P k (F(Ω)) for k ∈ N 0 .An equivalent formulation to the HDG FEM from [35] Here, R : V h → P k+1 (T ) is defined as in (17) and for any v h = (v T , v F ), w h = (w T , w F ) ∈ V h .This method is also known under the label of weak Galerkin FEM [39].It is straightforward to verify that ∇ pw Ru 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 The efficiency η res (T ) |||u−Ru h ||| pw follows from the arguments in the proofs of Lemma 7-8.

Equilibrium-based a posteriori HHO error analysis
The residual-based guaranteed upper bound (GUB) of the error |||u − Ru h ||| pw from Subsection 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,30,6,5] to establish, up to the well-known Poincaré constant C P ≤ 1/π, a constant-free guaranteed upper bounds for a tight error control.

Guaranteed error control
The guaranteed upper bounds of this section involves two post-processings of the potential reconstruction Ru h ∈ P k+1 (T ) of the discrete solution u h to (18).First, the patch-wise design of a flux reconstruction Q p ∈ RT k+p (T ) with p ∈ N 0 from [3,10,31] provides an H(div, Ω)-conforming approximation to ∇ pw Ru h with the equilibrium Π r f + div Q p = 0 in Ω and r from (35) below.Second, the nodal average ARu h ∈ S k+1 0 (T ) ⊂ V results in an Vconforming approximation of Ru h by averaging all values of the discontinuous function Ru h at each Lagrange point of S k+1 0 (T ).This, the split (7), and the solution property (20) give rise to the guaranteed upper bound (GUB) with r ∈ N 0 defined by The main result of this section states the reliability and efficiency (up to data oscillations) of η eq,p for all parameters p ∈ N 0 .
Theorem 3 (equilibrium-based GUB for HHO) Let u ∈ V resp.u h ∈ V h solve (3) resp.(18).Given a parameter p ∈ N 0 , there exists Q p ∈ RT k+p (T ) such that the error estimator η eq,p (T ) from ( 34) is an efficient GUB for any q ∈ N 0 and osc k−1 (f, T ) ≤ C 9 η eq,p (T ).The constants C 8 and C 9 exclusively depend on the polynomial degree k ∈ N 0 , the parameter q ∈ N 0 , and the shape-regularity of 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 sh (v h , v h ) := T ∈T sT (v h , v h ) from [25] defined, for 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, The subsequent subsection continues with some explanations on the flux reconstruction Q p ∈ RT k+p (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 − ∇ pw Ru h in 2D.The efficiency of the averaging |||(1 − A)Ru h ||| pw up to data oscillations follows in Subsections 4.3-4.4.Subsection 4.6 concludes with the proof of Theorem 3.

Construction of equilibrated flux
This subsection defines the post-processed H(div, Ω)-conforming equilibrated flux Q p ∈ RT k+p (T ) that enters the GUB η eq,p from (34) based on local patch-wise minimization problems in the spirit of [10,30,31].
Consider the shape-regular vertex-patch ω(z) := int( {T ∈ T (z)}) covered by the neighbouring simplices T (z) := {T ∈ T : z ∈ T } sharing a given vertex z ∈ V with the facet spider F(z) := {F ∈ F : z ∈ F }. Recall the space of piecewise Raviart-Thomas functions RT pw k (T ) from Subsection 1.4 and define Throughout the remaining parts of this section, abbreviate G h := ∇ pw Ru h ∈ P k (T ; R n ).Given a vertex z ∈ V with the P 1 -conforming nodal basis function ϕ z ∈ S 1 (T ), the property (20) provides compatible data such that the discrete affine space is not empty.Consequently, is well defined as the with the piecewise Raviart-Thomas interpolation ) is a piecewise Raviart-Thomas function of degree k + p. Hence I RT (ϕ z G h ) = ϕ z G h and I RT could be omitted in the formula (41).The partition of unity z∈V ϕ z ≡ 1 and 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(div) minimization by residual [10,31]) Given any vertex z ∈ V, a piecewise Raviart-Thomas function σ z ∈ RT pw q (T (z)) and a piecewise polynomial r z ∈ P q (T (z)) of degree q ∈ N 0 , define the residual for all v ∈ H 1 (Ω).If z ∈ V(Ω) is an interior vertex, then suppose additionally that Res z (1) = 0. Then holds for a constant C s that exclusively depends on the shape-regularity (and is in particular independent of the polynomial degree q).

Remark 1
The patch-wise construction of Q p in Subsection 4.2 typically generates local data oscillation osc k+p (ϕ z f, T (z)) in the error analysis as in the proof of Theorem 3 in Subsection 4.6 below or, e.g., [30,Theorem 3.17].A straightforward computation osc k+p (ϕ z f, T (z)) ≤ osc k+p−1 (f, T (z)) apparently leads to a loss of one degree in the data oscillation but Lemma 8 verifies for any p, q ∈ 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].

Local equivalence of stabilizations
The first improvement to the current HHO literature is the local equivalence of the two stabilizations sh from ( 37) and s h from (2).The authors of this paper could not find any motivation for the alternative stabilization sh 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

Theorem 4 (local equivalence of stabilizations)
The constants C 10 and C 11 exclusively depend on the polynomial degree k and the shape regularity of 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 , and the shape-regularity h F ≈ h T for all F ∈ F(T ) reveal Since Π 0 ϕ k = 0 (from the design of Rv h ), a Poincaré inequality shows On the one hand, an integration by parts provides On the other hand, an integration by parts and the definition of R imply Since ∆ϕ k ∈ P k (T ), the L 2 projection Π k on the right-hand side of ( 50) is redundant.Hence, the combination of ( 50)-( 51) with A Cauchy inequality on the right-hand side of (52), a discrete trace inequality, and , the combination of ( 48)-( 49) with (53) concludes the proof of (47).

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 ∈ N 0 , the solution u to (3) and the discrete solution u h to (18) The constant C 12 exclusively depends on k, p, and the shape regularity of T .
Proof Given k, p ∈ N 0 , let u ∈ V solve the Poisson model problem −∆ u = Π k+p f with the right-hand side Π k+p f .Subsection 4.3 in [32] constructs a stable enriching operator J : V h → 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 J : V h → V of the interpolation I : Consequently, the smoother J leads to a discrete solution u h = (u T , u F ) ∈ V h in the modified HHO discretization of [32] as a quasi-best approximation of the above solution ũ.The point is that u h ∈ V h coincides with the original HHO solution u h from (18).The arguments from the proof of [32,Theorem 4.10]

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 ∈ V resp.u h ∈ V h solve (3) resp.(18).Then Ru h and ARu h satisfy, for any p ∈ N 0 , that The constant C 13 exclusively depends on k, p, and the shape regularity of the triangulation T .
The proof can follow the proof of [25,Theorem 4.7] but additionally utilizes the two significant improvements from Subsections 4.3-4.4 that allow a stabilization-free efficiency in (54).

Remark 2 (p-robustness)
The H 1 (Ω)-conforming post-processing of Ru 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 It remains an open question whether this term can be controlled p-robustly by |||u−Ru h ||| pw +osc k (f, T ) (with a multiplicative constant independent of the polynomial degree k).

Proof of Theorem 3
Let p ∈ N 0 and r = 0 if k = 0 and r = k + p if k ≥ 1 as in (35) be given.Recall the abbreviation G h = ∇ pw Ru h ∈ P k (T ) with the discrete solution u h ∈ V h to (18) and let Q p = z∈V Q z,h ∈ H(div, Ω) denote the flux reconstruction from Subsection 4.2.The proof establishes (36) in five steps.
Step 1 provides the GUB |||u−Ru h ||| pw ≤ η eq,p (T ).This can follow from the paradigm of [1,2,30] as outlined below.The choice G := G h and w := ARu h in (8) and a triangle inequality lead to and the previously displayed formula result in the reliability |||u − Ru h ||| pw ≤ η eq,p (T ).
Step 3 reveals, for any q ∈ N 0 , the efficiency of the flux reconstruction for any polynomial degree k ≥ 1.The case k = 0 follows in Step 4 below.
Recall (39) for any vertex z ∈ V in the construction of Q p from Subsection 4.2 and set σ z := I RT (ϕ z G h ) ∈ RT pw k+p (T (z)).The piecewise Raviart-Thomas interpolation for any facet F ∈ F(T ) of a simplex T ∈ T and the normal trace γ T : This and elementary algebra with the product rule div pw (ϕ Recall the residual Res z (v) with v ∈ H 1 (Ω) for the vertex z ∈ V from (44).The commuting diagram property (57)-(58) establish the identity This, a piecewise integration by parts, and the property (20 = 0 for any interior vertex z ∈ V(Ω).Hence, Lemma 9 applies for any vertex z ∈ V and provides for the local contributions Q z,h of Q p = z∈V Q z,h from (41).The identity ) allows for a kand p-robust efficiency control of the flux reconstruction error with a constant C 14 that solely depends on the shape regularity of T .This is deemed noteworthy and motivates two different approaches for the bound of the residual Res z (v) on the right-hand side of (59) for p = 0 and p ≥ 1 below.
Step 3.1 provides (56 of the patch ω(z).This, the commuting diagram property (57), and (58) with ϕ z div pw G h ∈ P k (T (z)) in the definition of the residual Res z (v) from (44) verify The shape regularity of T , a Poincaré inequality for interior vertices z ∈ V(Ω), and a Friechrichs inequality for boundary vertices z ∈ V(∂Ω) provide Given any facet F ∈ F(z) ∩ F(Ω) in the facet spider F ∈ F(z) with facet patch ω(F ), a trace inequality thus shows Cauchy inequalities on the right-hand side of (61), the stability of the L 2 projection, ϕ z L ∞ (ωz) = 1, and (62)-(63) prove that Res z (v) is controlled by .
The efficiency of those residual terms follows from the proof of Theorem 2 in Section 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 ≥ 1.Given any normalized v ∈ H 1 * (ω(z)) with ∇v L 2 (ω(z)) = 1, (58) and the identity I RT (ϕ z G h ) = ϕ z G h in the definition of the residual Res z (v) from (44) verify that Res z (v) is equal to This, a piecewise integration by parts, and the product rule The weak formulation (3) with the test function This, a Cauchy-Schwarz, and a piecewise Poincaré inequality with the normalization ∇v L 2 (ω(z)) = 1 provide Since ∇ϕ z L ∞ (ω(z)) ≈ h −1 , the Leibniz rule, a triangle inequality, and (62) show that can be bounded by a constant independent of h.Thus, the combination of (43) with (59) and (65) results in the efficiency of Q p − G h in (60) with a kand p-robust constant C 14 .Since ϕ z Π k−1 f ∈ P k (T ) ⊂ P k+p (T ), the Pythagoras theorem and This and Lemma 8 implies the efficiency of the local oscillations from (65) for any q ∈ N 0 as in Remark 1.This and (60) prove (56) for p ≥ 1 and conclude the proof of Step 3.2.
Step 4 affirms the efficiency of the flux reconstruction for the polynomial degree k = 0 and any q ∈ N 0 .Let u ∈ V solve the Poisson model problem −∆ u = Π 0 f with piecewise constant right-hand side Π 0 f ∈ P 0 (T ).A careful inspection reveals that all arguments from Step 3 apply to the case k = 0 for f replaced by Π 0 f .This leads to Q p − G h ≤ C 15 ||| u − Ru h ||| pw with a constant C 15 that solely depends on the shape of T (because the data oscillations on the right-hand side of (60) vanish).Therefore, the triangle inequality ||| u − Ru h ||| pw ≤ |||u − u||| + |||u − Ru h ||| pw , the standard bound |||u − u||| ≤ C P osc 0 (f, T ), and C P = 1/π < 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 − Ru h ||| pw + osc k−1 (f, T ) η 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 ≥ 1 and in (66) for k = 0 imply the efficiency η eq,p |||u − Ru h ||| pw + osc q (f, T ) for any q ∈ N 0 .This concludes the proof of Theorem 3.

Numerical experiments
This section provides numerical evidence for optimal convergence and a comparison of the stabilization-free GUBs η res and η eq,p from the Sections 3 and 4 with the original error estimator η HHO from [25,Theorem 4.3] for the HHO method in three 2D benchmarks.The equilibrated GUB from Theorem 3, depends on the post-processed quantity Q ∆ p := Q p −∇ pw Ru h with the Raviart-Thomas function Q p ∈ RT k+p (T ) of degree k + p for p ∈ N 0 from Subsection 4.2.An algorithmic description of the computation of Q ∆ p for arbitrary p follows in Appendix A. Theorem 3 shows that η eq,p is efficient and reliable for all p ∈ N 0 .The residual-based error estimator η HHO from [25,Theorem 4.3] reads Proof Abbreviate ˜ (F ) := (n + 1)h F h 2 T |T | −1 for the facet F ∈ F(T ) of T and apply Lemma 3 to the singleton triangulation {T } and f −Π 0 f ∈ H 1 (T ).This and the Poincaré inequality The assertion (67) follows from the observation that ˜ (F ) is maximized on the facet F ∈ F(T ) with h F = h T .

Implementation and adaptive algorithm
Our implementation of the HHO method in MATLAB uses nodal bases for the spaces P k (F), P k (T ), and P k+1 (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 ∈ T , with the modified jump [•] F := • × n F along a boundary side F ∈ F(∂Ω) and the newest-vertex-bisection (NVB).The sum of this over all triangles is, up to some multiplicative constants, equivalent to η 2 res (T ).

Algorithm 1 AFEM algorithm
Input: Initial regular triangulation T 0 and polynomial degree k ∈ N 0 of the HHO method for levels := 0, 1, 2... do Solve (18) for discrete solution u ∈ V exactly on T and compute Ru Compute (refinement indicators) η 2 res (T ) for all T ∈ T Mark minimal subset M ⊂ T with 1 Figure 1 displays the energy norm |||e||| 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 η res , η eq,p , and η HHO are efficient and therefore equivalent to the energy error |||e||| pw , see Figure 1 on the right for k = 0, 2.

Corner singularity in the L-shaped domain
The third benchmark problem is set in the L-shaped domain Ω = (−1, 1) 2 \ [0, 1) 2 with constant right-hand side f ≡ 1 with an exact solution u ∈ H 1+s (Ω) for all 0 ≤ 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 for T ∈ T that are induced from the GUB η res , η HHO , and η eq,0 .Here, the norm |||e||| pw of the distance e from the discrete solution Ru u ∈ P k+1 ( T ) is the HHO approximation of u on an adaptive refinement T of T with at least 2|T | ≤ | T | elements.The three adaptive schemes (Algorithm 1, driven by η res (T ), η HHO (T ), or η eq,0 (T )) recover optimal rates of convergence and lead to similar local refinement of the adaptive mesh sequences as in Figure 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 η res , η HHO , and η 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 Figure 6.All three benchmarks verify that the considered error estimators are GUB with reliability constant 1, while the post-processing in the equilibrated GUB η eq,p produces minimal overestimation with significant additional computational costs.A.2 Degrees of freedom for RT pw q (T ) This subsection introduces a basis for the Raviart-Thomas finite element on T ∈ T (z) and starts with the definition of some linear functionals on H(div, T ).For any σ ∈ H(div, T ), set λ ,m T,div (σ) := T div σ x 1 x m 2 dx 0 ≤ + m ≤ q, λ j T,E (σ) := E σ |T • ν T s j ds 0 ≤ j ≤ q, E ∈ F (T ).
(The set Λ 0 T itself is linear independent [37, Lemma 3.1].)Given any Λ T with (72), denote the remaining Nq = q(q − 1)/2 degrees of freedom Λ T \ Λ 0 T by λ r T for r = 1, . . ., Nq.Let B RT,T = {ϕ j T,E , ϕ ,m T,div , ϕ r T } be the unique basis of RTq(T ) dual to Λ T with inferred indices from Λ T .Then, the collection B RT := T ∈T (z) B RT,T is a basis of RT pw q (T (z)) and any function σz ∈ RT pw q (T (z)) has the representation with coefficients c j T,E = λ j T,E (σz), c ,m T,div = λ ,m T,div (σz), and c r T = λ r T (σz) for all T ∈ T (z), E ∈ F (T ), and 0 ≤ j ≤ q, 1 ≤ + m ≤ q, 1 ≤ r ≤ Nq.By duality, the coefficients c j T,E with 0 ≤ j ≤ q uniquely determine the normal trace (σ z|T ) |E • ν E ∈ Pq(E) on the edge E ∈ F (T ) of T ∈ T (z).Any set of degrees of freedom Λ T with (72) works with the equilibration algorithm in A.3.
Example 3 (Construction of Λ T ) This example presents a generic procedure to obtain such a set from Λ T .The integration by parts formula shows that the lowest-order divergence moment λ 0,0 T, div = E∈F (T ) λ 0 T,E depends linearly on the lowest-order edge moments and, similarly, the sums λ ,m T,div + λ −1,m T,1 + mλ ,m−1 T,2 ∈ (P +m (F (T ))) * (74) are functionals on P +m (F (T )) (summands with negative indices are understood as zero).This relation allows for the substitution of volume moments in Λ T for divergence moments λ ,m T,div , 1 ≤ + m ≤ q, and leads to Λ T with (72).The remaining degrees of freedom Λ T \ Λ 0 T are volume moments of the form λ r T = λ ,m T,α for a fixed α ∈ {1, 2}, e.g.,
1: Initialize all coefficients c j T,E , c ,m T,div , c r T in (73) with zero.2: for a := 1 : N do 3: c 0 Ta,E a−1 Ta,E a−1
This subsection presents the equilibration procedure, starting with Algorithm 2, that computes an admissible function σ∆ z ∈ S(z) in terms of the representation (73).Enumerate the N := |T (z)| triangles T ∈ T (z) from 1 to N as in Figure 7. Any two neighbouring triangles Ta, T a+1 share an edge Ea := Ta ∩ T a+1 for a = 1, ..., N − 1.If z ∈ V(Ω) is an interior vertex, T 1 and T N share an additional edge E 0 := E N := T 1 ∩ T N .For a boundary vertex z ∈ V(∂Ω), T 1 , T N have the distinct boundary edges E 0 , E N ∈ F (z) ∩ F (∂Ω).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 from (69).Both proofs are provided in A.4.
Lemma 11 Given z ∈ V, let {ϕ j T,E , ϕ ,m T,div , ϕ r T } be the basis of RTq(T ) dual to Λ T with (72) for all T ∈ T (z).Suppose f ∈ L 2 (ω(z)) and G h ∈ Pq(T (z); R 2 ) satisfy Then the output of Algorithm 2 with input f and G h defines a function σ∆ z ∈ S(z).

Let b T ∈ S n+1 0 (
T ) with 0 ≤ b T ≤ 1 = max b T denote the volume bubblefunction on T ∈ T .The equivalence of norms in the finite-dimensional space T ).Hence, the decomposition of Lemma 1 and Lemma 5 result in |||u − Ru h ||| pw ≤ η res (T ).

Table 1
Explicit constants C 1 , Cst, C 2 for right-isosceles triangles with respect to the maximum interior angle ωmax of the polygonal domain Ω.