Guaranteed lower bounds on eigenvalues of elliptic operators with a hybrid high-order method

This paper introduces a novel hybrid high-order (HHO) method to approximate the eigenvalues of a symmetric compact differential operator. The HHO method combines two gradient reconstruction operators by means of a parameter 0<α<1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0<\alpha <~1$$\end{document} and introduces a novel cell-based stabilization operator weighted by a parameter 0<β<∞\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0<\beta <\infty $$\end{document}. Sufficient conditions on the parameters α\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} and β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document} are identified leading to a guaranteed lower bound property for the discrete eigenvalues. Moreover optimal convergence rates are established. Numerical studies for the Dirichlet eigenvalue problem of the Laplacian provide evidence for the superiority of the new lower eigenvalue bounds compared to previously available bounds.


Introduction
The eigenvalue problem for symmetric compact differential operators is a fundamental task in the numerical analysis with a well-understood a priori error analysis for conforming finite element methods (FEM) leading to optimal asymptotic convergence rates [4,5]. The Rayleigh-Ritz min-max principle shows that the discrete FEM eigenvalues are also guaranteed upper bounds of the exact eigenvalues, even in the pre-asymptotic range of coarse triangulations. In practice, guaranteed lower bounds (GLBs) can be even more important in a safety analysis in computational mechanics or for the detection of spectral gaps. The computation of lower eigenvalue bounds has been achieved based on the solution of nonconforming finite element schemes followed by a simple post-processing in [15,16]. In particular, letting κ 2 := π −2 + (2n(n + 1)(n + 2)) −1 in nD, and if λ CR ( j) is the j-th discrete eigenvalue computed with the Crouzeix-Raviart FEM, [16] proves (without extra conditions; the linear independency condition in [15,16] can be neglected [19,27]) that thereby delivering a GLB on the j-th continuous Dirichlet eigenvalue λ( j) for the Laplacian. Several other contributions [10,11,[27][28][29]34,37,38] derived GLBs using the maximal mesh-size h max as a global parameter as in (1.1). This results in a possible underestimation for locally refined triangulations. This motivated the design of novel schemes for the direct computation of GLBs without a global post-processing. The first example can be found in [20] for a hybridizable discontinuous Galerkin (HDG) method with Lehrenfeld-Schöberl stabilization, also studied under the label weak Galerkin scheme. The scheme proposed in [20] requires the fine tuning of skeletal-based stabilization terms. The main contribution of the present work is to introduce a new scheme which leads to GLBs on the eigenvalues under a simplified tuning of the stabilization terms. To achieve our goal, we rely on the framework of hybrid high-order (HHO) methods. HHO methods have been introduced in [23,24] for linear diffusion and locking-free linear elasticity and have been bridged to HDG and nonconforming virtual element methods (ncVEM) in [12]. Recall that in HHO methods the discrete unknowns are polynomials of degree k ≥ 0 attached to the mesh faces and polynomials of degree ∈ {k − 1, k, k + 1}, ≥ 0, attached to the mesh cells. In the present setting, we consider the degree = k + 1 for the cell unknowns. Moreover the two key ingredients in HHO methods are a local gradient reconstruction operator and a local stabilization operator. The HHO method devised herein introduces two novelties with respect to the literature. The first novelty is that the gradient reconstruction combines an operator mapping to piecewise Raviart-Thomas functions of degree k and an operator mapping to the piecewise gradient of piecewise polynomials of degree at most (k + 1). Although Raviart-Thomas reconstructions were considered previously in [1,21], it is the first time that they are combined with another reconstruction. Note that the present analysis cannot employ the single Raviart-Thomas reconstruction. The second novelty is that the stabilization operator is not skeletal-based but cell-based.
The weak formulation of the continuous Laplace eigenvalue problem seeks (λ, u) ∈ R + × H 1 0 ( ) such that a(u, v) = λb(u, v) for all v ∈ H 1 0 ( ) and b(u, u) = 1. (1. 2) The discrete eigenvalue problem seeks (λ h , u h ) ∈ R + × V h with While the bilinear form b h represents the L 2 ( ) scalar product b, the gradient-like approximations in the bilinear form a h involve two reconstructions R and G RT of the discrete unknowns in V h := P k+1 (T )×P k (F) with the space of piecewise polynomials of degree at most (k +1) on each simplex in the triangulation P k+1 (T ) and the space of piecewise polynomials of degree at most k on each face P k (F). The precise definition of the linear maps R : V h → L 2 ( ; R) and G RT : V h → L 2 ( ; R n ) can be found below in Sect. 3.1. Given two positive parameters 0 < α < 1 and 0 < β< ∞, for any u h := (u T , u F ) and v h := (v T , v F ) ∈ V h , the energy scalar product reads (1.4) Section 7.1.2 presents an algorithm for an effective parameter selection in the lowestorder case. The Raviart-Thomas reconstruction G RT of the gradient does not require a stabilization in the source problem [1,21], but the second term on the right-hand side in (1.4) has a negative sign and another stabilization with the reconstruction R in P k+1 (T ) is added. This paper investigates the a priori error analysis for discrete eigenvalues and eigenvectors and confirms optimal convergence rates. Moreover it shows that the following GLB property holds for the j-th discrete eigenvalue λ h ( j) of (1.3) and the j-th eigenvalue λ( j) of (1.2): with δ := κh max and σ , κ related to the stability of L 2 -projections onto piecewise polynomial spaces. (Notice that (1.5) means that each of the conditions (i) σ 2 β + δ 2 λ( j) ≤ α < 1 (a priori) and (ii) σ 2 β + δ 2 λ h ( j) ≤ α < 1 (a posteriori) implies the GLB property.) Numerical examples study the feasibility of the condition identified in the GLB (1.5) and the relation to (1.1) and the bounds in [20]. The remaining parts of this paper are organised as follows. After a short summary of the notation in Sect. 2, Sect. 3 introduces the new method (1.3) with all the necessary operators and the discrete bilinear forms a h and b h . Theorem 4.1 in Sect. 4 establishes (1.5). Section 5 contains the a priori error analysis which hinges on the Babuška-Osborn theory [4] and is inspired by [9] for the eigenvalue approximation by means of the standard HHO method (which does not have the GLB property). Section 6 concentrates on an alternative formulation of the lowest-order version for comparison with the Crouzeix-Raviart method. The numerical experiments in the final Sect. 7 illustrate the advantage of the direct lower bounds delivered by the present HHO method in the case of non-convex domains where adaptive mesh-refinement is necessary for optimal convergence rates. These results also provide numerical evidence for the superiority of the new GLBs compared to the aforementioned methods.

Triangulations
Let T denote a shape-regular triangulation of a bounded polyhedral Lipschitz domain ⊂ R n into closed n-simplices in the sense of Ciarlet [2,6,7,26]. For any simplex T ∈ T , let F(T ) denote the set of its (n + 1) sides and let N (T ) denote the set of its (n + 1) vertices. The intersection T 1 ∩ T 2 of two distinct, non-disjoint simplices T 1 and T 2 in T is the shared sub-simplex conv{N (T 1 ) ∩ N (T 2 )} = ∂ T 1 ∩ ∂ T 2 of their shared vertices. Furthermore, F := T ∈T F(T ) (resp. F( ) or F(∂ )) denotes the set of all (resp. interior or boundary) sides. For any simplex or sub-simplex K , let h K := diam(K ) denote its diameter. The piecewise constant function h T ∈ P 0 (T ) takes the value h T | T = h T on each simplex T ∈ T and h max := max T ∈T h T denotes the maximal mesh-size. Throughout this paper, ν T is the piecewise constant function which denotes for each simplex T ∈ T the outer unit normal vector ν T | T = ν T .

Discrete function spaces and L 2 -projections
For any M ∈ T or M ∈ F, ∈ N 0 , m ∈ N, let P (M; R m ) denote the set of polynomials of total degree at most in each component regarded as functions in L ∞ (M; R m ) and set and we omit R m whenever m = 1. The piecewise Raviart-Thomas space is RT pw (T ) := P (T )x + P (T ; R n ). The associated L 2 -projections are denoted The following two properties of the L 2 -projections are useful in the analysis below. These properties are classical (see, e.g., [

1)
and for all φ ∈ H m (T ; R n ), The constant C apx depends on the shape-regularity of T and on the polynomial degree , but is independent of the cell diameter h T .
The following refined stability estimates for L 2 -projections play an important role in the devising of guaranteed lower bounds on the eigenvalues.
Moreover there is κ > 0 such that for all T ∈ T and all f ∈ H 1 (T ), The constants C st and κ depend on the shape-regularity of T and on the polynomial degree , but are independent of the cell diameter h T . [20] where it is shown that, for any T ∈ T , the conditions

Remark 2.4 (Constants C st and κ) The stability estimate (2.3) is established in
are equivalent to the existence of an h T -independent constant C st (T ) > 0 such that (2.3) holds on T ∈ T , and one then sets 3) lead to (2.4) with κ ≤ C P C st . The Poincaré constant reads C P := 1/ j 11 for n = 2 with the first root of the first Bessel function j 11 [30] and C P ≤ 1/π for n ≥ 3 [3,32]. For = 0, an upper bound on the constant κ was first computed in [16] and improved in [15] for n = 2. The appendix of [20] proves κ 2 ≤ π −2 + (2n(n + 1)(n + 2)) −1 for any space dimension n.

Vector and matrix notation
For a, b ∈ R m×k , let a · b = a b ∈ R k×k and a ⊗ b = ab ∈ R m×m . The notation |•| depends on the context and denotes the Euclidean length of a vector, the cardinality of a finite set, the n-or (n − 1)-dimensional Lebesgue measure of a subset of R n . Furthermore a b abbreviates that there exists a generic constant C (independent of the mesh-size) with a ≤ Cb, whereas a ≈ b abbreviates a b a.

The modified HHO method
This section is devoted to the reconstruction and stability operators of the new method (1.3) and their properties. Moreover, the discrete bilinear forms are discussed and the abstract matrix eigenvalue problem is introduced.

Operators of interest
Let k ≥ 0 be the polynomial degree and set Define the reconstruction operator R : V h → P k+1 (T ) such that for any v h := Let G denote the Galerkin projection onto P k+1 (T ), i.e., for any φ ∈ H 1 (T ), we have Lemma 3.1 [23,24] The Raviart-Thomas reconstruction G RT : The comparison of (3 The stabilization operator S : (3.4) Lemma 3.3 Given any φ ∈ V = H 1 0 ( ) and any simplex T ∈ T , the stability term fulfils with σ 2 := C 2 st − 1 and the constant C st from Theorem 2.3.

Proof The definition (3.4) and Lemma 3.1 imply (3.5). Since
, the Pythagoras theorem proves for any T ∈ T , The first term is estimated in (2.3) by C 2 st ∇φ − k (∇φ) 2 L 2 (T ) , whereas the best approximation property of k proves that The combination finishes the proof of (3.6), and (3.7) follows by summation over T ∈ T .

Remark 3.4 (Piecewise evaluation) For all
, so that all these quantities can be computed piecewise and in parallel.

Discrete bilinear forms
Given global constants 0 < α < 1 and 0 < β< ∞, the reconstructions in (3.2)-(3.3), and the stabilization (3.4), the bilinear form a h : The L 2 -projection property of k and the definition of the piecewise bilinear form a pw imply that the bilinear form a h can be rewritten as follows for any The definiteness of • h is known and that of • a,h follows from Lemma 3.5 below. This guarantees that a h (•, •) is a scalar product on V h × V h which induces a norm on V h . Consequently the source problem associated with a h is well posed owing to the Lax-Milgram lemma.

Lemma 3.5
There exist constants γ , γ > 0, independent of the mesh-size, such that Proof The Pythagoras theorem and the best approximation property of F,k show for The trace inequality from [8,13] and [22,Eq. (1.42)] together with the Poincaré inequality and the shape-regularity of T lead to It is shown in [1,21] (see, e.g., Lemma 1 in [1]) that there is C 1 such that for any The combination of the preceding three displayed inequalities leads to an estimator of This and the triangle inequality show that The combination of (3.11)-(3.12) and the last identity from (3.10) shows that with the constant γ −1 := max{2(1 + C 2 (n + 1))/β, On the other hand, for any T ∈ T , (3.3) with q RT = G RT (v h ), an integration by parts, and the Cauchy-Schwarz inequality lead to A combination of a trace and an inverse estimate for Raviart-Thomas functions shows that h The sum over all T ∈ T reads The boundedness of the stability contribution follows from the triangle inequality, the estimate |R(v h ) | pw ≤ G RT (v h ) L 2 ( ) shown above and the last bound, leading to The second identity in (3.10) shows pw . This concludes the proof with γ := 2β + (1 + 2β)2 max{1, C 2 3 }.

Matrix eigenvalue problem
The algebraic realization of the eigenvalue problem (1.3) leads to a matrix eigenvalue problem with coefficient matrices and vector ( (3.13) The bilinear form b h solely depends on the volume components. The elimination of the face unknowns leads to the Schur complement and the equivalent matrix eigenvalue problem 14) The mass matrix B T T ∈ R dimP k+1 (T )×dimP k+1 (T ) is positive definite and allows the approximation of dimP k+1 (T ) eigenvalues and the application of the min-max principle (e.g. [33,Chapter 6]). The alternative formulation (3.14) will be exploited in Sect. 5 below.

Lower eigenvalue bounds
This section proves the main theorem, namely that the modified HHO method (1.3) provides guaranteed lower bounds for the continuous eigenvalues. Recall the constants C st and κ from Theorem 2.3, set σ 2 := C 2 st − 1, and δ := κh max . Theorem 4.1 Let λ( j) denote the j-th continuous eigenvalue of (1.2) and λ h ( j) the j-th discrete eigenvalue of (1.

3). Then each of the conditions
Proof To alleviate the notation we simply write λ h and λ.

Step 4. Lower bound for b h
Step 3, the estimate (2.4) and the Pythagoras theorem show that Since Step 5. Upper bound for a h (I (φ), I (φ)). Given φ ∈ H 1 0 ( ) from Step 3, the alternative form of a h in (3.8) and Lemma 3.2 prove The commuting property from Lemma 2.1, the Pythagoras theorem, and (4.5) The combination of (4.3)-(4.5) with (3.7) proves (for 0 < α < 1) that Step 6. Finish of the proof. The combination of Step 3-Step 5 shows that The pre-factor on the left-hand side is non-negative in the case where the assumption (ii) holds and (4.6) proves the claim λ h ≤ λ. In the case where the assumption (i) holds, (4.6) implies that For contradiction assume λ − λ h < 0 and divide the previous inequality by this difference so that 1 ≤ δ 2 (1 − k )(∇φ) 2 L 2 ( ) . This is smaller than or equal to This contradiction concludes the proof.

Convergence analysis
This section proves that the discrete eigenvalues converge with the expected optimal rates. The analysis hinges on the Babuška-Osborn theory for the spectral approximation of compact selfadjoint operators and adapts the arguments devised in [9] to analyze the spectral approximation using the standard HHO method.

Babuška-Osborn theory
Let H denote a Hilbert space with inner product ( If w n, j ∈ ker(μ n, j I − T n ) is a unit vector in the eigenspace of μ n, j for 1 ≤ j ≤ m, then there exists u ∈ E μ = ker(μI − T ) such that for all n ∈ N H ) . The constants C a and C b may depend on μ but are independent of n.

The source problem and relevant solution operators
Given a right-hand side f ∈ L 2 ( ), the weak formulation of the Poisson model problem and its solution is associated with the solution operator T : The source problem for the new method (1.3) seeks u h ∈ V h such that is associated with the solution operator T h : L 2 ( ) → V h with T h ( f ) := u h , and well-posed by Lemma 3.5. Using Lemma 3.5 and proceeding as in [9,Lemma 3.2] shows that the operator T h is bounded uniformly with respect to the mesh-size. The where the above right-hand side represents the consistency error. The solution property − u = φ a.e. in , and a piecewise integration by parts show, The last equality holds since v F = (v F ) F∈F ∈ P k (F( )) is single-valued for F ∈ F( ) and vanishes along F ⊂ ∂ , and v T = (v T ) T ∈T . The combination of the alternative form of a h from (3.8), Lemma 3.2, and Lemma 2.1 proves that

3) of G RT (v h ) and an integration by parts show that
The last three displayed equalities lead to For any T ∈ T , the combination of (2.1)-(2.2) and 0 < α < 1 proves that Hence, the Cauchy-Schwarz inequality shows that The term S 2 is controlled similarly, and the term S 3 is controlled via Lemma 3.3 and (2.1). The bounds on S 1 , S 2 , S 3 combined with the estimate (5.4) conclude the proof.

Error analysis for the eigenvalue problem
Based on Theorem 5.2, the difference T − T h in the summands in Theorem 5.1 are estimated as in [9] resulting in Theorem 5.6 below. Let us briefly outline the main steps. Recall the index s ∈ (1/2, 1] of the (reduced) elliptic regularity on the polyhedral domain ⊂ R n .
Proof This is [9, Lemma 4.1] for a symmetric a h (•, •) in combination with Theorem 5.2.
For a spectral value μ ∈ σ (T )\{0}, the smoothness of the functions in the eigenspaces E μ is quantified by t ∈ [s, k + 1] (depending on μ) and a constant C t such that Since possibly E μ ⊂ H s+ε ( ) for some ε > 0, we can have t > s > 1/2.
On E μ × E μ this bound can be improved.

Lemma 5.5
For any s ≤ t ≤ k + 1 and μ ∈ σ (T )\{0} verifying (5.6), any φ, ψ ∈ E μ satisfy Proof We detail the proof since it is slightly different from [9, Lemma 4.3]. For any φ, ψ ∈ E μ , the definition of T h and T h show that where we used that T h is the solution operator in (5.2). Using the symmetry of a h and again that T h is the solution operator in (5.2), we infer that To exploit the additional smoothness (5.6) of ψ ∈ H 1+t ( ) for S 4 , the projection property is combined with the Cauchy-Schwarz inequality and (2.1) to obtain The term S 5 needs a bit more algebraic manipulation. We first write The Cauchy-Schwarz inequality for a h , Theorem 5.2, and (5.6) prove for the first part that

With the abbreviation S 6 := a h (I T (φ), I T (ψ)) − a(T (φ), T (ψ)) the solution properties of T and T h show for the second summand of S 5 that
where the last inequality follows as in (5.7). The alternative form of a h from (3.8), the definition of a, as well as a combination of Lemma 3.2 and Lemma 2.1 lead to The Cauchy-Schwarz inequality, Lemma 3.3, and (2.1)-(2.2) prove for all φ, ψ ∈ E μ verifying (5.6) that The combination of (5.8)-(5.10) proves |S 5 | h 2(1+t) max φ L 2 ( ) ψ L 2 ( ) . This and (5.7) conclude the proof.

Given a discrete eigenfunction u
The constants C a , C b , C c > 0 may depend on μ, the polynomial degree k, the domain , and the mesh regularity but are independent of the mesh-size.
Proof The combination of Lemma 5.4-5.5 with the Babuška-Osborn theory in Theorem 5.1 provides the first two assertions. The proof for the final claim is analogue to Corollary 4.6 in [9]. In fact Lemma 5.5 provides the bound (5.10) for δ u := a h (I u, I u) − a(u, u) = S 6 .

Remark 5.7
The eigenvalues λ h in (1.3) and λ in (1.2) are associated with μ h = λ −1 h and μ = λ −1 , respectively, which leads to the same estimates. If the smoothness is optimal, i.e., t = k + 1, the proven convergence is of order h 2k+2 for the eigenvalues and h k+1 for the eigenvectors in the H 1 -seminorm.

Lowest-order case
This section is devoted to the analysis of the lowest-order case and provides a comparison to the Crouzeix-Raviart method. Throughout this section we set k := 0 so that V h := P 1 (T ) × P 0 (F( )).
Let mid(K ) := 1 m m j=1 P j ∈ K denote the barycenter of K := conv{P 1 , . . . , P m }. The associated piecewise constant function mid(T ) ∈ P 0 (T ; R n ) takes on each simplex T ∈ T the value mid(T )| T := mid(T ). Let us set

Comparison with the Crouzeix-Raviart method
The Crouzeix-Raviart (CR) finite element space reads The vector spaces P 0 (F( )) and CR 1 0 (T ) can be identified as follows. The extension operator I CR : In other words we have I NC = I CR • F ,0 . Recall the operators R from (3.2), G RT from (3.3), and S from (3.4).

Proof (a) For any p ∈ P 1 (T ) the definition (3.2) of R and an integration by parts
show that p).
concludes the proof of (a). (b) For any q RT ∈ RT 0 (T ), the definition (3.3) of G RT and an integration by parts show that On the other hand, since any q RT ∈ RT pw 0 (T ) can be written as the Pythagoras theorem for 0 and the definition of s(T ) show that The comparison of the last two displayed equalities shows that That completes the proof of (b). (c) This follows directly from the definition (3.4) and (a).

Proposition 6.2 (a) The bilinear forms in the lowest-order case read for all u h
The discrete solution u h := (u T , u F ) ∈ V h of the lowest-order EVP satisfies Proof (a) Using the alternative form for a h from (3.8) for the lowest-order case k = 0, the substitution of the operators with Lemma 6.1 proves the claim. 3) The choice 0 (v T ) ∈ P 0 (T ) in (6.3) shows that 1−α s(T ) 0 (u T − I CR (u F )) = λ h 0 (u T ) and equivalently Since for any v 1 ∈ P 1 (T ), we have v 1 = 0 (v 1 ) + ∇ pw v 1 · (• − mid(T )), the last displayed identity concludes the proof of (6.1) with Since the operator I CR : Proof The proof is similar to [20, Thm.6.2]. Since I CR F ,0 (v CR ) = v CR for any Crouzeix-Raviart function v CR ∈ CR 1 0 (T ), Proposition 6.2 shows that for v CR,h := The discrete min-max principles for λ h ( j) and λ CR ( j) conclude the proof.
The constant σ is more difficult to estimate (recall that σ := C 2 st − 1 but the proof of the existence of C st in [20, Thm 3.1] is non-constructive). Actually the constant σ is only needed in Lemma 3.3 to bound |S(I φ) | pw for all φ ∈ V := H 1 0 ( ). To bound this quantity in the lowest-order case, we can use an inverse inequality. Let us set for all T ∈ T , c inv (T ) := sup v 1 ∈P 1 (T ) ∇v 1 L 2 (T ) / v 1 L 2 (T ) . Notice that c inv (T ) is the square root of the maximal eigenvalue of the generalized matrix eigenvalue problem with the stiffness and mass matrix of P 1 (T ). Let us set c inv := max T ∈T c inv (T ). Recall from Remark 2.4 the Poincaré constant C P such that For any T ∈ T , the projection properties of 0 and 1 and the inverse estimate for affine functions in P 1 (T ) show that guaranteed lower bounds of [16,20] for the lowest-order variants on regular triangulations of the polygonal domains in Fig. 1 in 2D.

Implementation
The implementation is realized in MATLAB based on the data structure and assembling from [8,Section 7.8]. The resulting algebraic eigenvalue problem is solved with the MATLAB routine eigs exactly; the termination and round-off errors are neglected for simplicity. (An algebraic bound [31, 15.9.1] could partly circumvent the issue of inexact solve as in [16] -however, this bound solely clarifies the existence of an eigenvalue but gives no information on its numbering.) Given a regular triangulation T of the bounded polygonal Lipschitz domain ⊂ R 2 into triangles, the Crouzeix-Raviart eigenpairs (λ CR ( j), u CR ( j)) of (6.4) and the post-processed GLBCR( j) (1.1) from [16] are computed, together with the lowest-order version of the skeletal method (SM) [20] defined in (7.1) below. All these methods deliver guaranteed lower bounds which are compared to those delivered by the present modified HHO method. The SM method features the same discrete space V h = P 1 (T ) × P 0 (F( )) as the modified HHO method, whereas the discrete bilinear forms in the 2D case read (with κ defined in Remark 2.4 and I CR from Sect. 6.1) These quantities are computed and compared with the discrete eigenvalues λ h ( j) of the lowest-order modified HHO method (1.3) with the bilinear forms of Proposition 6.2.

Setting the parameters of the modified HHO method
The computable condition σ 2 β + δ 2 λ h ( j) ≤ α < 1 shows in Theorem 4.1 that the discrete eigenvalue λ h ( j) of the new method (1.3) is a lower bound for the exact eigenvalue λ( j) of (1.2). This condition restricts the choice of the parameters 0 < α < 1 and 0 < β < ∞. For right-isosceles triangles Lemma 6.4 shows that σ ≤ √ 72/ j 11 ≤ 2.2145 in (3.7) and δ ≤ κ h max in (2.4). The numerical bound κ ≤ 0.1893 from [27] slightly improves the analytical bound κ ≤ 0.2983 from [15] in the numerical experiments below. If λ CR ( j) denotes the j-th Crouzeix-Raviart eigenvalue and λ h ( j) the j-th eigenvalue of (1.3), Corollary 6.3 proves λ h ( j) ≤ λ CR ( j). This means that on a given triangulation T with known λ CR ( j), the parameter choice 0 < α < 1 and β = (α − δ 2 λ CR ( j))/σ 2 > 0 leads to a guaranteed lower bound λ h ( j) of the j-th exact eigenvalue λ( j). If the parameters are chosen beforehand, the condition σ 2 β +δ 2 λ h ( j) ≤ α may not be satisfied on a coarser mesh due to the impact of h max . In this case the computed value λ h ( j) is replaced by zero (which is an obvious guaranteed lower bound).

Adaptive mesh refinement
Adaptive mesh refinement may recover optimal convergence rates. For the related Crouzeix-Raviart adaptive finite element method (AFEM) driven by the estimator η, whose local contributions for any T ∈ T of area |T | read  [17] proves optimality for the principal eigenvalue of the CR-EVP. Since the a posteriori error analysis for the new modified HHO method and the skeletal method [20] is left open, the refinement indicator (7.2) drives adaptive mesh-refinement in the AFEM algorithm [14, Algorithm 2.2] with Dörfler marking for bulk parameter θ = 0.5 (and θ = 1 for uniform refinement) and newest-vertex bisection. This refinement preserves the interior angles in the triangulation.

Experiments on the L-shaped domain
On the non-convex L-shaped domain := (−1, 1) 2 \[0, 1) × (−1, 0], the principal eigenvalue λ(1) = 9.6397238389738806 is computed with a P 2 finite element method on uniformly refined triangulations with Aitken extrapolation. The associated eigenvector is apparently in H 1 0 ( ) \ H 2 ( ) resulting in the reduced convergence rate 0.8 for uniform mesh-refinement in Fig. 2a. As soon as the initial triangulation from Fig. 1a is refined three times, λ h (1) slightly improves the known bound λ SM (1) = GLBCR(1) on the uniform meshes. The adaptive mesh-refinement driven by the estimator (7.2) allows to recover the optimal convergence rate with the skeletal method in [20] and the modified HHO method. Remarkably, the modified HHO method with the parameter choice α = 0.4 and β = 0.07 convinces with sharper bounds. In [36] the multiple eigenvalue λ(103) = λ(104) = λ(105) = 50π 2 and the associated eigenfunction in C ∞ ( ) are presented. Figure 2b shows for these eigenvalues the optimal convergence rate of one with uniform mesh-refinement. The similar plots obtained with adaptive mesh-refinement are omitted for brevity. The parameter choice α = 0.4 and β = 0.06 guarantees the condition σ 2 β + δ 2 λ h ( j) ≤ α for λ(103) = λ(104) = λ(105) on coarser meshes. Figure 3 illustrates that the new method improves all known guaranteed lower bounds for the principal eigenvalue on the L-shaped domain on at least three times uniform refined meshes for an appropriate parameter choice and that the parameter range that leads to improvement grows with mesh-refinement, but the improvement is more impressive on the coarser triangulation. For higher eigenvalues the parameter studies show similar results.

Experiments on the slit domain
On the non-convex slit domain := (−1, 1) 2 \ ([0, 1) ×{0}), the principal eigenvalue λ(1) = 8.371330522443726 and the sixth λ(6) = 30.535991049204789 are approximated with a P 2 finite element method on uniformly refined meshes and Aitken extrapolation for comparison. The first and sixth eigenfunction on the non-convex slit domain are obviously in H 1 0 ( )\H 2 ( ). For uniform mesh-refinement, this leads to the reduced convergence rates 0.4 for the first in Fig. 4a and 0.6 for the sixth in Fig. 4b. The AFEM algorithm with bulk parameter θ = 0.5 driven by the estimator (7.2) allows to recover the optimal rates for both λ SM ( j) and λ h ( j), j ∈ {1, 6}, in Fig. 4. The parameter choice α = 0.4 and β = 0.07 allows to compute sharper bounds with the new method on finer meshes. The orange line illustrates that with the parameter choice α = 0.4 and β = (α − δ 2 λ CR ( j))/σ 2 , the discrete eigenvalue λ h ( j) is a guaranteed lower bound on each triangulation. For the slit domain, Fig. 5 illustrates that the new method improves all known guaranteed lower bounds on the moderately (three resp. four times for λ(1) and λ(6)) uniformly refined triangulation for an appropriate parameter choice with a wider range of appropriate parameters on a finer mesh.

Experiments on the isospectral domains
The isospectral drums with the initial triangulation of Fig. 1c, d have the same eigenvalues. The paper [25] displays approximations for the first 25 identical eigenvalues on these domains and [36] gives the approximation λ(50) = 54.187936. Figure 6 presents convergence plots for the error in the principal and the fiftieth eigenvalue. The first eigenfunction to the principal eigenvalue λ(1) = 2.53794399980 is in H 1 0 ( )\H 2 ( ) and leads to the reduced convergence rate 0.8 for uniform meshrefinement. The AFEM algorithm driven by (7.2) (after three uniform refinements to guarantee σ 2 β + δ 2 λ h ( j) ≤ α) recovers the optimal convergence rate for the direct lower bounds in Fig. 6a. In [36] are no remarks on the smoothness of the fiftieth eigenfunction. The numerical results displayed in Fig. 6b suggest that this eigenfunction is indeed in H 2 ( ). Uniform and adaptive mesh-refinement lead to optimal convergence rates. For brevity the adaptive results are not displayed. The parameter studies in Fig. 7 indicate that a parameter 0.1 ≤ α ≤ 0.5 and an appropriate β from Sect. 7.1.2 improve all known guaranteed lower bounds on refined meshes.

Conclusions
This subsection summarizes the empirical observations of the numerical experiments in Sects. 7.2-7.4.
(i) All experiments confirm the a priori convergence rates of Theorem 5.1. The convergence rate depends only on the smoothness of the approximated eigenfunction. For instance in Figs. 2b and 6b, the optimal convergence rate is one for uniform refinement despite the reduced convergence rate in Figs. 2a and 6a for the principal eigenvalue in H 1 0 ( )\H 2 ( ). (ii) Theorem 5.1 predicts a convergence provided the initial mesh is sufficiently fine.
In all examples the convergence rate is visible for moderate triangulations, so this restriction does not affect the numerical examples too much. illustrates that the scheme can be competitive to other methodologies for the computation of guaranteed lower eigenvalue bounds. For the appropriate parameter selection, the scheme can provide the sharpest bounds in comparison to [16,20]. Numerical benchmarks with the higher-order versions of the method suggested in the paper are even more promising provided the mesh is adapted appropriately and will be investigated in future research.
in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.