Adaptive hybrid high-order method for guaranteed lower eigenvalue bounds

The higher-order guaranteed lower eigenvalue bounds of the Laplacian in the recent work by Carstensen, Ern, and Puttkammer [Numer. Math. 149, 2021] require a parameter $C_{\mathrm{st},1}$ that is found $\textit{not}$ robust as the polynomial degree $p$ increases. This is related to the $H^1$ stability bound of the $L^2$ projection onto polynomials of degree at most $p$ and its growth $C_{\rm st, 1}\propto (p+1)^{1/2}$ as $p \to \infty$. A similar estimate for the Galerkin projection holds with a $p$-robust constant $C_{\mathrm{st},2}$ and $C_{\mathrm{st},2} \le 2$ for right-isosceles triangles. This paper utilizes the new inequality with the constant $C_{\mathrm{st},2}$ to design a modified hybrid high-order (HHO) eigensolver that directly computes guaranteed lower eigenvalue bounds under the idealized hypothesis of exact solve of the generalized algebraic eigenvalue problem and a mild explicit condition on the maximal mesh-size in the simplicial mesh. A key advance is a $p$-robust parameter selection. The analysis of the new method with a different fine-tuned volume stabilization allows for a priori quasi-best approximation and improved $L^2$ error estimates as well as a stabilization-free reliable and efficient a posteriori error control. The associated adaptive mesh-refining algorithm performs superior in computer benchmarks with striking numerical evidence for optimal higher empirical convergence rates.


Introduction
This paper proposes and analyzes a new hybrid high-order (HHO) eigensolver for the direct computation of guaranteed lower eigenvalue bounds (GLB) for the Laplacian.
1.1.Three categories of GLB.The min-max principle enables guaranteed upper eigenvalue bounds but prevents a direct computation of a GLB by a conforming approximation in a Rayleigh quotient.So GLB shall be based on nonconforming finite element methods (FEM), on modified mass and/or stiffness matrices (with reduced integration or fine-tuned stabilization terms), or on further post-processing.The last decade has seen a few GLB we group in three categories (i)-(iii).
(ii) Classical nonconforming FEM [20,22,45] and mixed FEM [39] allow for the GLB λ h /(1 + δλ h ) ≤ λ with the discrete eigenvalue λ h and a known parameter δ ∝ h 2 max in terms of the maximal mesh-size h max .On the positive side, the GLB provides unconditional information on the exact eigenvalue λ from the computed discrete eigenvalue λ h .On the negative side, the global parameter h max can spoil a very accurate approximation λ h in this GLB and is of lowest-order only.A finetuned stabilization of the classical nonconforming FEM in [23], however, provides a first (but low-order) remedy of the third category.
1.2.Motivation and outline of Section 2. The constants σ 2 1 := C 2 st,1 − 1 and κ := C P C st,1 in (1.1) depend on the Poincaré constant C P ≤ 1/π and a stability constant C st,1 .The latter has to be contrasted with the constant C st,2 , where C st,1 and C st,2 are the best possible constants in the stability estimates for all f ∈ H 1 (T ) (1.3) in a given simplex T ⊂ R n with the (component-wise) L 2 projection Π m and the Galerkin projection G m onto polynomials of total degree at most m ∈ N 0 .The two constants C st,1 and C st,2 are independent of the diameter h T := diam(T ) of T , but might depend on the shape of T and the polynomial degree p. Figure 1 illustrates the behaviour of C st,1 and C st,2 for different triangular shapes and various polynomial degrees p. Section 2 investigates the p-robustness of C st,2 and reveals that C st,2 ≤ C st,1 ∝ √ p + 1 tends to infinity as p → ∞, while we conjecture C st,2 ≤ √ 2 for triangles T with maximum interior angle ω ≤ π/2.Notice that a large constant C st,1 leads to a large σ 1 in (1.1) and so, α < 1 enforces small β and restricts the GLB to very fine meshes.The main motivation of this work arises from the convenient bound C st,2 ≤ √ 2: Can we design a discretization method of the third category (iii) based on σ 2 := C P C st,2 ≤ √ 2/π in (1.1)?
1.3.A modified HHO method and outline of Section 3.This paper provides an affirmative answer with a new fine-tuned stabilization in a modified HHO scheme in Section 3 and a new criterion sufficient for the GLB λ h ≤ λ.One advantage of (1.4) over (1.1) is the straightforward and p-robust parameter selection β := α/σ 2  2 .It turns out that σ 2 ≤ κ and so (1.4) improves on (1.1) in the sense that σ 2 2 h 2 max λ ≤ α holds on much coarser triangulations for higher polynomial degrees p.
Given a bounded polyhedral Lipschitz domain Ω ⊂ R n , let V := H 1 0 (Ω) denote the Sobolev space endowed with the energy scalar product a(u, v) := (∇u, ∇v) L 2 (Ω) and the L 2 scalar product b(u, v) := (u, v) L 2 (Ω) for all u, v ∈ V .This paper considers the model problem that seeks an eigenpair (λ, u) ∈ R + × (V \ {0}) such that a(u, v) = λb(u, v) for any v ∈ V. (1.5)The HHO methodology has been introduced in [32,31] and is related to HDG and nonconforming virtual element methods [27].Given a regular triangulation T into simplices, the ansatz space V h = P p+1 (T ) × P p (F(Ω)) consists of piecewise polynomials of (total) degree at most p + 1 attached to the simplices and piecewise polynomials of degree at most p attached to the interior faces.Two reconstruction operators link the two components of v h ∈ V h : The potential reconstruction Rv h ∈ P p+1 (T ) provides a discrete approximation to v in the space of piecewise polynomials P p+1 (T ) of degree at most p + 1.The gradient reconstruction Gv h ∈ RT pw p (T ) approximates the gradient ∇v in the space of piecewise Raviart-Thomas functions RT pw p (T ) [1,30].Let Sv h := v T − Rv h ∈ P p+1 (T ) for any v h = (v T , v F ) ∈ V h denote the additional cell-based stabilization.Given positive parameters 0 < α < 1 and 0 < β < ∞, the bilinear forms a h : The definitions of R, G, and further details follow in Section 3 below.1.4.GLB with p-robust parameters and outline of Section 4. The discrete bilinear form a h from [19] 6).The two stabilizations are locally equivalent, but the innovative difference is that the parameter selection in the new scheme circumvents an inverse inequality, and rather builds it into the scheme.Section 4 verifies the sufficient condition (1.4) for exact GLB under the assumption of exact solve.

1.5.
A priori error analysis of the new scheme and outline of Section 5. A quasi-best approximation for the source problem [38] allows for quasi-best approximation results in Theorem 5.1 for a simple eigenvalue λ, namely with a positive constant C 1 and the minimum 0 < s ≤ 1 of the index of elliptic regularity and one; the HHO interpolation I : V → V h is recalled in Subsection 3.3 below.Compared to earlier results in [12,19], (1.9) provides an additional positive power s of h max in the L 2 error.This is important as it eventually enables the absorption of higher-order terms in the a posteriori error analysis.
1.6.Stabilization-free a posteriori error analysis and outline of Section 6.
Let p h := Π p Gu h ∈ P p (T ; R n ) denote the L 2 projection of the gradient reconstruction Gu h ∈ RT pw p (T ) onto the space of vector-valued piecewise polynomials with the normal jump [p h • ν F ] F and the tangential jump [p h × ν F ] F of p h across a side F of T .Theorem 6.1 asserts reliability and efficiency of the error estimator η 2 := T ∈T η 2 (T ) for sufficiently small mesh-sizes h max in the sense that 1.7.Adaptive mesh-refining algorithm and outline of Section 7. Three 2D computer experiments in Section 7 provide striking numerical evidence that the criterion (1.4) indeed leads to confirmed lower eigenvalue bounds.The adaptive mesh-refining algorithm driven by the refinement indicator η from (1.10) recovers the optimal convergence rates of the eigenvalue error λ − λ h in all numerical benchmarks with singular eigenfunctions.This is the first time that p-robust higher-order GLB of the third category are displayed.
1.8.General notation.Standard notation for Lebesgue and Sobolev function spaces applies throughout this paper.In particular, (•, •) L 2 (ω) denotes the L 2 scalar product and H(div, ω) is the space of Sobolev functions with weak divergence in L 2 (ω) for a domain ω ⊂ R n .Recall the abbreviation V := H 1 0 (Ω) for the space of Sobolev functions, endowed with the energy scalar product a(u, v) := (∇u, ∇v) L 2 (Ω) and the L 2 scalar product b(u, v) := (u, v) L 2 (Ω) for all u, v ∈ V .
For a subset M ⊂ R n of diameter h M , let P p (M ) denote the space of polynomials of maximal (total) degree p regarded as functions defined in M .Given a simplex T ⊂ R n , the space of Raviart-Thomas finite element functions reads The Galerkin projection G := G p+1 : H 1 (T ) → P p+1 (T ) maps f ∈ H 1 (T ) to the unique solution Gf to Π 0 Gf = Π 0 f and (∇Gf, ∇p p+1 ) L 2 (T ) = (∇f, ∇p p+1 ) L 2 (T ) for all p p+1 ∈ P p+1 (T ) (1.12) with the convention In 2D, C P ≤ 1/j 11 = 0.260980 with the first positive root of the Bessel function J 1 [44] and C P ≤ 1/π in any space dimension [49,5].The context-sensitive notation | • | may denote the absolute value of a scalar, the Euclidean norm of a vector, the length of a side, or the volume of a simplex.The notation A ≲ B abbreviates A ≤ CB for a generic constant C independent of the mesh-size and A ≈ B abbreviates A ≲ B ≲ A. Throughout this paper, C 1 , . . ., C 14 denote positive constants independent of the mesh-size.

Stability estimates
This section discusses the behaviour of the constants C st,1 , C st,2 from (1.2)-(1.3)as p → ∞ and the computation of σ 2 := C P C st,2 with the Poincaré constant C P in (1.4) that arises from the stability estimates in Lemma 2.2 below.

2.1.
Stability constants and estimates.The following theorem asserts that C st,2 is p-robust (and small in general, see Figure 1) whereas C st,1 → ∞ as p → ∞.
Theorem 2.1.For any simplex T , there exist positive constants Proof.The existence of the constants 1 ≤ C st,1 ≤ C st,2 follows from [25, Theorem 3.1]; cf.Appendix A for further details.The technical proof of the p-robustness of C st,2 involves a linear bounded operator R curl : H −1 (T ; R 3 ) → L 2 (T ; R 3 ) from [42,28,46] and is carried out in Appendix B. The robustness holds for n = 2 with a simpler and hence omitted proof.The remaining parts of this proof concern the growth of C st,1 .Let X := H 1 (T )/R denote the Hilbert space with inner product (∇•, ∇•) L 2 (T ) and note that ker(∇(1 ≤ |||ϕ||| for every ϕ ∈ X, the definition of the operator norm of the oblique projection 1 − Π p+1 ∈ L(X; X) provides Kato's oblique projection lemma [52] for the Hilbert space X leads to t and Curl w := curl w for n = 3.For any g ∈ H −1 (T ; R 2n−3 ) in the dual space of H 1 0 (T ; R 2n−3 ) endowed with the operator norm The gradients ∇P p+1 (T ) of polynomials P p+1 (T ) of degree at most p + 1 form a subspace of P p (T ; R n ) and give rise to the L 2 orthogonal decomposition P p (T ; are defined, for any q p , r p ∈ Q p , by a(q p , r p ) := (q p , r p ) L 2 (T ) and b(q p , r p ) := ((−∆) −1 curl q p , curl r p ) L 2 (T ) .
(2.2) Theorem 2.3 (stability constant).The maximal eigenvalue of the eigenvalue problem a(q p , r p ) = λb(q p , r p ) for all r p ∈ Q p (2.4) Notice that (2.4) is a finite-dimensional eigenvalue problem and (−∆) −1 q p in b(q p , r p ) can be approximated by, e.g., a conforming FEM.Numerical experiments below even suggest that the bound C st,2 = m p is exact in n = 2 dimensions.
The remaining parts of the proof therefore assume p ≥ 1.Given f ∈ H 1 (T ), assume without loss of generality that ∇f ⊥ ∇P p+1 (T ) in ).Throughout this proof, abbreviate q p := Π p ∇f ∈ P p (T ; R n ).A Helmholtz decomposition leads to q p = ∇a + Curl b with a ∈ H 1 (T ) and b ∈ H 1 0 (T ; R 2n−3 ).For any v ∈ H 1 0 (T ; R 2n−3 ), the L 2 orthogonality Curl v ⊥ ∇a in L 2 (T ; R n ), an integration by parts, and a Cauchy inequality prove (The proof solely involves elementary algebra and is therefore omitted.)Hence, (2.5) implies ).Since ∇P p+1 (T ) ⊂ P p (T ; R n ), the best approximation of q p in ∇P p+1 (T ) satisfies the L 2 orthogonality q p ⊥ ∇P p+1 (T ).This and the Pythagoras theorem provide min vp+1∈Pp+1(T ) On the other hand, the constant m p from (2.3) satisfies min vp+1∈Pp+1(T ) Hence, (2.6) implies , a triangle inequality, the estimate 2(∇a, ∇(f − a)) L 2 (T ) ≤ δ|||a||| 2 + |||f − a||| 2 /δ, and (2.7) show, for all positive parameters δ > 0, that b(q p , q p ) and the orthogonal decomposition P p (T ; R n ) = Q p ⊕ ∇P p+1 (T ) with curl ∇P p+1 (T ) ≡ 0 reveal a(q p , q p )/b(q p , q p ) (2.9) with the bilinear forms a and b from (2.2).The min-max principle [4,Sec. 8] and (2.9) show that m 2 p is the maximal eigenvalue of (2.4).This concludes the proof.□ Example 2.4 (numerical example).Table 1 displays the computed maximal eigenvalue m 2 p ≥ C 2 st,2 of the eigenvalue problem (2.4) for the right-isosceles triangle T .The right-hand side is approximated by the Courant FEM of polynomial degree 10 on a uniform triangulation of T with 50721 degrees of freedom.The lower bounds sup for C st,2 and C st,1 from (1.2) are computable Rayleigh quotients and displayed in Figure 1.Computer experiments provide numerical evidence for the convergence of the lower bounds of C st,2 to m p as N → ∞ and, hence, for C st,2 = m p .The lower bound of C st,1 ∝ √ p + 1 displays the expected growth.
Undisplayed numerical experiments suggest that a small minimal interior angle does not affect the asymptotic bound of C st,2 , but leads to increased growth of C st,1 as p → ∞.We observed C st,2 = m p and the convergence C 2 st,2 → 2 as p → ∞ for different isosceles and various right triangles, whereas an interior angle ω > π/2 has a mild influence on the maximal value of C st,2 as shown for isosceles triangles

The modified HHO method
This section introduces the HHO method and the discrete eigenvalue problem.
, there exists a unique T ∈ T with F ∈ F(T ).Then ω F := int(T ), ν F := ν T is the exterior unit vector of F ∈ F(T ), and [v] F := v| F .The triangulation T gives rise to the space H 1 (T ) := {v ∈ L 2 (Ω) : v| T ∈ H 1 (T )} of piecewise Sobolev functions.The differential operators div pw , ∇ pw , and ∆ pw denote the piecewise applications of div, ∇, and ∆ without explicit reference to the triangulation T .
3.2.Discrete spaces.Let P p (T ), P p (F), and RT pw p (T ) denote the space of piecewise functions with restrictions to T ∈ T or F ∈ F in P p (T ), P p (F ), and RT p (T ).The local mesh-sizes give rise to the piecewise constant function h T ∈ P 0 (T ) with T ) onto P p (T ), P p (F), and RT pw p (T ) are computed cell-wise.For vector-valued functions τ ∈ L 1 (Ω; R n ), the L 2 projection Π p onto P p (T ; R n ) = P p (T ) n applies componentwise.The Pythagoras theorem implies the stability of L 2 projections, for any τ ∈ L 2 (Ω; R n ) and v ∈ L 2 (Ω), The Galerkin projection Gf of f ∈ H 1 (T ) is computed cell-wise by (1.12) with The inclusion denote the ansatz space of the HHO method for p ∈ N 0 .The interior sides F(Ω) give rise to the subspace P p (F(Ω)) of all (v F ) F ∈F ∈ P p (F) with the convention v F ≡ 0 on any boundary side F ∈ F(∂Ω) for homogenous boundary conditions.In other words, the notation [32, Eq. ( 28)] or [31, Eq. ( 41)] reads

Potential reconstruction. The potential reconstruction Rv
The bilinear form (∇ pw •, ∇ pw •) L 2 (Ω) on the left-hand side of (3.5) defines a scalar product and the right-hand side of (3.5) is a linear functional in the quotient space P p+1 (T )/P 0 (T ).The Riesz representation Rv h ∈ P p+1 (T ) of this linear functional in P p+1 (T )/P 0 (T ) is selected by The unique solution Rv h ∈ P p+1 (T ) to (3.5)-(3.6)defines the potential reconstruction operator R : V h → P p+1 (T ).

Gradient reconstruction. The gradient reconstruction Gv
In other words, Gv h is the Riesz representation of the linear functional on the righthand side of (3.7) in the Hilbert space RT pw p (T ) endowed with the L 2 scalar product.Since ∇ pw P p+1 (T ) ⊂ RT pw p (T ), (3.7) implies the L 2 orthogonality Gv h − Rv h ⊥ ∇ pw P p+1 (T ).The following lemma recalls the commutativity of G and R [32,31,1,33].The Galerkin projection G is defined in (1.12).
in the definition (1.6) of a h .Since the two stabilizations are locally equivalent, this leads to the assertion.□ The discrete eigenvalue problem (3.9) gives rise to the symmetric generalized algebraic eigenvalue problem The application of the Schur complement as in [19,Section 3.3]

A priori error analysis
The Babuška-Osborn theory [4] for the spectral approximation of compact selfadjoint operators leads to a priori convergence rates for the approximation of λ and of u in the energy norm [12,19].This section establishes the quasi-best approximation estimate (1.9) for a simple eigenvalue λ, that eventually allows for the absorption of higher-order terms in the a posteriori error analysis of Section 6.
Theorem 5.1 (a priori).If h max is sufficiently small, then (1.9) holds.The constant C 1 exclusively depends on p, n, Ω, and the shape regularity of T .

Lemma 5.2 (enriching operator).
There exists a linear bounded operator J : V h → V that is a right-inverse of I, i.e., v h = IJv h = (Π p+1 Jv h , Π p F Jv h ) for all v h ∈ V h , and stable in the sense that The constants ∥J∥ and C 4 solely depend on p, n, and the shape regularity of T .
Proof.The construction of the enriching operator J : V h → V in spirit of [53] involves standard averaging and bubble-function techniques from [54] and is explained in [38,Section 4.3] for a related HHO method without the proof of (5.1).Notice that J from [38] (called stabilized bubble smoother E H therein) only satisfies Jv h − v T ⊥ P p−1 (T ) for any given v h = (v T , v F ) ∈ V h .However, a straight-forward modification of [38,Eq. (4.16)] (in the notation of [38], B K v M ∈ P p+1 (K) should be defined by equation (4.16) therein for all q ∈ P p+1 (K)) immediately provides a right-inverse J of I.The arguments from [38,Propositions 4.5 and 4.7] apply and lead to the stability of J with respect to the equivalent discrete norm It remains to prove (5.1) which is well-known for the Crouzeix-Raviart finite element method with an appropriate interpolation I and the conforming companion J from [ with the L 2 projection Π p F onto P p (F ) for all faces F ∈ F. A triangle inequality, the stability of Π p F on a face F , and a discrete trace inequality show for all F ∈ F(T ) and T ∈ T .This, a triangle inequality for J := A + B(1 − A), (5.2), and the second inequality on [38, p. 2180] result in 3), the aforementioned inequalities, the trace inequality, and the piecewise application of the Poincaré inequality imply . This, the triangle inequality and the L 2 orthogonality ∇ pw (v − RIv) ⊥ ∇ pw P p+1 (T ) conclude the proof of (5.1).□ The second lemma proves quasi-best approximation estimates for a source problem.
and the data oscillation osc(f, The commutativity Π RT ∇v = GIv for v ∈ V from Lemma 3.1 enters this proof in two ways.First, it verifies Π p GI u = Π p ∇ u with v := u so that (3.8) reads Second, for v := J e h , the resulting L 2 orthogonality ∇J e h − G e h ⊥ RT pw p (T ) to the piecewise Raviart-Thomas functions RT pw p (T ) provides Since u ∈ V solves −∆ u = f in Ω, this and (5.6)-(5.7)verify with the Galerkin projection G from (1.12).Hence, the Poincaré inequality shows A Cauchy and a piecewise application of the Poincaré inequality reveal (5.10) The combination of (5.8)-(5.10)with a Cauchy inequality provides . This, (3.2)-(3.3), the stability ∥∇J e h ∥ L 2 (Ω) ≤ ∥J∥∥ e h ∥ a,h from Lemma 5.2, a Cauchy inequality, and The final lemma links (3.9) to (5.4).Recall the simple eigenpair (λ, u) of (1.5) and the associated discrete eigenpair (λ h , u h ) of (3.9) with with the constant Proof.This follows as in [21, Lem.2.4] with straight-forward modifications and is hence omitted.□ Proof of Theorem 5.1.The proof of (1.9) is split into three steps.
Step 1 provides the L 2 error estimate
Proof.This follows immediately from Theorem 5.1, the stability (1.3), and standard approximation properties of piecewise polynomials [11,Lemma 4.3.8]. □ The techniques of this section also apply to the HHO method of [19] and lead to the optimal rate s + t for the L 2 error towards a simple eigenvalue therein.
Theorem 6.2 (reliability and efficiency).For sufficiently small mesh-sizes h max , p h := Π p Gu h ∈ P p (T ; R n ) and η from (1.10) satisfy (1.11).The constants C eff and C rel exclusively depend on p, n, Ω, and the shape regularity of T .
Proof.The first part of the proof verifies that p h = Π p Gu h satisfies (A1)-(A2).

Numerical examples
The section presents three numerical benchmarks for the approximation of Dirichlet eigenvalues of the Laplacian on nonconvex domains Ω ⊂ R 2 .7.1.Parameter selection.For right-isosceles triangles, recall C st,2 ≤ √ 2 from Example 2.4 and C P = 1/( √ 2π) from [44].Throughout this section, let α = 0.5 and β := α/σ 2 2 = 4.934802 with σ 2 2 = C 2 P C 2 st,2 = 1/π 2 = 0.101321.The computable (a posteriori) condition σ 2 2 max{β, h 2 max λ h (j)} ≤ α from Theorem 4.1 leads to GLB(j) := λ h (j) ≤ λ(j).Since the parameters are chosen before-hand, the condition h 2 max λ h ≤ α/σ 2 2 = 4.934802 may not be satisfied on a coarse mesh with large h max and j.In this case, GLB(j) := 0 (which is a guaranteed lower eigenvalue bound), so only GLB are displayed in this section.7.2.Numerical realization.The algebraic eigenvalue problem (3.10) is realized with the iterative solver eigs from the MATLAB standard library in an extension of the data structures and short MATLAB programs in [3,17]; the termination and round-off errors are expected to be very small and neglected for simplicity.The a posteriori estimate from Theorem 6.1 motivates the refinement indicator η 2 (T ) from (1.10) with η 2 = T ∈T η 2 (T ).The standard adaptive algorithm [18, Algorithm 2.2] is modified in that, if h 2 max λ h ≤ α/σ 2 2 is not satisfied, the mesh is uniformly refined.It runs with the initial triangulations from Figure 3, the default bulk parameter θ = 0.5, and polynomial degrees p displayed in Figure 4.
The uniform and adaptive mesh-refinements lead to convergence history plots of the eigenvalue error λ(j) − GLB(j) and the a posteriori estimate η 2 plotted against the number of degrees of freedom of V h (ndof) in log-log plots below; dashed (resp.solid) lines represent uniform (resp.adaptive) mesh-refinements.7.3.L-shaped domain.The first example concerns the principle Dirichlet eigenvalue on the domain Ω := (−1, 1) 2 \ ([0, 1) × [0, −1)) with a re-entering corner at (0, 0) and the reference value λ(1) = 9.6397238440219410 from [9].This leads to the suboptimal convergence rate 2/3 for λ(1)−GLB(1) and η 2 (for all p) on uniform triangulations in Figure 5.The adaptive mesh-refining algorithm refines towards the origin as displayed in Figure 6 and recovers the optimal convergence rates p + 1 for λ(1) − GLB(1) and η 2 .7.4.Isospectral domain.The isospectral drums are pairs of non-isometric domains with identical spectrum of the Laplace operator.This subsection considers the domain Ω shown in Figure 3.b from [41]; the reference values λ(1) = 2.53794399980 and λ(25) = 29.5697729132are from [9] and [34].Figure 7 shows the suboptimal convergence rate 2/3 for λ(1)−GLB(1) and η 2 for the approximation of the principle eigenvalue λ(1) on uniformly refined triangulations.The adaptive mesh-refining algorithm refines towards four singular corners (for p = 3) as depicted in Figure 9 and recovers the optimal convergence rates p + 1 for λ(1) − GLB(1) and η 2 .Figure 8 displays the empirical convergence rate 1 for both λ(25) − GLB(25) and η 2 in case p = 0, while it is the expected rate 2/3 for p ≥ 1 in the presence of a typical corner singularity in the eigenfunction.We conjecture that the singular contribution to the corresponding eigenfunction in this particular example has a very small coefficient and the reduced asymptotic convergence rate 2/3 is therefore barely visible unless a very high accuracy is (e.g., absolute error in the eigenvalues much smaller than 5 × 10 −4 ).The adaptive mesh-refining algorithm refines towards four re-entering corners and recovers the optimal convergence rates p + 1 for λ(25) − GLB(25) and η 2 .There are two reasons for the plateau observed in the convergence history plot of λ(25) − GLB (25) in Figure 8.a.First, a larger pre-asymptotic range is expected and observed for the approximation of larger eigenvalues.Second, the condition h 2 max λ h ≤ α is not satisfied for the first triangulations, whence GLB is set to zero.An asymptotic behaviour is observed beyond 30,000 degrees of freedom for all displayed polynomial degrees p = 0, . . ., 4.   for p ≥ 1.Hence, there may be no reduction of the maximal mesh-size h max .Figure 11 displays suboptimal convergence rate 0.5 for the errors λ(1) − GLB(1) and η 2 for p = 0, . . ., 4. The adaptive mesh-refining recovers the optimal convergence rates p + 1.
7.6.Conclusions.The computer experiments provide empirical evidence for optimal convergence rates of the adaptive mesh-refining algorithm.The ad hoc choice α = 1/2 is robust in all computer experiments.For β = α/σ 2 2 , the computable condition σ 2 2 h 2 max λ h (j) ≤ α leads to confirmed lower eigenvalue bounds and holds on triangulations into right-isosceles triangles, whenever the maximal mesh-size h max satisfies λ h h 2 max ≤ απ 2 .In all displayed numerical benchmarks, λ h is a lower eigenvalue bound of λ even for λ h h 2 max > απ 2 .The computed (but otherwise undisplayed) efficiency indices 7×10 −2 ≤ I := |λ − λ h |η −2 ≤ 4×10 −3 range in the numerical examples from .07 to 4×10 −3 for an asymptotic range 2×10 4 ≤ ndof ≤ 10 5 ; the quotient I decreases for larger polynomial degree p.The overall numerical experience provides convincing evidence for the efficiency and reliability of the stabilization-free a posteriori error estimates of this paper.Higher polynomial degrees p lead to significantly more accurate lower bounds and clearly outperform lowest-order discretizations.
[55] T Wurzer.Stability of the trace of the polynomial L2-projection on triangles.Technical Report 36, Institute for Analysis and Scientific Computing, Vienna (2010).
Appendix: On p-robustness of constants in refined H 1 stability estimates This appendix provides details of the proof of Theorem 2.1 in the paper with focus on the constants C st,1 , C st,2 and their dependence on the polynomial degree p ∈ N 0 in three space dimensions.
Overview.Let |||•||| := ∥∇ • ∥ L 2 (T ) abbreviate the seminorm in the Sobolev space H 1 (T ) := H 1 (int(T )) and let Π p denote the L 2 projection onto the space P p (T ) of polynomials of total degree at most p ∈ N 0 for a fixed tetrahedron T ⊂ R 3 .For any Sobolev function f ∈ H 1 (T ), the Galerkin projection Gf ∈ P p+1 (T ) is the unique polynomial of degree at most p + 1 with the prescribed integral mean Π 0 Gf = Π 0 f and the orthogonality ∇(f − Gf ) ⊥ ∇P p+1 (T ) in L 2 (T ; R 3 ).The constants C st,1 and C st,2 are the best possible constants in the stability estimates Theorem 2.1 asserts the following properties of C st,1 and C st,2 .

3. 1 .
Triangulation.Let T be a regular triangulation of Ω into simplices in the sense of Ciarlet such that ∪ T ∈T T = Ω.Given a simplex T ∈ T of positive volume |T | > 0, let F(T ) denote the set of the n+1 hyperfaces of T , called sides of T .Define the set of all sides F = ∪ T ∈T F(T ) and the set of interior sides F(Ω) = F \{F ∈ F : F ⊂ ∂Ω} in T .For any interior side F ∈ F(Ω), there exist exactly two simplices T + , T − ∈ T such that ∂T + ∩ ∂T − = F .The orientation of the outer normal unit

1 Figure 3 . 4 Figure 4 .
Figure 3. Initial triangulations of the L-shaped domain in Subsection 7.3, the isospectral drum in Subsection 7.4, and the dumbbell-slit domain in Subsection 7.5.

Table 1 .
The constant C st,2 = m p on right-isosceles triangles.2.2.Numerical comparison and conjecture.The following theorem considers the computation of guaranteed upper bounds of C st,2 in n = 2, 3 space dimensions for a control of σ 2 in (2.1).Given and 2p U p 2for all p ≥ 2, whence √ p ≲ ∥Π p ∥ on the reference tetrahedron T .This and a scaling argument with an affine transformation concludes the proof for a general tetrahedron.□