Benchmark results for testing adaptive (cid:12)nite element eigenvalue procedures part 2 (conforming eigenvector and eigenvalue estimates)

We present an hp-adaptive continuous Galerkin (hp-CG) method for approximating eigenvalues of elliptic operators, and demonstrate its utility on a collection of benchmark problems having features seen in many important practical applications—for example, high-contrast discontinuous coeﬃcients giving rise to eigenfunctions with reduced regularity. In this continuation of our benchmark study, we concentrate on providing reliability estimates for assessing eigenfunction/invariant subspace error. In particular, we use these estimates to justify the observed robustness of eigenvalue error estimates in the presence of repeated or clustered eigenvalues. We also indicate a means for obtaining eﬃciency estimates from the available eﬃciency estimates for the associated boundary value (source) problem. As in the ﬁrst part of the paper we provide extensive numerical tests for comparison with other high-order methods and also extend the list of analyzed benchmark problems.


Introduction
In Part 1 of this work [15], we presented an hp-adaptive discontinuous Galerkin (hp-DG) method for approximating eigenvalues of elliptic operators, and demonstrated its utility on a collection of benchmark problems having features seen in many important practical applications-for example, high-contrast discontinuous coefficients giving rise to eigenfunctions with reduced regularity. This approach was shown to be highly efficient in terms of cost per correct digit, and provided computed error estimates which were asymptotically identical to the actual errors. Our DG work did not, however, provide a means of assessing eigenfunction/invariant subspace error; nor did it justify the observed robustness of eigenvalue error estimates in the presence of repeated eigenvalues. In this continuation of our benchmark study, we present a continuous Galerkin hp-adaptive (hp-CG) method which addresses these issues. As in Part 1, we provide extensive numerical tests for comparison with other high-order methods. Our CG work builds on the abstract framework of [21,7] for eigenvalue/invariant subspace error estimation, as outlined in Section 2, providing a more robust theory than was possible in our previous approach. Concerning hp-adaptive approximation methods for eigenvalue problems, much more has been written about their a priori error analysis than their a posteriori error analysis, and we mention [5,26,31], as well as the more comprehensive general survey [11], as recent references in this vein. A recent contribution to a posteriori analysis which is most readily compared with our own is that of [4], and we make relevant comparisons and contrasts with this work in Sections 3 and 4 .
Our model problem is as follows. Let Ω ⊂ R 2 be a bounded polygonal domain, and let ∂Ω D ⊂ ∂Ω have positive (1D) Lebesgue measure. We define the space H = {v ∈ H 1 (Ω) : v = 0 on ∂Ω D }, where these boundary values are understood in the sense of trace. We are interested in the eigenvalue problem: where We assume that the diffusion matrix A is piecewise constant and uniformly positive definite a.e., and the scalar c is also piecewise constant and non-negative. The assumption that A and c are piecewise constant, as opposed to just piecewise smooth, is a theortical convenience, as is the assumption that Ω is a polygon. These are not practical limitations to our approach. The discontinuities in these coefficients permit investigation of many interesting eigenvalue model problems for composite materials, such as those which are of interest for methods of nondestructive sensing (cf. [1, 2]).
Our interest in problems of this sort is motivated by considerations of photonic crystals and related problems, (cf. [3,14]). Such applications are not directly addressed in this work, but many of their inherent computational challenges are captured in our model problem.
The assumptions on the coefficients A and c are sufficient to guarantee that the bilinear form is both bounded and coercive on H 1 (Ω), B(v, w) ≤ β 1 ∥v∥ 1 ∥w∥ 1 , B(v, v) ≥ β 0 ∥v∥ 2 1 for all v, w ∈ H , for some constants β 0 , β 1 > 0. So B(·, ·) is an inner product on H, whose induced energy-norm ||| · ||| is equivalent to ∥ · ∥ 1 . Here and elsewhere, | · | 1 and ∥ · ∥ 1 denote the usual semi-norm and norm on H 1 (Ω), and ∥ · ∥ 0 denotes the norm on L 2 (Ω). For some results, it is convenient to restrict the norms/semi-norms to a subset S ⊂ Ω (in the obvious way), and we denote these restrictions by ∥ · ∥ 0,S , | · | 1,S , ∥ · ∥ 1,S and ||| · ||| S . The we point out that ||| · ||| S may actually be a semi-norm, though our assumptions on the coefficients A and c guarantee that there are local constants β 0S , β 1S > 0 such that β 0S |v| 2 1,S ≤ |||v||| 2 S ≤ β 1S ∥v∥ 2 1,S , and the seminorm in the lower bound can be replaced with the full norm (after modifying The variational eigenvalue problem (1) is satisfied by the positive sequence of eigenvalues and a sequence of eigenvectors (ψ i ) i∈N , Here we have counted the eigenvalues according to their multiplicity. Furthermore, the sequence (λ i ) i∈N has no finite accumulation points. For path-wise connected domains Ω, we know via the Peron-Frobenius Theorem, that λ 1 < λ 2 holds and that ψ 1 can be chosen to be continuous and strictly positive in Ω. We will also use the notation to denote the spectrum of the variational eigenvalue problem, and the spectral subspace (invariant subspace) associated to λ ∈ Spec B , noting that these subspaces are finite dimensional in our setting. Furthermore, let E λ be the L 2 -orthogonal projection onto M(λ). Then ∑ λ∈Spec B E λ = I and the spaces M(λ) = Ran E λ and M(µ) = Ran E µ are mutually orthogonal for distinct λ, µ ∈ Spec B . We finally note that and so we obtain an alternative representation of the energy norm The rest of this paper is organized as follows. In Section 2 we introduce basic notation related to the hp-discretization, and discuss approximation defects as ideal error estimates, with key results from [21,7] extended for use in the present context. These extensions make possible the incorporation of results from [28, Section 3] to obtain efficient and reliable estimates of eigenvalue and eigenvector approximations. This development is given in Section 3, where our key eigenvalue and eigenvector results, Theorem 3.4 and Theorem 3.8, are presented. Section 4, which constitutes the bulk of the paper, is devoted to numerical experiments on a variety of different kinds of problems to assess the practical behavior of the proposed approach.

Discretization and ideal error estimates
We discretize (1) using hp-finite element spaces, which we now briefly describe. Let T = T h be a triangulation of Ω with the piecewise-constant mesh function h : T h → (0, 1), h(K) = diam(K) for K ∈ T h . Throughout we implicitly assume that all meshes -even on the coarsest level -are aligned with all discontinuities of the data A and c, as well as any locations where the (homogeneous) boundary conditions change between Dirichlet and Neumann. Given a piecewiseconstant distribution of polynomial degrees, p : T h → N, we define the space where P j is the collection of polynomials of total degree no greater than j on a given set. Suppressing the mesh parameter h for convenience, we also define the set of edges E in T , and distinguish interior edges E I , and edges on the Neumann boundary E N (if there are any). Additionally, we let T (e) denote the one or two triangles having e ∈ E as an edge, and we extend p to E by p(e) = max K∈T (e) p(K). As is standard, we assume that the family of spaces satisfy the following regularity properties on T h and p: There is a constant γ > 0 for which It is really just a matter of notational convenience that a single constant γ is used for all of these upper and lower bounds. The shape regularity assumption (C1) implies that the diameters of adjacent elements are comparable.
In what follows we consider the discrete versions of (1): We also assume, without further comment, that the solutions are ordered and indexed as in (3), with (ψ i ,ψ j ) = δ ij . We are interested in assessing approximation errors in collections of computed eigenvalues and associated invariant subspaces. Let s m = {µ k } m k=1 ⊂ (a, b) be the set of all eigenvalues of B, counting multiplicities, in the interval (a, b), and let S m = span{ϕ k } m k=1 be the associated invariant subspace, with (ϕ i , ϕ j ) = δ ij . The discrete problem (6) is used to compute corresponding Given computed approximationsŝ m andŜ m as described above, we define the corresponding approximation defects as: where u(f ) andû(f ) satisfy the boundary value problems: In Theorems 2.2 and 2.4 below, we state key theorems from [20] and its published version [17], which were used for finite element computations in [21,7]. The results of [20,17] show that these approximation defects would yield ideal error estimates for eigenvalue and eigenvector computation if they could be computed. Note that for the results of [20,17] to hold we do not haave to make any assumptions save that B is a positive definite quadratic form and that V is from its associated form domain.
The constant C m is given by an explicit formula which is a reasonable practical overestimate, see [21,7,20,17,16] for details. In particular the requirement ηm(Ŝm) is according to [16,17] a reasonable restriction to guarantee a second order -in thew norm of the residual -pairing between Ritz values and eigenvalues. A similar result holds for the eigenvectors. We point the interested reader to [21,Theorem 4.1 and equation (3.10)] and [7, Theorem 3.10] for the use of the abstract results from [20,17,16,18] in the finite element setting. Specifically, let us emphasize that neither B has to be a divergence type quadratic form nor does V does have to be a finite element space for the results to hold. Subsequently all of the constants which appear are independent of the mesh parameters and are solely problem dependent. Remark 2.3. Although λ 1 < λ 2 for the particular problems we consider numerically in the present work, much of the theory carries over to problems where Ω is not pathwise connected, or the boundary conditions are periodic (as examples). In these cases the Peron-Frobenius theorem does not apply, and it is quite possible that the smallest eigenvalue is degenerate. If this is the case, and λ 1 = λ m , then the constantλ 1 /2λ m in (10) can be replaced by 1.
Another notable feature of these ideal estimates is that they are asymptotically exact, both as eigenvalue and eigenvector error estimators, as the following theorem indicates in the case of a single degenerate eigenvalue and its corresponding invariant subspace.
be the computed approximation of the invariant subspace corresponding to λ q . Then, taking the pairing of eigenvectors ϕ i and Ritz vectorsφ i as in [21], we have Remark 2.5. We emphasize that max-min formulation for the approximation defects has at it heart relative errors for source problems (8)-(9), for which a posteriori error estimation techniques are better developed. It will be clear from the development below that our computable estimates of the approximation defects will inherit the strengths and weaknesses, both theoretical and practical, of the underlying a posteriori techniques for source problems. For a theoretical study of the convergence rates as well as comparisons with other approaches see [17].

Error estimates in a Riesz vector basis
We now turn to practical estimates of the approximation defects, and we must first address their max-min variational definition. Recall that we use u(·) andû(·) to denote the solution operators from (8) and (9). Takingφ 1 , . . . ,φ m to be the Ritz vectors fromŜ m as described above, we define the matrices E, G ∈ R m×m by It was shown in [21] that are the eigenvalues of the generalized eigenproblem for the matrix pair (E, G). Furthermore, since G is positive definite, these eigenvalues are precisely those of At this stage we see that max-min problem is reduced to a small (generalized) eigenvalue problem, but we must still approximate the entries of E and G in a useful way. In [21,7], which concerned low-order h-methods, we used hierarchical basis techniques to approximate the error µ jφ j , not just some norm of them, and these were used to approximate each entry in E and G. We cannot take the analogous hierarchical basis approach for our high-order work, because key theoretical results (e.g. reliability) in the hp-setting are not yet available. We instead employ the residual-based error estimates of [28], for which efficiency and reliability are proven, and these results are adapted for use in our setting in Section 3. The residual-based approach is designed to estimate the norm errors, |||u(φ j )−û(φ j )|||, so we can no longer approximate the off-diagonal entries of E and G in a useful way. To address this, we "diagonalize" the max-min problem as follows, to close out this section. Letting Because it is a relative estimate, it is expected that D l < 1 even for finite element spaces V of fairly small dimension. In any case, we have D ≤ G ≤ (1 + D l )D. Here we use ≤ to denote the Löwner order relation between Hermitian matrices, for further properties and definition see [8,23]. Therefore, or, equivalently, we obtain the following result.
Lemma 2.6. It holds that

Error estimation for hp-approximations of eigenvalues and eigenvectors
Using Lemma 2.6, we have reduced the problem of estimating the approximation defects, and hence the error in our eigenvalue/eigenvector computations, to that of estimating error in associated boundary value problems. In particular, we must estimate |||u(μ iφi )−û(μ iφi )||| 2 = |||u(μ iφi )−φ i ||| 2 for each Ritz vector, whereŜ m = span{φ 1 , . . . ,φ m } is our approximation of S m = span{ϕ 1 , . . . , ϕ m }. We modify key results from [28], which were stated only for the Laplacian, to our context. We define the element residuals R i for K ∈ T , and the edge (jump) residuals r i for e ∈ E, by For interior edges e ∈ E I , K and K ′ are the two adjacent elements, having outward unit normals n K and n K ′ , respectively; and for Neumann boundary edges e ∈ E N (if there are any), K is the single adjacent element, having outward unit normal n K . We note that R is a polynomial of degree no greater than p(K) on K, and r is a polynomial of degree no greater than p(e) on e, because of our assumption that A and c are piecewise constant.
where E I (K) and E N (K) denote the interior edges and Neumann boundary edges of K, respectively. An inspection the proof of [28, Lemma 3.1] (which was stated for the Laplacian) makes the following assertion clear.
Lemma 3.1. There is a constant C > 0 depending only on the hp-constant γ and the coercivity A few remarks are in order concerning the lemma above and how it relates to [28, Lemma 3.1]. First, the bound in [28, Lemma 3.1] includes an additional term involving the difference between the righthand side (in our caseμ i ϕ i ) and its projection on K into a space of polynomials. This additional term only arises in their result because they have chosen to use the projection of the righthand side, instead of the righthand side itself, to define the element residual (here called R i ). They do this in order to employ certain polynomial inverse estimates, which hold in our case outright because our righthand sides are piecewise polynomial. Their result also involves a parameter α ∈ [0, 1], which we have taken to be 0. The result [28, Lemma 3.1] is based on Scott-Zhang type quasi-interpolation, which naturally gives rise to errors measured in H 1 . Mimicking their arguments with our indicator, one would arrive at a result of the form Here, ω K is the patch of elements which share an edge with K. The global continuity constant β 1 could be replaced in Lemma 3.2 by a local continuity constant β 1ω K if desired.
Remark 3.3. The p-dependence in local efficiency bound of Lemma 3.2 is unfortunately unavoidable in the proof, and would suggest decreased efficiency of the estimator as p K is increased if this estimate were sharp. Our numerical experiments do seem to indicate that there may be a modest decrease in the efficiency of the estimator under hp-refinement in practical computations.
With these results we now state the main theorem concerning relative eigenvalue error estimates The constant C 1 depends solely on the ratioλ 1 /(2λ 2 ), the hp-regularity constant γ, the continuity constant β 1 , and the maximal polynomial degreep = max K∈T p(K). The constant C 2 depends solely on the relative distance to the unwanted component of the spectrum, the hp-regularity constant γ and the coercivity constant β 0 .
Proof. LetŜ m = span{φ 1 , . . . ,φ m } be our approximation of S m = span{ϕ 1 , . . . , ϕ m }. We start from the inequality (10) Using (16), we obtain the estimatê From (21) and Lemmas 3.1 and 3.2 the assertions of the theorem now follow directly. In particular due to the uniform estimate D l ≤ 1 we havê From this we readily deduce the claims on the constants C 1 and C 2 .
Remark 3.5. It is relative local indicatorsλ −1 i ε 2 i (K) which will be used to mark elements for refinement, as described in Section 4.
Remark 3.6. The dependence of the efficiency constant C 1 on the maximal polynomial degree is solely due to the conclusion of Lemma 3.2. If one were to use instead of (19) a different boundary value estimator for which the efficiency estimate, similar to Lemma 3.2, involves constants which do not depend on the polynomial degree, then the constant C 1 from Theorem 3.4 would directly inherit such a property. One possible class of estimators which could exhibit such a feature would be a generalization of hierarchical error estimators, cf. [7]. This will be a topic for further research.
A similar result holds for the eigenvectors and eigenspaces. We let be the L 2 orthogonal projection onto the eigenspace belonging to the first m eigenvalues of the form B as given in Theorem 3.4. We also take ∥ · ∥ S 2 to be the Hilbert-Schmidt norm on the ideal of all Hilbert-Schmidt operators (see [32]), which is the natural extension of the matrix Frobenius norm, . Note that we are considering subsets of the space of bounded (and compact) operators from L 2 to L 2 . Our first approximation result is The constant C m,T depends solely on the relative distance to the unwanted component of the spectrum (e.g. λm−λ m+1 λm+λ m+1 ), the hp-regularity constant γ and the continuity constant β 1 .
This result is a robust reliability estimate which ensures the convergence of invariant subspaces. With additional information on eigenvalue separation we present more detailed efficiency and reliability estimate in an eigenvector setting. Let λ 1 = λ s 1 < · · · < λ s k be all elements of the set {λ i : i = 1, · · · , m}. We define the following gap measure

m be eigenvectors and Ritz vectors which satisfy both the assumptions of Theorem 2.2 and the paring of Theorem 2.4, as well as
The constant c V depends solely on the ratioλ 1 /(2λ 2 ), the hp-regularity constant γ, the continuity constant β 1 , and the maximal polynomial degreep = max K∈T p(K The conclusion of the theorem now follows from Lemmas 2.6, 3.1 and 3.2. We close this section with a few remarks concerning [4], which uses a similar eigenvalue/eigenvector approximation error estimator, also based on the estimator from [28]. Although the final practical estimates are similar, the theoretical developments are not. Our approach provides a more robust analysis of the reliability of the eigenvalue/eigenvector estimator for multiple and clustered eigenvalues, whereas theoretical results from [4] were for simple eigenvalues (cf. [4, Proposition 3.1]). Note here that although ϵ i do depend on the chosen Riesz vector basis ofŜ m , the sum does not. According to (14) it is trace type invariant and so it is meaningful to reduce it as a basis independent characteristic of the Riesz spaceŜ m . This is the basis of our basis independent claim for the resolution of eigenvalue multiplicity and for computing a subspace dependent mesh refinement. We note that multiple or clustered eigenvalues appear as result of symmetries or near symmetries of the problem and are a typical feature of 2D or 3D eigenvalue problems. In a later work, [30, Section 7.1], these authors suggest that one might work around the issue of multiple eigenvalues arising due to symmetries by exploiting these symmetries to reduce the problem to a different domain. We see in our work that this is not necessary in theory or practice.

Experiments
In the numerical experiments we illustrate the efficiency of the estimator (20) on several problems of the general form for a second-order, linear elliptic operator L, where homogeneous Dirichlet or Neumann conditions are imposed on the boundary. Following [6], we assume an error model of the form √ DOFs for problems whose eigenvectors are expected to be smooth, and for problems such as those on non-convex polygonal domains and/or discontinuous coefficients, whose eigenvectors are expected to have isolated singularities. We use DOFs = dim(V p k ) to denote the size of the discrete problem. The constants C and α are determined by least-squares fitting, and α is reported for each problem. Plots are given of the total relative error, its a posteriori estimate, and the associated effectivity index, shown, respectively, below: In the case of a single eigenvalue λ i the effectivity index reduces (λ i − λ i )/ε 2 i , and we make the following comparison with what is presented in [4], in which hp-adaptivity is also used for eigenvalue problems. The effectivities reported in [4] are in terms of eigenfunction error, which corresponds closely with the square root of the effectivities reported here. This difference should be taken into consideration when comparing the effectivities reported here with those in [4] or other similar contributions. For problems in which the exact eigenvalues are known, we use these values in our error analysis. For most problems, we use highly accurate computations on very large problems to produce "exact eigenvalues" for our comparisons, as discussed in the introduction.
In all simulations we used an hp-adaptive algorithm in order to get the best convergence possible. To drive the hp-adaptivity we use the element-wise contributions to the quantity ∑ m i=1λ −1 i ε 2 i , to provide local error indicators. Then, we apply a simple fixed-fraction strategy to mark the elements to adapt. For each marked element, the choice of whether to locally refine it or vary its approximation order is made by estimating the local analyticity of the computed eigenvectors in the interior of the element by computing the coefficients of the L 2 -orthogonal polynomial expansion (cf. [13,22]). In contrast, the authors of [4] make this decision by comparing the local indicator on a given marked triangle with a prediction of that error derived from the indicator of its parent.
Our algorithm, presented as Algorithm 1, has a very simple structure that consists of a repeatuntil loop. During each iteration of the loop a new approximation of the eigenpair(s) of interest is computed, then the (group) error estimator is calculated and, if the global group error estimator ∑ m i=1λ −1 i ε 2 i is smaller than the prescribed tolerance tol the algorithm stops; otherwise the mesh T and the space V p k are refined and another iteration follows. The algorithm is also dependent on two parameters θ,θ ∈ ⟨0, 1]. The parameter θ controls the fixed fraction element marking strategy. The parameterθ controls the convergence rate by coupling the p refinement with the h-refinement strategy, for details see [12,13,22]. Our p refinement strategy is, as in [13,22], based on testing for local analyticity.

Algorithm 1 Cluster oriented hp-adaptive algorithm
The function Refine does three things: marks the elements for refinement, decides the refinement pattern for each marked element choosing either h-or p-adaptivity, and refines the finite element space. When h-adaptivity is applied to marked elements, those elements are refined using redrefinement, and red-refinement is also applied to the neighbours to avoid any constraints on the newly created degrees of freedom on the marked elements. This will typically introduce "hanging nodes" on the neighbors of marked elements, which is why we consider one-irregular meshes. On the other hand, when p-adaptivity is applied to some elements, the orders of such elements are increased by one and also the orders of the neighbourhood elements are accordingly increased to avoid any new constrained degrees of freedom.

Dirichlet Laplacian on the unit triangle
As a simple problem for which the eigenvalues and eigenvectors are explicitly known (cf. [27]), we consider the problem where: L = −∆, Ω is equilateral triangle of having unit edge-length, and ψ = 0 on ∂Ω. The eigenvalues can be indexed as λ mn = 16π 2 9 (m 2 + mn + n 2 ) , and we refer interested readers to [27] for explicit descriptions of the eigenvectors. In Figure 2(a) we plot the total relative error for the first four eigenvalues, together with the associated error estimate; and in Figure 2(b) we plot the effectivity quotient. In this case we have obtained α = 0.5070. It is clear that the convergence is exponential in this case, and that the effectivity undergoes a mild degradation as the problem size increases. This modest decrease in effectivity is in line with Remark 3.3, and it is also seen in several of our remaining experiments.

Dirichlet Laplacian on the unit triangle with a hole
We now consider the problem where L = −∆, Ω is the equilateral triangle having edge-length 2 with an equilateral triangle having edge-length 1/2 removed from its center (see Figure 1), and ψ = 0 on ∂Ω. For such a problem, it is expected that some of the eigenvectors will have an r 3/5type singularity at each of the three interior corners, where r is the distance to the nearest corner. In this case, the exact eigenvalues are unknown, so we computed the following reference values of In Figure 3(a) we plot the relative error and error estimates together, for the first three eigenvalues, and in Figures 3(b) we plot the corresponding values of the effectivity quotient. We again see exponential convergence with α = 0.2190 and a modest deterioration of effectivity.

Square domain with discontinuous reaction term
For this pair of problems we take Ω = (0, 1) 2 , ∇ψ · n = 0 on ∂Ω, and L ψ = −∆ψ + κV M D · ψ, where V M D is the characteristic function of the touching squares labelled M 1 in Figure 4. We consider two values of the constant parameter, κ = 10, 100. It is straightforward to see that the corresponding bilinear form is an inner-product in this case (no zero eigenvalues), and that all eigenvectors are at least in H 2 .
For κ = 10, we have in Figure 5(a) the total relative error and error estimates for the first four eigenvalues; and the effectivity quotient is given in Figure 5(b). For these simulations we used the following reference values for the first four eigenvalues, which are 1e-8 accurate: 4.150242455, 10.706070962, 18.779725462, 25.150325247. The analogous plots for the first four eigenvalues in the case κ = 100 are given in Figure 6(a) and Figure 6(b). For these simulations, we used the following reference values for the first four eigenvalues, which are 1e-8 accurate: 13.210576406, 13.990033964, 60.294151672, 64.840268299. In both cases we see apparent exponential convergence with α = 0.2495 and α = 0.1827 respectively, and reasonable effectivity behavior. It is clear from the error plots that for both values of κ the convergence is exponential.

Square domain with discontinuous diffusion term
Using the domain Ω = (0, 1) 2 , partitioned into regions M 1 and M 2 as in Figure 4, and homogeneous Dirichlet conditions ψ = 0 on ∂Ω, we consider the operator L = −∇ · (a∇), where a = 1 in    Figure 9 we reported the final mesh and the final distribution of polynomials orders for κ = 100.

A Kellogg problem
We here consider a variant of the previous problem type for which we can give more specific information about the kinds of singularities which can be expected in terms of the size of the jump    discontinuity. More specifically, we consider problems of the form ∫ where a has jump discontinuities across certain internal interfaces. We refer to this type of problem (24) as a Kellogg eigenvalue problem, in reference to the work of that author on boundary value problems of this sort-although an argument could be made for calling the previous problem type by this name.
If Ω is the unit disk and a = κ = β 2 in the first and third quadrants, and a = 1 in the second and fourth quadrants (see Figure 1), we can describe the eigenpairs explicitly. We assume that β > 1. The eigenvalues and functions can be split into three different classes, which we now describe.

Obviously, Ψ
(1) km is discarded when k = 0. We see that the eigenvalues of this type are independent of β, as are the eigenvectors ψ , θ ∈ [3π/2, 2π) when k is even , when k is odd .
Class 3. For k ≥ 1 and m ≥ 1, let z when k is even , when k is odd .
It is clear from these expressions that singularities of type r γ for any γ ∈ (0, 1) may be achieved by choosing β large enough-these may be obtained by Class 2 eigenvectors when k = 0, for example. If we choose κ = β 2 = 10 for the circle domain, the eigenvectors associated with the smallest three eigenvalues are The second of these has an r σ 0 -type singularity at the origin, where σ 0 ≈ 0.389964; the third of these has an r ρ 1 -type singularity at the origin, where ρ 1 ≈ 1.61004. So it is clear that the second eigenvector is the most singular. We compute eigenvalues on the analogous square domain (Figure 4), with a = 1 in M 2 and a = κ = β 2 = 10 in M 1 . The singular behavior of the eigenvectors near the cross point will be the same as for the circular domain. In Figures 10(a)-10(b) we report the total relative error and error estimates for the first three eigenvalues, and the effectivity index. For these simulations we used the following reference values for the first three eigenvalues: 19.739208802 (1e-8), 30.264820 (1e-5), 70.310149038 (1e-8). Again we see apparent exponential convergence, with α = 0.2624. √ DOF s , where α = 0.2642. We also include the analogous data for pure h-adaptive refinement using quadratic elements to illustrate the difference in performance from the hp version.
To give some indication of the nature of the eigenvectors in the interior, we briefly consider a related problem where Ω is the unit disk with a slit along the positive x-axis, as pictured in Figure 1, with the same boundary conditions. In this case, the eigenvalues and eigenvectors are known explicitly. For k ≥ 0 and m ≥ 1, let z km be the m th positive root of the first-kind Bessel function J k/2 . It is straightforward to verify that, up to renormalization of eigenvectors, the eigenpairs can be indexed by λ km = z 2 km , ψ km = J k/2 (z km r) cos(kθ/2) , k ≥ 0 , m ≥ 1 .
We see that ψ km ∼ cos(kθ/2) ( z km r 2 ) k/2 as r → 0, so singularities of type r k/2 occur infinitely many times in the spectrum. The strongest of these singularities is of type r 1/2 , and it occurs in the eigenvector associated with the second eigenvalue, for example. The same asymptotic behavior of the eigenvectors near the crack tip is expected for the square and circular domains, and in Figure 11 we show a contour plot of the second eigenvalue for the square domain.
In Figure 12 we plot the total relative errors and error estimates for the first four eigenvalues with α = 0.3314, and in Figure 13 the individual eigenvalue errors are shown. It is clear from the second of these figures that the second eigenvalue, which corresponds to the most singular eigenvector, clearly has the worst convergence rate (as expected), and that this is what "spoils" the convergence of the cluster of the first four eigenvalues. This becomes even more apparent when Figure 14 (with α = 0.3121), which corresponds to the second eigenvalue alone, is compared with Figure 12-they are nearly identical. Moreover in Figures 15 and 16 we report the final mesh and the final distribution of polynomials orders for the second eigenvalue. As can be seen, the adaptive procedure has automatically heavily refined in the center, where the singularity is located.    The polynomial degree distribution is not symmetric. The symmetry could be enforced using the symmetry group of the domain, but both the theory and the experiments confirm that we achieve exponential convergence without imposing this further restriction.