Variational eigenvalue approximation of non-coercive operators with application to mixed formulations in elasticity

We present an abstract framework for the eigenvalue approximation of a class of non-coercive operators. We provide sufficient conditions to guarantee the spectral correctness of the Galerkin scheme and to obtain optimal rates of convergence. The theory is applied to the convergence analysis of mixed finite element approximations of the elasticity and Stokes eigensystems.


Introduction
In many common applications of solid mechanics, mixed formulations derived from the Hellinger-Reissner variational principle perform better than the standard displacement-based formulation. They deliver direct and accurate approximations of the stress tensor and they are free from the locking phenomenon in the nearly incompressible case [11].
The symmetry constraint on the Cauchy stress tensor has been the main difficulty in the construction of stable conforming discretizations of stress-displacement mixed formulations. The first important progress in this direction is due to Arnold and Winther [3]. This work led to further developments in conforming mixed finite elements on simplicial and rectangular meshes for both 2D and 3D; see [1,6,23] and the references therein. However, these mixed finite elements require the simultaneous imposition of H(div)-conformity and strong symmetry, which entails too many degrees of freedom and complicates the implementation of the corresponding Galerkin schemes. Moreover, they are not amenable to hybridization. To overcome this difficulty one can either consider non-conforming or DG approximations Dedicated to the memory of Francisco-Javier Sayas.
B Salim Meddahi salim@uniovi.es [7,20,32] or relax the symmetry constraint as in [2,4,5,14,21,31]. We point out that the latter alternative, which is the option of choice in this paper, incorporates a further variable (a Lagrange multiplier called the rotor that approximates the skew-symmetric part of the displacement gradient) to enforce weakly the symmetry restriction at the discrete level.
The approximation of eigenvalue problems in mixed form has been the object of several papers; see part 3 of [8] and the references therein. In particular, it is known from [9] that the usual stability conditions for discrete mixed source problems (namely the 'coercivity in the kernel' and the inf-sup conditions) are not sufficient to ensure correct spectral approximations. Recently, a dual-mixed eigenvalue formulation of the elasticity problem with reduced symmetry has been considered in [26]. The eigenproblem resulting from this approach doesn't fit into any of the previously existing theories for mixed spectral problems. Nevertheless, the abstract spectral approximation theory of Descloux-Nassif-Rappaz [16] could be successfully adapted in [26] to show that the Galerkin method based on the first order Arnold-Falk-Winther element [5] is free from spurious modes and converges at optimal rates for eigenvalues and eigenvectors. The same strategy has been applied to a pseudostress formulation of the Stokes eigenproblem [28] and to a stress-pressure formulation of a fluid-structure interaction spectral problem [27].
The aim of this paper is to provide a general theory for the spectral approximation of a class of symmetric and noncoercive operators, so that the studies carried out in [26][27][28] all fit into the same framework. The analysis given here is performed according to the ideas in [26] and builds on the theory developed in [16,17]. The resulting unified approach reveals a new criterion (see Assumption 5 below) to determine the spectral correctness of a given Galerkin approximation. This allows to validate more families of mixed finite elements for the approximation of the elasticity eigenproblem as mentioned in Remark 6.1.
We also highlight that the analysis considered in [26] relies on the regularity of an auxiliary elasticity source problem. Here, we can circumvent the use of this property, which allows us to treat the important case of heterogeneous material coefficients.
The paper is organized as follows: In Sect. 2 we set out the abstract spectral problem and we describe its continuous Galerkin approximation in Sect. 3. In Sect. 4 we provide sufficient conditions ensuring the spectral correctness of the approximation in the sense of [16]. In Sect. 5, we establish rates of convergence for eigenvalues and eigenfunctions. Section 6 is devoted to applications. We show that the abstract framework can be applied to the stress formulation with weak symmetry of the elasticity and Stokes eigenproblems. We present numerical results for the latter example that confirm the theoretical convergence rates.
with the a priori estimate (see [18,Theorem 2.6] where B stands for the norm of the bilinear form B. It follows that the direct decomposition u = u 0 + u − u 0 into components u 0 ∈ K and u − u 0 ∈ K ⊥ B is stable. As a consequence of Proposition 2.1, there exists a unique continuous projector P : X → X with range K ⊥ B and kernel K . We are now going to provide a description of the spectrum of T under the following conditions. Assumption 2 We assume that (i) the inclusion P(X ) → H is compact, (ii) and the inclusion P(X ) ∩ T (X ) → X is compact.
We notice that, as P(X ) = K ⊥ B is T -invariant, the inclusion T (P(X )) ⊂ P(X ) ∩ T (X ) holds true and Assumption 2 (ii) implies that T : K ⊥ B → K ⊥ B is compact. The following result is then a consequence of the spectral characterization of compact operators.
is a sequence of finite multiplicity eigenvalues of T that converges to 0 and the corresponding eigenspaces lie in K ⊥ B ; iii) if T is non-injective, η = 0 is an eigenvalue of T with associated eigenspace ker(T ).

Remark 2.1 If we assume that
, then it can also be shown that the ascent of each eigenvalue η k ∈ (0, 1) is 1, c.f. [26,Proposition A.2].

A continuous Galerkin discretization
We introduce a family {X h } h≥0 ⊂ X of finite dimensional subspaces of X . The continuous Galerkin discretization of the variational eigenproblem 2.1 reads as follows: find 0 = u h ∈ X h and κ h ∈ R such that (3.1) We will use the notation for the distance in X between an element u and a closed subspace W ⊂ X .
To proceed with the analysis of problem (3.1) we need the following discrete inf-sup conditions.

Assumption 4
We assume that Under Assumption 3 and Assumption 4 (ii), we can prove (as in Proposition 2.1) that the splitting is direct and uniformly stable with respect to h. We can also associate to this direct decomposition a unique projector P h : X h → X h with range K ⊥ B h and kernel K h , which is uniformly bounded with respect to h.
Moreover, thanks to Assumption 4 (i), the linear operatorT h : X → X h defined, for all u ∈ X , by We point out that T h :=T h | X h reduces to the identity on K h , which means that 1 is an eigenvalue of T h with associated eigenspace K h . Moreover, κ h = 0 is an eigenvalue of Problem (3.1) if and only if η h = 1/κ h is an eigenvalue of T h and the corresponding eigenspaces are the same. Finally, here again, the symmetry of We are then in a position to provide the following spectral decomposition of T h .
Proof The result follows from the decomposition then, the eigenvalues η hk ∈ (0, 1) are non-defective.

Correctness of the spectral approximation
Henceforth, given any positive functions F h and G h depending on the parameter h, the abbreviation F h G h means that F h ≤ C G h with a constant C > 0 independent h. Moreover, the norm of a linear and continuous operator L : V 1 → V 2 between two Hilbert spaces V 1 and V 2 is denoted The spectral approximation theory developed in [16] for non-compact operators relies essentially on the condition to prove that T h : X h → X h provides a correct spectral approximation of T (in a sense that will be precised in Theorem 4.1 below). The aim of this section is to show that the following key assumption guarantees (4.1).

Assumption 5
There exists a linear operator h : Proof Taking into account that T − T h vanishes identically on K h ⊂ K we obtain, Next, we deduce from the triangle inequality and Céa estimate (3.3) that and the uniform boundedness of T h with respect to h gives the result.
To achieve (4.1), let us first prove the following auxiliary result.

Lemma 4.2
Under Assumptions 1 (ii), 3, 4 (ii) and 5 (iii) it holds, Proof Let us first notice that, by virtue of Assumption 5 (iii), The triangle inequality yields where the last estimate is a consequence of (4.2), Assumption 5 (iii) and the fact that Next, we use the inf-sup condition provided by Assumption 4 (ii) to deduce from (4.3) and Assumption 3 that and the result follows.

Lemma 4.3 Under Assumptions 1-5 it holds
Proof Let us first notice that Combining Lemma 4.1 and Lemma 4.2 with the last estimate yields Now, by virtue of Assumption 2 (ii) and Assumption 5 (i)-(ii), T P : X → X is compact and the operator I − h : (K ⊥ B , · X ) → X is uniformly bounded and converges pointwise to zero. Hence, (I − h )T P : X → X converges uniformly to zero; namely, For the sake of completeness, we finalize this section by adapting the results of [16] (see also [26]) to show that Assumptions 1-5 are sufficient to ensure the correctness of the spectral approximation. Let us first recall that the resolvent the operator of T is given by The mapping z → R z (T ) L(X ) is continuous for all z / ∈ sp(T ) and goes to zero as |z| → ∞. Consequently, it is bounded on any compact subset F ⊂ C satisfying F ∩sp(T ) = ∅, namely, there exists a constant C F > 0 such that where C F > 0 is the constant appearing in (4.8). For E and F closed subspaces of X , we set the latter being the so called gap between subspaces E and F. Let F ⊂ C\{0, 1} be a compact set whose boundary is a smooth Jordan curve not intersecting sp(T ). It is well known [24] that the linear and bounded operator is a projector onto the finite dimensional space E(X ) spanned by the generalized eigenfunctions associated with the finite set of eigenvalues of T contained in . It follows from Lemma 4.4 that, for h small enough, the linear operator The following auxiliary result is essential for this purpose.
Proof We reproduce here the proof given in [26,Lemma 2]. By virtue of Lemma 4.4, there exists h 0 > 0 such that the resolvent identity is satisfied. Hence, for any v h ∈ X h , and the result follows from (4.8) and Lemma 4.4.

Lemma 4.6
Assume that Assumptions 1-5 are satisfied. There exists h 0 > 0 such that Combining the last estimate with (4.9) gives

Consequently, by virtue of the uniform boundedness of
It follows that and the result is a consequence of the last estimate, Lemma 4.5 and (4.12).
We are now in a position to establish the convergence properties of the eigenvalues and eigenfunctions.

Moreover, if E(X ) is the T -invariant subspace of X spanned by the generalized eigenfunctions corresponding to the set of eigenvalues {η
Proof We deduce from Lemma 4.6, Lemma 4.3 and Assumption 5 (ii) and from the fact that As a consequence, E(X ) and E h (X h ) have the same dimension provided h is sufficiently small, c.f. [24]. Finally, as the eigenvalues {η F 1 , . . . , η F m } are isolated, for a sufficiently small > 0, we can consider

Asymptotic estimates for the eigenvalue and eigenfunction error
We proved in Sect. 4 that, under Assumptions 1-5, the Galerkin scheme (3.1) does not pollute the spectrum of T with spurious modes. Moreover, we established the convergence of eigenvalues and eigenfunctions with correct multiplicity. However, in practice the space E η (X ) of generalized eigenfunctions corresponding to a given isolated eigenvalue η = 1 enjoys individual smoothness properties. Therefore, in order to be able to claim that the Galerkin method (3.1) has optimal convergence rates we need to estimate the error for a particular eigenvalue η = 1 and for the corresponding eigenspace E η (X ) only in terms of δ(E η (X ), X h ). This question has been addressed in [17] for noncompact operators under the condition of coercivity for the bilinear form A. In this section we extend the results to the abstract framework we are considering here. Hereafter, we focus on a particular isolated eigenvalue η = 1 of T of algebraic multiplicity m and let D η ⊂ C be a closed disk centered at η with boundary γ such that D η ∩sp(T ) = {η}. We denote by E η := 1 2πi γ R z (T ) dz : X → X the projector onto the eigenspace E η (X ) of η and we define, for h small enough, the projector by E η,h : We begin our analysis by proving an analogue of Lemma 4.
Lemma 5.1 Assume that Assumptions 1-5 are satisfied. Let F ⊂ C be an arbitrary compact subset such that F ∩ sp(T ) = ∅. There exist C F > 0 and h 0 > 0 such that, The last estimate and the triangle inequality yield and the result follows from the uniform boundedness of T h L(X ) .
It follows from Lemma 5.1 that, for h small enough, R z (T h ) : X −→ X is linear and bounded uniformly in h for all z ∈ γ . Hence, the linear operator is the eigenspace corresponding to the eigenvalues of T h contained in γ .

Theorem 5.1 Under Assumptions 1-5 and for h small enough, it holds
Now, due to the fact that E η (X ) ⊂ K ⊥ B is finite dimensional and T -invariant, we deduce from (5.2), Céa estimate (3.3) and Assumption 5 (ii) that It follows that On the other hand, we also deduce from (5.3) that the operatorẼ η,h : E η (X ) → E η,h (X h ) converges uniformly to the identity, which proves that it is invertible, for h small enough. We denote its inverse η,h : ,X ) ≤ 2 and, again by (5.3), The result follows from the last estimate and (5.4).
We recall that κ = 1/η an eigenvalue of Problem 2.1 with the same m-dimensional eigenspace E η (X ). Analogously, if η i,h , i = 1, . . . , m, are the eigenvalues of T h (repeated accordingly to their respective algebraic multiplicities) that converge to η then, κ i,h = 1/η i,h are the eigenvalues of Problem 3.1 converging to κ and the corresponding generalized eigenfunctions span E η,h (X h ). The last step of this section is the following theorem, in which we establish a double order of convergence for the eigenvalues. To this end we need the following assumption.
Proof We denote by u i,h an eigenfunction corresponding to κ i,h satisfying u i,h X = 1. There exists an eigenfunction u ∈ E η (X ) satisfying It follows that, for h small enough, u X is bounded from below and above by a constant independent of h. Furthermore, Assumption 6 and the fact that E η (X ) is finite-dimensional imply the existence of c > 0, independent of h, such that B(u, u) ≥ c u X for all u ∈ E η (X ). Using (5.5) and the uniform boundedness of u X , it is straightforward deduce that B(u ih , u ih ) ≥ c 2 for h sufficiently small. We can now use the identity to obtain the estimate and the result follows.

Applications
We present two applications of the abstract theory developed in the previous sections. They concern dual mixed formulations for the elasticity and Stokes eigensystems. Along this paper we convene to apply all differential operators row-wise. Hence, given regular tensors σ : → M and vector fields u : → R d , we set the divergence div σ : → R d , the gradient ∇u : → M, and the linearized strain tensor ε(u) : → S as

The continuous problem
Let ⊂ R d (d = 2, 3) be a bounded Lipschitz polygon/polyhedron representing a linearly elastic body with mass density ∈ L ∞ ( ) satisfying (x) ≥ 0 > 0 a.e. in . For simplicity, we assume that the structure is fixed at the boundary ∂ . We denote by A(x) : M → M the symmetric and positive-definite 4th-order tensor (known as the compliance tensor) that relates the Cauchy stress tensor σ to the strain tensor through the linear material law A(x)σ = ε(u).
Our aim is to find natural frequencies ω ∈ R such that div σ + ω 2 (x)u = 0 in . Here, we opt for combining this equilibrium equation with the constitutive law to eliminate the displacement field u and impose σ as a primary variable. This procedure leads to the following grad-div eigensystem: Find 0 = σ : → S, 0 = r : → K and eigenmodes ω ∈ R such that, We notice that we introduced above the skew symmetric tensor r := 1 2 ∇u − (∇u) t (the rotation) by writing Hooke's law Aσ = ∇u − r. This additional unknown will act as a Lagrange multiplier for the symmetry restriction.
We point out that the continuous inclusion X → H is not compact.
Thanks to Corollary 6.1, we can define the source operator T : X → X in terms of problem (2.2). We recall that ker(I − T ) = K and that the symmetry of T with respect to B(·, ·) yields T (K ⊥ B ) ⊂ K ⊥ B . In addition, the direct and stable splitting X = K ⊕ K ⊥ B holds true. Our aim now is to characterize the unique projector P : X → X with range K ⊥ B and kernel K associated to this splitting.
For any (σ , r) ∈ X , we consider P(σ , r) := ( σ , r) with σ = A −1 ε( u) and r := 1 2 ∇ u − (∇ u) t , where u is the unique solution of the classical displacement based variational formulation of the elasticity problem in with volume load div σ , namely, We point out that div σ = div σ by construction. In addition, Korn's inequality provides the stability estimate which ensures the continuity of P : X → X . Now, it is clear that P • P = P and ker P = K . Besides, for any (σ , r) ∈ X , , r), (τ , s)) = B (( σ , r), (τ , s) which proves that P(X ) ⊂ K ⊥ B . Finally, we notice that (I − P)X ⊂ K , and hence, We conclude that P : X → X is indeed the unique continuous projector corresponding to the direct and stable decomposition X = K ⊕ K ⊥ B .

Lemma 6.2 Assumption 2 is satisfied.
Proof Let {(σ n , r n )} n be a weakly convergent sequence in X . As P ∈ L(X ), the sequence {( σ n , r n )} n := {P((σ n , r n ))} n is also weakly convergent in X . By definition, σ n = A −1 ε( u n ) and r n := 1 2 ∇ u n − (∇ u n ) t , where u n ∈ H 1 0 ( , R d ) solves (6.4) with righthand side div σ n . It follows from (6.5) that u n is bounded in H 1 0 ( , R d ) and the compactness of the embedding H 1 ( , R d ) → L 2 ( , R d ) implies that { u n } n admits a subsequence (denoted again { u n } n ) that converges strongly in L 2 ( , R d ). Next, we deduce from the Green identity that { σ n } n is a Cauchy sequence in L 2 ( , M). Moreover, the identity (A( σ p − σ q ), τ ) + (τ , r p − r q ) = 0, ∀τ ∈ H (div 0 , , M) and the inf-sup condition (6.3) yield The last estimate ensures that { r n } n is also a Cauchy sequence in L 2 ( , M). We then come to the conclusion that the image under P of any bounded sequence in X contains a converging subsequence in H , which proves that Assumption 2 (i) is satisfieded. On the other hand, testing (2.2) with (τ , 0) and choosing the components of τ : → M indefinitely differentiable and compactly supported in , we readily obtain that, if (σ * , r * ) := T ((σ , r)) then

Consequently,
and the compactness of the embedding T (X ) ∩ P(X ) → X follows. We conclude that Assumption 2 (ii) is fulfilled.
Finally, we notice that for all 0 = (σ , r) ∈ P(X ), , r), (σ , r)) ≥ B((σ , r), (σ , r)) = (Aσ , σ ) > 0, (6.6) where we used that σ is symmetric and that {0} × L 2 ( , K) ⊂ K . We can invoke now Theorem 6.1 and Remark 2.1 to conclude that we have the following spectral characterization for the source operator T . 1) is a sequence of finite-multiplicity eigenvalues of T that converges to 0. The ascent of each of these eigenvalues is 1 and the corresponding eigenfunctions lie in P(X ). Moreover, η = 1 is an infinite-multiplicity eigenvalue of T with associated eigenspace K and η = 0 is not an eigenvalue.

The discrete problem
We consider a family {T h } h of shape regular simplicial meshes of¯ satisfying the standard finite element conformity assumptions. We denote by h K the diameter of triangles/tetrahedra K ∈ T h and let the parameter h : Hereafter, given an integer m ≥ 0 and D ⊂ R d , P k (D, E) is the space of functions with domain D and values in E, where E is either R d , M or K, and whose scalar components are polynomials of degree at most m. Likewise, the spaces of E-valued functions with piecewise polynomial scalar components of degree ≤ m relatively to T h are defined by H (div, , M). We point out that the set {W h , P k (T h , R d ), P k (T h , K)} constitutes the mixed finite element of Arnold-Falk-Winther [5] for linear elasticity. The key property ensuring the stability of this triplet of spaces is given by the following result, c.f. [5,Theorem 7.1].

that is uniformly bounded with respect to h and that satisfies
We point out that H (div 0 , , M). The same procedure used in the proof of Corollary 6.1 can be used verbatim to deduce from Lemma 6.3 that Assumption 4 is satisfied. Moreover, as Consequently, the discrete counterpart of (6.6) is satisfied. Indeed, for all 0 = (σ h , r h ) ∈ P h (X h ) we have that A ((σ h , r h ), (σ h , r h )) ≥ B((σ h , r h ), (σ h , r h ) We deduce from Theorem 3.1 and Remark 3.1 the following result.

Proposition 6.2 The spectrum of T h consists of m := dim(X h ) eigenvalues, repeated accordingly to their respective multiplicities. It holds sp(T
The eigenspace associated to η h = 1 is K h . The real numbers η hk ∈ (0, 1), k = 1, . . . , m 0 , are non-defective eigenvalues with eigenspaces lying in K ⊥ B h and η h = 0 is not an eigenvalue.
Let us now recall some well-known approximation properties of the finite element spaces introduced above. Given s > 1/2 the tensorial version of the canonical interpolation operator h : H s ( , M) → W h associated with the Brezzi-Douglas-Marini (BDM) mixed finite element [12], satisfies the following classical error estimate, see [11,Proposition 2.5.4], Moreover, we have the well-known commutativity property,  H s ( , M). In such a case, we can directly define the operator linear operator h : K ⊥ B → X h by h (σ , r) := ( h σ , S h r) and deduce (as shown below) that Assumption 5 is satisfied. However, instead of relying on regularity results that are difficult to establish for the elasticity system in the case of general domains, boundary conditions and material properties, we resort to the quasi-interpolation operator constructed in [13,19,25] by combining the BDM interpolation operator h with a mollification technique. The resulting projector has domain H (div, , M) and range W h and preserves all the properties of h listed above. More precisely, we will use tensorial version of the following result [19, Theorem 6.5] (see also [13,25]). σ for all σ ∈ H (div, , M).  K). Indeed, by virtue of property (iii) of Theorem 6.1 and because div σ = div σ h , we have that

Proof For an arbitrary
This proves Assumption 5 (iii). On the other hand, using this time property (ii) of Theorem 6.1, we deduce that, if we let ( σ , r) = P(σ , r) ∈ K ⊥ B for an arbitrary (σ , r) ∈ X , we obtain and again by property (iii) of Theorem 6.1 It follows immediately from (6.12)-(6.13) and the triangle inequality that Assumption 5 (i) is satisfied. Moreover, the error estimates (6.8), (6.10) and (6.11) and classical density results ensure that Assumption 5 (ii) is a consequence of (6.13).
We conclude that the Galerkin method (3.1) provides a correct spectral approximation of the eigenproblem (6.1) in the sense of Theorem 4.1.
We recall that η ∈ (0, 1) is an eigenvalue of T with multiplicity m if and only if κ = 1/η is an eigenvalue of Problem 2.1 with the same multiplicity and the corresponding eigenfunctions coincide. Analogously, η i,h , i = 1, . . . , m, are the eigenvalues of T h (repeated accordingly to their respective multiplicities) that converge to η if and only if κ i,h = 1/η i,h are the eigenvalues of Problem 3.1 converging to η. Moreover, the corresponding eigenfunctions coincide. Taking into account that Assumption 6 is satisfied because of (6.6), the following rates of convergence for the eigenfunctions and eigenvalues are a direct consequence of Theorem 5.1, (6.6) and Theorem 5.2, together with the interpolation error estimates (6.8), (6.10) and (6.11).

Theorem 6.2 Let η = 1 be an eigenvalue of T of algebraic multiplicity m and let E η (X )
be the corresponding eigenspace. There exists h 0 > 0 such that for all h ≤ h 0 , T h admits exactly m eigenvalues η i,h , i = 1, . . . , m, repeated according to their respective multiplicity, such that Remark 6.1 For the sake of brevity and clarity of exposition, we only considered here an approximation based on the Arnold-Falk-Winther element [5]. However, we could equally have defined the Galerkin method in base of the families of mixed finite elements introduced by Cockburn-Gopalakrishnan-Guzmán (CGG(k)) in [14] and by Gopalakrishnan-Guzmán (GG(k)) in [21] for k ≥ 1. Indeed, the same BDM quasi-interpolation operator J h can be used in the proof of Assumption 5 for GG(k) and the quasi-interpolation operator corresponding to the Raviart-Thomas element can be used instead for CGG(k).
The following Poincaré-Friedrichs inequality is essential in our analysis.

Corollary 6.2 The bilinear form
is coercive on H (div, , M).
We point out that Corollary 6.2 also implies the coerciveness of the bilinear form on H (div 0 , , M). Consequently, following the same steps given in the proof of Corollary 6.1 we deduce that Assumption 1 is satisfied. Hence, the source operator T : X → X defined by problem (2.2) is bounded and symmetric with respect to B(·, ·). In addition, we have the direct and stable splitting X = K ⊕ K ⊥ B , where K := H (div 0 , , M) × L 2 ( , K) is the eigenspace corresponding to the essential eigenvalue η = 1 of (2.1). We also have that , r), (τ , s)) = 0 for all (τ , s) ∈ X then σ = σ t , (σ , 1) = 0 and Taking τ = σ D in (6.18) we deduce that σ D = 0. This condition holds true if and only if σ = pI , with p = 1 d tr(σ ) ∈ L 2 ( ) and ∇ p = div σ ∈ L 2 ( , R d ). Testing now (6.18) with skew symmetric tensors and using a density argument we deduce that r = 0. It follows that ker(T ) = {q I ; q ∈ H 1 ( ) ∩ L 2 0 ( )} × {0} is the eigenspace corresponding to the eigenvalue η = 0 of T .
The characterization of the projector P : X → X associated to this decomposition follows the same pattern of last section. For any (σ , r) ∈ X , we consider P(σ , r) = ( σ , r) ∈ X with σ := ε(ũ) −pI and r := 1 2 ∇ũ − (∇ũ) t , where (ũ,p) ∈ H 1 0 ( , R d ) × L 2 0 ( ) is the unique solution of the classical velocity-pressure variational formulation of Stokes problem in with load − div σ , i.e., (6.19) We know [11] that there exists a constant C > 0 such that moreover, div σ = div σ by construction. This ensures the continuity of P : X → X . Besides, it is clear that P • P = P, ker P = K , and for any (σ , r) ∈ X , , r), (τ , s) which proves that P(X ) ⊂ K ⊥ B . The inclusion K ⊥ B ⊂ P(X ) is a consequence of (I − P)X ⊂ K . We conclude that P : X → X is the unique continuous projector associated to the direct and stable decomposition X = K ⊕ K ⊥ B .

Lemma 6.6 Assumption 2 is satisfied.
Proof Let {(σ n , r n )} n be a weakly convergent sequence in X . The continuity of P : X → X implies that {( σ n , r n )} n := {P(σ n , r n )} n = { ∇ũ n −p n I , 1 2 ∇ũ n − (∇ũ n ) t } n also converges weakly in X , where (ũ n ,p n ) ∈ H 1 0 ( , R d ) × L 2 0 ( ) solves (6.19) with datum − div σ n . As a consequence of (6.20), {ũ n } n is bounded in H 1 ( , R d ). In turn, we deduce from the identities −∇p n = div σ n = div σ n that {p n } n is also bounded in H 1 ( ). The compactness of the embedding that {σ n ,p n } n admits a subsequence (also denotes {σ n ,p n } n ) that converges strongly in L 2 ( , R d ) × L 2 ( ). Next, it follows from (6.19) and a Green formula that Taking v =ũ p −ũ q in the last identity yields which proves that {ε(ũ n )} n is a Cauchy sequence in L 2 ( , M). Moreover, testing the identity We can now use the inf-sup condition (6.3) as in Lemma 6.2 to deduce that { r n } n is also a Cauchy sequence in L 2 ( , M). We conclude that {( σ n , r n )} n := {P(σ n , r n )} n admits a subsequence that converges strongly in H and Assumption 2 (i) follows.
Let us denote by (σ * , r * ) := T (σ , r) the image of any (σ , r) ∈ X by the source operator T . Testing (2.2) with (I , 0) ∈ X we deduce that (σ * , 1) = (σ , 1) . Testing now the same equation with (τ , 0), and choosing the entries of τ : → M indefinitely differentiable and compactly supported in , we deduce that Therefore, the subspace is compactly embedded in X and Assumption 2 (ii) is also satisfied.

The discrete problem
The Galerkin scheme (3.1) is based on the same finite element spaces used in the last section. Namely, for k ≥ 0, we take X h := W h × P k (T h , K), with W h := P k+1 (T h , M) ∩ H (div, , M) and here again K h = W 0 h ×P k (T h , K) ⊂ K . It follows easily from Lemma 6.3 that Assumption 4 holds true. We define h : X → X h as in Lemma 6.4, which ensures Assumption 5 and guaranties by the way that the Galerkin method (3.1) provides a correct spectral approximation of the stress formulation of the Stokes eigenproblem (2.1) in the sense of Theorem 4.1. Moreover, Assumption 6 is satisfied thanks to (6.21). We can then rely on Theorem 5.2 to obtain the following rates of convergence for the eigenfunctions and eigenvalues.

Remark 6.3
The analysis given in this section can be adapted to deal with Dirichlet-Neumann boundary conditions for eigenproblems (6.1) and (6.15) by defining the space X as in [26] and by employing the quasi-interpolation operator with partial boundary conditions recently introduced by Licht in [25, Theorem 6.3].

Numerical results
On the unit disk, the eigenvalues of the Stokes eigenproblem (6.15) are given by the sequence { 1 2 j 2 n } n≥1, ≥1 , where j nk is the -th positive zero of the Bessel function J n of the first kind of order n. Accurate approximations of the first 4 eigenvalues are given by We use the open-source finite element software Netgen/NGSolve [30] to implement the Galerkin method (3.1) of the Stokes eigenproblem (6.15). We take advantage of the Netgen/NGSolve support for curved mixed finite elements of arbitrary order to base the construction of X h on an H (div)-conforming BDM-parametric element associated to an exact triangulation T h of the unit disk. We notice that this leads to a conforming finite element approximation X h of X .
We denote by λ jh the approximation of λ j computed by solving problem 3.1 with A and B given by (6.17). We introduce the experimental rates of convergence where h andĥ are two consecutive mesh sizes. We present in Table 1 the first four eigenvalues computed on a series of exact partitions of the unit disk¯ with decreasing mesh sizes h, and for polynomial degrees k + 1 = 1, 2, 3, 4. We also report in the same table the arithmetic mean of the experimental rates of convergence obtained for each eigenvalue via (6.24). We observe that a convergence of order 2(k + 1) is attained for each eigenvalue, as predicted by the error estimate given in (6.22).