An efficient finite element method and error analysis for eigenvalue problem of Schrödinger equation with an inverse square potential on spherical domain

We provide an effective finite element method to solve the Schrödinger eigenvalue problem with an inverse potential on a spherical domain. To overcome the difficulties caused by the singularities of coefficients, we introduce spherical coordinate transformation and transfer the singularities from the interior of the domain to its boundary. Then by using orthogonal properties of spherical harmonic functions and variable separation technique we transform the original problem into a series of one-dimensional eigenvalue problems. We further introduce some suitable Sobolev spaces and derive the weak form and an efficient discrete scheme. Combining with the spectral theory of Babuška and Osborn for self-adjoint positive definite eigenvalue problems, we obtain error estimates of approximation eigenvalues and eigenvectors. Finally, we provide some numerical examples to show the efficiency and accuracy of the algorithm.


Introduction
The Schrödinger eigenvalue problem with the inverse-square (IS) or centrifugal potential is widely used in nuclear physics, quantum physics, nonlinear optics, and so on [1][2][3][4][5][6]. The potential in many electronic equations produces singularity and can describe the attraction or repulsion between objects, which usually leads to strong singularities of the eigenfunctions, and this cannot be simply regarded as a perturbation term [7][8][9][10][11]. Therefore new tools and techniques, different from linear elliptic operators with bounded coefficients, are urgently needed for such an eigenvalue problem with the IS potential, both in analysis and in numerics. Ghanbari et al. [12][13][14][15] and Khater et al. [16][17][18][19][20] discussed some effective numerical methods to remove the singularity of nonlocal operators. In addition, some other works [21][22][23][24][25][26] mainly focus on studying exact solitary wave solutions.
In the last decade, increasing attention is paid to the numerical approximation of Schrödinger operators with similar singular potentials [1,[27][28][29][30]. If we solve the Schrödinger eigenvalue problem directly in a three-dimensional domain, it will take a lot of computing time and memory capacity to obtain high-precision numerical solutions [11,31,32]. In practice, we usually need to solve the Schrödinger eigenvalue problem on a three-dimensional spherical domain. As far as we know, there is little work discussing the eigenvalue problem with IS potential on a spherical domain.
The purpose of this work is developing an effective finite element method for the eigen- The rest of the paper is organized as follows. In Sect. 2, we obtain a dimensional reduction scheme for the Schrödinger eigenvalue problem with IS potential. In Sect. 3, we construct a weak form and numerical discretization scheme. In Sect. 4, we obtain error estimates of approximation eigenvalues and eigenfunctions. In Sect. 5, we study the implementation details of the algorithm. In Sect. 6, we test the accuracy and convergence of the numerical algorithm. Finally, in Sect. 7, we provide some concluding remarks.

Dimension reduction scheme
In this paper, we consider the following eigenvalue problem with Dirichlet boundary conditions: where c is a bounded constant, and = {(x, y, z) ∈ R 3 : 0 ≤ a < r < b} with r = x 2 + y 2 + z 2 .
Let x = r sin θ cos φ, y = r sin θ sin φ, z = r cos θ , ψ(r, θ , φ) = u(x, y, z). We take the Laplacian in spherical coordinates Then (1)-(2) can be rewritten as the eigenvalue problem of the Schrödinger equation in spherical coordinates: We recall that an important property of the spherical harmonics {Y l m } (as normalized in [33]) is that they are eigenfunctions of the Laplace-Beltrami operator S . More precisely, and from the definition of L 2 (S) we find Using spherical harmonic expansion, we obtain that First of all, we take into account the case a = 0. Define the differential operatorL l ψ m l = 1 r 2 (∂ r (r 2 ∂ r ψ m l )l(l + 1)ψ m l ). By substituting expression (10) into (5) and (6) we can obtain a series of equivalent one-dimensional eigenvalue problems: Let Then (11)-(12) can be rewritten as Analogously, for the case a > 0, inserting expression (10) into (5) and (7), we can derive the following one-dimensional eigenvalue problem:

Weak form and discrete scheme
When a = 0, we can define the space The inner product and norm can be defined as follows: Thus the weak form of (13)- (14) is: where a l u m l , v m l = We denote byṼ h (l) =P 1h ∩ H 1 0,l (I) the approximation space, whereP 1h is a piecewise linear interpolation polynomial space. Thus the corresponding numerical scheme of (19) is: Find When a > 0, the usual space can be denoted as We define the inner product and norm as Thus the weak form of (17)-(18) is: where We denote by V h = P 1h ∩ H 1 0 (I) the approximation space, where P 1h is a piecewise linear interpolation polynomial space. Then the corresponding numerical scheme of (21) is: Find

Error estimation of approximation solutions
In this section, we prove error estimates of approximate eigenvalues and eigenfunctions. Without loss of generality, we only consider the case a > 0.
Use the technique of [34], we deduce the following results. Proof Using the Cauchy-Schwarz inequality, we find Let V (λ l ) and V (λ lh ) be the eigenfunction spaces corresponding to the eigenvalues λ l and λ lh , respectively. Let where u m l a l = a l (u m l , u m l ). From Theorems 1 and 2 we know that a l (u m l , v m l ) (resp., b l (u m l , v m l )) is a symmetric, continuous, and coercive bilinear form on H 1 0 (I) × H 1 0 (I) (resp., L 2 (I) × L 2 (I)).
By the spectral theory of eigenvalue problems [34] we have the following theorem.
Theorem 3 Let (λ l , u m l ) and (λ lh , u m lh ) be respectively the eigenpairs of (21) and (22). Then the following inequalities hold: Define the piecewise linear interpolation operator I h : where p li (t) denotes the linear interpolation polynomial of u m l on the interval Then from an error formula of linear interpolating remainder term we derive that where ξ i (t) ∈ I i is a function depending on t.
For u m l , we have the following error results.
, u m l ∈ H 1 0 (I). Suppose that u m l is smooth enough such that |∂ Thus we obtain From the Poincaré inequality we derive that This ends the proof.
We can get the following standard error estimation results.
Theorem 5 Let (λ l , u m l ) and (λ lh , u m lh ) be the eigenpairs of (21) and (22), respectively. If u m l ∈ H 1 0 (I) satisfy the condition of Theorem 4, then the following inequalities hold: Proof Since Combining this with Theorem 3, we obtain (28).

Implementation of the numerical scheme
In this section, we present the algorithm to solve problems (20) and (22). First, we define some basis functions. Let where i = 1, . . . , N -1. We find that Substituting expression (29) into (20) and noticing the v m lh , we can observe the linear system where A = (a ij ), From the properties of the basis functions we know that the stiff matrices and mass matrix in (30) are all tridiagonal sparse matrices.
Case 2. a > 0. Let Substituting expression (31) into (22) and taking v m lh in V h , we obtain the linear system where Similarly, from the properties of the basis functions we know that the stiff matrices and mass matrix in (32) are all tridiagonal sparse matrices.
Remark 1 Our numerical method can transform three-dimensional problems into a series of eigenvalue problems. By constructing appropriate basis functions these onedimensional value problems will be discretized into asparse stiffness matrix and mass matrix, which can be efficiently solved.

Numerical experiments
In this section, we present several numerical examples to check the convergence and accuracy of the numerical algorithm.
Example 1 We take c = 0, a = 0, b = 1, and l = 0, 1, 2 as our example. In Tables 1-3, we provide the first four eigenvalues with different l and h of Example 1.
Example 2 We take c = 1 3 , a = 0, b = 1, and l = 0, 1, 2 as our example. In Tables 1-3, we give the first four eigenvalues with different l and h of Example 2.
We know from Tables 1-6 that the approximation eigenvalues achieve at least threedigit accuracy with h ≤ 1 512 for l = 0, 1, 2. To further show the convergence of approximation eigenvalues, we let the numerical solution of h = 1 1024 be the reference solution. The error of the approximate eigenvalues with different h are provided in Figs. 1-6.
We observe from Figs. 1-6 that the numerical eigenvalues are convergent.      Example 3 We take c = 1 2 , a = 1, b = 2, and l = 0, 1, 2 as our example. The first four eigenvalues for different l and h are listed in Tables 7-9.
We know from Tables 7-9 that the numerical eigenvalues can achieve six-digit accuracy at least for h ≤ 1 512 for l = 0, 1, 2. Similarly, we select the solutions with h = 1 1024 as the reference solutions. In Figs. 7-8, we plot the error of the approximate eigenvalues. We observe from Figs. 7-8 that the numerical eigenvalues are also convergent.

Conclusions
In this work, an efficient finite element method is constructed to solve the Schrödinger eigenvalue problem with the IS potential on a spherical domain. By using spherical coordinate transformation and variable separation technique, we reduce the original problem into a series of equivalent and independent eigenvalue problems, which not only overcomes the difficulty brought by the singular coefficient, but also reduces the degrees of freedom greatly by dimension reduction. Thus we can spend less computing time and memory capacity to obtain high-precision numerical solutions. Numerical results show that our algorithm is very effective. We believe that this method can be extended to more complex practical problems.