Eric Cancès
email: cances@cermics.enpc.fr

Geneviève Dusson
email: dusson@ann.jussieu.fr

Yvon Maday
email: maday@ann.jussieu.fr

Benjamin Stamm
email: stamm@mathcces.rwth-aachen.de

Martin Vohralík
email: martin.vohralik@inria.fr

Guaranteed and robust a posteriori bounds for Laplace eigenvalues and eigenvectors: conforming approximations *

Keywords: Laplace eigenvalue problem, a posteriori estimate, eigenvalue error, eigenvector error, guaranteed bound, conforming finite element method, inexact solver, stopping criteria

This paper derives a posteriori error estimates for conforming numerical approximations of the Laplace eigenvalue problem with a homogeneous Dirichlet boundary condition. In particular, upper and lower bounds for an arbitrary simple eigenvalue are given. These bounds are guaranteed, fully computable, and converge with optimal speed to the given exact eigenvalue. They are valid without restrictions on the computational mesh or on the approximate eigenvector; we only need to assume that the approximate eigenvalue is separated from the surrounding smaller and larger exact ones, which can be checked in practice. Guaranteed, fully computable, optimally convergent, and polynomial-degree robust bounds on the energy error in the approximation of the associated eigenvector are derived as well, under the same hypotheses. Remarkably, there appears no unknown (solution-, regularity-, or polynomial-degree-dependent) constant in our theory, and no convexity/regularity assumption on the computational domain/exact eigenvector(s) is needed. Two improvements of the multiplicative constant appearing in our estimates are presented. First, it is reduced by a fixed factor under an explicit, a posteriori calculable condition on the mesh and on the approximate eigenvector-eigenvalue pair. Second, when an elliptic regularity assumption on the corresponding source problem is satisfied with known constants, the multiplicative factor can be brought to the optimal value of one. Inexact algebraic solvers are taken into account; the estimates are valid on each iteration and can serve for the design of adaptive stopping criteria. The application of our framework to conforming finite element approximations of arbitrary polynomial degree is provided, along with a numerical illustration on a set of test problems.

Introduction

Precise numerical approximation of eigenvalues and eigenvectors is crucial in countless applications. Thus, there has been a long-standing interest in answering the question: what is the size of the errors in computed eigenvalues and eigenvectors? This question is usually tackled via a posteriori error estimates. For elliptic source problems such as the Laplace one, conclusive answers are today given by, in particular, the theory of equilibrated fluxes following Prager and Synge [START_REF] Prager | Approximations in elasticity based on the concept of function space[END_REF], see Destuynder and Métivet [START_REF] Destuynder | Explicit error bounds in a conforming finite element method[END_REF], Braess et al. [START_REF] Braess | Equilibrated residual error estimates are p-robust[END_REF], Ern and Vohralík [START_REF]Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations[END_REF], and the references therein. The structure of the Laplace eigenvalue problem appears rather richer in comparison with the elliptic source case.

Recently, though, there has been an important progress in obtaining guaranteed lower bounds for the eigenvalues, especially for the first one: Luo et al. [START_REF] Luo | Computing the lower and upper bounds of Laplace eigenvalue problem: by combining conforming and nonconforming finite element methods[END_REF], Hu et al. [START_REF] Hu | Lower bounds for eigenvalues of elliptic operators: by nonconforming finite element methods[END_REF][START_REF] Hu | The lower/upper bound property of approximate eigenvalues by nonconforming finite element methods for elliptic operators[END_REF], Carstensen and Gedicke [16], Yang et al. [START_REF] Yang | The lower/upper bound property of the Crouzeix-Raviart element eigenvalues on adaptive meshes[END_REF], or Liu [START_REF] Liu | A framework of verified eigenvalue bounds for self-adjoint differential operators[END_REF] achieve so via the lowest-order nonconforming finite element method, Kuznetsov and Repin [START_REF] Kuznetsov | Guaranteed lower bounds of the smallest eigenvalues of elliptic differential operators[END_REF] and Šebestová and Vejchodský [START_REF] Šebestová | Two-sided bounds for eigenvalues of differential operators with applications to Friedrichs, Poincaré, trace, and similar constants[END_REF][START_REF]Two-sided bounds of eigenvalues -local efficiency and convergence of adaptive algorithm[END_REF] give numerical-method-independent estimates based on flux (functional) estimates, Liu and Oishi [START_REF] Liu | Verified eigenvalue evaluation for the Laplacian over polygonal domains of arbitrary shape[END_REF] elaborate fine a priori approximation estimates for lowestorder conforming finite elements, and, most recently, Xie et al. [START_REF] Xie | Fully computable error bounds for eigenvalue problem[END_REF] also rely on fluxes. Earlier work comprises Kato [START_REF] Kato | On the upper and lower bounds of eigenvalues[END_REF], Forsythe [START_REF] Forsythe | Asymptotic lower bounds for the fundamental frequency of convex membranes[END_REF], Weinberger [START_REF] Weinberger | Upper and lower bounds for eigenvalues by finite difference methods[END_REF], Bazley and Fox [START_REF] Bazley | Lower bounds for eigenvalues of Schrödinger's equation[END_REF], Fox and Rheinboldt [START_REF] Fox | Computational methods for determining lower bounds for eigenvalues of operators in Hilbert space[END_REF], Moler and Payne [START_REF] Moler | Bounds for eigenvalues and eigenvectors of symmetric operators[END_REF], Kuttler and Sigillito [START_REF] Kuttler | Bounding eigenvalues of elliptic operators[END_REF][START_REF]Estimating eigenvalues with a posteriori/a priori inequalities[END_REF], Still [START_REF] Still | Computable bounds for eigenvalues and eigenfunctions of elliptic differential operators[END_REF], Goerisch and He [START_REF] Goerisch | The determination of guaranteed bounds to eigenvalues with the use of variational methods. I[END_REF], Plum [START_REF] Plum | Guaranteed numerical bounds for eigenvalues[END_REF], Behnke et al. [START_REF] Behnke | Eigenvalue inclusions via domain decomposition[END_REF], and Armentano and Durán [START_REF] Armentano | Asymptotic lower bounds for eigenvalues by nonconforming finite element methods[END_REF], see also the references therein. Sometimes, though, restrictions may apply. A condition on relative closeness to the (first) eigenvalue is necessary in [START_REF] Kuznetsov | Guaranteed lower bounds of the smallest eigenvalues of elliptic differential operators[END_REF]Remark 3.2], [54, condition (3.6)], and [55, condition (5.23)] (in these references, the bounds actually do not converge with the correct speed); solution of an auxiliary eigenvalue problem for nonconvex domains is requested [START_REF] Liu | Verified eigenvalue evaluation for the Laplacian over polygonal domains of arbitrary shape[END_REF]; potential overestimation on adaptively generated meshes may hamper the bounds of [START_REF] Liu | Verified eigenvalue evaluation for the Laplacian over polygonal domains of arbitrary shape[END_REF]16,[START_REF] Liu | A framework of verified eigenvalue bounds for self-adjoint differential operators[END_REF], relying on a priori estimates and employing the largest mesh element diameter; an auxiliary global flux problem needs to be solved in [START_REF] Xie | Fully computable error bounds for eigenvalue problem[END_REF]; a saturation assumption may be necessary, see the discussion in [START_REF] Hu | Lower bounds for eigenvalues of elliptic operators: by nonconforming finite element methods[END_REF].

The question of precision for both eigenvalues and eigenvectors has also been investigated previously. For conforming finite elements, relying on the a priori error estimates resumed in Babuška and Osborn [START_REF] Babuška | Eigenvalue problems[END_REF] and Boffi [START_REF] Boffi | Finite element approximation of eigenvalue problems[END_REF], see also the references therein, a posteriori error estimates have been obtained by Verfürth [START_REF] Verfürth | Finite element discretizations of elliptic equations[END_REF], Maday and Patera [START_REF] Maday | Numerical analysis of a posteriori finite element bounds for linear functional outputs[END_REF], Larson [START_REF] Larson | A posteriori and a priori error analysis for finite element approximations of self-adjoint elliptic eigenvalue problems[END_REF], Heuveline and Rannacher [START_REF] Heuveline | A posteriori error control for finite approximations of elliptic eigenvalue problems[END_REF], Durán et al. [START_REF] Durán | A posteriori error estimates for the finite element[END_REF], Grubišić and Ovall [START_REF] Grubišić | On estimators for eigenvalue/eigenvector approximations[END_REF], Rannacher et al. [START_REF] Rannacher | Adaptive finite element solution of eigenvalue problems: balancing of discretization and iteration error[END_REF], and Šolín and Giani [START_REF] Solin | An iterative adaptive finite element method for elliptic eigenvalue problems[END_REF], see also the references therein. These estimates, though, systematically contain uncomputable terms, typically higher order on fine enough meshes.

Let Ω ⊂ R d , d = 2, 3, be a polygonal/polyhedral domain with a Lipschitz boundary, and let λ i , u i be the eigenvalues and associated eigenvectors of the Laplace operator -∆ on Ω with Dirichlet boundary conditions. The purpose of the present paper is to derive guaranteed and optimally convergent a posteriori bounds on both an arbitrary separated Laplace eigenvalue and the associated eigenvector for conforming (variational) methods. Nonconforming methods including nonconforming, discontinuous Galerkin, or mixed finite elements are treated in Cancès et al. [START_REF]Guaranteed and robust a posteriori bounds for Laplace eigenvalues and eigenvectors: a unified framework[END_REF]. We describe the setting in details in Section 2. Section 3 and Section 4 then contain a collection of equivalence inequalities between respectively the i-th eigenvalue error and the square of the i-th eigenvector energy error, the i-th eigenvector energy error and dual norm of the residual, and between the dual norm of the residual and its computable estimates. These results are valid under one key assumption: λ ih needs to be confined like λ i-1 < λ ih < λ i+1 (the left inequality of course only needs to hold when i > 1), see (5.2) below. This can be guaranteed in many cases of practical interest by a domain inclusion argument Ω -⊆ Ω ⊆ Ω + with known smaller and larger eigenvalues λ i-1 ≤ λ i-1 (Ω -) and λ i+1 (Ω + ) ≤ λ i+1 and by requesting λ i-1 =: λ i-1 (Ω -) < λ ih < λ i+1 := λ i+1 (Ω + ). Numerical bounds λ i-1 ≥ λ i-1 (typically available during the calculation) and λ i+1 ≤ λ i+1 (obtained a on a coarse mesh by the approach of [START_REF] Liu | Verified eigenvalue evaluation for the Laplacian over polygonal domains of arbitrary shape[END_REF]16,[START_REF] Liu | A framework of verified eigenvalue bounds for self-adjoint differential operators[END_REF]) can also be used, see Remarks 5.3 and 5.4 below. We also suppose that the approximation spaces consist of appropriate piecewise polynomials. For improved versions of our bounds, we additionally need to check the smallness of the L 2 (Ω)-norm of the Riesz representation of the residual, see the a posteriori calculable conditions (5.6) and (5.9) below. These can be always satisfied by refining the computational mesh/increasing the polynomial degree of the approximate solution. Note that no condition of Galerkin orthogonality of the residual to the finite element hat functions needs to be satisfied: the entire analysis is presented in the context of inexact algebraic solvers. Our estimates are valid on each iteration subject to the above inclusion of λ ih and can be used for efficient adaptive stopping criteria of iterative eigenvalue solvers, as promoted in, e.g., Mehrmann and Miedlar [START_REF] Mehrmann | Adaptive computation of smallest eigenvalues of self-adjoint elliptic partial differential equations[END_REF] or Carstensen and Gedicke [16].

In Section 5, the results of Section 3-Section 4 are turned into actual a posteriori bounds. First, upper and lower bounds for the i-th eigenvalue are given in Theorems 5.1 and 5.2. For a finite element approximation with an exact algebraic solver for simplicity, we obtain

λ ih -η 2 i ≤ λ i ≤ λ ih -η2 i , (1.1a)
with η i = m ih ∇u ih + σ ih,dis and ηi = ηi (r ih ) being fully computable quantities. Here u ih is the approxima-tion of the i-th exact eigenvector u i , • is the L 2 (Ω)-norm, σ ih,dis is an equilibrated flux reconstruction by mixed finite element local residual problems, and r ih is formed by conforming finite element local residual liftings. The associated eigenvector energy estimates are given next, with Theorem 5.6 revealing

∇(u i -u ih ) ≤ η i , η i ≤ C i ∇(u i -u ih ) , (1.1b)
where C i is a constant that only depends on λ 1 , λ i-1 , λ ih , λ i+1 , on the space dimension d, and on some Poincaré-Friedrichs-type constant C cont,PF together with a discrete stability constant C st , both only depending on the shape regularity of the mesh. In particular, C i is independent of the polynomial degree of u ih , leading to the polynomial-degree robustness. Moreover, a computable bound on C i is given. The constant C i , however, deteriorates for increasing eigenvalues. We distinguish three different cases. In Cases A and B of Theorems 5.1, 5.2, and 5.6, the multiplicative factor m ih of the estimator η i contains the factor max{(

λ ih λi-1 -1) -1 , (1 -λ ih λ i+1
) -1 } and similarly for ηi ; Case B improves the overall size of m ih under the fine-enough-mesh condition (5.6). The results of these two cases hold without any assumption on the convexity of the computational domain Ω and on the regularity of the weak solutions. If, additionally, elliptic regularity of the corresponding source problem is known, the interpolation and stability constants are computable (typically when d = 2 and Ω is convex), and the condition (5.9) holds, the factor m ih in front of the principal term ∇u ih + σ ih,dis has the optimal behavior 1 + O(h 2 ), as summarized in Case C of Theorems 5.1, 5.2, and 5.6.

We show how to apply the above general results to conforming finite elements of arbitrary order in Section 6. Numerical experiments presented in Section 7 fully support the theoretical findings; in particular the necessary conditions hold from quite coarse meshes. We only treat here simple eigenvalues and associated eigenvectors; clustered and multiple eigenvalues will be dealt with in a forthcoming contribution. Finally, building on these results, guaranteed error bounds and fully adaptive strategies with dynamic stopping criteria may become possible for nonlinear eigenvalue problems; some of our first results in this direction are summarized in [START_REF] Cancès | A perturbation-method-based a posteriori estimator for the planewave discretization of nonlinear Schrödinger equations[END_REF].

Setting

We denote by H 1 (Ω) the Sobolev space of L 2 (Ω) functions with weak gradients in [L 2 (Ω)] d and by V := H 1 0 (Ω) its zero-trace subspace. Similarly, H(div, Ω) stands for the space of [L 2 (Ω)] d functions with weak divergences in L 2 (Ω). The notations ∇ and ∇• are used respectively for the weak gradient and divergence. Moreover, for ω ⊂ Ω, (∇u, ∇v) ω stands for ω ∇u•∇v dx and (u, v) ω for ω uv dx; we also denote ∇v 2 ω := ω |∇v| 2 dx and v 2 ω := ω v 2 dx and drop the index whenever ω = Ω.

The Laplace eigenvalue problem

We consider here the problem: find eigenvector and eigenvalue pairs (u k , λ k ), with u k satisfying a homogeneous Dirichlet boundary condition over ∂Ω and subject to the constraint

u k = 1, such that -∆u k = λ k u k in Ω. In a weak form, (u k , λ k ) ∈ V × R + with u k = 1 and (∇u k , ∇v) = λ k (u k , v) ∀v ∈ V. (2.1)
Actually, cf. Gilbarg and Trudinger [START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF], Babuška and Osborn [START_REF] Babuška | Eigenvalue problems[END_REF], Boffi [START_REF] Boffi | Finite element approximation of eigenvalue problems[END_REF], or Strang and Fix [START_REF] Strang | An analysis of the finite element method[END_REF], u k , k ≥ 1, form a countable orthonormal basis of L 2 (Ω) consisting of vectors from V , whereas 0 < λ 1 < λ 2 ≤ λ 3 ≤ . . . going to +∞. The smallest eigenvalue λ 1 is positive and simple and the associated eigenvector u k to each simple λ k is unique up to the sign that we fix here by the condition (u k , χ k ) > 0, where χ k ∈ L 2 (Ω) is typically a characteristic function of Ω (for k = 1) or of its subdomain (for k > 1). Note that it follows from (2.1) and the scaling u k = 1 that ∇u k 2 = λ k . Below, we shall often employ the Parseval identity, giving for any v ∈ L 2 (Ω)

v 2 = k≥1 (v, u k ) 2 .
(2.2)

As (u k / √ λ k ) k≥1 form an orthonormal basis of V , for which one in particular uses that (∇u k , ∇u l ) = λ k (u k , u l ) = 0 for k = l, for any v ∈ V , we also obtain

∇v 2 = k≥1 (∇v, ∇u k ) 2 λ k = k≥1 λ k (v, u k ) 2 . (2.3)

Residual and its dual norm

The derivation of a posteriori error estimates usually exploits the concept of the residual and of its dual norm. We will proceed in this way as well. Let V stand for the dual of V .

Definition 2.1 (Residual and its dual norm). For any pair

(u ih , λ ih ) ∈ V × R, define the residual Res(u ih , λ ih ) ∈ V by Res(u ih , λ ih ), v V ,V := λ ih (u ih , v) -(∇u ih , ∇v) ∀v ∈ V.
(2.4a)

Its dual norm is then Res(u ih , λ ih ) -1 := sup v∈V, ∇v =1 Res(u ih , λ ih ), v V ,V . (2.4b)
We will also often work with the Riesz representation of the residual r (ih) ∈ V ,

(∇r (ih) , ∇v) = Res(u ih , λ ih ), v V ,V ∀v ∈ V, (2.5a
)

∇r (ih) = Res(u ih , λ ih ) -1 .
(2.5b)

Generic equivalences

In extension of some classical results, see [START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF][START_REF] Babuška | Eigenvalue problems[END_REF][START_REF] Boffi | Finite element approximation of eigenvalue problems[END_REF][START_REF] Strang | An analysis of the finite element method[END_REF], we establish in this section generic equivalence results between the following three quantities: the i-th eigenvalue error ∇u ih 2 -λ i , which can potentially be negative, the square of the i-th eigenvector energy error ∇(u i -u ih ) 2 , and the square of the dual norm of the residual Res(u ih , λ ih ) 2 -1 . These equivalences may for the moment contain uncomputable terms like the eigenvalues λ i-1 , λ i , λ i+1 or the Riesz representation norm r (ih) , but all such terms will be removed later. To proceed in an abstract way allowing for inexact algebraic solvers, we rather work with the eigenvalue error given by ∇u ih 2 -λ i instead of λ ih -λ i ; of course these coincide when the discrete Rayleigh quotient link ∇u ih 2 = λ ih holds, typically upon solver convergence. A generalization to any self-adjoint operator with compact resolvent can be found in Cancès et al. [START_REF]Guaranteed and robust a posteriori bounds for Laplace eigenvalues and eigenvectors: a unified framework[END_REF].

Our first two lemmas are similar in parts to the developments in [START_REF] Kuttler | Bounding eigenvalues of elliptic operators[END_REF][START_REF] Kuznetsov | Guaranteed lower bounds of the smallest eigenvalues of elliptic differential operators[END_REF][START_REF] Šebestová | Two-sided bounds for eigenvalues of differential operators with applications to Friedrichs, Poincaré, trace, and similar constants[END_REF][START_REF]Two-sided bounds of eigenvalues -local efficiency and convergence of adaptive algorithm[END_REF], giving a computable bound on the L 2 (Ω) error u i -u ih . Let i ≥ 1 and define

C ih := min 1 - λ ih λ i-1 2 , 1 - λ ih λ i+1 2 . ( 3.1)
The left term needs to be disregarded for i = 1.

Lemma 3.1 (L 2 (Ω) bound via a quadratic residual inequality). Let (u ih , λ ih ) ∈ V × R + with u ih = 1 and (u i , u ih ) ≥ 0 be the i-th approximate eigenvector-eigenvalue pair, i ≥ 1. Let λ i be simple and let λ i-1 < λ ih when i > 1 and λ ih < λ i+1 . Then

u i -u ih ≤ α ih := √ 2C -1 2 ih r (ih) . (3.2)
Proof. Characterizations (2.1), (2.4a), and (2.5a) give

(r (ih) , u k ) = (∇u k , ∇r (ih) ) λ k = (λ ih (u ih , u k ) -(∇u ih , ∇u k )) λ k = λ ih λ k -1 (u ih , u k ). (3.3)
Consequently, the Parseval equality (2.2)

with v = r (ih) yields r (ih) 2 = k≥1 (r (ih) , u k ) 2 = k≥1 1 - λ ih λ k 2 (u ih , u k ) 2 . (3.4)
Observe that the function x ∈ R + → 1 -λ ih x 2 reaches its minimum at x = λ ih and is decreasing on (0, λ ih ] and increasing on [λ ih , ∞). Thus the constant C ih in (3.1) equals min k≥1, k =i 1 -λ ih λ k

2

. Further, employing the scalings u i = 1 and u ih = 1,

(u ih -u i , u i ) = (u ih , u i ) -u i 2 = (u ih , u i ) - u i 2 2 - u ih 2 2 = - 1 2 u i -u ih 2 . (3.5)
As (u i , u k ) = 0 for k ≥ 1, k = i from the orthogonality of u k , elaborating (3.4) further while adding and subtracting C ih (u ih -u i , u i ) 2 and using (3.1) and (3.5) gives

r (ih) 2 = λ ih λ i -1 2 (u ih , u i ) 2 + k≥1, k =i 1 - λ ih λ k 2 (u ih -u i , u k ) 2 ≥ λ ih λ i -1 2 (u ih , u i ) 2 + C ih k≥1 (u ih -u i , u k ) 2 -C ih (u ih -u i , u i ) 2 (3.6) = λ ih λ i -1 2 (u ih , u i ) 2 + C ih u i -u ih 2 - C ih 4 u i -u ih 4 ,
where we have also employed (2.2) with v = u ih -u i . Dropping the first (nonnegative and presumably small) term on the right-hand side and denoting e ih := u i -u ih 2 , we conclude the validity of the quadratic residual inequality in e ih C ih 4 e 2 ih -C ih e ih + r (ih) 2 ≥ 0.

(

From the sign assumption (u i , u ih ) ≥ 0, employing

u i = u ih = 1, e ih = u i -u ih 2 = 2 -2(u i , u ih ) ≤ 2, (3.8)
so that C ih e ih ≤ 2 r (ih) 2 , i.e., (3.2). Note that inspecting more closely the quadratic inequality (3.7), the improved bound e ih ≤ 2 -4 -2α 2 ih ( √ 2-times better for e ih approaching zero) follows under condition r (ih) 2 < C ih that we prefer to avoid.

In addition to (3.1), define also (disregarding again the left term for i = 1)

C ih := min λ i-1 1 - λ ih λ i-1 2 , λ i+1 1 - λ ih λ i+1 2 .
(3.9) Lemma 3.2 (L 2 (Ω) bound with respect to ∇r (ih) ). Under the assumptions of Lemma 3.1, there also holds Recall the sign characterization (u i , χ i ) > 0 with χ i ∈ L 2 (Ω), i ≥ 1. The sign condition (u i , u ih ) ≥ 0 necessary in Lemmas 3.1 and 3.2 is typically always satisfied; the following lemma can be used for its rigorous verification: Lemma 3.3 (Sign verification). Let (u ih , λ ih ) ∈ V × R + satisfy u ih = 1, (u ih , χ i ) > 0, λ i-1 < λ ih when i > 1 and λ ih < λ i+1 , and α ih ≤ χ i -1 (u ih , χ i ) for α ih given by (3.2) or (3.10). Then the sign condition (u i , u ih ) ≥ 0 is satisfied.

u i -u ih ≤ α ih := √ 2 C -1 2 ih ∇r (ih) . (3.10) Proof. Developing (2.3) for v = r (ih) via (3.3) gives ∇r (ih) 2 = k≥1 λ k (r (ih) , u k ) 2 = k≥1 λ k 1 - λ ih λ k 2 (u ih , u k ) 2 . (3.11) Next, min k≥1, k =i λ k 1 -λ ih
Proof. Suppose -(u i , u ih ) > 0. Then the bounds of Lemmas 3.1 and 3.2 hold for -u ih in place of u ih , i.e., u i + u ih ≤ α ih . Consequently, a contradiction follows,

(u ih , χ i ) = -(u i , χ i ) + (u i + u ih , χ i ) < (u i + u ih , χ i ) ≤ u i + u ih χ i ≤ (u ih , χ i ).

i-th eigenvalue error equivalences

We first show how to exploit the L 2 (Ω) bound for equivalence between the eigenvalue error and the eigenvector error.

Theorem 3.4 (Eigenvalue bounds). Let u ih ∈ V with u ih = 1, i ≥ 1, be arbitrary subject to u i -u ih ≤ α ih for some α ih ∈ R + . Then ∇(u i -u ih ) 2 -λ i α 2 ih ≤ ∇u ih 2 -λ i ≤ ∇(u i -u ih ) 2 .
(3.12)

Under the additional assumption α 2 1h ≤ 2, there also holds, for the first eigenpair,

1 2 1 - λ 1 λ 2 1 - α 2 1h 4 ∇(u 1 -u 1h ) 2 ≤ ∇u 1h 2 -λ 1 . (3.13)
Proof. Using the weak solution characterization (2.1) and (3.5),

∇u ih 2 -λ i = ∇(u ih -u i ) 2 + 2(∇(u ih -u i ), ∇u i ) = ∇(u ih -u i ) 2 + 2λ i (u i , u ih -u i ) = ∇(u ih -u i ) 2 -λ i u i -u ih 2 .
(3.14)

Dropping the (nonpositive and presumably small) last term, the upper bound in (3.12) follows; estimating it using u i -u ih ≤ α ih , we arrive at the lower bound in (3.12). The bound (3.13) only seems to hold for the first eigenpair. To prove it, we use (2.2)

-(2.3) for v = u 1 -u 1h . First, ∇(u 1 -u 1h ) 2 -λ 1 u 1 -u 1h 2 = k≥1 (λ k -λ 1 )(u 1 -u 1h , u k ) 2 = k≥2 (λ k -λ 1 )(u 1 -u 1h , u k ) 2 . (3.15) Using λ k ≥ λ 2 for k ≥ 2, λ 2 > λ 1 , (3.5
) for i = 1, and the Cauchy-Schwarz inequality,

∇(u 1 -u 1h ) 2 -λ 1 u 1 -u 1h 2 ≥ (λ 2 -λ 1 ) k≥1 (u 1 -u 1h , u k ) 2 -(λ 2 -λ 1 )(u 1 -u 1h , u 1 ) 2 = (λ 2 -λ 1 ) u 1 -u 1h 2 - λ 2 -λ 1 4 u 1 -u 1h 4 .
Using u 1 -u 1h ≤ α 1h and reemploying (2.2) for v = u 1 -u 1h , we arrive at, second,

∇(u 1 -u 1h ) 2 -λ 1 u 1 -u 1h 2 ≥ (λ 2 -λ 1 ) u 1 -u 1h 2 -α 2 1h λ 2 -λ 1 4 u 1 -u 1h 2 = k≥1 (λ 2 -λ 1 ) 1 - α 2 1h 4 (u 1 -u 1h , u k ) 2 .
Summing this with (3.15) with weights 1 2 yields

∇(u 1 -u 1h ) 2 -λ 1 u 1 -u 1h 2 ≥ k≥1 λ k -λ 1 2 + λ 2 -λ 1 2 1 - α 2 1h 4 (u 1 -u 1h , u k ) 2 . Now notice that, using (2.3) for v = u 1 -u 1h , 1 2 1 - λ 1 λ 2 1 - α 2 1h 4 ∇(u 1 -u 1h ) 2 = k≥1 λ k 2 1 - λ 1 λ 2 1 - α 2 1h 4 (u 1 -u 1h , u k ) 2 .
A simple calculation (note 1 2 ≤ 1 -

α 2 1h 4 ≤ 1) shows that λ k -λ 1 2 + λ 2 -λ 1 2 1 - α 2 1h 4 ≥ λ k 2 1 - λ 1 λ 2 1 - α 2 1h 4 k ≥ 1.
Thus

∇(u 1 -u 1h ) 2 -λ 1 u 1 -u 1h 2 ≥ 1 2 1 - λ 1 λ 2 1 - α 2 1h 4 ∇(u 1 -u 1h ) 2 ,
and (3.13) follows using (3.14).

i-th eigenvector error equivalences

We next investigate the equivalence between the eigenvector error ∇(u i -u ih ) and the dual norm of the residual Res(u ih , λ ih ) -1 . Recall the definition (3.1) and also set

C ih := 1 if i = 1, C ih := max λ ih λ 1 -1 2 , 1 if i > 1. (3.16)
Furthermore, let

γ ih := ∇(u i -u ih ) 2 if λ i ≤ ∇u ih 2 is known to hold, max{ ∇(u i -u ih ) 2 , λ i α 2 ih } otherwise;
(3.17) we refer to Remark 5.4 below for the discussion when λ i ≤ ∇u ih 2 .

Theorem 3.5 (Eigenvector bounds).

Let (u ih , λ ih ) ∈ V × R + with u ih = 1, i ≥ 1, be arbitrary subject to u i -u ih ≤ α ih for some α ih ∈ R + . Then ∇(u i -u ih ) 2 ≤ Res(u ih , λ ih ) 2 -1 + (λ ih + λ i )α 2 ih , (3.18a
)

Res(u ih , λ ih ) 2 -1 ≤ λ ih -∇u ih 2 + γ ih 2 λ i + C ih ∇(u i -u ih ) 2 . (3.18b) Let λ i-1 < λ ih when i > 1, λ ih < λ i+1
, and α 2 ih ≤ 2 λ1 λi . Then there also holds

∇(u i -u ih ) 2 ≤ C -1 ih 1 - λ i λ 1 α 2 ih 4 -1 Res(u ih , λ ih ) 2 -1 . (3.19)
Proof. Starting from (3.11), adding and subtracting

C ih λ i (u ih -u i , u i ) 2 , using (u i , u k ) = 0 for k ≥ 1, k = i, (3.5)
, and the Cauchy-Schwarz inequality, we observe

∇r (ih) 2 ≥ λ i λ ih λ i -1 2 (u ih , u i ) 2 + C ih k≥1 λ k (u ih -u i , u k ) 2 -C ih λ i (u ih -u i , u i ) 2 = λ i λ ih λ i -1 2 (u ih , u i ) 2 + C ih ∇(u i -u ih ) 2 - C ih 4 λ i u i -u ih 4 ≥ C ih ∇(u i -u ih ) 2 - C ih 4 λ i u i -u ih 4 .
Using the Poincaré-Friedrichs inequality

u i -u ih 2 ≤ 1 λ1 ∇(u i -u ih ) 2 , ∇r (ih) 2 ≥ C ih ∇(u i -u ih ) 2 - C ih 4 λ i λ 1 ∇(u i -u ih ) 2 α 2 ih ,
where we have also employed u i -u ih ≤ α ih . Thus (3.19) follows via (2.5b).

The proof of Lemma 3. (3.16). Thus, (3.11) together with the Cauchy-Schwarz inequality and u i = u ih = 1 give

1 gives sup k≥1, k =i 1 -λ ih λ k 2 = C ih , recalling
∇r (ih) 2 ≤ λ i λ ih λ i -1 2 + C ih k≥1, k =i λ k (u ih -u i , u k ) 2 ≤ (λ ih -λ i ) 2 λ i + C ih ∇(u i -u ih ) 2 .
Using the inequalities (3.12) and the definition (3.17) of γ ih ,

|λ ih -λ i | ≤ λ ih -∇u ih 2 + | ∇u ih 2 -λ i | ≤ λ ih -∇u ih 2 + γ ih , so that (3.18b) is proven.
Finally, (3.18a) can be seen as in, e.g., Carstensen and Gedicke [15, Lemma 3.1] combined with u i -

u ih ≤ α ih [-].

Dual norm of the residual equivalences

We now estimate the dual residual norm Res(u ih , λ ih ) -1 , for u ih ∈ V a piecewise polynomial of degree p ≥ 1 and λ ih ∈ R. For the upper bound, following [START_REF] Prager | Approximations in elasticity based on the concept of function space[END_REF][START_REF] Destuynder | Explicit error bounds in a conforming finite element method[END_REF][START_REF] Braess | Equilibrated residual error estimates are p-robust[END_REF][START_REF]Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations[END_REF] and [21, [START_REF] Papež | Estimating and localizing the algebraic and total numerical errors using flux reconstructions[END_REF][START_REF] Papež | Sharp algebraic and total a posteriori error bounds via a multilevel approach[END_REF] for inexact solvers, see also the references therein, we introduce an equilibrated flux reconstruction. This is a vector field σ ih constructed from the local residual of (u ih , λ ih ) by solving patchwise mixed finite element problems such that

σ ih ∈ V h ⊂ H(div, Ω), (4.1a)
∇•σ ih = λ ih u ih -ρ ih , λ -1 2 1
ρ ih sufficiently small.

Meshes and discrete spaces

We first introduce some more notation. Let henceforth {T h } h be a family of matching simplicial partitions of the domain Ω, shape regular in the sense that the ratio of each element diameter to the diameter of its largest inscribed ball is uniformly bounded by a constant κ T > 0. We denote by K a generic element of T h . The set of vertices is denoted by V h , with interior vertices V int h , vertices located on the boundary V ext h , and a generic vertex a. We call T a the patch of elements of T h which share the vertex a ∈ V h , ω a the corresponding subdomain, and n ωa its outward unit normal. We often tacitly extend functions defined on ω a by zero outside of ω a , whereas V h (ω a ) stands for the restriction of the space V h to ω a . Next, ψ a for a ∈ V h stands for the piecewise affine "hat" function taking value 1 at the vertex a and zero at the other vertices. Remarkably, (ψ a ) a∈V h form a partition of unity via a∈V h ψ a = 1| Ω .

Let P s (K), s ≥ 0, stand for polynomials of total degree at most s on K ∈ T h , and P s (T h ) for piecewise polynomials on T h , without any continuity requirement. Let also Brezzi and Fortin [10] or Roberts and Thomas [START_REF] Roberts | Mixed and hybrid methods[END_REF]. We also denote by Π Q h the L 2 (Ω)-orthogonal projection onto Q h .

V h × Q h ⊂ H(div, Ω) × L 2 (Ω) stand for the Raviart-Thomas-Nédélec (RTN) mixed finite element spaces of for degree p + 1, i.e., V h := {v h ∈ H(div, Ω); v h | K ∈ [P p+1 (K)] d + P p+1 (K)x} and Q h := P p+1 (T h ), see

Equilibrated flux reconstruction for inexact solvers

Let r ih ∈ P p (T h ) be a discontinuous piecewise p-degree polynomial that lifts the misfit in the Galerkin orthogonality of the residual Res(u ih , λ ih ), i.e.

Res(u

ih , λ ih ), v h V ,V = λ ih (u ih , v h ) -(∇u ih , ∇v h ) = (r ih , v h ) ∀v h ∈ P p (T h ) ∩ V. (4.2)
A simple elementwise construction of r ih is proposed in [47, equation (5.2)]. Typically, r ih = 0 for an "exact" discrete algebraic solve that we do not suppose here.

We construct σ ih in two steps. First, solve the following homogeneous local Neumann (Neumann-Dirichlet close to the boundary) discrete problems on patches ω a :

Definition 4.1 (Equilibrated flux reconstruction). [-] For a ∈ V h , set V a h := {v h ∈ V h (ω a ); v h •n ωa = 0 on ∂ω a }, Q a h := {q h ∈ Q h (ω a ); (q h , 1) ωa = 0}, a ∈ V int h , V a h := {v h ∈ V h (ω a ); v h •n ωa = 0 on ∂ω a \ ∂Ω}, Q a h := Q h (ω a ), a ∈ V ext h . Then define σ ih,dis := a∈V h σ a ih,dis ∈ V h , where σ a ih,dis ∈ V a h solve σ a ih,dis := arg min v h ∈V a h , ∇•v h =Π Q h (λ ih u ih ψa-∇u ih •∇ψa-r ih ψa) ψ a ∇u ih + v h ωa ∀a ∈ V h . (4.3)
Note that the Euler-Lagrange equations for (4.3) give the standard mixed finite element formulation, cf. [START_REF]Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations[END_REF]Remark 3.7

]: find σ a ih,dis ∈ V a h and p a h ∈ Q a h such that (σ a ih,dis , v h ) ωa -(p a h , ∇•v h ) ωa = -(ψ a ∇u ih , v h ) ωa ∀v h ∈ V a h , (4.4a)
(∇•σ a ih,dis , q h ) ωa = (λ ih u ih ψ a -∇u ih •∇ψ a -r ih ψ a , q h ) ωa ∀q h ∈ Q a h . (4.4b)
Consequently, ∇•σ ih,dis = λ ih u ih -r ih , cf., e.g., [START_REF] Papež | Sharp algebraic and total a posteriori error bounds via a multilevel approach[END_REF]Lemma 3.6]. Now, proceeding as in [46, Section 3.2], one can construct in a multilevel way a second flux reconstruction σ ih,alg ∈ V h such that ∇•σ ih,alg = r ih . Consequently, setting σ ih := σ ih,dis + σ ih,alg , (4.1b) follows with ρ ih = 0. Other strategies are pursued in [21,[START_REF] Papež | Estimating and localizing the algebraic and total numerical errors using flux reconstructions[END_REF], where the algebraic residual is included differently into (4.3)/(4.4b). These approaches yield

∇•σ ih,alg = r ih -ρ ih (4.5)
with ρ ih = 0 and are based on precomputing some algebraic solver iterations in order to ensure that ρ ih is sufficiently small with respect to the two other contributions in (4.9a) below, more precisely verifying (4.9b).

Conforming local residual liftings

To estimate Res(u ih , λ ih ) -1 from below, we solve conforming primal counterparts of problems (4.4), without the term with r ih . On each patch ω a around the vertex a ∈ V h , define

H 1 * (ω a ) := {v ∈ H 1 (ω a ); (v, 1) ωa = 0}, a ∈ V int h , (4.6a)
H 1 * (ω a ) := {v ∈ H 1 (ω a ); v = 0 on ∂ω a ∩ ∂Ω}, a ∈ V ext h , (4.6b)
and let X a h be an arbitrary discrete subspace of

H 1 * (ω a ), typically P p+1 (T a ) ∩ H 1 * (ω a ). Definition 4.2 (Conforming local Neumann problems). Define r a ih ∈ X a h by (∇r a ih , ∇v h ) ωa = Res(u ih , λ ih ), ψ a v h V ,V ∀v h ∈ X a h
for each a ∈ V h . Then set r ih := a∈V h ψ a r a ih . The functions r a ih are discrete Riesz projections of the local residual with hat-weighted test functions. As all ψ a r a ih ∈ H 1 0 (ω a ), r ih ∈ V , though r a ih ∈ V .

Dual norm of the residual equivalences

Following Carstensen and Funken [13, Theorem

ψ a ∇u ih + σ a ih,dis ωa ≤ C st sup v∈H 1 * (ωa); ∇v ωa =1 { Res(u ih , λ ih ), ψ a v V ,V -(r ih , ψ a v) ωa }. (4.8)
The constant C st > 0 again only depends on κ T , and a computable upper bound on C st is given in [START_REF]Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations[END_REF]Lemma 3.23]. We can summarize the main result of this section:

Theorem 4.3 (Residual equivalences). Let (u ih , λ ih ) ∈ P p (T h ) ∩ V × R
Res(u ih , λ ih ) -1 ≤ ∇u ih + σ ih,dis + σ ih,alg + λ -1 2 1 ρ ih , (4.9a
)

∇u ih + σ ih,dis + σ ih,alg + λ -1 2 1 ρ ih ≤ 3(d + 1)C st C cont,PF Res(u ih , λ ih ) -1
when σ ih,alg +λ

-1 2 1 ρ ih ≤ (2(d + 1)C st C cont,PF ) -1 ∇u ih + σ ih,dis , (4.9b
)

Res(u ih , λ ih ), r ih V ,V ∇r ih ≤ Res(u ih , λ ih ) -1 . (4.9c) Proof. Fix v ∈ V with ∇v = 1.
Using definition (2.4a), adding and subtracting (σ ih , ∇v), and employing the Green theorem and the equilibrium (4.1b) yields

Res(u ih , λ ih ), v V ,V = λ ih (u ih , v) -(∇u ih , ∇v) = (ρ ih , v) -(∇u ih + σ ih , ∇v).
Thus, definition (2.4b) of the dual norm of the residual and the Cauchy-Schwarz, Poincaré-Friedrichs, and triangle inequalities yield the bound (4.9a). This actually also holds for V h being the cheaper RTN space of order p and not p + 1, as (4.1b) still holds. To prove (4.9b), we proceed as in [

Guaranteed and fully computable upper and lower bounds

We combine here the different results of the previous sections to derive the actual guaranteed and fully computable bounds for eigenvalues (in Section 5.1) and eigenvectors (in Section 5.2). A discussion of the results is provided in Section 5.3. We will sometimes use ζ (ih) ∈ V , the solution of the Laplace source problem -∆ζ

(ih) = r (ih) in Ω, ζ (ih) = 0 on ∂Ω, i.e., (∇ζ (ih) , ∇v) = (r (ih) , v) ∀v ∈ V. (5.1)
We also denote by V h := P 1 (T h ) ∩ V the lowest-order conforming finite element space, i.e., the span of ψ a over all a ∈ V int h , and by h the maximal diameter of all K ∈ T h .

Eigenvalues

We first tackle the upper and lower bounds for the i-th eigenvalue λ i . We discuss the necessary auxiliary bounds below in Remark 5.3.

Theorem 5.1 (Guaranteed lower bounds for the i-th eigenvalue). Let the i-th eigenvalue, i ≥ 1, be simple and suppose the auxiliary bounds

λ 1 ≤ λ 1 , λ i ≤ λ i , λ i+1 ≤ λ i+1 , as well as λ i-1 ≤ λ i-1 when i > 1, for λ 1 , λ i , λ i+1 , λ i-1 > 0. Let (u ih , λ ih ) be any element of P p (T h ) ∩ V × R + verifying [-] u ih = 1 and the inequalities λ i-1 < λ ih when i > 1, λ ih < λ i+1 . (5.2)
Let next σ ih,dis and r ih be respectively constructed following Definitions 4.1 and 4.2, let σ ih,alg ∈ V h verify (4.5) for an inexact solver, and define

η i,res := ∇u ih + σ ih,dis + σ ih,alg + λ -1 2 1 ρ ih .
Set

c ih := max λ ih λ i-1 -1 -1 , 1 - λ ih λ i+1 -1 , (5.3a
)

cih := max λ i-1 -1 2 λ ih λ i-1 -1 -1 , λ -1 2 i+1 1 - λ ih λ i+1 -1 , (5.3b)
with the left terms in the max disregarded for i = 1. Then

∇u ih 2 -η 2 i ≤ λ i , (5.4)
where we distinguish the following three cases:

Case A (No smallness assumption) If (u i , u ih ) ≥ 0 is known to hold, define α ih := √ 2c ih η i,res ; if only (u ih , χ i ) > 0 holds, set α ih := √ 2(1-u ih -Π i u ih ) -1 2 cih η i,res
, where Π i u ih stands for the L 2 (Ω)-orthogonal projection of u ih on the span of χ i . Then (5.4) holds with

η 2 i := η 2 i,res + λ ih + λ i α 2 ih . (5.5)
Case B (Improved estimates under a smallness assumption) Let (u ih , χ i ) > 0, define α ih := √ 2c ih η i,res and request

α ih ≤ min 2λ 1 λ i 1 2 , χ i -1 (u ih , χ i ) . (5.6)
Then, (5.4) holds with

η 2 i := c 2 ih 1 - λ i λ 1 α 2 ih 4 -1 η 2 i,res . (5.7)
Case C (Optimal estimates under elliptic regularity assumption) Let (u ih , χ i ) > 0 and assume that the solution ζ (ih) of problem (5.1) belongs to the space H 1+δ (Ω), 0 < δ ≤ 1, so that the approximation and stability estimates

min v h ∈V h ∇(ζ (ih) -v h ) ≤ C I h δ |ζ (ih) | H 1+δ (Ω) , (5.8a
)

|ζ (ih) | H 1+δ (Ω) ≤ C S r (ih) (5.8b) are satisfied. Define α ih := √ 2c ih C I C S h δ η i,res +λ -1 2 1 σ ih,alg + λ -1 2 1
ρ ih and let

α ih ≤ χ i -1 (u ih , χ i ).
(5.9)

Then (5.4) holds with η 2 i given by (5.5).

Theorem 5.2 (Improved guaranteed upper bounds for the i-th eigenvalue). Let the assumptions of Theorem 5.1 be satisfied, with the auxiliary bounds

λ 1 ≤ λ 1 , λ i ≤ λ i ≤ λ i , for λ 1 , λ i , λ i > 0. Let also λ i ≤ ∇u ih 2 , see Remark 5.4. Set c ih := 1 if i = 1, c ih := max λ ih λ 1 -1 2 , 1 if i > 1, (5.10a
)

d ih := λ 2 i c 2 ih + 4λ i Res(u ih , λ ih ), r ih 2 V ,V ∇r ih 2 + 4λ i c ih λ ih -∇u ih 2 .
(5.10b)

Then λ i ≤ ∇u ih 2 -η2 i , (5.11)
with, in Cases A and C,

η2 i := max -λ i α 2 ih + 1 2 d ih -λ i c ih + 2 λ ih -∇u ih 2 , 0 , (5.12)
and, in Case B, for i = 1 only,

η2 1 := max 1 4 1 - ∇u 1h 2 λ 2 1 - α 2 1h 4 d 1h -λ 1 + 2 λ 1h -∇u 1h 2
, 0 .

(5.13) Remark 5.3 (Auxiliary bounds λ 1 , λ i , and λ i+1 ). A straightforward consequence of the min-max principle for self-adjoint operators, see, e.g., Gilbarg and Trudinger [START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF], is that

Ω ⊆ Ω + ⇒ λ k (Ω + ) ≤ λ k and Ω -⊆ Ω ⇒ λ k ≤ λ k (Ω -) for all k ≥ 1, where λ k (Ω ±
) is the k-th eigenvalue on Ω ± . We can then obtain all λ 1 , λ i , and λ i+1 necessary in Theorem 5.1 by this domain inclusion for Ω + with known exact eigenvalues (typically rectangular d-parallelepipeds or d-spheres, cf. [START_REF] Trefethen | Computed eigenmodes of planar regions[END_REF]). In what concerns λ i , a very precise choice is to use

λ i := ∇u ih 2 -η 2 i
, where η 2 i was first computed with a rather rough bound λ i . For λ i+1 , if the analytic bounds are too rough to be useful, guaranteed and easily computable numerical bounds can be used from Liu and Oishi [START_REF] Liu | Verified eigenvalue evaluation for the Laplacian over polygonal domains of arbitrary shape[END_REF] (on convex domains for d = 2), Carstensen and Gedicke [16], or Liu [START_REF] Liu | A framework of verified eigenvalue bounds for self-adjoint differential operators[END_REF], typically on a quite coarse mesh. Finally, as a "practical gratis" strategy for λ i+1 , one may simply use λ (i+1)h computed by the linear algebra toolbox when solving for (λ ih , u ih ), see, e.g., Saad [START_REF] Saad | Numerical methods for large eigenvalue problems[END_REF] and the references therein. Then Theorem 5.1 (and 5.6 below) may no longer hold stricto sensu, but sharp upper and lower bounds are still expected/observed in practice.

Remark 5.4 (Auxiliary bounds λ i-1 and λ i ). When (u ih , λ ih ) is given by the conforming finite element method of Section 6 below, with an exact solver leading to satisfaction of (6.1) with (u ih , u jh ) = δ ij , there holds λ i ≤ λ ih = ∇u ih 2 and similarly λ i-1 ≤ λ (i-1)h = ∇u (i-1)h 2 , leading to rather precise auxiliary bounds λ i and λ i-1 . For the first eigenvalue, there holds λ 1 ≤ ∇u 1h 2 for an arbitrary u 1h ∈ H 1 0 (Ω). For the i-th eigenvalue, i > 1, we in general need to resort to the min-max principle giving

λ i ≤ max ξ∈R i , ξ =1 ∇ i k=1 ξ k u kh 2 i k=1 ξ k u kh 2
for an arbitrary linearly independent i-tuple (u 1h , . . . , u ih ), where ξ 2 = i k=1 ξ 2 k .

Remark 5.5 (Constants C I and C S ). Let Ω be a convex polygon in R 2 . Then it is classical that the solution ζ (ih) of (5.1) belongs to H 2 (Ω) and |ζ (ih) | H 2 (Ω) = ∆ζ (ih) = r (ih) , so that δ = 1 and C S = 1, see Grisvard [START_REF] Grisvard | Elliptic problems in nonsmooth domains[END_REF]Theorem 4.3.1.4]. In this situation, calculable bounds on C I can be found in Liu and Kikuchi [START_REF] Liu | Analysis and estimation of error constants for P 0 and P 1 interpolations over triangular finite elements[END_REF] and Carstensen et al. [17], see also Liu and Oishi [41, Section 2] and the references therein; in particular, for a mesh formed by isosceles right-angled triangles, C I ≤ 0.493 √ 2 . We now prove Theorems 5.1 and 5.2, separately for each case: Proof (Case A). 1) Lower bound of Theorem 5.1. If (u i , u ih ) ≥ 0 is known to hold, we can start from the L 2 (Ω) bound (3.10). If this is not the case but (u ih , χ i ) > 0 holds, we first inspect the proof of Lemma 3.2 to obtain an alternative L 2 (Ω) estimate. We have -2(u i , u ih ) = -2(u i , u ih -Π i u ih ) -2(u i , Π i u ih ). Note that the second term is negative by the sign assumption (u i , χ i ) > 0 on u i . So, instead of (3.8), as

u i = 1 and u ih -Π i u ih < 1, u i -u ih 2 ≤ 2 + 2 u ih -Π i u ih =: δ ih < 4.
Consequently, the quadratic inequality in the proof of Lemma 3.2 implies u i -u ih 2 ≤ ∇r (ih) 2 C -1 ih (1δ ih /4) -1 . Thus, the bound (4.9a) and assumption (5.2) enable us to give a computable upper bound on the L 2 (Ω) error by the estimator α ih ; note that min{a, b} -1 2 = max{a -1 2 , b -1 2 }, linking the constant C ih of (3.9) with cih of (5.3b). Consequently, the bound in (5.4) follows by combining the upper bounds in (3.12), (3.18a), and once again (4.9a).

2) Upper bound of Theorem 5.2. We start from the lower bound in (3.12). We then need to bound ∇(u i -u ih ) 2 from below, for which we use (3.18b). Relying on the simplifying assumption λ i ≤ ∇u ih 2 , satisfied namely in cases discussed in Remark 5.4, γ ih of (3.17) simplifies to ∇(u i -u ih ) 2 . Thus (3.18b) forms a quadratic inequality for ∇(u i -u ih ) 2 , yielding, in combination with (4.9c),

∇(u i -u ih ) 2 ≥ 1 2 d ih -λ i c ih + 2 λ ih -∇u ih 2 .
(5.14) Thus (5.11) with the estimator (5.12) follows.

Proof (Case B). The proof proceeds as above. Note that conditions in (5.6) imply that α ih ≤ 2 λ1 λi and α ih ≤ χ i -1 (u ih , χ i ) for α ih of (3.10). We can thus use Lemma 3.3 to find that (u i , u ih ) is indeed nonnegative, Lemma 3.2 for the L 2 (Ω) bound, and the improved estimates (3.19) of Theorem 3.5 and (3.13) of Theorem 3.4. For the latter, that seems only to hold for the first eigenpair, we also employ the inequality

r (ih) 2 = (∇ζ (ih) , ∇r (ih) ) = (∇(ζ (ih) -ζ ih ), ∇r (ih) ) + (r ih , ζ ih ),
where ζ ih ∈ V h is the minimizer in (5.8a). Employing (4.5), the Green theorem, the Poincaré-Friedrichs

inequality ζ ih ≤ λ -1 2 1
∇ζ ih , and the stability ∇ζ ih ≤ ∇ζ (ih)

(r ih , ζ ih ) = -(σ ih,alg , ∇ζ ih ) + (ρ ih , ζ ih ) ≤ σ ih,alg + λ -1 2 1 ρ ih ∇ζ (ih) .
Noting that (5.1) gives ∇ζ (ih) ≤ λ

Eigenvectors

We now summarize our estimate on the energy error in the approximation of the i-th eigenvector, as well as its efficiency and robustness:

Theorem 5.6 (Guaranteed and robust bound for the i-th eigenvector error). Let the assumptions of Theorem 5.1 be verified. Then the energy error can be bounded via

∇(u i -u ih ) ≤ η i , (5.15)
where η i is defined in the Cases A and C by (5.5) and in Case B by (5.7), with appropriate α ih . Under condition (4.9b), all these estimators η i are efficient as

η 2 i,res ≤ 3 2 (d + 1) 2 C 2 st C 2 cont,PF λ ih -∇u ih 2 + γ ih 2 λ i + C ih ∇(u i -u ih ) 2 .
(5.16)

Proof. The guaranteed error bound (5.15) follows as in Theorem 5.1 upon combining the upper bounds in estimates (3.18) or (3.19) together with (4.9a). The efficiency (5.16) is a consequence of (4.9b) and of (3.18b).

Comments

We collect here comments about Theorems 5.1, 5.2, and 5.6.

Remark 5.7 (Inexact solvers). The results of Theorems 5.1, 5.2, and 5.6 are presented in a general context of an inexact algebraic solver. Shall the algebraic residual representer r ih in (4.2) be zero and the Rayleigh quotient link ∇u ih 2 = λ ih hold, many simplifications apply. Namely, σ ih,alg = 0, ρ ih = 0, and ∇u ih 2 can be systematically replaced by λ ih , making in particular disappear numerous terms in the estimators η i and ηi . The polynomial-degree-robust efficiency (5.16) holds under the condition (4.9b) only, which is a typical inexactness (stopping) criterion. For the elliptic regularity Case C, though, it appears wise to only stop the iterations when λ

-1 2 1 σ ih,alg + λ -1 2
1
ρ ih is comparable to the first term in α ih .

Remark 5.8 (Sharpness and comparison of the different bounds of Theorems 5.1 and 5.6). The advantage of Case A is that it holds on an arbitrarily coarse mesh, provided that only the structural assumption (5.2) holds. It may, however, lead to a larger overestimation of the error. Case B, under the "fine enough mesh" condition (5.6), then significantly improves the multiplicative factor in front of the central term

η i,res = ∇u ih + σ ih,dis + σ ih,alg + λ -1 2 1
ρ ih , in limit to the factor c ih given by (5.3a). The bound of Case B still holds without any regularity/convexity/dimension assumption and all the quantities appearing are known. Finally, also the factor c ih is asymptotically removed in Case C, when δ > 0 and h → 0. Here, however, elliptic regularity is needed, see Remark 5.5.

Remark 5.9 (Dependence on the maximal element diameter h). The maximal element diameter h is not present at all in Cases A and B of Theorem 5.1 and it does not necessarily need to tend to zero in Case C: it only appears as a multiplicative factor of the principal estimator η i,res . [-] This stands in contrast to previous guaranteed results like [41,

(u i -u ih ) 2 C ih + ∇(ui-u ih ) 2 λi
for an exact algebraic solver in the context of the finite element method (6.1) below. Noting that ∇(ui-u ih ) 2 λi ≤ 2(λi+λ ih ) λi (in fact this term becomes negligible with mesh refinement/increasing the polynomial degree), we conclude that the result of Theorem 5.6 is fully robust with respect to the polynomial degree p of u ih : all the constants in the comparison between the error ∇(u i -u ih ) and the estimate featuring ∇u ih + σ ih,dis are independent of p. Note, though, that the factor C ih given by (3.16) deteriorates for higher eigenvalues.

Remark 5.11 (Error localization and mesh adaptivity). Since there holds η 2 i,res ≤ 3

K∈T h ∇u ih + σ ih,dis 2 K + σ ih,alg 2 K + λ -1 1 ρ ih 2 K
, a crucial advantage of the estimators of Theorems 5.1 and 5.6 is that these local contributions can be used in adaptive mesh refinement based on marking strategies. This is tightly linked to Remark 5.9.

Application to conforming finite elements

We verify in this section the conditions of the application of our results to the conforming finite element method.

Let V h := P p (T h ) ∩ V for a given polynomial degree p ≥ 1. In the finite element method, the exact i-th

eigenpair (u ih , λ ih ) ∈ V h × R + is such that (u ih , u jh ) = δ ij , 1 ≤ i, j ≤ dim V h , and (∇u ih , ∇v h ) = λ ih (u ih , v h ) ∀v h ∈ V h , (6.1)
with the signs ideally fixed by (u i , u ih ) ≥ 0, practically by (u ih , χ i ) > 0. Thus, upon verifying (5.2) and possibly checking (5.6) or (5.9), all the results of Theorems 5.1, 5.2, and 5.6 hold for any p ≥ 1. Note that an inexact solution of (6.1) in the form (4.2) is taken into account. Shall (6.1) hold, r ih in (4.2) vanishes and, moreover, choosing v h = u ih in (6.1) yields

∇u ih 2 = λ ih . N h ndof λ 2 -λ 1h (5.2) χ 1 -1 (u 1h , χ 1 ) -α 1h (5.9) λ 1 =1.5π 2 λ 2 =4.5π

Numerical experiments

We finally numerically illustrate the estimates of Theorems 5.1, 5.2, and 5.6 on three test cases in R 2 , for conforming finite elements (6.1) of order p = 1. We actually only use the cheaper Raviart-Thomas-Nédélec space of degree p = 1 for the flux equilibration instead of p + 1. This still gives guaranteed bounds, see the proof of Theorem 4.3, and we do not observe any asymptotic loss of efficiency. The implementation was done in the FreeFem++ code [START_REF] Hecht | New development in FreeFem++[END_REF]. When we only consider one eigenvalue, it is implicitly assumed that we have chosen χ 1 = 1 for the sign characterization. We consider five test settings with an exact solver and illustrate the use of an inexact solver in a sixth one.

First eigenvalue on the unit square

We start by testing the framework on a unit square Ω = (0, 1) 2 and focus on the first eigenvalue. The eigenvalues on a square of size H being π 2 (k 2 + l 2 )/H 2 , k, l = 1, . . . , ∞, the first and second eigenvalues are λ 1 = 2π 2 and λ 2 = 5π 2 , respectively. In consequence, we can easily choose different λ 1 ≤ λ 1 and λ 2 ≤ λ 2 for the auxiliary eigenvalue bounds and analyze the sensitivity of our results with respect to these choices. The first eigenfunction is given by u 1 (x, y) = sin(πx) sin(πy). We focus here on the refined elliptic regularity of Case C, since d = 2 and the domain is convex, with constants C S = 1 and δ = 1 given in Remark 5.5.

Structured mesh

We first illustrate in Table 1 how quickly the computable conditions (5.2) and (5.9) are satisfied under a uniform refinement of a structured mesh. We take C I = 0.493 √ 2 following Remark 5.5 and consider N = 3, 4, 5 subdivisions of each boundary of Ω for the two choices λ 1 = 1.5π 2 , λ 2 = 4.5π 2 and λ 1 = 0.5π 2 , λ 2 = 3π 2 , respectively. Note that the finite element space on the coarsest mesh such that all conditions are satisfied contains 25, respectively 36, degrees of freedom only. Indeed, it turns out that our conditions are rather mild.

Next, Figure 1 (left) illustrates the convergence of the error λ 1h -λ 1 as well as of its lower and upper bounds η2

1 , η 2 1 given by Case C of Theorems 5.1 and 5.2. We also plot the energy error in the eigenfunction ∇(u 1 -u 1h ) and its upper bound η 1 of Theorem 5.6, Case C. The convergence rates are optimal as expected from the theory.

We present in Table 2 precise numbers of the lower and upper bounds λ 1h -η 2 1 ≤ λ 1 ≤ λ 1h -η2 1 on the exact eigenvalue λ 1 , the effectivity indices of the lower and upper bounds η2

1 ≤ λ 1h -λ 1 ≤ η 2 1 of the error λ 1h -λ 1 , and the effectivity index of the upper bound ∇(u 1 -u 1h ) ≤ η 1 , given respectively by

I lb λ,eff := λ 1h -λ 1 η2 1 , I ub λ,eff := η 2 1 λ 1h -λ 1 , I ub u,eff := η 1 ∇(u 1 -u 1h ) . (7.1)
We observe rather sharp results, and this also for the relative size of the first eigenvalue confidence interval

E λ,rel := 2 (λ 1h -η2 1 ) -(λ 1h -η 2 1 ) (λ 1h -η2 1 ) + (λ 1h -η 2 1 ) . ( 7
h η 2 1 λ 1h -λ1 η2 1 η1 ∇(u1 -u 1h )
Figure 1: [Unit square] Error in the eigenvalue and eigenvector approximation, its lower bound (eigenvalue only), and its upper bound for the choice λ 1 = 1.5π 2 , λ 2 = 4.5π 2 ; sequence of structured (left) and unstructured but quasi-uniform (right) meshes; Case C Table 2: [Unit square, structured mesh] Lower and upper bounds on the exact eigenvalue λ 1 , the effectivity indices, and size of the relative λ 1 confidence interval; λ 1 = 1.5π 2 , λ 2 = 4.5π 2 ; Case C

N h ndof λ 1 λ 1h λ 1h -η 2 1 λ 1h -η2 1 I lb λ,eff I ub λ,eff E λ,rel I ub u,eff 10

Unstructured mesh

Consider now a sequence of unstructured quasi-uniform meshes, obtained by an initial partition of each boundary edge into N intervals. Conditions (5.2) and (5.9) turn here to be satisfied similarly as in Table 1.

The convergence plots for this case are presented in Figure 1 (right), showing a similar behavior as for the structured meshes. This time, we use the upper bound on C I according to [START_REF] Liu | Analysis and estimation of error constants for P 0 and P 1 interpolations over triangular finite elements[END_REF]Eqn. (46)]:

C I = 0.493 max K∈T h 1+| cos(θ K )| sin(θ K ) ν+(α K ,θ K ) 2 h [40] K h K .
We refer to [START_REF] Liu | Analysis and estimation of error constants for P 0 and P 1 interpolations over triangular finite elements[END_REF] for the definition of h [START_REF] Liu | Analysis and estimation of error constants for P 0 and P 1 interpolations over triangular finite elements[END_REF] K and other notation. We observe in Table 3 that the results are similar to structured meshes; in particular the case of λ 1 = 0.5π 2 , λ 2 = 3π 2 is less sensitive to the unstructured mesh (not presented).

First eigenvalue on an L-shaped domain: mesh adaptivity

We next consider the L-shaped domain Ω := (-1, 1) 2 \ [0, 1] × [-1, 0], where λ 1 ≈ 9.6397238440 is known to high accuracy [START_REF] Trefethen | Computed eigenmodes of planar regions[END_REF]. Including Ω into the square Ω + = (-1, 1) 2 , cf. Remark 5.3, we take λ 1 = λ 1 (Ω + ) = π 2 /2, whereas λ 2 = 15.1753 from Table 1 of [START_REF] Liu | A framework of verified eigenvalue bounds for self-adjoint differential operators[END_REF] is employed. We test here the Cases A and B within an adaptive refinement strategy. To do so, we use the local character of our estimators, see Remark 5.11. We employ the Dörfler marking with θ = 0.6 and the newest vertex bisection mesh refinement.

Table 4 illustrates whether the conditions (5.2) and (5.6) are satisfied under this adaptive refinement. Figure 2 (right) illustrates the error in the eigenvalue and the eigenvector and their bounds (5.4), (5.11), and (5.15). Optimal convergence rates are indicated by dashed lines. The initial mesh is structured with 22 degrees of freedom and the conditions (5.2) and (5.6) are all satisfied starting from 140 degrees of freedom. The transition from Case A to Case B in Theorems 5.1, 5.2, and 5.6 is marked by a dotted line. Figure 2 (left) then depicts an adaptively refined mesh and We next consider a domain with a polygonal hole, see Figure 3 (left) illustrating the mesh used at iteration 20 of our adaptive mesh refinement strategy. The lower bounds λ 1 and λ 2 on the first and second eigenvalue have been obtained once and for all before starting the adaptive algorithm following the estimates derived in [START_REF] Liu | A framework of verified eigenvalue bounds for self-adjoint differential operators[END_REF], on a uniform mesh with 1143 nodes. Figure 3 (right) shows the interval between our lower (λ 1h -η 2 1 ) and upper (λ 1h -η2

1 ) bounds on the first eigenvalue, relying on Case B of Theorems 5.1 and 5.2, whose assumptions hold starting from 2494 degrees of freedom; Table 6 states the numbers. Note that the interval size (λ 1h -η2 1 ) -(λ 1h -η 2 1 ) = η 2 1 -η2 1 behaves like 1/ndof.

Higher eigenvalues

We now test the upper and lower bounds for higher eigenvalues. First we consider the unit triangle with vertices (0, 0), (1, 0), (0, 1) and a family of structured meshes. The auxiliary lower bounds are obtained by a computation on a fixed coarse mesh with 2145 triangles following [START_REF] Liu | A framework of verified eigenvalue bounds for self-adjoint differential operators[END_REF], which results in Figure 4 gives the convergence plots for the first four eigenvalues and Table 7 provides more details on absolute numbers and efficiency. As the domain is convex (case C), we obtain excellent upper bounds for the error in all four eigenvalue/eigenvector pairs. The lower bound of the eigenvalue error (the improved eigenvalue upper bound of Theorem 5.2) is, however, degrading for higher eigenvalues. We now apply the same setting to the L-shaped domain where we obtain again the auxiliary lower bounds by the method presented in [START_REF] Liu | A framework of verified eigenvalue bounds for self-adjoint differential operators[END_REF] for a coarse structured mesh with 3201 triangles resulting in 8 provides more details on the efficiency. We now observe that the efficiency also degrades for the upper bound of the eigenvalue and eigenvector error. Further, improved lower bounds of the eigenvalue error are not available for the considered meshes for i > 1. This appears as the resulting ηi are all equal to zero, see (5.12), respectively (5.13), so that our eigenvalue upper bound stays that of the finite element method. For all meshes and all considered eigenvalues, though, our estimates still give a rather tight guaranteed eigenvalue confidence interval and quite reasonable eigenvector effectivity indices. We can also observe by a jump of the blue curve (η 2 i ) the change between the cases A and B. The critical mesh size where this change occurs seems to degrade with increasing eigenvalues.

Inexact algebraic eigenvalue solvers

We finally consider inexact eigenvalue solvers. Since we are using FreeFem++, we rely on an algebraic eigenvalue solver based on the ARPACK package that is built in FreeFem++. Here a user-specified tolerance can be provided and we choose it in a mesh-dependent way as tol(h) = h 2 to materialize an inexact solver. We set σ ih,dis following Definition 4.1. In order to compute σ ih,alg in (4.5), we proceed as in [START_REF] Papež | Estimating and localizing the algebraic and total numerical errors using flux reconstructions[END_REF] and the references therein and first compute a second reconstructed flux σih,dis corresponding to some additional algebraic iterations (here corresponding to the tolerance h 2 /100 in ARPACK); then σ ih,alg := σih,dis -σ ih,dis . Figure 6 demonstrates that we still obtain excellent lower and upper bounds. Adaptive stopping criteria of the form (4.9b), leading to savings in algebraic solver iterations, are not investigated here.

Comparison with existing results

We finally compare our results with some existing ones from [16,[START_REF] Liu | Verified eigenvalue evaluation for the Laplacian over polygonal domains of arbitrary shape[END_REF][START_REF] Liu | A framework of verified eigenvalue bounds for self-adjoint differential operators[END_REF]. In what concerns the unit square and the first eigenvalue of Section 7.1, our estimates appear sharper while comparing Table 2 with the estimates presented in [16], see Figure 6.2 therein. For the L-shaped domain and uniformly refined meshes of Section 7.4 for the first eigenvalue, we also obtain better results than those presented in [16, Figure 6.4], where an efficiency issue appears; compared to the results presented in [41, Table 5.5], we observe that our lower bound λ 1h -η 2 1 of the exact eigenvalue is a little less sharp, whereas the upper bound λ 1h -η2 1 is not present in [START_REF] Liu | Verified eigenvalue evaluation for the Laplacian over polygonal domains of arbitrary shape[END_REF]. Recall also from Section 1 that our estimates are much cheaper here than those of [START_REF] Liu | Verified eigenvalue evaluation for the Laplacian over polygonal domains of arbitrary shape[END_REF] (there is no auxiliary eigenvalue problem to solve). For adaptive meshes, we observe that our efficiency of the confidence interval for the first eigenvalue as measured in [16] by or [START_REF] Liu | A framework of verified eigenvalue bounds for self-adjoint differential operators[END_REF]. For the eigenvalue upper bounds in the CR case, we evaluate the Rayleigh quotient on the P 1 conforming nodal averaging of the original eigenvectors. On the convex triangle, the present approach seems to give the sharpest results, whereas on the L-shaped domain, the method based on the CR finite elements with the constant from [START_REF] Liu | A framework of verified eigenvalue bounds for self-adjoint differential operators[END_REF] is better for the lower bound. Recall, though, that important advantages of the present theory are that it additionally gives a guaranteed control of the eigenvector error by the same estimators, is not specific to a particular scheme but yields general results that are here applied to any order conforming finite element method and extended in [START_REF]Guaranteed and robust a posteriori bounds for Laplace eigenvalues and eigenvectors: a unified framework[END_REF] to basically any numerical scheme, and achieves polynomial-degree robustness. It can also be noted that the present estimators take elementwise form immediately suitable for adaptive mesh refinement. Table 8: [L-shaped domain, uniform mesh refinement] Lower and upper bounds on the first four exact eigenvalues λ i , the effectivity indices, and the sizes of the relative λ i confidence intervals

h ndof λ i λ ih λ ih -η 2 i λ ih -η2 i I lb λ,eff I ub λ,eff E λ,rel I ub u,eff 40
[13] C. Carstensen and S. A. Funken, Fully reliable localized error control in the FEM, SIAM J. Sci. Comput., 21 (1999/00), pp. 1465-1484.

[14] C. Carstensen and D. Gallistl, Guaranteed lower eigenvalue bounds for the biharmonic equation, Numer. Math., 126 (2014), pp. 33-51.

[15] C. Carstensen and J. Gedicke, An oscillation-free adaptive FEM for symmetric eigenvalue problems, Numer. Math., 118 (2011), pp. 401-427.

[16] , Guaranteed lower bounds for eigenvalues, Math. Comp., 83 (2014), pp. 2605-2629.

[17] C. Carstensen, J. Gedicke, and D. Rim, Explicit error estimates for Courant, Crouzeix-Raviart and Raviart-Thomas finite element methods, J. Comput. Math., 30 (2012), pp. 337-353.

λ k 2 =

 2 C ih . Thus, similarly to (3.6)-(3.7), with e ih := u i -u ih 2 , C ih 4 e 2 ih -C ih e ih + ∇r (ih) 2 ≥ 0. One then concludes as in Lemma 3.1.

  conforming residual liftings following [3, Section 5.1], [51, Section 4.1.1], [22, Section 3.3] will allow us to construct r ih ∈ X h ⊂ V leading to a lower bound on Res(u ih , λ ih ) -1 .

-1 2 1r

 2 (ih) , the Cauchy-Schwarz inequality, estimates (5.8), and the characterization (2.5b) altogether give r (ih) ≤ C I C S h δ Res(u ih , λ ih ) -1 +λ -

  Theorem 4.3], [16, Theorem 3.2], or [39, Theorem 2.1]. Remark 5.10 (Polynomial-degree robustness). The multiplicative factor in the parenthesis in (5.16) takes the form ∇

λ 1 =

 1 49.2883, λ 2 = 98.4296, λ 3 = 127.937, λ 4 = 166.975, λ 5 = 196.439.

λ 1 = 9 .

 19 60692, λ 2 = 15.1695, λ 3 = 19.6932, λ 4 = 29.4166, λ 5 = 31.7363.

Figure 5

 5 Figure 5 plots the convergence of the errors the estimators whereas Table8provides more details on the efficiency. We now observe that the efficiency also degrades for the upper bound of the eigenvalue and eigenvector error. Further, improved lower bounds of the eigenvalue error are not available for the considered meshes for i > 1. This appears as the resulting ηi are all equal to zero, see (5.12), respectively (5.13), so that our eigenvalue upper bound stays that of the finite element method. For all meshes and all considered eigenvalues, though, our estimates still give a rather tight guaranteed eigenvalue confidence interval and quite reasonable eigenvector effectivity indices. We can also observe by a jump of the blue curve (η 2 i ) the change between the cases A and B. The critical mesh size where this change occurs seems to degrade with increasing eigenvalues.

4 ∇(u 4 -Figure 4 :

 444 Figure 4: [Triangular domain, structured meshes] Errors in the first four eigenvalue and eigenvector approximations, their lower bounds (eigenvalues only), and their upper bounds

N

λ 2h -λ 2 η 2 ∇(u 2 -λ 3h -λ 3 η 3 ∇(u 3 -λ 4h -λ 4 η 4 ∇(u 4 -Figure 5 :

 2233445 Figure 5: [L-shaped domain, unstructured meshes] Error in the first four eigenvalue and eigenvector approximations, their lower bounds (eigenvalues only), and their upper bounds

  Moreover, the key result of Braess et al.[START_REF] Braess | Equilibrated residual error estimates are p-robust[END_REF] Theorem 7], see[START_REF]Stable broken H 1 and H(div) polynomial extensions for polynomial-degree-robust potential and flux reconstruction in three space dimensions[END_REF] Theorem 2.2] for three space dimensions, states that the reconstructions of Definition 4.1 satisfy the following stability property,

3.1], Braess et al. [8, Section 3], or [22, Lemma 3.12], there exists a constant C cont,PF only depending on the mesh regularity parameter κ T such that ∇(ψ a v) ωa ≤ C cont,PF ∇v ωa ∀v ∈ H 1 * (ω a ), ∀a ∈ V h . (4.7)

  be arbitrary. Then, for σ ih,dis of Definition 4.1 and r ih of Definition 4.2,

Table 1 :

 1 [Unit square, structured mesh] Validation of assumptions (5.2) and (5.9)

		3 0.4714	16	19.04 ( )	-0.64 (×)
	2	4 0.3536	25	21.55 ( )	0.12 ( )
		5 0.2828	36	22.69 ( )	0.40 ( )
	λ 1 =0.5π 2 λ 2 =3π 2	3 0.4714 4 0.3536 5 0.2828	16 25 36	4.233 ( ) 6.743 ( ) 7.887 ( )	-3.49 (×) -0.66 (×) 0.02 ( )

  Table 5 presents more details on the errors and efficiencies.

	N	h	ndof	λ 1	λ 1h	λ 1h -η 2 1 λ 1h -η2 1 I lb λ,eff I ub λ,eff	E λ,rel	I ub u,eff
	10 0.1698	143 19.7392 20.0336	17.9458	20.6491	-	7.09 1.40e-01	2.65
	20 0.0776	523 19.7392 19.8139	19.6366	19.7909	3.24	2.37 7.83e-03	1.54
	40 0.0413	1975 19.7392 19.7573	19.7307	19.7434	1.30	1.47 6.42e-04	1.21
	80 0.0230	7704 19.7392 19.7436	19.7383	19.7396	1.10	1.20 6.41e-05	1.09
	160 0.0126 30666 19.7392 19.7403	19.7391	19.7393	1.07	1.12 1.04e-05	1.06

Table 3 :

 3 [Unit square, unstructured mesh] Lower and upper bounds on the exact eigenvalue λ 1 , the effectivity indices, and size of the relative λ 1 confidence interval; λ 1 = 1.5π 2 , λ 2 = 4.5π 2 ; Case C

	Level	h	ndof λ 2 -λ 1h (5.2)	α 1h λ 1h /2λ 1 (5.6)	χ 1 -1 (u 1h , χ 1 ) -α 1h (5.6)
	1 0.7500	22	1.8223 (×)	2.97 (×)	-6.17 (×)
	4 0.7071	34	3.8799 ( )	0.94 ( )	-1.27 (×)
	10 0.5000	140	5.2053 ( )	0.33 ( )	0.13 ( )

Table 4 :

 4 [L-shaped domain, adaptive mesh refinement] Validation of the assumptions (5.2) and (5.6) for λ 1 = π 2 /2 and λ 2 = 15.1753 7.3 First eigenvalue on a domain with a hole: mesh adaptivity

Table 5 :

 5 [L-shaped domain, adaptive mesh refinement] Lower and upper bounds on the first exact eigenvalue λ 1 , the effectivity indices, and the size of the relative λ 1 confidence interval; λ 1 = π 2 /2 and λ 2 = 15.1753, Case B

Table 6 :

 6 Domain with a hole, adaptive mesh refinement] Mesh of the adaptive algorithm at iteration 20 (left) and the lower and upper bounds for the exact eigenvalue λ 1 (right); Case B 60.541 60.494 60.455 60.422 60.401 60.387 60.376 60.367 60.359 60.354 λ 1h -η 2 1 37.223 43.710 48.428 52.058 54.275 55.680 56.799 57.719 58.436 58.910 [Domain with a hole, adaptive mesh refinement] Lower and upper bounds on the exact eigenvalue λ 1 as a function of the degrees of freedom; Case B

	1 2 (η 2 1 -η2 1 )/|λ 1 -λ 1h + 1 2 (η 2 1 + η 2 1 )| is

To facilitate the comparisons, we finally present in Tables 9 and 10 several methods for the tests of

[START_REF] Liu | Verified eigenvalue evaluation for the Laplacian over polygonal domains of arbitrary shape[END_REF][START_REF] Armentano | Asymptotic lower bounds for eigenvalues by nonconforming finite element methods[END_REF]

Table 5 .

 5

2 (h = 1/64) and Table

5

.3 (h = 1/32)]. We compare in particular the approach presented in this article, lowest-order conforming finite elements from

[START_REF] Liu | Verified eigenvalue evaluation for the Laplacian over polygonal domains of arbitrary shape[END_REF]

, and the lowest-order Crouzeix-Raviart (CR) method presented in [16], with explicit upper bound of the interpolation constants derived in either

[14]

Table 7 :

 7 [Triangular domain, uniform mesh refinement] Lower and upper bounds on the first four exact eigenvalues λ i , the effectivity indices, and the sizes of the relative λ i confidence intervals

	N	h	ndof	λ i	λ ih	λ ih -η 2 i	λ ih -η2 i	I lb λ,eff	I ub λ,eff	E λ,rel	I ub u,eff
	40 0.0839	1437	9.6397	9.6955	9.1450	9.6870	6.59	9.87 5.76e-02	3.16
	80 0.0459	5674	9.6397	9.6588	9.4719	9.6559	6.52	9.78 1.92e-02	3.15
	160 0.0234	21878	9.6397	9.6467	9.5779	9.6456	6.58	9.86 7.04e-03	3.16
	320 0.0125	86810	9.6397	9.6423	9.6167	9.6419	6.63	9.93 2.62e-03	3.14
	640 0.0059 352256	9.6397	9.6407	9.6310	9.6406	6.73	9.98 9.94e-04	2.77
	40 0.0839	1437 15.1973 15.2440	14.2080 15.2440	-	22.17 1.64e-01	4.70
	80 0.0459	5674 15.1973 15.2092	14.9577 15.2092	-	21.11 4.09e-02	4.60
	160 0.0234	21878 15.1973 15.2002	15.1378 15.2002	-	20.87 1.02e-02	4.57
	320 0.0125	86810 15.1973 15.1980	15.1825 15.1980	-	20.81 2.55e-03	4.55
	640 0.0059 352256 15.1973 15.1974	15.1936 15.1974	-	20.81 6.36e-04	4.09
	40 0.0839	1437 19.7392 19.8216	18.7524 19.8216	-	12.97 1.75e-01	3.59
	80 0.0459	5674 19.7392 19.7597	19.5056 19.7597	-	12.38 4.44e-02	3.52
	160 0.0234	21878 19.7392 19.7444	19.6805 19.7444	-	12.23 1.14e-02	3.50
	320 0.0125	86810 19.7392 19.7405	19.7246 19.7405	-	12.19 2.84e-03	3.48
	640 0.0059 352256 19.7392 19.7395	19.7356 19.7395	-	12.20 7.01e-04	3.07
	40 0.0839	1437 29.5215 29.7057 -154.6818 29.7057	-1000.68	-31.53
	80 0.0459	5674 29.5215 29.5675	-10.5379 29.5675	-	871.81 3.08e+00 29.51
	160 0.0234	21878 29.5215 29.5331	26.5255 29.5331	-	258.67 2.59e-01 16.08
	320 0.0125	86810 29.5215 29.5244	28.8467 29.5244	-	231.37 6.37e-02 15.16
	640 0.0059 352256 29.5215 29.5222	29.3595 29.5222	-	225.32 1.56e-02 13.45

Table 9 :

 9 [Triangular domain, structured meshes] Comparison of different methods; CR is the Crouzeix-Raviart method based approach presented in [16] and the constants indicated in the reference Figure6: [Triangular domain, structured meshes, inexact solver] Error in the first eigenvalue and eigenvector approximation, its lower bound (eigenvalues only), and its upper bound for a uniform refinement; the convergence plots for an exact solver are indicated in gray λ 1 = 9.6380 in this work Liu&Oishi[START_REF] Liu | Verified eigenvalue evaluation for the Laplacian over polygonal domains of arbitrary shape[END_REF] CR with[START_REF] Liu | A framework of verified eigenvalue bounds for self-adjoint differential operators[END_REF] CR with[17]

	λ 1 = 49.348	in this work Liu&Oishi [41] CR with [39] CR with [14]
	Lower bound: 49.341	49.254	49.288	49.225
	Upper bound: 49.351	49.400	49.402
	λ 2 = 98.696	in this work Liu&Oishi [41] CR with [39] CR with [14]
	Lower bound: 98.562	98.352	98.430	98.179
	Upper bound: 98.762	98.931	98.944

Table 10 :

 10 [L-shaped domain, structured meshes] Comparison of different methods; CR is the Crouzeix-Raviart method based approach presented in [16] and the constants indicated in the reference approximation of eigenvalue problems, Math. Models Methods Appl. Sci., 13 (2003), pp. 1219-1229.

* This work was supported by the ANR project MANIF "Mathematical and numerical issues in first-principle molecular simulation". The last author has also received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 647134 GATIPOR).