A Posteriori Eigenvalue Error Estimation for the Schrödinger Operator with the Inverse Square Potential

We develop an a posteriori error estimate of hierarchical type for Dirichlet eigenvalue problems of the form (−∆ + (c/r))ψ = λψ on bounded domains Ω, where r is the distance to the origin, which is assumed to be in Ω. This error estimate is proven to be asymptotically identical to the eigenvalue approximation error on a family of geometrically-graded meshes. Numerical experiments demonstrate this asymptotic exactness in practice.


Introduction
We consider the eigenvalue problem in a bounded domain Ω ⊂ R 2 with the Dirichlet boundary condition φ| ∂Ω = 0, where c > 0 and r = |x| is the distance to the origin, which is assumed to be in Ω. This eigenvalue problem is associated with the Schrödinger equation with the Dirichlet boundary condition u| ∂Ω = 0. For simplicity of the theory, we assume that c is constant, and that Ω is a polygon. However, we will also consider non-polygonal domains in the examples as well as in the experiments throughout the paper, to which our results extend. In fact, our analysis can also deal with certain type of non-constant functions c possessing multiple inverse-square singularities, as in [15].
The eigenvalue problem (1) with the inverse-square, or centrifugal, potential (c/r) 2 is of importance in quantum mechanics (for example cf. [23,12,11,28]). This potential presents the same "differential order" as the Laplacian near the origin, as is apparent when the Laplacian is expressed in polar coordinates. The strong singularity r −2 in the potential generally causes singular behavior (unbounded gradient) in the solution of (2) as well as in some of the eigenfunctions of (1). In addition to the singular potential, the geometry (smoothness) of the domain and boundary conditions may also play a critical role in determining the regularity of the solution. Therefore, new analytical tools, different from techniques for standard elliptic operators with bounded coefficients, are needed to develop well-posedness and regularity results, as well as effective numerical algorithms for (2) and (1). For Schrödinger operators with similar singular potentials, the analysis is generally carried out in Sobolev spaces with special weights, instead of in the usual Sobolev space H m (for example, see [10,11,15,21,22] and references therein). In particular, based on the weighted estimates, effective finite element methods associated with a class of graded meshes were proposed in [21] to approximate singular solutions of the Schrödinger equation at the optimal rate. An a posteriori error estimate of hierarchical type for these optimal finite element algorithms was developed in [22], and it provides a practical stopping criterion for approximating the solution of (2). The present paper builds on this work, adapting it to eigenvalue problems. We prove, and then numerically demonstrate, that our cheaply-computable error estimate is asymptotically identical to the error in our eigenvalue approximation, independent of singularities present in the eigenfunctions or whether the eigenvalues are degenerate.
Finite element methods for elliptic eigenvalue problems are nearly as old those for the associated boundary value problems, so there is a rich literature, and basic analysis is well-developed, at least for standard  second-order elliptic operators. We do not attempt a comprehensive overview of the relevant literature, but merely cite two classic references for the basic theory, [3,25], and mention some three recent papers concerning a posteriori error estimates for lower-order methods which might be most readily compared to our own, [7,24,6]. In both [24] and [6], eigenvalue error estimates are developed for standard elliptic operators. These are proven to asymptotically exact under certain assumptions on mesh structure and smoothness of the eigenfunctions. Non-self-adjoint problems having real eigenvalues are also considered in [24]. The work [6] employs hierarchical bases for error estimation, in the same manner as we do here, but the effectivity analysis for boundary value problems is quite different from ours, as is the theoretical bridge between boundary value problems and eigenvalue problems-which is done here via a key identity (Lemma 3.1). A certain class of non-linear eigenvalue problems, also relevant in certain quantum physical applications, is considered in [7]. Asymptotic exactness of the eigenvalue errors is not considered in [7], and cannot be achieved for the type of error estimates used, but the important issue of proving convergence of the associated adaptive method is addressed.
The rest of the paper is organized as follows. In Section 2, we define the weighted Sobolev spaces used in the analysis of [21,22], state key regularity results, and present basic eigenvalue theory for (1). Two examples are presented to provide some intuition about these eigenvalue problems, and one of these is revisited explicitly in the experiments. In Section 3, we first formulate the finite element approximation of the eigenvalue problem on graded meshes (Definition 3.2). Then, using finite element analysis in weighted spaces, we prove the exactness of our a posteriori error estimate in Theorem 3.6, our main result. In Section 4, we report numerical tests for different domains with different singular eigenfunctions. These tests confirm our theoretical prediction on the effectivity of the a posteriori estimate.

Basic Definitions and Results
Throughout, we use the following notation for the L 2 -inner-product and norm, For multi-indices α = (α 1 , α 2 ) ∈ N 2 0 , we employ the standard conventions |α| = α 1 +α 2 , and for . Let Q consist of the origin and the corners of Ω. These are the points at which one might expect an eigenfunction of (1) to have an unbounded gradient (cf. [1,2,13,14,27,17,18,16,20,8,4]). For x ∈ Ω, let ρ(x) be the distance between x and Q. We define the following weighted Sobolev spaces and their corresponding norms We note that K 0 0 = L 2 (Ω). Letting we define the following bilinear form on H, and note that it is, in fact, an inner-product. We denote the induced norm by ||| · |||. It can be shown (cf. [22]) that Lemma 2.1. The norms ||| · ||| and · K 1 1 are equivalent on H. With these definitions in hand, the variational form of our eigenvalue problem is given by We will refer to a solution (λ, φ) of (8) as an eigenpair of B on H, with eigenvalue λ and eigenfunction φ. Before stating a few basic facts about the eigenvalue problem (8), we introduce a related family of boundary value problems and remark on their well-posedness. Lemma 2.1 leads to the well-posedness of (9) in H by the Riesz Representation Theorem.
A stronger regularity result is proven for the boundary value problem (9) in [21,Theorem 3.3]: There is a constant η > 0 depending only on Ω and the constant c ≥ 0 in (7) such that, for any f ∈ K m−1 a−1 , where m ∈ N 0 and |a| < η, we have u(f ) ∈ K m+1 a+1 . More specifically, it holds that where C depends on m and a, but not on f .
Let K : L 2 (Ω) → L 2 (Ω) be the solution operator associated with (9), i.e. Kf = u(f ). Theorem 2.2 implies that K is bounded. The fact that H 1 0 (Ω) is compactly embedded in L 2 (Ω) easily implies that H ⊂ H 1 0 (Ω) is compactly embedded in L 2 (Ω), because any bounded sequence in H is clearly bounded in H 1 0 (Ω) as well. Therefore, K is a self-adjoint compact operator on L 2 (Ω), and we have the following basic results for the eigenvalue problem associated with its inverse -the operator defined by the bilinear form B on H: (1) The eigenvalues of (8) form a sequence of positive numbers with no finite accumulation points. We will assume that they are ordered (2) It is possible to choose a corresponding sequence of eigenfunctions φ n , i.e.
Furthermore these eigenfunctions form a Hilbert basis for L 2 (Ω). (3) There is a min-max variational characterization of the eigenvalues Such results are standard for symmetric elliptic problems with L ∞ coefficients. Given an eigenvalue λ, we denote its invariant subspace by In other words, E(λ) consists of all the eigenfunctions associated with λ, as well as the zero function. As indicated above by the non-strict inequalities λ i ≤ λ i+1 , eigenvalues may be degenerate, having geometric multiplicity dim E(λ) > 1.
For bounded domains, it is straight-forward to see that K m a ⊂ K m b if b < a. Let η = η(Ω, c) be as in Theorem 2.2, and choose |a| < η. Let (λ, φ) be an eigenpair of (8) for any j ≤ 2n. We therefore have the following corollary. Corollary 2.3. There is a constant η > 0 depending only on Ω and the constant c ≥ 0 in (7) such that, for any |a| < η and any eigenfunction φ, it holds that φ ∈ K n a+1 for all n ≥ 0; more briefly, φ ∈ K ∞ a+1 . We close this section with two examples which help provide some intuition about these eigenvalue problems-particularly the types of singularities which can occur in the eigenfunctions.
Example 2.4. Suppose Ω is the unit disk, r < 1. Expressing the eigenvalue problem in polar coordinates and using separation-of-variables, we find the eigenvalues λ mn and corresponding invariant subspaces E(λ mn ) for n ≥ 0 and m ≥ 1, where j m (ν) is the m th positive root of the first-kind Bessel function J ν (z), and σ n = √ n 2 + c 2 . When n = 0 these subspaces are one-dimensional. These formulas hold for c ≥ 0, but we will primarily be interested in the case c ∈ (0, 1).  Table 2. The smallest eight eigenvalues for the unit sector problem, Example 2.5, listed together with their indices and multiplicities for c = 0; α = 1/2 (left) and α = 2/3 (right).
have asymptotic behavior r c near the origin for this problem, so the gradient will be unbounded at the origin if c ∈ (0, 1). If n ≥ 1, the gradient of an eigenfunction in E(λ mn ) vanishes at the origin. Determining the location in the spectrum of all eigenfunctions having a specific regularity is untenable as it would require knowledge of the interlacing of roots of the Bessel functions J σn . However, a couple of specific instances will shed light on typical behavior. Table 1 gives the smallest eight eigenvalues (counting multiplicities) for this problem when c = 1/2 and c = 2/3. These eigenvalues are correct in all digits shown, up to rounding in the last digit.
Example 2.5. Fixing α ≥ 1/2, suppose Ω is the sector of the unit disk, with r < 1 and 0 < θ < π/α, where θ is the opening angle of the sector. The limiting case α = 1/2 represents the unit disk with the positive x-axis removed. As before, we find the eigenvalues λ mn and corresponding invariant subspaces E(λ mn ) for m, n ≥ 1, where and σ n = (nα) 2 + c 2 . Again, these formulas hold for c ≥ 0, and the case c = 0 (the Laplace eigenvalue problem) illustrates the type of singularities which can occur solely because of re-entrant corners, i.e. α < 1. We provide the first eight eigenvalues for c = 0, with α = 1/2 (slit disk) and α = 2/3 (L-shape) in Table 2.

Discretization and Error Estimation
In this section, we consider the finite element approximation of solutions to the eigenvalue problem (8), with focus on the estimation of error in the computed eigenvalue approximations. Before getting into the details of our finite element discretization, we make a few relevant claims which hold more generally. We restrict our attention to finite dimensional subspaces V ⊂ H. The natural analogues of (8) and (9) are As before, we will refer to a solution (λ,φ) of (12) as an eigenpair of B on V , with eigenvalueλ and eigenfunctionφ. These discrete problems are well-posed by basic linear algebra and by the coercivity of the bilinear form on V (Lemma 2.1). More specifically, if {v 1 , v 2 , . . . , v N } is a basis for V , then (12) is equivalent to the generalized eigenvalue problem where the matrices A = (a ij ) and M = (m ij ). The following analogues from the continuous eigenvalue problem apply: (1) There are precisely N = dim(V ) eigenvalues for the system (14), which we take to be ordered as 0 <λ 1 ≤λ 2 ≤ · · · ≤λ N .
(2) It is possible to choose corresponding eigenfunctionsφ n , i.e. B(φ n , v) =λ n (φ n , v) for all v ∈ V , such that (φ i ,φ j ) = δ ij . These eigenfunctions clearly form a Hilbert basis for V . (3) There is a min-max variational characterization of the eigenvalueŝ This characterization implies thatλ n ≥ λ n for n ≤ N .
Proof. Using the fact that φ = φ = 1, we first have the following identities, Subtracting these identities, we obtain the well-known error formula, Using the fact that B(φ, v) =λ(φ, v) for all v ∈ H, we further manipulate (16), from which (15) follows directly.
Our computed estimate of λ −λ will be based on approximating |||φ −φ||| 2 , treating the rest of the bound in (15) as higher-order terms. To make this more precise, we now shift to definitions of our finite element spaces, and a few key results.
Given a triangulation T of Ω, let V be the vertex set (the vertices of all triangles), which we assume includes all singular points Q. We define the two spaces The space P k consists of polynomials of total degree k or less. We note that it is necessary that v(0) = 0 for v ∈ V . We will approximate the solution of (8) in the space V and assess the error of this approximation in the space W .  . Let T be a triangulation of Ω whose vertices include Q, such that no triangle in T has more than one of its vertices in Q. For κ ∈ (0, 1/2], a κ refinement of T , denoted by κ(T ), is obtained by dividing each edge AB of T in two parts as follows: • If neither A nor B is in Q, then we divide AB into two equal parts.
• Otherwise, if say A is in Q, we divide AB into AC and CB such that |AC| = κ|AB|. This will divide each triangle of T into four triangles. Figure 1 shows a triangle having a singular vertex (the vertex on the top), together with three subsequent κ-refinements, with κ = 1/4. Given an initial triangulation T 0 , the associated family of graded triangulations {T n : n ≥ 0} is defined recursively, T n+1 = κ(T n ).
Remark 3.3. Although it may be useful in practice to have a different grading ratio κ q for each q ∈ Q, which is not difficult to implement, we do not pursue that theoretical generality here.
Given a family {T n } of κ-refined triangulations, we set V n = V (T n ) and W n = W (T n ), and u n (f ) ∈ V n is the solution of (13) on V n . We see that dim(V n ) ∼ dim(W n ) ∼ 4 n , because each refinement increases the number of triangles by precisely a factor of 4. We also define ε n (f ) ∈ W n by We collect two key results from [21,22]. Theorem 3.4. Let η be as in Theorem 2.2, and for 0 < a < min(η, 1) choose κ = 2 −1/a . There is a constant C which is independent of f and n, such that where the related numbers σ > 1 and a < ξ < min(2a, η, 1) are also independent of f and n.
Although we generally think of f as remaining fixed, these results allow for f to vary with n, and we exploit this fact below. In what follows, we let {(λ k,n , φ k,n ) : 1 ≤ k ≤ N = dim(V n )}, be eigenpairs for B on V n with (φ i,n , φ j,n ) = δ ij , as discussed at the beginning of this section; therefore, λ k,n =λ k , φ k,n =φ k and V = V n , for example. Corollary 3.5. Setting ψ k,n = u(λ k,n φ k,n ) and ε k,n = ε n (λ k,n φ k,n ), and employing the assumptions of Theorem 3.4, we have where C is independent of k and n.
We emphasize that the computation of ε k,n involves solving the problem Using the standard bases for V n and W n , it is shown in [22,Theorem 3.6] that the condition number of the matrix associated with (24) is well-conditioned independent of n. In fact, it is spectrally equivalent to its own diagonal. This makes (24) very inexpensive to solve, particularly when compared with computing solutions to the eigenvalue problem Ax =λM x on V n , where the stiffness matrix A is known to have a condition number which grows like 4 n .
Suppose we fix k and consider the sequence of discrete eigenpairs {(λ k,n , φ k,n ) : n ≥ 0} for the Schrödinger operator with meshes appropriately graded near the set Q (see Theorem 3.4). We do not rehearse standard finite element convergence theory for eigenvalue problems (cf. [9, Section 3.3], [3]), but by the approximation property given in Theorem 3.4, the following results for our eigenvalue problem can be derived in a similar fashion: • The approximate eigenvalues λ k,n converge (down) to λ k quadratically, • The distance between φ k,n and the invariant subspace E(λ k ) associated with λ k decreases linearly in the energy norm, We emphasize that λ k may be a degenerate eigenvalue (repeated in the sequence of eigenvalues), so E(λ k ) may have dimension greater than one. The analogous statements on V n hold as well. In light of this, it does not necessarily make sense to say that {φ k,n : n ≥ 0} converges, even up to sign. Nevertheless, we do have convergence in the sense of (26), and we refer to this (loosely) as eigenvector "convergence". We also remark that, although the eigenfunction v ∈ E(λ k ) which is nearest to φ k,n may not be of unit length, but we do not lose (26) if we add this restriction.
In practice, the eigenvalue convergence is precisely quadratic, and the eigenvector "convergence" is precisely linear on these properly graded meshes.
Combining these pieces, we arrive at our key eigenvalue error theorem.

Numerical Experiments
In this section we report the outcome of several numerical experiments, to demonstrate how well the theory of previous sections-particularly Theorem 3.6-are realized in practice. The data of interest are the eigenvalue errors λ k,n − λ k , their computed estimates |||ε k,n ||| 2 , and the associated effectivities The software package PLTMG [5] was used for these experiments, with suitable modifications for employing hierarchical basis error estimation and graded mesh refinement, with ARPACK [19] in shift-and-invert mode as the algebraic eigenvalue solver. In order reduce the width of tables of numerical data, we use the following abbreviation of scientific notation, a × 10 m ↔ a m . For example, We first revisit the unit disk problem of Example 2.4, considering case c = 1/2, for which we know that the eigenfunctions associated with λ 1 ≈ 9.86960440 and λ 6 ≈ 39.4784176 have an r 1/2 singularity (see Table 1). The grading ratio κ = 0.2 was used for refinement. Note that by Theorem 3.4, the upper bound of the grading parameter κ near the origin is 2 −1/(1/2) = 0.25 to achieve the optimal convergence rate. Therefore, we have chosen an appropriate grading ratio here. The data for these experiments are in Table 3. The eigenvalue convergence is seen to be quadratic, i.e. linear in N = dim(V n ), and the effectivities are very close 1. The top row of data is absent for λ 6 because, on this coarse mesh, the approximate eigenvalue 33.2876671 was actually (slightly) nearer to λ 4 = λ 5 ≈ 27.1817273 than to λ 6 . The effectivity of the error estimate when this was taken into account was 1.010.
We now consider the degenerate eigenvalue λ = λ 2 = λ 3 ≈ 15.9205134. The corresponding invariant subspace (eigenspace) is spanned by Of course the ordering of φ 2 and φ 3 is arbitrary, as is that particular choice of basis for this invariant subspace. These functions are smooth enough to be optimally approximated on a sequence of uniformly refined meshes, κ = 0.5; for comparison, grading ratio κ = 0.4 and κ = 0.2 were used as well. On each mesh, two approximate eigenvalues and eigenvectors were computed and error estimates for both were computed.
The results indicate that it really is irrelevant which of the approximate eigenpairs is used to estimate error in the eigenvalue approximation, as indicated by our theory. Since the code (PLTMG+ARPACK) assigns an order to the approximate eigenpairs, we employ this order as well, (λ 2,n , φ 2,n ), (λ 3,n , φ 3,n ). The computed eigenvalues λ 2,n and λ 3,n agreed with each other to far more digits than they agreed with λ 2 = λ 3 , so the reported errors are identical. It is only the error estimates, and hence effectivities, which are slightly different.
In terms of the grading, all three grading choices gave optimal order convergence, as the theory predicts, with uniform refinement (κ = 0.5) yielding the smallest errors and κ = 0.2 yielding the largest errors. In terms of effectivities, uniform refinement was the worst, followed in order by κ = 0.4 and κ = 0.2, though all were close to 1. To save space, only the data for κ = 0.5 and κ = 0.2 are reported in Table 4. In order to n N λ 2,n − λ 2 |||ε 2,n ||| 2 EF F λ 3,n − λ 3 |||ε 3,n ||| 2 EF F 0 48 1.019 +0 9.  Table 4. Data for the Unit Disk problem, corresponding to approximations of λ 2 = λ 3 . Uniform refinement (top) and κ = 0.2 graded refinement (bottom), c = 1/2. demonstrate the "drift" in approximate eigenfunctions associated with degenerate eigenvalues, we provide a sequence of contour plots for φ 2,n , n = 1, 2, 3, 4 in Figure 2. The contour plots of φ 3,n are essentially obtained by rotating the given plots by 90 degrees. These plots illustrate the assertion in Section 3 that the sequence {φ k,n } may not converge, though the terms are getting successively closer to E(λ k ).
Finally, we turn to the L-shape domain Ω = (−3, 3) 2 \ [1, 3) 3 . We consider the case c = 1/2, for which there will be eigenfunctions having an r 1/2 -singularity at the origin and an r 2/3 -singularity at the point (1, 1). We use the grading ratios κ = 0.2 for triangles touching the origin, and κ = 0.3 for triangles touching (1, 1). Table 5 contains our approximations and error estimates for the first four eigenvalues. Contour plots of the first four eigenfunctions are given in Figure 3. As an interesting comparison, we also consider the case c = 0, for which no singularity is present at the origin, and grading is only needed near the point (1, 1). The eigenvalues in this case have been obtained elsewhere to very high accuracy [26] using a computational very well-suited to the Laplacian, and we report their values here, rescaling them by a factor of four due to the fact that our domain has four times the area of theirs: λ 1 ≈ 2.4099310 , λ 2 ≈ 3.7993130 , λ 3 = π 2 2 ≈ 4.9348022 , λ 4 ≈ 7.3803703 .